Expectation of product of quadratic forms in Gaussian distribution

Question

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)$.

Answer 1

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.

Answer 2

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}