State-space Representation of Matérn and Damped Simple Harmonic Oscillator Gaussian Processes

Gaussian processes (GPs) are now widely used in the analysis of astronomical time series. Efficient algorithms have been devised that allow likelihood evaluation in ${ \mathcal O }(n)$ operations, where n is the number of observations (see Foreman-Mackey et al. 2017, and references therein). Linear Gaussian state space models allow likelihood evaluation in ${ \mathcal O }(n)$ operations via the Kalman filter (Kalman 1960). In this note we present analytic linear Gaussian state space model representations for the popular Matérn and Damped Simple Harmonic Oscillator (DSHO) GPs, offering an alternative route to efficient likelihood evaluation for these processes. This representation can offer advantages in terms of numerical stability for the Matérn case (D. Foreman-Mackey 2021, personal communication).

There is a correspondence between a stationary GP over ${\mathbb{R}}$ with kernel $k(\cdot )$ and power spectral density $S(\omega )$ that is a rational function of ω, i.e., $S(\omega )\propto p(\omega )/q(\omega )$ with p and q polynomials, and the prior implied by an order p linear stochastic differential equation given by

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{p}f}{{{dt}}^{p}}+{a}_{1}\displaystyle \frac{{d}^{p-1}f}{{{dt}}^{p-1}}+\cdots +{a}_{p}\displaystyle \frac{{df}}{{dt}}\\ \qquad ={b}_{0}\displaystyle \frac{{dW}}{{dt}}+{b}_{1}\displaystyle \frac{{d}^{2}W}{{{dt}}^{2}}+\cdots +{b}_{q}\displaystyle \frac{{d}^{q+1}W}{{{dt}}^{q+1}},\end{array}\end{eqnarray} \tag{ 1 }$$

where W is a Wiener process (Brownian motion) with variance ${\sigma }_{w}^{2}$ (Hartikainen & Sarkka 2010; Saatçi 2012; Brockwell & Davis 2013). Equation (1) is equivalent to a state space model with the following state equation

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\boldsymbol{Ax}}(t)+{\boldsymbol{L}}{dW}(t)\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{x}}=(f(t),\tfrac{{df}}{{dt}},\ldots ,\tfrac{{d}^{p-1}f}{{{dt}}^{p-1}})$, ${\boldsymbol{L}}={\left(\mathrm{0,0},\ldots ,1\right)}^{\top }$ and

$$\begin{eqnarray*}{\boldsymbol{A}}=\left(\begin{array}{ccccc}0 & 1 & \cdots & 0 & 0\\ \vdots & \ddots & \cdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 & 0\\ -{a}_{p} & -{a}_{p-1} & \cdots & -{a}_{2} & -{a}_{1}\end{array}\right),\end{eqnarray*}$$

and the following measurement equation

$$\begin{eqnarray*}&&y={{\boldsymbol{b}}}^{\top }{\boldsymbol{x}}\end{eqnarray*}$$

where ${\boldsymbol{b}}={({b}_{0},{b}_{1},\ldots ,{b}_{p-1})}^{\top }$ and b_j =0 for $j\gt q$. The observed process y would in the astronomical context used to model a time series y(t) observed at times ${\{{t}_{i}\}}_{i=1}^{n}$ with associated measurement errors ${\{{\sigma }_{i}\}}_{i=1}^{n}$.

If ${\boldsymbol{x}}(t)$ is stationary the process y defined by the above state space model can be shown (e.g., Brockwell 2009) to be a zero-mean Gaussian CARMA($p,q$) model with power spectral density

$$\begin{eqnarray}&&S(\omega )=\displaystyle \frac{{\sigma }_{w}^{2}}{2\pi }{\left|\displaystyle \frac{b(i\omega )}{a(i\omega )}\right|}^{2}\end{eqnarray} \tag{ 3 }$$

and variance ${\boldsymbol{\Sigma }}$ that satisfies

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{\Sigma }}(t)}{{dt}}={\boldsymbol{A}}{\boldsymbol{\Sigma }}(t)+{\boldsymbol{\Sigma }}{{\boldsymbol{A}}}^{\top }+{\sigma }_{w}^{2}{\boldsymbol{L}}{{\boldsymbol{L}}}^{\top }.\end{eqnarray} \tag{ 4 }$$

We can express the process given by Equation (2) as a linear Gaussian state space model observed at a set of input points $\{{\boldsymbol{x}}({t}_{i})\}$

$$\begin{eqnarray}\begin{array}{rcl}p({\boldsymbol{x}}({t}_{i})| {\boldsymbol{x}}({t}_{i-1})) & = & N({{\boldsymbol{\Phi }}}_{i-1}{\boldsymbol{x}}({t}_{i-1}),{{\boldsymbol{Q}}}_{i-1})\\ p(y({t}_{i})| {\boldsymbol{x}}({t}_{i})) & = & N({{\boldsymbol{b}}}^{\top }{\boldsymbol{x}}({t}_{i}),{\sigma }_{n}^{2})\end{array}\end{eqnarray} \tag{ 5 }$$

where ${\delta }_{i}\equiv {t}_{i}-{t}_{i-1}$, ${{\boldsymbol{\Phi }}}_{i-1}\equiv \exp ({\boldsymbol{A}}{\delta }_{i})$ and

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i-1}={\sigma }_{w}^{2}{\int }_{0}^{{\delta }_{i}}{\left({\delta }_{i}-h\right)}^{2}{{\boldsymbol{c}}}_{i-1}{{\boldsymbol{c}}}_{i-1}^{\top }\,{dh},\end{eqnarray} \tag{ 6 }$$

where ${\sigma }_{n}^{2}$ is the measurement error, ${{\boldsymbol{c}}}_{i-1}$ is the last column of the matrix $\exp ({\boldsymbol{A}}({\delta }_{i}-h))$ and ${\boldsymbol{H}}$ is a projection matrix that selects the first component of the vector, which is the process in our state space model. With all the conditional probabilities of our process determined, we can perform inference using a filtering algorithm.

A filtering algorithm computes the conditional probabilities of the states ${{\boldsymbol{x}}}_{k}$ given all previous measurements $\{{{\boldsymbol{y}}}_{1},{{\boldsymbol{y}}}_{2},\ldots ,{{\boldsymbol{y}}}_{k}\}\equiv {{\boldsymbol{y}}}_{1:k}$ and a set of parameters ${\boldsymbol{\theta }}$, i.e., $p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{y}}}_{1:k},{\boldsymbol{\theta }})$. For the linear Gaussian case this can done in ${ \mathcal O }(n)$ operations with a set of recursive equations known as the Kalman filter, which are given by ⁴

$$\begin{eqnarray}\begin{array}{rcl}p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{y}}}_{1:k-1},{\boldsymbol{\theta }}) & = & N({{\boldsymbol{\mu }}}_{k}^{-},{{\boldsymbol{P}}}_{k}^{-})\\ p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{y}}}_{1:k},{\boldsymbol{\theta }}) & = & N({{\boldsymbol{\mu }}}_{k},{{\boldsymbol{P}}}_{k})\\ p({{\boldsymbol{y}}}_{k}| {{\boldsymbol{y}}}_{1:k-1},{\boldsymbol{\theta }}) & = & N({{\boldsymbol{H}}}_{k}{{\boldsymbol{\mu }}}_{k}^{-},{{\boldsymbol{S}}}_{k}).\end{array}\end{eqnarray} \tag{ 7 }$$

The mean and covariance matrices of the distributions above are obtained with the Kalman filter prediction and update steps (note that the matrices can depend on ${\boldsymbol{\theta }}$ although we have not written explicitly that dependence). The prediction step is given by

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\mu }}}_{k}^{-} & = & {{\boldsymbol{\Phi }}}_{k-1}{{\boldsymbol{\mu }}}_{k-1}\\ {{\boldsymbol{P}}}_{k}^{-} & = & {{\boldsymbol{\Phi }}}_{k-1}{{\boldsymbol{P}}}_{k-1}{{\boldsymbol{\Phi }}}_{k-1}^{\top }+{{\boldsymbol{Q}}}_{k-1},\end{array}\end{eqnarray} \tag{ 8 }$$

and the update step by

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{v}}}_{k} & = & {{\boldsymbol{y}}}_{k}-{{\boldsymbol{H}}}_{k}{{\boldsymbol{\mu }}}_{k}^{-}\\ {{\boldsymbol{S}}}_{k} & = & {{\boldsymbol{H}}}_{k}{{\boldsymbol{P}}}_{k}^{-}{{\boldsymbol{H}}}_{k}^{\top }+{{\boldsymbol{R}}}_{k}\\ {{\boldsymbol{K}}}_{k} & = & {{\boldsymbol{P}}}_{k}^{-}{{\boldsymbol{H}}}_{k}^{\top }{{\boldsymbol{S}}}_{k}^{-1}\\ {{\boldsymbol{\mu }}}_{k} & = & {{\boldsymbol{\mu }}}_{k}^{-}+{{\boldsymbol{K}}}_{k}{{\boldsymbol{v}}}_{{\boldsymbol{k}}}\\ {{\boldsymbol{P}}}_{k} & = & {{\boldsymbol{P}}}_{k}^{-}-{{\boldsymbol{K}}}_{k}{{\boldsymbol{S}}}_{k}{{\boldsymbol{K}}}_{k}^{\top }.\end{array}\end{eqnarray} \tag{ 9 }$$

The recursion is initialized with a prior mean ${{\boldsymbol{\mu }}}_{0}$ and variance ${{\boldsymbol{P}}}_{0}$. To learn the parameters of a GP that has a corresponding state space representation we are interested in the marginal likelihood which can be computed recursively as

$$\begin{eqnarray}&&p({{\boldsymbol{y}}}_{1:n}| {\boldsymbol{\theta }})=\displaystyle \prod _{k=1}^{n}p({{\boldsymbol{y}}}_{k}| {{\boldsymbol{y}}}_{1:k-1},{\boldsymbol{\theta }}),\end{eqnarray} \tag{ 10 }$$

where we define $p({{\boldsymbol{y}}}_{1}| {{\boldsymbol{y}}}_{0:1})\equiv p({{\boldsymbol{y}}}_{1})$. The terms in the product are part of the output of the Kalman filter (see Equation (7)). Thus, the marginal likelihood can be calculated in ${ \mathcal O }(n)$ operations. Armed with the likelihood, we can sample the posterior using a variety of techniques. We note that the use of the Kalman filter to evaluate the likelihood of a CARMA($p,q$) model in the astronomical context has been presented before in Kelly et al. (2014). In their work the likelihood for a CARMA($p,q$) process is evaluated numerically using the algorithm of Jones & Ackerson (1990) which relies on the diagonalization of the ${\boldsymbol{A}}$ matrix using the roots of the a(z) polynomial. In this note we provide analytic expressions for the ${\boldsymbol{\Phi }}$ and ${\boldsymbol{Q}}$ matrices for some GP kernels of wide astronomical application.

2.1. The Matérn Family

One of the most widely used kernels is the Matérn family, which are indexed by ν and have an associated spectral density for a one-dimensional input given by Rasmussen & Williams (2006)

$$\begin{eqnarray}&&S(\omega )={\sigma }_{f}^{2}\displaystyle \frac{2{\pi }^{1/2}{\rm{\Gamma }}(\nu +1/2){\lambda }^{2\nu }}{{\rm{\Gamma }}(\nu ){\left({\lambda }^{2}+{\omega }^{2}\right)}^{\nu +1/2}}\end{eqnarray} \tag{ 11 }$$

and is thus of the form of Equation (3) for $\nu =n+1/2$, $n\in {\mathbb{N}}$. More precisely, we have $a{(z)=(\lambda +z)}^{n+1}$ and $b(z)=1$. As an example, for the popular $\nu =3/2$ kernel the corresponding stochastic differential equation is

$$\begin{eqnarray}\displaystyle \frac{d{\boldsymbol{x}}(t)}{{dt}}=\left(\begin{array}{cc}0 & 1\\ -{\lambda }^{2} & -2\lambda \end{array}\right){\boldsymbol{x}}(t)+\left(\begin{array}{c}0\\ 1\end{array}\right){dW}(t),\end{eqnarray} \tag{ 12 }$$

where ${\boldsymbol{x}}({\boldsymbol{t}})\equiv (x(t),\tfrac{{dx}(t)}{{dt}})$ and the variance of dW is given in terms of ${\sigma }_{f}^{2}$ and λ by

$$\begin{eqnarray*}&&{\sigma }_{w}^{2}={\sigma }_{f}^{2}\displaystyle \frac{4{\pi }^{3/2}{\rm{\Gamma }}(2){\lambda }^{3}}{{\rm{\Gamma }}(3/2)}\end{eqnarray*}$$

in order that the power spectra density of the process corresponds to that of Equation (11) for $\nu =3/2$.

For the Matérn family with $\nu =n+1/2$ the matrices ${{\boldsymbol{\Phi }}}_{i-1}$ and ${{\boldsymbol{Q}}}_{i-1}$ defined in Equation (6) can be computed analytically using the Laplace transform and a symbolic mathematics package. We provide the expressions of these matrices for a set of values of ν below.

2.2. Stochastically Driven Damped Simple Harmonic Oscillator

Foreman-Mackey et al. (2017) proposed the stochastic process given by the following stochastic differential equation (SDE) as a physically motivated process to describe stellar photometric variations

$$\begin{eqnarray}&&\displaystyle \frac{d2f}{{dt}}+\displaystyle \frac{{\omega }_{0}}{Q}\displaystyle \frac{{df}}{{dt}}+{\omega }_{0}^{2}f(t)={dW}(t).\end{eqnarray} \tag{ 13 }$$

The corresponding power spectral density is given by

$$\begin{eqnarray}&&S(\omega )=\displaystyle \frac{{\sigma }_{w}^{2}}{{\left({\omega }^{2}-{\omega }_{0}^{2}\right)}^{2}+{\omega }^{2}{\omega }_{0}^{2}/{Q}^{2}}\end{eqnarray} \tag{ 14 }$$

where ${\sigma }_{w}^{2}$ is the spectral density of the white noise process dW. In terms of the power at ${\omega }_{0}$ we have that ${\sigma }_{w}^{2}=S({\omega }_{0}){\omega }_{0}^{4}/{Q}^{2}$, and in terms of the total power ${S}_{\mathrm{tot}}$ we have ⁵ that ${\sigma }_{w}^{2}=2{S}_{\mathrm{tot}}{\omega }_{0}^{3}/Q$.

From Equation (13) we see that the DSHO is equivalent to a CARMA(2,0) process with $a(z)={w}_{0}^{2}+({w}_{0}/Q)z+{z}^{2}$ and $b(z)=1$. The corresponding vectorial Markovian SDE is

$$\begin{eqnarray}\displaystyle \frac{d{\boldsymbol{x}}(t)}{{dt}}=\left(\begin{array}{cc}0 & 1\\ -{\omega }_{0}^{2} & -{\omega }_{0}/Q\end{array}\right){\boldsymbol{x}}(t)+\left(\begin{array}{c}0\\ 1\end{array}\right){dW}(t).\end{eqnarray} \tag{ 15 }$$

where ${\boldsymbol{x}}({\boldsymbol{t}})\equiv (x(t),\tfrac{{dx}(t)}{{dt}})$. Just as the case of the Matérn family, the matrices ${{\boldsymbol{\Phi }}}_{i-1}$ and ${{\boldsymbol{Q}}}_{i-1}$ defined in Equation (6) can be computed analytically. Note that when $Q=1/2$ this process reduces to a Matérn 3/2 process with $\lambda ={\omega }_{0}$.

We follow the procedure described in Saatçi (2012) to calculate analytic expressions for the matrices needed for the likelihood calculation for the Matérn and DSHO processes. The matrices ${\boldsymbol{\Phi }}$ are calculated using the identity

$$\begin{eqnarray}&&\exp ({\boldsymbol{A}}t)={{ \mathcal L }}^{-1}\{{\left(s{\boldsymbol{I}}-{\boldsymbol{A}}\right)}^{-1}\},\end{eqnarray} \tag{ 16 }$$

where ${{ \mathcal L }}^{-1}$ denotes the inverse Laplace transform and ${\boldsymbol{I}}$ is the identity matrix. In the case of the Matérn family of processes, the ${\boldsymbol{Q}}$ matrices can be computed analytically via the integral given in Equation (6). In the case of the DSHO, we use the matrix fraction decomposition. If we define matrices ${\boldsymbol{C}}$ and ${\boldsymbol{D}}$ such that ${\boldsymbol{\Sigma }}={\boldsymbol{C}}{{\boldsymbol{D}}}^{-1}$, then Equation (4) is satisfied when ${\boldsymbol{C}}$ and ${\boldsymbol{D}}$ take the form

$$\begin{eqnarray}&&\left(\begin{array}{c}{\boldsymbol{C}}(t)\\ {\boldsymbol{D}}(t)\end{array}\right)=\exp \left\{\left(\begin{array}{cc}{\boldsymbol{A}} & {\boldsymbol{L}}{{\boldsymbol{L}}}^{\top }\\ 0 & -{{\boldsymbol{A}}}^{\top }\end{array}\right)t\right\}\left(\begin{array}{c}{\boldsymbol{C}}(0)\\ {\boldsymbol{D}}(0)\end{array}\right),\end{eqnarray} \tag{ 17 }$$

and ${\boldsymbol{\Sigma }}{(t)={\boldsymbol{C}}(t){\boldsymbol{D}}(t)}^{-1}$. To calculate the ${\boldsymbol{Q}}$ matrices for the DSHO we use Equation (17) with ${\boldsymbol{C}}(0)=0$ and ${\boldsymbol{D}}(0)={\boldsymbol{I}}$, calculate the matrix exponential using the inverse Laplace transform, and then set ${{\boldsymbol{Q}}}_{i-1}={\boldsymbol{\Sigma }}({\delta }_{i})$. Note that we set ${\sigma }_{w}^{2}=1$ to reduce the notational clutter, ${\boldsymbol{Q}}$ matrices presented here should be multiplied by the actual value of this variance when using them. We also denote ${\omega }_{0}$ for the damped simple harmonic oscillator simply by ω.

3.1. Matérn $\nu =1/2$ (Exponential or Ohrstein–Uhlenbeck Process)

$$\begin{eqnarray}&&{{\boldsymbol{\Phi }}}_{i-1}=\exp (-\lambda {\delta }_{i})\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i-1}=-\displaystyle \frac{\exp (-2{\delta }_{i}\lambda )-1}{2\lambda }\end{eqnarray} \tag{ 19 }$$

3.2. Matérn $\nu =3/2$

$$\begin{eqnarray}{{\boldsymbol{\Phi }}}_{i-1}=\exp (-\lambda {\delta }_{i})\left(\begin{array}{cc}\lambda {\delta }_{i}+1 & {\delta }_{i}\\ -{\lambda }^{2}{\delta }_{i} & (1-\lambda {\delta }_{i})\end{array}\right)\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{i-1}=\displaystyle \frac{1}{4{\lambda }^{2}}\left(\begin{array}{cc}{\lambda }^{-1}[1-\exp (-2\delta \lambda )(2{\delta }^{2}{\lambda }^{2}+2\delta \lambda +1)] & 2{\delta }^{2}{\lambda }^{2}\exp (-2\delta \lambda )\\ 2{\delta }^{2}{\lambda }^{2}\exp (-2\delta \lambda ) & \lambda [1-\exp (-2\delta \lambda )(2{\delta }^{2}{\lambda }^{2}-2\delta \lambda +1)]\end{array}\right)\end{eqnarray} \tag{ 21 }$$

3.3. Matérn $\nu =5/2$

$$\begin{eqnarray}{{\boldsymbol{\Phi }}}_{i-1}=\exp (-\lambda {\delta }_{i})\left(\begin{array}{ccc}\tfrac{\left({\lambda }^{2}{\delta }^{2}+2\lambda \delta +2\right)}{2} & \delta \left(\lambda \delta +1\right) & \tfrac{{\delta }^{2}}{2}\\ -\tfrac{{\lambda }^{3}{\delta }^{2}}{2} & \left(-{\lambda }^{2}{\delta }^{2}+\lambda \delta +1\right) & -\tfrac{\delta \left(\lambda \delta -2\right)}{2}\\ \tfrac{{\lambda }^{3}\delta \left(\lambda \delta -2\right)}{2} & {\lambda }^{2}\delta \left(\lambda \delta -3\right) & \tfrac{\left({\lambda }^{2}{\delta }^{2}-4\lambda \delta +2\right)}{2}\end{array}\right)\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{i-1}=\left(\begin{array}{ccc}\tfrac{3}{16\,{\lambda }^{5}} & 0 & \tfrac{1}{16\,{\lambda }^{3}}\\ 0 & \tfrac{1}{16\,{\lambda }^{3}} & 0\\ \tfrac{1}{16\,{\lambda }^{3}} & 0 & \tfrac{3}{16\,\lambda }\end{array}\right)+\exp (-\lambda {\delta }_{i})\times \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}\left(\begin{array}{ccc}-\tfrac{\left(2\,{\delta }^{4}\,{\lambda }^{4}+4\,{\delta }^{3}\,{\lambda }^{3}+6\,{\delta }^{2}\,{\lambda }^{2}+6\,\delta \,\lambda +3\right)}{16\,{\lambda }^{5}} & \tfrac{{\delta }^{4}}{8} & \tfrac{\left(2\,\delta \,\lambda +2\,{\delta }^{2}\,{\lambda }^{2}+4\,{\delta }^{3}\,{\lambda }^{3}-2\,{\delta }^{4}\,{\lambda }^{4}+1\right)}{16\,{\lambda }^{3}}\\ \tfrac{{\delta }^{4}}{8} & -\tfrac{\left(2\,{\delta }^{4}\,{\lambda }^{4}-4\,{\delta }^{3}\,{\lambda }^{3}+2\,{\delta }^{2}\,{\lambda }^{2}+2\,\delta \,\lambda +1\right)}{16\,{\lambda }^{3}} & \tfrac{{\delta }^{2}{\left(\delta \lambda -2\right)}^{2}}{8}\\ \tfrac{\left(2\,\delta \,\lambda +2\,{\delta }^{2}\,{\lambda }^{2}+4\,{\delta }^{3}\,{\lambda }^{3}-2\,{\delta }^{4}\,{\lambda }^{4}+1\right)}{16\,{\lambda }^{3}} & \tfrac{{\delta }^{2}{\left(\delta \lambda -2\right)}^{2}}{8} & -\tfrac{\left(2\,{\delta }^{4}\,{\lambda }^{4}-12\,{\delta }^{3}\,{\lambda }^{3}+22\,{\delta }^{2}\,{\lambda }^{2}-10\,\delta \,\lambda +3\right)}{16\,\lambda }\end{array}\right)\end{eqnarray} \tag{ 24 }$$

3.4. DSHO with Parameters ω and Q. Case $Q\geqslant 1/2$.

Define $\beta =\sqrt{{Q}^{2}-1/4}$.

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{\Phi }}}_{i-1}=\exp (-\omega {\delta }_{i}/2Q)\\ \times \,\left(\begin{array}{cc}\cos (\omega \beta {\delta }_{i}/Q)+\sin (\omega \beta {\delta }_{i}/Q)/2\beta & Q\sin (\omega \beta {\delta }_{i}/Q)/\omega \beta \\ -Q\omega \sin (\omega {\delta }_{i}\beta /Q)/\beta & \cos (\omega {\delta }_{i}\beta /Q)-\sin (\omega {\delta }_{i}\beta /Q)/2\beta \end{array}\right)\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i-1}=Q\exp (-\omega {\delta }_{i}/Q)\times \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\left(\begin{array}{cc}\tfrac{(\cos (2\omega {\delta }_{i}\beta /Q)-1)-2\beta \sin (2\omega {\delta }_{i}\beta /Q)+4{\beta }^{2}(\exp (\omega {\delta }_{i}/Q)-1)}{8{\omega }^{3}{\beta }^{2}} & \tfrac{Q{\sin }^{2}(\omega {\delta }_{i}\beta /Q)}{2{\omega }^{2}{\beta }^{2}}\\ \tfrac{Q{\sin }^{2}(\omega {\delta }_{i}\beta /Q)}{2{\omega }^{2}{\beta }^{2}} & \tfrac{(\cos (2\omega {\delta }_{i}\beta /Q)-1)+2\beta \sin (2\omega {\delta }_{i}\beta /Q)+4{\beta }^{2}(\exp (\omega {\delta }_{i}/Q)-1)}{8\omega {\beta }^{2}}\end{array}\right)\end{eqnarray} \tag{ 27 }$$

3.5. DSHO with Parameters ω and Q. Case $Q\lt 1/2$.

Define $\beta =\sqrt{1/4-{Q}^{2}}$. The expressions for ${\boldsymbol{\Phi }}$ and ${\boldsymbol{Q}}$ correspond to those of the case $Q\geqslant 1/2$ using the fact that $\cos ({ix})=\cosh (x)$ and $\sin ({ix})=i\sinh (x)$.

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{\Phi }}}_{i-1}=\exp (-\omega {\delta }_{i}/2Q)\\ \ \times \,\left(\begin{array}{cc}\cosh (\omega \beta {\delta }_{i}/Q)+\sinh (\omega \beta {\delta }_{i}/Q)/2\beta & Q\sinh (\omega \beta {\delta }_{i}/Q)/\omega \beta \\ -Q\omega \sinh (\omega {\delta }_{i}\beta /Q)/\beta & \cosh (\omega {\delta }_{i}\beta /Q)-\sinh (\omega {\delta }_{i}\beta /Q)/2\beta \end{array}\right)\end{array}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{i-1}=Q\exp (-\omega {\delta }_{i}/Q)\times \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}\left(\begin{array}{cc}\tfrac{(1-\cosh (2\omega {\delta }_{i}\beta /Q))-2\beta \sinh (2\omega {\delta }_{i}\beta /Q)+4{\beta }^{2}(\exp (\omega {\delta }_{i}/Q)-1)}{8{\omega }^{3}{\beta }^{2}} & \tfrac{Q{\sinh }^{2}(\omega {\delta }_{i}\beta /Q)}{2{\omega }^{2}{\beta }^{2}}\\ \tfrac{Q{\sinh }^{2}(\omega {\delta }_{i}\beta /Q)}{2{\omega }^{2}{\beta }^{2}} & \tfrac{(1-\cosh (2\omega {\delta }_{i}\beta /Q))+2\beta \sinh (2\omega {\delta }_{i}\beta /Q)+4{\beta }^{2}(\exp (\omega {\delta }_{i}/Q)-1)}{8\omega {\beta }^{2}}\end{array}\right)\end{eqnarray} \tag{ 30 }$$

Note that the ${\boldsymbol{\Phi }}$ and ${\boldsymbol{Q}}$ matrices for the DSHO process have well-defined limits as $Q\to {0.5}^{+}$ and $Q\to {0.5}^{-}$, and are given by the corresponding matrices for the Matérn $\nu =3/2$ process. This follows from the fact that ${\mathrm{lim}}_{x\to 0}\sin (x)/x={\mathrm{lim}}_{x\to 0}\sinh (x)/x\,=\,1$.