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

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]


• I wonder if first solving for $E[X^T A X B X^T X]$ and $E[X^T X B X^T CX]$ would help lead one to a general formula for $E[X^T A X B X^T CX]$.
– JimB
Sep 24, 2022 at 22:51
• Your question involves a multivariate normal with a diagonal covariance matrix (rather than the general covariance matrix in the link you presented) and so using Moment is pretty fast for your example.
– JimB
Sep 24, 2022 at 23:55
• Some speculation: If $X^T.A.X.B.X^T.C.X==X^T.C.X.B.X^T.A.X$, then one could target solving $X^T.D.X.B.X^T.D.X$ as $X^T.A.X.B.X^T.C.X = (X^T.(A+C).X.B.X^T.(A+C).X - X^T.A.X.B.X^T.A.X - X^T.C.X.B.X^T.C.X)/2$.
– JimB
Sep 25, 2022 at 23:04

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 :)

• Does your formula work? Isserlis theorem is for centered variables. I haven't seen the non-centered version, this is why I used MomentConvert in the notebook linked Sep 26, 2022 at 8:09
• edited the question to clarify that I'm looking for generic version, not necessarily centered Sep 26, 2022 at 8:12
• +1 Excellent! I think you have a typo in $E[Z^TAZBZ^TCZ]=(\text{tr} A)(\text{tr} C) HBH + HB^T H \text{tr} (A^T C) + \text{tr} (BH) \text{tr} (A^TC) H$ in the the last $\text{tr} (A^TC)$ should be $\text{tr} (A C)$ (i.e., no transpose for $A$). The Mathematica code would be Tr[A] Tr[CC] H . B . H + H . Transpose[B] . H Tr[Transpose[A] . CC] + Tr[B . H] Tr[A . CC] H (where I've used CC instead of C because the symbol C is reserved).
– JimB
Sep 26, 2022 at 16:13
• Two items: (1) This matrix solution is tens (or hundreds) of thousands times faster than brute force approach using replacement rules, and (2) you should post (or would it be cross-post?) your answer for Yaroslav's bounty on CrossValidated.
– JimB
Sep 26, 2022 at 16:26
• Thanks @JimB, I edited. I think I will stick to this friendly SE for the time being :) Sep 26, 2022 at 16:29

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

• The "how to get matrix form of moment" seems similar to "how to get matrix form of derivative". Here Sjoerd does this using replacement rules Sep 24, 2022 at 23:37