[NumPy-Tickets] [NumPy] #1842: inconsistency in RandomState

[NumPy Trac]

#1842: inconsistency in RandomState
--------------------------+-------------------------------------------------
 Reporter:  gtakacs       |       Owner:  somebody   
     Type:  defect        |      Status:  new        
 Priority:  normal        |   Milestone:  Unscheduled
Component:  numpy.random  |     Version:  1.6.0      
 Keywords:                |  
--------------------------+-------------------------------------------------

Comment(by josefpktd):

 see discussion at http://mail.scipy.org/pipermail/numpy-
 discussion/2012-February/060575.html

 multivariate_normal uses svd to transform the normal random variables

 In the case of multiplicity of singular values, as in this case, svd
 doesn't guarantee a deterministic order, with ATLAS and MingW we get
 different results, with MKL it looks we always get the same results.
 scipy.linalg.svd also produced identical (deterministic) results in the
 runs I tried.

 Pauli's variation of the example

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

 }}}

 pauli's conclusion :

 "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.

 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.
 (end quote)

 It's not clear to me how to fix this, v is not just a reordering, but a
 different decomposition of the same subspace?

-- 
Ticket URL: <http://projects.scipy.org/numpy/ticket/1842#comment:1>
NumPy <http://projects.scipy.org/numpy>
My example project