A QMCPy Quick Start

Posted byAleksei Sorokin and Sou-Cheng ChoiPosted inBlog

We have created QMCPy, an open-source, object-oriented quasi- Monte Carlo (QMC) software framework in Python 3, which contains standardized parent classes for modeling integrands, measure, discrete distribution, and stopping criteria.  We hope QMCPy could enable researchers and users to quickly extend and experiment with novel algorithmic components that validate their theories, or to simply leverage or compare proven ones as they develop their applications.

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,\boldsymbol{I}_d/2),\] \[\begin{aligned} \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{1}{2} \sum_{j=1}^d\Phi^{-1}(x_j)}\right) \, \rm d \boldsymbol{x} \end{aligned}\]

where $||\boldsymbol{x}||$ is the Euclidean norm, $\boldsymbol{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 with n samples d dimensions 
    returns n-vector of the Kesiter function evaluations
    """
    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
lattice = qmcpy.Lattice(dimension = 2)
gaussian = qmcpy.Gaussian(lattice, mean = 0, covariance = 1/2)
cf_keister = qmcpy.CustomFun(gaussian, custom_fun = keister)
stopping_criterion = qmcpy.CubQMCLatticeG(cf_keister, abs_tol = 1e-4)

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)
Solution: 1.8082         
CustomFun (Integrand Object)
Lattice (DiscreteDistribution Object)
    dimension       2^(1)
    randomize       1
    seed            None
    backend         gail
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      2^(-1)
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         1.00e-04
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(16)
    solution        1.808
    r_lag           2^(2)
    time_integrate  0.070

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. J., 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. J., 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. J., 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.

Aleksei Sorokin
Website | + posts

Aleksei Sorokin is an Applied Mathematics and Data Science student at the Illinois Institute of Technology.

Sou-Cheng Choi
+ posts

Dr. Sou-Cheng T. Choi is Chief Data Scientist at the Kamakura Corporation and Research Associate Professor of Applied Mathematics at the Illinois Institute of Technology.

%d bloggers like this: