# [SciPy-user] slow integrals

David M. Cooke cookedm at physics.mcmaster.ca
Wed Jun 7 15:40:03 CDT 2006

```On Wed, Jun 07, 2006 at 12:18:29PM +0200, Eric Emsellem wrote:
> Hi!
>
> I am performing a double integration using scipy.integrate.quad and it
> seems rather slow in fact so I would like to know if there is any way to
> make it more efficient (I have a pure C version of that integral and it
> is much faster, maybe by a factor of more than 10!).

doesn't take advantage of speedup from using arrays.

> I am doing this integral as:
>
> result = scipy.integrate.quad(IntlosMu1, -Inf, Inf,...)
>
> so an integral between - and + infinity.

That's always fun. Note that quadrature() will use a Fourier integral
for infinite limits. I can suggest the usual tricks:

- use a weighting function. The quad_explain() function has the info
on what's usuable
- map to a finite interval using some transform (arctan, say)

> The Integrand (IntlosMu1) is itself an integral which I am computing
> using a direct Quadrature (it is an integral on an adimensional variable
> T which varies between 0 and 1).
> So I compute once and for all the abs and weight for the quadrature using
>
>
> and then I pass on Xquad and Wquad as parameters to
> scipy.integrate.quad, to be used in the integrand.

You might want to try dblquad; it does essentially what you're doing above,
but does the inner integral adaptively (AFAIK), so it may use fewer points.

> (One last note: I am computing this double integral many many times on
> different points so I really need efficiency here...)
>
> So the questions I have are:
> - can you already see something wrong with this in terms of efficiency?
> (I doubt it since I don't provide much info, but just in case)
> - are there other integration scheme I could use to do that ?

I would suggest integrate.quadrature, which does a simple Gaussian
integration scheme (you'll need to do it on a finite interval).

Looking at quadrature(), though, it looks like it wraps the passed function
in an inefficient value-at-a-time wrapper (vec_func). Try removing that
wrapper from the source, so it will call your function with an array of
points to evaluate. Or raise a TypeError in your function if passed a scalar.

> - how should I try to test things and see where the bottleneck is?

Profiling, I guess. It may not pick up the callbacks to your function from

the stuff in there.

--
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke                      http://arbutus.physics.mcmaster.ca/dmc/
|cookedm at physics.mcmaster.ca

```