# Numerical evaluation/approximation of non-central high-order moments of high-dimensional Gaussian measures?

I need to numerically evaluate/approximate non-central high-order moments of high-dimensional Gaussian measures/distributions with given mathematical expectations and covariance matrices. The Gaussian measures have dimension $d$ with say $d>1000$ and the moments have order $k$ with say $k>>10000$.

More precisely, I "only" need to evaluate the ratios of successive moments

$\frac{{\mathbb{E}X_1^{{k_1}}...X_i^{{k_i} + 1}...X_d^{{k_d}}}}{{\mathbb{E}X_1^{{k_1}}...X_i^{{k_i}}...X_d^{{k_d}}}}\quad k = \sum\limits_{i = 1}^d {{k_i}}$

The classical Wick-Isserlis theorem/formula and similar symbolic formulae look helpless because they involve an astronomically large number of terms $k!!$.

As suggested, I tried a Laplace's approximation/expansion but that does not work well: unless I'm mistaken, in standard form, it gives a wrong approximation that does not depend on ${{k_i}}$ and in fully exponential form, we fall on a multivariate system of quadratic equations of dimension $d$ for the maximum of the function...

I'm unable to find a single reference on this problem, that's quite unexpected.

Of course, Monte-Carlo methods is always an option but it is the last one: fast and deterministic algorithms are highly preferred.

So, please, what are my best options in order to numerically evaluate those moments in general? Computational complexity theoretic results welcome too.

Thanks.

PS.: I finally found a nice, ad hoc, deterministic method allowing to evaluate the ratios of successive moments instantaneously. Unfortunately I can't disclose it right now. I really don't understand why I can't find a single paper and why I don't get any answer about such a basic and fundamental problem??????

• Have you tried using the moment generating function for a multivariate Gaussian? Jan 11, 2018 at 13:19
• @NawafBou-Rabee The moment generating function of a Gaussian is a Gaussian. Which leads to Isserlis formula. However one can probably obtain a recursive version (integrating by parts).
– lcv
Jan 11, 2018 at 13:58
• @NawafBou-Rabee One nice reference is arxiv.org/pdf/1310.2559.pdf Jan 11, 2018 at 14:20
• Thanks, but where in that reference do they state or give the Wick-Isserlis theorem or formula? Jan 11, 2018 at 14:27
• @NawafBou-Rabee Another one thphys.uni-heidelberg.de/~amendola/teaching/compstat-hd.pdf, page 28, but you can find plenty of them by yourself. Unfortunately, the Wick-Isserlis formula for moments of order $k$ involve $k!!$ terms. Hence, with $k>>1000$, a different method is required for evaluating my moments. Astonishingly, for the time being, I'm unable to find numerical methods and algorithms for high-order moments of Gaussian r.v.. Jan 11, 2018 at 16:03

It seems you may want to turn the large $k$ property to your advantage and use a Laplace-type asymptotic method with controlled bounds. For fixed $d$ you may find some useful tools in the book "Analytic Combinatorics in Several Variables" by Pemantle and Wilson.
If also $d$ becomes large, then this starts looking like a problem in statistical field theory. There is a method people used for this kind of thing called the Lipatov argument (see for instance this article by Spencer).
• Thanks for your input. $d$ is fixed, $k$ is slightly variable but no asymptotics. For sure, good old Laplace is always an option, but I guess there exist dedicated methods and algorithms for those fundamental integrals. Jan 12, 2018 at 8:11