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).
beta_true = np.random.multivariate_normal(mean_true, cov_true, size=1).T
--- 15.3049829006 seconds ---
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. @Yugip=5000
and these three lines.