How to evaluate $\mathrm{Cov}(\vec x^T A \vec x,\vec x^T B \vec x )$ for $\vec x\sim N(\vec 0, \Sigma)$ and square matrices $A,B$?

Question

I need to evaluate $\mathrm{Cov}(\vec x^T A \vec x,\vec x^T B \vec x )$ for $\vec x\sim N(\vec 0, \Sigma)$ and square (nilpotent) matrices $A,B$ (not symmetric, non-random). $\Sigma$ is known. It becomes incredibly messy if I write it out, but I believe the final result will be quite tidy because of numerics. This makes me speculate that there exists some identity that makes it easier to evaluate, i.e. by somehow factoring the $A$ and $B$ matrices out of the covariance in exchange for more terms or something similar?

Answer 1

Let $\tilde{A}=\frac12(A+A^T)$, then $\tilde{A}$ is symmetric and $\vec{x}^T\tilde{A}\vec{x}=\vec{x}^TA\vec{x}$.
For symmetric $\tilde{A},\tilde{B}$ and $\vec{x}\sim N(\vec{0},\Sigma)$, $ \mathsf{Cov}(\vec{x}^T\tilde{A}\vec{x},\vec{x}^T\tilde{B}\vec{x})=2\,\text{tr} (\tilde{A}\Sigma \tilde{B}\Sigma)$.