celerite.solver.LinAlgError: failed to factorize or solve matrix about celerite HOT 15 CLOSED

Thanks for you help! I will try new alternatives as it seems to be some "internal" issue with the algebra.
I thought that tau was |x -x'|^2, but it is not.

from celerite.

I'm not sure that "internal" issue is quite the right take home here. You'll find the same thing with every other linear algebra system (including industry standards like LAPACK)!

from celerite.

This means that your parameters went somewhere crazy and you're getting an invalid kernel. You should catch this error and return -np.inf from your log probability function. Better yet would be to choose prior bounds on the parameters that don't include crazy values :) but that can sometimes be tricky.

from celerite.

@dfm @huhell There is a discussion of choosing a positive definite kernel in Appendix A of the celerite paper, which amounts to making sure that the power spectrum (which is the Fourier transform of the kernel) does not go negative. Sturm's theorem can be utilized to accomplish this, and as Dan mentions, if this is not satisfied, the likelihood can be set to zero.

from celerite.

That's not what is happening here - Sturm's theorem is already being applied. You can still get numerically unstable solves when the kernel is formally positive.

from celerite.

I have well defined priors on my parameters, and my log probability function returns -np.inf in case of issues. The only parameters without priors are the kernel parameters, but this is very hard to prioritize to my opinion, no ?

from celerite.

It depends on the kernel and what you're trying to do with it. But, you should always put priors on all your parameters. For example, it probably doesn't make sense for one of the hyperparameters to be 1e15 (or if it does, you should change units so that it doesn't). If you put reasonable prior bounds on the parameters then this probably won't happen (especially if you use a shoterm instead of a complexterm directly) but you can still do as I suggested and just catch the exception and return -np.inf because that means that it's an invalid model.

from celerite.

Alright, thank you very much !

from celerite.

Hi,
I have been trying celerite and I encountered the same issue mentioned above (LinAlgError: failed to factorize or solve matrix). My log-likelihood function does return -np.inf with invalid inputs and I also added bounds to the GP hyperparameters (e.g., terms.ComplexTerm(..., bounds=bounds)). The weird thing is that the same piece of code works with some data (fitting a lightcurve), but not with another set of values. I am not sure what else I could try. I am normalising my values. I have tried similar things with george and it works, but not with celerite.
Thanks!

from celerite.

@temuller The power spectrum needs to be positive everywhere so that the matrix can be factorized; this corresponds to having a positive definite covariance matrix, which is needed to carry out Cholesky decomposition. What type of kernel are you creating? If you allow the term to have arbitrary complex coefficients, then positive definiteness is not guaranteed. However, you can check the coefficients beforehand using Sturm's theorem, as described in the appendix A of the celerite paper. This is implemented through the function check_parameters() described here: https://celerite.readthedocs.io/en/stable/python/kernel/?highlight=sturm If you are already using check_parameters(), then I would be interested to know more about their values.

from celerite.

I am creating a simple kernel with log_a=log_b and log_c=log_d. The bounds are such that a,b,c and d are all positive. I will use check_parameters() and see what I get. I will check the appendix in the paper as well. I will post an update, thanks!

from celerite.

If I use check_parameters() I get True which means that the parameters are valid. I noticed something weird in the process. If I add bounds to my mean function, when I call gp.get_parameter_bounds() these bounds do not appear, but they do appear inside terms.ComplexTerm(), why is that?

from celerite.

@temuller I'm not sure about the parameter bounds; @dfm can address that.

Sorry, I should have looked at the discussion above your comment before mentioning Sturm's theorem again!

Setting a=b and c=d leads to a power spectrum of the form S(\omega) = \sqrt{2/\pi} 4ac^3/(\omega^4 + 4c^4). So, yes, this ought to be positive definite as long as ac > 0. What values of a and c are you evaluating this for?

from celerite.

@ericagol so, for a I am using 0.02 and for c I am using 0.67.

I have a related question as well. I have tried using terms.RealTerm() with ac > 0, which should give me an squared-exponential kernel. This works, but the resulting fit does not look like an squared-exponential kernel (i.e., it does not result in a smooth fit), it looks more like straight lines in between data points (this depends on the value of c).

from celerite.

@temuller: There's a lot going on here :)

1. Poorly conditioned matrices can still be numerically negative even if they are formally positive definite. That's what you're seeing when you get the LinAlgError. As I mentioned about, you should just catch that exception and label that as an invalid model. Alternatively, you can pass quiet=True to the log_likelihood function to have it return -np.inf instead of throwing the error.

2. A RealTerm will never be similar to an exponential squared: exp(-|dt|) vs. exp(-|dt|^2). I'm not sure what you mean by "straight line between data points", but I generally don't recommend using a plain RealTerm. It's always better to use a SHOTerm because of its high-frequency behavior. If you believe that there is a bug, please open a separate issue since this is unrelated.

3. I also don't understand your question about bounds. If there is a bug or something that isn't clear, please open a new issue with sample code that demonstrates the problem.

from celerite.