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$.

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

In [4]:

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

    returns n-vector of the Keister 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

import numpy as np def keister(x): """ x: nxd numpy ndarray n samples d dimensions returns n-vector of the Keister 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.

In [5]:

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)

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/.

In [6]:

solution, data = stopping_criterion.integrate()
print(data)

solution, data = stopping_criterion.integrate() print(data)

Data (Data)
    solution        1.808
    comb_bound_low  1.808
    comb_bound_high 1.809
    comb_bound_diff 0.001
    comb_flags      1
    n_total         2^(13)
    n               2^(13)
    time_integrate  0.013
CubQMCLatticeG (AbstractStoppingCriterion)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_limit         2^(30)
CustomFun (AbstractIntegrand)
Gaussian (AbstractTrueMeasure)
    mean            0
    covariance      2^(-1)
    decomp_type     PCA
Lattice (AbstractLDDiscreteDistribution)
    d               2^(1)
    replications    1
    randomize       SHIFT
    gen_vec_source  kuo.lattice-33002-1024-1048576.9125.txt
    order           NATURAL
    n_limit         2^(20)
    entropy         331737392988868539347762820736945757498

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 Integration 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.