[NumPy-Tickets] [NumPy] #1842: inconsistency in RandomState
NumPy Trac
numpy-tickets@scipy....
Thu Feb 16 09:01:11 CST 2012
#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
More information about the NumPy-Tickets
mailing list