CODEWITHSUNDEEP
pythonmatrixscipynumerical-integration
cody_tastic
cody_tastic

Reputation: 31

Numerical integration for matrix values in Python

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

Answers (4)

divenex
divenex

Reputation: 17266

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.

Upvotes: 1

Nico Schlömer
Nico Schlömer

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

KMB
KMB

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

Gianluca Micchi
Gianluca Micchi

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

Related Questions