On Thu, Feb 16, 2012 at 9:08 AM, Pauli Virtanen <p...@iki.fi> wrote: > 16.02.2012 14:54, josef.p...@gmail.com kirjoitti: > [clip] >> If I interpret you correctly, this should be a svd ticket, or an svd >> ticket as "duplicate" ? > > I think it should be a multivariate normal ticket. > > "Fixing" SVD is in my opinion not sensible: its only guarantee is that A > = U S V^H down to numerical precision and S are sorted. If the algorithm > assumes something extra, it is wrong. This sort of reproducibility > issues affect potentially all code (depends on the compiler and > libraries used), and trying to combat it at the linalg level is IMHO not > our business --- if someone really wants it, they should tell their C > compiler and all libraries to use a reproducible FP model.
I agree, I added the comments to the ticket. > > However, we should ensure the algorithms we provide are stable against > rounding error. In this case, the random number generation is not, so it > should be fixed. storing the last column of v vli = [] for i in range(10): (u,s,v) = svd(R) print('v[:,-1]') print(v[:,-4:]) vli.append(v[:, -1]) >>> np.unique([tuple(vv.tolist()) for vv in vli]) array([[-0.31622777, -0.11785113, 0.08706383, 0.42953906, 0.75736963, -0.31048693, -0.01693654, 0.10328164, -0.04417299, -0.10540926], [-0.31622777, -0.03661979, 0.61237244, -0.15302481, 0.0664198 , 0.11341968, 0.38265194, 0.51112292, -0.10540926, 0.25335061]]) The different v are not just a reordering of each other. If my linear algebra is correct, then the algorithm provides different basis vectors for the subspace with identical singular values. I don't see any way to fix multivariate_normal for this case, except for dropping svd or for random perturbing a covariance matrix with multiplicity of singular values. 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