Expectation of product of quadratic forms in Gaussian distribution Let $x \in \mathbb{R}^n$ be a random vector with i.i.d. entries distributed as $\mathcal{N}(0, 1)$.  Let $A, B$ be two $n \times n$ symmetric matrices.  I would like to find $\mathbb{E} (x^T A x) (x^T B x)$.  
 A: That’s the expectation of a polynomial of degree $4$ in the $x_i$. Only the terms with even exponents for all $x_i$ survive. There are three types of those: $A_{ii}B_{ii}x_i^4$, $A_{ii}B_{jj}x_i^2x_j^2$ and $A_{ij}B_{ji}x_i^2x_j^2$ (where $i\ne j$). The moments are $\mathbb E\left[x_i^4\right]=3$ and $\mathbb E\left[x_i^2\right]=1$. Thus
\begin{eqnarray}
\mathbb E\left[(x^\top Ax)(x^\top Bx)\right]
&=&
\sum_{ijkl}x_iA_{ij}x_jx_kB_{kl}x_l
\\
&=&
3\sum_iA_{ii}B_{ii}+\sum_{i\ne j}A_{ii}B_{jj}+2\sum_{i\ne j}A_{ij}B_{ji}
\\
&=&
\sum_{ij}A_{ii}B_{jj}+2\sum_{ij}A_{ij}B_{ji}
\\
&=&
\operatorname{Tr}A\operatorname{Tr}B+2\operatorname{Tr}AB\;.
\end{eqnarray}
Note that the result is basis-independent, as it must be, since applying a unitary transform to $x$ again yields a vector of components with independent standard normal distributions.
A: Fact 1: Let $S$ be a (real) symmetric matrix. Then $\mathbb{E}[(x^TSx)^2] = 2\mathrm{Tr}(S^2) + (\mathrm{Tr}(S))^2$.
From Fact 1, by using $4ab = (a+b)^2 - (a-b)^2$, we have
\begin{align}
\mathbb{E}[(x^TAx)(x^TBx)]
&= \frac{1}{4}\mathbb{E}[(x^T(A+B)x)^2] - \frac{1}{4}\mathbb{E}[(x^T(A-B)x)^2]\\
&= \frac{1}{2}\mathrm{Tr}((A+B)^2) + \frac{1}{4}(\mathrm{Tr}(A+B))^2 \\
&\quad - \frac{1}{2}\mathrm{Tr}((A-B)^2)
- \frac{1}{4}(\mathrm{Tr}(A-B))^2\\
&= 2\mathrm{Tr}(AB) + \mathrm{Tr}(A)\mathrm{Tr}(B).
\end{align}
$\phantom{2}$
Proof of Fact 1: Let $S = U\mathrm{diag}(d_1,d_2, \cdots, d_n)U^T$ be the eigendecomposition.
Let $y = U^Tx$. Then $y$ is a random vector with i.i.d. entries distributed as $\mathcal{N}(0, 1)$.
Note that $\mathbb{E}[y_i^4] = 3$ and $\mathbb{E}[y_i^2]=1$.
We have
\begin{align}
\mathbb{E}[(x^TSx)^2] &= \mathbb{E}[(y^T\mathrm{diag}(d_1,d_2, \cdots, d_n)y)^2]\\
&= \mathbb{E}\Big[\sum_{i=1}^n d_i^2 y_i^4 + \sum_{1\le i < j\le n} 2d_id_jy_i^2y_j^2\Big]\\
&= 3\sum_{i=1}^n d_i^2 + \sum_{1\le i < j\le n} 2d_id_j\\
&= 2\mathrm{Tr}(S^2) + (\mathrm{Tr}(S))^2.
\end{align}
