MultivariateNormal can't handle semi-positive definite covariance matrix

[xieliaing]

What would you like to submit? (put an 'x' inside the bracket that applies)

• question
• bug report
• feature request

Robust Multivariate Normal Sampling with semi-PD cov matrix

generate method in MultivariateNormalDistribution class requires strictly positive definite covariance matrix but in reality, we could facing cases where a semi-PD covariance matrix is present and we still want to do sampling which is 100% valid. Think of a case where you have two columns that are 100% correlated.

Current implementation leverages on Cholesky Decomposition for the following formula:
Cov = L'L
z~N(0, I)
sample = \mu + L'z
Cholesky Decomposition requires strictly PD cov matrix but its variant LDL' decomposition is able to handle semi-PD cov matrix and is considered as a robust alternative.

I can contribute the LDL' decomposition

Note: If you would like to support the development for this feature or resolution of this bug, consider adding a bounty to it later in https://www.bountysource.com/teams/accord-net/issues

[xieliaing]

Alternative to LDL' decomposition, we can also use SVD for a robust implementation of simulation from multivariate normal with semi-PD cov matrix. It is a little more intensive in computation but all tools are readily available.

[cesarsouza]

Hi @xieliaing!

Thanks a lot again for opening an issue!

The framework currently has an implementation for the LDLt decomposition that can be obtained by passing robust = true to the CholeskyDecomposition constructor. However, the MultivariateNormalDistribution class is not leveraging this functionality right now.

Instead, the way the framework currently offers to overcome the issue where a covariance matrix might be non-positive definite is by computing a SVD decomposition instead of a (robust) Cholesky. However, I hadn't researched quite thoroughly yet how to sample values from a Multivariate Normal whose covariance matrix had been decomposed as an SVD, and for this reason, the Generate method is not supporting this case yet. If you know about the proper way to do it, please let me know!

Otherwise, another possibility would be to update the MultivariateNormal to consider different levels of robustness when decomposing the covariance matrix, such that it could be possible to use a normal Cholesky for when covariance matrices are known to be positive-definite, a robust Cholesky for when they are known to be possibly positive semi-definite, and a SVD for the completely degenerate/singular cases (I have to admit I am not completely sure right now whether the LDL' would have been able to handle singular covariance matrices, or even if it could, whether it would be able to be used for sampling in those cases).

If you would like, please let me know of some examples of cases that you think should have been handled by the framework, posting the given inputs and the expected outputs (even if they would result in non-positive definite matrices), and I might be able to take a look on how to wire the Generate method to work in those cases.

Thanks!

Regards,
Cesar

[cesarsouza]

Oh actually I had posted my last answer before seeing your update about the SVD. If you happen to know how to sample values from a Multivariate Normal distribution whose covariance matrix has been decomposed using an SVD instead of a (robust) Cholesky, please let me know - I hadn't have the time to look into it yet!

Thanks!

Regards,
Cesar

[xieliaing]

@cesarsouza Sure thing. I will check it out and post here, may be some python sample code if the approach works.

[xieliaing]

Looks like we can just follow the svd approach used in multivariate_normal method from Numpy:
https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx

They modified multivariate_normal method specifically for the same reason listed here:
numpy/numpy#5726

[xieliaing]

@cesarsouza pushed an implementation using SVD.

[cesarsouza]

Hi @xieliaing,

I've manually created a pull request from your fork: #1086

Regards,
Cesar