Needs a test. Would be worth checking why CI is failing too.

@eric-wieser

Because the generators may change, the simplest test might be some rough statistics on the result. What do you suggest?

Re the test: We can just use the singular matrix from my example and then check that the first two samples are approx. equal for svd and eigh. Not sure what other things are tested already, but we should probably at least do some checks on the (marginal) variances and means if that is not done yet.

I do think that we should be consistent about the factor - typically the correlation is done as y = F @ x s.t. E[y y.T] = E[F @ x @ x.T @ F.T] = F @ E[x @ x.T] @ F.T = F @ F.T, so we should have F @ F.T = cov.

Note also that if many evals are zero, we could remove these evecs (or columns of U,V for svd) and work with sth like the following:

def eigh_sample(mean, cov, cutoff=1e-10):
    evals, evecs = np.linalg.eigh(cov)
    pos = evals > cutoff
    evecs_pos = evecs[:, pos]
    base_samples = rg.normal((np.sum(pos), 1))
    scale = np.sqrt(evals[pos])
    return mean + (scale * evecs_pos) @ base_samples

Not sure if it's worth doing this if we re-compute the eigendecomposition/svd every time though.

Not sure if it's worth doing this if we re-compute the eigendecomposition/svd every time though.

There was some discussion of adding the option of passing the left-hand factor, but that actually seems more appropriate to a class based implementation, that is, one could ask for a generator that returned multivariate normals with given covariance. If I needed a bunch of variates, that is how I would do it, but that might not be appropriate for NumPy.

I agree this probably doesn't belong into numpy. I'm surprised scipy's mvn for "frozen" rvs doesn't deal with this, looks like the matrix decomposition is recomputed for each sample here: https://github.com/scipy/scipy/blob/v1.4.1/scipy/stats/_multivariate.py#L660

Torch's MVN does store the factor (though it doesn't allow singular matrices): https://github.com/pytorch/pytorch/blob/master/torch/distributions/multivariate_normal.py#L146-L15

I'm surprised scipy's mvn for "frozen" rvs doesn't deal with this

Might be worth opening a SciPy issue.

but we should probably at least do some checks on the (marginal) variances and means if that is not done yet.

That is what I was thinking. It shouldn't require too large a sample to get sufficient accuracy to catch big errors like the current one. I would suggest starting with an L and passing L @ L.T

EDIT: Or even just generating a stream with L, then resetting the generator state and passing L @ L.T. The resulting streams should be pretty similar and it will avoid problems with changing generators.

Might be worth opening a SciPy issue.

scipy/scipy#11772

Changed to transpose on trailing dims. Added generator-agnostic test for singular covariances. Maybe someone else can take on additional tests checking some basic statistics on the generated samples in a separate PR.

I think the test is too special cased. I think that a rough check the statistics all methods would be useful. For instance, looks like cholesky also yields wrong results.

In [36]: b = default_rng(123).multivariate_normal([0, 0], cov, size=(1000,), method='svd')                                              

In [37]: (b.T @ b)/1000                                                                                                                 
Out[37]: 
array([[1.06677611, 1.05240103],
       [1.05240103, 2.04555217]])

In [38]: b = default_rng(123).multivariate_normal([0, 0], cov, size=(1000,), method='eigh')                                             

In [39]: (b.T @ b)/1000                                                                                                                 
Out[39]: 
array([[ 0.38497443, -0.0150749 ],
       [-0.0150749 ,  2.79194006]])

In [40]: b = default_rng(123).multivariate_normal([0, 0], cov, size=(1000,), method='cholesky')                                         

In [41]: (b.T @ b)/1000                                                                                                                 
Out[41]: 
array([[2.14015452, 1.06968605],
       [1.06968605, 1.03675997]])

In [42]: cov                                                                                                                            
Out[42]: 
array([[1., 1.],
       [1., 2.]])

In [43]: numpy.__version__                                                                                                              
Out[43]: '1.19.0.dev0+f1a247f'

Oh so the Cholesky samples are broken too. We may as well just go all in on the _factor @ x representation then and only transpose the svd factor

We will need a release note to go with this. You can see how to do this by looking at the current release note files in doc/release/upcoming_changes/.

Changed representation internally. Was easiest to use einsum for the output, hopefully that's ok. Added a test that is pretty loose and should only catch significant errors. At a later stage one could do some proper hypothesis testing with an extremely low significance level to avoid flaky tests.
Added upcoming_changes.

Should the tests cover the new numpy.random.default_rng as well, or are there already similar tests elsewhere?

not sure about other rngs, maybe someone who knows the codebase better can clarify?

Using x @ _factor.T, splitting out the basic stats test & relaxing tolerances.

Looks good, let's see if the tests pass.

Thanks @Balandat. It is very nice to have those errors fixed and tested.

@charris It would be a good idea. @WarrenWeckesser discovered a casting issue using some Monte Carlo methods IIRC. It would make sense to start with the things that are "newer", e.g., the Ziggurat exponentials, gamma, normals, and the bounded integers. These are all widely used in other distributions.