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] \E[B] \E[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 all variables are unbiased. 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] & = \Cov(A, B) \Cov(C, D) \\ & + \Cov(A, C) \Cov(B, D) \\ & + \Cov(A, D) \Cov(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 - 1)!!}{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_{\pi(g(1)), \, \pi(g(2))} \times \dots \times \Sigma_{\pi(g(k - 1)), \, \pi(g(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 now a combinatorics problem. For each sequence $(\pi(1), \dots, \pi(k))$ involving $\ell$ distinct incides, let's consider the partition $p(1) + \dots + p(\ell) = k$ that describes the number of times $\pi$ visits each index. We'll write $$M_p = (\E X_i^{p(1)}) \dots(\E X_i^{p(\ell)})$$ for the contribution to the sum $\sum_\pi \E\bigl[ X_{\pi(1)} \dots X_{\pi(k)}\bigr]$ made by this sequence $\pi.$ Let's also classify sequences by the ordered sets $I$ describing which of the $N$ indices they involve, and the order they introduce them.

For a given partition $p(1) + \dots + p(\ell)$ and an ordered set $I$ of size $\ell,$ there are a certain number of sequences $\pi$ associated with both. However, this number only depends on the partition $p.$ (Specifically, this is the number of groupings of a $k$-element set into subsets with sizes agreeing with the partition $p.$) Meanwhile, the number of ordered sets of size $\ell$ grows with $N$; this is $N (N - 1) \dots (N - \ell + 1) = (N)_\ell$. Overall, we get formulas like $$\begin{align*} \E Z^3 & = N M_{3} \\ \E Z^4 & = 3 M_{2, 2} (N)_2 + M_4 N \\ \E Z^5 & = 10 M_{2, 3} (N)_2 + M_5 N \\ \E Z^6 & = 15 M_{2, 2, 2} (N)_3 + \bigl(15 M_{2, 4} + 10 M_{3, 3} \bigr) (N)_2 + M_6 N. \end{align*}$$ If we divide $Z$ by $\sqrt N,$ the formulas for even moments converge to quantities like $3 M_{2, 2}$ and $15 M_{2, 2, 2},$ and formulas for odd moments converge to zero at rate $O(1/\sqrt N).$

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 large parts of this contribution turn out to depend only on some restricted properties of the parts?