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.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
cody_tastic
  • 31
  • 3
  • 1
    Why is `integrand` a function of `k`, but doesn't use it? But in any case `quad` integrates a scalar function. – hpaulj Jul 19 '19 at 02: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. – cody_tastic Jul 19 '19 at 08: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) – Gianluca Micchi Jul 19 '19 at 09: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. – cody_tastic Jul 19 '19 at 09:23
  • Actually, take a look here: https://stackoverflow.com/a/41226624/5048010 – Gianluca Micchi Jul 19 '19 at 09:47

4 Answers4

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
KMB
  • 11
  • 2
  • 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 :) – Elias Jul 20 '19 at 11:05
  • thanks for your feedback. I wrote a little snippet for you :). I don't see how the error relates to the problem. – KMB Jul 20 '19 at 15:33
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)
Gianluca Micchi
  • 1,584
  • 15
  • 32
  • 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. – cody_tastic Jul 19 '19 at 16:07
  • 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. – Gianluca Micchi Jul 23 '19 at 08:16
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.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249
0

2003 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.

divenex
  • 15,176
  • 9
  • 55
  • 55