Repository navigation
random.multivariate_normal is extremely slow #14193
Description
It seems as though using np.random.multivariate_normal to generate a random vector of a fairly moderate size (1881) is very slow. generating the random variables via cholesky decomposition is much faster.
Reproducing code example:
import numpy as np
mean = np.random.rand(1881)
a = np.random.rand(1881, 1881)
cov = a @ a.T / mean.size
# here I define a custom function to generate the random vector using Cholesky decompostion.
def chol_sample(mean, cov):
return mean + np.linalg.cholesky(cov) @ np.random.standard_normal(mean.size)
Timing the speed in an ipython session gives:
In [7]: %timeit chol_sample(mean, cov)
184 ms +- 4.78 ms per loop (mean +- std. dev. of 7 runs, 10 loops each)
In [8]: %timeit np.random.multivariate_normal(mean, cov)
6.45 s +- 32.3 ms per loop (mean +- std. dev. of 7 runs, 1 loop each
The Cholesky factorization method is 35 times faster. I think this is a serious problem when one needs to call this function repeatedly for an MCMC sampling problem.
Numpy/Python version information:
('1.17.0', '3.6.8 (default, Apr 28 2019, 00:51:53) \n[GCC 8.2.1 20181127]')
Activity
zoj613 commentedon Aug 3, 2019
I looked at the source code for
np.random.multiariate_normaland it seems like the line slowing the generation of random values is line3756numpy/numpy/random/mtrand.pyx
Line 3756 in 71c8a10
Is it possible to replace that line with
numpy.linalg.cholesky(cov)instead since it seems faster?charris commentedon Aug 3, 2019
Yep, but we are stuck with it for the legacy rng because of backwards compatibility. There was extensive discussion of it some years back. Looks like the new rng also uses that method, we should fix that. @bashtage.
charris commentedon Aug 3, 2019
What does the timing give if you use
eigh?bashtage commentedon Aug 3, 2019
I think cholesky is less robust than SVD, in particular if the covariance is indefinite.
zoj613 commentedon Aug 3, 2019
Its faster than the current method, but slower than Cholesky.
Timing shows that it is 2.15 times faster than using svd
charris commentedon Aug 3, 2019
SVD is a bad choice because the result hides negative eigenvalues and the resulting factorization is not symmetric in that case. Using
eighwould be a better choice for robustness, I believe R uses that method, and it allows easy checking for negative eigenvalues and dealing with them, either raising an error or setting them to zero if sufficiently small. Cholesky is probably the fastest and suitable for many applications, could maybe add a keyword selecting the method.bashtage commentedon Aug 3, 2019
zoj613 commentedon Aug 3, 2019
Yeah, I would be much happier even with the
eighmethod as the default since it at least faster than the current one. I think implementing it with a symmetry check and an option to use Cholesky would not be difficult.charris commentedon Aug 3, 2019
Cholesky assumes symmetry also, but non-symmetric matrices are a user error rather than an implementation problem. We try to check that now with the SVD by checking the symmetry of the factorization, but it is a hack (I added that bit).
zoj613 commentedon Aug 3, 2019
I can submit a PR for this if we agree that using
eighand adding ause_choleskyis a reasonable workaround to the speed issue.charris commentedon Aug 3, 2019
I was thinking of a
methodkeyword with possible values{"fast", "default", "legacy"}. With theeighmethod as the default. Not sure of the best way to handle negative eigen values when they occur, maybe set them to zero, and check if the resulting covariance is close to the original. I can see some cases where a reduced form of the covariance might be desirable, but expect that most inputs will be small enough that it doesn't really matter. Note that the original SVD method did not handle negative eigenvalues correctly, nor did it check for symmetry IIRC.bashtage commentedon Aug 3, 2019
zoj613 commentedon Aug 3, 2019
Just a note, SVD has a boolean keyword
hermitianthat makes the factorization more efficient when set to True for symmetric input. Here is the timing difference:I'm not sure if this should be set to True by default inside the function given that a covariance matrix is supposed to symmetric. I aslo noticed that setting
hermitianto True gives comparable speed to theeighmethod. What do you think @bashtage @charrischarris commentedon Aug 3, 2019
@zoj613 It still has the same problem with negative eigenvalues as the svd does, they are hidden away in a difference of sign between left/right eigenvectors. What happens when a non-symmetric input is passed? I suspect
hermitian=Truejust callseighand adjusts the signs of the eigenvectors associated with the negative eigenvalues.zoj613 commentedon Aug 3, 2019
you're right. I looked into the code and it does make a call to
eighand then modify the sign of singular values afterwards.charris commentedon Aug 3, 2019
eighshould be more numerically accurate than the genericsvd. I think Robert's original motivation was to deal with near, or completely singular covariances. In those cases it generally won't matter much if small perturbations are added to the result, it's not much worse than just ignoring the negative eigenvalues in general. It is, however, more difficult to determine that adjustments have been made.I like Cholesky because of its speed, although for generating long streams, that too does not matter much. The ability to reuse a factorization, as suggested by @bashtage would be useful as the factorization would only need to be done once. @zoj613 Does your use case use a fixed covariance in multiple calls?
zoj613 commentedon Aug 3, 2019
@charris my use case requires me to calculate a new covariance matrix in each iteration of the MCMC sampler. The newly calculated covariance matrix in each iteration only varies in the diagonal entries. The off-diagonal entries never change.
bashtage commentedon Aug 4, 2019
@zoj613 Sounds like an option that allows you to directly supply the factored covariance would be far and away the fastest then.
zoj613 commentedon Aug 4, 2019
For my use case, I found the best solution was to use the cholesky module in the
scikit-sparsepackage since my covariances were sparse. The speed ups were incredible. I wrote a function that basically combines the cases where a dense or sparse covariance is used. It can be found here:https://github.com/zoj613/OccuSpytial/blob/36f00459a70995806886f52a6520618d6eb06c3b/occuspytial/utils/utils.py#L116-L151
However for general usage I do think the
multivariate_normalfunction needs improvement. I wrote some code to incorporate the suggestions made above, but I would first need to build the cloned repo so I can test the changes and then figure out the rules to contributing on this repo. That is if nobody intends on working on this already.zoj613 commentedon Aug 4, 2019
the unittests pass but there seems to be 2 failing regressions tests. Can someone please help with dealing with this issue? It seems to be caused by an unhandled
RuntimeErrorregarding a non-positive semidefinite input. I've attached the logfile of that error message.errormsg.log
@charris @bashtage
charris commentedon Aug 4, 2019
The test covariance is all zero, i.e., semidefinite. Cholesky will not handle that.