1

I am working on a homework assignment and I noticed that when the dimension of mean and covariance is very high, multivariate_normal will occupy all CPU forever, without generating any results.

An example code snippet,

cov_true  = np.eye(p)
mean_true = np.zeros(p)
beta_true = multivariate_normal(mean_true, cov_true, size=1).T

when p=5000, this will run forever. environment, python3.4 and python3.5, numpy 1.11.0

Is it really a bug or did I miss something?

13
  • It works for me. Same versions. change the 3rd line to this see if it works: beta_true = np.random.multivariate_normal(mean_true, cov_true, size=1).T Oct 18, 2016 at 7:20
  • yes, just different import, how long does it take for you to run this line?
    – 1a1a11a
    Oct 18, 2016 at 22:27
  • Does "yes" means it worked and didn't take 100% cpu? mine: --- 15.3049829006 seconds --- Oct 19, 2016 at 10:42
  • Oh, I used a different import, from np.random import multivariate_normal, yes, I mean using 100% CPU, on my 13' MacBook pro, real 0m53.319s, user 1m40.845s, sys 0m2.128s, and on a modern workstation, it is slightly better, but it uses all 48 cores, I can't understand why. @Yugi
    – 1a1a11a
    Oct 19, 2016 at 19:01
  • I think there must be something wrong if your time is 15 seconds, my time is 50 seconds, are you using p=5000? and my test program is just import, p=5000 and these three lines.
    – 1a1a11a
    Oct 19, 2016 at 19:04

1 Answer 1

4

What takes so much time?

To account for relations between variables NumPy computes the singular value decomposition of your covariance matrix and this takes the majority of the time (the underlying GESDD is in general Θ(n3), and 50003 is already a bit).

How can things be sped up?

In the simplest case with all variables independent, you could just use random.normal:

from numpy.random import normal

sample = normal(means, deviations, len(means))

Otherwise, if your covariance matrix happens to be full rank (hence positive-definite), supplant svd with cholesky (still Θ(n3) in general, but with a smaller constant):

from numpy.random import standard_normal
from scipy.linalg import cholesky

l = cholesky(covariances, check_finite=False, overwrite_a=True)
sample = means + l.dot(standard_normal(len(means)))

If the matrix may be singular (as is sometimes the case), then either wrap SPSTRF or consider helping with scipy#6202.

Cholesky will likely be noticeably faster, but if that's not sufficient, then further you could research if if it wouldn't be possible to decompose the matrix analytically, or try using a different base library (such as ACML, MKL, or cuSOLVER).

3
  • It reports error, sample = normal(mean_true, cov_true, len(mean_true)) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "mtrand.pyx", line 1495, in mtrand.RandomState.normal (numpy/random/mtrand/mtrand.c:10068) ValueError: scale <= 0
    – 1a1a11a
    Oct 23, 2016 at 1:23
  • You can use something like np.sqrt(cov_true.diagonal()) to pick standard deviations from the covariance matrix. The first method is only correct if all off-diagonal entries are zero.
    – wrwrwr
    Oct 23, 2016 at 13:54
  • @wrwrwr does that mean that the complexity is constant in the number of samples, i.e. np.random.multivariate_normal(size=100) isn't too much slower than np.random.multivariate_normal(size=10)? Or is it linear in the size? May 1, 2017 at 19:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.