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

Daniel Foreman-Mackey11affiliation: Sagan Fellow 22affiliation: Astronomy
Department, University of Washington, Seattle, WA , Eric Agol22affiliation:
Astronomy Department, University of Washington, Seattle, WA 33affiliation:
Guggenheim Fellow , Sivaram Ambikasaran44affiliation: Department of
Computational and Data Sciences, Indian Institute of Science, Bangalore, India
, and Ruth Angus55affiliation: Simons Fellow 66affiliation: Department of
Astronomy, Columbia University, New York, NY

###### 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.

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

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

at coordinates

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

is

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

where

$\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}}$
(7)

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

$\displaystyle
k_{\boldsymbol{\alpha}}(\tau_{nm})=\sigma_{n}^{2}\,\delta_{nm}+a\,\exp(-c\,\tau_{nm})$
(8)

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)

$\displaystyle
k_{\boldsymbol{\alpha}}(\tau_{nm})=\sigma_{n}^{2}\,\delta_{nm}+\sum_{j=1}^{J}a_{j}\,\exp(-c_{j}\,\tau_{nm})\quad.$
(9)

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

$\displaystyle
k_{\boldsymbol{\alpha}}(\tau_{nm})=\sigma_{n}^{2}\,\delta_{nm}+\sum_{j=1}^{J}$
$\displaystyle\left[\frac{1}{2}(a_{j}+i\,b_{j})\,\exp\left(-(c_{j}+i\,d_{j})\,\tau_{nm}\right)\right.$
(10)
$\displaystyle+\left.\frac{1}{2}(a_{j}-i\,b_{j})\,\exp\left(-(c_{j}-i\,d_{j})\,\tau_{nm}\right)\right]$

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

$\displaystyle
k_{\boldsymbol{\alpha}}(\tau_{nm})=\sigma_{n}^{2}\,\delta_{nm}+\sum_{j=1}^{J}$
$\displaystyle\left[a_{j}\,\exp\left(-c_{j}\,\tau_{nm}\right)\,\cos\left(d_{j}\,\tau_{nm}\right)\right.$
(11)
$\displaystyle+\left.b_{j}\,\exp\left(-c_{j}\,\tau_{nm}\right)\,\sin\left(d_{j}\,\tau_{nm}\right)\right]\quad.$

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

$\displaystyle
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}}\quad.$
(12)

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

$\displaystyle
k_{j}(\tau_{nm})=a_{j}\,\exp\left(-c_{j}\,\tau_{nm}\right)\,\cos\left(d_{j}\,\tau_{nm}\right)$
(13)

and the PSD for this component is

$\displaystyle
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]\quad.$
(14)

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

$\displaystyle k_{j}(\tau_{nm})=a_{j}\,\exp\left(-c_{j}\,\tau_{nm}\right)$
(15)

with the PSD

$\displaystyle
S_{j}(\omega)=\sqrt{\frac{2}{\pi}}\,\frac{a_{j}}{c_{j}}\,\frac{1}{1+{\left(\frac{\omega}{c_{j}}\right)}^{2}}\quad.$
(16)

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

$\displaystyle
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)]$
(17)

where

$\displaystyle\tilde{a}_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,(a_{j}\,a_{k}\pm b_{j}\,b_{k})$ (18)
$\displaystyle\tilde{b}_{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,(b_{j}\,a_{k}\mp a_{j}\,b_{k})$ (19)
$\displaystyle\tilde{c}$ $\displaystyle=$ $\displaystyle c_{j}+c_{k}$ (20)
$\displaystyle\tilde{d}_{\pm}$ $\displaystyle=$ $\displaystyle d_{j}\mp
d_{k}\quad.$ (21)

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

$\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)=\epsilon(t)$
(22)

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)

$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}}$
(23)

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

$\displaystyle a_{j}$ $\displaystyle=$ $\displaystyle S_{0}\,\omega_{0}\,Q$
(24) $\displaystyle b_{j}$ $\displaystyle=$
$\displaystyle\frac{S_{0}\,\omega_{0}\,Q}{\sqrt{4\,Q^{2}-1}}$ $\displaystyle
c_{j}$ $\displaystyle=$ $\displaystyle\frac{\omega_{0}}{2\,Q}$ $\displaystyle
d_{j}$ $\displaystyle=$
$\displaystyle\frac{\omega_{0}}{2\,Q}\sqrt{4\,Q^{2}-1}\quad,$

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

$\displaystyle a_{j\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,S_{0}\,\omega_{0}\,Q\,\left[1\pm\frac{1}{\sqrt{1-4\,Q^{2}}}\right]$
(25) $\displaystyle b_{j\pm}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle c_{j\pm}$ $\displaystyle=$
$\displaystyle\frac{\omega_{0}}{2\,Q}\,\left[1\mp\sqrt{1-4\,Q^{2}}\right]$
$\displaystyle d_{j\pm}$ $\displaystyle=$ $\displaystyle 0\quad.$

For these definitions, the kernel is

$k_{\mathrm{SHO}}(\tau;\,S_{0},\,Q,\,\omega_{0})=S_{0}\,\omega_{0}\,Q\,e^{-\frac{\omega_{0}\,\tau}{2Q}}\,\begin{cases}\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{cases}$ (26)

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

$\displaystyle
S(\omega)=\sqrt{\frac{2}{\pi}}\,\frac{S_{0}}{(\omega/\omega_{0})^{4}+1}\quad.$
(27)

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

$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)}\quad.$
(28)

  * •

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

$\displaystyle
S(\omega)=\sqrt{\frac{2}{\pi}}\,\frac{S_{0}}{\left[(\omega/\omega_{0})^{2}+1\right]^{2}}$
(29)

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

$\displaystyle k(\tau)$ $\displaystyle=$ $\displaystyle\lim_{f\to
0}\,\frac{1}{2}\,S_{0}\,\omega_{0}\,\left[\left(1+1/f\right)\,e^{-\omega_{0}\,(1-f)\,\tau}+\left(1-1/f\right)\,e^{-\omega_{0}\,(1+f)\,\tau}\right]$
(30) $\displaystyle=$ $\displaystyle
S_{0}\,\omega_{0}\,e^{-\omega_{0}\,\tau}\,[1+\omega_{0}\,\tau]$ (31)

or, equivalently (using Equation 11 and Equation 24)

$\displaystyle k(\tau)$ $\displaystyle=$ $\displaystyle\lim_{f\to
0}\,S_{0}\,\omega_{0}\,e^{-\omega_{0}\,\tau}\,\left[\cos(f\,\tau)+\frac{\omega_{0}}{f}\,\sin(f\,\tau)\right]$
(32) $\displaystyle=$ $\displaystyle
S_{0}\,\omega_{0}\,e^{-\omega_{0}\,\tau}\,[1+\omega_{0}\,\tau]\quad.$ (33)

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

$\displaystyle k(\tau)\approx
S_{0}\,\omega_{0}\,Q\,\exp\left(-\frac{\omega_{0}\,\tau}{2\,Q}\right)\,\cos\left(\omega_{0}\,\tau\right)\quad.$
(34)

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.

Figure 1: (left) The power spectrum of a stochastically-driven simple harmonic
oscillator (Equation 23) plotted for several values of the quality factor $Q$.
For comparison, the dashed line shows the Lorentzian function from Equation
(14) with $c_{j}=\omega_{0}/(2\,Q)=1/20$ and normalized so that
$S(d_{j})/S(0)=100$. (middle) The corresponding autocorrelation functions with
the same colors. (right) Three realizations in arbitrary units from each model
in the same colors.

##  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

$\displaystyle
K_{n,m}=\left\\{\begin{array}[]{ll}\sum_{r=1}^{R}U_{n,r}\,V_{m,r}&\mathrm{if}\,m<n\\\
A_{n,n}&\mathrm{if}\,m=n\\\
\sum_{r=1}^{R}U_{m,r}\,V_{n,r}&\mathrm{otherwise}\end{array}\right.$ (38)

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

$\displaystyle
K=A+\mathrm{tril}(U\,V^{\mathrm{T}})+\mathrm{triu}(V\,U^{\mathrm{T}})$ (39)

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$:

$\displaystyle K=L\,D\,L^{\mathrm{T}}$ (40)

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

$\displaystyle L=I+\mathrm{tril}(U\,W^{T})$ (41)

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

$\displaystyle
A+\mathrm{tril}(U\,V^{\mathrm{T}})+\mathrm{triu}(V\,U^{\mathrm{T}})$
$\displaystyle=$
$\displaystyle\left[I+\mathrm{tril}(U\,W^{T})\right]\,D\,{\left[I+\mathrm{tril}(U\,W^{T})\right]}^{\mathrm{T}}$
$\displaystyle=$ $\displaystyle
D+\mathrm{tril}(U\,W^{\mathrm{T}})\,D+D\,\mathrm{triu}(W\,U^{\mathrm{T}})$
$\displaystyle+\,\mathrm{tril}(U\,W^{\mathrm{T}})\,D\,\mathrm{triu}(W\,U^{\mathrm{T}})$

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

$\displaystyle S_{n,j,k}$ $\displaystyle=$ $\displaystyle
S_{n-1,j,k}+D_{n-1,n-1}\,W_{n-1,j}\,W_{n-1,k}$ $\displaystyle D_{n,n}$
$\displaystyle=$ $\displaystyle
A_{n,n}-\sum_{j=1}^{R}\sum_{k=1}^{R}U_{n,j}\,S_{n,j,k}\,U_{n,k}$
$\displaystyle W_{n,j}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{n,n}}\left[V_{n,j}-\sum_{k=1}^{R}U_{n,k}\,S_{n,j,k}\right]$
(43)

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

$\displaystyle\ln\det K=\sum_{n=1}^{N}\ln D_{n,n}$ (44)

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

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

$\displaystyle{\boldsymbol{z}}=K^{-1}\,{\boldsymbol{y}}={(L^{\mathrm{T}})}^{-1}\,D^{-1}\,L^{-1}\,{\boldsymbol{y}}\quad.$
(45)

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

$\displaystyle f_{n,j}$ $\displaystyle=$ $\displaystyle
f_{n-1,j}+W_{n-1,j}\,z_{n-1}^{\prime}$ $\displaystyle z_{n}^{\prime}$
$\displaystyle=$ $\displaystyle y_{n}-\sum_{j=1}^{R}U_{n,j}\,f_{n,j}$ (46)

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

$\displaystyle g_{n,j}$ $\displaystyle=$ $\displaystyle
g_{n+1,j}+U_{n+1,j}\,z_{n+1}$ $\displaystyle z_{n}$ $\displaystyle=$
$\displaystyle\frac{z_{n}^{\prime}}{D_{n,n}}-\sum_{j=1}^{R}W_{n,j}\,g_{n,j}$
(47)

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

$\displaystyle U_{n,2\,j-1}$ $\displaystyle=$ $\displaystyle
a_{j}\,e^{-c_{j}\,t_{n}}\,\cos(d_{j}\,t_{n})+b_{j}\,e^{-c_{j}\,t_{n}}\,\sin(d_{j}\,t_{n})$
$\displaystyle U_{n,2\,j}$ $\displaystyle=$ $\displaystyle
a_{j}\,e^{-c_{j}\,t_{n}}\,\sin(d_{j}\,t_{n})-b_{j}\,e^{-c_{j}\,t_{n}}\,\cos(d_{j}\,t_{n})$
$\displaystyle V_{m,2\,j-1}$ $\displaystyle=$ $\displaystyle
e^{c_{j}\,t_{m}}\,\cos(d_{j}\,t_{m})$ $\displaystyle V_{m,2\,j}$
$\displaystyle=$ $\displaystyle e^{c_{j}\,t_{m}}\,\sin(d_{j}\,t_{m})$
$\displaystyle A_{n,n}$ $\displaystyle=$
$\displaystyle\sigma_{n}^{2}+\sum_{j=1}^{J}a_{j}\quad.$ (48)

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

$\displaystyle\tilde{U}_{n,2\,j-1}$ $\displaystyle=$ $\displaystyle
a_{j}\,\cos(d_{j}\,t_{n})+b_{j}\,\sin(d_{j}\,t_{n})$
$\displaystyle\tilde{U}_{n,2\,j}$ $\displaystyle=$ $\displaystyle
a_{j}\,\sin(d_{j}\,t_{n})-b_{j}\,\cos(d_{j}\,t_{n})$
$\displaystyle\tilde{V}_{m,2\,j-1}$ $\displaystyle=$
$\displaystyle\cos(d_{j}\,t_{m})$ $\displaystyle\tilde{V}_{m,2\,j}$
$\displaystyle=$ $\displaystyle\sin(d_{j}\,t_{m})$
$\displaystyle\tilde{W}_{n,2\,j-1}$ $\displaystyle=$ $\displaystyle
e^{-c_{j}\,t_{n}}\,W_{n,2\,j-1}$ $\displaystyle\tilde{W}_{n,2\,j}$
$\displaystyle=$ $\displaystyle e^{-c_{j}\,t_{n}}\,W_{n,2\,j}$ (49)

and the new set of variables

$\displaystyle\phi_{n,2\,j-1}=\phi_{n,2\,j}=e^{-c_{j}\,(t_{n}-t_{n-1})}$ (50)

with the constraint

$\displaystyle\phi_{1,2\,j-1}=\phi_{1,2\,j}=0\quad.$ (51)

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

$\displaystyle S_{n,j,k}$ $\displaystyle=$
$\displaystyle\phi_{n,j}\,\phi_{n,k}\,\left[S_{n-1,j,k}+D_{n-1,n-1}\,\tilde{W}_{n-1,j}\,\tilde{W}_{n-1,k}\right]$
$\displaystyle D_{n,n}$ $\displaystyle=$ $\displaystyle
A_{n,n}-\sum_{j=1}^{2\,J}\sum_{k=1}^{2\,J}\tilde{U}_{n,j}\,S_{n,j,k}\,\tilde{U}_{n,k}$
$\displaystyle\tilde{W}_{n,j}$ $\displaystyle=$
$\displaystyle\frac{1}{D_{n,n}}\left[\tilde{V}_{n,j}-\sum_{k=1}^{2\,J}\tilde{U}_{n,k}\,S_{n,j,k}\right]\quad.$
(52)

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

$\displaystyle f_{n,j}$ $\displaystyle=$
$\displaystyle\phi_{n,j}\,\left[f_{n-1,j}+\tilde{W}_{n-1,j}\,z_{n-1}^{\prime}\right]$
$\displaystyle z_{n}^{\prime}$ $\displaystyle=$ $\displaystyle
y_{n}-\sum_{j=1}^{2\,J}\tilde{U}_{n,j}\,f_{n,j}$ (53)

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

$\displaystyle g_{n,j}$ $\displaystyle=$
$\displaystyle\phi_{n+1,j}\,\left[g_{n+1,j}+\tilde{U}_{n+1,j}\,z_{n+1}\right]$
$\displaystyle z_{n}$ $\displaystyle=$
$\displaystyle\frac{z_{n}^{\prime}}{D_{n,n}}-\sum_{j=1}^{2\,J}\tilde{W}_{n,j}\,g_{n,j}$
(54)

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$.

Figure 2: _(left)_ The maximum numerical error introduced by the algorithm
described in the text as a function of the number of _celerite_ terms $J$ and
data points $N$. Each line corresponds to a different $J$ and the error
increases with $J$ as shown in the legend. Each point in this figure is the
average result across 10 systems with randomly sampled parameters. _(right)_
The distribution of fractional error on the log determinant of $K$ computed
using the algorithm described in the text compared to the algorithm
implemented in NumPy (Van Der Walt et al., 2011) for the same systems shown in
the left-hand figure. This histogram only includes comparisons for $N\leq
2048$ because standard methods become computationally expensive for larger
systems.  Figure 3: A benchmark showing the computational scaling of celerite
using the Cholesky factorization routine for rank-$2\,J$ semiseparable
matrices as described in the text. _(left)_ The cost of computing Equation (5)
with a covariance matrix given by Equation (11) as a function of the number of
data points $N$. The different lines show the cost for different numbers of
terms $J$ increasing from bottom to top. To guide the eye, the straight black
line without points shows linear scaling in $N$. For comparison, the
computational cost for the same model using a general Cholesky factorization
routine implemented in the Intel MKL linear algebra package is shown as the
dashed black line. _(right)_ The same information plotted as a function of $J$
for different values of $N$. Each line shows the scaling for a specific value
of $N$ increasing from bottom to top. The black line shows quadratic scaling
in $J$.

###  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

$\displaystyle\left.\frac{\,\mathrm{d}k_{\mathrm{SHO}}(\tau)}{\,\mathrm{d}\tau}\right|_{\tau=0}=0\quad,$
(55)

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

$\displaystyle
p({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}}\,|\,{\boldsymbol{y}},\,X)\propto
p({\boldsymbol{y}}\,|\,{X,\,{\boldsymbol{\theta}}},\,{\boldsymbol{\alpha}})\,p({\boldsymbol{\theta}},\,{\boldsymbol{\alpha}})$
(56)

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.

Table 1: The parameters and priors for Example 1.  parameter | prior
---|---
$\ln(S_{0})$ | $\mathcal{U}(-10,\,10)$
$\ln(Q)$ | $\mathcal{U}(-10,\,10)$
$\ln(\omega_{0})$ | $\mathcal{U}(-10,\,10)$
Figure 4: (top left) A simulated dataset (black error bars), and MAP model
(blue contours). (bottom left) The residuals between the mean predictive model
and the data shown in the top left figure. (right) The inferred PSD – the blue
contours encompass 68% of the posterior mass – compared to the true PSD
(dashed black line).

###  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

$\displaystyle
k_{\mathrm{true}}(\tau)=A\,\exp\left(-\frac{\tau^{2}}{2\,\lambda^{2}}\right)\,\cos\left(\frac{2\,\pi\,\tau}{P_{\mathrm{true}}}\right)$
(57)

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

$\displaystyle
S_{\mathrm{true}}(\omega)=\frac{\lambda\,A}{2}\,\left[\exp\left(-\frac{\lambda^{2}}{2}\,\left(\omega-\frac{2\,\pi}{P_{\mathrm{true}}}\right)^{2}\right)+\exp\left(-\frac{\lambda^{2}}{2}\,\left(\omega+\frac{2\,\pi}{P_{\mathrm{true}}}\right)^{2}\right)\right]$
(58)

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

$\displaystyle
k(\tau)=k_{\mathrm{SHO}}(\tau;\,S_{0},\,Q,\,\omega_{1})\,k_{\mathrm{SHO}}(\tau;\,S_{0}=1,\,Q=1/\sqrt{2},\,\omega_{2})$
(59)

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)

$\displaystyle\mathrm{BIC}=-2\,\mathcal{L}^{*}+K\,\log N$ (60)

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.

Table 2: The parameters and priors for Example 2.  parameter | prior
---|---
$\ln(S_{0})$ | $\mathcal{U}(-15,\,5)$
$\ln(Q)$ | $\mathcal{U}(-10,\,10)$
$\ln(\omega_{1})$ | $\mathcal{U}(-10,\,10)$
$\ln(\omega_{2})$ | $\mathcal{U}(-5,\,5)$
Figure 5: (top left) A simulated dataset (black error bars), and MAP model
(blue contours). (bottom left) The residuals between the mean predictive model
and the data shown in the top left figure. (right) The inferred period of the
process. The true period is indicated by the vertical orange line, the
posterior inference using the correct model is shown as the blue dashed
histogram, and the inference made using the “wrong” effective model is shown
as the black histogram.  Figure 6:  The Bayesian information criterion (BIC)
evaluated for the data from Figure 5 and different kernel choices. The
$x$-axis shows the number of terms included in each model. The blue circles
show models where the terms have been summed and the orange squares indicate
that the model is a product of terms. In the product models, the parameters
$S_{0}$ and $Q$ are held fixed at $1$ and $1/\sqrt{2}$ respectively for all
but the first term. The dashed black line shows the minium value of the BIC
and this corresponds to the model chosen in the text.

###  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

$\displaystyle
k(\tau)=A\,\exp\left(-\frac{\tau^{2}}{2\,\ell^{2}}-\Gamma\,\sin^{2}\left(\frac{\pi\,\tau}{P_{\mathrm{rot}}}\right)\right)$
(61)

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

$\displaystyle
k(\tau)=\frac{B}{2+C}\,e^{-\tau/L}\,\left[\cos\left(\frac{2\,\pi\,\tau}{P_{\mathrm{rot}}}\right)+(1+C)\right]$
(62)

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.

Table 3: The parameters and priors for Example 3.  parameter | prior
---|---
$\ln(B/\mathrm{ppt}^{2})$ | $\mathcal{U}(-10.0,\,0.0)$
$\ln(L/\mathrm{day})$ | $\mathcal{U}(1.5,\,5.0)$
$\ln(P/\mathrm{day})$ | $\mathcal{U}(-3.0,\,5.0)$
$\ln(C)$ | $\mathcal{U}(-5.0,\,5.0)$
Figure 7: Inferred constraints on a quasiperiodic GP model using the
covariance function in Equation (62) and two quarters of Kepler data. (top
left) The Kepler data (black points) and the MAP model prediction (blue curve)
for a 60 day subset of the data used. The solid blue line shows the predictive
mean and the blue contours show the predictive standard deviation. (bottom
left) The residuals between the mean predictive model and the data shown in
the top left figure. (right) The posterior constraint on the rotation period
of KIC 1430163 using the dataset and model from Figure 7. The period is the
parameter $P_{\mathrm{rot}}$ in Equation (62) and this figure shows the
posterior distribution marginalized over all other nuisance parameters in
Equation (62). This is consistent with the published rotation period made
using the full Kepler baseline shown as a vertical gray line (Mathur et al.,
2014).

###  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

$\displaystyle\omega_{0,\,j}=2\,\pi\,(\nu_{\mathrm{max}}+j\,\Delta\nu+\epsilon)$
(63)

and the amplitudes are given by

$\displaystyle
S_{0,\,j}=\frac{A}{Q^{2}}\,\exp\left(-\frac{[j\,\Delta\nu+\epsilon]^{2}}{2\,W^{2}}\right)$
(64)

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)

$\displaystyle\nu_{\mathrm{max}}=171.94\pm
3.62\,\mu\mathrm{Hz}\quad\mathrm{and}\quad\Delta\nu=13.28\pm
0.29\,\mu\mathrm{Hz}\quad.$ (65)

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.

Table 4: The parameters and priors for Example 4.  parameter | prior
---|---
$\ln(S_{g}/\mathrm{ppm}^{2})$ | $\mathcal{U}(-15,\,15)$
$\ln(\omega_{g}/\mathrm{day}^{-1})$ | $\mathcal{U}(-15,\,15)$
$\ln(\nu_{\mathrm{max}}/\mu\mathrm{Hz})$ | $\mathcal{U}(\ln(130),\,\ln(190))$
$\ln(\Delta\nu/\mu\mathrm{Hz})$ | $\mathcal{U}(\ln(12.5),\,\ln(13.5))$
$\epsilon/\mathrm{day}^{-1}$ | $\mathcal{N}(0,\,1)$
$\ln(A/\mathrm{ppm}^{2}\,\mathrm{day})$ | $\mathcal{U}(-15,\,15)$
$\ln(Q)$ | $\mathcal{U}(-15,\,15)$
$\ln(W/\mathrm{day}^{-1})$ | $\mathcal{U}(-3,\,3)$
$\ln(\sigma/\mathrm{ppm})$ | $\mathcal{U}(-15,\,15)$
Figure 8:  (top) The Kepler data (black points) and the MAP model prediction
(blue curve) for a 10 day subset of the month-long dataset that was used for
the fit. The solid blue line shows the predictive mean and the blue contours
show the predictive standard deviation. (bottom) The residuals between the
mean predictive model and the data shown in the top figure.  Figure 9: The
probabilistic constraints on $\nu_{\mathrm{max}}$ and $\Delta\nu$ from the
inference shown in Figure 8 compared to the published value (error bar) based
on several years of Kepler observations (Pinsonneault et al., 2014). The two-
dimensional contours show the 0.5-, 1-, 1.5, and 2-sigma credible regions in
the marginalized planes and the histograms along the diagonal show the
marginalized posterior for each parameter.  Figure 10: A comparison between
the Lomb-Scargle estimator of the PSD and the posterior inference of the PSD
as a mixture of stochastically-driven simple harmonic oscillators. (top) The
periodogram of the Kepler light curve for KIC 11615890 computed on the full
four year baseline of the mission. The orange line shows a smoothed
periodogram and the blue line indicates the level of the measurement
uncertainties. (middle) The same periodogram computed using about a month of
data. (bottom) The power spectrum inferred using the mixture of SHOs model
described in the text and only one month of Kepler data. The black line shows
the median of posterior PSD and the gray contours show the 68% credible
region.

###  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.

Table 5: The parameters, priors, and (if known) the true values used for the simulation in Example 5.  | parameter | prior | true value
---|---|---|---
kernel: ${\boldsymbol{\alpha}}\quad$ | $\ln(B/\mathrm{ppt}^{2})$ | $\mathcal{U}(-10.0,\,0.0)$ | —
$\ln(L/\mathrm{day})$ | $\mathcal{U}(1.5,\,5.0)$ | —
$\ln(P_{\mathrm{rot}}/\mathrm{day})$ | $\mathcal{U}(-3.0,\,5.0)$ | —
$\ln(C)$ | $\mathcal{U}(-5.0,\,5.0)$ | —
$\ln(\sigma/\mathrm{ppt})$ | $\mathcal{U}(-5.0,\,0.0)$ | —
mean: ${\boldsymbol{\theta}}\quad$ | $f_{0}/\mathrm{ppt}$ | $\mathcal{U}(-0.5,\,0.5)$ | 0
$\ln(P_{\mathrm{orb}}/\mathrm{day})$ | $\mathcal{U}(\ln(7.9),\,\ln(8.1))$ | $\ln(8)$
$\ln(R_{p}/R_{\star})$ | $\mathcal{U}(\ln(0.005),\,\ln(0.1))$ | $\ln(0.015)$
$\ln(T/\mathrm{day})$ | $\mathcal{U}(\ln(0.4),\,\ln(0.6))$ | $\ln(0.5)$
$t_{0}/\mathrm{day}$ | $\mathcal{U}(-0.1,\,0.1)$ | $0$
$b$ | $\mathcal{U}(0,\,1.0)$ | $0.5$
$q_{1}$ | $\mathcal{U}(0,\,1)$ | $0.5$
$q_{2}$ | $\mathcal{U}(0,\,1)$ | $0.5$
Figure 11: _(top)_ A month-long segment of Kepler light curve for KIC 1430163
with a synthetic transit model injected (black points) and the MAP model for
the stellar variability (blue line). _(bottom)_ The MAP “de-trending” of the
data in the top panel. In this panel, the MAP model for the stellar
variability has been subtracted to leave only the transits. The de-trended
light curve is shown by black error bars and the MAP transit model is shown as
a blue line.  Figure 12: The marginalized posterior constraints on the
physical parameters of the planet transit in the light curve shown in the top
panel of Figure 11. The two-dimensional contours show the 0.5-, 1-, 1.5, and
2-sigma credible regions in the marginalized planes and the histograms along
the diagonal show the marginalized posterior for each parameter. The true
values used in the simulation are indicated by blue lines. For each parameter,
the inference is consistent with the true value.

###  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.

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

##  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

$\displaystyle
k_{\mathrm{CARMA}}(\tau)=\sum_{j=1}^{J}A_{j}\,\exp\,(r_{j}\,\tau)$ (66)

where

$\displaystyle
A_{j}=\sigma^{2}\,\frac{\left[\sum_{k=0}^{K}\beta_{k}\,{(r_{j})}^{k}\right]\,\left[\sum_{k=0}^{K}\beta_{k}\,{(-r_{j})}^{k}\right]}{-2\,\mathrm{Re}(r_{j})\,\prod_{k=1,\,k\neq
j}^{J}(r_{k}-r_{j})({r_{k}}^{*}+r_{j})}$ (67)

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).

Figure 13: The same as Figure 3, but using an optimized C++ implementation of
the Kalman filter method developed by Kelly et al. (2014). Comparing this
figure to Figure 3, we see that celerite is about an order of magnitude faster
in all cases. _(left)_ The cost of computing Equation (5) with a covariance
matrix given by Equation (11) as a function of the number of data points $N$.
The different lines show the cost for different numbers of terms $J$
increasing from bottom to top. To guide the eye, the straight black line
without points shows linear scaling in $N$. _(right)_ The same information
plotted as a function of $J$ for different values of $N$. Each line shows the
scaling for a specific value of $N$ increasing from bottom to top. The black
line shows quadratic scaling in $J$.

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

$\displaystyle\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}}>0\quad.$
(A1)

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

$\displaystyle a_{j}\,c_{j}>-b_{j}\,d_{j}$ (A2) $\displaystyle
a_{j}\,c_{j}>b_{j}\,d_{j}\quad.$ (A3)

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

$\displaystyle|b_{j}\,d_{j}|<a_{j}\,c_{j}\quad.$ (A4)

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}$:

$S(\omega)=\sum_{j=1}^{J}\frac{q_{j}\,z+r_{j}}{z^{2}+s_{j}\,z+t_{j}}=0$ (A5)

where

$\displaystyle q_{j}$ $\displaystyle=$ $\displaystyle
a_{j}\,c_{j}-b_{j}\,d_{j}$ (A6) $\displaystyle r_{j}$ $\displaystyle=$
$\displaystyle(d_{j}^{2}+c_{j}^{2})(b_{j}\,d_{j}+a_{j}\,c_{j})$ (A7)
$\displaystyle s_{j}$ $\displaystyle=$ $\displaystyle 2(c_{j}^{2}-d_{j}^{2})$
(A8) $\displaystyle t_{j}$ $\displaystyle=$
$\displaystyle(c_{j}^{2}+d_{j}^{2})^{2}.$ (A9)

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

$Q_{0}(z)=\sum_{j=1}^{J}(q_{j}z+r_{j})\prod_{k\neq
j}\left(z^{2}+s_{k}z+t_{k}\right)=0\quad,$ (A10)

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

$\displaystyle Q_{k}(z)=-{\rm rem}(Q_{k-2},\,Q_{k-1})$ (A11)

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

$\displaystyle{\boldsymbol{y}}=K\,{\boldsymbol{z}}$ (B1)

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

$\displaystyle f^{+}_{n,j}$ $\displaystyle=$
$\displaystyle\phi_{n+1,j}\,\left[f^{+}_{n+1,j}+\tilde{U}_{n,j}\,z_{n+1}\right]$
(B2) $\displaystyle f^{-}_{n,j}$ $\displaystyle=$
$\displaystyle\phi_{n,j}\,\left[f^{-}_{n-1,j}+\tilde{V}_{n-1,j}\,z_{n-1}\right]$
(B3) $\displaystyle y_{n}$ $\displaystyle=$ $\displaystyle
A_{n,n}\,z_{n}+\sum_{j=1}^{2\,J}\left[\tilde{V}_{n,j}f^{+}_{n,j}+\tilde{U}_{n-1,j}f^{-}_{n,j}\right]\quad,$
(B4)

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

${\boldsymbol{y}}={\boldsymbol{\mu}}_{\boldsymbol{\theta}}({\boldsymbol{t}})+K_{\boldsymbol{\alpha}}^{1/2}\,{\boldsymbol{q}}$
(B5)

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

$\displaystyle
K_{\boldsymbol{\alpha}}=K_{\boldsymbol{\alpha}}^{1/2}\,{\left(K_{\boldsymbol{\alpha}}^{1/2}\right)}^{T}\quad.$
(B6)

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

$\displaystyle f_{n,j}$ $\displaystyle=$
$\displaystyle\phi_{n,j}\,\left[f_{n-1,j}+\tilde{W}_{n-1,j}\,\sqrt{D_{n-1,n-1}}\,q_{n-1}\right]$
(B7) $\displaystyle y_{n}$ $\displaystyle=$
$\displaystyle\sqrt{D_{n,n}}\,q_{n}+\sum_{j=1}^{2\,J}\tilde{U}_{n,j}\,f_{n,j}$
(B8)

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

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

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

$\displaystyle{\boldsymbol{\mu}}^{*}$ $\displaystyle=$
$\displaystyle{\boldsymbol{\mu}}_{\boldsymbol{\theta}}({\boldsymbol{t}}^{*})+K({\boldsymbol{t}}^{*},\,{\boldsymbol{t}})\,{K({\boldsymbol{t}},\,{\boldsymbol{t}})}^{-1}\,\left[{\boldsymbol{y}}-{\boldsymbol{\mu}}_{\boldsymbol{\theta}}({\boldsymbol{t}})\right]$
(B11)

and covariance

$\displaystyle K^{*}$ $\displaystyle=$ $\displaystyle
K({\boldsymbol{t}}^{*},\,{\boldsymbol{t}}^{*})-K({\boldsymbol{t}}^{*},\,{\boldsymbol{t}})\,{K({\boldsymbol{t}},\,{\boldsymbol{t}})}^{-1}\,K({\boldsymbol{t}},\,{\boldsymbol{t}}^{*})$
(B12)

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)

$\displaystyle\mu^{*}_{m}$ $\displaystyle=$ $\displaystyle
K(t^{*}_{m},{\boldsymbol{t}}){\boldsymbol{z}}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{j=1}^{J}e^{-c_{j}|t^{*}_{m}-t_{n}|}\,\left[a_{j}\,\cos{(d_{j}\,|t^{*}_{m}-t_{n}|)}+b_{j}\,\sin{(d_{j}\,|t^{*}_{m}-t_{n}|)}\right]\,z_{n}\quad.$

where

$\displaystyle{\boldsymbol{z}}={K({\boldsymbol{t}},\,{\boldsymbol{t}})}^{-1}\,\left[{\boldsymbol{y}}-{\boldsymbol{\mu}}_{\boldsymbol{\theta}}({\boldsymbol{t}})\right]\quad.$
(B14)

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

$\displaystyle\mu^{*}_{m}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{J}\sum_{n=1}^{n_{0}}e^{-c_{j}\,(t^{*}_{m}-t_{n})}\,\left[a_{j}\,\cos{(d_{j}\,(t^{*}_{m}-t_{n}))}+b_{j}\,\sin{(d_{j}\,(t^{*}_{m}-t_{n}))}\right]\,z_{n}$
$\displaystyle+$
$\displaystyle\sum_{j=1}^{J}\sum_{n=n_{0}+1}^{N}e^{-c_{j}\,(t_{n}-t^{*}_{m})}\,\left[a_{j}\,\cos{(d_{j}\,(t_{n}-t^{*}_{m}))}+b_{j}\,\sin{(d_{j}\,(t_{n}-t^{*}_{m}))}\right]\,z_{n}\quad.$

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

$\displaystyle Q^{-}_{n,k}$ $\displaystyle=$
$\displaystyle\left[Q^{-}_{n-1,k}+z_{n}\,\tilde{V}_{n,k}\right]\,e^{-c_{k//2}\,(t_{n+1}-t_{n})},$
(B16) $\displaystyle X^{-}_{m,n,k}$ $\displaystyle=$ $\displaystyle
e^{-c_{k//2}\,(t_{m}^{*}-t_{n+1})}\,\tilde{U}_{m,k}^{*},$ (B17) $\displaystyle
Q^{+}_{n,k}$ $\displaystyle=$
$\displaystyle\left[Q^{+}_{n+1,k}+z_{n}\,\tilde{U}_{n,k}\right]\,e^{-c_{k//2}(t_{n}-t_{n-1})},$
(B18) $\displaystyle X^{+}_{m,n,k}$ $\displaystyle=$ $\displaystyle
e^{-c_{k//2}\,(t_{n-1}-t_{m}^{*})}\,\tilde{V}_{m,k}^{*}\quad,$ (B19)

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

$\mu^{*}_{m}=\sum_{k=1}^{2\,J}\left[Q^{-}_{n_{0},k}X^{-}_{m,n_{0},k}+Q^{+}_{n_{0}+1,k}X^{+}_{m,n_{0}+1,k}\right]\quad.$
(B20)

## References

  * Agol et al. (2005) Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, Monthly Notices of the Royal Astronomical Society, 359, 567
  * Aigrain et al. (2016) Aigrain, S., Parviainen, H., & Pope, B. J. S. 2016, Monthly Notices of the Royal Astronomical Society, 706
  * Aigrain et al. (2015) Aigrain, S., Llama, J., Ceillier, T., et al. 2015, Monthly Notices of the Royal Astronomical Society, 450, 3211
  * Ambikasaran (2015) Ambikasaran, S. 2015, Numer. Linear Algebra Appl., 22, 1102
  * Ambikasaran et al. (2016) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2016, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bischof, C., et al. 1999, LAPACK Users’ Guide, 3rd edn. (Philadelphia, PA: Society for Industrial and Applied Mathematics)
  * Anderson et al. (1990) Anderson, E. R., Duvall, Jr., T. L., & Jefferies, S. M. 1990, ApJ, 364, 699
  * Angus et al. (2017) Angus, R., Morton, T., Aigrain, S., Foreman-Mackey, D., & Rajpaul, V. 2017, ArXiv e-prints, arXiv:1706.05459
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Barclay et al. (2015) Barclay, T., Endl, M., Huber, D., et al. 2015, ApJ, 800, 46
  * Bezanzon et al. (2012) Bezanzon, J., Karpinski, S., Shah, V., & Edelman, A. 2012, in Lang.NEXT
  * Bond et al. (1999) Bond, J. R., Crittenden, R. G., Jaffe, A. H., & Knox, L. 1999, Comput. Sci. Eng., Vol. 1, No. 2, p. 21 - 35, 1, 21
  * Bond & Efstathiou (1987) Bond, J. R., & Efstathiou, G. 1987, Monthly Notices of the Royal Astronomical Society, 226, 655
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
  * Brewer & Stello (2009) Brewer, B. J., & Stello, D. 2009, Monthly Notices of the Royal Astronomical Society, 395, 2226
  * Brown (2001) Brown, T. M. 2001, The Astrophysical Journal, 553, 1006
  * Byrd et al. (1995) Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM Journal on Scientific Computing, 16, 1190
  * Campante et al. (2016) Campante, T. L., Schofield, M., Kuszlewicz, J. S., et al. 2016, The Astrophysical Journal, 830, 138
  * Carter & Winn (2009) Carter, J. A., & Winn, J. N. 2009, The Astrophysical Journal, 704, 51
  * Chaplin et al. (2011) Chaplin, W. J., Kjeldsen, H., Christensen-Dalsgaard, J., et al. 2011, Science, 332, 213
  * Chaplin et al. (2013) Chaplin, W. J., Basu, S., Huber, D., et al. 2013, The Astrophysical Journal Supplement Series, 210, 1
  * Claret & Bloemen (2011) Claret, A., & Bloemen, S. 2011, Astronomy & Astrophysics, 529, A75
  * Corsaro & Ridder (2014) Corsaro, E., & Ridder, J. D. 2014, Astronomy & Astrophysics, 571, A71
  * Czekala et al. (2017) Czekala, I., Mandel, K. S., Andrews, S. M., et al. 2017, ApJ, 840, 49
  * Dawson et al. (2014) Dawson, R. I., Johnson, J. A., Fabrycky, D. C., et al. 2014, ApJ, 791, 89
  * Dillon et al. (2013) Dillon, J. S., Liu, A., & Tegmark, M. 2013, Phys. Rev. D, 87, 043005
  * Dörrie (1965) Dörrie, H. 1965, 100 Great Problems of Elementary Mathematics: Their History and Solution, Dover Books on Mathematics Series, §24, (Dover Publications), 112–116
  * Dumusque et al. (2014) Dumusque, X., Boisse, I., & Santos, N. C. 2014, The Astrophysical Journal, 796, 132
  * Evans et al. (2015) Evans, T. M., Aigrain, S., Gibson, N., et al. 2015, Monthly Notices of the Royal Astronomical Society, 451, 680
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 24, doi:10.21105/joss.00024
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, dfm/celerite: celerite v0.2.1, , , doi:10.5281/zenodo.806847
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hoyer, S., Bernhard, J., & Angus, R. 2014, george: George (v0.2.0), , , doi:10.5281/zenodo.11989
  * Foreman-Mackey & Morton (2016) Foreman-Mackey, D., & Morton, T. 2016, dfm/transit: v0.3.0, , , doi:10.5281/zenodo.159478
  * Foreman-Mackey et al. (2016) Foreman-Mackey, D., Morton, T. D., Hogg, D. W., Agol, E., & Schölkopf, B. 2016, AJ, 152, 206
  * Gibson et al. (2012) Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012, MNRAS, 419, 2683
  * Gilliland et al. (2010) Gilliland, R. L., Brown, T. M., Christensen-Dalsgaard, J., et al. 2010, Publications of the Astronomical Society of the Pacific, 122, 131
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Gould et al. (2015) Gould, A., Huber, D., Penny, M., & Stello, D. 2015, Journal of The Korean Astronomical Society, 48, 93
  * Grunblatt et al. (2016) Grunblatt, S. K., Huber, D., Gaidos, E. J., et al. 2016, AJ, 152, 185
  * Guennebaud et al. (2010) Guennebaud, G., Jacob, B., et al. 2010, Eigen v3, http://eigen.tuxfamily.org, ,
  * Harvey (1985) Harvey, J. 1985, in ESA Special Publication, Vol. 235, Future Missions in Solar, Heliospheric & Space Plasma Physics, ed. E. Rolfe & B. Battrick
  * Haywood et al. (2014) Haywood, R. D., Cameron, A. C., Queloz, D., et al. 2014, Monthly Notices of the Royal Astronomical Society, 443, 2517
  * Holman (2005) Holman, M. J. 2005, Science, 307, 1288
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398
  * Huber et al. (2009) Huber, D., Stello, D., Bedding, T. R., et al. 2009, Communications in Asteroseismology, 160, 74
  * Huber et al. (2011) Huber, D., Bedding, T. R., Stello, D., et al. 2011, ApJ, 743, 143
  * Hunter et al. (2007) Hunter, J. D., et al. 2007, Computing in science and engineering, 9, 90
  * Ivezić et al. (2008) Ivezić, v., Tyson, J. A., Acosta, E., et al. 2008, arXiv:0805.2366v4
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python, ,
  * Kallinger et al. (2014) Kallinger, T., De Ridder, J., Hekker, S., et al. 2014, A&A, 570, A41
  * Kelly et al. (2014) Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., & Uttley, P. 2014, ApJ, 788, 33
  * Kelly et al. (2011) Kelly, B. C., Sobolewska, M., & Siemiginowska, A. 2011, ApJ, 730, 52
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152
  * Knutson et al. (2007) Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183
  * Littlefair et al. (2017) Littlefair, S. P., Burningham, B., & Helling, C. 2017, MNRAS, 466, 4250
  * Loeb & Gaudi (2003) Loeb, A., & Gaudi, B. S. 2003, The Astrophysical Journal, 588, L117
  * Luger et al. (2016) Luger, R., Agol, E., Kruse, E., et al. 2016, AJ, 152, 100
  * Luger et al. (2017) Luger, R., Kruse, E., Foreman-Mackey, D., Agol, E., & Saunders, N. 2017, ArXiv e-prints, arXiv:1702.05488
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, ApJ, 721, 1014
  * Mandel & Agol (2002) Mandel, K., & Agol, E. 2002, The Astrophysical Journal, 580, L171
  * Mathur et al. (2014) Mathur, S., García, R. A., Ballot, J., et al. 2014, A&A, 562, A124
  * McAllister et al. (2016) McAllister, M. J., Littlefair, S. P., Dhillon, V. S., et al. 2016, Monthly Notices of the Royal Astronomical Society, 464, 1353
  * McQuillan et al. (2013) McQuillan, A., Aigrain, S., & Mazeh, T. 2013, MNRAS, 432, 1203
  * McQuillan et al. (2014) McQuillan, A., Mazeh, T., & Aigrain, S. 2014, ApJS, 211, 24
  * Michel et al. (2009) Michel, E., Samadi, R., Baudin, F., et al. 2009, A&A, 495, 979
  * Nocedal & Wright (2006) Nocedal, J., & Wright, S. J. 2006, Numerical Optimization (Springer)
  * Pancoast et al. (2014) Pancoast, A., Brewer, B. J., Treu, T., et al. 2014, Monthly Notices of the Royal Astronomical Society, 445, 3073
  * Pinsonneault et al. (2014) Pinsonneault, M. H., Elsworth, Y., Epstein, C., et al. 2014, ApJS, 215, 19
  * Pope et al. (2016) Pope, B. J. S., Parviainen, H., & Aigrain, S. 2016, MNRAS, 461, 3399
  * Press & Rybicki (1998) Press, W. H., & Rybicki, G. B. 1998, ApJ, 507, 108
  * Rajpaul et al. (2015) Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S., & Roberts, S. 2015, Monthly Notices of the Royal Astronomical Society, 452, 2269
  * Rasmussen & Williams (2006) Rasmussen, C. E., & Williams, K. I. 2006, Gaussian Processes for Machine Learning (MIT Press)
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. M. Oschmann, M. Clampin, G. G. Fazio, & H. A. MacEwen (SPIE-Intl Soc Optical Eng)
  * Rybicki & Press (1992) Rybicki, G. B., & Press, W. H. 1992, ApJ, 398, 169
  * Rybicki & Press (1995) —. 1995, Physical Review Letters, 74, 1060
  * Schwarz et al. (1978) Schwarz, G., et al. 1978, The annals of statistics, 6, 461
  * Smith et al. (2012) Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000
  * Sokal (1989) Sokal, A. D. 1989, Monte Carlo methods in statistical mechanics: foundations and new algorithms, Troisieme cycle de la physique en Suisse Romande
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints, arXiv:1503.03757
  * Stello et al. (2016) Stello, D., Cantiello, M., Fuller, J., et al. 2016, Nature, 529, 364
  * Stello et al. (2009) Stello, D., Chaplin, W. J., Basu, S., Elsworth, Y., & Bedding, T. R. 2009, MNRAS, 400, L80
  * Stello et al. (2013) Stello, D., Huber, D., Bedding, T. R., et al. 2013, The Astrophysical Journal, 765, L41
  * Stumpe et al. (2012) Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985
  * Uhlenbeck & Ornstein (1930) Uhlenbeck, G. E., & Ornstein, L. S. 1930, Physical review, 36, 823
  * Uttley et al. (2005) Uttley, P., McHardy, I. M., & Vaughan, S. 2005, Monthly Notices of the Royal Astronomical Society, 359, 345
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Vandebril et al. (2007) Vandebril, R., Van Barel, M., & Mastronardi, N. 2007, Matrix computations and semiseparable matrices: linear systems, Vol. 1 (JHU Press)
  * VanderPlas (2017) VanderPlas, J. T. 2017, ArXiv e-prints, arXiv:1703.09824
  * Verner et al. (2011) Verner, G. A., Elsworth, Y., Chaplin, W. J., et al. 2011, MNRAS, 415, 3539
  * Wandelt & Hansen (2003) Wandelt, B. D., & Hansen, F. K. 2003, Phys. Rev. D, 67, 023001
  * Wang et al. (2012) Wang, Y., Khardon, R., & Protopapas, P. 2012, ApJ, 756, 67
  * Wilson (2014) Wilson, A. G. 2014, PhD thesis, University of Cambridge
  * Wilson & Adams (2013) Wilson, A. G., & Adams, R. P. 2013, in ICML (3), 1067–1075
  * Wilson et al. (2015) Wilson, A. G., Dann, C., & Nickisch, H. 2015, arXiv preprint arXiv:1511.01870, http://arxiv.org/abs/1511.01870
  * Wilson & Nicisch (2015) Wilson, A. G., & Nicisch, H. 2015, Proceedings of the 32nd International Conference on Machine Learning, 1775
  * Zhu et al. (1997) Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Transactions on Mathematical Software (TOMS), 23, 550
  * Zinn et al. (2016) Zinn, J. C., Kochanek, C. S., Kozłowski, S., et al. 2016, ArXiv e-prints, arXiv:1612.04834
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80
  * Zucker et al. (2007) Zucker, S., Mazeh, T., & Alexander, T. 2007, The Astrophysical Journal, 670, 1326