# 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 '20 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 '20 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... – Yaroslav Bulatov Sep 13 '20 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 '20 at 0:38
• Might want to change TensorProduct to CircleTimes here? – Carl Woll Sep 18 '20 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)
$$$$