Prove that $\mathrm{Cov}(x^TAx,x^TBx) = 2 \mathrm{Tr}(A \Sigma B \Sigma) + 4 \mu^TA \Sigma B \mu$

Question

Suppose $\vec x \sim N(\vec \mu, \Sigma)$ is multivariate normal. I want to see that $$\mathrm{Cov}(\vec x^TA\vec x,\vec x^TB\vec x) = 2 \mathrm{Tr}(A \Sigma B \Sigma) + 4 \vec \mu^TA \Sigma B \vec \mu$$ I have been searching the internet for a while, and found multiple sources confirming it, but no proof, for example here: $Var(Q)=2\ tr(A\Sigma A\Sigma)+4\mu^TA\Sigma A\mu$

Do you know how to show it?

Answer 1

We have

\begin{align} \operatorname{Var}(x^T(A+B)x)&=\operatorname{Var}(x^TAx+x^TBx) \\&=\operatorname{Var}(x^TAx)+\operatorname{Var}(x^TBx)+2\operatorname{Cov}(x^TAx,x^TBx) \end{align}

So, $$\operatorname{Cov}(x^TAx,x^TBx)=\frac12\left[\operatorname{Var}(x^T(A+B)x)-\operatorname{Var}(x^TAx)-\operatorname{Var}(x^TBx)\right]\tag{1}$$

Note that for the calculation of variance and covariance when $x$ is multivariate normal, it is assumed that $A$ and $B$ are symmetric, which also makes $A+B$ symmetric.

For $y\sim N(0,I)$, it is shown here that $$\operatorname{Var}(y^TAy)=2\operatorname{tr}(A^2)\tag{2}$$

Now suppose $x\sim N(\mu,\Sigma)$ where $\Sigma$ is positive definite.

Then there exists a nonsingular matrix $C$ such that $\Sigma=CC^T$.

And $$x\sim N(\mu,\Sigma)\implies y=C^{-1}(X-\mu)\sim N(0,I)$$

So, $$x^TAx=(\mu+Cy)^TA(\mu+Cy)=\mu^TA\mu+2\mu^TACy+y^T(C^TAC)y$$

Therefore noting that $y$ and $y^T(C^TAC)y$ are uncorrelated,

\begin{align} \operatorname{Var}(x^TAx)&=4\operatorname{Var}(\mu^TACy)+\operatorname{Var}\left(y^T(C^TAC)y\right) \\&=4\mu^T(AC)(AC)^T\mu + 2\operatorname{tr}((C^TAC)(C^TAC))&\small\left[\text{ using }(2)\right] \\&=4\mu^T A\Sigma A\mu + 2\operatorname{tr}(C^TA\Sigma AC) \\&=4\mu^T A\Sigma A\mu + 2\operatorname{tr}(A\Sigma A\Sigma) &\small\left[\because \,\operatorname{tr}(AB)=\operatorname{tr}(BA)\right] \tag{3} \end{align}

Now it follows from $(1)$ and $(3)$ that

$$\operatorname{Cov}(x^TAx,x^TBx)=4\mu^TA\Sigma B\mu + 2\operatorname{tr}(A\Sigma B\Sigma)$$

For a direct calculation of the covariance and hence the variance, one may refer to Graybill's Matrices with Applications in Statistics, 2nd edition (1983).

Answer 2

Just trying to provide a sketch for the proof, the key here is to show the following equations: $$\begin{align}\operatorname{E}[x'Ax]&=\operatorname{tr}(A\Sigma)+\mu'A\mu, \label{1}\tag{1}\\\operatorname{E}[x'Axx'Bx]&=2\operatorname{tr}(A\Sigma B\Sigma)+4\mu'A\Sigma B\mu+(\operatorname{tr}(A\Sigma)+\mu'A\mu)(\operatorname{tr}(B\Sigma)+\mu'B\mu). \label{2}\tag{2}\end{align}$$ The desired result follows immediately since by definition $$\begin{align}\operatorname{Cov}(x'Ax, x'Bx)&=\operatorname{E}[(x'Ax-E[x'Ax])(x'Bx-\operatorname{E}[x'Bx])']\\&=\operatorname{E}[x'Axx'Bx]-\operatorname{E}[x'Ax]\operatorname{E}[x'Bx].\end{align}$$

The proof to equation ($\ref{1}$) is simple and can be found in many introductory text. Equation ($\ref{2}$) is the real deal here, but fortunately a proof can be found in Proofs Section 5 of the Matrix Reference Manual. Check out 5.18 and 5.19 for Isserlis' theorem, and finally 5.28 where they derived an expression for a much more general form: $$\operatorname{E}[(Ax-a)'(Bx-b)(Cx-c)'(Dx-d)].$$