contribution: code for multi-dimensional adaptive integration

[Steven G. Johnson]

Dear GSL folks,

For some time now, I've felt that GSL has cried out for a multidimensional 
cubature routine, similar to its 1d adaptive quadrature routines.  Now, 
I've written one (~700 lines in C) based on the well-known Genz-Malik 
embedded cubature rule (with help from a GPL'ed library called HIntLib, 
see below), and I hope that you're interested to have it included in GSL. 
I've attached the unmodified code; I'm willing to do the work to modify it 
for your coding style, etcetera.

Adaptive cubature is a huge win over competing methods in GSL for moderate 
dimensionality (dimensions 2-6), as I demonstrate with a couple benchmarks 
below, and these kinds of dimensionalities are obviously important in lots 
of scientific applications.

My Genz-Malik cubature is based in part on code from a GPL'ed library 
called HIntLib (http://www.cosy.sbg.ac.at/~rschuer/hintlib/).  HIntLib is 
a much more sophisticated, full-featured package that is focused on 
extremely high-dimensional integration, including many other cubature and 
Monte-Carlo methods.  Unfortunately, this means it is a large and complex 
library that is overkill for low-dimensional integrals, and is not 
suitable for inclusion in GSL (also, it is in C++).  My implementation 
strips out all the dozens of C++ classes defined in HIntLib and implements 
just the serial adaptive cubature with just the degree-7 Genz-Malik rule 
for dimensions > 1 (although it would be easy to add more embedded 
cubature rules), and I think it should be easy to include in GSL.

Actually, Genz and Malik wrote their own implementation of their adaptive 
cubature rule as a Fortran code adapt.f which can be found as part of the 
public-domain CMLIB distributed by NIST.  Interestingly, however, I found 
that HIntLib's (and my) implementation of the same cubature rule to often 
be more robust -- the Genz Malik code seems to be easily "fooled" 
(converges to the wrong result) for discontinuous and/or strongly peaked 
objectives, at least in my experience.  (These problems were the whole 
reason I switched from adapt.f to HIntLib in my own applications.)

The only other alternative that I know of is CUBPACK, which is non-free.

Cordially,
Steven G. Johnson

PS. Note that my attached code includes a 1d quadrature based on GSL's 
15-point QAGS.  This for completeness, since the Genz-Malik rule does not 
work in 1d.  Obviously, a GSL version of my routine would just call the 
your 1d integration routine when dimensions = 1.

------------------------------benchmarks-------------------------------- 
Currently, to perform multi-dimensional integration in GSL you have to 
either use Monte-Carlo integration or recursively call 1d adaptive 
integrations.  Both of these are extremely inefficient in many cases 
compared to adaptive cubature, and to show this I've benchmarked my 
routine against GSL for two simple test functions (whose integrals I know 
analytically):

cos: a simple smooth integrand = product of cos(x[i]) for dimensions i

gaussian: a gaussian integral that tests facility with strongly peaked
           integrands, with the axes rescaled via x -> (1-x)/x so that
           the integration goes from 0 to infinity.  The integrand
 	  (before coordinate rescaling) is:
 	        product of 2/sqrt(pi) * exp(-x[i]^2)
           (which integrates to 1).

The integration volume was from x[i] = 0 to 1+0.4*sin(i) for the cos 
integrand, and from 0 to 1 for the gaussian integrand (rescaled to 
0..infinity).

In each case, I integrated to a relative error tolerance of 1e-3, and 
counted the function evaluations. (For the Monte-Carlo algorithms, I 
repeatedly doubled the number of iterations until the estimated error 
achieved the given relative error.)  I used the Mersenne Twister RNG. 
For recursive quadrature, I used QAG with GSL_INTEG_GAUSS15.  The 
resulting counts for each method were as follows:

                ***** cos integrand (smooth) *****
#dimensions   cubature  recursive-QAG  MC-plain  MC-miser  MC-vegas
      2           17          225        163840    20462     2890
      3           33          3375       327680    40936     5120
      4           57         50625       327680    81885     12500
      5           93        759375       327680   163788     10240
      6          447       11390625         -     163794     10935
      7          1205     170859375      655360   163801     12800
      8          5213         -             -     327616     25600
      9         31185         -             -     327620     25600
     10        160605         -         1310720   655268     25600
     11        603693         -         1310720   655274     20480

(A "-" indicates that the "converged" value had an error bigger than the 
specified tolerance, or that, in the case of QAG, it was taking longer 
than I cared to wait.)

***** gaussian integrand (smooth, peaked, boundary singularity) *****
#dimensions   cubature  recursive-QAG  MC-plain  MC-miser  MC-vegas
      2          255        5085           -       40923      6250
      3         1353       337635       5242880    163759    49130
      4         8721     22314225      10485760    327564      -
      5        58869         -         20971520    655194    168070
      6       327055         -         20971520    1310458   390625
      7      2079589         -            -        2621013   409600
-------------- next part --------------
A non-text attachment was scrubbed...
Name: integrator.c
Type: text/x-csrc
Size: 26177 bytes
Desc: 
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20050610/e484581f/attachment.bin>