The blog Why Add Q to MC?, explained the advantages of carefully chosen, low discrepancy sampling sites for approximating multivariate integrals, or equivalently, expectations of functions of multivariate random variables. This blog post explains at a Bayesian approach to determining the sample size required to satisfy the user’s error tolerance.
Recall that the problem of interest takes the following form:
\begin{align}
\mu: &= \int_{[0,1]^d} f(\boldsymbol{x}) \, \text{d}{\boldsymbol{x}} = {\mathbb{E}}[f(\boldsymbol{x})], \; \text{where } \boldsymbol{x} \sim \mathcal{U}[0,1]^d, \\
\mu \approx \widehat{\mu}_n &: = \frac 1n \sum_{i=1}^n f(\boldsymbol{x}_i),
\text{where } \boldsymbol{x}_1, \boldsymbol{x}_2, \ldots \text{are the sampling sites or nodes}.
\end{align}
The question of when enough samples have been taken to satisfy the specified error criterion, $\left \lvert{\mu – \widehat{\mu}_n}\right \rvert \le \varepsilon$, is an important one. Stopping criteria based on Bayesian credible intervals are one good way of answering this question.
Traditionally, Bayesian credible intervals take $\mathcal{O}(n^3)$ operations to construct, where $n$ is the sample size. This is due to the need to invert an $n \times n$ (covariance) matrix. In contrast, computing $\widehat{\mu}$ requires only $\mathcal{O}(n)$ operations. Fortuitously, there is a way to reduce the computational cost of constructing the Bayesian credible interval to $\mathcal{O}(n \log (n))$ when using lattice or digital net low discrepancy samples. These Bayesian stopping criteria are implemented as CubBayesLatticeG and CubBayesNetG. This blog explains the key points.
The Bayesian Approach
The Bayesian approach to approximating integrals assumes that the integrand, $f:[0,1]^d \to \mathbb{R}$, is drawn from a Gaussian stochastic process parameterized by a (constant) mean and a covariance function defined by scale and shape parameters. This is denoted $f \sim \mathcal{GP}(m,s^2 C_{\boldsymbol{\theta}})$, and means that
\begin{align*}
{\mathbb{E}}[f(\boldsymbol{x})] = m, \; {\textup{cov}}\bigl(f(\boldsymbol{t}),f(\boldsymbol{x}) \bigr) = {\mathbb{E}}[(f(\boldsymbol{t}) – m)(f(\boldsymbol{x}) – m)] = s^2 C_{\boldsymbol{\theta}}(\boldsymbol{t},\boldsymbol{x}), \; \forall \boldsymbol{t}, \boldsymbol{x} \in [0,1]^d.
\end{align*}
Here, $m$, $s$, and $C_{\boldsymbol{\theta}}$ must be specified or estimated. Because $f$ is a Gaussian process, the multivariate integral of $f$, i.e., $\mu$, has a Gaussian distribution. Automatic Bayesian cubature uses these assumptions to increase $n$ until we reach
$
\mathbb{P}_f [\left\lvert { \mu – \widehat{\mu}_n} \right\rvert \le \varepsilon ] \ge 99\%.
$
The Bayesian credible interval, which depends on the sampling nodes, and the parameters defining the Gaussian process is
\begin{equation*}
\mathbb{P}_f \left[
|\mu-\widehat{\mu}_n| \leq \text{err}_\text{CI}
\right] = 99\%, \;
{\text{err}_{\text{CI}}} := 2.58 \; s \; \sqrt{c_{0,{\boldsymbol{\theta}}} – {\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X})}^T \mathsf{C}_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X})^{-1}\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X})}
\end{equation*}
assuming that the mean of the Gaussian process, $m$, is zero. The formulae for the case of general $m$ are more complicated, but similar. Here,
\begin{gather*}
\label{eqn:fGaussDist_c0}
\mathsf{X} = ({\boldsymbol{x}}_1, \ldots, {\boldsymbol{x}}_n)^T, \qquad c_{0,{\boldsymbol{\theta}}} := \int_{[0,1]^{d} \times [0,1]^{d}} C_{\boldsymbol{\theta}} (\boldsymbol{t},\boldsymbol{x}) \, \text{d}{\boldsymbol{t}} \, \text{d}{\boldsymbol{x}} , \\
\boldsymbol{c}_{\boldsymbol{\theta}}(\mathsf{X}) := \left( \int_{[0,1]^d} C_{\boldsymbol{\theta}} (\boldsymbol{t},\boldsymbol{x}_i) \, \text{d}{\boldsymbol{t}} \right)_{i=1}^n, \qquad \mathsf{C}_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X}) = \left( C_{\boldsymbol{\theta}}({\boldsymbol{x}}_i,{\boldsymbol{x}}_j) \right)_{i,j=1}^n,
\end{gather*}
and $2.58$ represents $99.5\%$ quantile of the standard Gaussian distribution.
Bayesian cubature proceeds by increasing $n$ until $\text{err}_{\text{CI}}$ is no greater than the error tolerance.
The Bayesian credible interval as a stopping criterion assumes that the integrand, $f$, is not an outlier with respect to this Gaussian process. To improve the chances that this is the case, the hyperparameters $m$, $s$, and $\boldsymbol{\theta}$ should be tuned using the function data. In particular, $s$ should depend on the magnitude of the fluctuations of $f$. The standard approaches for tuning the hyperparameters require optimization, for which each iteration involves an eigendecomposition of $ C_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X})$.
Fast Bayesian Cubature
The computation of $\text{err}_{\text{CI}}$ for a single $\mathsf{X}$ and $\boldsymbol{\theta}$ requires inversion of the covariance matrix $ C_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X})$, which typically costs $\mathcal{O}(n^3)$. An eigendecompostion of $ C_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X})$ costs at least this much. The way to overcome this computational burden is to use families of covariance kernels, $C_{\boldsymbol{\theta}}$, that match the design, $\mathsf{X}$, in a way that the eigenvectors and eigenvalues of the matrix $C_{\boldsymbol{\theta}}(\mathsf{X},\mathsf{X})$ can be obtained via fast transforms. For integration lattices and kernels $C_{\boldsymbol{\theta}}$ such as
\begin{align*}
%\label{eq:shInv}
C_\boldsymbol{\theta}(\boldsymbol{t},\boldsymbol{x}) = K_\boldsymbol{\theta}(\boldsymbol{t} – \boldsymbol{x} \bmod \boldsymbol{1}),
\end{align*} where
\begin{align*}
\label{the_kernel_eqn_bernoulli}
K_\boldsymbol{\theta}(\boldsymbol{x}) =
\prod_{l=1}^d \biggl[
1 – (-1)^{r} \eta B_{2r}( x_l) \biggr], \;
\forall \boldsymbol{t},\boldsymbol{x} \in [0,1]^d, \ \boldsymbol{\theta} = (r,\eta), \ r \in \mathbb{N}, \ \eta > 0,
\end{align*}
and $B_{2r}$ is an even degree Bernoulli polynomial,
the fast transform corresponds to a fast Fourier transform. For digital nets (see our blog Digital Sequences, the Niederreiter Construction) and $C_\boldsymbol{\theta}$ that are digitally-shift invariant, the fast transform corresponds to a fast Walsh transform. In both cases, the computational burden attributable to tuning the hyperparameters and computing the width of the credible interval is a reasonable $\mathcal{O}(n \log(n))$.
The following code shows how to use a Bayesian stopping criterion for Keister’s example.
>>> from qmcpy import *
>>> from numpy import Keister, Lattice, CubBayesLatticeG
>>>
>>> tol = .005
>>>
>>> integrand = Keister(Lattice(dimension=2, order='linear'))
>>> keister_2d_exact = integrand.exact_integ()
>>> solution, data = CubBayesLatticeG(integrand, abs_tol=tol, n_init=2 ** 5).integrate()
>>> print('Integration error: ', abs(solution - keister_2d_exact) < tol)
Listing 1: Example usage of Lattice Bayesian Cubature algorithm.
This example can be run in Google Colab without any installation using this notebook.
References
- Hickernell, F. J. Blog: Why Add Q to MC? https://qmcpy.org/2020/06/25/why_add_q_to_mc/ (2020).
- Choi, S.-C. T., Hickernell, F., McCourt, M. & Sorokin, A. QMCPy: A quasi-Monte Carlo Python Library https://qmcsoftware.github.io/QMCSoftware/. 2020.
- Jagadeeswaran, R. & Hickernell, F. Fast Automatic Bayesian Cubature Using Lattice Sampling. Stat. Comput. 29, 1215–1229 (2019).
- Sloan, I. H. & Joe, S. Lattice Methods for Multiple Integration (Oxford University Press, Oxford, 1994).
- Hickernell, F. & Niederreiter, H. The Existence of Good Extensible Rank-1 Lattices. J. Complexity 19, 286–300 (2003).
- Cooley, J. W. & Tukey, J. W. An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation 19, 297–301 (1965).
- Ebert, A. Blog: Digital Sequences, the Niederreiter Construction. https://qmcpy.org/2021/06/04/digital-sequences-the-niederreiter-construction/ (2021).
- Jagadeeswaran, R. Fast Automatic Bayesian Cubature Using Matching Kernels and Designs PhD thesis (Illinois Institute of Technology, 2019).
- Fino, B. J. & Algazi, V. R. Unified Matrix Treatment of the Fast Walsh-Hadamard Transform. IEEE Transactions on Computers C-25, 1142–1146 (1976).
Fred Hickernell is Professor of Applied Mathematics and Vice Provost for Research at Illinois Institute of Technology. His research spans computational mathematics and statistics.