[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!).
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
More information about the SciPy-user
mailing list