Kernel Building

The class of kernels that can be used in celerite are sums of terms of the following form:

\[\begin{split}k(\tau) &= \sum_{j=1}^J \frac{1}{2}\left[ (a_j + i\,b_j)\,e^{-(c_j+i\,d_j)\,\tau} + (a_j - i\,b_j)\,e^{-(c_j-i\,d_j)\,\tau} \right] \\ &= \sum_{j=1}^J \frac{1}{2}\left[ \alpha_j\,e^{-\beta_j\,\tau} + {\alpha_j}^*\,e^{-{\beta_j}^*\,\tau} \right]\end{split}\]

At a basic level, the interface involves providing the solver with lists of \(\alpha_j\) and \(\beta_j\). It’s useful to know if \(b_j = 0\) and \(d_j = 0\) for any terms because we can use that fact to gain computational efficiency. Therefore, there are 6 arrays that should be provided to the code:

  • a_real: an array of \(a_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),
  • c_real: an array of \(c_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),
  • a_comp: an array of \(a_j\) for the other terms,
  • b_comp: an array of \(b_j\) for the other terms,
  • c_comp: an array of \(c_j\) for the other terms, and
  • d_comp: an array of \(d_j\) for the other terms.

In practice, this is further abstracted and users should use the provided term objects described below or subclass Term to implement a custom kernel. Subclasses should overload the Term.get_real_coefficients() and Term.get_complex_coefficients() methods to provide custom \(\alpha\) and \(\beta\) arrays.

The following gives a simple example for a custom term of the form:

\[k(\tau) = \frac{a}{2+b}\,e^{-c\,\tau}\left[ \cos\left(\frac{2\,\pi\,\tau}{P}\right) + (1+b) \right]\]

with the parameters log_a, log_b, log_c, and log_P.

# If we use the autograd wrappers around numpy then the gradients
# of the term will be computed automatically.
import autograd.numpy as np

class CustomTerm(terms.Term):
    parameter_names = ("log_a", "log_b", "log_c", "log_P")

    def get_real_coefficients(self, params):
        log_a, log_b, log_c, log_P = params
        b = np.exp(log_b)
        return (
            np.exp(log_a) * (1.0 + b) / (2.0 + b), np.exp(log_c),
        )

    def get_complex_coefficients(self, params):
        log_a, log_b, log_c, log_P = params
        b = np.exp(log_b)
        return (
            np.exp(log_a) / (2.0 + b), 0.0,
            np.exp(log_c), 2*np.pi*np.exp(-log_P),
        )

In this example, all of the returned coefficients are scalars but they can also be returned as arrays. To use this kernel, we would use something like the following:

bounds = dict(log_a=(None, None), log_b=(None, 5.0), log_c=(-1.0, 1.0),
              log_P=(-0.5, 0.5))
kernel = CustomTerm(log_a=0.1, log_b=0.5, log_c=-0.01, log_P=0.0,
                    bounds=bounds)

where bounds is optional.

Finally, terms are combined by adding or multiplying them:

term1 = terms.RealTerm(0.5, -5.0)
term2 = terms.ComplexTerm(1.0, 0.1, 0.6, -0.3)
term3 = terms.SHOTerm(1.0, 0.1, 0.6)

# Adding kernels works
kernel = term1 + term2

# So does multiplying them
kernel = term1 * term2

# Or, any combination of these operations
kernel = term1 * (term2 + term3)

These combined terms are also Term objects and they provide the same interface, concatenating the parameters and coefficients in the correct order.

Since Term objects implement the Modeling Protocol, for the kernel described above, we can do things like the following:

print(kernel.get_parameter_dict())

With the output:

OrderedDict([('terms[0]:log_a', 0.5),
             ('terms[0]:log_c', -5.0),
             ('terms[1]:log_a', 1.0),
             ('terms[1]:log_b', 0.1),
             ('terms[1]:log_c', 0.6),
             ('terms[1]:log_d', -0.2)])

The following are the Term objects provided by celerite.

class celerite.terms.Term(*args, **kwargs)

The abstract base “term” that is the superclass of all other terms

Subclasses should overload the terms.Term.get_real_coefficients() and terms.Term.get_complex_coefficients() methods.

check_parameters()

Check for negative power in the PSD using Sturm’s theorem

Returns:True for valid parameters.
coefficients

All of the coefficient arrays

This property is the concatenation of the results from terms.Term.get_real_coefficients() and terms.Term.get_complex_coefficients() but it will always return a tuple of length 6, even if alpha_complex_imag was omitted from get_complex_coefficients.

Returns:(array[j_real], array[j_real], array[j_complex], array[j_complex], array[j_complex], array[j_complex]): alpha_real, beta_real, alpha_complex_real, alpha_complex_imag, beta_complex_real, and beta_complex_imag as described above.
Raises:ValueError – For invalid dimensions for the coefficients.
get_complex_coefficients(params)

Get the arrays alpha_complex_* and beta_complex_*

This method should be overloaded by subclasses to return the arrays alpha_complex_real, alpha_complex_imag, beta_complex_real, and beta_complex_imag given the current parameter settings. By default, this term is empty.

Returns:(array[j_complex], array[j_complex], array[j_complex], array[j_complex]): alpha_complex_real, alpha_complex_imag, beta_complex_real, and beta_complex_imag as described above. alpha_complex_imag can be omitted and it will be assumed to be zero.
get_psd(omega)

Compute the PSD of the term for an array of angular frequencies

Parameters:omega (array[..]) – An array of frequencies where the PSD should be evaluated.
Returns:The value of the PSD for each omega. This will have the same shape as omega.
get_real_coefficients(params)

Get the arrays alpha_real and beta_real

This method should be overloaded by subclasses to return the arrays alpha_real and beta_real given the current parameter settings. By default, this term is empty.

Returns:alpha_real and beta_real as described above.
Return type:(array[j_real], array[j_real])
get_value(tau)

Compute the value of the term for an array of lags

Parameters:tau (array[..]) – An array of lags where the term should be evaluated.
Returns:The value of the term for each tau. This will have the same shape as tau.
terms

A list of all the terms included in a sum of terms

class celerite.terms.RealTerm(*args, **kwargs)

The simplest celerite term

This term has the form

\[k(\tau) = a_j\,e^{-c_j\,\tau}\]

with the parameters log_a and log_c.

Strictly speaking, for a sum of terms, the parameter a could be allowed to go negative but since it is somewhat subtle to ensure positive definiteness, we require that the amplitude be positive through this interface. Advanced users can build a custom term that has negative coefficients but care should be taken to ensure positivity.

Parameters:
  • log_a (float) – The log of the amplitude of the term.
  • log_c (float) – The log of the exponent of the term.
class celerite.terms.ComplexTerm(*args, **kwargs)

A general celerite term

This term has the form

\[k(\tau) = \frac{1}{2}\,\left[(a_j + b_j)\,e^{-(c_j+d_j)\,\tau} + (a_j - b_j)\,e^{-(c_j-d_j)\,\tau}\right]\]

with the parameters log_a, log_b, log_c, and log_d. The parameter log_b can be omitted and it will be assumed to be zero.

This term will only correspond to a positive definite kernel (on its own) if \(a_j\,c_j \ge b_j\,d_j\) and the log_prior method checks for this constraint.

Parameters:
  • log_a (float) – The log of the real part of amplitude.
  • log_b (float) – The log of the imaginary part of amplitude.
  • log_c (float) – The log of the real part of the exponent.
  • log_d (float) – The log of the imaginary part of exponent.
class celerite.terms.SHOTerm(*args, **kwargs)

A term representing a stochastically-driven, damped harmonic oscillator

The PSD of this term is

\[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}\]

with the parameters log_S0, log_Q, and log_omega0.

Parameters:
  • log_S0 (float) – The log of the parameter \(S_0\).
  • log_Q (float) – The log of the parameter \(Q\).
  • log_omega0 (float) – The log of the parameter \(\omega_0\).
class celerite.terms.Matern32Term(*args, **kwargs)

A term that approximates a Matern-3/2 function

This term is defined as

\[k(\tau) = \sigma^2\,\left[ \left(1+1/\epsilon\right)\,e^{-(1-\epsilon)\sqrt{3}\,\tau/\rho} \left(1-1/\epsilon\right)\,e^{-(1+\epsilon)\sqrt{3}\,\tau/\rho} \right]\]

with the parameters log_sigma and log_rho. The parameter eps controls the quality of the approximation since, in the limit \(\epsilon \to 0\) this becomes the Matern-3/2 function

\[\lim_{\epsilon \to 0} k(\tau) = \sigma^2\,\left(1+ \frac{\sqrt{3}\,\tau}{\rho}\right)\, \exp\left(-\frac{\sqrt{3}\,\tau}{\rho}\right)\]
Parameters:
  • log_sigma (float) – The log of the parameter \(\sigma\).
  • log_rho (float) – The log of the parameter \(\rho\).
  • eps (Optional[float]) – The value of the parameter \(\epsilon\). (default: 0.01)
class celerite.terms.JitterTerm(*args, **kwargs)

A diagonal jitter or “white noise” term

This term has the form

\[k(\tau_{n,m}) = \sigma^2\,\delta_{n,m}\]

with the parameter log_sigma.

Note

This kernel has the special property that it is never used for celerite.GP.predict().

Parameters:log_sigma (float) – The log of the amplitude of the white noise.