[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842
Thu Feb 16 07:54:50 CST 2012
On Thu, Feb 16, 2012 at 8:45 AM, Pauli Virtanen <email@example.com> wrote:
> 16.02.2012 14:14, firstname.lastname@example.org 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()
> 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
> Please file a bug ticket, so this doesn't get forgotten...
the multivariate normal case is already
I can add the diagnosis.
If I interpret you correctly, this should be a svd ticket, or an svd
ticket as "duplicate" ?
> NumPy-Discussion mailing list
More information about the NumPy-Discussion