GP implementation of a precomputed covariance matrix

[gjgilbert]

PyMC3 supports multiplying a GP covariance function with a precomputed covariance matrix, as seen with the example below copied from the PyMC3 docs (https://docs.pymc.io/notebooks/GP-MeansAndCovs.html). Is it possible to do something similar directly with exoplanet?

Thanks!

import matplotlib.pyplot as plt
import pymc3 as pm
import theano
import numpy as np

# data vector
X = np.linspace(0, 2, 200)[:,None]

# first evaluate a covariance function into a matrix
period = 0.2
cov_cos = pm.gp.cov.Cosine(1, period)
K_cos = theano.function([], cov_cos(X))()

# now multiply it with a covariance *function*
cov = pm.gp.cov.Matern32(1, 0.5) * K_cos

K = cov(X).eval()

plt.figure(figsize=(14,4))
plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K).random(size=3).T);
plt.title("Samples from the GP prior");
plt.ylabel("y");
plt.xlabel("X");

[dfm]

This is really only supported for products of celerite kernels (in which case, you can just take the product of the Term objects) because the product of a celerite matrix with a general matrix will be a general dense matrix and it won't be possible to get scalable solves. As far as I can tell, there aren't any general results for the inverse and determinant of element-wise products of matrices that will help in this case.

What's your use case?

[dfm]

PS. If you don't need fast solves, everything is easy, but presumably you want to use exoplanet for the fast GPs?

[gjgilbert]

Yes, exactly, I want to take advantage of the fast GPs.

In practice, I might be okay with the slow case. For my specific use here, I've computed the covariance matrix based on the autocorrelation function (ACF) of the out-of-transit flux from a Kepler light curve. I plan to multiply this matrix by a constant term GP. So, I'm not actually fitting a computationally expensive GP; in the future I might want to, which is why I'd prefer to build this on celerite's fast GPs.

It sounds like the most straightforward solution might be for me to use a fixed kernel rather than a covariance matrix. I think it might work if I model the ACF as a product/sum of celerite Term objects. If I'm not varying the inputs to the kernel during a transit fit, even a model with fairly complicated kernel should run quickly. Does this sound reasonable to you? Or did you have something else in mind when you suggested everything is easy if I don't need fast solves?

[gjgilbert]

I should also mention that my covariance matrix is a sparse block diagonal matrix because I'm only concerned with correlations on the timescale of transits.

[dfm]

It sounds like the most straightforward solution might be for me to use a fixed kernel rather than a covariance matrix. I think it might work if I model the ACF as a product/sum of celerite Term objects. If I'm not varying the inputs to the kernel during a transit fit, even a model with fairly complicated kernel should run quickly.

You definitely don't want to fit the ACF - just fit the GP directly to the out of transit points then you can fix the hyperparameters. In general, I find that this doesn't make much of a difference compared to just marginalizing over the hyperparameters simultaneously with the transit fit (especially if the correlations are on a transit timescale), but maybe it'll matter more in your case!

Or did you have something else in mind when you suggested everything is easy if I don't need fast solves?

What I meant by that is just that PyMC3's GP implementation is just using a standard Cholesky factorization to compute the GP so it doesn't care about the details of the kernel! If PyMC3's GP is fast enough then you can just use that.

I should also mention that my covariance matrix is a sparse block diagonal matrix because I'm only concerned with correlations on the timescale of transits.

I think you probably mean "banded" rather than "block diagonal" here. In that case, there are fast things that you can do, but I'm not sure that those are implemented in Theano.

Since this is specific enough to your use case, I'm going to close this and perhaps you can email me directly if you have more questions about this one? I'll set up a mailing list one of these days :)

[gjgilbert]

Thanks, for your help, Dan. I agree - this is a pretty specific case. I'm going to keep working on this and I'll email you with any more questions I have.

[ericagol]

@gjgilbert I've done something similar before in which I computed the ACF of a Kepler light curve, and then fit the ACF with the sum of celerite terms. This was with an older, slower version of celerite, but it worked well for that purpose! https://github.com/ericagol/tyc3559_julia

[gjgilbert]

Thanks, @ericagol. I'll take a look at this.

[dfm]

Let it be known that I will never condone fitting a model to an ACF estimator for any reason. The whole point of celerite is that that isn't necessary anymore!

[ericagol]

Indeed; I did that before Sivaram developed the faster method, and I didn’t think I could optimize the kernel parameters since I was applying it to an entire Kepler light curve, so the likelihood computation was quite slow.

Eric Agol Astronomy Professor University of Washington

On Dec 13, 2019, at 5:38 AM, Dan Foreman-Mackey notifications@github.com wrote:

Let it be known that I will never condone fitting a model to an ACF estimator for any reason. The whole point of celerite is that that isn't necessary anymore!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.