Allowing to use covariance matrices with no Cholesky decomposition

[fjaviersanchez]

The current code breaks if the covariance matrix Cholesky decomposition fails. However, there are cases where np.linalg.inv works, but the Cholesky decomposition fails (and I've found several instances where covariance matrices are like that).

It is true that if the Cholesky decomposition fails, it is likely that the covariance is near singular or the condition number is very high, and any results using those matrices will not be very trustworthy. The current approach places a nice safeguard against very naïve usage. However, I think that for some debugging or advanced usage it would be convenient to just raise a warning, and allow the log-likelihood to be computed as -0.5*np.einsum('i, ij, j', data-theory, inv_cov, data-theory) instead of the current np.dot(x, x) using the Cholesky decomposition.

I have to say that this feature is not 100% necessary for any of the projects that I am currently involved with, but would make things a bit easier (Some covariance matrices that I've used needed additional tlc in order for them to work with firecrown).

Anyway, thoughts? It would be great to hear from people that have run many different chains with different covariance matrices.

[JaimeRZP]

Hi All!

I am in favor of continuing to raise the exception and put the burden on the user to regularize their covariance if it is singular.
Otherwise, users might not even realize their covariance is weird.

[marcpaterno]

I think the arguments in favor of using Cholesky decomposition are:

1. it is fast, because it takes advantage of the expectation that the matrix is symmetric positive-definite, and
2. it fails if the matrix is not symmetric positive definite, thus providing a warning to the user that the covariance matrix is "defective".

Because this is an operation done only once by Firecrown, I do not think the speed argument is very important. But I agree with @JaimeRZP that the exception thrown makes it clear that the covariance matrix needs to be fixed.

A alternative would be to use SVD as a fallback, printing a warning that the covariance matrix is suspect. It is quite easy to ignore a warning in an often large output from a real program. The risk of users not noticing the warning, and thus proceeding to do questionable analysis with an covariance matrix that is suspect but not noticed to be suspect is considerable.

[fjaviersanchez]

Yes, I agree with your points above, the risk of not noticing the warnings is high. I am happy to continue with the current approach. Please, feel free to close this discussion :).

[elisachisari]

I would prefer the exception. But what is the different technique you have in mind? Could this SVD option in section 4 of https://arxiv.org/pdf/astro-ph/9912149 be a fall back? Does anyone have thoughts on using this?

[marcpaterno]

Here was my argument of the dangerous feature of having the poorly-conditioned covariant matrix yield a result (rather than an exception):

A alternative would be to use SVD as a fallback, printing a warning that the covariance matrix is suspect. It is quite easy to ignore a warning in an often large output from a real program. The risk of users not noticing the warning, and thus proceeding to do questionable analysis with an covariance matrix that is suspect but not noticed to be suspect is considerable.

[tilmantroester]

Rather than implement a fallback within firecrown, another option might be to allow providing the covariance in other forms (such as providing the inverse directly). That way it's up to the user to deal with ill-conditioned covariances but if they somehow can invert it, firecrown can just use the inverse directly without having to implement these more advanced methods itself.

[marcpaterno]

@bclevine is working on a new likelihood class that saves an inverse covariance matrix, which is different from what ConstGaussian does.

@tilmantroester would this satisfy the use you have in mind?