Moment calculations seem incorrect

[arthur-lin1027]

The compute_moments_inefficient_implementation function should return consistent results with increasing maxdeg argument. However, it does not. Below I show a case with a multivariate gaussian centered at (1,0,0) with unit covariances.

Example:

A = np.eye(3,3)
a = np.array([1, 0, 0])
maxdeg = 4
moments = compute_moments_inefficient_implementation(A, a, maxdeg)
np.round(moments,4)

Output:

array([[[ 15.7496,   0.    ,  15.7496,   0.    ,  47.2488],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [ 15.7496,   0.    ,  15.7496,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [ 47.2488,   0.    ,   0.    ,   0.    ,   0.    ]],

       [[ 15.7496,   0.    ,  15.7496,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [ 15.7496,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ]],

       [[ 31.4992,   0.    ,  31.4992,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [ 31.4992,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ]],

       [[ 62.9984,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ]],

       [[157.4961,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ],
        [  0.    ,   0.    ,   0.    ,   0.    ,   0.    ]]])

If I repeat this same code with maxdeg=3, I obtain the following results.

A = np.eye(3,3)
a = np.array([1, 0, 0])
maxdeg = 3
moments_small = compute_moments_inefficient_implementation(A, a, maxdeg)
np.round(moments_small,4)

Output:

array([[[15.7496,  0.    , 15.7496,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [15.7496,  0.    ,  0    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ]],

       [[15.7496,  0.    , 15.7496,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [15.7496,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ]],

       [[31.4992,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ]],

       [[62.9984,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ],
        [ 0.    ,  0.    ,  0.    ,  0.    ]]])

(moments[:4, :4, :4] == moments_small).all()

Output:

np._False

You can see that the subset of moments don't match with moments_small, i.e. moments[:4, :4, :4] != moments_small, for example, moments[0,2,2] = 15.7496 and moments_small[0,2,2] = 0. They should be the exact same value.

In general, these moments matrices are looking way too sparse, for example, moments[2,:,:] being all zeros except for one value doesn't make sense. For example, moments[2,2,2] should definitely be a nonzero value. For the above case with $A=I^{3\times 3}$, $r_{ij} = [1,0,0]^T$, we should have moments[2,2,2] be:

$$ \begin{align*} &\int d^3r e^{-\frac{1}{2}(\mathbf{r-r_{ij}})\mathbf{A}(\mathbf{r-r_{ij}})} x^2 y^2 z^2 \\ =&\int dx e^{-\frac{1}{2}(x-1)^2} x^2 \int dy e^{-\frac{1}{2}y^2} y^2 \int dz e^{-\frac{1}{2}z^2} z^2 = (2*\sqrt{2\pi}) (\sqrt{2\pi}) (\sqrt{2\pi}) = 31.4992 \neq 0 \end{align*} $$

There is likely an issue with the way we are iterating that is failing to update some of these values properly.

[arthur-lin1027]

Ok... upon further thought, I understand that the intention of this is that moments[2,2,2] actually has degree 6, so it never gets calculated if we input a maxdeg < 6. However, I am curious to see if the zeros in the [2,2,2] position are ever used when we use a lower maximum degree, and if that is a correct, or an incorrect calculation.