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??????