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:

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

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.
Keister, B. D. Multidimensional Quadrature Algorithms. Computers in Physics 10, 119–122 (1996).
Oliphant, T., Guide to NumPy https://ecs.wgtn.ac.nz/foswiki/pub/Support/ManualPagesAndDocumentation/numpybook.pdf (Trelgol Publishing USA, 2006).
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).
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.