Efficient simulation of high dimensional Gaussian vectors

TL;DR

This paper introduces a linear-storage Markov chain Monte Carlo method for efficiently simulating high-dimensional Gaussian vectors, significantly reducing computational costs compared to traditional Cholesky-based methods.

Contribution

The authors develop a novel MCMC algorithm that approximates Gaussian vectors with linear storage and provide theoretical bounds on its accuracy and efficiency.

Findings

Linear storage cost in dimension d

Faster than standard methods by nearly a factor of d

Provides bounds on Wasserstein distance

Abstract

We describe a Markov chain Monte Carlo method to approximately simulate a centered d-dimensional Gaussian vector X with given covariance matrix. The standard Monte Carlo method is based on the Cholesky decomposition, which takes cubic time and has quadratic storage cost in d. In contrast, the storage cost of our algorithm is linear in d. We give a bound on the quadractic Wasserstein distance between the distribution of our sample and the target distribution. Our method can be used to estimate the expectation of h(X), where h is a real-valued function of d variables. Under certain conditions, we show that the mean square error of our method is inversely proportional to its running time. We also prove that, under suitable conditions, our method is faster than the standard Monte Carlo method by a factor nearly proportional to d. A numerical example is given.

Tables

Table 1. Table 1: Running times of the MC method. The second column gives the time to perform the Cholesky decomposition. The third column gives the time to simulate Z 1 , … , Z n ′ subscript 𝑍 1 … subscript 𝑍 superscript 𝑛 ′ {Z}_{1},\dots,{Z}_{n^{\prime}} , and calculating h ​ ( A ​ Z i ) ℎ 𝐴 subscript 𝑍 𝑖 h({A}{Z}_{i}) , 1 ≤ i ≤ n ′ 1 𝑖 superscript 𝑛 ′ 1\leq i\leq n^{\prime} , where n ′ = 10 4 superscript 𝑛 ′ superscript 10 4 n^{\prime}=10^{4} , and does not incorporate the Cholesky decomposition running time.

$d$	Cholesky	Simulations
	decomposition
10	0.000	0.007
${10}^2$	0.001	0.19
${10}^3$	2.3	17.3
${10}^4$	2306	1702

Table 2. Table 2: The MCMC and MC methods for estimating E ​ ( max 1 ≤ i ≤ d ⁡ Y ​ ( s i ) ) , E subscript 1 𝑖 𝑑 𝑌 subscript 𝑠 𝑖 {\rm E}(\max_{1\leq i\leq d}{Y}({s_{i}})), with n = 100 ​ d 𝑛 100 𝑑 n=100d and burnin = n / 2 burnin 𝑛 2 \text{burnin}=n/2 . The MCMC RMSE is an estimate of MSE ​ ( n ; n / 2 ) MSE 𝑛 𝑛 2 \sqrt{\text{MSE}(n;n/2)} , which is calculated as explained in Section  F . The standard deviation Σ Σ \Sigma is estimated using the MC method with n ′ = 10 4 superscript 𝑛 ′ superscript 10 4 n^{\prime}=10^{4} . The last column gives τ MC ​ ( ϵ ) / τ MCMC ​ ( ϵ , 1 / 2 ) subscript 𝜏 MC italic-ϵ subscript 𝜏 MCMC italic-ϵ 1 2 \tau_{\text{MC}}(\epsilon)/\tau_{\text{MCMC}}(\epsilon,1/2) , where ϵ = MSE ​ ( n ; n / 2 ) italic-ϵ MSE 𝑛 𝑛 2 \epsilon=\sqrt{\text{MSE}(n;n/2)} , and τ MC ​ ( ϵ ) subscript 𝜏 MC italic-ϵ \tau_{\text{MC}}(\epsilon) and τ MCMC ​ ( ϵ , 1 / 2 ) subscript 𝜏 MCMC italic-ϵ 1 2 \tau_{\text{MCMC}}(\epsilon,1/2) are calculated in seconds.

$d$	MCMC	MCMC	MCMC	$\Sigma$	${\tau }_{\text{MC}}/{\tau }_{\text{MCMC}}$
	Average	RMSE	comp. time
${10}^2$	2.38	0.119	$0.002$	2.7	4
${10}^3$	3.02	0.069	$0.19$	2.7	26
${10}^4$	3.57	0.049	$16$	2.7	173

Equations

$$A{A}^T=V.$$

$$X_{n+1}=X_n+(g_n-e_{i_n}^T{X}_n)(V{e}_{i_n}),$$

$${\mathcal{W}}_2(\mu ,{\mu }^{\prime })=(\underset{Y\sim \mu ,Y^{\prime }\sim {\mu }^{\prime }}{\inf }\mathrm{E}(||Y-Y^{\prime }|{|}^2){)}^{1/2}.$$

$$\text{MSE}(n)=\mathrm{E}((\frac{{\sum }_{j=0}^{n-1}h(X_j)}{n}-m{)}^2)$$

$$X^{\prime }=X+(g-e_1^{T}X)(V{e}_1),$$

$$X_{n+1}^{\text{HR}}=X_n^{\text{HR}}+g_n^{\text{HR}}{u}_n,$$

$$||A||=||A|{|}_V=\sqrt{\mathrm{tr}(A^{T}VA)}$$

$$Y_{n+1}=P_{i_n}{Y}_n+g_n{f}_{i_n}.$$

$$X_{n+1}=(I-V{e}_{i_n}{e}_{i_n}^T)X_n+g_n(V{e}_{i_n}),$$

$$Z_n=Y_n+M_n{Z}_0.$$

$${\lambda }_{\max }(V-\mathrm{Cov}(X_n))\le ||M_n|{|}^2.$$

$$\mathrm{E}(Z_n{Z}_m^T)=M_{m,n}.$$

$$\mathrm{Cov}(X_n)=V-\sqrt{V}{M}_n{M}_n^T\sqrt{V},$$

$$\mathrm{tr}(V-\mathrm{Cov}(X_n))=\mathrm{E}(||X_n-\sqrt{V}{Z}_n|{|}^2)=||M_n|{|}^2.$$

$$\mathrm{tr}(V-\mathrm{Cov}(X_n))=\mathrm{E}(||M_n|{|}^2).$$

$$\mathrm{E}(||X_n-\sqrt{V}{Z}_n|{|}^2|i_0,\dots ,i_n)=||M_n|{|}^2,$$

$$\mathrm{E}(||X_n-\sqrt{V}{Z}_n|{|}^2)=\mathrm{E}(||M_n|{|}^2).$$

$$\mathrm{E}(X_n{X}_n^T|i_0,\dots ,i_n)\le V.$$

$$\sum _{j=0}^n\mathrm{E}(||M_j|{|}^2)\le d^2,$$

$$\mathrm{tr}(V-\mathrm{Cov}(X_n))\le \frac{d^2}{n}.$$

$$\mathrm{E}(e_{i_j}{e}_{i_j}^T)$$

$$\mathrm{E}(f_{i_j}{f}_{i_j}^T)$$

$$\mathrm{E}(P_{i_j})=I-d^{-1}V.$$

$$\mathrm{E}(||v_{j+1}|{|}^2)$$

$$d^{-1}\mathrm{E}(v_j^{T}V{v}_j)=\mathrm{E}(||v_j|{|}^2)-\mathrm{E}(||v_{j+1}|{|}^2).$$

$$\sum _{j=0}^n\mathrm{E}(v_j^{T}V{v}_j)\le d.$$

$$v_j^{T}V{v}_j=v^T{M}_j^{T}V{M}_{j}v,$$

$$v^T(\mathrm{E}(\sum _{j=0}^n{M}_j^{T}V{M}_j))v\le d.$$

$$\mathrm{tr}(\mathrm{E}(\sum _{j=0}^n{M}_j^{T}V{M}_j))\le d^2,$$

$$\mathrm{E}((h(X)-h(X^{\prime }){)}^2)\le {\kappa }^2(\mathrm{E}(||X-X^{\prime }|{|}^2){)}^{\gamma }$$

Full text

Keywords: Cholesky factorisation, Gaussian vectors, Markov chains, Monte Carlo simulation

1 Introduction

Monte Carlo simulation of Gaussian vectors is commonly used in a variety of fields such as weather prediction [gel2004], finance [Hull12, Chap. 13], and machine learning [russo2014learning, russoBVR2016]. This paper considers the problem of efficiently sampling a $d$-dimensional Gaussian vector ${X}$ with a given mean and a given $d\times d$ covariance matrix ${V}$. Since any Gaussian random variable is an affine function of a standard Gaussian random variable, we assume throughout the paper that the components of ${X}$ are standard Gaussian random variables, and so the diagonal elements of ${V}$ are $1$. Then ${X}$ can be simulated [glasserman2004Monte, Subsection 2.3.3] as follows. Let ${Z}$ be a $d$-dimensional vector of independent standard Gaussian random variables, and let ${A}$ be a $d\times d$ matrix such that

[TABLE]

Then ${A}{Z}\sim N({0},{V})$, i.e. ${A}{Z}$ is a $d$-dimensional Gaussian vector with covariance matrix ${V}$.

Such a matrix ${A}$ can be computed in $O(d^{{3}})$ time and $O(d^{{2}})$ space using Cholesky factorization or one of its variants [golub2013matrix, Subsections 4.2.5 and 4.2.8]. Once ${A}$ is calculated, ${A}{Z}$ can be computed in $O(d^{2})$ time. But, in several applications (see e.g. [gel2004]), $d$ is in the tens of thousands or more, and so the calculation of a Cholesky factorization on a standard computer may not be possible in practice, due to the high running time and/or storage cost. Alternative methods for generating Gaussian vectors have been developed in special cases. For instance, exact and efficient simulation of Gaussian processes on a regular grid in $\mathbb{R}^{q}$, $q\geq 1$, can be performed [WoodChan94, dietrich1997fast] using Fast Fourier transforms if the covariance matrix is stationary with respect to translations. Similar Fast Fourier Transform methods can be used for exact simulation of fractional Brownian surfaces on a regular mesh [stein2002fast]. Sparse Cholesky decomposition [rue2001fast] and iterative methods [aune2013iterative] have been proposed to generate efficiently Gaussian vectors when the precision matrix ${V}^{-1}$ is sparse.

This paper develops a new Markov Chain Monte Carlo method for approximate generation of a Gaussian vector ${X}$ with correlation matrix ${V}$. Our method is straightforward to implement and can be applied to any correlation matrix ${V}$ whose elements are known or easy to compute. It has a total storage cost of $O(d)$. At iteration $n$, it produces a $d$-dimensional vector ${X}_{n}$ whose distribution converges (according to the quadratic Wasserstein distance) to $N({0},{V})$ as $n$ goes to infinity. Assuming each element of ${V}$ can be computed in $O(1)$ time, the running time of each iteration is $O(d)$. Our method can for instance be used to approximately simulate spatial Gaussian processes of various types such as Matérn, powered exponential, and spherical on any subset of size $d$ of $\mathbb{R}^{2}$ with $O(d)$ storage cost (background on spatial statistics can be found in [diggle2003introduction]). While FFT methods can simulate such processes on regular grids, certain applications (e.g. [gel2004]) require the simulation of spatial Gaussian processes on non-regular subsets of $\mathbb{R}^{2}$.

We now describe our method in more detail. Let $(i_{n}$), $n\geq 0$, be a deterministic or a random sequence in $\{1,\ldots,d\}$, and let $(g_{n}$), $n\geq 0$, be a sequence of independent standard Gaussian random variables, independent of $(i_{n})$, $n\geq 0$. Define the Markov chain of $d$-dimensional column vectors ${X}_{n}$, $n\geq 0$, as follows. Let ${X}_{0}=0$ and, for $n\geq 0$, let

[TABLE]

where ${e}_{i}$ is the $d$-dimensional column vector whose $i$-th coordinate is $1$, and remaining coordinates are [math] (if $t\in\mathbb{R}$ and $u$ is a vector, $tu$ is the scalar product of $t$ and $u$). Since ${V}{e}_{i_{n}}$ is the $i_{n}$-th column of ${V}$, ${X}_{n+1}$ can be calculated from ${X}_{n}$ in $O(d)$ time, with storage cost $O(d)$. The motivation behind (1.2) is explained in Section 2, where we also show that (1.2) is a variant of the hit-and-run algorithm. A general description of the hit-and-run algorithm can be found in [smith1984].

Section 3 shows that, if $i_{n}$ are independent random variables uniformly distributed over $\{1,\dots,d\}$, then the quadratic Wasserstein distance between the distribution of ${X}_{n}$ and $N({0},{V})$ is at most $d/\sqrt{n}$. The quadratic Wasserstein distance between two probability distributions $\mu$ and $\mu^{\prime}$ over $\mathbb{R}^{d}$ is defined as

[TABLE]

To put this result into perspective, denote by $\mu_{\epsilon}$ the distribution of $N({0},(1-\epsilon){V})$, for $0\leq\epsilon\leq 1$. Then $\mathcal{W}_{2}(\mu_{\epsilon},\mu_{0})=(1-\sqrt{1-\epsilon})\sqrt{d}$, by [dowson1982frechet, Eq. 16]. Thus, after $n=O(d\epsilon^{-2})$ steps, which can be performed in $O(d^{2}\epsilon^{-2})$ total time, the quadratic Wasserstein distance between the distribution of ${X}_{n}$ and $\mu_{0}$ is at most $\mathcal{W}_{2}(\mu_{\epsilon},\mu_{0})$. Section 4 shows that, if $h$ is a real-valued function on $\mathbb{R}^{d}$ satisfying certain conditions, and $i_{0},\dots,i_{n-1}$ are independent random variables uniformly distributed over $\{1,\dots,d\}$, then $m={\rm E}(h({X}))$, where ${X}\sim N({0},{V})$, is well approximated by $n^{-1}\sum^{n-1}_{j=0}h({X}_{j})$. More precisely, Theorem 4.1 give explicit bounds on the mean square error

[TABLE]

of this estimator. For instance, if $h$ is $\kappa$-Lipschitz, Theorem 4.1 implies that $n\text{MSE}(n)\leq 18\kappa^{2}d^{2}$. We give an example with $n=\Theta(d)$ where this bound is tight, up to a constant. To our knowledge, for general $V$, no previous methods achieve a similar tradeoff between the running time and the Wasserstein distance, or between the running time and the mean square error, when $n=\Theta(d)$. Section 5 assumes that ${V}$ is positive definite and shows that, under suitable conditions, $\text{MSE}(n)\sim cn^{-1}$ as $n$ goes to infinity, where $c$ is a constant. It also gives an explicit geometric bound on the Wasserstein distance between the distribution of ${X}_{n}$ and $N({0},{V})$, and an explicit bound on the mean square error of a related estimator of $m$. Section 6 gives examples and a numerical simulation, and shows that, under certain conditions, the total time needed by our method to achieve a given standarized mean square error is ${O}^{*}(d^{2})$. Concluding remarks are given in a closing section.

An introduction to MCMC methods can be found in [dellaportas2003introduction]. Our proof-techniques are based on coupling arguments. Conductance techniques can also be used to analyse mixing properties of Markov chains (see e.g. [sinclair1992improved, kahale1997semidefinite, diaconis2009markov]). Chernoff bounds for reversible discrete Markov chains in terms of the spectral gap have been established in [kahale1997large, gillman1998chernoff]. Previous theoretical results on the performance of hit-and-run algorithms have focused on their mixing properties (see [VempalaCousins2016, BRS1993hit] and references therein). For instance, after appropriate preprocessing, the hit-and-run algorithm for sampling from a convex body [lovasz1999hit] produces an approximately uniformly distributed sample point after $O^{*}(d^{3})$ steps. For general log-concave functions, after appropriate preprocessing [lovaszVempala2006], the hit-and-run algorithm mixes in $O^{*}(d^{4})$ steps. Note that, while the algorithms in [lovasz1999hit, lovaszVempala2006] require a pre-processing phase to make the target distribution ”well-rounded”, our method does not. When ${V}$ is positive definite, the Metropolis and Gibbs algorithms, and an algorithm for sampling from general log-concave functions using a Langevin stochastic differential equation [Moulines2016sampling], could be used to approximately sample from $N(0,V)$, but these algorithms require the calculation of ${V}^{-1}$. Standard algorithms for inverting a matrix take $\Theta(d^{3})$ time and $\Theta(d^{2})$ space, however, and so the pre-processing cost of these algorithms is as high as the Cholesky decomposition cost. Omitted proofs are in the appendix.

2 Motivation, notation and general properties

We motivate (1.2) by assuming that ${V}$ is positive definite, which implies the existence of a lower-triangular matrix ${A}$ satisfying (1.1) and such that ${A}$ and ${V}$ have the same first column [glasserman2004Monte, Subsection 2.3.3]. Let ${Z}\sim N(0,{I})$, where ${I}$ is the $d\times d$ identity matrix, and $g$ a standard Gaussian random variable independent of ${Z}$. Set ${Z}^{\prime}={Z}+(g-{e}_{1}^{T}{Z}){e}_{1}$, and let ${X}={A}{Z}$ and ${X}^{\prime}={A}{Z}^{\prime}={X}+(g-{e}_{1}^{T}{Z})({A}{e}_{1}$). Note that ${Z}^{\prime}$ is obtained from ${Z}$ by replacing its first component ${e}_{1}^{T}{Z}$ with $g$, and so ${Z}^{\prime}\sim N(0,{I})$. Hence ${X}\sim{X}^{\prime}\sim N({0},{V})$. But, since ${A}{e}_{1}$ (resp. ${V}{e}_{1}$) is the first column of ${A}$ (resp. ${V}$), ${A}{e}_{1}={V}{e}_{1}$. Furthermore, ${e}_{1}^{T}{A}={e}_{1}^{T}$since the first line of ${A}$ is ${e}_{1}^{T}$, and so ${e}_{1}^{T}{X}={e}_{1}^{T}{Z}$. Thus

[TABLE]

and so the RHS of (2.1) is a centered Gaussian vector with covariance matrix ${V}$. (1.2) is obtained from (2.1) by replacing ${e}_{1}$, $g$, ${X}$ and ${X}^{\prime}$ with ${e}_{i_{n}}$, $g_{n}$, ${X}_{n}$ and ${X}_{n+1}$, respectively.

We now describe a generic standard hit-and-run algorithm to approximately sample from a real-valued density function $f$ on $\mathbb{R}^{d}$. First, choose ${X}^{\text{HR}}_{0}$ from a certain distribution. If we are currently at point ${X}^{\text{HR}}_{n}$, we first choose a random vector ${u}_{n}\in\mathbb{R}^{d}$ according to a certain distribution, and then set

[TABLE]

where $g_{n}^{\text{HR}}$ is a random variable whose density at $t\in\mathbb{R}$ is proportional to $f({X}^{\text{HR}}_{n}+t{u}_{n})$.

If ${V}$ is positive definite, the density $f({x})$ of $N(0,V)$ at ${x}\in\mathbb{R}^{d}$ is $\exp(-{x}^{T}{V}^{-1}{x}/2)$, up to a multiplicative constant. In the standard hit and run algorithm, ${u}_{n}$ is a uniformly distributed unit vector. However, if we set ${X}^{\text{HR}}_{0}=0$ and ${u}_{n}={V}{e}_{i_{n}}$, it follows after some calculations that $g^{\text{HR}}_{n}\sim N(-{e}_{i_{n}}^{T}{X}^{\text{HR}}_{n},1)$. Thus, we can choose $g^{\text{HR}}_{n}=g_{n}-{e}_{i_{n}}^{T}{X}^{\text{HR}}_{n}$, which implies by induction that ${X}_{n}^{\text{HR}}={X}_{n}$.

If ${x}$ is a $d$-dimensional vector, denote by $||{x}||$ the $l_{2}$-norm of ${x}$. For any $d\times d$ matrix ${A}$, the matrix ${A}^{T}{V}{A}$ is positive semi-definite. Let

[TABLE]

be the Frobenius norm of the matrix $\sqrt{{V}}{A}$. If ${A}$ and ${B}$ are symmetric $d\times d$ matrices, we say that ${A}\leq{B}$ if ${B}-{A}$ is positive semi-definite, and we denote by $\lambda_{\max}({A})$ the largest eigenvalue of ${A}$. If ${Z}$ is a centered $d$-dimensional random vector such that ${\rm E}(||{Z}||^{2})$ is finite, let ${\rm Cov}({Z})={\rm E}({Z}{Z}^{T})$ denote the covariance matrix of ${Z}$.

For $1\leq i\leq d$, let ${f}_{i}=\sqrt{{V}}{e}_{i}$ and ${P}_{i}={I}-{f}_{i}{f}_{i}^{T}$. Note that $||{f}_{i}||^{2}={e}_{i}^{T}{V}{e}_{i}=1$. Thus ${f}_{i}$ is a unit vector and ${P}_{i}$ is a projection matrix, i.e. ${P}_{i}^{2}={P}_{i}$. Define the random sequence of $d$-dimensional vectors ${Y}_{n}$, $n\geq 0$, as follows: ${Y}_{0}=0$ and

[TABLE]

By rewriting (1.2) as

[TABLE]

it can be shown by induction that ${X}_{n}=\sqrt{{V}}{Y}_{n}$.

For $0\leq m\leq n$, let ${M}_{m,n}={P}_{i_{n-1}}{P}_{i_{n-2}}\cdots{P}_{i_{m}}$, with ${M}_{n,n}={I}$, and let ${M}_{n}={M}_{0,n}$. Let ${Z}_{0}$ be a $d$-dimensional vector of independent standard Gaussian random variables which is independent of the sequence $(g_{n},i_{n})$, $n\geq 0$. For $n\geq 1$, let

[TABLE]

Since $\lambda_{\max}({A})\leq{\rm tr}({A})$ for a positive semi-definite matrix ${A}$, the following lemma implies that, if the sequence $(i_{k})$, $k\geq 0$, is deterministic, then ${X}_{n}$ is centered Gaussian and

[TABLE]

As a consequence, any entry of ${V}-{\rm Cov}({X}_{n})$ is upper-bounded, in absolute value, by $||{M}_{n}||^{2}$.

Lemma 2.1**.**

If the sequence $(i_{k})$, $0\leq k\leq n-1$, is deterministic, then, for $0\leq m\leq n$, ${X}_{n}$ and ${Y}_{n}$ are centered Gaussian vectors, ${Z}_{n}\sim N(0,{I})$, and

[TABLE]

Furthermore,

[TABLE]

${\rm Cov}({X}_{n})\leq{V}$, and

[TABLE]

Lemma 2.1 forms the basis for the proofs of our main results. Indeed, if ${M}_{m,n}$ goes to [math] as $n-m$ goes to infinity then, by (2.4), ${\rm Cov}({X}_{n})$ converges to ${V}$ as $n$ goes to infinity. Furthermore, if both $m$ and $n-m$ are sufficiently large then, by (2.3), ${Z}_{n}$ and ${Z}_{m}$ are nearly independent and, by (2.2), ${Y}_{m}$ (resp. ${Y}_{n}$) is close to ${Z}_{m}$ (resp. ${Z}_{n}$). Thus, ${Y}_{m}$ and ${Y}_{n}$ are nearly independent, as well, and so are ${X}_{m}$ and ${X}_{n}$. These arguments are informal since we have not defined the terms ”nearly independent” and ”close”, but give intuition behind the proofs of Theorems 3.1 and 4.1.

Lemma 2.2 below generalizes some results of Lemma 2.1 when the sequence $(i_{k})$, $k\geq 0$, is random.

Lemma 2.2**.**

If the sequence $(i_{k})$, $k\geq 0$, is deterministic or random, the quadratic Wasserstein distance between the distribution of ${X}_{n}$ and $N({0},{V})$ is at most $\sqrt{{\rm E}(||{M}_{n}||^{2})}$. Furthermore, ${Z}_{n}\sim N(0,{I})$, ${X}_{n}$ is centered, ${\rm Cov}({X}_{n})\leq{V}$, and

[TABLE]

Proof.

By Lemma 2.1, conditioning on $i_{0},\dots,i_{n}$, ${Z}_{n}\sim N({0},{I})$. Thus, the unconditional distribution of ${Z}_{n}$ is $N({0},{I})$, and $\sqrt{{V}}{Z}_{n}\sim N({0},{V})$. On the other hand, by (2.5),

[TABLE]

and so, by the tower law,

[TABLE]

By (1.3), it follows that the quadratic Wasserstein distance between the distribution of ${X}_{n}$ and $N(0,V)$ is at most $\sqrt{{\rm E}(||{M}_{n}||^{2})}$. On the other hand, it follows from Lemma 2.1 that ${\rm E}({X}_{n}|i_{0},\dots,i_{n})=0$. By the tower law, we infer that ${X}_{n}$ is centered. Similarly, by Lemma 2.1,

[TABLE]

Hence, by the tower law, ${\rm E}({X}_{n}{X}_{n}^{T})\leq{V}$, and so ${\rm Cov}({X}_{n})\leq{V}$. Once again, (2.6) follows from (2.5) by the tower law. ∎

3 Upper bound on the Wasserstein distance

We first show the following lemma.

Lemma 3.1**.**

If ${P}$ is a $d\times d$ projection matrix and ${A}$ is a $d\times d$ matrix, then $||{A}{P}||\leq||{A}||$.

Proof.

Let ${H}={A}^{T}{V}{A}$. Since ${\rm tr}(BC)={\rm tr}(CB)$, ${\rm tr}(PHP)={\rm tr}(HP)={\rm tr}(PH)$, and so ${\rm tr}(H)-{\rm tr}(PHP)={\rm tr}(({I}-{P})H({I}-{P}))$. Since $H$ is positive semi-definite, so is $({I}-{P})H({I}-{P})$, and so ${\rm tr}(PHP)\leq{\rm tr}(H)$. Equivalently, $||{A}{P}||^{2}\leq||{A}||^{2}$. This concludes the proof. ∎

Under the conditions stated in Theorem 3.1 below, by an argument similar to that surrounding Lemma 2.1, it follows from (3.2) that each entry of the matrix ${V}-{\rm Cov}({X}_{n})$ is at most $d^{2}/n$ in absolute value.

Theorem 3.1**.**

Assume that $i_{n}$, $n\geq 0$, are independent random variables uniformly distributed over $\{1,\dots,d\}$. For $n\geq 1$, the quadratic Wasserstein distance between the distribution of ${X}_{n}$ and $N({0},{V})$ is at most $d/\sqrt{n}$,

[TABLE]

and the sequence ${\rm E}(||{M}_{j}||^{2})$ is decreasing. Furthermore, for $n\geq 1$, ${X}_{n}$ is centered, ${\rm Cov}({X}_{n})\leq{V}$, and

[TABLE]

Proof.

For any non-negative integer $j$,

[TABLE]

Thus,

[TABLE]

and so

[TABLE]

Let ${v}$ be a unit $d$-dimensional vector. For $j\geq 0$, set ${v}_{j}={M}_{j}{v}$. Since ${P}_{i}$ is a projection and ${v}_{j+1}={P}_{i_{j}}{v}_{j}$ for $j\geq 0$, it follows that

[TABLE]

Hence

[TABLE]

As ${v}_{0}={v}$, we conclude that

[TABLE]

But

[TABLE]

and so, for any unit vector ${v}$,

[TABLE]

Thus, any diagonal entry of the matrix ${\rm E}(\sum^{n}_{j=0}{M}_{j}^{T}{V}{M}_{j})$ is at most $d$. Hence,

[TABLE]

which implies (3.1). On the other hand, since ${M}_{j+1}^{T}={M}_{j}^{T}{P}_{i_{j}}$, Lemma 3.1 shows that $||{M}_{j+1}^{T}||\leq||{M}_{j}^{T}||$. Hence the sequence ${\rm E}(||{M}_{j}^{T}||^{2})$ is decreasing. Since ${M}_{j}$ and ${M}^{T}_{j}$ have the same distribution, ${\rm E}(||{M}_{j}^{T}||^{2})={\rm E}(||{M}_{j}||^{2})$. Thus, the sequence ${\rm E}(||{M}_{j}||^{2})$ is decreasing as well and so, by (3.1), $nE(||{M}_{n}||^{2})\leq d^{2}$. We conclude the proof using Lemma 2.2. ∎

4 Bounding the mean square error

We now define the class of $(\kappa,\gamma,W)$-Lipschitz functions, with $\kappa>0$ and $\gamma\in(0,1]$.

Definition 4.1**.**

Let $W$ be a $d\times d$ positive semi-definite matrix. We say that a real-valued Borel function $h$ of $d$ variables is $(\kappa,\gamma,W)$-Lipschitz if

[TABLE]

for any centered Gaussian column vector $\begin{pmatrix}{X}\\ {X}^{\prime}\\ \end{pmatrix}$ with ${\rm Cov}({X})\leq W$ and ${\rm Cov}({X}^{\prime})\leq W$, where ${X}$ and ${X}^{\prime}$ are $d$-dimensional.

We say that a function $h$ is $(\kappa,W)$-Lipschitz if it is $(\kappa,1,W)$-Lipschitz. For instance, if $h$ a real-valued $\kappa$-Lipschitz function on $\mathbb{R}^{d}$, i.e. $|h({x})-h({x}^{\prime})|\leq\kappa||{x}-{x}^{\prime}||$ for ${x}$, ${x}^{\prime}$ in $\mathbb{R}^{d}$, then $h$ is $(\kappa,W)$-Lipschitz for any $d\times d$ positive semi-definite matrix $W$. The following lemma gives an example of a $(\kappa,W)$-Lipschitz function in $\mathbb{R}$ which is not $\kappa^{\prime}$-Lipschitz for any $\kappa^{\prime}>0$.

Lemma 4.1**.**

Let $f(z)=e^{z}$. Then $f$ is $(e^{\nu}\sqrt{4\nu+1},\nu)$-Lipschitz for $\nu\geq 0$.

Let $h$ be a $(\kappa,\gamma,{V})$-Lipschitz function. Set $m={\rm E}(h({X}))$ and $\Sigma^{2}={\rm Var}(h({X}))$, where ${X}\sim N({0},{V})$, and denote by $\hat{h}$ the real-valued function on $\mathbb{R}^{d}$ defined by $\hat{h}({x})=h(\sqrt{{V}}{x})-m$. Note that ${\rm E}(\hat{h}({Z}))=0$ if ${Z}\sim N({0},{I})$, since $\sqrt{{V}}{Z}\sim N({0},{V})$. In particular, by Lemma 2.2, ${\rm E}(\hat{h}({Z}_{j}))=0$ for $j\geq 0$. In other words, ${\rm E}(h(\sqrt{V}{Z}_{j}))=m$, and so the average of $h(\sqrt{V}{Z}_{j})$, $b\leq j\leq n-1$, where $b$ is a burn-in period, is an unbiased estimator of $m$. The variance of this estimator equals $(n-b)^{-2}{\rm E}((\sum^{n-1}_{j=b}\hat{h}({Z}_{j}))^{2})$, which we bound using Lemma 4.2 below. Choices for the parameters $b$ and $\delta$ will be given in the sequel.

Lemma 4.2**.**

Let $b$, $n$ and $\delta$ be integers, with $0\leq\delta\leq b<n$. Then

[TABLE]

For $j\geq 0$, let

[TABLE]

The second moments of $\beta_{j}$ and of $\beta$ can be bounded as follows.

Lemma 4.3**.**

Let $b$ and $n$ be integers, with $0\leq b<n$. Then

[TABLE]

and

[TABLE]

Proof.

Assume first that the sequence $i_{0},\dots,i_{n-1}$, is deterministic. For $0\leq j\leq n-1$,

[TABLE]

The second equation follows from the relations ${\rm Cov}(\sqrt{{V}}{Z}_{j})={V}$ and ${\rm Cov}({X}_{j})\leq{V}$, and the last equation follows from (2.5). Thus, for any random sequence $i_{0},\dots,i_{n-1}$,

[TABLE]

The first inequality in the lemma then follows by taking expectations and using the tower law. The second inequality follows from the first one and the Cauchy-Schwartz inequality. ∎

Combining Lemmas 4.2 and 4.3 yields the following.

Lemma 4.4**.**

Let $b$, $n$ and $\delta$ be integers, with $0\leq\delta\leq b<n$. If $i_{0},\dots,i_{n-1}$ are independent random variables uniformly distributed over $\{1,\dots,d\}$, then

[TABLE]

Proof.

Since ${M}_{j,l}\sim{M}_{l-j}\sim{M}_{l-j}^{T}$, for any fixed $j\geq 0$,

[TABLE]

and

[TABLE]

Hence, by Lemma 4.2,

[TABLE]

As

[TABLE]

it follows by Lemma 4.3 that

[TABLE]

Since $\hat{h}({Y}_{j})=h({X}_{j})-m$, this concludes the proof. ∎

By applying Lemma 4.4 with $b=\delta=0$ and using (3.1), we get the following upper bound on $\text{MSE}(n)$.

Theorem 4.1**.**

Let $h$ be a $(\kappa,{V})$-Lipschitz function on $\mathbb{R}^{d}$, with $m={\rm E}(h({X}))$, where ${X}\sim N({0},{V})$. If $i_{0},\dots,i_{n-1}$ are independent random variables uniformly distributed over $\{1,\dots,d\}$, then

[TABLE]

4.1 Tightness of MSE bound

We now give an example where the bound on the mean square error in Theorem 4.1 is optimal, up to a multiplicative constant. Let ${V}={I}$ and let $h({x})=||{x}||$ for ${x}\in\mathbb{R}^{d}$. Thus $h$ is a $1$-Lipschitz function on $\mathbb{R}^{d}$. By [forbes2011statistical, Sec. 11.3],

[TABLE]

which implies by induction that $m\geq\sqrt{d/2}$. Furthermore, it follows from (1.2) and by induction on $n$ that $X_{n}$ has at most $n$ non-zero components, and that the non-zero components of $X_{n}$ are independent standard Gaussian random variables. Thus ${\rm E}(||X_{n}||^{2})\leq n$, and so ${\rm E}(||X_{n}||)\leq\sqrt{n}$ for $n\geq 0$. Thus,

[TABLE]

Hence, for $n=d/4$,

[TABLE]

and so

[TABLE]

Thus, the LHS of (4.3) is within an absolute constant from its RHS.

5 The positive semi-definite case

Let $h$ be a $(\kappa,\gamma,{V})$-Lipschitz function. Define $m$, $\Sigma$ and $\hat{h}$ as in Section 4. This section assumes that $V$ is positive definite and that, $i_{n}$, $n\geq 0$, are independent random variables uniformly distributed over $\{1,\dots,d\}$. Denote by $\lambda$ the smallest eigenvalue of ${V}$, and set

[TABLE]

The following lemma, combined with Lemma 2.2, implies a geometric bound on the Wasserstein distance between the distribution of ${X}_{n}$ and $N({0},{V})$ if ${V}$ is positive semi-definite.

Lemma 5.1**.**

For $j\geq 0$, ${\rm E}(||{M}_{j}||^{2})\leq d^{2}(1-\lambda d^{-1})^{j}$.

Proof.

We use the same notation as in the proof of Theorem 3.1. Since the largest eigenvalue of ${V}$ is at most ${\rm tr}({V})=d$, it follows from (3.4) that

[TABLE]

On the other hand, (3.3) implies that

[TABLE]

and so ${\rm E}(||{v}_{j}||^{2})\leq(1-\lambda d^{-1})^{j}$. Hence, for any unit-vector ${v}$,

[TABLE]

Thus each diagonal element of ${\rm E}({M}_{j}^{T}{V}{M}_{j})$ is at most $d(1-\lambda d^{-1})^{j}$, and so ${\rm tr}({\rm E}({M}_{j}^{T}{V}{M}_{j}))\leq d^{2}(1-\lambda d^{-1})^{j}$. This concludes the proof. ∎

As noted before, $n^{-1}{\sum^{n-1}_{j=0}h(\sqrt{{V}}{Z}_{j})}$ is an unbiased estimator of $m$. The following lemma implies that, if $c>0$, the variance of this estimator is $\Theta(n^{-1})$ as $n$ goes to infinity.

Lemma 5.2**.**

As $n$ goes to infinity, $n^{-1}{\rm E}((\sum^{n-1}_{j=0}\hat{h}({Z}_{j}))^{2})$ converges to $c$, where

[TABLE]

Theorem 5.1 below implies that, if $c>0$, $\text{MSE}(n)\sim cn^{-1}$ as $n$ goes to infinity.

Theorem 5.1**.**

As $n$ goes to infinity, $n{\rm E}((\frac{\sum^{n-1}_{j=0}h({X}_{j})}{n}-m)^{2})$ converges to $c$ .

Proof.

Define $\beta_{j}$ and $\beta$ via 4.2, with $b=0$. Let $\theta=(1-\lambda d^{-1})^{\gamma}$. By Lemma 4.3,

[TABLE]

The second inequality follows from Jensen’s inequality, and the last one from Lemma 5.1. Hence, by the Cauchy-Schwartz inequality,

[TABLE]

The last inequality follows from the relation $\theta^{1/2}\leq 1-\lambda\gamma d^{-1}/2$ (which is a consequence of Taylor’s formula with Lagrange remainder). On the other hand,

[TABLE]

But, by the Cauchy-Schwartz inequality and Lemma 5.2, for sufficiently large $n$,

[TABLE]

Using Lemma 5.2 once again, it follows that $n^{-1}{\rm E}((\sum^{n-1}_{j=0}\hat{h}({Y}_{j}))^{2})$ converges to $c$ as $n$ goes to infinity. This concludes the proof. ∎

We now show that the estimator ${(2/n)}{\sum^{n-1}_{j=n/2}h({X}_{j})}$ of $m$ has an exponentially decreasing bias. We also give a bound on the mean square error of this estimator which, in certain cases, is smaller than the RHS of (4.3).

Theorem 5.2**.**

Set $\delta=\left\lceil 4(\lambda\gamma)^{-1}d\ln(\kappa d/\Sigma)\right\rceil$. For $d\geq 3$ and even $n>0$,

[TABLE]

Furthermore, if $n>2\delta$,

[TABLE]

Proof.

Set $b=n/2$ and define $\beta$ via (4.2). Using calculations similar to the proof of Theorem 5.1, it follows that

[TABLE]

Since $1+x\leq e^{x}$ for $x\in\mathbb{R}$, it follows that

[TABLE]

This implies (5.1) since the $\hat{h}({Z}_{j})$’s are centered.

We now prove (5.2). We first note that, by applying (4.1) with ${X}\sim N(0,V)$, ${X}^{\prime}\sim N(0,V)$, ${X}$ and ${X}^{\prime}$ independent, if follows after some calculations that

[TABLE]

and so

[TABLE]

On the other hand, by Lemma (5.1) and Jensen’s inequality,

[TABLE]

The second equation follows from the inequality $(1-\lambda d^{-1})^{\gamma}\leq 1-\lambda\gamma d^{-1}$. Thus,

[TABLE]

Thus, by applying Lemma 4.4, it follows that

[TABLE]

By (5.4), this implies (5.2). ∎

6 Examples

Let $h$ be a real-valued Borel function of $d$ variables that can be calculated at any point in $O(d)$ time. Assume that ${V}$ is positive definite, and that both $m={\rm E}(h({X}))$ and $\Sigma^{2}={\rm Var}(h({X}))$ exist and are finite, where ${X}\sim N({0},{V})$. Denote by MCMC the algorithm that generates ${X}_{0},\dots,{X}_{n-1}$ via (1.2), where the $i_{j}$’s are independent and identically distributed over $\{1,\dots,d\}$, and estimates $m$ via

[TABLE]

where $b$ is a burn-in period. The standard Monte Carlo algorithm, referred to later as MC, first calculates a lower-triangular matrix ${A}$ satisfying (1.1) in $\Theta(d^{3})$ time via the procedure described in [glasserman2004Monte, Subsection 2.3.3], then generates $n^{\prime}$ independent $d$-dimensional vectors of independent standard Gaussian random variables ${Z}_{1},\dots,{Z}_{n^{\prime}}$, and estimates $m$ by taking the average of $h({A}{Z}_{j})$, $1\leq j\leq n^{\prime}$. The variance of this estimator is $V_{\text{MC}}(n^{\prime})=\Sigma^{2}/n^{\prime}$.

6.1 Comparison of the MC and MCMC methods

The mean square error of the $h_{n,b}$ estimator of $m$ is defined as

[TABLE]

Given $\epsilon\in(0,\Sigma)$, $n^{\prime}=\Sigma^{2}/\epsilon^{2}$ samples of the MC algorithm are needed to ensure that $V_{\text{MC}}(n^{\prime})=\epsilon^{2}$ (ignoring rounding issues). Calculating the Cholesky decomposition and $h({A}{Z}_{j})$, $1\leq j\leq n^{\prime}$, takes

[TABLE]

time. On the other hand, for $\epsilon>0$, if $h$ is $(\kappa,\gamma,{V})$-Lipschitz and $\xi\in\{0,1/2\}$, denote by $\tau_{\text{MCMC}}(\epsilon,\xi)$ the running time of the MCMC algorithm needed to ensure that $\text{MSE}(n;b)\leq\epsilon^{2}$, using burn-in period $b=\xi n$. If $\gamma=1$, by Theorem 4.1, after $n=18\kappa^{2}d^{2}/\epsilon^{2}$ steps of the MCMC algorithm, $\text{MSE}(n)\leq\epsilon^{2}$. Thus,

[TABLE]

Thus, if there is a constant $\phi\geq 1$ independent of $d$ such that $\kappa^{2}d\leq\phi\Sigma^{2}$ (this is equivalent to saying that (5.3) is tight, up to a constant), then, for fixed $\epsilon/\Sigma<1$,

[TABLE]

Similarly, under the assumptions of Theorem 5.2,

[TABLE]

Hence, if there are positive constants $\phi$ and $\phi^{\prime}$ independent of $d$ such that $\kappa\leq d^{\phi}\Sigma$ and $\lambda\gamma\geq\phi^{\prime}$, then, for fixed $\epsilon/\Sigma<1$,

[TABLE]

Examples where (6.1) or (6.2) hold are given below.

6.2 A Basket option

Consider a set of $d$ stocks $S_{1},\dots,S_{d}$. For $t\geq 0$, denote by $S_{i}(t)$ the price of $S_{i}$ at time $t$. Assume that $S_{1}(0)=S_{2}(0)=\dots=S_{d}(0)=1$. A Basket call option with maturity $T$ and strike $K$ is a financial derivative that pays the amount $((S_{1}(T)+\cdots+S_{d}(T))/d-K)^{+}$ at time $T$. Under a standard pricing model [glasserman2004Monte, Subsection 3.2.3], the price of a basket option is ${\rm E}(h({U}))$, where ${U}$ is a centered Gaussian vector with covariance matrix ${V}$ given by $V_{ij}=\text{Correl}(\ln(S_{i}(T)),\ln(S_{j}(T)))$ for $1\leq i\leq j\leq d$, and, for $x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$,

[TABLE]

where $r$ is the risk-free rate, and $\sigma_{i}$ is the volatility of $S_{i}$. Assume that the $\sigma_{i}$’s are bounded by a constant independent of $d$. It follows from Lemma 6.1 below that $h$ is $(\kappa,{V})$-Lipschitz, where $\kappa=O(d^{-1/2})$ as $d$ goes to infinity.

Lemma 6.1**.**

Let $g(x_{1},\dots,x_{d})=\max(\sum^{d}_{i=1}w_{i}e^{\sigma_{i}x_{i}}-K,0)$, where $w_{i}\geq 0$ for $1\leq i\leq d$. Then $g$ is $(\kappa,{V})$-Lipschitz, where $\kappa=\sqrt{\sum^{d}_{i=1}w_{i}^{2}e^{2\sigma_{i}^{2}}(4\sigma_{i}^{2}+1)}$.

Thus, by Theorem 4.1, $n=O(d/\epsilon^{2})$ steps of the MCMC algorithm are sufficient to ensure that $\text{MSE}(n)\leq\epsilon^{2}$, and so $\tau_{\text{MCMC}}(\epsilon,0)=O(d^{2}/\epsilon^{2})$. On the other hand, if $\Sigma=\Theta(1)$ (which is the case [Hull12, Sec. 25.14] if the volatilities and correlations are lower-bounded by a constant and $K=0$, for instance), then $\tau_{\text{MC}}(\epsilon)=\Theta(d^{3}+d^{2}/\epsilon^{2})$. In practice, though, $d$ is quite small.

6.3 The multivariate normal function

Let $a=(a_{1},\dots,a_{d})\in\mathbb{R}^{d}$, and $\hat{a}=\min_{1\leq i\leq d}|a_{i}|$. Set $h({x})=1_{{x}\leq a}$ for ${x}=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$, where ${x}\leq a$ if and only if $x_{i}\leq a_{i}$ for $1\leq i\leq d$. The following lemma and the analysis in Subsection 6.1 show that, if there are positive constants $\phi$ and $\phi^{\prime}$ independent of $d$ such that $\hat{a}\geq d^{-\phi}$, $\Sigma\geq d^{-\phi}$, and $\lambda\geq\phi^{\prime}$, then (6.2) holds for fixed $\epsilon/\Sigma<1.$

Lemma 6.2**.**

The function $h$ is $(3(d/\hat{a})^{1/3},1/3,W)$-Lipschitz for any $d\times d$ positive semi-definite matrix $W$.

Proof.

Let $\nu>0$. It is easy to see that if ${||{X}-{Y}||}<\nu$ and, for $1\leq i\leq d$, $|X_{i}|\not\in[|a_{i}|,|a_{i}|+\nu]$ and $|Y_{i}|\not\in[|a_{i}|,|a_{i}|+\nu]$, then $h({X})=h({Y})$. Hence

[TABLE]

By Chebyshev’s inequality, $\Pr(\nu\leq||{X}-{Y}||)\leq\nu^{-2}E(||{X}-{Y}||^{2}).$ On the other hand, a simple calculation shows that, for $z>0$, the density of any centered Gaussian random variable at $z$ is at most $1/z$. Hence,

[TABLE]

and a similar relation holds for ${Y}_{i}$. Thus,

[TABLE]

Minimizing over $\nu$ implies that

[TABLE]

∎

6.4 The maximum function

For $x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$, let $h(x)=\max_{1\leq i\leq d}x_{i}$. Then $h$ is $1$-Lipschitz. Let $X\sim N(0,V)$, where $V$ is a correlation matrix. Standard calculations show that $\Pr(h(X)>z)\leq de^{-z^{2}/2}$ for $z>0$. Thus, it follows after some calculations that $\Pr(h(X)^{2}>z\ln(d))\leq e^{-z/4}$ for $z\geq 2$ and $d\geq 3$, and so ${\rm E}(h(X)^{2})\leq 6\ln d$. On the other hand, since

[TABLE]

where $X_{1}$ is the first coordinate of $X$,

[TABLE]

Since $E(h(X))\leq 3\sqrt{\ln d}$, we conclude that $\Sigma^{2}\geq d^{-15}/\sqrt{2\pi}$. Hence (6.2) holds for fixed $\epsilon/\Sigma<1$ if there is a positive constant $\phi^{\prime}$ independent of $d$ such that $\lambda\geq\phi^{\prime}$.

6.5 A numerical example

In [gel2004], the temperatures $Y({s}_{1}),\cdots,Y({s}_{d})$ at a given set of $d$ locations ${s}_{1},\cdots,{s}_{d}$ in $\mathbb{R}^{2}$ at a given future time are modelled as a Gaussian vector where ${\rm E}(Y({s}_{i}))$ is a known function of ${s}_{i}$, ${\rm Var}(Y({s}_{i}))=\varrho$ and, for two different locations ${s}_{i}$ and ${s}_{j}$,

[TABLE]

where $\varrho$, $\sigma$ and $r$ are positive constants, with $\sigma^{2}<\varrho$. By simulating the vector $(Y({s}_{1}),\cdots,Y({s}_{d}))$, we can estimate the expected maximum temperature at these $d$ locations.

For simplicity of presentation, we assume thereafter that $Y({s}_{i})$ is centered for $1\leq i\leq d$, and so ${X}_{i}=\varrho^{-1/2}Y({s}_{i})$ is a standard Gaussian random variable. The correlation matrix ${V}$ of the Gaussian vector ${X}=(X_{1},\dots,X_{d})^{T}$ is given by

[TABLE]

for $i\neq j$. Since the matrix $(\exp(-||{s}_{i}-{s}_{j}||/r))_{1\leq i,j\leq d}$ is positive semi-definite [cressie2015statistics, Section 2.5], $\lambda\geq 1-(\sigma^{2}/\varrho)$. We use the MC and MCMC algorithms to estimate ${\rm E}(\max_{1\leq i\leq d}{Y}({s_{i}})).$ Note that $\max_{1\leq i\leq d}{Y}({s_{i}})=h(X)$, where $X=(X_{1},\dots,X_{d})$, and $h(x)=\sqrt{\varrho}\max(x_{1},\dots,x_{d})$ for $x=(x_{1},\dots,x_{d})\in\mathbb{R}^{d}$. The analysis in Subsection 6.4 shows that (6.2) holds for fixed $\epsilon/\Sigma<1$.

Our numerical simulations assume that $r=10$, $\varrho=8$, $\sigma^{2}=7.44$, and that, for $1\leq i\leq d$, the first (resp. second) coordinate of ${s}_{i}$ equals $\lfloor i/d^{\prime}\rfloor/d^{\prime}$ (resp. ($i\mod d^{\prime})/d^{\prime})$, where $d^{\prime}=\lceil\sqrt{d}\rceil$. After scaling, the $r$, $\varrho$, and $\sigma^{2}$ parameters are close to those estimated in [gel2004]. Our experiments were performed on a desktop PC with an Intel Pentium 2.90 GHz processor and 4 Go of RAM, running Windows 7 Professional. The codes were written in the C++ programming language, and the compiler used was Microsoft Visual C++ 2013. Computing times are given in seconds. Table 1 gives computing times of the MC method. Running the Cholesky factorization for $d=10^{5}$ without external storage causes memory overflow. Extrapolating the results in Table 1 shows that the Cholesky factorization would take a few weeks for $d=10^{5}$ if enough internal memory were available. Table 2 compares the MC and MCMC methods for $d$ up to $10^{4}$. For the tested parameters, the MCMC method is more efficient than the MC method, and its efficiency increases with $d$. For $d=10^{5}$, $n=100d$ and $b=n/2$, the MCMC average is $4.01$, and is calculated in $8969$ seconds. Lemma 6.2 shows that the MCMC method can also be used to calculate the probability that the maximum temperature over a set of $d$ points exceeds a certain level.

6.6 Other examples

Spatial Gaussian processes of various types such as Matérn, powered exponential, and spherical, restricted to any subset of size $d$ of $\mathbb{R}^{2}$, are centered Gaussian vectors with a covariance matrix whose entries can be calculated in $O(1)$ time. For instance, the covariance matrix $V$ of a powered exponential process restricted to a subset $\{{s}_{1},\dots,{s}_{d}\}$ of $\mathbb{R}^{2}$ is given by [diggle2003introduction]

[TABLE]

where $r>0$ and $0<\theta\leq 2$. The techniques of our paper can be used to simulate the restriction of such processes to any finite subset of $\mathbb{R}^{2}$.

7 Conclusion

We have shown how to simulate a Markov chain ${X}_{n}$, $n\geq 0$, such that the Wasserstein distance between the distribution of ${X}_{n}$ and $N(0,V)$ is at most $d/\sqrt{n}$. It takes $O(d)$ time to generate each step of the chain. Whereas the standard Monte Carlo simulation method has $\Theta(d^{2})$ storage cost, the storage cost of our method is $\Theta(d)$. Furthermore, by running the chain $n$ steps, our method can estimate ${\rm E}(h({X}))$, where ${X}$ is a centered Gaussian vector with covariance matrix ${V}$ and $h$ is a real-valued function of $d$ variables. Under certain conditions, we give an explicit upper bound on the mean square error of our estimate, and show that it is inversely proportional to the running time. We also prove that, in certain cases, the total time needed by our method to obtain a given standarized mean square error is ${O}^{*}(d^{2})$ time, whereas the standard Monte Carlo method takes $\Theta(d^{3})$ time.

Appendix A Proof of Lemma 2.1

Recall first that if ${Z}$ and ${Z}^{\prime}$ are independent centered $d$-dimensional random vectors such that ${\rm E}(||{Z}||^{2})$ and ${\rm E}(||{Z}^{\prime}||^{2})$ are finite, then ${\rm E}(||{Z}||^{2})={\rm tr}({\rm Cov}({Z}))$, ${\rm Cov}({A}{Z})={A}{\rm Cov}({Z}){A}^{T}$, and ${\rm Cov}({Z}+{Z}^{\prime})={\rm Cov}({Z})+{\rm Cov}({Z}^{\prime})$.

It can be shown by induction that ${Y}_{n}$ is a linear combination of $g_{{0}},\dots,g_{{n-1}}$, and so ${Y}_{n}$ is a centered Gaussian vector. Hence ${X}_{n}$ is also centered and Gaussian. Furthermore, ${Z}_{n}$ is centered Gaussian since it is a linear combination of $g_{{0}},\dots,g_{{n-1}}$ and of ${Z}_{0}$. Thus, ${Z}_{n}$ and $g_{n}$ are independent. For $0\leq l\leq n$,

[TABLE]

Thus, since ${Z}_{l}$ and $g_{l}$ are independent,

[TABLE]

It follows by induction that ${\rm Cov}({Z}_{l})={I}$, and so ${Z}_{l}\sim N(0,{I})$. Thus, (2.3) holds when $m=n$. Furthermore, since $g_{l}$ and ${Z}_{m}$ are independent for $0\leq m\leq l\leq\ n-1$, it follows from (A.1) that ${\rm E}({Z}_{l+1}{Z}^{T}_{m})={P}_{i_{l}}{\rm E}({Z}_{l}{Z}^{T}_{m})$. It follows by induction on $l$ that ${\rm E}({Z}_{l}{Z}^{T}_{m})={M}_{m,l}$ for $0\leq m\leq l\leq\ n$, Hence (2.3).

Since ${Y}_{n}$ is a linear combination of $g_{{0}},\dots,g_{{n-1}}$, the vectors ${Z}_{0}$ and ${Y}_{n}$ are independent. Thus, as ${\rm Cov}({Z}_{n})={I}$ and ${\rm Cov}({M}_{n}{Z}_{0})={M}_{n}{M}_{n}^{T}$, it follows from (2.2) that

[TABLE]

Hence (2.4), which implies that ${\rm Cov}({X}_{n})\leq{V}$. On the other hand, by (2.2), $\sqrt{{V}}{Z}_{n}-{X}_{n}=\sqrt{{V}}{M}_{n}{Z}_{0}$ and so,

[TABLE]

Hence,

[TABLE]

The second equation follows from the relation ${\rm tr}({A}{B})={\rm tr}({B}{A})$. Using (2.4), this implies (2.5). This concludes the proof. ∎

Appendix B Proof of Lemma 4.1

Let $(X,X^{\prime})$ be a Gaussian vector in $\mathbb{R}^{2}$, with $X\sim N({0},\nu)$, $X^{\prime}\sim N(0,\nu^{\prime})$ and $\nu^{\prime}\leq\nu$. We show the following, which immediately implies Lemma 4.1:

[TABLE]

Let

[TABLE]

Since ${\rm E}(e^{Z})=e^{\frac{1}{2}{\rm Var}({Z})}$ for any centered Gaussian random variable $Z$,

[TABLE]

The last equation follows from the inequalities $1-x\leq e^{-x}$ and $\cosh(x)\leq 1+x^{2}\cosh(x)/2$ (which is a consequence of Taylor’s formula with Lagrange remainder) for any real number $x$. On the other hand, $(\sqrt{\nu}-\sqrt{\nu^{\prime}})^{2}\leq\rho$ since ${\rm Cov}(X,X^{\prime})\leq\sqrt{\nu\nu^{\prime}}$, and so

[TABLE]

Hence, as $1\leq\cosh(x)$ for any real number $x$,

[TABLE]

This implies (B.1).

Appendix C Proof of Lemma 4.2

We first prove the following lemma which implies that, under certain conditions, the covariances between the components of ${Y}$ and ${Y}^{\prime}$ can be used to bound the covariance between $\hat{h}_{1}({Y})$ and $\hat{h}_{2}({Y}^{\prime})$.

Lemma C.1**.**

Let $\begin{pmatrix}{Y}\\ {Y}^{\prime}\\ \end{pmatrix}$ be a Gaussian column vector in $\mathbb{R}^{2d}$, with ${Y}\sim{Y}^{\prime}\sim N({0},{I})$ and ${\rm E}({Y}{Y}^{\prime T})={A}{B}^{T}$, where ${A}$ and ${B}$ are $d\times d$ matrices, with ${A}{A}^{T}\leq{I}$ and ${B}{B}^{T}\leq{I}$. Let $h_{1}$ and $h_{2}$ be two $(\kappa,\gamma,{V})$-Lipschitz functions on $\mathbb{R}^{d}$. Then

[TABLE]

Proof.

We first note that, if ${Z}\sim{Z}^{\prime}\sim N({0},{I})$ and ${Z}$, ${Z}^{\prime}$ are independent, then

[TABLE]

Furthermore, if $({Y},{Y}^{\prime})$ is a centered Gaussian vector in $\mathbb{R}^{2d}$ with ${\rm Cov}({Y})\leq{I}$ and ${\rm Cov}({Y}^{\prime})\leq{I}$, and $h$ is a $(\kappa,\gamma,{V})$-Lipschitz function on $\mathbb{R}^{d}$, then

[TABLE]

The second equation follows from the fact that $\sqrt{{V}}({Y}-{Y}^{\prime})$ is a centered Gaussian vector with covariance matrix $\sqrt{{V}}{\rm Cov}({Y}-{Y}^{\prime})\sqrt{{V}}$.

We now prove the lemma. Let $G$, ${G^{\prime}}$, ${G^{\prime\prime}}$, ${G}_{1}$ and ${G}_{2}$ be independent $d$-dimensional Gaussian vectors such that ${G}\sim{G^{\prime}}\sim{G^{\prime\prime}}\sim N({0},{I})$, ${G}_{1}\sim N({0},{I}-{A}{A}^{T})$, and ${G}_{2}\sim N({0},{I}-{B}{B}^{T})$. Note that ${G}_{1}$ and ${G}_{2}$ exist since ${I}-{A}{A}^{T}$ and ${I}-{B}{B}^{T}$ are positive semi-definite. Since ${{AG}}\sim N({0},{A}{A}^{T})$ and is independent of ${G}_{1}$, ${{AG}}+{G}_{1}\sim N({0},{I})$. Similarly, ${BG}+{G}_{2}\sim N({0},{I})$. Also, since ${G}$, ${G}_{1}$ and ${G}_{2}$ are independent and centered,

[TABLE]

Thus, the covariance matrix of the Gaussian vector $\begin{pmatrix}{{AG}}+{G}_{1}\\ {BG}+{G}_{2}\\ \end{pmatrix}$ is \left(\begin{array}[]{cc}{I}&{A}{B}^{T}\\ {B}{A}^{T}&{I}\\ \end{array}\right). Hence the centered Gaussian vectors $\begin{pmatrix}{{AG}}+{G}_{1}\\ {BG}+{G}_{2}\\ \end{pmatrix}$ and $\begin{pmatrix}{Y}\\ {Y}^{\prime}\\ \end{pmatrix}$ have the same covariance matrix, and so they have the same distribution. Thus,

[TABLE]

The second equation follows by applying (C.1) to each of the pairs $({{AG}}+{G}_{1},{BG}^{\prime\prime}+{G}_{2})$, $({{AG}}^{\prime}+{G}_{1},{BG}^{\prime\prime}+{G}_{2})$, and $({{AG}}^{\prime}+{G}_{1},{BG}+{G}_{2})$. The fourth equation follows from (C.2) and the relations ${\rm Cov}({{AG}}-{{AG}}^{\prime})=2{A}{A}^{T}$ and ${\rm Cov}({BG}-{BG}^{\prime\prime})=2{B}{B}^{T}$.

Hence

[TABLE]

Replacing $h_{2}$ with $-h_{2}$ concludes the proof.∎

We now prove Lemma 4.2. Assume first that the sequence $i_{0},\dots,i_{n-1}$, is deterministic. By Lemma 2.1, ${Z}_{j}\sim{Z}_{l}\sim N({0},{I})$ for $0\leq j\leq l\leq n$, and

[TABLE]

Let ${A}={M}_{j^{\prime},l}$ and ${B}={M}_{j,j^{\prime}}^{T}$, with $j^{\prime}=\lfloor(j+l)/2\rfloor$. Since ${A}^{T}$ is the product of $l-j^{\prime}$ projection matrices, it follows by induction on $l$ that $||{A}^{T}{x}||\leq||{x}||$ for $x\in\mathbb{R}^{d}$, and so ${A}{A}^{T}\leq{I}$. Similarly, ${B}{B}^{T}\leq{I}$. Since ${M}_{j,l}={A}{B}^{T}$, it follows from Lemma C.1 that

[TABLE]

Thus,

[TABLE]

The third equation follows from observing that $\hat{h}({Z}_{j})$ and $\hat{h}({Z}_{l})$ are centered and have standard deviation $\Sigma$. The last equation follows from the fact that each term ${M}_{j^{\prime},l}$ (resp. ${M}_{j,j^{\prime}}$) occurs at most twice in the last line.

Assume now that the sequence $i_{0},\dots,i_{n-1}$, is random. By (C.4),

[TABLE]

We conclude the proof by taking expectations and using the tower law. ∎

Appendix D Proof of Lemma 5.2

Let $a_{0}={\rm E}((\hat{h}({Z}_{0}))^{2})$ and, for $j>0$, let $a_{j}=2{\rm E}(\hat{h}({Z}_{0})\hat{h}({Z}_{j}))$. Since $({Z}_{j})$, $j\geq 0$, is a time-homogeneous Markov Chain and ${Z}_{j}\sim N(0,I)$, $2{\rm E}(\hat{h}({Z}_{k})\hat{h}({Z}_{k+j}))=a_{j}$ for $j>0$. Hence

[TABLE]

and so

[TABLE]

On the other hand, by applying (C.3) with $j=0$ and $l=j$, it follows that

[TABLE]

with $j^{\prime}=\lfloor j/2\rfloor$. The second equation follows from Jensen’s inequality. But, since ${M}_{j^{\prime},j}\sim{M}_{j-j^{\prime}}$,

[TABLE]

Furthermore, as $j^{\prime}\leq j-j^{\prime}$, it follows from Theorem 3.1 that

[TABLE]

Since ${M}_{j^{\prime}}\sim{M}_{j^{\prime}}^{T}$, we conclude that $|a_{j}|\leq 2\kappa^{2}({\rm E}(||{M}_{j^{\prime}}||^{2}))^{\gamma}.$ Hence, by Lemma 5.1, the series $\sum_{j=0}^{\infty}a_{j}$ is absolutely convergent. Thus, the RHS of (D.1) converges to $\sum^{\infty}_{j=0}a_{j}$ as $n$ goes to infinity. But ${\rm E}(\hat{h}({Z}_{0})\hat{h}({Z}_{j}))={\rm Cov}(h(\sqrt{{V}}{Z}_{0}),h(\sqrt{{V}}{Z}_{j}))$ since ${\rm E}(h(\sqrt{{V}}{Z}_{0}))={\rm E}(h(\sqrt{{V}}{Z}_{j}))=m$, and so $\sum^{\infty}_{j=0}a_{j}=c$. This concludes the proof.∎

Appendix E Proof of Lemma 6.1

Since the function $\max(z-K,0)$ is $1$-Lipschitz with respect to $z$, we can assume without loss of generality that $K=0$. Let $\begin{pmatrix}{U}\\ {U}^{\prime}\\ \end{pmatrix}$ be a centered Gaussian vector with ${\rm Cov}({U})\leq{V}$ and ${\rm Cov}({U}^{\prime})\leq{V}$, where ${U}$ and ${U}^{\prime}$ have dimension $d$. Let $U_{i}$ (resp. $U^{\prime}_{i}$) be the $i$-th component of ${U}$ (resp. ${U}^{\prime}$) and $\lambda_{i}=(4\sigma_{i}^{2}+1)e^{2\sigma_{i}^{2}}$. It follows from Lemma 4.1 that

[TABLE]

On the other hand, by the Cauchy-Schwartz inequality,

[TABLE]

Taking expectations, it follows that

[TABLE]

as desired. ∎

Appendix F Estimating the mean square error

This section describes a numerical method to estimate the mean square error of the MCMC method for moderate values of $d$. The method first calculates a matrix ${A}$ satisfying (1.1) using Cholesky factorization. An unbiased estimator of $m$ is then computed as follows. Define the Markov chain $({X}^{\prime}_{n})$, $n\geq 0$, by setting ${X}^{\prime}_{0}={A}{Z}_{0}$ and, for $n\geq 0$, let

[TABLE]

Thus ${X}^{\prime}_{n}$ satisfies the same recursion as ${X}_{n}$. By rewriting (F.1) as

[TABLE]

and using calculations similar to those in the proof of Lemma 2.1, it can be shown by induction that ${X}^{\prime}_{n}\sim N(0,V)$. Thus,

[TABLE]

is an unbiased estimator of $m$. Since $\text{MSE}(n;b)$ equals the variance of $h_{n,b}$ plus its square bias,

[TABLE]

In Subsection 6.5, $\text{MSE}(n;b)$ is estimated for $d\leq 10^{4}$ using (F.2), each term in the RHS of (F.2) being calculated via $100$ independent simulations of $h_{n,b}$ and $h^{\prime}_{n,b}$. For any $j\in\{0,\dots,n-1\}$, the same random variables $g_{j}$ and $i_{j}$ are used to calculate $h_{n,b}$ and $h^{\prime}_{n,b}$ via (1.2) and (F.1).

Appendix G Acknowledgments

This research has been presented at the Paris Bachelier Seminar, November 2016, and at the 2nd IMA Conference on the Mathematical Challenges of Big Data, London, December 2016. The author thanks Nicolas Chopin, Petros Dellaportas, Peter Glynn, Emmanuel Gobet, Benjamin Jourdain, and Didier Marteau for helpful conversations. This work was achieved through the Laboratory of Excellence on Financial Regulation (Labex ReFi) under the reference ANR-10-LABX-0095. It benefitted from a French government support managed by the National Research Agency (ANR).