Fourth order moments of the sum of multivariate normal distribution

Question

Suppose $X\sim \mathcal{N}_d(\boldsymbol{0},\boldsymbol{\Sigma})$ follows a $d$-dimensional multivariate normal distribution. Let $\text{A}_i$ and $\text{A}_j$ be two arbitrary symmetric $d\times d$ matrices. What is the expectation $\mathbf{E}\left[X^T\text{A}_i\boldsymbol{\Sigma}^{-1}XX^T\text{A}_j\boldsymbol{\Sigma}^{-1}X\right]$? I expect the answer is $3\,\text{tr}(\text{A}_i\text{A}_j)$, but how to prove it?

Answer 1

Let's solve the problem first with the identity matrix $\Sigma.$ To reduce subscripts, replace $A_i$ and $A_j$ by $d\times d$ symmetric matrices $A=(a_{ij})$ and $B=(b_{kl})$. Bearing in mind that for any subscript $i=1,2,\ldots, d$ the expectations of these Normal variables are

$$0 = E[x_i] = E[x_i^3];\ 1 = E[x_i^2];\ 3 = E[x_i^4]\tag{*}$$

we may compute (from linearity of expectation) that

$$\eqalign{ E[x^\prime A x\, x^\prime B x] &= E[(x^\prime A x)\, (x^\prime B x)] \\ &= E\left[\sum_{i,j} x_i a_{ij} x_j\, \sum_{k,l} x_k b_{kl} x_l\right] \\ &= \sum_{i,j,k,l} a_{ij}b_{kl}\,E\left[x_ix_jx_kx_l\right]. }$$

From $(*)$ and the independence of the $x_i$ (because $\Sigma$ is the identity) the only nonzero terms in these expectations have only even powers of individual $x_i,x_j,x_k,x_l.$ One way to work out the result is to return to the middle expression and extract the squares, writing

$$\sum_{i,j} x_i a_{ij} x_j = \sum_i a_{ii}x_i^2 + 2\sum_{i\lt j} a_{ij} x_i x_j$$

and likewise for the second factor. Expanding the products gives monomials in $x_i^2 x_k^2,$ $x_ix_j x_k^2,$, $x_i^2x_lx_k,$ and $x_ix_jx_kx_l$ where $i\lt j$ and $k\lt l$. Upon taking expectations the middle two will become zero. The last will be nonzero only when $i=k$ and $j=l.$ The first will equal $1$ except when $i=k,$ in which case it will equal $3:$ that is, $2$ is added to what we might otherwise think the expectation is. What remains is

$$\eqalign{ E\left[\sum_{i,j} x_i a_{ij} x_j\, \sum_{k,l} x_k b_{kl} x_l\right] &= E\left[\sum_{i,k} a_{ii}b_{kk}\,x_i^2 x_k^2 + 4 \sum_{i\lt j} a_{ij} b_{ij}\, x_i^2 x_j^2\right]\\ &= \sum_{i,k} a_{ii}b_{kk} + 2\sum_{i} a_{ii}b_{ii} + 2\left(\sum_{i,j} a_{ij}b_{ij} - \sum_{i} a_{ii}b_{ii}\right)\\ &= \sum_i a_{ii}\,\sum_k b_{kk} + 2\sum_{i,j} a_{ij}b_{ij} \\ &= {\operatorname{Tr}(A)\operatorname{Tr}(B) + 2 \operatorname{Tr}(A B).\tag{**}} }$$

(The last term is obviously $\operatorname{Tr}(A^\prime B),$ but the symmetry permits us to replace $A^\prime$ by $A.$ This will be useful in the last step below.)

---

When $\Sigma$ is not the identity matrix, it has a square root $\Lambda:$ that is, $$\Lambda\Lambda^\prime = \Sigma,$$ and $x$ has the distribution of $\Lambda z$ for standard multivariate Normal $z$ because (since all expectations of the $z_i$ and $x_i$ are zero and $E[zz^\prime]$ is the identity matrix)

$$\operatorname{var}{\Lambda z} = E[(\Lambda z)(\Lambda z)^\prime] = E[\Lambda\,zz^\prime\,\Lambda^\prime] = \Lambda\Lambda^\prime=\Sigma.$$

Thus, addressing the expression in the question and noting that

$$\Lambda^\prime \Sigma^{-1} \Lambda = 1_d,$$

$$\eqalign{E[x^\prime A \Sigma^{-1} x x^\prime B \Sigma^{-1} x] &= E[z^\prime \Lambda^\prime\, A \Sigma^{-1}\,\Lambda z\, z^\prime \Lambda^\prime\, B \Sigma^{-1}\, \Lambda z] \\ &= E[z^\prime (\Lambda^\prime\ A (\Lambda^\prime)^{-1})(\Lambda^\prime\ \Sigma^{-1}\Lambda) z\, z^\prime )(\Lambda^\prime\, B (\Lambda^\prime)^{-1}) (\Lambda^\prime \Sigma^{-1}\Lambda) z] \\ &=E[z^\prime A^{\Lambda} z\, z^\prime B^{\Lambda} z] }$$

where

$$A^{\Lambda} = \Lambda^\prime\ A (\Lambda^\prime)^{-1}$$

is the $\Lambda^\prime$ conjugate of $A$ (and likewise with $B$). Since the trace is invariate under conjugation, and $(AB)^{\Lambda}=A^\Lambda B^\Lambda,$ the original formula $(**)$ remains unchanged.