A Linear Time and Constant Space Algorithm to Compute the Mixed Moments of the Multivariate Normal Distributions

Shalosh B. EKHAD and Doron ZEILBERGER

Abstract: Using recurrences gotten from the Apagodu-Zeilberger Multivariate Almkvist-Zeilberger algorithm, we present a linear-time and constant-space algorithm to compute the general mixed moments of the k-variate general normal distribution, with any covariance matrix, for any specific k. Besides their obvious importance in statistics, they are also very significant in enumerative combinatorics, since, when the entries of the covariance matrix remain symbolic, they enable us to count in how many ways, in a species with k different genders, a bunch of individuals can all get married, keeping track of the different kinds of the k(k − 1)/2 possible heterosexual marriages, and the k possible same-sex marriages. We completely implement our algorithm (with an accompanying Maple package, MVNM.txt) for the bivariate and trivariate cases (and hence taking care of our own 2-sex society and a putative 3-sex society), but alas, the actual recurrences for larger k took too long for us to compute. We leave them as computational challenges.

Maple Package

This article is accompanied by a Maple package, MVNM.txt, available from

https://sites.math.rutgers.edu/~zeilberg/tokhniot/MVNM.txt .

The web-page of this article,

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/mvnm.html ,

contains input and output files, referred to in this paper.

The multivariate Normal Distribution

Recall that the probability density function (see [T] and [Wik]) of the multivariate normal distribution with mean 0 and (symmetric) covariance matrix C = (c_{ij})_{1≤i,j≤k} is

f_{C}(x) := \frac{e^{− \frac{1}{2} x^{T} C^{−1}x}}{\sqrt{(2π)^{k} det C}} .

By simple rescaling we can always assume that all the variances are 1, in other words, that the entries of the main diagonal of C are all 1.

We are interested in fast computation of the mixed moments

M_{C}(m_{1}, ⋯, m_{k}) := ∫_{R^{k}} x_{1}^{m_{1}} ⋯ x_{k}^{m_{k}} f_{C}(x_{1}, …, x_{k}) dx_{1} ⋯ dx_{k} .

One way (not a good one!) to compute these moments, for any specific (m_{1}, …, m_{k}) is to diagonalize C, make a change of variables and compute an integral of the form

∫_{R^{k}} ∏_{i=1}^{k} (∑_{j=0}^{k} b_{ij}x_{j})^{m_{i}} e^{− \frac{1}{2} (x_{1}^{2}+...+x_{k}^{2})} dx_{1} ⋯ dx_{k} .

Then expand ∏_{i=1}^{k} (∑_{j=0}^{k} b_{ij}x_{j})^{m_{i}} and use the fact that ∫_{−∞}^{∞} e^{−x^{2}/2}x^{r} dx is 0 if r is odd and \sqrt{2π} · \frac{r!}{2^{r/2}(r/2)!)} if r is even.

A much better way is via the moment generating function ([Wik][T])

∑_{0≤m_{1},...,m_{k}<∞} M_{C}(m_{1}, …, m_{k})\frac{t_{1}^{m_{1}} ⋯ t_{k}^{m_{k}}}{m_{1}! ⋯ m_{k}!} = e^{\frac{1}{2} (∑_{1≤i,j≤k} t_{i}c_{ij}t_{j})} .

This is implemented in procedure MOMd in the Maple package MVNM.txt mentioned above. For example to get the (3, 3, 3, 3)-mixed moment for the generic four-variate normal distribution, with a general (symbolic) covariance matrix

( \begin{matrix}1 & c12 & c13 & c14 \\ c12 & 1 & c23 & c24 \\ c13 & c23 & 1 & c34 \\ c14 & c24 & c34 & 1 \end{matrix} ) ,

enter

lprint(MOMd([[1,c12,c13,c14],[c12,1,c23,c24],[c13,c23,1,c34],[c14,c24,c34,1]],[3,3,3,3]));

Defining M := max(m_{1}, …, m_{k}), this requires O(M^{k}) time and memory.

Another way is to to use the fact that

∫_{R^{k}} \frac{∂}{∂x_{1}} (x_{1}^{m_{1}} ⋯ x_{k}^{m_{k}} f_{C}(x_{1}, …, x_{k})) dx_{1} ⋯ dx_{k} = 0 .

Using the product and the chain rule, and expanding, one gets a certain mixed recurrence, requiring to compute all the (up to) m_{1} ⋯ m_{k} ‘previous’ values, requiring, again O(M^{k}) memory and time.

But thanks to the Apagodu-Zeilberger [ApZ] multivariate extension of the Almkvist-Zeilberger

[AlZ] algorithm there exist pure recurrences, with polynomial coefficients in m_{1}, …, m_{k}, in each of the discrete coordinate directions. The ones for k = 2 are fairly simple (they are essentially secondorder), but the ones for k = 3 are already very complicated. But once found (and we did find them!) this enables a linear-time and constant-space algorithm for computing any (m_{1}, m_{2}, m_{3})- mixed moment. The recurrences are too complicated to be typeset here, but can be read from the Maple source-code of procedure MOM3 in our Maple package.

The syntax is

MOM3(c12,c13,c23,[m1,m2,m3]); .

For example to get the (10, 10, 10) mixed moment as a polynomial in the symbols c12, c13, c23, type:

MOM3(c12,c13,c23,[10,10,10]); .

This should (and does!) give the same answer as

MOMd([[1,c12,c13],[c12,1,c23],[c13,c23,1]],[10,10,10]); .

To really appreciate the superiority of our algorithm, using MOM3, over the straightforward MOM3d, try, for example

restart:  read `MVNM.txt`:  t0:=time():lu1:=MOM3(1/2,1/3,1/4,[570,560,750]); time()-
t0; ,

that would give you the very complicated lu1 in 2.56 seconds, while

t0:=time():  lu2:=MOMd([[1,1/2,1/3],[1/2,1,1/4],[1/3,1/4,1]],[570,560,750]); time()-
t0; ,

would confirm that lu2 and lu1 are the same (good check!), but it takes 631.007 seconds. Warning: Don’t even try to use floating-point! You would get garbage, due to the complexity of the calculations that accumulate the round-off errors. Both ways would give you erroneous answers unless Digits is set very high. If you keep c12, c12, c23 symbolic, the superiority of MOM3 over MOMd is even more apparent.

restart:  read `MVNM.txt`:  time(MOM3(c12,c13,c23,[100,50,40])); , is less than 12

seconds, while doing the same things with MOMd takes 100 times longer! Why this is also Important in Enumerative Combinatorics? Using what Herb Wilf [Wil] used to call generatiningfunctionlogy it is easy to see that, when the entries of the covariance matrix C are kept symbolic, then for the bivariate case, the coefficient of c^{r} in M_{[[1,c],[c,1]]}(m1, m2) is the exact number of ways that m1 men and m2 women can all get married and there are exactly r heterosexual marriages. The coefficiet of

c_{12}^{a_{12}} c_{13}^{a_{23}} c_{23}^{a_{23}} ,

in the polynomial

M_{[[1,c_{12},c_{13}],[c_{12},1,c_{23}],[c_{13},c_{23},1]]}(m_{1}, m_{2}, m_{3}) ,

is the exact number of ways that, in a 3-gender society, with genders S_{1}, S_{2}, S_{3}, that m_{1} individuals of gender S_{1}, m_{2} individuals of gender S_{2}, and m_{3} individuals of gender S_{3}, can all get married (note that you need their total number, m_{1} + m_{2} + m_{3} to be even, or else it is not possible) where there were exactly a_{12} {S1, S2} marriages, a_{13} {S1, S3} marriages, and a_{23} {S2, S3} marriages. For example, if you want to know the number of ways 300 men and 200 women can get married where there were exactly 100 heterosexual weddings (and hence 150 same-sex marriages), type:

coeff(MOM2(c,[300,200]),c,100); ,

to get a certain 564-digit integer. If you want to know, in a 3-gender society, the exact number of ways that 20 individuals of gender S1, 20 individuals of gender S2, and 20 individuals of gender S3, can get married (so altogether there are 30 weddings) with 9 {S1, S2} weddings, 7 {S1, S3} weddings, and 5 {S2, S3} weddings (and hence 30 − 9 − 7 − 5 = 9 same-sex marriages), type:

coeff(coeff(coeff(MOM3(c12,c13,c23,[20,20,20]),c12,9),c13,7),c23,5); ,

getting, in 0.533 seconds, that the number is:

444975998773143505634352562176000000000 .

Sample Data To see the list of lists of lists of polynomials in c12,c13,c23, let’s call it L, such that

L[m1][m2][m3]

is the (m1,m2,m3)-mixed moment of the trivariate normal distribution with covariance matrix

[[1,c12,c13],[c12,1,c13],[c13,c23,1]] for 1 ≤ m1, m2, m3 ≤ 20 look at the output file
https://sites.math.rutgers.edu/~zeilberg/tokhniot/oMVNM1.txt .

To see the first 35 diagonal mixed moments (i.e. up to the (70, 70, 70) mixed moment), see

https://sites.math.rutgers.edu/~zeilberg/tokhniot/oMVNM2.txt .

Enjoy! The recurrences for four dimensions took too long for us, and we leave them as computational challenges. Perhaps they can be done with Christoph Koutschan’s [K] very powerful Mathematica package? References

[AlZ] Gert Almkvist and Doron Zeilberger, The method of differentiating Under The integral sign, J. Symbolic Computation 10 (1990), 571-591.

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/duis.pdf .

[ApZ] Moa Apagodu and Doron Zeilberger, Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory, Adv. Appl. Math. 37 (2006), 139-152. [Special issue in honor of Amitai Regev]

https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/multiZ.html .

[K] Christoph Koutschan, Advanced applications of the holonomic systems approach, PhD thesis, Research Institute for Symbolic Computation (RISC), Johannes Kepler University, Linz, Austria, 2009.

http://www.koutschan.de/publ/Koutschan09/thesisKoutschan.pdf,
http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/ .

[T] Y.L. Tong, “The multivariate normal distribution. Springer Series in Statistics”. New York: Springer-Verlag, 1990.

[Wik] Wikipedia contributors. Multivariate normal distribution, Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, 5 Feb. 2022. Web. 8 Feb. 2022.

[Wil] Herbert S. Wilf, “generatingfunctionology, Academic Press, 1990. Second Edition: 1994; Third edition : 2005 (CRC Press). Freely downloadable from:

https://www2.math.upenn.edu/~wilf/gfologyLinked2.pdf

Shalosh B. Ekhad and Doron Zeilberger, Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

Email: [ShaloshBEkhad, DoronZeil] at gmail dot com .

Exclusively published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger and arxiv.org

Feb. 20, 2022