Subject: contribution: code for multi-dimensional adaptive integration
Date: Thu, 9 Jun 2005 21:13:55 -0400 (EDT)
From: "Steven G. Johnson" <stevenj@fftw.org>
Reply-To: stevenj@math.mit.edu
To: gsl-discuss@sources.redhat.com, help-gsl@gnu.org

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

    ---------------------------------------------------------------------
                   Name: integrator.c
   integrator.c    Type: NOTEPAD File
                         (application/x-unknown-content-type-cfile)
               Encoding: BASE64
