Cuba - a library for multidimensional numerical integration

The Cuba library offers a choice of four independent routines for multidimensional numerical integration: Vegas, Suave, Divonne, and Cuhre. They work by very different methods, summarized in the following table:

Routine Basic integration method Algorithm type Variance reduction
Vegas Sobol quasi-random sample
or Mersenne Twister pseudo-random sample		 Monte Carlo
Monte Carlo		 importance sampling
Suave Sobol quasi-random sample
or Mersenne Twister pseudo-random sample		 Monte Carlo
Monte Carlo		 globally adaptive subdivision
    + importance sampling
Divonne Korobov quasi-random sample
or Sobol quasi-random sample
or Mersenne Twister pseudo-random sample
or cubature rules		 Monte Carlo
Monte Carlo
Monte Carlo
deterministic		 stratified sampling,
    aided by methods from
    numerical optimization
Cuhre cubature rules deterministic globally adaptive subdivision

All four have a C/C++, Fortran, and Mathematica interface and can integrate vector integrands. Their invocation is very similar, so it is easy to substitute one method by another for cross-checking. For further safeguarding, the output is supplemented by a chi-square probability which quantifies the reliability of the error estimate.

The source code compiles with gcc, the GNU C compiler. The C functions can be called from Fortran directly, so there is no need for adapter code. Similarly, linking Fortran code with the library is straightforward and requires no extra tools.

Download:

Windows users: Cuba 3 uses fork(2) to parallelize the execution threads. This POSIX function is not part of the Windows API, however, and is furthermore used in an essential way such that it cannot be worked around simply with CreateProcess etc. The only feasible emulation seems to be available through Cygwin.

Ready-made MathLink executables (Version 2.1, statically linked as far as possible):

Linux x86 (32-bit):

Mac OS X (universal binaries):

Windows (32-bit):

Vegas is the simplest of the four. It uses importance sampling for variance reduction, but is only in some cases competitive in terms of the number of samples needed to reach a prescribed accuracy. Nevertheless, it has a few improvements over the original algorithm and comes in handy for cross-checking the results of other methods.

Suave is a new algorithm which combines the advantages of two popular methods: importance sampling as done by Vegas and subregion sampling in a manner similar to Miser. By dividing into subregions, Suave manages to a certain extent to get around Vegas' difficulty to adapt its weight function to structures not aligned with the coordinate axes.

Divonne is a further development of the CERNLIB routine D151. Divonne works by stratified sampling, where the partitioning of the integration region is aided by methods from numerical optimization. A number of improvements have been added to this algorithm, the most significant being the possibility to supply knowledge about the integrand. Narrow peaks in particular are difficult to find without sampling very many points, especially in high dimensions. Often the exact or approximate location of such peaks is known from analytic considerations, however, and with such hints the desired accuracy can be reached with far fewer points.

Cuhre employs a cubature rule for subregion estimation in a globally adaptive subdivision scheme. It is hence a deterministic, not a Monte Carlo method. In each iteration, the subregion with the largest error is halved along the axis where the integrand has the largest fourth difference. Cuhre is quite powerful in moderate dimensions, and is usually the only viable method to obtain high precision, say relative accuracies much below 1e-3.

Upward of 75% of all questions regarding Cuba have to do with how to choose bounds different from the unit hypercube in Fortran, C, and C++. The solution is not to choose bounds but to scale the integrand. For the mathematically challenged, the explicit transformation is (in one dimension)

∫ab dx f[x] → ∫01 dy f[a + (b - a) y] (b - a),

where the final (b - a) is the one-dimensional version of the Jacobian. This generalizes straightforwardly to more than one dimension.

For constant integration bounds, this transformation might be implemented in Fortran as 

        subroutine ScaledIntegrand(ndim, x, ncomp, result)
        implicit none
        integer ndim, ncomp
        double precision x(ndim), result(ncomp)

        integer maxdim
        parameter (maxdim = 16)

        double precision upper(maxdim)
        common /ubound/ upper

        double precision lower(maxdim)
        common /lbound/ lower

        integer dim, comp
        double precision range, jacobian, scaledx(maxdim)

        jacobian = 1
        do dim = 1, ndim
          range = upper(dim) - lower(dim)
          jacobian = jacobian*range
          scaledx(dim) = lower(dim) + x(dim)*range
        enddo

        call Integrand(ndim, scaledx, ncomp, result)

        do comp = 1, ncomp
          result(comp) = result(comp)*jacobian
        enddo
        end

This site and the programs offered here are not commercial. Cuba is an open-source package and free of charge. If you want to use Cuba in a commercial application, make sure you understand the GNU Lesser General Public License under which Cuba is distributed. Cuba is being developed at the Max Planck Institute for Physics in Munich.

Last update: 4 May 13 Thomas Hahn