Accelerare le aspettative gaussiane

Question

Ho scoperto che l'operatore Expectation di Mathematica è lento per le variabili casuali gaussiane, mi piacerebbe scalare a 100 dimensioni, ma attualmente ci vogliono 10 secondi per 10 dimensioni, in attesa di suggerimenti su come accelerare questo processo.

Problema: per una variabile casuale gaussiana $x$ nel $d$ dimensioni, calcolare le seguenti quantità $$E[xx'],E[xx'xx'],E[xx'\otimes xx']$$

$\otimes$ si riferisce a 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

Le identità consentono di calcolare queste matrici centinaia di volte più velocemente $$ \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) 
```