[Paper Review] GAIL---Guaranteed Automatic Integration Library in MATLAB: Documentation for Version 2.1

Sou‐Cheng T. Choi, Yuhan Ding+6 more|arXiv (Cornell University)|Mar 23, 2015

Numerical Methods and Algorithms7 references5 citations

TL;DR

GAIL v2.1 is an open-source MATLAB library providing guaranteed, adaptive numerical integration for univariate and multivariate functions using Monte Carlo and Quasi-Monte Carlo methods with rigorous error bounds. It ensures accuracy by leveraging cone conditions on Walsh-Fourier coefficients, offering guaranteed convergence under specified tolerance criteria for integrands in a defined function class.

ABSTRACT

Automatic and adaptive approximation, optimization, or integration of functions in a cone with guarantee of accuracy is a relatively new paradigm. Our purpose is to create an open-source MATLAB package, Guaranteed Automatic Integration Library (GAIL), following the philosophy of reproducible research and sustainable practices of robust scientific software development. For our conviction that true scholarship in computational sciences are characterized by reliable reproducibility, we employ the best practices in mathematical research and software engineering known to us and available in MATLAB. This document describes the key features of functions in GAIL, which includes one-dimensional function approximation and minimization using linear splines, one-dimensional numerical integration using trapezoidal rule, and last but not least, mean estimation and multidimensional integration by Monte Carlo methods or Quasi Monte Carlo methods.

Motivation & Objective

• To develop a reproducible, open-source MATLAB library for guaranteed numerical integration with rigorous error control.
• To implement adaptive integration algorithms that automatically adjust sampling based on function smoothness and error tolerance.
• To support both standard Monte Carlo and quasi-Monte Carlo methods (Sobol’ and rank-1 lattice) for multidimensional integration.
• To ensure mathematical guarantees on integration accuracy under cone conditions on Fourier coefficients.
• To promote sustainable scientific software practices through documented, versioned, and reusable code for computational mathematics.

Proposed method

• Employs a cone condition on the Walsh-Fourier coefficients of the integrand to enable error estimation and guaranteed convergence.
• Uses adaptive refinement with Sobol’ sequences for quasi-Monte Carlo integration, dynamically increasing sample points until error tolerance is met.
• Implements generalized tolerance control via 'max' or 'comb' types, combining absolute and relative error tolerances with configurable weights.
• Applies fast Walsh transforms to estimate coefficient decay rates and validate cone membership, ensuring theoretical error bounds.
• Supports multiple integration domains including bounded intervals, hypercubes, and infinite domains with uniform or normal measures.
• Integrates with MATLAB’s built-in random number generators and provides a modular, extensible architecture for algorithm plug-in.

Experimental results

Research questions

• RQ1 How can numerical integration be made adaptive and guaranteed in accuracy without prior knowledge of the integrand’s smoothness?
• RQ2 What conditions on Fourier coefficients ensure rigorous error bounds in quasi-Monte Carlo integration?
• RQ3 How can a balance between absolute and relative error tolerances be systematically managed in adaptive integration?
• RQ4 To what extent can Sobol’ sequences and rank-1 lattices improve convergence rates for high-dimensional integrals with guaranteed error?
• RQ5 Can a software library be designed to ensure reproducibility and reliability in scientific computing through formal error guarantees?

Key findings

• The cubSobol_g function successfully estimated the integral of f(x) = x1*x2 over [0,1]^2 with a result of 0.2500, matching the analytical value within the specified tolerance.
• For f(x) = x1^2*x2^2*x3^2 over R^3 with normal measure, the method returned q = 1.0004, closely approximating the true value of 1.
• The algorithm achieved a guaranteed error bound for f(x) = exp(-x1^2 - x2^2) on [-1,2]^2, yielding q = 2.6532 with a tolerance of 1e-3.
• In a financial option pricing example, the expected payoff of a European call was estimated at q = 2.0552 with a tolerance of 1e-4, demonstrating robustness in real-world applications.
• The library correctly identified when the integrand lay outside the cone condition (via exitflag = 2), providing clear warnings when guarantees could not be assured.
• The use of fudge parameter @(m) 5*2.^-m in the cone condition enabled stable convergence across diverse test functions, including high-dimensional and non-smooth cases.