# Gaussian fourth-moment formulas?

Suppose $$X$$ is an $$m\times n$$ random matrix where rows are I.I.D. samples of $$n$$-dimensional Gaussian, and $$A$$ is an $$n\times n$$ matrix. I'm looking for the value of the following quantities: $$E[X^T X A X^T X]$$

For zero-centered Gaussian with covariance $$\Sigma$$ and $$m=1$$, the following holds

$$E[X^T X A X^T X]=\\ \Sigma A \Sigma + \Sigma A^T \Sigma + \Sigma \operatorname{Tr}(\Sigma A)$$

What is the expression for more general case?

Edit The following formula from the answer has been confirmed to work:

$$$$\frac{1}{m^2}E[X^TXAX^TX] = \\\frac{(m-1)}{m}HAH + \frac{1}{m} \left( H(A+A^T)H - \mu\mu^T(A+A^T)\mu\mu^T + H\operatorname{Tr}(HA) \right),$$$$ where $$H=\Sigma+\mu\mu^T$$

• Maybe this section in Wikipedia's page for the Multivariate Normal can help.
– Yves
Sep 22 at 7:14
• Just diagonalize $A$ and apply the usual formulas for the second and fourth univariate moments.
– whuber
Sep 22 at 11:39
• I need a formula for symbolic $A$, analogous to $m=1$ case Sep 22 at 14:02

I have never seen an closed form expression for this. Probably because it is quite ugly. I have worked with a similar expression before, and I'd be happy to see if my expression is stands up to yours. So here is one way to compute the matrix expectation.

I will work with $$\frac{X^TX}{m}$$ rather than $$X^TX$$, since that is the empirical 2nd order moment matrix. It is a quantity that will converge.

Lemma First some tricks that need index manipulations for me to see. If $$y\sim \mathcal{N}(0,\Sigma)$$, $$\mu$$ a vector and $$A$$ is a matrix, $$E[y\mu{}^T A y \mu^T] = \Sigma{}A^T\mu\mu^T$$ $$E[\mu{}y^TA\mu{}y^T] = \mu\mu^T A^T \Sigma$$

Proof: Follow the indices! I use Einstein summation convention to be able declutter notation. $$E[y\mu{}^T A y \mu]_{ij} = E[y_i\mu_kA_{kl}y_l\mu_j] = E[y_iy_l]A^T_{lk} \mu_k \mu_j= \Sigma_{il}A^T_{lk}\mu_k\mu_j = [\Sigma{}A^T\mu\mu^T]_{ij}$$ The second equation is similarly proven. QED

Next, we will look at the zero-mean case as a warm-up exercise.

Lemma Assume $$m$$ independent gaussian random vectors $$x^{(i)} \overset{iid}{\sim} \mathcal{N}(0,\Sigma)$$. Row-stack them in the matrix $$X$$. Then $$$$\begin{split} E\left[\frac{X^TX}{m}A\frac{X^TX}{m}\right] &= \frac{m-1}{m} \Sigma{}A\Sigma{} \\ & + \frac{1}{m}\left( \Sigma{}(A+A^T)\Sigma{} + \operatorname{Tr}(\Sigma{}A)\Sigma{} \right) \end{split}$$$$

Proof:

$$\frac{1}{m^2}E[X^TXAX^TX] = \sum_{i,j}E[x^{(i)}{x^{(i)}}^TAx^{(j)}{x^{(j)}}^T]$$ There are $$m(m-1)$$ terms where $$i\neq{}j$$, and $$m$$ terms where they agree. Since $$x^{(i)}$$ and $$x^{(j)}$$ are iid, we can let $$x$$ denote either of them, and state $$\frac{1}{m^2}E[X^TXAX^TX] = \frac{1}{m^2}( m(m-1) E[xx^T]AE[xx^T] + m E[xx^TAxx^T])$$ The first term is just the variances. The second term needs isserlis theorem. It gives (again, einstein summation convention) \begin{align} E[xx^TAxx^T]_{ij} &= E[x_ix_kA_{kl}x_lx_j] \\&= A_{kl}(E[x_ix_k]E[x_lx_j]+E[x_ix_l]E[x_kx_j]+E[x_ix_j]E[x_lx_k] \\&= \left\{ \Sigma(A+A^T)\Sigma + \Sigma \operatorname{Tr}(A\Sigma)\right\}_{ij} \end{align} Armed with this, we can simply state $$\frac{1}{m^2}E[X^TXAX^TX] = \frac{m-1}{m^2} \Sigma{}A\Sigma{} + \frac{1}{m}(\Sigma(A+A^T)\Sigma + \Sigma \operatorname{Tr}(A\Sigma))$$ QED

Now, we look at the general case.

Lemma: Assume $$m$$ independent gaussian random vectors $$x^{(i)} \overset{iid}{\sim} \mathcal{N}(\mu,\Sigma)$$. Row-stack them in the matrix $$X$$. Let $$H=\Sigma+\mu\mu^T$$.

Then $$$$E\left[\frac{X^TX}{m}A\frac{X^TX}{m}\right] = \frac{(m-1)}{m}HAH + \frac{1}{m} \left( H(A+A^T)H - \mu\mu^T(A+A^T)\mu\mu^T + H\operatorname{Tr}(HA)) \right),$$$$

Proof: Let $$x \sim \mathcal{N}(\mu,\Sigma)$$. Define $$y=x-\mu$$.

As before, $$$$\frac{1}{m^2}E[X^TXAX^TX] = \frac{1}{m^2}( m(m-1) E[xx^T]AE[xx^T] + m E[xx^TAxx^T]) \label{eq:one}$$$$

The first term is simply $$$$\label{eq:two} E[xx^T]AE[xx^T] = (\Sigma+\mu\mu^T)A(\Sigma+\mu\mu^T).$$$$

The second term needs more work. We expand it in $$y+\mu$$, and notice that all terms with 1 or 3 factors of $$y$$ must have expectation 0 due to either Isserlis theorem, or to the zero mean. There are in total 16 terms, but 8 gets zeroes out in that way. $$$$\begin{split} E[xx^TAxx^T] &= E[(y+\mu) (y+\mu)^TA(y+\mu)(y+\mu)^T] \\ & = E[yy^TAyy^T] + E[\mu\mu^TA\mu\mu^T]\\ & + E[yy^TA\mu\mu^T] + E[\mu\mu^TAyy^T]\\ & + E[y\mu^TAy\mu^T] + E[\mu{}y^TA\mu{}y^T]\\ & + E[y\mu^TA\mu{}y^T] + E[\mu{}y^TAy\mu{}^T] \end{split}$$$$ All these are simple expectations. The first one due to Isserlis' theorem. The second is the expectation of a constant. Term three and four are just covariances. Term five and six are computed by the initial Lemma that I needed some index manipulations for. Term seven is due to the inner three factors are scalar so the $$y$$'s can be combined. Term eight uses $$y^TAy = \operatorname{Tr}(yy^TA)$$ and that the trace is linear and thus commutes with expectations. So, we find $$$$\begin{split} E[xx^TAxx^T] &= (\Sigma(A+A^T)\Sigma + \Sigma \operatorname{Tr}(A\Sigma)) + \mu\mu^TA\mu\mu^T\\ & + \Sigma{}A\mu\mu^T + \mu\mu^TA\Sigma{}\\ & + \Sigma{}A^T\mu\mu^T + \mu\mu^TA^T\Sigma \\ & + \Sigma\mu^TA\mu + \mu\mu^T\operatorname{Tr}(\Sigma{}A) \end{split}$$$$ Some regrouping of terms yields $$$$\begin{split} E[xx^TAxx^T] &= (\Sigma+\mu\mu^T)(A+A^T)(\Sigma+\mu\mu^T) - \mu\mu^T(A+A^T)\mu\mu^T\\ & + (\Sigma+\mu\mu^T)\operatorname{Tr}((\Sigma+\mu\mu^T)A) \end{split}$$$$

Combining the found results, we get $$$$E\left[\frac{X^TX}{m}A\frac{X^TX}{m}\right] = \frac{(m-1)}{m}HAH + \frac{1}{m} \left( H(A+A^T)H - \mu\mu^T(A+A^T)\mu\mu^T + H\operatorname{Tr}(HA)) \right),$$$$ where $$H=\Sigma+\mu\mu^T$$

QED

• dividing by $m$ simplifies a little. Assuming symmetry of A will also simplify the answer. Otherwise, I think this is about what to expect from the problem... Sep 22 at 14:10
• Thanks for going through the derivation, will try this formula out. It is indeed natural to assume symmetry of $A$, and to divide by $m^2$, since that demonstrates how formula interpolates between the simple formulas that exist for $m=1$ and $m=\infty$ Sep 22 at 14:23
• In differential geometry, this is known as a "debauch of indices." Applying a little bit of linear algebra at the outset to simplify the situation would shorten this demonstration considerably.
– whuber
Sep 22 at 15:00
• Thank you. Seems I did f' up then. Will double check and come back. Sep 22 at 18:25
• Woo hoo, works! Sep 22 at 20:38