Fourth order moments of the sum of multivariate normal distribution 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?
 A: 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.
