Fast and scalable Gaussian process modeling with applications to astronomical time series

TL;DR

This paper introduces a new linear-scaling Gaussian Process method for one-dimensional data, enabling efficient probabilistic analysis of large astronomical time series without requiring evenly spaced observations.

Contribution

The paper presents a novel Gaussian Process modeling approach with linear computational complexity for 1D data, applicable to large astronomical datasets and based on covariance functions as mixtures of complex exponentials.

Findings

Method scales linearly with dataset size

Successfully applied to simulated and real astronomical data

Provides interpretable probabilistic inference

Abstract

The growing field of large-scale time domain astronomy requires methods for probabilistic data analysis that are computationally tractable, even with large datasets. Gaussian Processes are a popular class of models used for this purpose but, since the computational cost scales, in general, as the cube of the number of data points, their application has been limited to small datasets. In this paper, we present a novel method for Gaussian Process modeling in one-dimension where the computational requirements scale linearly with the size of the dataset. We demonstrate the method by applying it to simulated and real astronomical time series datasets. These demonstrations are examples of probabilistic inference of stellar rotation periods, asteroseismic oscillation spectra, and transiting planet parameters. The method exploits structure in the problem when the covariance function is expressed as a mixture of complex exponentials, without requiring evenly spaced observations or uniform noise. This form of covariance arises naturally when the process is a mixture of stochastically-driven damped harmonic oscillators -- providing a physical motivation for and interpretation of this choice -- but we also demonstrate that it can be a useful effective model in some other cases. We present a mathematical description of the method and compare it to existing scalable Gaussian Process methods. The method is fast and interpretable, with a range of potential applications within astronomical data analysis and beyond. We provide well-tested and documented open-source implementations of this method in C++, Python, and Julia.

Tables

Table 1. Table 1: The parameters and priors for Example 1.

parameter	prior
$\ln \left(S_0\right)$	$𝒰(-10,\mathrm{\hspace{0.17em}10})$
$\ln \left(Q\right)$	$𝒰(-10,\mathrm{\hspace{0.17em}10})$
$\ln \left({\omega }_0\right)$	$𝒰(-10,\mathrm{\hspace{0.17em}10})$

Table 2. Table 2: The parameters and priors for Example 2.

parameter	prior
$\ln \left(S_0\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}5})$
$\ln \left(Q\right)$	$𝒰(-10,\mathrm{\hspace{0.17em}10})$
$\ln \left({\omega }_1\right)$	$𝒰(-10,\mathrm{\hspace{0.17em}10})$
$\ln \left({\omega }_2\right)$	$𝒰(-5,\mathrm{\hspace{0.17em}5})$

Table 3. Table 3: The parameters and priors for Example 3.

parameter	prior
$\ln \left(B/{\mathrm{ppt}}^2\right)$	$𝒰(-10.0,\mathrm{\hspace{0.17em}0.0})$
$\ln \left(L/\mathrm{day}\right)$	$𝒰(1.5,\mathrm{\hspace{0.17em}5.0})$
$\ln \left(P/\mathrm{day}\right)$	$𝒰(-3.0,\mathrm{\hspace{0.17em}5.0})$
$\ln \left(C\right)$	$𝒰(-5.0,\mathrm{\hspace{0.17em}5.0})$

Table 4. Table 4: The parameters and priors for Example 4.

parameter	prior
$\ln \left(S_g/{\mathrm{ppm}}^2\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}15})$
$\ln \left({\omega }_g/{\mathrm{day}}^{-1}\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}15})$
$\ln \left({\nu }_{\max }/\mu \mathrm{Hz}\right)$	$𝒰(\ln \left(130\right),\ln \left(190\right))$
$\ln \left(\Delta \nu /\mu \mathrm{Hz}\right)$	$𝒰(\ln \left(12.5\right),\ln \left(13.5\right))$
$ϵ/{\mathrm{day}}^{-1}$	$𝒩(0,\mathrm{\hspace{0.17em}1})$
$\ln \left(A/{\mathrm{ppm}}^2\mathrm{day}\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}15})$
$\ln \left(Q\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}15})$
$\ln \left(W/{\mathrm{day}}^{-1}\right)$	$𝒰(-3,\mathrm{\hspace{0.17em}3})$
$\ln \left(\sigma /\mathrm{ppm}\right)$	$𝒰(-15,\mathrm{\hspace{0.17em}15})$

Table 5. Table 5: The parameters, priors, and (if known) the true values used for the simulation in Example 5.

	parameter	prior	true value
kernel: $𝜶$	$\ln \left(B/{\mathrm{ppt}}^2\right)$	$𝒰(-10.0,\mathrm{\hspace{0.17em}0.0})$	—
	$\ln \left(L/\mathrm{day}\right)$	$𝒰(1.5,\mathrm{\hspace{0.17em}5.0})$	—
	$\ln \left(P_{\mathrm{rot}}/\mathrm{day}\right)$	$𝒰(-3.0,\mathrm{\hspace{0.17em}5.0})$	—
	$\ln \left(C\right)$	$𝒰(-5.0,\mathrm{\hspace{0.17em}5.0})$	—
	$\ln \left(\sigma /\mathrm{ppt}\right)$	$𝒰(-5.0,\mathrm{\hspace{0.17em}0.0})$	—
mean: $𝜽$	$f_0/\mathrm{ppt}$	$𝒰(-0.5,\mathrm{\hspace{0.17em}0.5})$	0
	$\ln \left(P_{\mathrm{orb}}/\mathrm{day}\right)$	$𝒰(\ln \left(7.9\right),\ln \left(8.1\right))$	$\ln \left(8\right)$
	$\ln \left(R_p/R_{\star }\right)$	$𝒰(\ln \left(0.005\right),\ln \left(0.1\right))$	$\ln \left(0.015\right)$
	$\ln \left(T/\mathrm{day}\right)$	$𝒰(\ln \left(0.4\right),\ln \left(0.6\right))$	$\ln \left(0.5\right)$
	$t_0/\mathrm{day}$	$𝒰(-0.1,\mathrm{\hspace{0.17em}0.1})$	$0$
	$b$	$𝒰(0,\mathrm{\hspace{0.17em}1.0})$	$0.5$
	$q_1$	$𝒰(0,\mathrm{\hspace{0.17em}1})$	$0.5$
	$q_2$	$𝒰(0,\mathrm{\hspace{0.17em}1})$	$0.5$

Table 6. Table 6: The computational cost and convergence stats for each example.

	$N$	$J$	directaaThe computational cost of computing the GP model using the general Cholesky factorization routine implemented in the Intel MKL.	celeritebbThe computational cost of computing the GP model using celerite.	dimensionccThe total number of parameters in the model.	evaluationsddThe total number of evaluations of the model to run MCMC to convergence.	$N_{\mathrm{eff}}$eeThe effective number of independent samples estimated by computing the integrated autocorrelation time of the chain.
			ms	ms
1	200	1	2.85	0.26	3	80000	1737
2	100	2	0.67	0.36	4	80000	1134
3	6950	2	8119.11	1.47	4	176000	2900
4	1440	8	828.93	2.74	9	640000	1443
5	1000	2	88.30	0.96	13	1280000	1490

Equations

$$\displaystyle{\boldsymbol{y}}=\left(\begin{array}[]{ccccc}y_{1}&&\cdots&&y_{N}\end{array}\right)^{\mathrm{T}}$$

$$\displaystyle X=\left(\begin{array}[]{ccccc}{\boldsymbol{x}}_{1}&&\cdots&&{\boldsymbol{x}}_{N}\end{array}\right)^{\mathrm{T}}$$

$$\ln \mathcal{L}(\mathbf{\theta },\mathbf{\alpha })=\ln p(\mathbf{y}|X,\mathbf{\theta },\mathbf{\alpha })=-\frac{1}{2}{{\mathbf{r}}_{\mathbf{\theta }}}^{\mathrm{T}}{K_{\mathbf{\alpha }}}^{-1}{\mathbf{r}}_{\mathbf{\theta }}-\frac{1}{2}\ln \det K_{\mathbf{\alpha }}-\frac{N}{2}\ln (2\pi )$$

$$\displaystyle{\boldsymbol{r}}_{\boldsymbol{\theta}}=\left(\begin{array}[]{ccccc}y_{1}-\mu_{\boldsymbol{\theta}}({\boldsymbol{x}}_{1})&&\cdots&&y_{N}-\mu_{\boldsymbol{\theta}}({\boldsymbol{x}}_{N})\end{array}\right)^{\mathrm{T}}$$

$$k_{\mathbf{\alpha }}({\tau }_{nm})={\sigma }_n^2{\delta }_{nm}+a\exp (-c{\tau }_{nm})$$

$$k_{\mathbf{\alpha }}({\tau }_{nm})={\sigma }_n^2{\delta }_{nm}+\sum _{j=1}^J{a}_j\exp (-c_j{\tau }_{nm}).$$

$$k_{\mathbf{\alpha }}({\tau }_{nm})={\sigma }_n^2{\delta }_{nm}+\sum _{j=1}^J$$

$$S(\omega )=\sum _{j=1}^J\sqrt{\frac{2}{\pi }}\frac{(a_j{c}_j+b_j{d}_j)({c_j}^2+{d_j}^2)+(a_j{c}_j-b_j{d}_j){\omega }^2}{{\omega }^4+2({c_j}^2-{d_j}^2){\omega }^2+{({c_j}^2+{d_j}^2)}^2}.$$

$$k_j({\tau }_{nm})=a_j\exp \left(-c_j{\tau }_{nm}\right)\cos \left(d_j{\tau }_{nm}\right)$$

$$S_j(\omega )=\frac{1}{\sqrt{2\pi }}\frac{a_j}{c_j}\left[\frac{1}{1+{\left(\frac{\omega -d_j}{c_j}\right)}^2}+\frac{1}{1+{\left(\frac{\omega +d_j}{c_j}\right)}^2}\right].$$

$$k_j({\tau }_{nm})=a_j\exp \left(-c_j{\tau }_{nm}\right)$$

$$S_j(\omega )=\sqrt{\frac{2}{\pi }}\frac{a_j}{c_j}\frac{1}{1+{\left(\frac{\omega }{c_j}\right)}^2}.$$

$$k_j(\tau )k_k(\tau )=e^{-\tilde{c}\tau }[{\tilde{a}}_{+}\cos ({\tilde{d}}_{+}\tau )+{\tilde{b}}_{+}\sin ({\tilde{d}}_{+}\tau )+{\tilde{a}}_{-}\cos ({\tilde{d}}_{-}\tau )+{\tilde{b}}_{-}\sin ({\tilde{d}}_{-}\tau )]$$

$${\tilde{a}}_{\pm }$$

$${\tilde{b}}_{\pm }$$

$$\tilde{c}$$

$${\tilde{d}}_{\pm }$$

$$\left[\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}+\frac{{\omega }_0}{Q}\frac{\mathrm{d}}{\mathrm{d}t}+{\omega }_0^2\right]y(t)=ϵ(t)$$

$$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0{\omega }_0^4}{({\omega }^2-{\omega }_0^2{)}^2+{\omega }_0^2{\omega }^2/Q^2}$$

$$a_j$$

$$b_j$$

$$c_j$$

$$d_j$$

$$a_{j\pm }$$

$$b_{j\pm }$$

$$c_{j\pm }$$

$$d_{j\pm }$$

$$k_{\mathrm{SHO}}(\tau ;S_0,Q,{\omega }_0)=S_0{\omega }_{0}Q{e}^{-\frac{{\omega }_0\tau }{2Q}}\{\begin{array}{ll}\cosh (\eta {\omega }_0\tau )+\frac{1}{2\eta Q}\sinh (\eta {\omega }_0\tau ),& 0<Q<1/2\\ 2(1+{\omega }_0\tau ),& Q=1/2\\ \cos (\eta {\omega }_0\tau )+\frac{1}{2\eta Q}\sin (\eta {\omega }_0\tau ),& 1/2<Q\end{array}$$

$$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0}{(\omega /{\omega }_0{)}^4+1}.$$

$$k(\tau )=S_0{\omega }_0{e}^{-\frac{1}{\sqrt{2}}{\omega }_0\tau }\cos \left(\frac{{\omega }_0\tau }{\sqrt{2}}-\frac{\pi }{4}\right).$$

$$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0}{{\left[(\omega /{\omega }_0{)}^2+1\right]}^2}$$

$$k(\tau )$$

$$=$$

$$k(\tau )$$

Full text

methods: data analysis — methods: statistical — asteroseismology — stars: rotation — planetary systems

††software: AstroPy (Astropy Collaboration et al., 2013), corner.py (Foreman-Mackey, 2016), Eigen (Guennebaud et al., 2010), emcee (Foreman-Mackey et al., 2013), george (Ambikasaran et al., 2016), Julia (Bezanzon et al., 2012), LAPACK (Anderson et al., 1999), matplotlib (Hunter et al., 2007), numpy (Van Der Walt et al., 2011), transit (Foreman-Mackey & Morton, 2016), scipy (Jones et al., 2001).

Section 1 Introduction

In the astrophysical literature, Gaussian Processes (GPs; Rasmussen & Williams, 2006) have been used to model stochastic variability in light curves of stars (Brewer & Stello, 2009), active galactic nuclei (Kelly et al., 2014), and the logarithmic flux of X-ray binaries (Uttley et al., 2005). They have also been used as models for the cosmic microwave background (Bond & Efstathiou, 1987; Bond et al., 1999; Wandelt & Hansen, 2003), correlated instrumental noise (Gibson et al., 2012), and spectroscopic calibration (Czekala et al., 2017; Evans et al., 2015). While these models are widely applicable, their use has been limited, in practice, by the computational cost and scaling. The cost of computing a general GP likelihood scales as the cube of the number of data points $\mathcal{O}(N^{3})$ and in the current era of large time domain surveys – with as many as $\sim 10^{4-9}$ targets with $\sim 10^{3-5}$ observations each – this scaling is prohibitive.

Existing astronomical time series datasets have already reached the limit where naïve application GP models is no longer tractable. NASA’s Kepler Mission (Borucki et al., 2010), for example, measured light curves with more than 60,000 observations each for about 190,000 stars. Current and forthcoming surveys such as K2 (Howell et al., 2014), TESS (Ricker et al., 2014), LSST (Ivezić et al., 2008), WFIRST (Spergel et al., 2015), and PLATO (Rauer et al., 2014) will continue to produce similar or larger data volumes.

In this paper, we present a method for directly and exactly computing a class of GP models that scales linearly with the number of data points $\mathcal{O}(N)$ for one dimensional data sets. Unlike earlier linear methods using Kalman filters (for example, Kelly et al., 2014), this novel algorithm exploits the semiseparable structure of a specific class of covariance matrices to directly factorize and solve the system.

This method can only be used with one-dimensional datasets and the covariance function must be represented by a mixture of exponentials; we will return to a discussion of what this means in detail in the following sections. However, the measurements don’t need to be evenly spaced and the uncertainties can be heteroscedastic.

This method achieves linear scaling by exploiting structure in the covariance matrix when it is generated by a mixture of exponentials. The semiseparable nature of these matrices was first recognized by Ambikasaran (2015), building on intuition from a twenty year old paper (Rybicki & Press, 1995). As we will discuss in the following pages, this choice of kernel function arises naturally in physical systems and we demonstrate it can be used as an effective model in other cases. This method is especially appealing compared to other similar methods – we discuss these comparisons in Section 7 – because it is exact, simple, and fast.

Our main expertise lies in the field of exoplanet discovery and characterization where GPs have become a model of choice. We are confident that this method will benefit this field, but we also expect that there will be applications in other branches of astrophysics and beyond. In Section 6, we present applications of the method to research problems in stellar rotation (Section 6.3), asteroseismic analysis (Section 6.4), and exoplanet transit fitting (Section 6.5). Some readers might consider starting with these section for motivating examples before delving into the detailed derivations of the earlier sections.

In the following pages, we motivate the general problem of GP modeling, derive our novel direct solver, and demonstrate the model’s application on real and simulated data sets. Alongside this paper, we have released well-tested and documented open source implementations written in C++, Python, and Julia. These are available online at GitHub https://github.com/dfm/celerite and Zenodo (Foreman-Mackey et al., 2017).

Section 2 Gaussian processes

GPs are stochastic models consisting of a mean function $\mu_{\boldsymbol{\theta}}({\boldsymbol{x}})$ and a covariance, autocorrelation, or “kernel” function $k_{\boldsymbol{\alpha}}({\boldsymbol{x}}_{n},\,{\boldsymbol{x}}_{m})$ parameterized by ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$ respectively. Under this model, the log-likelihood function $\mathcal{L}({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}})$ with a dataset

[TABLE]

at coordinates

[TABLE]

is

[TABLE]

where

[TABLE]

is the vector of residuals and the elements of the covariance matrix $K$ are given by $[K_{\boldsymbol{\alpha}}]_{nm}=k_{\boldsymbol{\alpha}}({\boldsymbol{x}}_{n},\,{\boldsymbol{x}}_{m})$. Equation (5) is the equation of an $N$-dimensional Gaussian and it can be derived as the generalization of what we astronomers call “$\chi^{2}$” to include the effects of correlated noise.

A point estimate for the values of the parameters ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$ for a given dataset $({\boldsymbol{y}},\,X)$ can be found by maximizing Equation (5) with respect to ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$ using a non-linear optimization routine (Nocedal & Wright, 2006). Furthermore, the uncertainties on ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$ can be quantified by multiplying the likelihood by a prior $p({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}})$ and using a Markov chain Monte Carlo (MCMC) algorithm to sample the joint posterior probability density.

GP models have been widely used across the physical sciences but their application is generally limited to small datasets, $N\lesssim 10^{3}$, because the cost of computing the inverse and determinant of a general matrix $K_{\boldsymbol{\alpha}}$ is $\mathcal{O}(N^{3})$. In other words, this cost is proportional to the cube of the number of data points $N$. Furthermore, the storage requirements also scale as $\mathcal{O}(N^{2})$. This means that for large datasets, every evaluation of the likelihood quickly becomes computationally intractable, and representing the matrix in memory can become expensive. In this case, the use of non-linear optimization or MCMC is no longer practical.

This is not a purely academic concern. While the cost of directly evaluating the GP likelihood using a tuned linear algebra library111For comparisons throughout this paper, we use the Intel Math Kernel Library (https://software.intel.com/en-us/intel-mkl) running on a recent Apple MacBook Pro. with a dataset of $N=1024$ measurements is less than a tenth of a second, the same calculation with $N=8192$ takes over 8 seconds. Furthermore, most optimization and sampling routines are iterative, requiring many evaluations of this model to perform parameter estimation or inference. Existing datasets from astronomical surveys like Kepler and K2 include tens of thousands of observations for hundreds of thousands of targets, and upcoming surveys like LSST are expected to produce thousands or tens of thousands of measurements each for billions of astronomical sources (Ivezić et al., 2008). Scalable methods will be required if we ever hope to use GP models with these datasets.

In this paper, we present a method for improving the cubic scaling for some use cases. We call our method and its implementations celerite.222The name celerite comes from the French word célerité meaning the speed of light in a vacuum. Throughout the paper, when referring to the celerite model, we typeset celerite in italics and, in reference to the implementation of this model, we typeset celerite in sans-serif. The celerite method requires using a specific model for the covariance $k_{\boldsymbol{\alpha}}({\boldsymbol{x}}_{n},\,{\boldsymbol{x}}_{m})$ and, although it has some limitations, we demonstrate in subsequent sections that it can increase the computational efficiency of many astronomical data analysis problems. The main limitation of this method is that it can only be applied to one-dimensional datasets, where by “one-dimensional” we mean that the input coordinates ${\boldsymbol{x}}_{n}$ are scalar, ${\boldsymbol{x}}_{n}\equiv t_{n}$. We use $t$ as the input coordinate because one-dimensional GPs are often applied to time series data but this isn’t a real restriction and the celerite method can be applied to any one-dimensional dataset. Furthermore, the covariance function for the celerite method is “stationary”. In other words, $k_{\boldsymbol{\alpha}}(t_{n},\,t_{m})$ is only a function of $\tau_{nm}\equiv|t_{n}-t_{m}|$.

Section 3 The celerite model

To scale GP models to larger datasets, Rybicki & Press (1995) presented a method of computing the first term in Equation (5) in $\mathcal{O}(N)$ operations when the covariance function is given by

[TABLE]

where $\{{\sigma_{n}}^{2}\}_{n=1}^{N}$ are the measurement uncertainties, $\delta_{nm}$ is the Kronecker delta, and ${\boldsymbol{\alpha}}=(a,\,c)$. The intuition behind this method is that, for this choice of $k_{\boldsymbol{\alpha}}$, the inverse of $K_{\boldsymbol{\alpha}}$ is tridiagonal and can be computed with a small number of operations for each data point. This model has been generalized to arbitrary mixtures of exponentials (for example, Kelly et al., 2011)

[TABLE]

In this case, the inverse is dense but Equation (5) can still be evaluated with a scaling of $\mathcal{O}(N\,J^{2})$, where $J$ is the number of components in the mixture (Kelly et al., 2014; Ambikasaran, 2015). In Section 5, we build on previous work (Ambikasaran, 2015) to derive a faster method with the same $\mathcal{O}(N\,J^{2})$ scaling that can be used when $k_{\boldsymbol{\alpha}}(\tau_{nm})$ is positive definite.

This kernel function can be made even more general by introducing complex parameters $a_{j}\to a_{j}\pm i\,b_{j}$ and $c_{j}\to c_{j}\pm i\,d_{j}$. In this case, the covariance function becomes

[TABLE]

and, for this function, Equation (5) can still be evaluated with $\mathcal{O}(N\,J^{2})$ operations. The details of this method and a few implementation considerations are outlined in Section 5, but we first discuss some properties of this covariance function.

By rewriting the exponentials in Equation (10) as sums of sine and cosine functions, we can see the autocorrelation structure is defined by a mixture of quasiperiodic oscillators

[TABLE]

For clarity, we refer to the argument within the sum as a “celerite term” for the remainder of this paper. The Fourier transform333Here and throughout we have defined the Fourier transform of the function $f(t)$ as $F(\omega)={(2\,\pi)}^{-1/2}\,\int_{-\infty}^{\infty}f(t)\,e^{i\,\omega\,t}\,\mathrm{d}t$. of this covariance function is the power spectral density (PSD) of the process and it is given by

[TABLE]

The physical interpretation of this model isn’t immediately obvious and we return to a more general discussion shortly but we start by considering some special cases.

If we set the imaginary amplitude $b_{j}$ for some component $j$ to zero, that term of Equation (11) becomes

[TABLE]

and the PSD for this component is

[TABLE]

This is the sum of two Lorentzian or Cauchy distributions with width $c_{j}$ centered on $\omega=\pm d_{j}$. This model can be interpreted intuitively as a quasiperiodic oscillator with amplitude $A_{j}=a_{j}$, quality factor $Q_{j}=d_{j}\,{(2\,c_{j})}^{-1}$, and period $P_{j}=2\,\pi\,{d_{j}}^{-1}$.

Similarly, setting both $b_{j}$ and $d_{j}$ to zero, we get

[TABLE]

with the PSD

[TABLE]

This model is often called an Ornstein–Uhlenbeck process (in reference to the classic paper, Uhlenbeck & Ornstein, 1930) and this is the kernel that was studied by Rybicki & Press (Rybicki & Press, 1992, 1995).

Finally, we note that the product of two terms of the form found inside the sum in Equation (11) can also be re-written as a sum with updated parameters

[TABLE]

where

[TABLE]

Therefore, the method described in Section 5 can be used to perform scalable inference on large datasets for any model, where the kernel function is a sum or product of celerite terms.

Section 4 celerite as a model of stellar variations

We now turn to a discussion of celerite as a model of astrophysical variability. A common concern in the context of GP modeling in astrophysics is the lack of physical motivation for the choice of kernel functions. Kernels are often chosen simply because they are popular, with little consideration of the impact of this decision. In this section, we discuss an exact physical interpretation of the celerite kernel that is applicable to many astrophysical systems, but especially the time-variability of stars.

Many astronomical objects are variable on timescales determined by their physical structure. For phenomena such as stellar (asteroseismic) oscillations, variability is excited by noisy physical processes and grows most strongly at the characteristic timescale but is also damped due to dissipation in the system. These oscillations are strong at resonant frequencies determined by the internal stellar structure, which are both excited and damped by convective turbulence.

To relate this physical picture to celerite, we consider the dynamics of a stochastically-driven damped simple harmonic oscillator (SHO). The differential equation for this system is

[TABLE]

where $\omega_{0}$ is the frequency of the undamped oscillator, $Q$ is the quality factor of the oscillator, and $\epsilon(t)$ is a stochastic driving force. If $\epsilon(t)$ is white noise, the PSD of this process is (Anderson et al., 1990)

[TABLE]

where $S_{0}$ is proportional to the power at $\omega=\omega_{0}$, $S(\omega_{0})=\sqrt{2/\pi}\,S_{0}\,Q^{2}$. The power spectrum in Equation (23) matches Equation (12) if

[TABLE]

for $Q\geq\frac{1}{2}$. For $0<Q\leq\frac{1}{2}$, Equation (23) can be captured by a pair of celerite terms with parameters

[TABLE]

For these definitions, the kernel is

[TABLE]

where $\eta=|1-(4\,Q^{2})^{-1}|^{1/2}$. Because of the damping, the characteristic oscillation frequency in this model, $d_{j}$, for any finite quality factor $Q>1/2$, is not equal to the frequency of the undamped oscillator, $\omega_{0}$.

The power spectrum in Equation (23) has several limits of physical interest:

• •

For $Q=1/\sqrt{2}$, Equation (23) simplifies to

[TABLE]

This functional form is commonly used to model for the background granulation noise in asteoreseismic and helioseismic analyses (Harvey, 1985; Michel et al., 2009; Kallinger et al., 2014). The kernel for this value of $Q$ is

[TABLE]

• •

Substituting $Q=1/2$, Equation (23) becomes

[TABLE]

with the corresponding covariance function (using Equation 11 and Equation 25)

[TABLE]

or, equivalently (using Equation 11 and Equation 24)

[TABLE]

This covariance function is also known as the Matérn-3/2 function (Rasmussen & Williams, 2006). This model cannot be directly evaluated using celerite because, as we can see from Equation (24), the parameter $b_{j}$ is infinite for $Q=1/2$. The limit in Equation (32), however, suggests that an approximate Matérn-3/2 covariance can be computed with celerite by using a small value of $f$ in Equation (32). The required value of $f$ will depend on the specific dataset and the precision requirements for the inference, so we encourage users to test this approximation in the context of their research before using it.

• •

Finally, in the limit of large $Q$, the model approaches a high quality oscillation with frequency $\omega_{0}$ and covariance function

[TABLE]

Figure 1 shows a plot of the PSD for these limits and several other values of $Q$. From this figure, it is clear that for $Q\leq 1/2$, the model has no oscillatory behavior and that for large $Q$, the shape of the PSD near the peak frequency approaches a Lorentzian.

These special cases demonstrate that the stochastically-driven simple harmonic oscillator provides a physically motivated model that is flexible enough to describe a wide range of stellar variations and we return to give some specific examples in Section 6. Low $Q\approx 1$ can capture granulation noise and high $Q\gg 1$ is a good model for asteroseismic oscillations. In practice, we take a sum over oscillators with different values of $Q$, $S_{0}$, and $\omega_{0}$ to give a sufficient accounting of the power spectrum of stellar time series. Since this kernel is exactly described by the exponential kernel, the likelihood (Equation 5) can be evaluated for a time series with $N$ measurements in $\mathcal{O}(N)$ operations using the celerite method described in the next section.

Section 5 Semiseparable matrices & celerite

The previous two sections describe the celerite model and its physical interpretation. In this section, we derive a new scalable direct solver for covariance matrices generated by this model to compute the GP likelihood (Equation 5) for large datasets. This method exploits the semiseparable structure of these matrices to compute their Cholesky factorization in $\mathcal{O}(N\,J^{2})$ operations.

The relationship between the celerite model and semiseparable linear algebra was first recognized by Ambikasaran (Ambikasaran, 2015) building on earlier work by Rybicki & Press (Rybicki & Press, 1995). Ambikasaran derived a scalable direct solver for any general semiseparable matrix and applied this method to celerite models with $b_{j}=0$ and $d_{j}=0$. During the preparation of this paper, we generalized this earlier method to include non-zero $b_{j}$ and $d_{j}$, but it turns out that we can derive a higher performance algorithm by restricting our method to positive definite matrices. In Appendix A, we discuss methods for ensuring that a celerite model is positive definite.

To start, we define a rank-$R$ semiseparable matrix as a matrix $K$ where the elements can be written as

[TABLE]

for some $N\times R$-dimensional matrices $U$ and $V$, and an $N\times N$-dimensional diagonal matrix $A$. This definition can be equivalently written as

[TABLE]

where the $\mathrm{tril}$ function extracts the strict lower triangular component of its argument and $\mathrm{triu}$ is the equivalent strict upper triangular operation. In this representation, rather than storing the full $N\times N$ matrix $K$, it is sufficient to only store the $N\times R$ matrices $U$ and $V$, and the $N$ diagonal entries of $A$, for a total storage volume of $(2\,R+1)\,N$ floating point numbers. Given a semiseparable matrix of this form, the matrix-vector products and matrix-matrix products can be evaluated in $\mathcal{O}(N)$ operations. We include a summary of the key operations in Appendix B, but the interested reader is directed to Vandebril et al. (2007) for more details about semiseparable matrices in general.

5.1 Cholesky factorization of positive definite semiseparable matrices

For our purposes, we want to compute the log-determinant of $K$ and solve linear systems of the form $K\,{\boldsymbol{z}}={\boldsymbol{y}}$ for ${\boldsymbol{z}}$, where $K$ is a rank-$R$ semiseparable positive definite matrix. These can both be computed using the Cholesky factorization of $K$:

[TABLE]

where $L$ is a lower triangular matrix with unit diagonal and $D$ is a diagonal matrix. We use the ansatz that $L$ has the form

[TABLE]

where $I$ is the identity, $U$ is the matrix from above, and $W$ is an unknown $N\times R$ matrix. Combining this with Equation (39), we find the equation

[TABLE]

that we can solve for $W$. Setting each element on the left side of Equation (5.1) equal to the corresponding element on the right side, we derive the following recursion relations

[TABLE]

where $S_{n}$ is a symmetric $R\times R$ matrix and every element $S_{1,j,k}=0$. The computational cost of one iteration of the update in Equation (43) is $\mathcal{O}(R^{2})$ and this must be run $N$ times – once for each row – so the total cost of this factorization scales as $\mathcal{O}(N\,R^{2})$.

After computing the Cholesky factorization of this matrix $K$, we can also apply its inverse and compute its log-determinant in $\mathcal{O}(N\,R)$ and $\mathcal{O}(N)$ operations respectively. The log-determinant of $K$ is given by

[TABLE]

where, since $K$ is positive definite, $D_{n,n}>0$ for all $n$.

The inverse of $K$ can be applied using the relationship

[TABLE]

First, we solve for ${\boldsymbol{z}}^{\prime}=L^{-1}\,{\boldsymbol{y}}$ using forward substitution

[TABLE]

where $f_{0,j}=0$ for all $j$. Then, using this result and backward substitution, we compute ${\boldsymbol{z}}$

[TABLE]

where $g_{N+1,j}=0$ for all $j$. The total cost of the forward and backward substitution passes scales as $\mathcal{O}(N\,R)$.

The algorithm derived here is a general method for factorizing, solving, and computing the determinant of positive definite rank-$R$ semiseparable matrices. As we will show below, this method can be used to compute Equation (5) for celerite models. In practice, when $K$ is positive definite and this method is applicable, we find that this method is about a factor of 20 faster than an optimized implementation of the earlier, more general method from Ambikasaran (2015).

5.2 Scalable computation of celerite models

The celerite model described in Section 3 can be exactly represented as a rank $R=2\,J$ semiseparable matrix444We note that if some celerite terms are real ($b_{j}=0$ and $d_{j}=0$) then these terms only add one to the rank. Therefore, if a celerite model has $J_{\mathrm{real}}$ real and $J_{\mathrm{complex}}$ general terms, the semiseparable rank is only $R=J_{\mathrm{real}}+2\,J_{\mathrm{complex}}=2\,J-J_{real}$. where the components of the relevant matrices are

[TABLE]

A naïve implementation of the algorithm in Equation (43) for this system will result in numerical overflow and underflow issues in many realistic cases because both $U$ and $V$ have factors of the form $e^{\pm c_{j}t}$, where $t$ can be any arbitrarily large real number. This issue can be avoided by reparameterizing the semiseparable representation from Equation (5.2). We define the following pre-conditioned variables

[TABLE]

and the new set of variables

[TABLE]

with the constraint

[TABLE]

Using this parameterization, we find the following, numerically stable, algorithm to compute the Cholesky factorization of a celerite model

[TABLE]

As before, a single iteration of this algorithm requires $\mathcal{O}(J^{2})$ operations and it is repeated $N$ times for a total scaling of $\mathcal{O}(N\,J^{2})$. Furthermore, $\tilde{V}_{n,j}$ and $A_{n,n}$ can be updated in place and used to represent $\tilde{W}_{n,j}$ and $D_{n,n}$, giving a total memory footprint of $(6\,J+1)\,N+J\,(J-1)/2$ floating point numbers.

The log-determinant of $K$ can still be computed as in the previous section, but we must update the forward and backward substitution algorithms to account for the reparameterization in Equation (5.2). The forward substitution with the reparameterized variables is

[TABLE]

where $f_{0,j}=0$ for all $j$, and the backward substitution algorithm becomes

[TABLE]

where $g_{N+1,j}=0$ for all $j$. The only differences, when compared to Equations (5.1) and (5.1) is the multiplication by $\phi_{n,j}$ in the recursion relations for $f_{n,j}$ and $g_{n,j}$. As before, the computational cost of this solve scales as $\mathcal{O}(N\,J)$.

5.3 Performance

We have derived a novel and scalable method for computing the Cholesky factorization of the covariance matrices produced by celerite models. In this section, we demonstrate the numerical and computational performance of this algorithm.

First, we empirically confirm that the method solves systems of the form $K\,{\boldsymbol{z}}={\boldsymbol{y}}$ to high numerical accuracy. The left panel of Figure 2 shows the L-infinity norm for the residual ${\boldsymbol{z}}-K^{-1}\,{\boldsymbol{y}}={\boldsymbol{z}}-K^{-1}K{\boldsymbol{z}}$ for a range of dataset sizes $N$ and numbers of terms $J$. For each pair of $N$ and $J$ we repeated the experiment 10 times with different randomly generated datasets and celerite parameters, and averaged across experiments to remove some sampling noise from this figure. The numerical error introduced by our method increases weakly for larger datasets $\sim N^{0.15}$ and roughly linearly with $J$.

For the experiments with small $N\leq 2048$, we also compared the log determinant calculation to the result calculated using the general log determinant function provided by the NumPy project (Van Der Walt et al., 2011), itself built on the LAPACK (Anderson et al., 1999) implementation in the Intel Math Kernel Library555https://software.intel.com/en-us/intel-mkl. The right-hand panel of Figure 2 shows a histogram of the fractional error in the log determinant calculation for each of the 480 experiments with $N\leq 2048$. These results agree to about $10^{-15}$ and the error only weakly increases with $N$ and $J$.

We continue by measuring the real-world computational efficiency of this method as a function of $N$ and $J$. This experiment was run on a 2013 MacBook Pro with a dual 2.6 GHz Intel Core i5 processor, but we find similar performance on other platforms. The celerite implementation used for this test is the Python interface with a C++ back end, but there is little overhead introduced by Python so we find similar performance using the implementation in C++ and an independent implementation in Julia (Bezanzon et al., 2012). Figure 3 shows the computational cost of evaluating Equation (5) using celerite as a function of $N$ and $J$. The empirical scaling as a function of data size matches the theoretical scaling of $\mathcal{O}(N)$ for all values of $N$ and the scaling with the number of terms is the theoretical $\mathcal{O}(J^{2})$ for $J\gtrsim 10$, with some overhead for smaller models. For comparison, the computational cost of the same computation using the Intel Math Kernel Library, a highly optimized general dense linear algebra package, is shown as a dashed line, which is independent of $J$. celerite is faster than the reference implementation for all but the smallest datasets with many terms, $J\gtrsim 64$ for $N\approx 512$, and $J$ increasing with larger $N$.

5.4 How to choose a celerite model

In any application of celerite – or any other GP model, for that matter – we must choose a kernel. A detailed discussion of model selection is beyond the scope of this paper and the interested reader is directed to Chapter 5 of Rasmussen & Williams (2006) for an introduction in the context of GPs, but we include a summary and some suggestions specific to celerite here.

For many of the astrophysics problems that we will tackle using celerite, it is possible to motivate the kernel choice physically. We give three examples of this in the next section and, in these cases, we have a physical description of the system that generated our data, we represent that system as a celerite model, and use celerite for parameter estimation. For other projects, we might be interested in comparing different theories and, in those case, we would need to perform formal model selection.

Sometimes, it is not possible to choose a physically motivated kernel when we use a GP as an effective model. In these cases, we recommend starting with a mixture of stochastically-driven damped SHOs (as discussed in Section 4) and select the number of oscillators using Bayesian model comparison, cross-validation, or another model selection technique. Despite the fact that it is less flexible than a general celerite kernel, we recommend the SHO model in most cases because it is once mean square differentiable

[TABLE]

giving smoother functions than the general celerite model that is only mean square continuous.

For a similar class of models – we return to these below in Section 7 – Kelly et al. (2014) recommend using an “information criterion”, such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC), as a model selection metric. These criteria are easy to implement and we use the BIC in Section 6.2 but, since these are only formally valid under specific sets of assumptions, we cannot recommend their use in general.

Section 6 Examples

To demonstrate the application of celerite, we use it to perform posterior inference in two examples with realistic simulated data and three examples of real-world research problems. These examples have been chosen to justify the use of celerite for a range of research problems, but they are by no means exhaustive. These examples are framed in the time domain with a clear bias in favor of large homogeneous photometric surveys, but celerite can also be used for other one-dimensional problems like, for example, spectroscopy, where wavelength – instead of time – is the independent coordinate (see Czekala et al., 2017, for a potential application).

In the first example, we demonstrate that celerite can be used to infer the power spectrum of a process when the data are generated from a celerite model. In the second example, we demonstrate that celerite can be used as an effective model even if the true process cannot be represented in the space of allowed models. This is an interesting example because, when analyzing real data, we rarely have any fundamental reason to believe that the data were generated by a GP model with a specific kernel. Even in these cases, GPs can be useful effective models and celerite provides computational advantages over other GP methods.

The next three examples show celerite used to infer the properties of stars and exoplanets observed by NASA’s Kepler Mission. Each of these examples touches on an active area of research so we limit our examples to be qualitative in nature and do not claim that celerite is the optimal method, but we hope these examples encourage interested readers to investigate the applicability of celerite to their research.

In each example, the joint posterior probability density is given by

[TABLE]

where $p({\boldsymbol{y}}\,|\,{X,\,{\boldsymbol{\theta}}},\,{\boldsymbol{\alpha}})$ is the GP likelihood defined in Equation (5) and $p({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}})$ is the joint prior probability density for the parameters. We sample each posterior density using MCMC and investigate the performance of the model when used to interpret the physical processes that generated the data. The specific parameters and priors are discussed in detail in each section, but we generally assume separable uniform priors with plausibly broad support.

6.1 Example 1: Recovery of a celerite process

In this first example, we simulate a dataset using a known celerite process and fit it with celerite to show that valid inferences can be made in this idealized case. We simulate a small ($N=200$) dataset using a GP model with a SHO kernel (Equation 26) with parameters $S_{0}=1\,\mathrm{ppm}^{2}$, $\omega_{0}=e^{2}\,\mathrm{day}^{-1}$, and $Q=e^{2}$, where the units are arbitrarily chosen to simplify the discussion. We add further white noise with the amplitude $\sigma_{n}=2.5\,\mathrm{ppm}$ to each data point. The simulated data are plotted in the top left panel of Figure 4.

In this case, when simulating and fitting, we set the mean function $\mu_{\boldsymbol{\theta}}$ to zero – this means that the parameter vector ${\boldsymbol{\theta}}$ is empty – and the elements of the covariance matrix are given by Equation (26) with three parameters ${\boldsymbol{\alpha}}=(S_{0},\,\omega_{0},\,Q)$. We choose a proper separable prior for ${\boldsymbol{\alpha}}$ with log-uniform densities for each parameter as listed in Table 1.

To start, we estimate the maximum a posteriori (MAP) parameters using the L-BFGS-B non-linear optimization routine (Byrd et al., 1995; Zhu et al., 1997) implemented by the SciPy project (Jones et al., 2001). The top left panel of Figure 4 shows the conditional mean and standard deviation of the MAP model over-plotted on the simulated data and the bottom panel shows the residuals away from this model. Since we are using the correct model to fit the data, it is reassuring that the residuals appear qualitatively uncorrelated in this figure.

We then sample the joint posterior probability (Equation 56) using emcee (Goodman & Weare, 2010; Foreman-Mackey et al., 2013). We initialize 32 walkers by sampling from a three-dimensional isotropic Gaussian centered the MAP parameter vector and with a standard deviation of $10^{-4}$ in each dimension. We then run 500 steps of burn-in and 2000 steps of MCMC. To assess convergence, we estimate the mean integrated autocorrelation time of the chain across parameters (Sokal, 1989; Goodman & Weare, 2010) and find that the chain results in 1737 effective samples.

Each sample in the chain corresponds to a model PSD and we compare this posterior constraint on the PSD to the true spectral density in the right-hand panel of Figure 4. In this figure, the true PSD is plotted as a dashed black line and the numerical estimate of the posterior constraint on the PSD is shown as a blue contour indicating 68% of the MCMC samples. It is clear from this figure that, as expected, the inference correctly reproduces the true PSD.

{floattable}

6.2 Example 2: Inferences with the “wrong” model

In this example, we simulate a dataset using a known GP model with a kernel outside of the support of a celerite process. This means the true autocorrelation of the process can never be correctly represented by the model that we are using to fit, but we use this example to demonstrate that, at least in this case, valid inferences can still be made about the physical parameters of the model.

The data are simulated from a quasiperiodic GP with the kernel

[TABLE]

where $P_{\mathrm{true}}$ is the fundamental period of the process. This autocorrelation structure corresponds to the power spectrum (Wilson & Adams, 2013)

[TABLE]

which, for large values of $\omega$, falls off exponentially. When compared to Equation (12) – which, for large $\omega$, goes as $\omega^{-4}$ at most – it is clear that a celerite model can never exactly reproduce the structure of this process. That being said, we demonstrate that robust inferences can be made about $P_{\mathrm{true}}$ even with this effective model.

We generate a realization of a GP model with the kernel given in Equation (57) with $N=100$, $A=1\,\mathrm{ppm}^{2}$, $\lambda=5\,\mathrm{day}$, and $P_{\mathrm{true}}=1\,\mathrm{day}$. We also add white noise with amplitude $\sigma=0.5\,\mathrm{ppm}$ to each data point. The top left panel of Figure 5 shows this simulated dataset.

We then fit this simulated data using the product of two SHO terms (Equation 26) where one of the terms has $S_{0}=1$ and $Q=1/\sqrt{2}$ fixed. The kernel for this model is

[TABLE]

where $k_{\mathrm{SHO}}$ is defined in Equation (26) and the period of the process is $P=2\,\pi/\omega_{1}$. We note that using Equation (17), the product of two celerite terms can also be expressed using celerite.

We choose this functional form by comparing the Bayesian information criterion (BIC) for a set of different models. In each model, we take the sum or product of $J$ SHO terms, find the maximum likelihood model using L-BFGS-B and compute the BIC (Schwarz et al., 1978)

[TABLE]

where $\mathcal{L}^{*}$ is the value of the likelihood at maximum, $K$ is the number of parameters, and $N$ is the number of data points. Figure 6 shows the value of the BIC for a set of products and sums of SHO terms and the model that we choose has the lowest BIC.

As in the previous example, we set the mean function to zero and can, therefore, omit the parameters ${\boldsymbol{\theta}}$. Table 2 lists the proper log-uniform priors that we choose for each parameter in ${\boldsymbol{\alpha}}=(S_{0},\,Q,\,\omega_{1},\,\omega_{2})$. These priors, together with the GP likelihood (Equation 5) fully specify the posterior probability density.

As above, we estimate the MAP parameters using L-BFGS-B and sample the posterior probability density using emcee. The top left panel of Figure 5 shows the conditional mean and standard deviation of the MAP model. The bottom left panel shows the residuals between the data and this MAP model and, even though this GP model is formally “wrong”, there are no obvious correlations in these residuals.

To perform posterior inference, we initialize 32 walkers by sampling from the 4-dimensional Gaussian centered on the MAP parameters with an isotropic standard deviation of $10^{-4}$. We then run 500 steps of burn-in and 2000 steps of MCMC. To estimate the number of independent samples, we estimate the integrated autocorrelation time of the chain for the parameter $\omega_{1}$ – the parameter of primary interest – and find 1134 effective samples.

For comparison, we run the same number of steps of MCMC to sample the “correct” joint posterior density. For this reference inference, we use a GP likelihood with the kernel given by Equation (57) and choose log-uniform priors on each of the three parameters $\ln A/\mathrm{ppm}^{2}\sim\mathcal{U}(-10,\,10)$, $\ln\lambda/\mathrm{day}\sim\mathcal{U}(-10,\,10)$, and $\ln P_{\mathrm{true}}/\mathrm{day}\sim\mathcal{U}(-10,\,10)$.

The marginalized posterior inferences of the characteristic period of the process are shown in the right panel of Figure 5. The inference using the correct model is shown as a dashed blue histogram and the inference made using the effective model is shown as a solid black histogram. These inferences are consistent with each other and with the true period used for the simulation (shown as a vertical gray line). This demonstrates that, in this case, celerite can be used as a computationally efficient effective model and hints that this may be true in other problems as well.

{floattable}

6.3 Example 3: Stellar rotation

A source of variability that can be measured from time series measurements of stars is rotation. The inhomogeneous surface of the star (spots, plage, etc.) imprints itself as quasiperiodic variations in photometric or spectroscopic observations (Dumusque et al., 2014). It has been demonstrated that for light curves with nearly uniform sampling, the empirical autocorrelation function provides a reliable estimate of the rotation period of a star (McQuillan et al., 2013, 2014; Aigrain et al., 2015) and that a GP model with a quasiperiodic covariance function can be used to make probabilistic measurements even with sparsely sampled data (Angus et al., 2017). The covariance function used for this type of analysis has the form

[TABLE]

where $P_{\mathrm{rot}}$ is the rotation period of the star. GP modeling with the same kernel function has been proposed as a method of measuring the mean periodicity in quasiperiodic photometric time series in general (Wang et al., 2012). The key difference between Equation (61) and other quasiperiodic kernels is that it is positive for all values of $\tau$. We construct a simple celerite covariance function with similar properties as follows

[TABLE]

for $B>0$, $C>0$, and $L>0$. The covariance function in Equation (62) cannot exactly reproduce Equation (61) but, since Equation (61) is only an effective model, Equation (62) can be used as a drop-in replacement for a significant gain in computational efficiency.

GPs have been used to measure stellar rotation periods for individual datasets (for example Littlefair et al., 2017), but the computational cost of traditional GP methods has hindered the industrial application to existing surveys like Kepler with hundreds of thousands of targets. The increase in computational efficiency and scalability provided by celerite opens the possibility of inferring rotation periods using GPs at scale of existing and forthcoming surveys like Kepler, TESS, and LSST.

As a demonstration, we fit a celerite model with a kernel given by Equation (62) to a Kepler light curve for the star KIC 1430163. This star has a published rotation period of $3.88\pm 0.58\,\mathrm{days}$, measured using traditional periodogram and autocorrelation function approaches applied to Kepler data from Quarters 0–16 (Mathur et al., 2014), covering about four years.

We select about 180 days of contiguous observations of KIC 1430163 from Kepler. This dataset has 6950 measurements and, using a tuned linear algebra implementation666We use the Intel Math Kernel Library https://software.intel.com/en-us/intel-mkl, a single evaluation of the likelihood requires over 8 seconds on a modern Intel CPU. This calculation, using celerite with the model in Equation (62), only takes $\sim 1.5\,\mathrm{ms}$ – a speed-up of more than three orders of magnitude per model evaluation.

We set the mean function $\mu_{\boldsymbol{\theta}}$ to zero and the remaining parameters ${\boldsymbol{\alpha}}$ and their priors are listed in Table 3. As with the earlier examples, we start by estimating the MAP parameters using L-BFGS-B and initialize 32 walkers by sampling from an isotropic Gaussian with a standard deviation of $10^{-5}$ centered on the MAP parameters. The left panels of Figure 7 show a subset of the data used in this example and the residuals away from the MAP predictive mean.

We run 500 steps of burn-in, followed by 5000 steps of MCMC using emcee. We estimate the integrated autocorrelation time of the chain for $\ln P_{\mathrm{rot}}$ and estimate that we have 2900 independent samples. These samples give a posterior constraint on the period of $P_{\mathrm{rot}}={3.80\pm 0.15}$ days and the marginalized posterior distribution for $P$ is shown in the right panel of Figure 7. This result is in good agreement with the literature value with smaller uncertainties. A detailed comparison of GP rotation period measurements and the traditional methods is beyond the scope of this paper, but Angus et al. (2017) demonstrate that GP inferences are, at a population level, more reliable than other methods.

The total computational cost for this inference using celerite is about 4 CPU-minutes. By contrast, the same inference using a general but optimized Cholesky factorization routine would require nearly 400 CPU-hours. This speed-up enables probabilistic measurement of rotation periods using existing data from K2 and forthcoming surveys like TESS and LSST where this inference will need to be run for at least hundreds of thousands of stars.

{floattable}

6.4 Example 4: Asteroseismic oscillations

The asteroseismic oscillations of thousands of stars were measured using light curves from the Kepler Mission (Gilliland et al., 2010; Huber et al., 2011; Chaplin et al., 2011, 2013; Stello et al., 2013) and asteroseismology is a key science driver for many of the upcoming large scale photometric surveys (Campante et al., 2016; Rauer et al., 2014; Gould et al., 2015). Most asteroseismic analyses have been limited to high signal-to-noise oscillations because the standard methods use statistics of the empirical periodogram. These methods cannot formally propagate the measurement uncertainties to the constraints on physical parameters and they instead rely on population-level bootstrap uncertainty estimates (Huber et al., 2009). More sophisticated methods that compute the likelihood function in the time domain scale poorly to large survey datasets (Brewer & Stello, 2009; Corsaro & Ridder, 2014).

celerite alleviates these problems by providing a physically motivated probabilistic model that can be evaluated efficiently even for large datasets. In practice, we model the star as a mixture of stochastically-driven simple harmonic oscillators where the amplitudes and frequencies of the oscillations are computed using a physical model, and evaluate the probability of the observed time series using a GP where the PSD is a sum of terms given by Equation (23). This gives us a method for computing the likelihood function for the parameters of the physical model (for example, $\nu_{\mathrm{max}}$ and $\Delta\nu$, or other more fundamental parameters) conditioned on the observed time series in $\mathcal{O}(N)$ operations. In other words, celerite provides a computationally efficient framework that can be combined with physically-motivated models of stars and numerical inference methods to make rigorous probabilistic measurements of asteroseismic parameters in the time domain. This has the potential to push asteroseismic analysis to lower signal-to-noise datasets and we hope to revisit this idea in a future paper.

To demonstrate the method, we use a simple heuristic model where the PSD is given by a mixture of 8 components with amplitudes and frequencies specified by $\nu_{\mathrm{max}}$, $\Delta\nu$, and some nuisance parameters. The first term is used to capture the granulation “background” (Kallinger et al., 2014) using Equation (27) with two free parameters $S_{g}$ and $\omega_{g}$. The remaining 7 terms are given by Equation (23) where $Q$ is a nuisance parameter shared between terms and the frequencies are given by

[TABLE]

and the amplitudes are given by

[TABLE]

where $j$ is an integer running from $-3$ to 3 and $\epsilon$, $A$, and $W$ are shared nuisance parameters. Finally, we also fit for the amplitude of the white noise by adding a parameter $\sigma$ in quadrature with the uncertainties given for the Kepler observations. All of these parameters and their chosen priors are listed in Table 4. As before, these priors are all log-uniform except for $\epsilon$ where we use a zero-mean normal prior with a broad variance of $1~{}\mathrm{day}^{2}$ to break the degeneracy between $\nu_{\mathrm{max}}$ and $\epsilon$.

To build a more realistic model, this prescription could be extended to include more angular modes, or $\nu_{\mathrm{max}}$ and $\Delta\nu$ could be replaced by the fundamental physical parameters of the star.

To demonstrate the applicability of this model, we apply it to infer the asteroseismic parameters of the giant star KIC 11615890, observed by the Kepler Mission. The goal of this example is to show that, even for a low signal-to-noise dataset with a short baseline, it is possible to infer asteroseismic parameters with formal uncertainties that are consistent with the parameters inferred with a much larger dataset. Looking forward to TESS (Ricker et al., 2014; Campante et al., 2016), we measure $\nu_{\mathrm{max}}$ and $\Delta\nu$ using only one month of Kepler data and compare our results to the results inferred from the full 4 year baseline of the Kepler Mission. For KIC 11615890, the published asteroseismic parameters measured using several years of Kepler observations are (Pinsonneault et al., 2014)

[TABLE]

Unlike typical giants, KIC 11615890 is a member of a class of stars where the dipole ($\ell=1$) oscillation modes are suppressed by strong magnetic fields in the core (Stello et al., 2016). This makes this target simpler for the purposes of this demonstration because we can neglect the $\ell=1$ modes and the model proposed above will be an effective model for the combined signal from the $\ell=0$ and 2 modes.

For this demonstration, we randomly select a month-long segment of PDC Kepler data (Stumpe et al., 2012; Smith et al., 2012). Unlike the previous examples, the posterior distribution is sharply multimodal and naïvely maximizing the posterior using L-BFGS-B is not practical. Instead, we start by estimating the initial values for $\nu_{\mathrm{max}}$ and $\Delta\nu$ using only the month-long subset of data and following the standard procedure described by Huber et al. (2009). We then run a set of L-BFGS-B optimizations with values of $\nu_{\mathrm{max}}$ selected in logarithmic grid centered on our initial estimate of $\nu_{\mathrm{max}}$ and initial values of $\Delta\nu$ computed using the empirical relationship between these two quantities (Stello et al., 2009). This initialization procedure is sufficient for this example, but the general application of this method will require a more sophisticated prescription.

Figure 8 shows a 10 day subset of the dataset used for this example. The MAP model is overplotted on these data and the residuals away from the mean prediction of this model are shown in the bottom panel of Figure 8. There is no obvious structure in these residuals, lending some credibility to the model specification.

We initialize 32 walkers by sampling from an isotropic Gaussian centered on the MAP parameters (the full set of parameters and their priors are listed in Table 4), run 5000 steps of burn-in, and run 15000 steps of MCMC using emcee. We estimate the mean autocorrelation time for the chains of $\ln\nu_{\mathrm{max}}$ and $\ln\Delta\nu$ and find 1443 effective samples from the marginalized posterior density. Figure 9 shows the marginalized density for $\nu_{\mathrm{max}}$ and $\Delta\nu$ compared to the results from the literature. This result is consistent within the published error bars and the posterior constraints are tighter than the published results. The top two panels of Figure 10 show the Lomb-Scargle periodogram (VanderPlas, 2017) estimator for the power spectrum based on the full 4 years of Kepler data and the month-long subset used in our analysis. The bottom panel of Figure 10 shows the posterior inference of the PSD using the celerite and the model described here. Despite only using one month of data, the celerite inference captures the salient features of the power spectrum and it is qualitatively consistent with the standard estimator applied to the full baseline. All asteroseismic analyses are known to have substantial systematic and method-dependent uncertainties (Verner et al., 2011) so further experiments would be needed to fully assess the reliability of this specific model.

This model requires about 10 CPU minutes to run the MCMC to convergence. This is more computationally intensive than traditional methods of measuring asteroseismic oscillations, but it is much cheaper than the same analysis using a general direct GP solver. For comparison, we estimate that repeating this analysis using a general Cholesky factorization implemented as part of a tuned linear algebra library777We use the Intel Math Kernel Library https://software.intel.com/en-us/intel-mkl would require about 15 CPU hours. An in-depth discussion of the benefits of rigorous probabilistic inference of asteroseismic parameters in the time domain is beyond the scope of this paper, but we hope to revisit this opportunity in the future.

{floattable}

6.5 Example 5: Exoplanet transit fitting

This example is different from all the previous examples – both simulated and real – because, in this case, do not set the deterministic mean function $\mu_{\boldsymbol{\theta}}$ to zero. Instead, we make inferences about $\mu_{\boldsymbol{\theta}}$ because it is a physically significant function and the parameters are fundamental properties of the system. This is an example using real data – because we use the light curve from Section 6.3 – but we multiply these data by a simulated transiting exoplanet model with known parameters. This allows us to show that we can recover the true parameters of the planet even when the transit signal is superimposed on the real variability of a star.

GP modeling has been used extensively for this purpose throughout the exoplanet literature (for example Dawson et al., 2014; Barclay et al., 2015; Evans et al., 2015; Foreman-Mackey et al., 2016; Grunblatt et al., 2016).

In Equation (5) the physical parameters of the exoplanet are called ${\boldsymbol{\theta}}$ and, in this example, the mean function $\mu_{\boldsymbol{\theta}}(t)$ is a limb-darkened transit light curve (Mandel & Agol, 2002; Foreman-Mackey & Morton, 2016) and the parameters ${\boldsymbol{\theta}}$ are the orbital period $P_{\mathrm{orb}}$, the transit duration $T$, the phase or epoch $t_{0}$, the impact parameter $b$, the radius of the planet in units of the stellar radius $R_{P}/R_{\star}$, the baseline relative flux of the light curve $f_{0}$, and two parameters describing the limb-darkening profile of the star (Claret & Bloemen, 2011; Kipping, 2013). As in Section 6.3, we model the stellar variability using a GP model with a kernel given by Equation (62) and fit for the parameters of the exoplanet ${\boldsymbol{\theta}}$ and the stellar variability ${\boldsymbol{\alpha}}$ simultaneously. The full set of parameters ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\theta}}$ are listed in Table 5 along with their priors and the true values for ${\boldsymbol{\theta}}$.

We take a 20 day segment of the Kepler light curve for KIC 1430163 ($N=1000$) and multiply it by a simulated transit model with the parameters listed in Table 5. Using these data, we maximize the joint posterior defined by the likelihood in Equation (5) and the priors in Table 5 using L-BFGS-B for all the parameters ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\theta}}$ simultaneously. The top panel of Figure 11 shows the data including the simulated transit as black points with the MAP model prediction over-plotted in blue. The bottom panel of Figure 11 shows the “de-trended” light curve where the MAP model has been subtracted and the MAP mean model $\mu_{\boldsymbol{\theta}}$ has been added back in. For comparison, the transit model is over-plotted in the bottom panel of Figure 11 and we see no obvious correlations in the residuals.

Sampling 32 walkers from an isotropic Gaussian centered on the MAP parameters, we run 10000 steps of burn-in and 30000 steps of production MCMC using emcee. We estimate the integrated autocorrelation time for the $\ln(P_{\mathrm{orb}})$ chain and find 1490 effective samples across the full chain. Figure 12 shows the marginalized posterior constraints on the physical properties of the planet compared to the true values. Even though the celerite representation of the stellar variability is only an effective model, the inferred distributions for the physical properties of the planet ${\boldsymbol{\theta}}$ are consistent with the true values. This promising result suggests that celerite can be used as an effective model for transit inference for a significant gain in computational tractability. Even for this example with small $N$, the cost of computing the likelihood using celerite is nearly two orders of magnitude faster than the same computation using a general direct solver.

{floattable}

6.6 Summary

In the previous subsections, we demonstrate five potential use cases for celerite. These examples span a range of data sizes and model complexities, and the last three are based on active areas of research with immediate real-world applicability. The sizes of these datasets are modest on the scale of some existing and forthcoming survey datasets because we designed the experiments to be easily reproducible but, even in these cases, the inferences made here would be intractable without substantial computational investment.

Table 6 lists the specifications of each example. In this table, we list the computational cost required for a single likelihood evaluation using a general Cholesky factorization and the cost of computing the same likelihood using celerite. Combining this with the total number of model evaluations required to run the MCMC algorithm to convergence (this number is also listed in Table 6), we can estimate the total CPU time required in each case. While the exact computational requirements for any problem also depend on the choice of inference algorithm and the specific precision requirements, this table provides a rough estimate of what to expect for projects at these scales.

{floattable}

Section 7 Comparisons to other methods

There are other methods of scaling GP models to large datasets and in this section we draw comparisons between celerite and other popular methods. Scalable GP methods tend to fall into two categories: approximate and restrictive. celerite falls into the latter category because, while the method is exact, it requires a specific choice of stationary kernel function and it can only be used in one-dimension.

The continuous autoregressive moving average (CARMA) models introduced to the astrophysics literature by Kelly et al. (2014) share many features with celerite. CARMA models are restricted to one-dimensional problems and the likelihood function for a CARMA model can be solved in $\mathcal{O}(N\,J^{2})$ using a Kalman filter (Kelly et al., 2014). The kernel function for a CARMA$(J,\,K)$ model is

[TABLE]

where

[TABLE]

and $\sigma$, $\{r_{j}\}_{j=1}^{J}$ and $\{\beta_{k}\}_{k=1}^{K}$ are parameters of the model; $\sigma$ is real, while the others are complex. Comparing Equation (66) to Equation (10), we can see that every CARMA model corresponds to an equivalent celerite model and the parameters $a_{j}$, $b_{j}$, $c_{j}$, and $d_{j}$ can be easily computed analytically. The inverse statement is not as simple. In practice, this means that celerite could be trivially used to compute any CARMA model. Using the CARMA solver to compute a celerite model, however, requires solving Equation (67) numerically for a given set of $\{A_{j}\}_{j=1}^{J}$.

The computational scaling of CARMA models is also $\mathcal{O}(N\,J^{2})$ using the Kalman filter method (Kelly et al., 2014), but we find that, in practice, this method is more computationally expensive than celerite, but it has a smaller memory footprint. Figure 13 shows the scaling of the Kalman filter solver for the systems shown in Figure 3. Comparing these figures, we find that the Kalman filter solver is slower than celerite for all the systems we tested with an average difference of about an order of magnitude.

Another benefit of celerite compared to Kalman filter solvers is the fact that it is an algebraic solver for the relevant matrix equations instead of an algorithm to directly compute Equation (5). This means that the inverse covariance matrix $K^{-1}$ can be applied to general vectors and matrices in $\mathcal{O}(N)$, while reusing the factorization. This is a crucial feature when GP models are combined with linear regression (for example Luger et al., 2017).

Another popular GP method uses the fact that, in the limit of evenly-spaced data and homoscedastic uncertainties, the covariance matrix is “Toeplitz” (for example Dillon et al., 2013). There are exact methods for solving Toeplitz matrix equations that scale as $\mathcal{O}(N\,\log N)$ and methods for computing determinants exactly in $\mathcal{O}(N^{2})$ or approximately in $\mathcal{O}(N\,\log N)$ (Wilson, 2014). The Toeplitz method is, in some ways, more flexible than celerite because it can be used with any stationary kernel, but it requires uniformly spaced data and the scaling is worse than celerite so it is less efficient when applied to large datasets.

Carter & Winn (2009) improved the scaling of Toeplitz methods by introducing a wavelet-based method for computing a GP likelihood with $\mathcal{O}(N)$ scaling. This method has been widely applied in the context of exoplanet transit characterization, but it requires evenly spaced observations and the power spectrum of the process must have the form $S(\omega)\propto\omega^{-1}$ to gain the computational advantage. This wavelet method has been demonstrated to improve parameter estimation for transiting exoplanets (Carter & Winn, 2009), but these strict requirements make this method applicable for only a limited set of use cases.

Another GP method that has been used extensively in astronomy is the hierarchical off-diagonal low rank (HODLR) solver (Ambikasaran et al., 2016). This method exploits the fact that many commonly used kernel functions produce “smooth” matrices to approximately compute the GP likelihood with the scaling $\mathcal{O}(N\,\log^{2}N)$. This method has the advantage that, unlike celerite, it can be used with any kernel function but, in practice, the cost can still prove to be prohibitively high for multi-dimensional inputs. The proportionality constant in the $N\,\log^{2}N$ scaling of the HODLR method is a function of the specific kernel and we find – using the george software package (Foreman-Mackey et al., 2014; Ambikasaran et al., 2016) – that this scales approximately linearly with $J$, but it requires substantial overhead making it much slower than celerite for all the models we tested. For large $J\gtrsim 256$ and small $N\lesssim 1000$, we find that george can approximately evaluate celerite models with the same or less computational cost than celerite, but that celerite is faster – and exact – in all other parts of parameter space.

The structured kernel interpolation (SKI/KISS-GP Wilson & Nicisch, 2015) framework is another approximate method that can be used to scale GPs to large datasets. This method uses fast interpolation of the kernel into a space where Toeplitz or Kronecker structure can be used to scale inference and prediction. The SKI/KISS-GP framework can be applied to scale GP inference with a wide range of kernels, but the computational cost will depend on the specific dataset, model, and precision requirements for the approximation. The SKI method has an impressive $\mathcal{O}(1)$ cost for test-time predictions and it is interesting to consider how this could be applied to celerite models.

Many other approximate methods for scaling GP inference exist (see, for example, Wilson et al., 2015, and references therein) and we make no attempt to make our discussion exhaustive. The key takeaway here is that celerite provides an exact method for GP inference for a specific class of one-dimensional kernel functions. Furthermore, since celerite models can be interpreted as a mixture of stochastically-driven, damped simple harmonic oscillators, they are a physically motivated choice of covariance function in many astronomical applications.

Section 8 Discussion

Gaussian Process models have been fruitfully applied to many problems in astronomical data analysis, but the fact that the computational cost generally scales as the cube of the number of data points has limited their use to small datasets with $N\lesssim 1000$. With the linear scaling of celerite  we envision that the application of Gaussian processes will be expanded to the existing and forthcoming large astronomical time domain surveys such as Kepler (Borucki et al., 2010), K2 (Howell et al., 2014), TESS (Ricker et al., 2014), LSST (Ivezić et al., 2008), WFIRST (Spergel et al., 2015), and PLATO (Rauer et al., 2014). Despite the restrictive form of the celerite kernel, with a sufficient number of components it is flexible enough to describe a wide range of astrophysical variability. In fact, the relation of the celerite kernel to the damped, stochastically-driven harmonic oscillator matches simple models of astrophysical variability, and makes the parameterization interpretable in terms of resonant frequency, amplitude, and quality factor.

Our background is in studying transiting exoplanets, a field which has only recently begun to adopt GP methods for analyzing the noise in stellar light curves and radial velocity datasets when detecting or characterizing transiting planets (for example, Carter & Winn, 2009; Gibson et al., 2012; Haywood et al., 2014; Barclay et al., 2015; Evans et al., 2015; Rajpaul et al., 2015; Aigrain et al., 2016; Foreman-Mackey et al., 2016; Grunblatt et al., 2016; Luger et al., 2016). All of these analyses have been limited to small datasets or restrictive kernel choices, but celerite weakens these requirements by providing a scalable method for computing the likelihood and a physical motivation for the choice of kernel. As higher signal-to-noise observations of transiting exoplanet systems are obtained, the effects of stellar variability will more dramatically impact the correct inference of planetary transit parameters. We predict that celerite will be important for scaling methods of transit detection (Pope et al., 2016; Foreman-Mackey et al., 2016), transit timing (Agol et al., 2005; Holman, 2005), transit spectroscopy (Brown, 2001), Doppler beaming (Loeb & Gaudi, 2003; Zucker et al., 2007), tidal distortion (Zucker et al., 2007), and phase functions (Knutson et al., 2007; Zucker et al., 2007) to the large astronomical time domain surveys of the future.

8.1 Other applications and limitations

Beyond these applications to model stellar variability, the method is generally applicable to other one-dimensional GP models. Accreting black holes show time series which may be modeled using a GP (Kelly et al., 2014); indeed, this was the motivation for the original technique developed by Rybicki & Press (Rybicki & Press, 1992, 1995). This approach may be broadly used for characterizing quasar variability (MacLeod et al., 2010), measuring time lags with reverberation mapping (Zu et al., 2011; Pancoast et al., 2014), modeling time delays in multiply-imaged gravitationally-lensed systems (Press & Rybicki, 1998), characterizing quasi-periodic variability in a high-energy source (McAllister et al., 2016), or classification of variable objects (Zinn et al., 2016). We expect that there are also applications beyond astronomy.

The celerite formalism can also be used for power spectrum estimation and quantification of its uncertainties. A mixture of celerite terms can be used to perform non-parametric probabilistic inference of the power spectrum despite unevenly-spaced data with heteroscedastic noise (see Wilson & Adams, 2013; Kelly et al., 2014, for examples of this procedure). This type of analysis will be limited by the quadratic scaling of celerite with the number of terms $J$, but this limits existing methods as well (Kelly et al., 2014).

There are many data analysis problems where celerite will not be immediately applicable. In particular, the restriction to one-dimensional problems is significant. There are many examples of multidimensional GP modeling the astrophysics literature (recent examples from the field of exoplanet characterization include Haywood et al., 2014; Rajpaul et al., 2015; Aigrain et al., 2016), where celerite cannot be used to speed up the analysis. It is plausible that an extension can be derived to tackle some multidimensional problems with the right structure – simultaneous parallel time series, for example – and we hope to revisit this possibility in future work.

8.2 Code availability

Alongside this paper, we have released a well-tested and documented open source software package that implements the method and all of the examples discussed in these pages. This software is available on GitHub https://github.com/dfm/celerite888This version of the paper was generated with git commit fd40d6a8 (2017-07-19). and Zenodo (Foreman-Mackey et al., 2017), and it is made available under the MIT license.

It is a pleasure to thank Megan Bedell, Ian Czekala, Will Farr, Sam Grunblatt, David W. Hogg, Dan Huber, Meredith Rawls, Dennis Stello, Jake VanderPlas, and Andrew Gordon Wilson for helpful discussions informing the ideas and code presented here. We would also like to thank the anonymous referee for thorough and constructive feedback that greatly improved the paper. This work was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. EA acknowledges support from NASA grants NNX13AF20G, NNX13A124G, NNX13AF62G, from National Science Foundation (NSF) grant AST-1615315, and from NASA Astrobiology Institute’s Virtual Planetary Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C. This research made use of the NASA Astrophysics Data System and the NASA Exoplanet Archive. The Exoplanet Archive is operated by the California Institute of Technology, under contract with NASA under the Exoplanet Exploration Program. This paper includes data collected by the Kepler Mission. Funding for the Kepler Mission is provided by the NASA Science Mission directorate. We are grateful to the entire Kepler team, past and present. These data were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013). Kepler

Appendix A Ensuring positive definiteness for celerite models

For a GP kernel to be valid, it must produce a positive definite covariance matrix for all input coordinates. For stationary kernels, this is equivalent – by Bochner’s theorem (see Section 4.2.1 in Rasmussen & Williams, 2006) – to requiring that the kernel be the Fourier transform of a positive finite measure. This means that the power spectrum of a positive definite kernel must be positive for all frequencies. This result is intuitive because, since the power spectrum of a process is defined as the expected squared amplitude of the Fourier transform of the time series, it must be non-negative.

Using Equation (12), we find that for a single celerite term, this requirement is met when

[TABLE]

The denominator is positive for all $c_{j}\neq 0$ and it can be shown that, when $c_{j}=0$, Equation (A1) is satisfied for all $\omega\neq d_{j}$ where the power is identically zero. Therefore, when $c_{j}\neq 0$, we require that the numerator is positive for all $\omega$. This requirement can also be written as

[TABLE]

Furthermore, we can see that $a_{j}$ must be positive since $k(0)=a_{j}$ should should be positive and, similarly, by requiring the covariance to be finite at infinite lag, we obtain the constraint $c_{j}\geq 0$. Combining these results, we find the constraint

[TABLE]

The constraint for $J>1$ is more complicated so, in most cases, we require that each term is positive definite using this relationship and use the fact that a product of sum of positive definite kernels will also be positive definite (Rasmussen & Williams, 2006) to show that the full model is positive. However, in the case of $J$ general celerite terms, we can check for negative values of the PSD by solving for the roots of the power spectrum; if there are any real, positive roots, then the power-spectrum goes negative (or zero), and thus does not represent a valid kernel. We rewrite the power spectrum, Equation 12), abbreviating with $z=\omega^{2}$:

[TABLE]

where

[TABLE]

The denominators of each term are positive, so we can multiply through by $\prod_{j}\left(z^{2}+s_{j}z+t_{j}\right)$ to find

[TABLE]

which is a polynomial with order $2\,(J-1)+1$. With $J=2$, this yields a cubic equation whose roots can be obtained exactly.

For arbitrary $J$, a procedure based upon Sturm’s theorem (Dörrie, 1965) allows one to determine whether there are any real roots within the range $(0,\,\infty]$. We first construct $Q_{0}(z)$ and its derivative $Q_{1}(z)={Q_{0}}^{\prime}(z)$, and then loop from $k=2$ to $k=2\,(J-1)+1$, computing

[TABLE]

where the function ${\rm rem}(p,\,q)$ is the remainder polynomial after dividing $p(z)$ by $q(z)$.

We evaluate the coefficients of each of the polynomials in the series by evaluating $f_{0}=\{Q_{0}(0),\,\ldots,\,Q_{2\,(J-1)+1}(0)\}$ to give us the signs of these polynomials evaluated at $z=0$. Likewise, we evaluate the coefficients of the largest order term in each polynomial that gives the sign of the polynomial as $z\rightarrow\infty$, $f_{\infty}$. With the sequence of coefficients $f_{0}$ and $f_{\infty}$, we then determine how many times the sign changes in each of these, where $\sigma(0)$ is the number of sign changes at $z=0$, and $\sigma(\infty)$ is the number of sign changes at $z\rightarrow\infty$. The total number of real roots in the range $(0,\,\infty]$ is given by $N_{+}=\sigma(0)-\sigma(\infty)$.

We have checked that this procedure works for a wide range of parameters, and we find that it robustly matches the number of positive real roots which we evaluated numerically. The advantage of this procedure is that it does not require computing the roots, but only carrying out algebraic manipulation of polynomials to determine the number of positive real roots. If a non-zero real root is found, the likelihood may be set to zero.

Appendix B Semiseparable matrix operations

In this appendix, we summarize some algorithms for manipulating semiseparable matrices and discuss how these can be used with GP models and celerite.

B.1 Multiplication

The product of a rank-$R$ semiseparable matrix with general vectors and matrices can be computed in $\mathcal{O}(N\,R)$ operations (Vandebril et al., 2007). In the context of celerite models, naïve application of these methods will result in numerical overflow and underflow issues, much like we found with the Cholesky factorization. To avoid these errors, we derive the following numerically stable algorithm for computing the dot product of a celerite covariance matrix

[TABLE]

using the semiseparable representation for $K$ defined in Equation 5.2:

[TABLE]

where $f^{+}_{N,j}=0$ and $f^{-}_{1,j}=0$. This algorithm requires two sweeps – $n=2,\,\ldots,\,N$ to compute $f^{-}$ and $n=N-1,\,\ldots,\,1$ to compute $f^{+}$ – with a scaling of $\mathcal{O}(N\,J)$.

B.2 Sampling data from a celerite process

To sample a dataset ${\boldsymbol{y}}$ from a GP model for fixed parameters ${\boldsymbol{\theta}}$ and ${\boldsymbol{\alpha}}$, we must compute

[TABLE]

where ${\boldsymbol{\mu}}_{\boldsymbol{\theta}}({\boldsymbol{t}})$ is the mean function evaluated at the input coordinates ${\boldsymbol{t}}$, ${\boldsymbol{q}}$ is a vector of draws from the unit normal $q_{i}\sim\mathcal{N}(0,\,1)$, and

[TABLE]

For a general matrix $K_{\boldsymbol{\alpha}}$, the cost of this operation scales as $\mathcal{O}(N^{3})$ to compute $K_{\boldsymbol{\alpha}}^{1/2}$ and $\mathcal{O}(N^{2})$ to compute the dot product. For celerite models, we use the semiseparable representation of the Cholesky factor that we derived in Section 5.1 to compute the dot product Equation (B5) in $\mathcal{O}(N\,J)$. Using the notation from Section 5.2, the algorithm to compute this dot product is

[TABLE]

where $f_{1,j}=0$ for all $j$. This scalable algorithm is useful for generating large simulated datasets for experiments with celerite and for performing posterior predictive checks.

B.3 Interpolation & extrapolation

The predictive distribution of a GP at a vector of $M$ input coordinates

[TABLE]

conditioned on a dataset $({\boldsymbol{t}},\,{\boldsymbol{y}})$ and GP parameters $({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}})$, is a normal ${\boldsymbol{y}}^{*}\sim\mathcal{N}({\boldsymbol{\mu}}^{*},\,K^{*})$ with mean

[TABLE]

and covariance

[TABLE]

where $K({\boldsymbol{v}},\,{\boldsymbol{w}})$ is the covariance matrix computed between $\boldsymbol{v}$ and $\boldsymbol{w}$ (Rasmussen & Williams, 2006).

Naïvely, the computational cost of computing Equation (B11) for a celerite model scales as $\mathcal{O}(N\,M)$ if we reuse the Cholesky factor derived in Section 5.1. It is, however, possible to improve this scaling to $\mathcal{O}(n\,N+m\,M)$ where $n$ and $m$ are integer constants. To derive this, we expand Equation (B11) using Equation (11)

[TABLE]

where

[TABLE]

Dividing the sum over $n$ into $\{t_{1},...,t_{n_{0}}\}<t^{*}_{m}$ and $\{t_{n0+1},...,t_{N}\}\geq t^{*}_{m}$, gives

[TABLE]

We compute this using two passes, a forward pass $n_{0}=1,\,\ldots,\,N$ and a backward pass $n_{0}=N,\,\ldots,\,1$. Defining

[TABLE]

where $t_{0}=t_{1}$, $t_{N+1}=t_{N}$, $Q^{-}_{0,k}=0$ and $Q^{+}_{N+1,k}=0$ for $k=1$ to $2\,J$, and in $X^{\pm}_{m,n,k}$ the expressions for $\tilde{U}_{m,i}^{*}$ and $\tilde{V}_{m,i}^{*}$ are evaluated at $t_{m}^{*}$. For each value of $m$, $Q^{\pm}$ are recursively updated over $n$ until $n_{0}$ is reached, at which point the predicted mean can be computed from

[TABLE]