PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Daniel Foreman-Mackey, Eric Agol, Sivaram Ambikasaran, Ruth Angus

arXiv:1703.09710v2astro-ph.IMastro-ph.EPastro-ph.SRphysics.data-anstat.AP

TL;DR

Large astronomical time-series datasets make general GP inference difficult because likelihood computation typically scales cubically with the number of observations. The paper introduces celerite, an exact linear-scaling method for one-dimensional GP models using structured exponential-mixture covariances, and demonstrates applications across astronomical time series. The method supports scalable inference while retaining a physical oscillator interpretation, but it is restricted to one-dimensional problems.

  • Problem

    General GP likelihood computation scales as O(N 3), limiting applications to large astronomical datasets.

  • Method

    Celerite exactly exploits semiseparable covariance structure for one-dimensional GP models, achieving O(N) scaling while allowing uneven sampling and heteroscedastic uncertainties.

  • Results

    The method is demonstrated on stellar rotation, asteroseismic oscillation, and exoplanet transit analyses, including examples that would otherwise be intractable without substantial computational investment.

  • Takeaways & Limitations

    Celerite can expand GP use to existing and forthcoming large astronomical time-domain surveys, and its kernel can describe a wide range of astrophysical variability.

  • Takeaways & Limitations

    Celerite is restricted to one-dimensional problems, so it cannot speed up many multidimensional astrophysical GP analyses.

Abstract

from arXiv · show

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.

1. INTRODUCTION

Gaussian Processes are widely used in astronomical data analysis, but cubic computational scaling limits their use on large time-series datasets. The paper introduces celerite, an exact linear-scaling method for one-dimensional data, and demonstrates applications across several astronomical problems.

  • Motivation: Gaussian Processes model stochastic variability and correlated noise across diverse astronomical applications.Examples include stellar, active galactic nucleus, and X-ray binary light curves, as well as cosmic microwave background and calibration problems.
  • Motivation: O(N 3) likelihood computation becomes prohibitive for surveys containing many targets and observations.Large time-domain surveys may contain approximately 10^4−10^9 targets with approximately 10^3−10^5 observations each.
  • Motivation: Existing and forthcoming surveys already produce datasets beyond the tractable range of naïve GP modeling.Kepler measured more than 60,000 observations for about 190,000 stars, with K2, TESS, LSST, WFIRST, and PLATO producing similar or larger volumes.
  • Contribution: Celerite directly and exactly computes a class of one-dimensional GP models in O(N) operations.It requires a covariance represented as a mixture of exponentials, while allowing uneven sampling and heteroscedastic uncertainties.
  • Applications: The method exploits covariance structure and is demonstrated on stellar rotation, asteroseismic analysis, and exoplanet transit fitting.The paper also provides well-tested, documented open-source implementations in C++, Python, and Julia.

2. GAUSSIAN PROCESSES

Gaussian Processes define probabilistic models through mean and covariance functions, but general likelihood evaluation requires cubic time and quadratic storage. Celerite targets this bottleneck with a specialized one-dimensional, stationary covariance model.

  • GP formulation: A GP combines a parameterized mean function with a parameterized covariance or kernel function.The covariance matrix contains kernel evaluations between pairs of input coordinates, and the likelihood supports parameter estimation and uncertainty quantification.
  • Computational challenge: O(N 3) computation and O(N 2) storage make general GP likelihood evaluations impractical for large datasets.Iterative optimization and MCMC become impractical when each likelihood evaluation requires expensive matrix inversion, determinant calculation, and storage.
  • Computational challenge: 8,192 measurements require over 8 seconds for a direct likelihood calculation, compared with less than a tenth of a second for 1,024 measurements.The paper notes that repeated evaluations make this scaling especially problematic for parameter estimation and inference.
  • Celerite scope: Celerite improves cubic scaling for one-dimensional datasets using a specific covariance model.Here, one-dimensional means scalar input coordinates, such as time, and the method can also apply to other scalar-coordinate datasets.
  • Celerite scope: The celerite covariance is stationary, depending on input pairs only through τ_nm = |t_n−t_m|.The method’s applicability is therefore defined by scalar coordinates and a covariance that depends only on coordinate separation.

3. THE CELERITE MODEL

The celerite model represents covariance functions with exponential components whose matrix structure enables fast GP calculations. This representation supports exact scalable inference and admits oscillator-based interpretations and algebraic closure under products.

  • Model structure: For a simple exponential covariance, the inverse matrix is tridiagonal and can be computed with a small number of operations per data point.The model includes measurement uncertainties and generalizes to arbitrary mixtures of exponentials.
  • Computational scaling: Mixtures of exponential components allow likelihood evaluation in O(N J2) operations, where J is the number of mixture components.The covariance inverse is dense in this formulation, but the likelihood remains scalable.
  • Model flexibility: Complex parameters extend the covariance family while preserving O(N J2) likelihood evaluation.The resulting covariance is expressed through complex-valued amplitudes and decay or frequency parameters.
  • Physical interpretation: Rewriting exponential terms with sine and cosine functions reveals a mixture of quasiperiodic oscillators.A celerite term has a power spectral density obtained from the Fourier transform of its covariance function.
  • Kernel composition: Products of celerite terms can be rewritten as sums with updated parameters, so kernels formed from sums or products remain amenable to scalable inference.This closure property extends the method beyond individual covariance components.

4. CELERITE AS A MODEL OF STELLAR VARIATIONS

Celerite has a physical interpretation as a mixture of stochastically driven, damped harmonic oscillators, providing models for diverse stellar variability. Quality factor controls the transition from non-oscillatory behavior to sharp oscillations, while some kernel limits require approximation.

  • Physical motivation: Celerite addresses concerns about physically motivated kernel choices in astrophysical GP modeling.The paper connects its kernel to physical systems relevant to stellar time variability.
  • Physical motivation: Stellar oscillations can be excited by noisy processes, strengthened near characteristic timescales, and damped by dissipation.For a stochastically driven damped SHO, the driving force is modeled as white noise and the resulting process has a specific power spectrum.
  • SHO representation: The SHO power spectrum matches the celerite form through parameters tied to oscillator frequency, quality factor, and power.For quality factors below or equal to one-half, the oscillator is represented by a pair of celerite terms with zero characteristic frequencies.
  • Limitations: The Matérn-3/2 limit cannot be evaluated directly because its celerite parameter b_j is infinite at Q = 1/2.An approximate implementation uses a small f, whose required value depends on the dataset and precision requirements.
  • Quality-factor limits: Q ≤ 1/2 produces no oscillatory behavior, whereas large Q yields a PSD whose peak approaches a Lorentzian shape.Figure 1 compares power spectra, autocorrelation functions, and simulated realizations across quality factors.
  • Astrophysical applications: Low Q ≈1 models granulation noise, while high Q ≫1 models asteroseismic oscillations; sums over oscillators broaden the stellar-variability model.The model uses different oscillator values of Q, S0, and ω0 to account for the power spectrum of stellar variations.

5. SEMISEPARABLE MATRICES & CELERITE

The paper derives a direct solver for celerite covariance matrices by exploiting semiseparable structure, enabling scalable factorization, solves, and log-determinant calculations. The method is numerically stable and empirically follows the predicted computational scaling.

  • Semiseparable matrices: Semiseparable representations reduce storage to (2R + 1)N values and support matrix-vector and matrix-matrix products in O(N) operations.The representation stores U, V, and the diagonal of A instead of the full covariance matrix.
  • Scalable computation: The derived factorization costs O(NR^2), while inverse application and triangular solves cost O(NR), and log-determinants cost O(N).For celerite models, the factorization therefore scales as O(NJ^2) and solves as O(NJ).
  • Performance: The method is about 20 times faster than an earlier general semiseparable solver when the covariance matrix is positive definite and applicable.The comparison concerns the factorization, solving, and determinant algorithm for positive definite rank-R semiseparable matrices.
  • Celerite representation: Celerite covariance matrices can be represented as rank R = 2J semiseparable matrices, with real terms reducing the rank contribution.For J_real real and J_complex complex terms, the rank is R = J_real + 2J_complex = 2J − J_real.
  • Numerical stability: Reparameterizing exponentially growing and decaying factors avoids numerical overflow and underflow in realistic datasets.The pre-conditioned variables yield a numerically stable Cholesky algorithm while preserving the required celerite structure.
  • Celerite solver: The celerite implementation uses O(NJ^2) operations and a memory footprint of (6J + 1)N + J(J − 1)/2 floating-point numbers.The implementation updates relevant quantities in place to reduce storage requirements.
  • Model selection: Bayesian information criteria such as AIC and BIC are not generally recommended because they are formally valid only under specific assumptions.The paper uses BIC in one application but cautions against general use.

6. EXAMPLES

The examples apply celerite to simulated and real astronomical time series, covering stellar rotation, asteroseismology, and exoplanet transit inference. Across these cases, the method supports probabilistic modeling with substantial computational savings, while the real-data demonstrations remain qualitative and non-exhaustive.

  • Simulated data: Celerite recovers the true power spectral density when a simulated process is generated and fit with the same model.Posterior sampling produced 1,737 effective samples, and the inferred PSD reproduced the true spectral density.
  • Simulated data: Even when the effective model is formally wrong, inferred periods remain consistent with the correct-model inference and the true simulated period.The result demonstrates computationally efficient effective modeling when the true process lies outside the allowed model space.
  • Stellar rotation: 1.5 ms versus over 8 seconds per likelihood evaluation demonstrates more than three orders of magnitude speed-up for a 6,950-measurement Kepler light curve.The comparison uses celerite against a tuned linear algebra implementation for the same stellar-rotation dataset.
  • Stellar rotation: The stellar-rotation inference measures Prot = 3.80 ± 0.15 days, agreeing well with the literature value and yielding smaller uncertainties.The authors note that detailed comparison with traditional methods is beyond the paper’s scope.
  • Asteroseismic oscillations and transiting planets: The demonstrations extend celerite to physically motivated stellar oscillation models and limb-darkened exoplanet transit light curves.The method evaluates the likelihood for physical parameters in O(N) operations and is applied across five use cases, including real and simulated data.
  • Asteroseismic oscillations: For asteroseismic data, posterior constraints on νmax and ∆ν are consistent with published error bars and tighter than the published constraints.The model requires about 10 CPU minutes for converged MCMC, compared with an estimated 15 CPU hours for a general Cholesky-based solver.

7. COMPARISONS TO OTHER METHODS

The section compares celerite with scalable GP alternatives, emphasizing its exact linear-in-N approach for a restricted one-dimensional kernel class and its practical speed advantage.

  • Scope and trade-offs: celerite is exact but requires a specific stationary kernel and one-dimensional inputs, unlike approximate methods that generally trade flexibility for approximation.Its covariance must be represented by a mixture of exponentials.
  • CARMA comparison: CARMA models share celerite’s one-dimensional setting but use Kalman filtering with O(N J^2) likelihood scaling.Each CARMA model corresponds to an equivalent celerite model, although the inverse parameter mapping is less simple.
  • CARMA comparison: About an order of magnitude: the Kalman filter solver is slower than celerite across all tested systems.The comparison used an optimized C++ implementation of the Kalman filter.
  • Algebraic advantages: celerite applies its factored inverse covariance matrix to general vectors and matrices in O(N), reusing the factorization.This is useful when GP models are combined with linear regression.
  • Other scalable methods: Toeplitz, wavelet, SKI/KISS-GP, and HODLR methods offer broader kernel flexibility or approximations but impose spacing, spectral, precision, or computational constraints.HODLR scales as O(N log^2 N), while Toeplitz methods require evenly spaced data.
  • Empirical comparison: For large J ≳256 and small N ≲1000, george can match or undercut celerite’s cost; elsewhere, celerite is faster and exact in the tested parameter space.This comparison concerns approximate evaluation of celerite models with george.

8. DISCUSSION

The discussion frames celerite as an exact, interpretable approach for extending GP analyses to large astronomical time-series datasets, while identifying important scope limits and future applications.

  • Discussion: For datasets previously limited to N ≲1000, celerite’s linear scaling is intended to extend GP modeling to large and forthcoming astronomical surveys.The method’s oscillator interpretation makes parameters interpretable as resonant frequency, amplitude, and quality factor.
  • Applications: The method has applications in stellar variability, exoplanet analysis, black-hole and quasar variability, reverberation mapping, gravitational lensing, and other one-dimensional problems.The paper also identifies spectroscopy as an example where wavelength replaces time as the independent coordinate.
  • Applications: Mixtures of celerite terms can support probabilistic power-spectrum inference with uneven sampling and heteroscedastic noise, but this application scales quadratically with the number of terms J.The J-scaling remains a limitation shared by existing methods discussed there.
  • Limitations: Celerite cannot speed up multidimensional GP analyses, although extensions for structured cases such as simultaneous parallel time series may be possible.The authors identify the one-dimensional restriction as significant.
  • Code availability: The authors provide a well-tested, documented open-source implementation, with availability under the MIT license alongside the paper’s examples.Implementations and examples are released through GitHub and Zenodo.

A. ENSURING POSITIVE DEFINITENESS FOR CELERITE MODELS

This appendix establishes conditions ensuring that celerite kernels produce positive definite covariance matrices, using power-spectrum positivity and algebraic root checks.

  • Positive definiteness: A stationary GP kernel is valid when its power spectrum is positive at every frequency, equivalently satisfying Bochner’s theorem.Positive definiteness is required for every set of input coordinates.
  • Single-term constraints: For a single celerite term, positive definiteness imposes |b_j d_j| < a_j c_j, with a_j positive and c_j nonnegative.The constraint follows from requiring the power-spectrum numerator to remain positive.
  • Multiple-term constraints: For multiple terms, checking each term separately is sufficient in common cases because products and sums of positive definite kernels remain positive definite.The full general constraint is more complicated than the single-term condition.
  • Root checking: A polynomial of order 2(J −1)+1 results in the general root-checking procedure; with J = 2, the problem reduces to a solvable cubic.Sturm’s theorem determines whether positive real roots exist without explicitly computing them.
  • Root checking: The Sturm-based procedure robustly matches numerically evaluated positive-root counts, and a non-zero real root allows the likelihood to be set to zero.The method uses sign changes at zero and infinity to count roots in (0, ∞].

B. SEMISEPARABLE MATRIX OPERATIONS

This appendix summarizes algorithms for manipulating semiseparable matrices and their use in Gaussian-process models with celerite.

  • Semiseparable matrix operations: The appendix presents semiseparable matrix operations and explains how they support celerite Gaussian-process models.It provides the linear-algebra background for the method’s scalable computations.

B.1. Multiplication

This section derives a numerically stable algorithm for multiplying a celerite covariance matrix using its semiseparable representation. The method uses two sweeps and scales linearly with dataset size for fixed rank.

  • B.1. Multiplication: Rank-R semiseparable matrix products can be computed in O(N R) operations, but naïve application to celerite models can overflow or underflow.The numerical problems resemble those encountered with Cholesky factorization.
  • B.1. Multiplication: The stable celerite multiplication algorithm is derived from the semiseparable representation of K.
  • B.1. Multiplication: Two sweeps compute the forward and backward quantities, with overall scaling O(N J).The sweeps proceed in opposite directions through the data.

B.2. Sampling data from a celerite process

This section develops scalable procedures for sampling data from a celerite Gaussian process. Semiseparable Cholesky structure reduces the relevant computations from general GP costs to O(N J).

  • B.2. Sampling data from a celerite process: Sampling a dataset from a GP with fixed parameters θ and α requires combining the mean function with draws from a unit normal vector.The mean function is evaluated at the input coordinates, and the draws satisfy q_i ∼ N(0, 1).
  • B.2. Sampling data from a celerite process: For a general covariance matrix, computing K_α^1/2 costs O(N^3), while its dot product costs O(N^2).
  • B.2. Sampling data from a celerite process: Celerite uses a semiseparable Cholesky factor to compute the sampling dot product in O(N J).
  • B.2. Sampling data from a celerite process: The scalable sampling algorithm supports generating large simulated datasets and performing posterior predictive checks.

B.3. Interpolation & extrapolation

This section describes scalable GP prediction at a vector of M input coordinates conditioned on observed data and model parameters. Reusing celerite structure improves naïve prediction scaling through forward and backward passes.

  • B.3. Interpolation & extrapolation: Conditioned on data and GP parameters, predictions at new coordinates follow a normal distribution with mean μ_* and covariance K_*.The covariance matrix is computed between the relevant input vectors.
  • B.3. Interpolation & extrapolation: Naïve celerite prediction scales as O(N M) when reusing the Cholesky factor.
  • B.3. Interpolation & extrapolation: The prediction scaling can be improved to O(n N + m M), where n and m are integer constants.The derivation expands the predictive expression using the celerite covariance structure.
  • B.3. Interpolation & extrapolation: The efficient computation uses forward and backward passes over the data and prediction coordinates.The method divides the sum around the prediction location and recursively updates intermediate quantities.
Loading 1703.09710v2…