Home
Reading
Searching
Subscribe
Sponsors
Statistics
Posting
Contact
Spam
Lists
Links
About
Hosting
Filtering
Features Download
Marketing
Archives
FAQ
Blog
 
Gmane
From: Pauli Virtanen <pav <at> iki.fi>
Subject: Re: strange behavior of numpy.random.multivariate_normal, ticket:1842
Newsgroups: gmane.comp.python.numeric.general
Date: Thursday 16th February 2012 13:45:39 UTC (over 5 years ago)
16.02.2012 14:14, [email protected] 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
 
CD: 2ms