Reputation: 31
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.
Upvotes: 3
Views: 2446
Reputation: 17266
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.
Upvotes: 1
Reputation: 58951
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.
Upvotes: 0
Reputation: 11
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
Upvotes: 1
Reputation: 1653
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)
Upvotes: 0