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 Input, typename Coeffs, typename Diag, typename LowRank, typename Dense>
void celerite2::core::to_dense(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, Eigen::MatrixBase<Dense> const &K_out)

Get the dense representation of a celerite matrix.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

  • a – (N,): The diagonal component

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

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

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

template<typename Input, typename Coeffs, typename Diag, typename LowRank, typename DiagOut, typename LowRankOut>
Eigen::Index celerite2::core::factor(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, 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

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

  • a – (N,): The diagonal component

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

  • V – (N, J): The second low rank 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 Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::solve_lower(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

Compute the solution of a lower triangular linear equation.

This computes Z such that:

Y = L * Y

where

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

This can be safely applied in place.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The right hand side

  • Z_out – (N, Nrhs): The solution of this equation

template<typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::solve_upper(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

Compute the solution of a upper triangular linear equation.

This computes Z such that:

Y = L^T * Y

where

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

This can be safely applied in place.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The right hand side

  • Z_out – (N, Nrhs): The solution of this equation

template<typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::matmul_lower(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

Apply a strictly lower matrix multiply.

This computes:

Z += tril(U * V^T) * Y

where tril is the strictly lower triangular function.

Note that this will update the value of Z.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The matrix to be multiplied

  • Z_out – (N, Nrhs): The matrix to be updated

template<typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::matmul_upper(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

Apply a strictly upper matrix multiply.

This computes:

Z += triu(V * U^T) * Y

where triu is the strictly lower triangular function.

Note that this will update the value of Z.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The matrix to be multiplied

  • Z_out – (N, Nrhs): The matrix to be updated

template<typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::general_matmul_lower(const Eigen::MatrixBase<Input> &t1, const Eigen::MatrixBase<Input> &t2, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

The general lower-triangular dot product of a rectangular celerite system.

Parameters

  • t1 – (N,): The left input coordinates (must be sorted)

  • t2 – (M,): The right input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (M, Nrhs): The matrix that will be multiplied

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

template<typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut>
void celerite2::core::general_matmul_upper(const Eigen::MatrixBase<Input> &t1, const Eigen::MatrixBase<Input> &t2, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out)

The general upper-triangular dot product of a rectangular celerite system.

Parameters

  • t1 – (N,): The left input coordinates (must be sorted)

  • t2 – (M,): The right input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (M, Nrhs): The matrix that will be multiplied

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

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 Input, typename Coeffs, typename Diag, typename LowRank, typename DiagOut, typename LowRankOut, typename Work>
Eigen::Index celerite2::core::factor(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<Diag> &a, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, 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

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

  • a – (N,): The diagonal component

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

  • V – (N, J): The second low rank 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 Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::solve_lower(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Apply a strictly lower matrix multiply.

This computes:

Z += tril(U * V^T) * Y

where tril is the strictly lower triangular function.

Note that this will update the value of Z.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The matrix to be multiplied

  • Z_out – (N, Nrhs): The matrix to be updated

  • F_out – (N, J*Nrhs): The workspace

template<bool update_workspace = true, typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::solve_upper(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &W, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Compute the solution of a upper triangular linear equation.

This computes Z such that:

Y = L^T * Y

where

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

This can be safely applied in place.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The right hand side

  • Z_out – (N, Nrhs): The solution of this equation

  • F_out – (N, J*Nrhs): The workspace

template<bool update_workspace = true, typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::matmul_lower(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Apply a strictly lower matrix multiply.

This computes:

Z += tril(U * V^T) * Y

where tril is the strictly lower triangular function.

Note that this will update the value of Z.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The matrix to be multiplied

  • Z_out – (N, Nrhs): The matrix to be updated

  • F_out – (N, J*Nrhs): The workspace

template<bool update_workspace = true, typename Input, typename Coeffs, typename LowRank, typename RightHandSide, typename RightHandSideOut, typename Work>
void celerite2::core::matmul_upper(const Eigen::MatrixBase<Input> &t, const Eigen::MatrixBase<Coeffs> &c, const Eigen::MatrixBase<LowRank> &U, const Eigen::MatrixBase<LowRank> &V, const Eigen::MatrixBase<RightHandSide> &Y, Eigen::MatrixBase<RightHandSideOut> const &Z_out, Eigen::MatrixBase<Work> const &F_out)

Apply a strictly upper matrix multiply.

This computes:

Z += triu(V * U^T) * Y

where triu is the strictly lower triangular function.

Note that this will update the value of Z.

Parameters

  • t – (N,): The input coordinates (must be sorted)

  • c – (J,): The transport coefficients

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

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

  • Y – (N, Nrhs): The matrix to be multiplied

  • Z_out – (N, Nrhs): The matrix to be updated

  • F_out – (N, J*Nrhs): The workspace

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 Eigen::Matrix<Scalar, Width, 1> CoeffVector

The Eigen type for a fixed width vector of coefficients

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

A tuple of vectors giving the coefficients for the celerite model

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

A tuple of matrices representing this celerite process

Public Functions

inline Term()
inline int get_width() const
inline 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.

inline Coeffs get_coefficients() const

Get the coefficients of the term as a tuple

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

Param S0

The power at omega = 0.

Param w0

The undamped angular frequency.

Param Q

The quality factor.

Param 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

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

Public Static Attributes

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

The simplest celerite model

Param a

The amplitude of the term.

Param 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

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

Public Static Attributes

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

A general celerite model

Param a

The real part of the amplitude.

Param b

The complex part of the amplitude.

Param c

The real part of the exponent.

Param 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

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

Public Static Attributes

static constexpr int Width = 2