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

enter image description here

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};

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_] := 
   expr, {x1 \[Distributed] dist, x2 \[Distributed] dist}];
m = Length[X];
{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

  • $\begingroup$ 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]$. $\endgroup$
    – JimB
    Commented Sep 24, 2022 at 22:51
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Commented Sep 24, 2022 at 23:55
  • 1
    $\begingroup$ 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$. $\endgroup$
    – JimB
    Commented Sep 25, 2022 at 23:04

3 Answers 3


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

  • $\begingroup$ 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 $\endgroup$ Commented Sep 26, 2022 at 8:09
  • $\begingroup$ edited the question to clarify that I'm looking for generic version, not necessarily centered $\endgroup$ Commented Sep 26, 2022 at 8:12
  • 1
    $\begingroup$ +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). $\endgroup$
    – JimB
    Commented Sep 26, 2022 at 16:13
  • $\begingroup$ 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. $\endgroup$
    – JimB
    Commented Sep 26, 2022 at 16:26
  • $\begingroup$ Thanks @JimB, I edited. I think I will stick to this friendly SE for the time being :) $\endgroup$
    – user293787
    Commented 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.)

  • 2
    $\begingroup$ 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 $\endgroup$ Commented Sep 24, 2022 at 23:37

I wrote a program (tensorgrad) that allows doing this kind of computation. It uses Gaussian integration by parts (basically the same as Isserlis' theorem), and allows easily computing:

\begin{align} E[X A X^T B X C X^T] &=B^T \cdot \mathrm{tr}( C) \cdot \mathrm{tr}(A ) \\&+ \mathrm{tr}(A C) \cdot B \\&+ \mathrm{tr}(A C^T) \cdot \mathrm{tr}( B) \cdot \mathbf{1} \end{align} for the mean 0, $\Sigma=I$ case, and

\begin{align} E[X A X^T B X C X^T] &=M A M^T B M C M^T \\&+ M A M^T B \cdot \mathrm{tr}(S C) \\&+ M A S C M^T \cdot \mathrm{tr}( B) \\&+ M A S C^T M^T B^T \\&+ B M C M^T \cdot \mathrm{tr}(A S) \\&+ B \cdot \mathrm{tr}(S C) \cdot \mathrm{tr}(A S) \\&+ B^T M A^T S C M^T \\&+ B^T \cdot \mathrm{tr}(S A^T S C) \\&+ \mathrm{tr}(C M^T B^T M A^T S) \cdot \mathbf{1} \\&+ \mathrm{tr}(S A^T S C^T) \cdot \mathrm{tr}( B) \cdot \mathbf{1} \end{align}

Where $E[X] = M$ and $S$ is the covariance matrix of each row.

I tested it on 1,000,000 samples to give the right expectation.

This shows that you don't need to parametrize in terms of $E[xx']$ to get a simple output. While using the covariance matrix, I believe this is as simple as user293787's solution.

Expressed in terms of tensor graphs, the formula looks like this:

enter image description here


Your Answer

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

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