# Speeding up Gaussian expectations

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 // Timing (* {10.3734, Null} *)
$$$$

• Are more general covariance matrices to be considered? What you've described are independent normal distributions with mean zero and different variances.
– JimB
Sep 12, 2020 at 23:24
• Just fiddling, I've gotten the time eater (X4) from 11.7 sec to 0.03 sec. Doing it for 100 takes about 5 minutes. Is that fast enough for your needs?
– ciao
Sep 12, 2020 at 23:53
• @jimb just zero centered Gaussians with arbitrary diagonal covariances (since any covariance matrix is diagonal in some coordinate system), ciao I was hoping for 5 seconds... Sep 13, 2020 at 0:06
• Just getting Outer[Times, xvec, xvec, xvec, xvec] probably takes 2 minutes with d=100. I suspect that @ciao is using the following: AbsoluteTiming[Outer[Times, xvec, xvec, xvec, xvec] /. {x[k_]^4 -> 3/k^2, x[k_]^3 -> 0, x[k_]^2 -> 1/k, x[k_] -> 0};].
– JimB
Sep 13, 2020 at 0:38
• Might want to change TensorProduct to CircleTimes here? Sep 18, 2020 at 16:26

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