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 6 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 |
|||