How to efficiently find moments of a multinormal distribution?

If $X = (X_1, ..., X_n)$ denotes a $n$-variate multivariate Normal random variable $N(\mu, \Sigma)$, then the moment generating function is given by:

$$ M(t)=E\left[e^{t'X}\right] = exp\left({t'\mu + \frac12 t' \Sigma t}\right)$$

The mgf can then be used to derive the raw moments you seek. For instance, in say a bivariate setting, the product raw moment $E\left[X_1^r X_2^s\right]$ can be obtained from the mgf $M(t)$ as follows:

$$E\left[X_1^r X_2^s\right] = \frac{\partial ^{r+s}M(t)}{\partial t_1^r\partial t_2^s}|_{\overset{\rightharpoonup }{t}=\overset{\rightharpoonup }{0}}$$

Advantages

This approach has 3 advantages:

FIRST, it is much faster than using Wolfram's in-built Moment function.

SECOND, the Isserlis theorem the OP has quoted above requires that the mean vector $\mu$ = 0. By contrast, the mgf approach given here does not require a zero mean vector, and is thus both more general and flexible.

THIRD, comparing to the solution given by @ybeltukov (using Isserlis), the LeafCount of the obtained solution is dramatically smaller using our mgf approach ... e.g. 261 using mgf here versus a leaf count of nearly 3 billion for ybeltukov (see the example below).

Implementation is easy:

For say a bivariate model:

tt = {Subscript[t, 1], Subscript[t, 2]};  
mu = {0, 0}; 
SS = Array[σ, {2, 2}];

Then, the mgf is;

mgf = Exp[tt.mu + 1/2 tt.SS.tt]

Write a little function to extract the raw moments from the mgf:

RawMoment[r_, s_] := D[mgf, {Subscript[t,1], r}, {Subscript[t,2], s}] /. Subscript[t,_] -> 0

All done.

Example: find $E\left[X_1^{40} X_2^{30}\right]$

Solution:

 RawMoment[40, 30];  // AbsoluteTiming

{0.130264, Null}

compared to:

Moment[MultinormalDistribution[{0, 0}, SS], {40, 30}]; // AbsoluteTiming

{1.938334, Null}

... which is about 8 times faster. The output is otherwise the same.

Comparing to ybeltukov

In another solution, ybeltukov suggested a function called moment (using the OP's Isserlis method). While this is fast, it unfortunately produces ungainly answers. To see this, compare the LeafCount:

RawMoment[40, 30] // LeafCount

261

versus

moment[{40, 30}] // LeafCount

2886864330

For more detail

For more detail, see, for instance, Chapter 6 of our book, "Mathematical Statistics with Mathematica". A free download of the chapter is now available here:

http://www.mathStatica.com/book/Rose_and_Smith_2002edition_Chapter6.pdf

[ We will be making the entire original first edition of the Springer book available as a free download in the next few days ... ]

Pairs is unordered and index in each pair is also unordered. Let us choose the order:

1. In each pair the first index is less then the second.

2. Pairs is sorted in ascending order by the first index. If the first indexes is equal then sort by the second.

Examples in the question have the right order.

Let we have powers $(0,\ldots,0,r_i,r_{i+1},\ldots,r_n)$. The first pair is always $\sigma_{i,j}$ where $i$ is the first nonzero power $r_i$. $j\ge i$ can take any index of nonzero powers in the residual list $(0,\ldots,0,r_i-1,r_{i+1},\ldots,r_n)$. So we sort out all the possible first pairs and recursively calculate moment with the residual powers $(0,\ldots,0,r_i-1,r_{i+1},\ldots,r_{j-1},r_j-1,r_{j+1},\ldots,r_n)$

ClearAll[moment]
moment[{0 ...}] = 1;
moment[x_List] /; Total[x] == 2 := 
   If[Length[#] == 1, σ[#[[1]], #[[1]]], σ[#[[1]], #[[2]]]] &@
      Position[x, y_?Positive][[All, 1]];
moment[x_List] := moment[x] = Module[{i = LengthWhile[x, # == 0 &] + 1, x2 = x}, x2[[i]]--; 
   Sum[If[x2[[j]] > 0, x2[[j]] σ[i, j] moment@MapAt[# - 1 &, x2, j], 0], 
      {j, i, Length[x]}]];

I use memoization here because it produces a great speedup.

moment[{2, 2}]

2 σ[1, 2]^2 + σ[1, 1] σ[2, 2]

m = moment[{400, 0}]; // AbsoluteTiming // First

0.025246

md[n_] := MultinormalDistribution[ConstantArray[0, n], 
     Table[σ[Min[i, j], Max[i, j]], {i, n}, {j, n}]];

m == Moment[md[2], {400, 0}] // AbsoluteTiming

{2.011387, True}

100x speedup!

 d = 2 RandomInteger[5, 5]

{4, 0, 2, 6, 10}

 m = moment[d]; // AbsoluteTiming // First

0.012637

 m == Moment[md[5], d] // AbsoluteTiming

{3.172868, True}

Even greater speedup!

Update: To decrease the leaf count add //Expand just after the MapAt[# - 1 &, x2, j]. It will be ~2 times slower but still faster then wolfies' RawMoment.

Explicit formula

The wolfies' answer gave me an idea that one can derive an explicit formula. Here it is!

$$ E(x_1^{r_1}x_2^{r_2}\cdots x_n^{r_n}) = \sum_{(p)}\prod_{i}\frac{r_i!\,\sigma_{ii}^{p_{ii}}}{(2p_{ii})!!}\prod_{i<j}\frac{\sigma_{ij}^{p_{ij}}}{p_{ij}!} $$

where sum is performed over all non-negative integer values of $p_{ii}$ and $p_{ij}$ with constrain

$$ \sum_{j=1}^{i-1}p_{ji}+2p_{ii}+\sum_{j=i+1}^{n}p_{ij}=r_i, \quad i=1,2,\ldots,n. $$

Implementation

moment2[x_List] := With[{n = Length[x]}, 
  With[{σii = Table[σ[i, i], {i, n}], σij = Join @@ Table[σ[i, j], {i, n}, {j, i + 1, n}],
    pij = Table[Unique[], {n (n - 1)/2}], 
    pos = n (n - 1)/2 - (n - Min[##] + 1) (n - Min[##])/2 + Abs[# - #2] &},
   With[{pii = Table[x[[i]] - Sum[If[i == j, 0, pij[[pos[i, j]]]], {j, n}], {i, n}]/2, 
     lim = Sequence @@ Join @@ Table[{pij[[pos[i, j]]], 
       x[[i]] - Sum[If[i == k, 0, pij[[pos[i, k]]]], {k, j - 1}], 
       0, -If[j == n, 2, 1]}, {i, n}, {j, i + 1, n}]},
    With[{arg = Times @@ ((x! σii^pii)/(2 pii)!!) Times @@ (σij^pij/pij!)},
      If[Length[{lim}] == 0, arg, Sum[arg, lim]]]]]]

moment2[{19, 20, 21}] // Hash // AbsoluteTiming
moment[{19, 20, 21}] // Hash // AbsoluteTiming

{0.036750, 4700900427412246901}

{2.762643, 4700900427412246901}

It is very fast and requred no memory for memoization (for large moments moment takes a huge amount of memory). The leaf count is small as in the wolfies' answer.

Derivation of the formula

One variable

At the beginning, let us consider the simplest case with one variable $E(x_1^{r_1})$. The moment generation function is

$$ m(t_1) = \exp\left(\frac{1}{2}\sigma_{11}t_1^2\right). $$

The first derivatives are

$$ \frac{\partial m}{\partial t_1}(t_1) = t_1 \sigma_{11}\exp\left(\frac{1}{2}\sigma_{11}t_1^2\right), $$

$$ \frac{\partial^2 m}{\partial t_1^2}(t_1) = (t_1^2 \sigma_{11}^2+\sigma_{11})\exp\left(\frac{1}{2}\sigma_{11}t_1^2\right), $$

$$ \frac{\partial^3 m}{\partial t_1^3}(t_1) = (t_1^3 \sigma_{11}^3+3t_1\sigma_{11}^2)\exp\left(\frac{1}{2}\sigma_{11}t_1^2\right), $$

$$ \frac{\partial^4 m}{\partial t_1^4}(t_1) = (t_1^4 \sigma_{11}^4+6t_1^2\sigma_{11}^3+3\sigma_{11}^2)\exp\left(\frac{1}{2}\sigma_{11}t_1^2\right) $$

and so on. Then to calculate the moment we need to put $t_1=0$. For the forth moment we have

$$ \frac{\partial^4 m}{\partial t_1^4}(0) = 3\sigma_{11}^2. $$

This process can be represented by the following scheme

We can take the derivative of:

1. The exponent. It increases the power of $t_{1}$ by 1 and multiply by $\sigma_{11}$. It is represented by the upward arrows (all of them has multiplication factor $\sigma_{11}$).

2. The pre-exponential polynomial. It decreases the power of $t_1$ by 1 and multiply by the current power of $t_1$. It is represented by the downward arrows (their multiplication factors correspond to the vertical position).

At the end we need to come to the zero vertical position. All other terms disappear after substitution $t_1=0$. To obtain the moment we need to sum products of factors of all possible paths.

Example for $r_1=10$:

One can check that we obtain the known result

$$ E(x_1^{r_1}) = \left\{\begin{array}{ll} (r_1-1)!!\sigma_{11}^{r_1/2} & \text{if }r_1\text{ is even},\\ 0 & \text{if }r_1\text{ is odd}. \end{array}\right. $$

For the odd $r_1$ there is simply no paths.

Two variables

With two variables the moment generating function is

$$ m(t_1,t_2) = \exp\left(\frac{1}{2}\sigma_{11}t_1^2+\frac{1}{2}\sigma_{22}t_2^2+\sigma_{12}t_1t_2\right). $$

The moment $E(x_1^{r_1}x_2^{r_2})$ can be calculated with

$$ E(x_1^{r_1}x_2^{r_2}) = \frac{\partial^{r_1+r_2} m}{\partial t_1^{r_1}\partial t_2^{r_2}}\Bigg|_{\substack{t_1=0,\\t_2=0}} $$

For definiteness we will:

$~~\rm A.$ take all derivatives with respect to $t_1$,

$~~\rm B.$ then take all derivatives with respect to $t_2$.

At the stage $\rm A$ there are three possibilities:

1. Increase the power of $t_{1}$ by 1 and multiply by $\sigma_{11}$.

2. Increase the power of $t_{2}$ by 1 and multiply by $\sigma_{12}$.

3. Decrease the power of $t_1$ by 1 and multiply by the current power of $t_1$.

After the stage $\rm A$ we substitute $t_1=0$.

At the stage $\rm B$ there are only to possibilities:

1. Increase the power of $t_{2}$ by 1 and multiply by $\sigma_{22}$.

2. Decrease the power of $t_2$ by 1 and multiply by the current power of $t_2$.

We don't consider production of powers of $t_1$ because they will be killed at the final substitution $t_1=0,t_2=0$.

Let us consider $E(x_1^6x_2^4)$. If at the stage $\rm A$ we always choose cases 1 or 3 (not 2) then it can be represented as multiplication of two full diagrams

If at the stage $\rm A$ we choose case 2 two times then it can be represented as multiplication of the diagrams

The first diagram is smaller by 2 because we consume 2 of 6 derivatives with respect to $t_1$ to produce $t_2$. The binomial coefficient $\binom{6}{2}$ is the number of possible choices of the case 2. The second diagram starts from the position 2 because we produce $t_2^2$ at the stage $\rm A$.

One can show that the sum of the diagram with $r_2$ derivatives and the initial position $k_2$ is $$ \left\{\begin{array}{ll} \frac{r_2!}{(r_2-k_2)!!}\sigma_{22}^{(r_2-k_2)/2} & \text{if }r_2-k_2\text{ is even},\\ 0 & \text{if }r_2-k_2\text{ is odd}. \end{array}\right. $$

Therefore, the term $\sigma_{11}^{p_{11}}\sigma_{12}^{p_{12}}\sigma_{22}^{p_{22}}$ in the moment $E(x_1^{r_1}x_2^{r_2})$ has the coefficient

$$ \frac{(2p_{11})!}{(2p_{11})!!}\frac{r_1!}{(p_{12})!(r_1-p_{12})!}\frac{r_2!}{(r_2-p_{12})!!} =\\ \frac{r_1!}{(2p_{11})!!}\frac{1}{p_{12}!}\frac{r_2!}{(2p_{22})!!} $$

where I use that $2p_{11}+p_{12}=r_1$ and $p_{12}+2p_{22}=r_2$. This formula tell us the form of the general formula which I wrote in the beginning. One can check that the formula has the same form for any number of variables. I didn't write it here because it is much more complicated. I just give an example for diagrams in the case of three variables.

Three variables

Let us consider $E(x_1^7x_2^6x_3^5)$ and coefficient before $\sigma_{11}^2\sigma_{12}^\vphantom{2}\sigma_{13}^2\sigma_{22}^2\sigma_{23}^\vphantom{2}\sigma_{33}^\vphantom{2}$.

The coefficient is

4!/4!! Multinomial[1, 2, 4] 5!/4!! Multinomial[1, 5] 5!/2!!

1701000

The general formula returns the same

moment2[{7, 6, 5}]

... + 1701000 σ[1, 1]^2 σ[1, 2] σ[1, 3]^2 σ[2, 2]^2 σ[2, 3] σ[3, 3] + ...

If anybody knows this formula please write where it is published!