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:
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}$).
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:
Increase the power of $t_{1}$ by 1 and multiply by $\sigma_{11}$.
Increase the power of $t_{2}$ by 1 and multiply by $\sigma_{12}$.
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:
Increase the power of $t_{2}$ by 1 and multiply by $\sigma_{22}$.
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!
AbsoluteTiming
of0.004
here. It looks like some serious improvements have happened in the new version. This apparently makes my answer to How to deal with complicated gaussian integrals in Mathematica roughly as fast as yours now... $\endgroup$