[Numpy-discussion] strange behavior of numpy.random.multivariate_normal, ticket:1842

[josef.pktd at gmail.com]

On Thu, Feb 16, 2012 at 9:08 AM, Pauli Virtanen <pav at iki.fi> wrote:
> 16.02.2012 14:54, josef.pktd at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion