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

josef.pktd@gmai... josef.pktd@gmai...
Thu Feb 16 07:54:50 CST 2012


On Thu, Feb 16, 2012 at 8:45 AM, Pauli Virtanen <pav@iki.fi> wrote:
> 16.02.2012 14:14, josef.pktd@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...

the multivariate normal case is already
http://projects.scipy.org/numpy/ticket/1842
I can add the diagnosis.

If I interpret you correctly, this should be a svd ticket, or an svd
ticket as "duplicate" ?

Thanks,

Josef


>
>        Pauli
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion


More information about the NumPy-Discussion mailing list