Pre-integration via active subspaces

For background on QMC see the monograph [7] and the survey article [6]. For a survey of RQMC see [21]. We will emphasize scrambled net integration with a recent description in [32]. Here we give a brief account.

QMC provides a way to estimate $\mu=\int_{[0,1]^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$ with greater accuracy than can be done by MC while sharing with MC the ability to handle larger dimensions $d$ than can be well handled by classical quadrature methods such as those in [5]. The QMC estimate, like the MC one takes the form $\hat{\mu}_{n}=(1/n)\sum_{i=0}^{n-1}f(\boldsymbol{x}_{i})$, except that instead of $\boldsymbol{x}_{i}\,{\buildrel\text{iid}\over{\sim}\,}\mathbb{U}[0,1]^{d}$ the sample points are chosen strategically to get a small value for

which is known as the star discrepancy. The Koksma-Hlawka inequality (see [15]) is

where $\|\cdot\|_{\mathrm{HK}}$ denotes total variation in the sense of Hardy and Krause. It is possible to construct points $\boldsymbol{x}_{i}$ with $D_{n}^{*}=O((\log n)^{d-1}/n)$ or to choose an infinite sequence of them along which $D_{n}^{*}=O((\log n)^{d}/n)$. Both of these are often written as $O(n^{-1+\epsilon})$ for any $\epsilon>0$ and this rate translates directly to $|\hat{\mu}_{n}-\mu|$ when $f\in\mathrm{BVHK}[0,1]^{d}$, the set of functions of bounded variation in the sense of Hardy and Krause. While the logarithmic powers are not negligible they correspond to worst case integrands that are not representative of integrands evaluated in practice, and as we see below, some RQMC methods provide a measure of control against them.

The QMC methods we study are digital nets and sequences. To define them, for an integer base $b\geqslant$ let

for $\boldsymbol{k}=(k_{1},\dots,k_{d})$ and $\boldsymbol{c}=(c_{1},\dots,c_{d})$ where $k_{j}\in\mathbb{N}_{0}$ and $0\leqslant a_{j}<b^{k_{j}}$. The sets $E(\boldsymbol{k},\boldsymbol{c})$ are called elementary intervals in base $b$. They have volume $b^{-|\boldsymbol{k}|}$ where $|\boldsymbol{k}|=\sum_{j=1}^{d}k_{j}$. For integers $m\geqslant t\geqslant 0$, the points $\boldsymbol{x}_{0},\dots,\boldsymbol{x}_{n-1}$ with $n=b^{m}$ are a $(t,m,d)$-net in base $b$ if every elementary interval $E(\boldsymbol{k},\boldsymbol{c})$ with $|\boldsymbol{k}|\leqslant m-t$ contains $b^{m-|\boldsymbol{k}|}$ of those points, which is exactly $n$ times the volume of $E(\boldsymbol{k},\boldsymbol{c})$.

For each $\boldsymbol{k}$ the sets $E(\boldsymbol{k},\boldsymbol{c})$ partition $[0,1)^{d}$ into $b^{|\boldsymbol{k}|}$ congruent half open sets. If $|\boldsymbol{k}|\leqslant m-t$, then the $n$ points $\boldsymbol{x}_{i}$ are perfectly stratified over that partition. The power of digital nets is that the points $\boldsymbol{x}_{i}$ satisfy ${m-t+d-1\choose d-1}$ such stratifications simultaneously. They attain $D_{n}^{*}=O((\log n)^{d-1}/n)$ after approximating each $[\boldsymbol{0},\boldsymbol{a})$ by sets $E(\boldsymbol{k},\boldsymbol{c})$ [24]. Given $b$ and $m$ and $d$, a smaller $t$ is the better. It is not always possible to get $t=0$.

For an integer $t\geqslant 0$, the infinite sequence $(\boldsymbol{x}_{i})_{i\geqslant 0}$ is a $(t,s)$-sequence in base $b$ if for all integers $m\geqslant t$ and $r\geqslant 0$, the points $\boldsymbol{x}_{rb^{m}},\dots,\boldsymbol{x}_{rb^{m}+b^{m}-1}$ are a $(t,m,d)$-net in base $b$. This means that the first $b^{m}$ points are a $(t,m,d)$-net as are the second $b^{m}$ points and if we take the first $b$ such nets together we get a $(t,m+1,d)$-net. Similarly the first $b$ such $(t,m+1,d)$-nets form a $(t,m+2,d)$-net and so on ad infinitum.

The nets and sequences in common use are called digital nets and sequences owing to a digital strategy used in their construction, where $x_{ij}$ is expanded in base $b$ and there are rules for constructing those base $b$ digits from the base $b$ representation of $i$. See [24]. The most commonly used nets and sequences are those of Sobol’ [39]. Sobol’ sequences are $(t,d)$-sequences in base $2$ and sometimes using base $2$ brings computational advantages for computers that work in base $2$. The value of $t$ is nondecreasing in $d$. The first $2^{m}$ points of a $(t,d)$-sequence can be a $(t^{\prime},m,d)$-net for some $t^{\prime}<t$.

While the Koksma-Hlawka inequality (1) shows that QMC is asymptotically better than MC for $f\in\mathrm{BVHK}[0,1]^{d}$ it is not usable for error estimation, and furthermore, many integrands of interest such as unbounded ones have $V_{\mathrm{HK}}(f)=\infty$. RQMC methods can address both problems. In RQMC the points $\boldsymbol{x}_{i}\sim\mathbb{U}[0,1]^{d}$ individually while collectively they have a small $D_{n}^{*}$. The uniformity property makes

an unbiased estimate of $\mu$ when $f\in L^{1}[0,1]^{d}$. If $f\in L^{2}[0,1]^{d}$ then $\mathrm{Var}(\hat{\mu}_{n})<\infty$ and it can be estimated by using independent repeated RQMC evaluations.

Scrambled $(t,m,d)$-nets have $\boldsymbol{x}_{i}\sim\mathbb{U}[0,1]^{d}$ individually and the collective condition is that they form a $(t,m,d)$-net with probability one. For an infinite sequence of $(t,m,d)$-nets one can use scrambled $(t,d)$-sequences. See [26] for both of these. The estimate $\hat{\mu}_{n}$ taken over a scrambled $(t,d)$-sequence satisfies a strong law of large numbers if $f\in L^{2}[0,1]^{1+\epsilon}$ [32]. If $f\in L^{2}[0,1]^{d}$, then $\mathrm{Var}(\hat{\mu}_{n})=o(1/n)$ giving the method asymptotically unbounded efficiency versus MC which has variance $\sigma^{2}/n$ for $\sigma^{2}=\mathrm{Var}(f(\boldsymbol{x}))$. For smooth enough $f$, $\mathrm{Var}(\hat{\mu}_{n})=O(n^{-3}(\log n)^{d-1})$ [28, 30] with the sharpest sufficient condition in [46]. Also there exists $\Gamma<\infty$ with $\mathrm{Var}(\hat{\mu}_{n})\leqslant\Gamma\sigma^{2}/n$ for all $f\in L^{2}[0,1]^{d}$ [29]. This bound involves no powers of $\log(n)$.

For $j\in 1{:}d$, $\int_{[0,1]^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\int_{[0,1]^{d-1}}\int_{0}^{1}f(\boldsymbol{x})\,\mathrm{d}x_{j}\,\mathrm{d}\boldsymbol{x}_{-j}$ which we can also write as $\mathbb{E}(f(\boldsymbol{x}))=\mathbb{E}(\mathbb{E}(f(\boldsymbol{x})\!\mid\!\boldsymbol{x}_{-j}))$ for $\boldsymbol{x}\sim\mathbb{U}[0,1]^{d}$. For $\boldsymbol{x}\in[0,1]^{d}$, define

It simplifies the presentation of some of our results, especially Theorem 3.2 on variance reduction, to keep $g$ defined as above on $[0,1]^{d}$ even though it does not depend at all on $x_{j}$ and could be written as a function of $\boldsymbol{x}_{-j}\in[0,1]^{d-1}$ instead. In pre-integration

which as we noted in the introduction is conditional MC except that we now use RQMC inputs.

Pre-integration can bring some advantages for RQMC. The plain MC variance of $g(\boldsymbol{x})$ is no larger than that of $f(\boldsymbol{x})$ and is generally smaller unless $f$ does not depend at all on $x_{j}$. Thus the bound $\Gamma\sigma^{2}/n$ is reduced. Next, the pre-integrated integrand $g$ can be much smoother than $f$ and (R)QMC improves on MC by exploiting smoothness. For example [12, 44] observe that for some option pricing integrands, pre-integrating certain variable can remove the discontinuities in the integrand or its gradient.

The integrands we consider here are defined with respect to a Gaussian random variable. We are interested in $\mu=\mathbb{E}(f(\boldsymbol{z}))$ for $\boldsymbol{z}\sim\mathcal{N}(\boldsymbol{0},\Sigma)$ and a nonsingular covariance $\Sigma\in^{d\times d}$. Letting $R_{0}\in^{d\times d}$ with $R_{0}R_{0}^{\mathsf{T}}=\Sigma$ we can write $\mu=\mathbb{E}(f(R_{0}\boldsymbol{z}))$ for $\boldsymbol{z}\sim\mathcal{N}(0,I)$. For an orthogonal matrix $Q\in^{d\times d}$ we also have $Q\boldsymbol{z}\sim\mathcal{N}(0,I)$. Then taking $\boldsymbol{z}=\Phi^{-1}(\boldsymbol{x})$ componentwise leads us to the estimate

for RQMC points $\boldsymbol{x}_{i}$. The choice of $Q$ or equivalently $R$ does not affect the MC variance of $\hat{\mu}$ but it can change the RQMC variance. We will consider some examples later. The mapping $\Phi^{-1}$ from $\mathbb{U}[0,1]^{d}$ to $\mathcal{N}(0,I)$ can be replaced by another one such as the Box-Muller transformation. The choice of transformation does not affect the MC variance but does affect the RQMC variance. Most researchers use $\Phi^{-1}$ but [25] advocates for Box-Muller.

When we are using pre-integration for a problem defined with respect to a $\mathcal{N}(0,\Sigma)$ random variable we must choose $R$ and then the coordinate $j$ over which to pre-integrate. Our approach is to choose $R$ to make coordinate $j=1$ as important as we can, using active subspaces.

For $f\in L^{2}[0,1]^{d}$ we can define an analysis of variance (ANOVA) decomposition from [16, 40, 8]. For details see [31, Appendix A.6]. This decomposition writes

where $f_{u}$ depends on $\boldsymbol{x}$ only through $x_{j}$ with $j\in u$ and also $\int_{0}^{1}f_{u}(\boldsymbol{x})\,\mathrm{d}x_{j}=0$ whenever $j\in u$. The decomposition is orthogonal in that $\int_{[0,1]^{d}}f_{u}(\boldsymbol{x})f_{v}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=0$ if $u\neq v$. The term $f_{\varnothing}$ is the constant function everywhere equal to $\mu=\int_{[0,1]^{d}}f(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$. To each effect $f_{u}$ there corresponds a variance component

For $|u|\geqslant 2$, the effect $f_{u}$ is called a $|u|$-fold interaction. The variance components sum to $\sigma^{2}=\mathrm{Var}(f(\boldsymbol{x}))$. We will use the ANOVA decomposition below when describing how to choose a pre-integration variable.

Sobol’ indices [41] are derived from the ANOVA decomposition. For $u\subseteq 1{:}d$ these are

They provide two ways to judge the importance of the set of variables $x_{j}$ for $j\in u$. They’re usually normalized by $\sigma^{2}$ to get an interpretation as a proportion of variance explained.

The mean dimension of $f$ is $\nu(f)=\sum_{u\subseteq 1:d}|u|\sigma^{2}_{u}/\sigma^{2}$. It satisfies $\nu(f)=\sum_{j=1}^{d}\overline{\tau}^{2}_{j}/\sigma^{2}$. A ridge function takes the form $f(\boldsymbol{x})=h(\Theta^{\mathsf{T}}\boldsymbol{x})$ for $\Theta\in^{d\times r}$ with $\Theta^{\mathsf{T}}\Theta=I$. For $\boldsymbol{x}\sim\mathcal{N}(0,I)$, the variance of $f$ does not depend on $d$ and the mean dimension is $O(1)$ as $d\to\infty$ [17] if $h$ is Lipschitz. If $h$ has a step discontinuity then it is possible to have $\nu(f)=\Omega(\sqrt{d})$ reduced to $O(1)$ by pre-integration over a component variable $x_{j}$ with $\theta_{j}$ bounded away from zero as $d\to\infty$ [17].

To get variance formulas for scrambled nets we follow Dick and Pillichshammer [7] who work with a Walsh function expansion of $L^{2}[0,1)^{d}$ for which they credit Pirsic [36]. Let $\omega_{b}=e^{2\pi i/b}$ with $i$ being the imaginary unit. For $k\in\mathbb{N}_{0}$ write $k=\sum_{j\geqslant 0}\kappa_{j}b^{j}$ for base $b$ digits $\kappa_{j}\in\{0,1,\dots,b-1\}$. For $x\in[0,1)$ write $x=\sum_{k\geqslant 1}\xi_{j}b^{-j}$ for base $b$ digits $\xi_{j}\in\{0,1,\dots,b-1\}$. Countably many $x$ have an expansion terminating in either infinitely many $0$s or in infinitely many $b-1$s. For those we always choose the expansion terminating in $0$s.

Using the above notation we can define the $k$-th $b$-adic Walsh function ${}_{b}\text{wal}_{k}:[0,1)\to{\mathbb{C}}$ as

The summation in the exponent is finite because $k<\infty$. Note that ${}_{b}\mathrm{wal}_{0}(x)=1$ for all $x\in[0,1)$. For $\boldsymbol{x}=(x_{1},\dots,x_{d})\in[0,1)^{d}$ and $\boldsymbol{k}=(k_{1},\ldots,k_{d})\in\mathbb{N}_{0}^{d}$, the $d$-dimensional Walsh functions are defined as

The Walsh series expansion of $f(x)$ is

The $d$-dimensional $b$-adic Walsh function system is a complete orthonormal basis in $L_{2}([0,1)^{d})$ [7, Theorem A.11] and the series expansion converges to $f$ in $L^{2}$.

While our integrand is real valued, it will also satisfy an expansion written in terms of complex numbers. For real valued $f$,

The variance under scrambled nets is different. To study it we group the Walsh coefficients. For $\boldsymbol{\ell}\in\mathbb{N}_{0}^{d}$ let

Then define

The functions $\beta_{\boldsymbol{\ell}}$ are orthogonal in that $\int_{[0,1)^{d}}\beta_{\boldsymbol{\ell}}(\boldsymbol{x})\macc@depth\char 1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{}\,\mathrm{d}\boldsymbol{x}=0$ when $\boldsymbol{\ell}^{\prime}\neq\boldsymbol{\ell}$. For $\boldsymbol{\ell}\neq\boldsymbol{0}$, $\beta_{\boldsymbol{\ell}}(\boldsymbol{x})$ has variance

If $\boldsymbol{x}_{i}$ are a scrambled version of original points $\boldsymbol{a}_{i}\in[0,1)^{d}$ then under this scrambling

for a collection of gain coefficients $\Gamma_{\boldsymbol{\ell}}\geqslant 0$ that depend on the $\boldsymbol{a}_{i}$. This expression can also be obtained through a base $b$ Haar wavelet decomposition [27]. Our $\Gamma_{\boldsymbol{\ell}}$ equals $nG_{\boldsymbol{\ell}}$ from [7]. The variance of $\hat{\mu}$ under IID MC sampling is $(1/n)\sum_{\boldsymbol{\ell}\in\mathbb{N}^{d}_{0}\setminus\{\boldsymbol{0}\}}\sigma^{2}_{\boldsymbol{\ell}}$ so $\Gamma_{\boldsymbol{\ell}}<1$ corresponds to integrating the term $\beta_{\boldsymbol{\ell}}(\boldsymbol{x})$ with less variance than MC does.

If scrambling of [26] or [22] is applied to $\boldsymbol{a}_{i}$ then

This holds for any $\boldsymbol{a}_{i}$ not just digital nets. When $\boldsymbol{a}_{i}$ are the first $b^{m}$ points of a $(t,d)$-sequence in base $b$ then $\Gamma=\sup_{\boldsymbol{\ell}}\Gamma_{\boldsymbol{\ell}}<\infty$ (uniformly in $m$) so that $\mathrm{Var}(\hat{\mu})\leqslant\Gamma\sigma^{2}/n$. Similarly for any $\boldsymbol{\ell}\in\mathbb{N}_{0}^{d}$ we have $\Gamma_{\boldsymbol{\ell}}\to 0$ as $n=b^{m}\to\infty$ in a $(t,d)$-sequence in base $b$ from which $\mathrm{Var}(\hat{\mu}_{n})=o(1/n)$. For a $(t,m,s)$-net in base $b$, one can show that the gain coefficients $\Gamma_{\boldsymbol{\ell}}=0$ for all $\boldsymbol{\ell}$ with $|\boldsymbol{\ell}|\leqslant m-t$.

Theorem 3.2 shows that pre-integration does not increase the variance under scrambling. This holds whether or not the underlying points are a digital net, though of course the main case of interest is for scrambling of digital nets and sequences.

Pre-integration has another benefit that is not captured by Theorem 3.2. By reducing the input dimension from $d$ to $d-1$ we might be able to find better points in $[0,1]^{d-1}$ than the $(t,m,d)$-net we would otherwise use in $[0,1]^{d}$. Those improved points might have some smaller gain coefficients or they might be a $(t^{\prime},m,s-1)$ net in base $b$ with $t^{\prime}<t$.

For scrambled net sampling, reducing the dimension reduces an upper bound on the variance. For any function $f\in L^{2}[0,1)^{d}$, the variance using a scrambled $(t,m,d)$-net in base $b$ is at most $b^{t+d}$ times the MC variance. Reducing the dimension reduces the bound to $b^{t+d-1}$ times the MC variance. For digital constructions in base $2$, there are sharper bounds on these ratios, $2^{t+d-1}$ and $2^{t+d-2}$, respectively [33] and for some nets described there, even lower bounds apply.

As remarked above, pre-integration over a variable $x_{j}$ that $f(\boldsymbol{x})$ uses will reduce the variance under scrambled net sampling. This reduction does not require $f$ to be monotone in $x_{j}$, though such cases have the potential to bring a greater improvement. Pre-integration can either increase or decrease the mean dimension because

and pre-integration could possibly reduce the denominator by a greater proportion than it reduces the numerator.

In order to choose $x_{j}$ to pre-integrate over, we can look at the variance reduction we get. Pre-integrating over $x_{j}$ reduces the scrambled net variance by

Evaluating this quantity for each $j\in 1{:}d$ might be more expensive than getting a good estimate of $\mu$. However we don’t need to find the best $j$. Any $j$ where $f$ depends on $x_{j}$ will bring some improvement. Below we develop a principled and computationally convenient choice by choosing the $j$ which is most important as measured by a Sobol’ index [41] from global sensitivity analysis [37].

A convenient proxy replacement for (7) is

where $\sigma^{2}_{u}$ is the ANOVA variance component for the set $u$. The equality above follows because the ANOVA can be defined, as Sobol’ [40] did, by collecting up the terms involving $x_{j}$ for $j\in u$ from the orthogonal decomposition. Sobol’ used Haar functions. The right hand side of (8) equals $\overline{\tau}^{2}_{j}/n$. It counts all the variance components in which variable $j$ participates. From the orthogonality properties of ANOVA effects it follows that

The Jansen estimator [19] is an estimate of the above integral that can be done by a $d+1$ dimensional MC or QMC or RQMC sampling algorithm. Our main interest in this Sobol’ index estimator is that we use it as a point of comparison to the use of active subspaces in choosing a projection of a Gaussian vector along which to pre-integrate.

The discussion in Section 3 motivates pre-integration of $f(\boldsymbol{x})$ for $\boldsymbol{x}\sim\mathcal{N}(0,I)$ over a linear combination $\theta^{\mathsf{T}}\boldsymbol{x}$ having the largest Sobol’ index over unit vectors $\theta$. For $\theta_{1}=\theta$, let $\Theta=(\theta_{1},\theta_{2},\dots,\theta_{d})\in^{d\times d}$ be an orthogonal matrix and write

Then we define $\overline{\tau}^{2}_{\theta}(f)$ to be $\overline{\tau}^{2}_{1}$ in the ANOVA of $f_{\Theta}(\boldsymbol{x})$. First we show that $\overline{\tau}^{2}_{\theta}$ does not depend on the last $d-1$ columns of $\Theta$.

Let $z$, $\tilde{z}$ and $\boldsymbol{y}$ be independent with distributions $\mathcal{N}(0,1)$, $\mathcal{N}(0,1)$ and $\mathcal{N}(0,I_{d-1})$ respectively. Let $\boldsymbol{x}=z\theta_{1}+\Theta_{-1}\boldsymbol{y}$ and $\tilde{\boldsymbol{x}}=\tilde{z}\theta_{1}+\Theta_{-1}\boldsymbol{y}$. Using the Jansen formula for $\theta=\theta_{1}$,

Now for an orthogonal matrix $Q\in^{(d-1)\times(d-1)}$, let

where $\tilde{\theta}_{1}=\theta_{1}$. In this parameterization we get