Speeding up Gaussian expectations

Question

I found that Mathematica's Expectation operator is slow for Gaussian random variables, I'd like to scale to 100 dimensions, but it currently takes 10 seconds for 10 dimensions, looking forward to tips on how to speed this up.

Problem: for a Gaussian random variable $x$ in $d$ dimensions, compute the following quantities $$E[xx'],E[xx'xx'],E[xx'\otimes xx']$$

$\otimes$ refers to KroneckerProduct

CircleTimes = KroneckerProduct;

problemSetup[d_] := (
   cov = DiagonalMatrix@Table[1/k, {k, 1, d}];
   dist = MultinormalDistribution[N@cov];
   Clear[x];
   xvec = Array[x, d];
   X2 = Expectation[xvec\[CircleTimes]xvec, 
     xvec \[Distributed] dist];
   X2X2 = 
    Expectation[(xvec\[CircleTimes]xvec).(xvec\[CircleTimes]xvec),
      xvec \[Distributed] dist]; 
   X4 = Expectation[Outer[Times, xvec, xvec, xvec, xvec], 
     xvec \[Distributed] dist];
   X4flat = Flatten[X4, {{1, 2}, {3, 4}}];
   );
problemSetup[10] // Timing (* {10.3734, Null} *)
```

Answer 1

Identities allow computing these matrices 100s of times faster $$ \begin{eqnarray} P_d&=&\text{nsymm[d], the symmetrizer matrix}\\ \mu&=&E[x]\\ X^2&=&E[xx']\\ E[(xx')(xx')]&=&2X ^2+X^2 \text{Tr} X^2 -2\|\mu\|^2 \mu \mu'\\ E[xx'\otimes xx']&=&2P_d( X^2\otimes X^2) + \text{vec} X^2 (\text{vec} X^2)'-2(\mu \mu')\otimes (\mu \mu') \\ \end{eqnarray} $$

(* Symmetrizer matrix, see \
https://mathematica.stackexchange.com/questions/230167/commutation-\
symmetrizer-and-duplication-matrices *)

getperm[perm_, n_] := Permute[IdentityMatrix[n*n], perm]
kcomm[n_] := 
 With[{mtx = ArrayReshape[Range[n*n], {n, n}]}, 
  getperm[FindPermutation[vec[Transpose[mtx]], vec[mtx]], Length[mtx]]]

nsymm[n_] := (kcomm[n] + IdentityMatrix[n^2])/2
vec[W_] := Join @@ Transpose[W];

CircleTimes = KroneckerProduct;
Clear[xx];
d = 4;
x = Array[xx, d];
sigma = RandomInteger[{-5, 5}, {d, d}];
sigma = sigma.Transpose[sigma] + 
  IdentityMatrix[
   d]; (* strictly posdef to avoid crash in Expectation *)
dist = 
 MultinormalDistribution[RandomInteger[{-5, 5}, {d}], sigma];
Ex[expr_] := 
 Expectation[expr, 
  x \[Distributed] dist]; (* Expectation with respect to x *)

X2 = Ex[x⊗x]; (* same as sigma+mu⊗mu *)
mu = Ex[x];
Ex[(x⊗x).(x⊗x)] == Tr[X2] X2 + 2 X2.X2 - 2 mu.mu (mu⊗mu)
Ex[(x⊗x)⊗(x⊗x)] == 2 nsymm[d].(X2⊗X2) + vec[X2]⊗vec[X2] - 2 (mu⊗mu)⊗(mu⊗mu) 
```