[SciPy-user] slow integrals

[David M. Cooke]

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!).

What do you use for your quadrature routine in your C version?

quad() uses QUADPACK, and calls your function once for each point, so it
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
> 
> [Xquad, Wquad] = real(scipy.special.orthogonal.ps_roots(Nquad))
> 
> 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
QUADPACK, though.

Also look at the infodict returned by quad(); quad_explain() describes
the stuff in there.

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