Obtaining formula of $E[X^T A X B X^T CX]$ for Gaussian $X$

Question

Suppose $X$ is an $m\times n$ random matrix where rows are I.I.D. samples $x$ of $n$-dimensional Gaussian. The code at the end of the post verifies a candidate (incorrect) formula for the following expression in terms of $H=E[xx']$ (this parameterization simplified related formula considerably compared to $\mu,\Sigma$ parameterization). What's the correct formula?

$$E[X^T A X B X^T CX]\approx m(m-1)HBH$$

There's this detailed answer on computing arbitrary moments of Gaussian distributions.

Alternatively there's this notebook which automatically derives non-central Wick's Theorem in summation form by doing the following:

1. MomentConvert[Moment[{1, 1, 1, 1}], "Cumulant"]
2. Set cumulants of order 3 and 4 to zero
3. MomentConvert back

A similar approach may be feasible for the question at hand, but a bit of a programming challenge!

d = 3;
x1 = Array[xx1, d];
x2 = Array[xx2, d];
X = {x1, x2};
SeedRandom[1];
sigma = DiagonalMatrix[RandomInteger[{1, 10}, {d}]];
mu = RandomInteger[{1, 10}, {d}];
mu = ConstantArray[0, {d}];
H = sigma + Outer[Times, mu, mu];
dist = MultinormalDistribution[mu, sigma];
EX[expr_] := 
  Expectation[
   expr, {x1 [Distributed] dist, x2 [Distributed] dist}];
m = Length[X];
Unprotect[C];
{A, C} = RandomInteger[{-5, 5}, {2, m, m}];
B = RandomInteger[{-5, 5}, {d, d}];
Print["Check passed: ", 
 EX[X[Transpose] . A . X . B . X[Transpose] . C . X] == 
  m (m - 1) H . B . H]

Related unanswered question on stats.SE

Answer 1

Note: Many thanks to @JimB for pointing out a mistake in an earlier version.

This question can be worked out by pencil and paper.

Special case. Let me first assume that $Y$ is a rectangular matrix whose entries are i.i.d. and standard normal, $N(0,1)$. Then: $$ E[Y^TAYBY^TCY] = (\text{tr} A)(\text{tr} C) B + B^T \text{tr} (A^T C) + \text{tr} (B) \text{tr} (AC) \mathbb{1} $$ The proof goes like this: $$ E[Y^TAYBY^TCY]_{ij} = \sum_{a,b,c,d,e,f} E[Y_{ai} A_{ab} Y_{bc} B_{cd} Y_{ed} C_{ef} Y_{fj}] = \sum_{a,b,c,d,e,f} A_{ab} B_{cd} C_{ef} E[Y_{ai} Y_{bc} Y_{ed} Y_{fj}] $$ Using Isserli's theorem this is equal to $$ = \sum_{a,b,c,d,e,f} A_{ab} B_{cd} C_{ef} (\delta_{ab} \delta_{ic} \delta_{ef} \delta_{dj} + \delta_{ae} \delta_{id} \delta_{bf} \delta_{cj} + \delta_{af} \delta_{ij} \delta_{be} \delta_{cd}) $$ where $\delta$ is the Kronecker symbol. It is now just a matter of checking what is contracted with what, for example in the first term we have $A_{ab} \delta_{ab}$ which gives $\text{tr} A$, and so on. The result follows.

By a similar calculation: $$ E[YAY] = E[Y^T A Y^T] = A^T \qquad E[YAY^T] = E[Y^T AY] = \text{tr}(A) \mathbb{1} $$

General case (but still zero mean). Now assume that $Z = YK$ where $K$ is some constant matrix, and denote $H = K^T K$. Then $$ E[Z^TAZBZ^TCZ] = (\text{tr} A)(\text{tr} C) HBH + HB^T H \text{tr} (A^T C) + \text{tr} (BH) \text{tr} (AC) H $$

The proof is, using the special case above: $$ E[Z^TAZBZ^TCZ] = E[K^TY^TAYKBK^T K^T Y^TCYK] = K^T E[Y^TAY(KBK^T) Y^TCY]K = (\text{tr} A)(\text{tr} C) K^TKBK^TK + K^T KB^T K^TK \text{tr} (A^T C) + \text{tr} (KBK^T) \text{tr} (AC) K^T K $$ and then replacing $K^T K = H$.

By a similar calculation: $$ E[ZAZ] = A^T H \qquad E[Z^T A Z^T] = H A^T$$ $$ E[ZAZ^T] = \text{tr}(AH) \mathbb{1} \qquad E[Z^T AZ] = \text{tr}(A) H $$

General case (non-zero mean). Suppose $M$ is the mean, $E[X] = M$. This is equivalent to $X = Z + M$ where $Z$ is as before. Then $$ E[X^TAXBX^TCX]\\ = (\text{tr} A)(\text{tr} C) HBH + HB^T H \text{tr} (A^T C) + \text{tr} (BH) \text{tr} (AC) H\\ + \text{tr}(A) HBM^TCM\\ + H (AMB)^T CM\\ + \text{tr}(BH) M^TACM\\ + \text{tr}(AMBM^TC) H \\ + M^TA (BM^TC)^T H\\ + \text{tr}(C) M^TAMB H \\ + M^TAMBM^TCM $$

The proof goes like this (expectation values of terms that are linear or cubic in $Z$ are equal to zero): $$ E[X^TAXBX^TCX]\\ = E[Z^TAZBZ^TCZ]\\ + E[Z^TAZBM^TCM]\\ + E[Z^TAMBZ^TCM]\\ + E[M^TAZBZ^TCM]\\ + E[Z^TAMBM^TCZ]\\ + E[M^TAZBM^TCZ]\\ + E[M^TAMBZ^TCZ]\\ + E[M^TAMBM^TCM]\\ = E[Z^TAZBZ^TCZ]\\ + E[Z^TAZ]BM^TCM\\ + E[Z^T(AMB)Z^T]CM\\ + M^TAE[ZBZ^T]CM\\ + E[Z^T(AMBM^TC)Z]\\ + M^TAE[Z(BM^TC)Z]\\ + M^TAMBE[Z^TCZ]\\ + M^TAMBM^TCM $$ Now it is just a matter of plugging in the results given for $Z$.

Note. It would be easy to make a mistake in such a calculation, so use at your own risk :)

Answer 2

This is just an extended comment in that one can find the expectation much quicker by using replacement rules rather than using Expectation. (I understand that the real answer must be in matrix terms.)

d = 3; (* Dimensions *)
m = 2; (* Number of independent samples *)
X = Table[x[i, j], {i, m}, {j, d}];
H = DiagonalMatrix[Range[d]];  (* Covariance matrix *)
dist = MultinormalDistribution[H];
A = Array[a, {d, d}];
W = X[Transpose] . X . A . X[Transpose] . X;
expectations = {x[i_, j_]^3 -> 0, x[i_, j_]^4 -> 3 j^2, x[i_, j_]^2 x[i_, k_]^2 -> j k, 
  x[i_, j_]^2 -> j, x[i_, j_] -> 0};
AbsoluteTiming[(X[Transpose] . X . A . X[Transpose] . X // Expand) /. expectations;]
(* {0.0028874, Null} *)

It takes 1.4 seconds using Expectation with your code.

The expectations in the replacement list are found using Moments and recognizing that each random term (after Expand) is a product of powers of 1, 2, or 3 elements of the $X$ matrix such as $x_{1,1}^2 x_{1,3}^2 and $x_{2,1} x_{2,3}^3. And the expectation of any of those products with any term with an odd power is zero.

So with $d=3$ we have the expectation of $x_{1,1}^2 x_{1,3}^2$ being

Moment[dist, {2, 0, 2}]
(* 3 *)

For $x_{2,2}^4$ the expectation is

Moment[dist, {0, 4, 0}]
(* 12 *)

And the replacement list was constructed by noticing patterns. (Unfortunately the order of the replacements is important as I need to replace x[i, j]^2 before replacing x[i, j] with zero.)