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

[josef.pktd at gmail.com]

On Wed, Feb 15, 2012 at 10:52 PM,  <josef.pktd at gmail.com> wrote:
> Doing a bit of browsing in the numpy tracker, I found this. From my
> search this was not discussed on the mailing list.
>
> http://projects.scipy.org/numpy/ticket/1842
>
> The multivariate normal random sample is not always the same, even
> though a seed is specified.
>
> It seems to alternate randomly between two patterns, the sum of the
> random numbers is always the same.
>
> Is there a randomness in the svd?  Given the constant sum, I guess the
> underlying random numbers are the same.
>
> (windows 7, python 32bit on 64 machine, numpy 1.5.1)
>
> import numpy
>
> d = 10
> alpha = 1 / d**0.5
> mu = numpy.ones(d)
> R = alpha * numpy.ones((d, d)) + (1 - alpha) * numpy.eye(d)
> rs = numpy.random.RandomState(587482)
> rv1 = rs.multivariate_normal(mu, R, 1)
> print rv1, rv1.sum()
>
> numpy.random.seed(587482)
> rv2 = numpy.random.multivariate_normal(mu, R, 1)
> print rv2[0][0], rv2.sum()
>
>
> running it a few times
>
>>>>
> [[ 0.43028555  1.06584226 -0.03496681 -0.31372591 -0.49716804 -1.50641838
>   0.99209124  0.57236839 -0.32107663  0.5865379 ]] 0.973769580498
> 0.0979690286727 0.973769580498
>>>>
> [[ 0.09796903  1.41010513 -1.10250773  0.71321445  0.09903517  0.36432555
>  -1.27590062  0.04533834  0.37426153  0.24792873]] 0.973769580498
> 0.430285553017 0.973769580498
>>>>
> [[ 0.09796903  1.41010513 -1.10250773  0.71321445  0.09903517  0.36432555
>  -1.27590062  0.04533834  0.37426153  0.24792873]] 0.973769580498
> 0.430285553017 0.973769580498
>>>>
> [[ 0.43028555  1.06584226 -0.03496681 -0.31372591 -0.49716804 -1.50641838
>   0.99209124  0.57236839 -0.32107663  0.5865379 ]] 0.973769580498
> 0.0979690286727 0.973769580498
>
> Josef


numpy linalg.svd doesn't produce always the same results

running this gives two different answers,
using scipy.linalg.svd I always get the same answer, which is one of
the numpy answers
(numpy random.multivariate_normal is collateral damage)

What I don't understand is that numpy.random uses numpy.dual.svd which
I thought is scipy.linalg if available, but it looks like it takes the
numpy svd.

--------
import numpy as np
#from numpy.dual import svd
from numpy.linalg import svd
#from scipy.linalg import svd

d = 10
alpha = 1 / d**0.5
mu = numpy.ones(d)
R = alpha * numpy.ones((d, d)) + (1 - alpha) * numpy.eye(d)

for i in range(10):
    (u,s,v) = svd(R)
    print 'v[-1]', v[-1]
-----------

Josef