3

I am trying to integrate over some matrix entries in Python. I want to avoid loops, because my tasks includes 1 Mio simulations. I am looking for a specification that will efficiently solve my problem.

I get the following error: only size-1 arrays can be converted to Python scalars

from scipy import integrate
import numpy.random as npr

n = 1000
m = 30

x = npr.standard_normal([n, m])


def integrand(k):
    return k * x ** 2


integrate.quad(integrand, 0, 100)

This is a simplied example of my case. I have multiple nested functions, such that I cannot simple put x infront of the integral.

Improve this question
5
  • 1
    Why is integrand a function of k, but doesn't use it? But in any case quad integrates a scalar function. Commented Jul 19, 2019 at 2:34
  • There was a typo. I corrected for k but my problem is the same. I integrate over k. I want to get 1 Mio integrals. Commented Jul 19, 2019 at 8:29
  • If I understand correctly, each element of the matrix needs to be integrated independently. You are just asking for a possibly vectorized version of integrand = lambda k: k * x[i, j]**2 \n for i in range(n): \n \t for j in range(m): \n \t \t integrate.quad(integrand, 0, 100) (sorry for poor formatting) Commented Jul 19, 2019 at 9:09
  • Yes, correct. I am looking for a package/function that will do that without loops, because I have n = 1 000 000. Say you want to take the ln of x, that would simply be sp.ln(x) and it would again return a matrix, where each entry is ln(x[i,j]). I am looking for the equivalent function for integration. Commented Jul 19, 2019 at 9:23
  • Actually, take a look here: stackoverflow.com/a/41226624/5048010 Commented Jul 19, 2019 at 9:47

4 Answers 4

Reset to default
1

Well you might want to use parallel execution for this. It should be quite easy as long as you just want to execute integrate.quad 30000000 times. Just split your workload in little packages and give it to a threadpool. Of course the speedup is limited to the number of cores you have in your pc. I'm not a python programer but this should be possible. You can also increase epsabs and epsrel parameters in the quad function, depending on the implemetation this should speed up the programm as well. Of course you'll get a less precise result but this might be ok depending on your problem.

import threading
from scipy import integrate
import numpy.random as npr

n = 2
m = 3
x = npr.standard_normal([n,m])

def f(a):
    for j in range(m):
        integrand = lambda k: k * x[a,j]**2
        i =integrate.quad(integrand, 0, 100)
        print(i) ##write it to result array

for i in range(n):
    threading.Thread(target=f(i)).start();
##better split it up even more and give it to a threadpool to avoid
##overhead because of thread init
Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Although your answer might be correct your answer is pretty dry. Try adding some code snippets as those are often way more helpful. Adding to that, this question is specific to something called an only size-1 arrays can be converted to Python scalars and you're saying how to do it which, again, might be true but doesn't help with the actual question. Anyways I see that you are a New Contributor, I welcome you to the StackOverflow community and encourage you to help (and get help) with hight quality/informational posts :)
thanks for your feedback. I wrote a little snippet for you :). I don't see how the error relates to the problem.
1

2023 Update: quad_vec

As of 2023, there is a Scipy function integrate.quad_vec for efficient quadrature of vector functions.

A solution to the question is the following highly-vectorized procedure.

from scipy import integrate
import numpy as np

x = np.random.standard_normal([1000, 30])

def integrand(k):
    return k * x**2

res = integrate.quad_vec(integrand, 0, 100)

The output res[0] contains a 1000x30 matrix with the numerical integrals for every parameter x.

This function uses standard adaptive quadrature, like QUADPACK, however, the interval subdivision is the same for all components of the vector function. The subdivision is chosen to ensure that every component of the vector function satisfies the selected convergence criteria. This implies that it makes sense to use quad_vec only when the different components of the vector function have a qualitatively similar behaviour.

Improve this answer

2 Comments

Very very nice! Is there a simple explanation of why this enables parallelization? Naively each element in the vector function can behave quite differently.. Also the title reads 2003 rather than 2023 :)
@Histoscienology I addressed your question in my revision.
0

This is maybe not the ideal solution but it should help a bit. You can use numpy.vectorize. Even the doc says: The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop. But still, a %timeit on the simple example you provided shows a 2.3x speedup.

The implementation is 

from scipy import integrate
from numpy import vectorize
import numpy.random as npr  

n = 1000
m = 30

x = npr.standard_normal([n,m])

def g(x):
    integrand = lambda k: k * x**2
    return integrate.quad(integrand, 0, 100)


vg = vectorize(g)
res = vg(x)
Improve this answer

2 Comments

Thank you very much. I tried the run with n = 5000 and I have been waiting for 2 hours already and then I accidentally stopped the run. The solution is sadly not optional. Any additional suggestions are welcome.
As explained in stackoverflow.com/a/41226624/5048010, quad doesn't support vectorization due to the nature of the algorithm itself. You can use instead scipy.integrate.simps and it should work automatically.
0

quadpy (a project of mine, commercial) does vectorized quadrature:

import numpy
import numpy.random as npr
import quadpy


x = npr.standard_normal([1000, 30])


def integrand(k):
    return numpy.multiply.outer(x ** 2, k)


scheme = quadpy.line_segment.gauss_legendre(10)
val = scheme.integrate(integrand, [0, 100])

This is much faster than all other answers.

Improve this answer

1 Comment

This solution requires buying a license

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.