Use Cholesky decomposition in random.multivariate_normal() - #3938

Is there really any circumstance where using SVD is better? Basically I'm
asking if we really need to add an option. How much faster is it for large
arrays?
On 17 Oct 2013 20:13, "Tal Yarkoni" [email protected] wrote:

In numpy.random.multivariate_normal(), decomposition of the covariance
matrix is currently accomplished using SVD. As discussed in #1821#1821 resolution), this is problematic since if a user passes a
non-positive definite covariance matrix, the covariance of the resulting
random variables may diverge considerably from what's expected, without the
user receiving any warning. Currently the only mention of this issue is at
the end of the docstring, which many users are unlikely to see. Other
packages invariably seem to use a decomposition method that raises an
exception if a non-positive definite matrix is provided, and this seems
like sensible behavior.

This patch changes the method's default behavior to use Cholesky
decomposition, which raises an exception when a non-positive definite
matrix is passed. Unless the size of the requested array is very large,
users are unlikely to perceive any difference in speed, and can always

explicitly select svd with the newly added 'decomposition' argument.

You can merge this Pull Request by running

git pull https://github.com/tyarkoni/numpy master

Or view, comment on, or merge it at:

#3938
Commit Summary

• random.multivariate_normal() now uses cholesky decomposition

File Changes

• M numpy/random/mtrand/mtrand.pyxhttps://github.com/Use Cholesky decomposition in random.multivariate_normal() #3938/files#diff-0(29)

Patch Links:

• https://github.com/numpy/numpy/pull/3938.patch
• https://github.com/numpy/numpy/pull/3938.diff

It isn't. IIRC, there was another discussion of this where @rkern said he used the svd because it gave better results on some corner cases. One problem that arises with an svd when there is an negative eigenvalue is that the decomposition it no longer symmetric. The requirement of non-negative singular values changes the sign of some of the orthogonal matrix entries. I'm in favor of Cholesky myself or maybe a straight eigenvalue decomposition with an error raised for negative eigen values and an option to set them to zero.

Thread at http://thread.gmane.org/gmane.comp.python.numeric.general/48163/focus=48191

I didn't have any particular basis in mind for keeping SVD other than the general concern that there might be some good reason for the precedent. R's mvnorm() doc converges with what you're saying, @charris: "The matrix decomposition is done via eigen; although a Choleski decomposition might be faster, the eigendecomposition is stabler." So maybe svd with explicit error raising is indeed the way to go.

I'd recommend eigh rather than SVD as it exposes negative eigenvalues. Setting the negative eigenvalues to zero then provides the closest L_2 approximation to a non-negative definite matrix.

So skimming that thread and the discussion here, it sounds like the current state is:

• An eigendecomposition is the way to go to get the best combination of numerical stability and safety
• Reproducibility of random sampling requires that we have a deterministic eigendecomposition routine, and we're not sure whether we do.
  ?

I'm pretty sure Cholesky is deterministic. The eigen decompostion always has the problem of being not unique for repeated or close eigenvalues and, as it is iterative, the termination condition may vary. Also the sign of the eigenvectors can be changed which will have an effect. The Cholesky can be indeterminate upto the choice of sign for the square root, but in practice I don't know why anyone would choose the negative square root.

We don't have reproducible eig. Lapack output can vary from platform to platform.

I think there are two issues here: (a) what decomposition to use for maximal stability and reproducibility, and (b) irrespective of approach, whether an exception should be raised when a non-positive definite matrix is passed. It sounds like there isn't a consensus on the first issue yet, but hopefully we can all agree that silently returning unexpected results is a bad idea. Given that Cholesky seems to be preferable to svd (though perhaps not to eigh), is deterministic, and will appropriately raise an exception, I'd suggest switching to Cholesky and doing away with the decomposition option in my commit. This seems like the simplest solution code-wise. If there's a broad consensus in favor of eigh in the future, it'll be easy enough to patch again. Sound good?

I'd guess the Cholesky decomposition could be improved with column and row scaling, and there are some references to that in the literature. LDL^T can be used, but requires pivoting for semidefinite matrices.

If SVD is numerically superior except for the part where it messes up with
non-definite matrices, then another possible compromise is to just check
the eigenvalues explicitly and then use SVD. Relatively expensive but
better than being untrustworthy...

On Thu, Oct 17, 2013 at 11:00 PM, Charles Harris
[email protected]:

I'd guess the Cholesky decomposition could be improved with column and row
scaling, and there are some references to that in the literature. LDL^T can
be used, but requires pivoting for semidefinite matrices.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/3938#issuecomment-26556793
.

The eigh decompostion can be turned into SVD by multiplying the appropriate columns in the left orthogonal factor by -1 for negative eigenvalues and using the absolute value of the eigenvalues as the singular values. The SVD isn't unique (and neither is the eigh).

Okay, right. And Wikipedia at least claims cholesky is unique
(mathematically, not numerically) for positive definite matrices (but not
positive-semidefinite matrices).

There is some value in being able to sample from multivariate normal with a
positive semidefinite covariance matrix, though it may not be critical...

On Thu, Oct 17, 2013 at 11:17 PM, Charles Harris
[email protected]:

The eigh decompostion can be turned into SVD by multiplying the
appropriate columns in the left orthogonal factor by -1 for negative
eigenvalues and using the absolute value of the eigenvalues as the singular
values. The SVD isn't unique (and neither is the eigh).

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/3938#issuecomment-26557894
.

General note, what is the policy here? It seems to me that ideally we don't just change the output of the RNG over versions. Maybe we need to set a flag whether the random state was set externally (for that instance of the RNG), and if so give a FutureWarning if None of the modes is specified and that flag is set? That way most people should not notice anything and get the new behaviour, but those who care about reproducibility are forced to consider it over versions.

The policy is to fix bugs in the distributions, but not to make other improvements that would change the output. For example, we would not change the normal distribution to use the more efficient ziggurat method. I don't think we try to preserve things to the last ULP, especially for distributions like this that have a lot of computation beyond that of the raw random bits, so if the new approach is at least roughly the same as the current one (likely), I think this would fall into "fixing bugs".

That said, these covariance matrices are usually small, so extra linear algebra computations are fairly cheap. Can we just check for indefiniteness via eigvals() and continue using the SVD to compute the matrix square root? The initial implementation did use the Cholesky decomposition until I ran into situations where we had semidefinite covariance matrices, and the Cholesky decomposition failed.

As an alternate, we could use eigh with a new keyword cutoff, where eigenvalues whose absolute value is < cutoff are set to 0 and remaining negative eigenvalues raise an error. Numerically, that should be equivalent to SVD but a bit faster and it would handle semidefinite matrices more gracefully. A default value of cutoff could be something like 1e-14, or we could make it dependent on matrix size. Scaling the matrix as d_C_d, where d is 1/sqrt(C.diagonal()) would probably help, although that is more a gut feeling than something I have experience with.

Cholesky with complete pivoting would is also a possibility, but I don't think we have that available, although it would be easy to make one. See http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.139.8064

Actually, any of these methods will produce a drastically different output stream (because the untransformed standard normal variates are isotropic there is a rotational freedom that can be folded into the matrix square root of the covariance). I support raising an error for indefinite inputs, but I think it should be done as a check in addition to the current algorithm. Users who are currently passing a correct covariance matrix should continue to get the same results in the future.

You can check for negative eigenvalues using the SVD

In [14]: a = eye(2)

In [15]: a[0, 1] = a[1, 0] = 1 + 1e-10

In [16]: eigvals(a)
Out[16]: array([  2.00000000e+00,  -1.00000008e-10])

In [17]: dot(v, u)
Out[17]: 
array([[  1.00000000e+00,   2.22044605e-16],
       [  2.22044605e-16,  -1.00000000e+00]])

Probably makes a good check for symmetry also; the off diagonal elements of the product should be small.

Or somewhat faster, O(n^2)

In [4]: sum(v*u, 1)
Out[4]: array([ 1., -1.])

Reusing the SVD seems like the best solution, as it would keep the number stream the same. $n^2$ is quite cheap.

Chiming here on the Cholseky with pivot method for semi-definite positive matrices, if you want to support Gaussian distributions with rank-deficient covariance matrices: it is available in Lapack, as discussed in google-deepmind/torch-randomkit#19 (a package that tries to port bits of mtrand.pyx to the Torch library)

I'm glad that it exists, but changing the algorithm would change the data stream. As a policy try to avoid that except in cases of outright bugs.

Fair enough. I'm somewhat lucky to be starting from scratch now in torch, in that regard: it was only added to Lapack in 2008, apparently.

OK, looks like reusing the SVD is the way to go. Remaining question is what to do when a negative eigenvalue is encountered. I'd suggest setting the corresponding singular value to zero and raising a warning. Thoughts?

Closing this and adding a task.