[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842
Thu Feb 16 07:45:39 CST 2012
16.02.2012 14:14, email@example.com kirjoitti:
> 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
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()
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
Please file a bug ticket, so this doesn't get forgotten...
More information about the NumPy-Discussion