A QMCPy Quick Start

In this tutorial, we introduce QMCPy [1] by an example. QMCPy can be installed with pip install qmcpy or cloned from the QMCSoftware GitHub repository.

Consider the problem of integrating the Keister function [2] with respect to a d-dimensional Gaussian measure:

f(\boldsymbol{x}) = \pi^{d/2} \cos(||\boldsymbol{x}||), \qquad \boldsymbol{x} \in \mathbb{R}^d, \qquad \boldsymbol{X} \sim \mathcal{N}(\boldsymbol{0}_d,\mathsf{I}_d/2),
\\ \mu  =  \mathbb{E}[f(\boldsymbol{X})] := \int_{\mathbb{R}^d} f(\boldsymbol{x}) \, \pi^{-d/2} \exp( - ||\boldsymbol{x}||^2) \,  \rm d \boldsymbol{x}
\\     =  \int_{[0,1]^d} \pi^{d/2}  \cos\left(\sqrt{ \frac 12 \sum_{j=1}^d\Phi^{-1}(x_j)}\right)  \, \rm d \boldsymbol{x},

where ||\boldsymbol{x}|| is the Euclidean norm, \mathsf{I}_d is the d-dimensional identity matrix, and \Phi denotes the standard normal cumulative distribution function. When d=2, \mu \approx 1.80819 and we can visualize the Keister function and realizations of the sampling points depending on the tolerance values, \varepsilon, in the following figure:

Keister Function

Keister Function

The Keister function is implemented below with help from NumPy [3] in the following code snippet:

import numpy as np
def keister(x):
    """
    x: nxd numpy ndarray
       n samples
       d dimensions

    returns n-vector of the Kesiter function
    evaluated at the n input samples
    """
    d = x.shape[1]
    norm_x = np.sqrt((x**2).sum(1))
    k = np.pi**(d/2) * np.cos(norm_x)
    return k # size n vector

In addition to our Keister integrand and Gaussian true measure, we must select a discrete distribution, and a stopping criterion [4]. The stopping criterion determines the number of points at which to evaluate the integrand in order for the mean approximation to be accurate within a user-specified error tolerance, \varepsilon. The discrete distribution determines the sites at which the integrand is evaluated.

For this Keister example, we select the lattice sequence as the discrete distribution and corresponding cubature-based stopping criterion [5]. The discrete distribution, true measure, integrand, and stopping criterion are then constructed within the QMCPy framework below.

import qmcpy
d = 2
discrete_distrib = qmcpy.Lattice(dimension = d)
true_measure = qmcpy.Gaussian(discrete_distrib, mean = 0, covariance = 1/2)
integrand = qmcpy.CustomFun(true_measure,keister)
stopping_criterion = qmcpy.CubQMCLatticeG(integrand = integrand, abs_tol = 1e-3)

Calling integrate on the stopping_criterion instance returns the numerical solution and a data object. Printing the data object will provide a neat summary of the integration problem. For details of the output fields, refer to the online, searchable QMCPy Documentation at https://qmcpy.readthedocs.io/.

solution, data = stopping_criterion.integrate()
print(data)
LDTransformData (AccumulateData Object)
    solution        1.808
    comb_bound_low  1.808
    comb_bound_high 1.809
    comb_flags      1
    n_total         2^(13)
    n               2^(13)
    time_integrate  0.017
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
CustomFun (Integrand Object)
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
Lattice (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       1
    order           natural
    entropy         273562359450377681412227949180408652150
    spawn_key       ()

This guide is not meant to be exhaustive but rather a quick introduction to the QMCPy framework and syntax. In an upcoming blog, we will take a closer look at low-discrepancy sequences such as the lattice sequence from the above example.

References

  1. Choi, S.-C. T., Hickernell, F., McCourt, M., Rathinavel J., & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library. https://qmcsoftware.github.io/QMCSoftware/. 2020.

  2. Keister, B. D. Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996).

  3. Oliphant, T., Guide to NumPy https://ecs.wgtn.ac.nz/foswiki/pub/Support/ManualPagesAndDocumentation/numpybook.pdf (Trelgol Publishing USA, 2006).

  4. Hickernell, F., Choi, S.-C. T., Jiang, L. & Jimenez Rugama, L. A. in WileyStatsRef-Statistics Reference Online (eds Davidian, M.et al.) (John Wiley & Sons Ltd., 2018).

  5. Jimenez Rugama, L. A. & Hickernell, F. Adaptive Multidimensional Inte-gration Based on Rank-1 Lattices in Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (eds Cools, R. & Nuyens, D.) 163.arXiv:1411.1966 (Springer-Verlag, Berlin, 2016), 407–422.