state space gaussian process

Gaussian processes (GPs) scale as $O(N^3)$, where $N$ is the number of data points, which is problematic for large datasets.

By parametrising the GP as a stochastic differential equation, we can reformulate the GP regression problem into a linear gaussian ssm, where it can be solved using kalman filtering and smoothing with linear time complexity ($O(N)$).

See the Temporal Gaussian Process Regression in Logarithmic Time and Kalman filtering and smoothing solutions to temporal Gaussian process regression modelspapers, as well as Adrien Corenflos’ implementation, the BayesNewton implementationand the GPy implementation. Perhaps could be implemented in dynamax by wrapping the filter step within a larger log-likelihood to optimise the lengthscale and variance, but no success so far.

Mátern-5/2 kernel

This is the case $p=2$, so $\nu =p+1/2=5/2$, a typical choice for spatial and temporal problems.

We need to train lengthscale ($\ell$) and variance ($\alpha$) of kernel during fitting, as well as the likelihood noise $\sigma$.

$$\begin{array}{r}\underset{z_t}{\underbrace{\left(\begin{array}{c}{z}_t^1\\ z_t^2\\ z_t^3\end{array}\right)}}=\underset{F}{\underbrace{\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -{\lambda }^3& -3{\lambda }^2& -3\lambda \end{array}\right)}}\underset{z_{t-1}}{\underbrace{\left(\begin{array}{c}{z}_{t-1}^1\\ z_{t-1}^2\\ z_{t-1}^3\end{array}\right)}}+q_t\end{array}$$

where

$$\begin{array}{r}\lambda =\sqrt{5}/\ell \\ q_t=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)w(t)\end{array}$$

and $w(t)$ is some white noise process.

The observation model is

$$\begin{array}{rl}{y}_t& =\underset{H}{\underbrace{\left(\begin{array}{ccc}1& 0& 0\end{array}\right)}}\underset{z_t}{\underbrace{\left(\begin{array}{c}{z}_t^1\\ z_t^2\\ z_t^3\end{array}\right)}}+r_t\end{array}$$

where $r_t\sim \mathcal{N}(0,{\sigma }^2)$.

This is a linear gaussian ssm.

Matérn-1/2 kernel

Matérn-1/2 is a less common covariance function, but it results in even simpler matrices:

$$\begin{array}{rl}\underset{z_t}{\underbrace{\left(\begin{array}{c}{z}_t\end{array}\right)}}& =\underset{F}{\underbrace{\left(\begin{array}{c}-\lambda \end{array}\right)}}\underset{z_{t-1}}{\underbrace{\left(\begin{array}{c}{z}_{t-1}\end{array}\right)}}+q_t\\ \lambda & =\sqrt{5}/\ell \\ q_t& =w(t)\end{array}$$

with the observation model

$$\begin{array}{rl}{y}_t& =\underset{H}{\underbrace{\left(\begin{array}{c}1\end{array}\right)}}\underset{z_t}{\underbrace{\left(\begin{array}{c}{z}_t\end{array}\right)}}+r_t\end{array}$$

where $r_t\sim \mathcal{N}(0,{\sigma }^2)$.