On Stein's method for multivariate normal approximation

Introduction

for all $f$ for which the left-hand side exists and is finite. The operator $T_{o}$ defined on $C^{1}$ functions by

is called the characterizing operator of the standard normal distribution. The left-inverse to $T_{o}$, denoted $U_{o}$, is defined by the equation

where $Z$ is a standard normal random variable; the boundedness properties of $U_{o}$ are an essential ingredient of Stein’s method.

with the random quantities $E_{1},E_{2},E_{3}$ being small compared to $\lambda$, then $W$ is indeed approximately Gaussian, and its distance to Gaussian (in some metric) can be bounded in terms of the $E_{i}$ and $\lambda$.

The addition of a random error to this equation was not needed in the applications in , but is a straightforward modification of the theorems proved there.

After the initial draft of appeared on the ArXiv, a preprint was posted by Reinert and Röllin which generalized one of the abstract normal approximation theorems of . Instead of condition (4) above, they required

where $\Lambda$ is a positive definite matrix and $E$ is a random error. This more general condition allowed them to estimate the distance to Gaussian random vectors with non-identity (even singular) covariance matrices. They then introduced an insightful new method, “the embedding method” for approximating real random variables by the normal distribution, by observing that in many cases in which the condition (1) does not hold, the random variable in question can be viewed as one component of a random vector which satisfies condition (5) with a non-diagonal $\Lambda$. Many examples are given, both of the embedding method and the multivariate normal approximation theorem directly, including applications to runs on the line, statistics of Bernoulli random graphs, U-statistics, and doubly-indexed permutation statistics.

The approach taken in the published version of is to give smoothness conditions instead by requiring bounds on the quantities

The Wasserstein distance $d_{W}(X,Y)$ between the random variables $X$ and $Y$ is defined by

where $M_{1}(g)=\sup_{x\neq y}\frac{|g(x)-g(y)|}{|x-y|}$ is the Lipschitz constant of $g$. On the space of probability distributions with finite absolute first moment, Wasserstein distance induces a stronger topology than the usual one described by weak convergence, but not as strong as the topology induced by the total variation distance. See for detailed discussion of the various notions of distance between probability distributions.

The $n\times n$ identity matrix is denoted $I_{n}$ and the $n\times n$ matrix of all zeros is denoted $0_{n}$.

that is, $M_{k}(f)$ is the Lipschitz constant of the $k-1$-st derivative of $f$.

This general definition of $M_{k}$ is a departure from what was done by Raič in ; there, smoothness conditions on functions are also given in coordinate-independent ways, and $M_{1}$ and $M_{2}$ are defined as they are here, but in case $k=3$, the quantity $M_{3}$ is defined as the Lipschitz constant of the Hessian with respect to the Hilbert-Schmidt norm as opposed to the operator norm.

Abstract Approximation Theorems

This section contains the basic lemmas giving the Stein characterization of the multivariate Gaussian distribution and bounds to the solution of the Stein equation, together with two multivariate abstract normal approximation theorems and their proofs. The first theorem is a reworking of the theorem of Reinert and Röllin on multivariate normal approximation with the method of exchangeable pairs for vectors with non-identity covariance. The second is an analogous result in the context of “continuous symmetries” of the underlying random variable, as has been previously studied by the author in , , and (jointly with S. Chatterjee) in .

is a solution to the differential equation

Part (1) follows from integration by parts.

and so $\mathcal{L}(Y)=\mathcal{L}(Z)$ since $C^{\infty}$ is dense in the class of bounded continuous functions, with respect to the supremum norm.

For part (3), first note that since $g$ is Lipschitz, if $t\in\left(0,1\right)$

which is integrable on $\left(0,1\right)$, so the integral exists by the dominated convergence theorem.

To show that $U_{o}g$ is indeed a solution to the differential equation (9), let

The next lemma gives useful bounds on $U_{o}g$ and its derivatives in terms of $g$ and its derivatives. As in , bounds are most naturally given in terms of the quantities $M_{i}(g)$ defined in the introduction.

If, in addition, $\Sigma$ is positive definite, then

Remark: Bounds (3), (4), and (5) are mainly of use when $\Sigma$ has a fairly simple form, since they require an estimate for $\|\Sigma^{-1/2}\|_{op}$. They are also of theoretical interest, since they show that if $\Sigma$ is non-singular, then the operator $U_{o}$ is smoothing; functions $U_{o}g$ are typically one order smoother than $g$. The bounds (1) and (2), while not showing the smoothing behavior of $U_{o}$, are useful when $\Sigma$ is complicated (or singular) and an estimate of $\|\Sigma^{-1/2}\|_{op}$ is infeasible or impossible.

Write $h(x)=U_{o}g(x)$ and $Z_{x,t}=\sqrt{t}x+\sqrt{1-t}\Sigma^{1/2}Z$. Note that by the formula for $U_{o}g$,

for unit vectors $u_{1},\ldots,u_{k}$, and part (1) follows immediately.

For the second part, note that (10) implies that

For part (3), note that it follows by integration by parts on the Gaussian expectation that

For part (4), again using integration by parts on the Gaussian expectation,

For notational convenience, let $f=U_{o}g$. By the exchangeability of $(X,X^{\prime})$,

where $R$ is the error in the Taylor approximation. By conditions (1) and (2), it follows that

that is (making use of the definition of $f$),

where the first line is by the Cauchy-Schwarz inequality, the second is by the standard bound $\|AB\|_{H.S.}\leq\|A\|_{op}\|B\|_{H.S.},$ and the third uses the bounds (2) and (4) from Lemma 2.

Finally, by Taylor’s theorem and Lemma 2,

The first bound of the theorem results from choosing the first term from each minimum; the second bound results from the second terms.

For notational convenience, let $f=U_{o}g$. Beginning as before,

where $R$ is the error in the Taylor approximation.

Now, by Taylor’s theorem, there exists a real number $K$ depending on $f$, such that

Breaking up the expectation over the sets on which $|X_{\epsilon}-X|^{2}$ is larger and smaller than a fixed $\rho>0$,

The second term tends to zero as $\epsilon\to 0$ by condition 3; condition 2 implies that the first is bounded by $CK\|\Lambda^{-1}\|_{op}\rho$ for a constant $C$ depending on the distribution of $X$. It follows that

is stronger than condition (3) of Theorem 4 and may be used instead; this is what is done in the application given in Section 3.

In , singular covariance matrices are treated by comparing to a nearby non-singular covariance matrix rather than directly. However, this is not necessary as all the proofs except those explicitly involving $\Sigma^{-1/2}$ go through for non-negative definite $\Sigma$.

Examples

The following example was treated by Reinert and Röllin as an example of the embedding method. It should be emphasized that showing that the number of $d$-runs on the line is asymptotically Gaussian seems infeasible with Stein’s original method of exchangeable pairs because of the failure of condition (1) from the introduction, but in , the random variable of interest is embedded in a random vector whose components can be shown to be jointly Gaussian by making use of the more general condition (5) of the introduction. The example is reworked here making use of the analysis of together with Theorem 3, yielding an improved rate of convergence.

assuming the torus convention, namely that $X_{n+k}=X_{k}$ for any $k$. For this example, we assume that $d<\frac{n}{2}$. To make an exchangeable pair, $d-1$ sequential elements of $X:=(X_{1},\ldots,X_{n})$ are resampled. That is, let $I$ be a uniformly distributed element of $\{1,\ldots,n\}$ and let $X_{1}^{\prime},\ldots,X_{n}^{\prime}$ be independent copies of the $X_{i}$. Let $X^{\prime}$ be constructed from $X$ by replacing $X_{I},\ldots,X_{I+d-2}$ with $X_{I}^{\prime},\ldots,X_{I+d-2}^{\prime}$. Then $(X,X^{\prime})$ is an exhangeable pair, and, defining $V_{i}^{\prime}:=V_{i}(X)$ for $i\geq 1$, it is easy to see that

where sums $\sum_{a}^{b}$ are taken to be zero if $a>b$. It follows that

Standard calculations show that, for $1\leq j\leq i\leq d,$

It then follows from (21) that, for $1\leq i,j\leq d$,

Condition (1) of Theorem 3 thus applies with $E=0$ and $\Lambda$ as above.

To apply Theorem 3, an estimate on $\|\Lambda^{-1}\|_{op}$ is needed. Following Reinert and Röllin, we make use of known estimates of condition numbers for triangular matrices (see, e.g., the survey of Higham ). First, write $\Lambda=:\Lambda_{E}\Lambda_{D}$, where $\Lambda_{D}$ is diagonal with the same diagonal entries as $\Lambda$ and $\Lambda_{E}$ is lower triangular with diagonal entries equal to one and $(\Lambda_{E})_{ij}=\frac{\Lambda_{ij}}{\Lambda_{jj}}$ for $i>j$. Note that all non-diagonal entries of $\Lambda_{E}$ are bounded in absolute value by $\frac{2\sqrt{p}}{d-1}$. From Lemeire , this implies the bounds

From Higham, $\|\Lambda_{E}^{-1}\|_{op}\leq\sqrt{\|\Lambda_{E}^{-1}\|_{1}\|\Lambda_{E}^{-1}\|_{\infty}},$ thus

Trivially, $\|\Lambda_{D}^{-1}\|_{op}=\frac{n}{d-1}$, and thus

It was determined by Reinert and Röllin that

Using these bounds in inequality (13) from Theorem 3 yields the following.

where $Z$ is a standard $d$-dimensional Gaussian random vector.

Remarks: Compare this result to that obtained in :

where $|h|_{2}=\sup_{i,j}\left\|\frac{\partial^{2}h}{\partial x_{i}\partial x_{j}}\right\|_{\infty}$ and $|h|_{3}=\sup_{i,j,k}\left\|\frac{\partial^{3}h}{\partial x_{i}\partial x_{j}\partial x_{k}}\right\|_{\infty}$.

2. Eigenfunctions of the Laplacian

Consider a compact Riemannian manifold $M$ with metric $g$. Integration with respect to the normalized volume measure is denoted $d\overline{\rm vol}$, thus $\int_{M}1d\overline{\rm vol}=1.$ For coordinates $\left\{\frac{\partial}{\partial x_{i}}\right\}_{i=1}^{n}$ on $M$, define

where the coefficient implicit in the $O(\epsilon^{3})$ depends on $f_{i}$ and $\gamma$ and $d_{x}f_{i}$ denotes the differential of $f_{i}$ at $x$. Recall that $d_{x}f_{i}(v)=\left\langle\nabla f_{i}(x),v\right\rangle$ for $v\in T_{x}M$ and the gradient $\nabla f_{i}(x)$ defined as above. Now, for $X$ fixed, $V$ is distributed according to normalized Lebesgue measure on $S_{X}M$ and $d_{X}f_{i}$ is a linear functional on $T_{X}M$. It follows that

exists and is finite; we will take $s(\epsilon)=\epsilon^{2}.$ Indeed, it is well-known (see, e.g., Theorem 11.12 of ) that

For the second condition of Theorem 4, it is necessary to determine

Choose coordinates $\left\{\frac{\partial}{\partial x_{i}}\right\}_{i=1}^{n}$ in a neighborhood of $X$ which are orthonormal at $X$. Then

(As before, the convergence requirement is satisfied since the $f_{i}$ are smooth and $M$ is compact.)

(where the implicit constants depend on the $f_{i}$ and on $k$), thus condition (3) of Theorem 4 is satisfied.

All together, we have proved the following.

Let $M$ be a compact Riemannian manifold and $f_{1},\ldots,f_{k}$ an orthonormal (in $L_{2}(M)$) sequence of eigenfunctions of the Laplacian on $M$, with corresponding eigenvalues $-\mu_{i}$. Let $X$ be a uniformly distributed random point of $M$. Then if $W:=(f_{1}(X),\ldots,f_{k}(X))$,

In this example, Theorem 6 is applied to the value distributions of eigenfunctions on flat tori. The class of functions considered here are random functions; that is, they are linear combinations of eigenfunctions with random coefficients.

Eigenfunctions of $\Delta_{B}$ are given by the real and imaginary parts of functions of the form

Consider a collection of $k$ random eigenfunctions $\{f_{j}\}_{j=1}^{k}$ of $\Delta_{B}$ on the torus which are linear combinations of eigenfunctions with random coefficients:

and, applying Theorem 6, we have shown that

Remarks: Note that if the elements of $\cup_{r=1}^{k}\mathcal{V}_{r}$ are mutually orthogonal, then the right-hand side becomes

Acknowledgements. The author thanks M. Meckes for many useful discussions. This research was supported by an American Institute of Mathematics five-year fellowship.

References