# 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? – Nawaf Bou-Rabee Jan 11 '18 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 '18 at 13:58
• @NawafBou-Rabee One nice reference is arxiv.org/pdf/1310.2559.pdf – Fabrice Pautot Jan 11 '18 at 14:20
• Thanks, but where in that reference do they state or give the Wick-Isserlis theorem or formula? – Nawaf Bou-Rabee Jan 11 '18 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.. – Fabrice Pautot Jan 11 '18 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. – Fabrice Pautot Jan 12 '18 at 8:11