[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842

[Pauli Virtanen]

16.02.2012 14:14, josef.pktd at gmail.com kirjoitti:
[clip]
> We had other cases of several patterns in quasi-deterministic linalg
> before, but as far as I remember only in the final digits of
> precision, where it didn't matter much except for reducing test
> precision in my cases.
> 
> In the random multivariate normal case in the ticket the differences
> are large, which makes them pretty unreliable and useless for
> reproducability.

Now that I read your mail more carefully, the following piece of code
indeed does not give reproducible results on Linux with ATLAS either:

--------
import numpy as np
from numpy.linalg import svd

d = 10
alpha = 1 / d**0.5
mu = np.ones(d)
R = alpha * np.ones((d, d)) + (1 - alpha) * np.eye(d)

for i in range(10):
    u, s, vH = svd(R)
    print vH[-1,1], abs(u.dot(np.diag(s)).dot(vH)-R).max()
print s
-----------

Of course, the returned SVD decomposition *is* correct in all cases.

The reason seems to be that the matrix has 9 coinciding singular values,
and the (alignment-dependent) rounding error is sufficient to perturb
the choice (or order?) of singular vectors.

So, the algorithm used to generate multivariate normal random numbers is
then actually numerically unstable, as it relies on the order of
singular vectors returned by SVD.

I'm not sure how to fix this. Maybe the vectors returned by SVD should
be sorted if there are numerically close singular values. Just ensuring
alignment of the input probably won't guarantee reproducibility across
platforms.

Please file a bug ticket, so this doesn't get forgotten...

	Pauli