random.multivariate_normal is extremely slow

[zoj613]

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]')

[zoj613]

I looked at the source code for np.random.multiariate_normal and it seems like the line slowing the generation of random values is line 3756

numpy/numpy/random/mtrand.pyx

Line 3756 in 71c8a10

(u, s, v) = svd(cov)

Is it possible to replace that line with numpy.linalg.cholesky(cov) instead since it seems faster?

[charris]

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]

What does the timing give if you use eigh?

[bashtage]

I think cholesky is less robust than SVD, in particular if the covariance is indefinite.

[zoj613]

What does the timing give if you use eigh?

Its faster than the current method, but slower than Cholesky.

import numpy as np

def eigh_sample(mean, cov): 
    s, v = np.linalg.eigh(cov) 
    return mean + v * np.sqrt(s) @ np.random.standard_normal(mean.size)

Timing shows that it is 2.15 times faster than using svd

In [30]: %timeit eigh_sample(mean, cov)                                                                                                                       
3 s +- 361 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)

[charris]

SVD is a bad choice because the result hides negative eigenvalues and the resulting factorization is not symmetric in that case. Using eigh would 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]

eigh will hide non symmetric in outs though, so might need a symmetry check.

…

On Sat, Aug 3, 2019, 17:51 Charles Harris ***@***.***> wrote: SVD is a bad choice because the result hides negative eigenvalues and the resulting factorization is not symmetric in that case. Using eigh would 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. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#14193?email_source=notifications&email_token=ABKTSRODONNUEROLXU2VGRDQCWZSJA5CNFSM4IJDMRX2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD3PR52A#issuecomment-517938920>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABKTSRLYFFZKVNGXKZRNDGTQCWZSJANCNFSM4IJDMRXQ> .

[zoj613]

eigh will hide non symmetric inputs though, so might need a symmetry check.

Yeah, I would be much happier even with the eigh method 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]

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]

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).

I can submit a PR for this if we agree that using eigh and adding a use_cholesky is a reasonable workaround to the speed issue.

[charris]

I was thinking of a method keyword with possible values {"fast", "default", "legacy"}. With the eigh method 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]

Another option would be to allow the user to provide the matrix square root themselves.

…

On Sat, Aug 3, 2019, 20:46 Charles Harris ***@***.***> wrote: I was thinking of a method keyword with possible values {"fast", "default", "legacy"}. With the eigh method 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. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#14193?email_source=notifications&email_token=ABKTSRJYXAAZDWZVL2QNS3LQCXOBVA5CNFSM4IJDMRX2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD3PUTWA#issuecomment-517949912>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABKTSRK42HVZUKX7HXCCAWLQCXOBVANCNFSM4IJDMRXQ> .

[zoj613]

Just a note, SVD has a boolean keyword hermitian that makes the factorization more efficient when set to True for symmetric input. Here is the timing difference:

In [119]: %timeit np.linalg.svd(cov, hermitian=True)                                                                                                          
2.67 s +- 41.8 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)

In [120]: %timeit np.linalg.svd(cov)                                                                                                                          
6.73 s +- 224 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)

In [121]: %timeit np.linalg.eigh(cov)                                                                                                                         
2.62 s +- 30.9 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)

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 hermitian to True gives comparable speed to the eigh method. What do you think @bashtage @charris

[charris]

@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=True just calls eigh and adjusts the signs of the eigenvectors associated with the negative eigenvalues.

[zoj613]

@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=True just calls eigh and adjusts the signs of the eigenvectors associated with the negative eigenvalues.

you're right. I looked into the code and it does make a call to eigh and then modify the sign of singular values afterwards.