On Thu, Feb 16, 2012 at 8:45 AM, Pauli Virtanen <p...@iki.fi> wrote: > 16.02.2012 14:14, josef.p...@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 _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
