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?

  • $\begingroup$ Are you sure you wrote the question correctly? A dimensional analysis easily shows the expectation must depend on $\Sigma,$ but $\Sigma$ doesn't appear in your expected answer. $\endgroup$
    – whuber
    Aug 15, 2018 at 15:18
  • $\begingroup$ Since the covariance matrix of $X$ is $\Sigma$, I expect the $\Sigma$ and $\Sigma^{-1}$ cancel out. So there is no $\Sigma$ in the answer. $\endgroup$
    – Eggplant
    Aug 19, 2018 at 11:11

1 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] }$$


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.