C++ interface

All of the computationally expensive parts of celerite2 are implemented in C++. This interface was not designed to be as user friendly as the Python interfaces and it isn’t as thoroughly documented, but these details should be enough to get you started if you’re interested in using or contributing to the low-level code directly. Please feel free to open an issue if you have questions!

One thing to note when using these functions is that the performance will be greatly improved (especially for small values of J) by using fixed size matrices because Eigen can work its magic. If you use the Terms interface defined below to construct your low rank matrices, the results will be right (using the Width parameter), but you might need to enforce this manually if you’re building your own matrices. Much of the code in the pybind11 interface is dedicated to making sure that small matrices correctly typed.

It is also worth noting that the notation in the equations below (and throughout these pages) is a little sloppy since it generally neglects the P matrix. Take a look at the original celerite paper for more details about the numerical stability introduced by using this matrix, and how it affects the algorithms.

Basic interface

If you don’t need access to the backpropagation functions, you can use the following functions to compute the standard celerite operations. If you need support for backpropagation, use the functions defined in the Autodiff Interface section, but know that they will require a larger memory footprint.

template<typename Diag, typename LowRank, typename Dense>
void celerite2::core::to_dense(const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<LowRank> &P, Eigen::MatrixBase<Dense> const &K_out)

Get the dense representation of a celerite matrix.

Parameters

  • a: (N,): The diagonal component

  • U: (N, J): The first low rank matrix

  • V: (N, J): The second low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • K_out: (N, N): The dense matrix

template<typename Diag, typename LowRank, typename DiagOut, typename LowRankOut>
Eigen::Index celerite2::core::factor(const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<LowRank> &P, Eigen::MatrixBase<DiagOut> const &d_out, Eigen::MatrixBase<LowRankOut> const &W_out)

Compute the Cholesky factorization of the system.

This computes d and W such that:

diag(a) + tril(U*V^T) + triu(V*U^T) = L*diag(d)*L^T

where

L = 1 + tril(U*W^T)

This can be safely applied in place: d_out can point to a and W_out can point to V, and the memory will be reused.

Parameters

  • a: (N,): The diagonal component

  • U: (N, J): The first low rank matrix

  • V: (N, J): The second low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d_out: (N,): The diagonal component of the Cholesky factor

  • W_out: (N, J): The second low rank component of the Cholesky factor

template<typename Diag, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::solve(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &X_out)

Solve a linear system using the Cholesky factorization.

This computes X in the following linear system:

K * X = Y

where K is the celerite matrix. This uses the results of the Cholesky factorization implemented by celerite2::core::factor.

This can be safely applied in place by setting X_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The right hand side vector or matrix

  • X_out: (N, Nrhs): The solution to the linear system

template<typename Diag, typename LowRank, typename RightHandSide, typename Norm, typename RightHandSideOut>
void celerite2::core::norm(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<Norm> const &norm_out, Eigen::MatrixBase<RightHandSideOut> const &X_out)

Compute the norm of vector or matrix under the celerite metric.

This computes Y^T * K^-1 * Y where K is the celerite matrix. This uses the results of the Cholesky factorization implemented by celerite2::core::factor.

This can be safely applied in place by setting X_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The right hand side vector or matrix

  • norm_out: (Nrhs, Nrhs): The norm of Y

  • X_out: (N, Nrhs): An intermediate result of the operation

template<typename Diag, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::dot_tril(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

Compute product of the Cholesky factor with a vector or matrix.

This computes L * Y where L is the Cholesky factor of a celerite system computed using celerite2::core::factor.

This can be safely applied in place by setting Z_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The target vector or matrix

  • Z_out: (N, Nrhs): The result of the operation

template<typename LowRank, typename RightHandSide, typename Indices, typename RightHandSideOut>
void celerite2::core::conditional_mean(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<RightHandSide> &z, const Eigen::MatrixBase<LowRank> &U_star, const Eigen::MatrixBase<LowRank> &V_star, const Eigen::MatrixBase<Indices> &inds, Eigen::MatrixBase<RightHandSideOut> const &mu_out)

Compute the conditional mean of a celerite process.

This computes K_star * K^-1 * y where K is the celerite matrix and K_star is the rectangular covariance between the training points and the test points.

See the original celerite paper (https://arxiv.org/abs/1703.09710) for the definitions of the _star matrices.

Note that this operation cannot be safely applied in place.

Parameters

  • U: (N, J): The first low rank matrix

  • V: (N, J): The second low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • z: (N,): The solution to the linear system x = K^-1 * y

  • U_star: (M, J): The first low rank matrix for K_star

  • V_star: (M, J): The second low rank matrix for K_star

  • inds: (M,): The indices where the test points would be inserted into the training data; this must be sorted

  • mu_out: (M,): The conditional mean

Autodiff Interface

If you don’t need access to the backpropagation functions, you should use the functions defined in the Basic interface section above. The functions defined in this section track the extra matrices that are required for the reverse pass in reverse-mode automatic differentiation. See Foreman-Mackey (2018) for more information about the reverse pass calculations.

Each of these functions follows the same argument naming conventions. As an example, imagine that we have a function called func that has two input arguments, x (lowercase variables are vectors) and Y (and uppercase variables are matrices), and and two outputs A and b. Let’s say that the reverse pass also requires a cached variable (or “workspace”) called S. In this case, the signature of the forward pass would be something like:

template <typename VectorIn, typename MatrixIn, typename MatrixOut,
          typename VectorOut, typename Workspace>
void func(
  const Eigen::MatrixBase<VectorIn> &x,
  const Eigen::MatrixBase<MatrixIn> &Y,
  Eigen::MatrixBase<MatrixOut> const &A_out,
  Eigen::MatrixBase<VectorOut> const &b_out,
  Eigen::MatrixBase<Workspace> const &S_out
)

Note that, in general, the _out parameters don’t need to have the right shape because they will be resized in place. The reverse pass of this function, would be implemented in another function func_rev with the following signature:

template <typename VectorIn, typename MatrixIn, typename Workspace,
          typename VectorOut, typename MatrixOut>
void func_rev(
  // Original inputs
  const Eigen::MatrixBase<VectorIn> &x,
  const Eigen::MatrixBase<MatrixIn> &Y,
  // Original outputs
  const Eigen::MatrixBase<MatrixIn> &A,
  const Eigen::MatrixBase<VectorIn> &b,
  const Eigen::MatrixBase<Workspace> &S,
  // The sensitivities of the outputs, note: S is not included
  const Eigen::MatrixBase<MatrixIn> &bA,
  const Eigen::MatrixBase<VectorIn> &bb,
  // The (resulting) sensitivities of the inputs
  Eigen::MatrixBase<VectorOut> const &bx_out,
  Eigen::MatrixBase<MatrixOut> const &bY_out
)

where the b prefix before a parameter indicates the overbar from the notation in Foreman-Mackey (2018):

\[\bar{x} = \frac{\partial \mathcal{L}}{\partial x}\]

The forward methods for each celerite operation are documented below, and the reverse passes are all implemented following the rules listed above. Take a look at the source code to see the signatures.

template<bool update_workspace = true, typename Diag, typename LowRank, typename DiagOut, typename LowRankOut, typename Work>
Eigen::Index celerite2::core::factor(const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<LowRank> &P, Eigen::MatrixBase<DiagOut> const &d_out, Eigen::MatrixBase<LowRankOut> const &W_out, Eigen::MatrixBase<Work> const &S_out)

Compute the Cholesky factorization of the system.

This computes d and W such that:

K = L*diag(d)*L^T

where K is the celerite matrix and

L = 1 + tril(U*W^T)

This can be safely applied in place: d_out can point to a and W_out can point to V, and the memory will be reused. In this particular case, the celerite2::core::factor_rev function doesn’t use a and V, but this won’t be true for all _rev functions.

Parameters

  • a: (N,): The diagonal component

  • U: (N, J): The first low rank matrix

  • V: (N, J): The second low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d_out: (N,): The diagonal component of the Cholesky factor

  • W_out: (N, J): The second low rank component of the Cholesky factor

  • S_out: (N, J*J): The cached value of the S matrix at each step

template<bool update_workspace = true, typename Diag, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::solve(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &X_out, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out, Eigen::MatrixBase<Work> const &G_out)

Solve a linear system using the Cholesky factorization.

This computes X in the following linear system:

K * X = Y

where K is the celerite matrix. This uses the results of the Cholesky factorization implemented by celerite2::core::factor.

This can be safely applied in place as long as you don’t need to compute the reverse pass. To compute the solve in place, set X_out = Y and Z_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The right hand side vector or matrix

  • X_out: (N, Nrhs): The final solution to the linear system

  • Z_out: (N, Nrhs): An intermediate result of the operation

  • F_out: (N, J*Nrhs): The workspace for the forward sweep

  • G_out: (N, J*Nrhs): The workspace for the backward sweep

template<bool update_workspace = true, typename Diag, typename LowRank, typename RightHandSide, typename Norm, typename RightHandSideOut, typename Work>
void celerite2::core::norm(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<Norm> const &X_out, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Compute the norm of vector or matrix under the celerite metric.

This computes Y^T * K^-1 * Y where K is the celerite matrix. This uses the results of the Cholesky factorization implemented by celerite2::core::factor.

This can be safely applied in place as long as you don’t need to compute the reverse pass. To compute the solve in place, set Z_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The target vector or matrix

  • X_out: (Nrhs, Nrhs): The norm of Y

  • Z_out: (N, Nrhs): An intermediate result of the operation

  • F_out: (N, J*Nrhs): The workspace for the forward sweep

template<bool update_workspace = true, typename Diag, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::dot_tril(const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<Diag> &d, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Compute product of the Cholesky factor with a vector or matrix.

This computes L * Y where L is the Cholesky factor of a celerite system computed using celerite2::core::factor.

This can be safely applied in place as long as you don’t need to compute the reverse pass. To compute the solve in place, set Z_out = Y.

Parameters

  • U: (N, J): The first low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • d: (N,): The diagonal component of the Cholesky factor

  • W: (N, J): The second low rank component of the Cholesky factor

  • Y: (N, Nrhs): The target vector or matrix

  • Z_out: (N, Nrhs): The result of the operation

  • F_out: (N, J*Nrhs): The workspace for the forward sweep

template<bool update_workspace = true, typename Diag, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::matmul(const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<LowRank> &P, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &X_out, Eigen::MatrixBase<RightHandSideOut> const &M_out, Eigen::MatrixBase<Work> const &F_out, Eigen::MatrixBase<Work> const &G_out)

Compute a matrix-vector or matrix-matrix product.

This computes X = K * Y where K is the celerite matrix.

Note that this operation cannot be safely applied in place.

Parameters

  • a: (N,): The diagonal component

  • U: (N, J): The first low rank matrix

  • V: (N, J): The second low rank matrix

  • P: (N-1, J): The exponential difference matrix

  • Y: (N, Nrhs): The target vector or matrix

  • X_out: (N, Nrhs): The result of the operation

  • M_out: (N, Nrhs): The intermediate state of the system

  • F_out: (N, J*Nrhs): The workspace for the forward sweep

  • G_out: (N, J*Nrhs): The workspace for the backward sweep

Terms

The only the most basic terms are implemented in C++ and they are mostly used for testing purposes, but it would be possible to generalize them to other use cases.

template<typename T, int J_ = Eigen::Dynamic>
class celerite2::Term

The abstract base class from which terms should inherit

Public Types

typedef T Scalar

The underlying scalar type of this Term (should probably always be double)

typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector

An Eigen vector with data type Scalar

typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Width, Order> LowRank

The Eigen type for the low-rank matrices used internally

typedef std::tuple<Vector, Vector, Vector, Vector, Vector, Vector> Coeffs

A tuple of vectors giving the coefficients for the celerite model

typedef std::tuple<Vector, LowRank, LowRank, LowRank> Matrices

A tuple of matrices representing this celerite process

Public Functions

Term()
int get_width() const
void set_coefficients(const Vector &ar, const Vector &cr, const Vector &ac, const Vector &bc, const Vector &cc, const Vector &dc)

Set the coefficients of the term

Parameters

  • ar: (J_real,): The real amplitudes.

  • cr: (J_real,): The real exponential.

  • ac: (J_comp,): The complex even amplitude.

  • bc: (J_comp,): The complex odd amplitude.

  • cc: (J_comp,): The complex exponential.

  • dc: (J_comp,): The complex frequency.

Coeffs get_coefficients() const

Get the coefficients of the term as a tuple

Matrices get_celerite_matrices(const Vector &x, const Vector &diag) const

Get the matrices required to represent the celerite process

Parameters

  • x: (N,): The independent coordinates of the data.

  • diag: (N,): The diagonal variance of the process.

template<typename Other>
Term<typename std::common_type<Scalar, typename Other::Scalar>::type, sum_width<Width, Other::Width>::value> operator+(const Other &other) const

Adding two terms builds a new term where the coefficients have been concatenated

Parameters

  • other: (Term): The term to add to this one.

template<typename T>
class celerite2::SHOTerm : public celerite2::Term<T, 2>

A term representing a stochastically-driven, damped harmonic oscillator

Parameters

  • S0: The power at omega = 0.

  • w0: The undamped angular frequency.

  • Q: The quality factor.

  • eps: A regularization parameter used for numerical stability.

Public Types

typedef T Scalar

The underlying scalar type of this Term (should probably always be double)

Public Functions

SHOTerm(const Scalar &S0, const Scalar &w0, const Scalar &Q, const Scalar &eps = 1e-5)

Public Static Attributes

constexpr int Width = 2
template<typename T>
class celerite2::RealTerm : public celerite2::Term<T, 1>

The simplest celerite model

Parameters

  • a: The amplitude of the term.

  • c: The exponent of the term.

Public Types

typedef T Scalar

The underlying scalar type of this Term (should probably always be double)

Public Functions

RealTerm(const Scalar &a, const Scalar &c)

Public Static Attributes

constexpr int Width = 1
template<typename T>
class celerite2::ComplexTerm : public celerite2::Term<T, 2>

A general celerite model

Parameters

  • a: The real part of the amplitude.

  • b: The complex part of the amplitude.

  • c: The real part of the exponent.

  • d: The complex part of the exponent.

Public Types

typedef T Scalar

The underlying scalar type of this Term (should probably always be double)

Public Functions

ComplexTerm(const Scalar &a, const Scalar &b, const Scalar &c, const Scalar &d)

Public Static Attributes

constexpr int Width = 2