[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842
josef.pktd@gmai...
josef.pktd@gmai...
Thu Feb 16 09:20:48 CST 2012
On Thu, Feb 16, 2012 at 9:08 AM, Pauli Virtanen <pav@iki.fi> wrote:
> 16.02.2012 14:54, josef.pktd@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
More information about the NumPy-Discussion
mailing list