Gaussian Moments

Any moment of a Gaussian distribution can be computed from its mean and covariance. For instance, $$\begin{align*} \E[A B C] & = \E[A] \Cov(B, C) \\ & + \E[B] \Cov(A, C) \\ & + \E[C] \Cov(A, B) \end{align*}$$ when $(A, B, C)$ are jointly Gaussian, which can be derived from the fact that $$\E[A B C] = 0$$ when $\E [A] = \E [B] = \E [C] = 0.$ However, the fourth moment may be a little less familiar. For unbiased and jointly Gaussian variables $(A, B, C, D),$ we can show that $$\begin{align*} \E[A B C D] & = \E[A B] \E[C D] \\ & + \E[A C] \E[B D] \\ & + \E[A D] \E[B C]. \end{align*}$$ Where does this expression come from?

The simplest way to compute moments is probably through the moment generating function $$m(\lambda) = \E_{X \sim p} [e^{\langle \lambda, X \rangle}].$$ For a centered Gaussian with covariance $\Sigma,$ a little algebra gets us $$m(\lambda) = \exp\left( \frac{1}{2} \lambda^T \Sigma \lambda \right).$$ In the case that $X$ is real-valued, so $\lambda$ and $\Sigma = \sigma^2$ are scalars, developing the series gives $$m(\lambda) = 1 + \frac{1}{2} (\sigma \lambda)^2 + \frac{3}{4!} (\sigma \lambda)^4 + \frac{5 \times 3}{6!} (\sigma \lambda)^6 + \dots.$$ Thus, the second moment is $\sigma^2,$ the fourth moment is $3 \sigma^4,$ the sixth moment is $15 \sigma^6,$ and so on. In the multivariable case, we can write $$m(\lambda) = \sum_r \frac{1}{2^r r!} (\lambda^T \Sigma \lambda)^{r} = \sum_r \frac{k!!}{k!} \langle \Sigma^{\otimes r}, \lambda^{\otimes k} \rangle$$ where $k = 2r$ and $\Sigma^{\otimes r}$ is a symmetric tensor power of the second moment. Explicitly, if we view $\Sigma^{\otimes r}$ as an ordinary tensor which happens to be symmetric, its value at an index $\pi = (\pi(1), \dots, \pi(k))$ is given by $$(\Sigma^{\otimes r})_\pi = \frac{1}{k!} \sum_{g} \Sigma_{g(\pi(1)), \, g(\pi(2))} \times \dots \times \Sigma_{g(\pi(k - 1)), \, g(\pi(k))}$$ where $g$ runs over all permutations of $k$ elements.

This general expression for the higher moments isn't too surprising; how could you get a symmetric tensor of order $k$ from the second moment besides by powering and symmetrizing? However, that explanation is a little vague, and the derivation with characteristic functions feels a little artificial. The following is a line of thought that helped clarified the story for me.

Let's think of a sum $$Z = \sum_{i = 1}^N X_i$$ of many unbiased i.i.d. samples $(X_1, X_2, \dots).$ Let's see if we can understand why the higher moments of $Z$ depend predominantly on the second moment of $X_i$ without using the CLT.

We'll work with ordinary polynomials first. Take four linear functions $a, b, c, d \colon \R^d \to \R$ of our vector-valued variables, and consider $$\E \bigl[a(Z) b(Z) c(Z) d(Z) \bigr] = \E \left[ \sum_i a(X_i) \sum_i b(X_i) \sum_i c(X_i) \sum_i d(X_i)\right].$$ This develops into expectation of products $a(X_{\pi(1)}) b(X_{\pi(2)}) c(X_{\pi(3)}) d(X_{\pi(4)})$ as $\pi$ runs over sequences of indices. When a sample $X_i$ appears alone in a product, the expectation of the whole product is zero. This leaves us with products involving either one or two distinct samples. Furthermore, the two-sample products occur more frequently: $$\begin{align*} \E \bigl[a(Z) b(Z) c(Z) d(Z) \bigr] & = N (N - 1) \E \bigl[ a(X_1) b(X_1) c(X_2) d(X_2) \bigr] \\ & + N(N - 1) \E \bigl[ a(X_1) b(X_2) c(X_1) d(X_2) \bigr] \\ & + N(N - 1) \E \bigl[ a(X_1) b(X_2) c(X_2) d(X_1) \bigr] \\ & + N \E \bigl[ a(X_1) b(X_1) c(X_1) d(X_1) \bigr]. \end{align*}$$ The first three terms can be evaluated in terms of the second moment of the population distribution, and correspond to our formula above for the fourth moment of a Gaussian. The last term, which depends on the fourth moment of the sample population, vanishes for large $N$ when we divide $Z$ by $\sqrt N.$

These calculations can be done more generally using tensor products. For convenience, let's write symmetric tensor products as scalar products. The $k$th moment of $Z$ expands as $$\E Z^k = \sum_\pi \E \bigl[ X_{\pi(1)} \dots X_{\pi(k)} \bigr],$$ where $\pi$ runs over all $k$-long sequences of indices. Since the variables $X_i$ are independent, each expectation in the sum above is a product of moments; for instance, $$\E \bigl[ X_1 X_1 \dots X_{k/2} X_{k/2} \bigr] = (\E X_i^2)^{k/2}.$$ Just as in our previous example, some products of moments contribute more substantially than others when $N$ is large. Specifically, a product will appear more times in our sum when it uses more distinct variables $X_i.$

Developing the exact expression for $\E Z^k$ is just a combinatorics problem. For each partition $p(1) + \dots + p(\ell) = k$ of $k,$ let's write $$M_p = \binom{k}{p(1), \dots, p(\ell)} (\E X_i^{p(1)}) \dots(\E X_i^{p(\ell)}).$$ This is the total contribution to the sum $\sum_\pi \E\bigl[ X_{\pi(1)} \dots X_{\pi(k)}\bigr]$ made by sequences $\pi$ which run over a given set of $\ell$ indices in a prescribed order and whose corresponding products have powers matching the partition $p(1) + \dots + p(\ell) = k.$ Where $(N)_\ell = N (N - 1) \dots (N - \ell + 1)$ counts the number of ordered sets of size $\ell,$ we conclude that $$\begin{align*} \E Z^3 & = N M_{3} \\ \E Z^4 & = (N)_2 M_{2, 2} + N M_4 \\ \E Z^5 & = (N)_2 M_{2, 3} + N M_5 \\ \E Z^6 & = (N)_3 M_{2, 2, 2} + (N)_2 \bigl( M_{2, 4} + M_{3, 3} \bigr) + N M_6 \end{align*}$$ and so on. Dividing $Z$ by $\sqrt N$ gives, as you like, either the Isserlis theorem or the central limit theorem.

In some sense, I suppose this is an algebraic form of CLT; that is, we have explained why the expectation of a polynomial in $Z / \sqrt{N}$ converges to an expression that only involves the second moment of $X_i.$ I wonder, are there other situations where a property of system with many "parts" can be seen as the contribution of many "pathways", and an overwhelming part of this contribution depends only on some restricted properties of the parts?