marginparsep has been altered.
topmargin has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Free-Form Variational Inference for Gaussian Process State-Space Models

Xuhui Fan 1  Edwin V. Bonilla 2  Terence J. O’Kane 3  Scott A. Sisson 4

††footnotetext:  1University of Newcastle, Australia 2CSIRO’s Data61,
Australia 3CSIRO’s Environment, Australia 4University of New South Wales,
Australia. Correspondence to: Xuhui Fan<EMAIL_ADDRESS>
Preprint. Copyright 2023 by the author(s).

###### Abstract

Gaussian process state-space models (GPSSMs) provide a principled and flexible
approach to modeling the dynamics of a latent state, which is observed at
discrete-time points via a likelihood model. However, inference in GPSSMs is
computationally and statistically challenging due to the large number of
latent variables in the model and the strong temporal dependencies between
them. In this paper, we propose a new method for inference in Bayesian GPSSMs,
which overcomes the drawbacks of previous approaches, namely over-simplified
assumptions, and high computational requirements. Our method is based on free-
form variational inference via stochastic gradient Hamiltonian Monte Carlo
within the inducing-variable formalism. Furthermore, by exploiting our
proposed variational distribution, we provide a collapsed extension of our
method where the inducing variables are marginalized analytically. We also
showcase results when combining our framework with particle MCMC methods. We
show that, on six real-world datasets, our approach can learn transition
dynamics and latent states more accurately than competing methods.

##  1 Introduction

State-space models (ssms; Murphy, 2023, Ch. 29) characterize the underlying
dynamics of a latent state given a set of observations via a transition (or
evolution) function and an observation model. As a modeling framework, they
provide a general approach for understanding time-series data (Kitagawa, 1987)
and for data assimilation problems (Katzfuss et al., 2016). Applications of
ssms abound and span diverse areas such as econometrics (Tsay, 2005),
meteorology (Hernandez et al., 1991), control engineering (Ogata et al., 2010)
and neuroscience (Brown et al., 1998).

In this paper, we focus on Bayesian ssms, where the transition function
describing the dynamics of the system is given a prior distribution. A
paradigmatic example of Bayesian ssms are Gaussian process state-space models
(gpssms; Frigola, 2015), where this prior distribution is a Gaussian process
(gp; Williams & Rasmussen, 2006). Due to their Bayesian non-parametric nature,
gpssms represent a principled and flexible approach to Bayesian ssms.

However, the flexibility of Gaussian processes (gps) adds significant
computational and statistical challenges to the already difficult problem of
inference in Bayesian ssms. Indeed, even for non-Bayesian ssms, standard
problems such as filtering, smoothing, and prediction are, in general,
analytically intractable111With the notable exception of linear state
transition and observation models with additive Gaussian noise, for which the
solution to the optimal filtering problem is given by the Kalman filter
(Kalman, 1960).. Having a gp prior over the transition function in ssms
increases the number of latent variables significantly; incorporates strong
(and potentially long-term) dependencies across states; and introduces a cubic
time complexity as a function of the number of observations.

Within the gp community, significant advances have been made addressing the
computational issues in gp regression and classification problems, most
notably using inducing-variable approximations Titsias (2009); Hensman et al.
(2015); Rossi et al. (2021) but also random-feature expansions Cutajar et al.
(2017); Marmin & Filippone (2022) and, more recent innovative approaches such
as the Vecchia approximation Sauer et al. (2022). Although these
approximations are applicable to gp-based dynamic models, the challenges above
remain prevalent within the context of gpssms.

Nevertheless, previous approaches have developed insightful and practical
algorithms for inference in gpssms, mainly based on inducing-variable
approximations, which is also our main underpinning methodology for scalable
gps. In outlining the most relevant approaches, our main object of interest is
the approximate joint posterior over state trajectories
${\boldsymbol{\mathbf{x}}}_{0:T}$ and inducing variables
${\boldsymbol{\mathbf{u}}}$,
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$, where $T+1$ is
the length of the trajectory. We consider two main aspects of this joint
distribution: (i) whether the dependencies between state trajectories and
inducing variables are captured and (ii) whether their corresponding
distributions are unconstrained, i.e., not restricted to a sub-optimal
parametric form. The seminal variational Gaussian process state space model
(vgpssm) proposed by Frigola et al. (2014) as well as the subsequent
identifiable Gaussian process state space model (igpssm) of Eleftheriadis et
al. (2017) use mean-field approaches, therefore, ignoring the posterior
dependencies between state trajectories and inducing variables. The more
recent methods, namely the probabilistic recurrent state space model (prssm)
of Doerr et al. (2018) and the variationally coupled dynamic trajectories
(vcdt) of Ialongo et al. (2019) introduce couplings across state trajectories
and inducing variables. However, their posteriors are constrained to be
Gaussians. Thus, these previous works have either assumed independence between
state trajectories and inducing variables or imposed strong parametric
constraints in their corresponding posteriors or both.

As shown by Ialongo et al. (2019), a mean-field posterior can yield poor
practical performance. Similarly, a Gaussian assumption on the state posterior
or the inducing variable posterior is also very strong and, by definition,
will not generally capture the true posterior even in the limit of infinite
computation. To address these issues, we propose a free-form variational
inference approach to posterior estimation in gpssms that models the full
joint distribution over states and inducing variables,
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$, without any
mean-field or parametric assumptions. We refer to our method as ffvd and
summarize the major differences between its posterior assumptions and those of
previous approaches in Table 1. Below we describe our contributions in more
detail.

(i) Flexible posterior: We develop an inference algorithm for gpssms based on
stochastic gradient Hamiltonian Monte Carlo (sghmc; Chen et al., 2014; Havasi
et al., 2018), which represents the posterior over states and inducing
variables using samples. Our approach naturally captures the posterior
correlations between states and inducing variables; does not constrain this
posterior to any parametric form, and, is scalable to a large number of
observations.

(ii) Collapsed inference that accelerates convergence: we show that (i) our
formulation allows us to collapse the inducing variables
${\boldsymbol{\mathbf{u}}}$, (ii) sample from the lower-dimensional marginal
$q({\boldsymbol{\mathbf{x}}}_{0:T})$ and, at the end of the sampling
procedure, (iii) obtain samples from the conditional
$q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$, for which
we derive a closed-form expression. We show that collapsing accelerates
convergence significantly.

(iii) Extensions with particle Markov chain Monte Carlo (pmcmc): we further
investigate whether more elaborate inference algorithms, such as pmcmc
(Andrieu et al., 2010), that account for the sequential nature of the problem
can provide more accurate posteriors.

(iv) State-of-the-art performance: we showcase the properties and benefits of
our approach compared to previous methods such as vgpssm, prssm and vcdt in a
synthetic example and six system identification benchmarks. Overall, our
method provides state-of-the-art performance when evaluated on these problems,
while having comparable computational requirements to previous approaches.

Table 1: Comparison across methods in terms of their assumptions on the variational distribution. The rows refer to whether the variational posterior captures the dependencies between the state trajectories and the inducing variables (coupled $q({{\boldsymbol{\mathbf{x}}}}_{0:T},{{\boldsymbol{\mathbf{u}}}})$); the distribution over states is unconstrained (unconstrained $q({{\boldsymbol{\mathbf{x}}}}_{0:T}\,|\,{{\boldsymbol{\mathbf{u}}}})$ or $q({{\boldsymbol{\mathbf{x}}}}_{0:T})$), i.e., not restricted to a sub-optimal parametric form; and whether the distribution over the inducing variables is also unconstrained (unconstrained $q({{\boldsymbol{\mathbf{u}}}})$). Our method is referred to as ffvd.  | vgpssm | igpssm | prssm | vcdt | ffvd
---|---|---|---|---|---
Coupled $q({{\boldsymbol{\mathbf{x}}}}_{0:T},{{\boldsymbol{\mathbf{u}}}})$ | ✗ | ✗ | ✓ | ✓ | ✓
Unconstrained $q({{\boldsymbol{\mathbf{x}}}}_{0:T}\,|\,{{\boldsymbol{\mathbf{u}}}})$ or $q({{\boldsymbol{\mathbf{x}}}}_{0:T}$) | ✓ | ✗ | ✗ | ✗ | ✓
Unconstrained $q({\boldsymbol{\mathbf{u}}})$ | ✓ | ✗ | ✗ | ✗ | ✓

##  2 Inference in Gaussian Process Models

Gaussian processes (gps) are priors over functions where every subset of
function values follow a Gaussian distribution. We use
$f({\boldsymbol{\mathbf{x}}})\sim{\mathcal{GP}}\left(m_{f}({\boldsymbol{\mathbf{x}}}),\kappa_{f}({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}^{\prime};{\boldsymbol{\mathbf{\theta}}})\right)$
to denote that $f$ is distributed according to a gp with mean function
$m_{f}(\cdot)$ and covariance function
$\kappa_{f}(\cdot,\cdot;{\boldsymbol{\mathbf{\theta}}})$, where
${\boldsymbol{\mathbf{\theta}}}$ are referred to as the gp hyper-
parameters222We note here that the mean function can also have additional
parameters but, in this work, we have simply assumed identity mean functions..
By definition, a gp prior over functions implies a finite prior over $T$
function values
${\boldsymbol{\mathbf{f}}}=[f({\boldsymbol{\mathbf{x}}}_{1}),\ldots,f({\boldsymbol{\mathbf{x}}}_{T})]^{\top}$,
i.e.,
${\boldsymbol{\mathbf{f}}}\sim\mathcal{N}({\boldsymbol{\mathbf{f}}};{\boldsymbol{\mathbf{m}}},{\boldsymbol{\mathbf{K}}})$,
where ${\boldsymbol{\mathbf{m}}}$ and ${\boldsymbol{\mathbf{K}}}$ are obtained
by evaluating the mean function and covariance function at all the inputs
${\boldsymbol{\mathbf{X}}}:=\\{{\boldsymbol{\mathbf{x}}}_{1},\ldots,{\boldsymbol{\mathbf{x}}}_{T}\\}$.

In supervised learning settings, we are given input-output observations
$\\{{\boldsymbol{\mathbf{x}}}_{t},y_{t}\\}_{t=1}^{T}$ and a conditional
likelihood model $p({\boldsymbol{\mathbf{y}}}\,|\,{\boldsymbol{\mathbf{f}}})$.
Inference involves estimating the posterior distribution
$p({\boldsymbol{\mathbf{f}}}\,|\,{\boldsymbol{\mathbf{y}}},{\boldsymbol{\mathbf{X}}})$
and the hyper-parameters ${\boldsymbol{\mathbf{\theta}}}$ from data.
Consequently, we can use these to estimate the posterior predictive
distribution at a new point ${\boldsymbol{\mathbf{x}}}_{\star}$, i.e.,
$p(f({\boldsymbol{\mathbf{x}}}_{\star})\,|\,{\boldsymbol{\mathbf{X}}},{\boldsymbol{\mathbf{y}}},{\boldsymbol{\mathbf{x}}}_{\star})$.
Notoriously, these tasks have cubic time complexity as a function of the
number of training observations, i.e., $\mathcal{O}(T^{3})$, arising from
algebraic operations involving the computation of the inverse covariance and
its log determinant. This motivates the need for sparse approximations.

###  2.1 Sparse gp Approximations via Inducing Variables

Various approximations have been proposed to deal with the cubic computational
complexity of inference in gp models. In this paper, we focus on inducing-
variable approximations based on variational inference, as originally proposed
by Titsias (2009) and made scalable to very large datasets by Hensman et al.
(2013). The main idea of these approximations is to augment the space of
function values with a set of $M$ inducing variables
${\boldsymbol{\mathbf{u}}}:=\\{u_{i}\\}$ and their corresponding inducing
inputs ${\boldsymbol{\mathbf{Z}}}:=\\{{\boldsymbol{\mathbf{z}}}_{i}\\}$. Thus,
inference involves estimating the posterior
$q({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}})\approx
p({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}\,|\,{\boldsymbol{\mathbf{X}}},{\boldsymbol{\mathbf{y}}})$
and the gp hyper-parameters ${\boldsymbol{\mathbf{\theta}}}$ via variational
inference. Under the assumption that
$q({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}):=q({\boldsymbol{\mathbf{u}}})p({\boldsymbol{\mathbf{f}}}\,|\,{\boldsymbol{\mathbf{u}}})$
where $p({\boldsymbol{\mathbf{f}}}\,|\,{\boldsymbol{\mathbf{u}}})$ is the
conditional prior, the variational objective, the so-called evidence lower
bound (elbo), decomposes over the observations, and inference can be carried
out with time complexity of $\mathcal{O}(M^{3})$, providing significant
advantages when $M\ll T$.

##  3 Gaussian Process State-Space Models

Let us assume we are given a time series of $T$ multi-dimensional observations
${\boldsymbol{\mathbf{y}}}_{1:T}$ and denote their corresponding latent states
with ${\boldsymbol{\mathbf{x}}}_{0:T}$, where ${\boldsymbol{\mathbf{x}}}_{0}$
is the initial state. Here we denote the time series
${\boldsymbol{\mathbf{x}}}_{t_{1}:t_{2}}:=\\{{\boldsymbol{\mathbf{x}}}_{t_{1}},\ldots,{\boldsymbol{\mathbf{x}}}_{t_{2}}\\}$,
and similarly for ${\boldsymbol{\mathbf{y}}}_{1:T}$. In general,
${\boldsymbol{\mathbf{x}}}_{t}\in\mathbb{R}^{d_{x}}$ and
${\boldsymbol{\mathbf{y}}}_{t}\in\mathbb{R}^{d_{y}}$. Gaussian process state-
space models (gpssms) formulate a discrete-time state-space model (ssm) where
the transition dynamics is given by a gp. The full generative process is:

$\displaystyle{\boldsymbol{\mathbf{x}}}_{0}$ $\displaystyle\sim
p({\boldsymbol{\mathbf{x}}}_{0}),\text{
}f({\boldsymbol{\mathbf{x}}})\sim{\mathcal{GP}}(m_{f}({\boldsymbol{\mathbf{x}}}),\kappa_{f}({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}^{\prime};{\boldsymbol{\mathbf{\theta}}})),$
(1) $\displaystyle{\boldsymbol{\mathbf{f}}}_{t}$
$\displaystyle:=f({\boldsymbol{\mathbf{x}}}_{t-1}),\quad{\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t}\sim\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{f}}}_{t},{\boldsymbol{\mathbf{Q}}}),$
(2)
$\displaystyle{\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t}$
$\displaystyle\sim
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t},{\boldsymbol{\mathbf{\phi}}}),$
(3)

where ${\boldsymbol{\mathbf{Q}}}$ is the transition process covariance, and
${\boldsymbol{\mathbf{\phi}}}$ is the vector of parameters of the conditional
likelihood
$p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t},{\boldsymbol{\mathbf{\phi}}})$.
Although our framework does not make any parametric assumptions about this
conditional likelihood, in our experiments in § 7, we adopt the same setting
as in previous works (Ialongo et al., 2019; Doerr et al., 2018) and set
$p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t},{\boldsymbol{\mathbf{\phi}}}):=\mathcal{N}({\boldsymbol{\mathbf{y}}}_{t};{\boldsymbol{\mathbf{C}}}{\boldsymbol{\mathbf{x}}}_{t}+{\boldsymbol{\mathbf{d}}},{\boldsymbol{\mathbf{R}}})$,
with
${\boldsymbol{\mathbf{\phi}}}=\\{{\boldsymbol{\mathbf{C}}},{\boldsymbol{\mathbf{d}}},{\boldsymbol{\mathbf{R}}}\\}$,
where ${\boldsymbol{\mathbf{C}}},{\boldsymbol{\mathbf{d}}}$ are the weights
and bias of the linear transformation on ${\boldsymbol{\mathbf{x}}}_{t}$ and
${\boldsymbol{\mathbf{R}}}$ is the observation covariance.

###  3.1 Joint Distribution

As shown by Frigola (2015), to sample ${\boldsymbol{\mathbf{f}}}_{t}$, instead
of conditioning on an infinite-dimensional function, we can condition only on
the transitions seen up to (but not including) time $t$, i.e.,
$\\{({\boldsymbol{\mathbf{x}}}_{i-1},{\boldsymbol{\mathbf{f}}}_{i})\\}_{i=1}^{t-1}$.
We can then write the joint distribution over latent variables and
observations as:

$p({\boldsymbol{\mathbf{y}}}_{1:T},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{f}}}_{1:T})=\\\
p({\boldsymbol{\mathbf{x}}}_{0})\prod_{t=1}^{T}p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1})p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t}),$
(4)

where
$p({\boldsymbol{\mathbf{x}}}_{0}),p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t}),p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})$
are defined as above and, for the edge case of $t=1$, we have:
$p({\boldsymbol{\mathbf{f}}}_{1}\,|\,{\boldsymbol{\mathbf{x}}}_{0})=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{1};m_{f}({\boldsymbol{\mathbf{x}}}_{0}),\kappa_{f}({\boldsymbol{\mathbf{x}}}_{0},{\boldsymbol{\mathbf{x}}}_{0};{\boldsymbol{\mathbf{\theta}}}))$.
Furthermore, we recognize each conditional distribution
$p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1})$
in Eq. 4 as the gp prediction at a single point
${\boldsymbol{\mathbf{x}}}_{t-1}$ using noiseless outputs
${\boldsymbol{\mathbf{f}}}_{1:t-1}$ and inputs
${\boldsymbol{\mathbf{x}}}_{0:t-2}$,

$\displaystyle
p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1})$
$\displaystyle=p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-2})$
(5)
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{f},{\boldsymbol{\mathbf{\Sigma}}}_{f}),$
(6)

with conditional mean and covariance given by

$\displaystyle{\boldsymbol{\mathbf{\mu}}}_{f}$
$\displaystyle={\boldsymbol{\mathbf{m}}}_{t-1}+{\boldsymbol{\mathbf{K}}}_{t-1,0:t-2}{\boldsymbol{\mathbf{K}}}_{0:t-2}^{-1}({\boldsymbol{\mathbf{f}}}_{1:t-1}-{\boldsymbol{\mathbf{m}}}_{0:T-2})$
(7) $\displaystyle{\boldsymbol{\mathbf{\Sigma}}}_{f}$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{t-1}-{\boldsymbol{\mathbf{K}}}_{t-1,0:t-2}{\boldsymbol{\mathbf{K}}}_{0:t-2}^{-1}{\boldsymbol{\mathbf{K}}}_{0:t-2,t-1}),$
(8)

where the subscript notation indicates the mean vectors and covariance
matrices obtained from evaluating the mean function and covariance function,
respectively, at the corresponding ranges,
${\boldsymbol{\mathbf{m}}}_{t_{1}:t_{2}}:=m_{f}({\boldsymbol{\mathbf{x}}}_{t_{1}:t_{2}})$,
${\boldsymbol{\mathbf{m}}}_{t}:={\boldsymbol{\mathbf{m}}}_{t:t}$,
${\boldsymbol{\mathbf{K}}}_{t_{1}:t_{2},t_{3}:t_{4}}:=\kappa_{f}({\boldsymbol{\mathbf{x}}}_{t_{1}:t_{2}},{\boldsymbol{\mathbf{x}}}_{t_{3}:t_{4}};{\boldsymbol{\mathbf{\theta}}})$,
${\boldsymbol{\mathbf{K}}}_{t_{1}:t_{2}}:={\boldsymbol{\mathbf{K}}}_{t_{1}:t_{2},t_{1}:t_{2}}$
and ${\boldsymbol{\mathbf{K}}}_{t}:={\boldsymbol{\mathbf{K}}}_{t:t}$.

###  3.2 Multidimensional Latent States & Control Inputs

In the case of multidimensional latent states, i.e., $d_{x}>1$, we assume
independent gps on each dimension, each with its own mean function and
covariance function. Because of this independence assumption, each gp only has
to condition on its own function evaluations. Therefore, for simplicity in the
notation, we do not index the means and covariances with respect to their
dimension $d$ and consider the underlying gps, their means, and covariance
functions as multi-dimensional. Furthermore, as we shall see in § 7, our
experiments consider additional control inputs
(${\boldsymbol{\mathbf{a}}}_{t}$) that affect the transitions in a Markovian
way. This is easy to incorporate in our framework by augmenting our input
space and, therefore, indexing the gps in the higher-dimensional space given
by the concatenation
$[{\boldsymbol{\mathbf{x}}}^{\top},{\boldsymbol{\mathbf{a}}}^{\top}]^{\top}$.
However, since the control inputs are fixed and deterministic, we do not need
to include them as part of out inference method and do not make them explicit
in our subsequent mathematical development.

###  3.3 Sparse gpssm

In order to make the gpssm inherently scalable, similarly to the standard
supervised regression setting with gps described in § 2, we augment the full
gpssm model with $M$ inducing variables and corresponding inducing inputs
${\boldsymbol{\mathbf{u}}}:=\\{u_{i}\\}_{i=1}^{M}$ and
${\boldsymbol{\mathbf{Z}}}:=\\{{\boldsymbol{\mathbf{z}}}_{i}\\}_{i=1}^{M}$.
This gives rise to the sparse gpssm model:

$p({\boldsymbol{\mathbf{y}}}_{1:T},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{f}}}_{1:T},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})=p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})p({\boldsymbol{\mathbf{x}}}_{0})\\\
\prod_{t=1}^{T}p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t}),$
(9)

where the prior over the inducing variables is determined by the gp prior,
i.e.,
$p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})=\mathcal{N}({\boldsymbol{\mathbf{m}}}_{Z},{\boldsymbol{\mathbf{K}}}_{Z})$
with ${\boldsymbol{\mathbf{m}}}_{Z}:=m_{f}({\boldsymbol{\mathbf{Z}}})$ and
${\boldsymbol{\mathbf{K}}}_{Z}:=\kappa_{f}({\boldsymbol{\mathbf{Z}}},{\boldsymbol{\mathbf{Z}}};{\boldsymbol{\mathbf{\theta}}})$.
Besides having a prior over the inducing variables
${\boldsymbol{\mathbf{u}}}$, the main difference with our previous full-model
formulation of § 3.1 is that the conditional distributions over
${\boldsymbol{\mathbf{f}}}_{t}$ have been augmented with inducing variables
${\boldsymbol{\mathbf{u}}}$ and corresponding inducing inputs
${\boldsymbol{\mathbf{Z}}}$. Each conditional,
$p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})=p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{0:t-2},{\boldsymbol{\mathbf{Z}}})$,
is the predictive (gp regression) distribution of
${\boldsymbol{\mathbf{f}}}_{t}$ at test input
${\boldsymbol{\mathbf{x}}}_{t-1}$ when observing
$({\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{u}}})$ at their
respective locations
$({\boldsymbol{\mathbf{x}}}_{0:t-2},{\boldsymbol{\mathbf{Z}}})$.

##  4 Free-Form Variational Inference

In this section, we develop a posterior estimation method using variational
inference (vi). We show that under a specific form of the approximate
posterior, we can estimate the joint posterior over
$\\{{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}\\}$ in free-
form, i.e., optimally, without making any assumptions such as independence
between state trajectories and inducing variables, typical of mean-field
approaches (Frigola et al., 2014; Eleftheriadis et al., 2017), or constraining
the form of the posterior to sub-optimal parametric forms (Ialongo et al.,
2019; Doerr et al., 2018).

###  4.1 Variational Family and Evidence Lower Bound

Variational inference is underpinned by the maximization of the evidence lower
bound (elbo), which is equivalent to minimizing the Kullback-Leibler (kl)
divergence between the approximate posterior and the true posterior. This
objective is given by

$\mathcal{L}_{\text{ELBO}}(q):=\mathbb{E}_{q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}_{1:T})}\left[\log{p({\boldsymbol{\mathbf{y}}}_{1:T},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{f}}}_{1:T},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})}\right.\\\
\left.-\log{q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}_{1:T})}\right],$
(10)

where
$p({\boldsymbol{\mathbf{y}}}_{1:T},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{f}}}_{1:T},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})$
is the joint distribution in Eq. 9 and
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}_{1:T})$
is our proposed approximate joint posterior, which we define as

$\displaystyle
q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$
$\displaystyle:=\prod_{t=1}^{T}p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}}),$
(11) $\displaystyle
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{f}}}_{1:T})$
$\displaystyle:=q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}).$
(12)

It is easy to show that the elbo, as defined in Eq. 10, is a lower bound on
the log marginal likelihood, i.e., $\mathcal{L}_{\text{ELBO}}(q)\leq\log
p({\boldsymbol{\mathbf{y}}}\,|\,{\boldsymbol{\mathbf{Z}}})$. Furthermore, the
definition of our variational distribution above, where we used the
conditional prior in Eq. 11, will yield an inference algorithm that alleviates
the cubic time complexity on the number of observations and, therefore, scales
up to large datasets. Crucially, our proposed joint variational distribution
in Eq. 12, will allow us to derive an _optimal_ variational distribution
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$ in free-form,
without imposing any parametric constraints over it and, instead, represent it
via samples. We will describe this in the next section.

###  4.2 Evidence Lower Bound Maximization

Our first step is to expand the expression for the elbo in Eq. 10 using our
joint model distribution in Eq. 9 and our proposed variational distribution in
Eqs. 11 and 12. We first note that our definition of the variational
distribution in Eq. 12 uses the same conditional prior as in the joint
distribution in Eq. 9. Therefore, this term cancels out, avoiding the
computation of operations on fully-coupled high-dimensional distributions over
latent functions. Thus, we have that:

$\displaystyle\mathcal{L}_{\text{ELBO}}(q)=\int
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\Bigg{\\{}-\log
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$
$\displaystyle+\log\Big{[}p({\boldsymbol{\mathbf{x}}}_{0})p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})\prod_{t=1}^{T}p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\Big{]}$
$\displaystyle+\mathbb{E}_{q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})}\Big{[}\log\prod_{t=1}^{T}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})\Big{]}\Bigg{\\}}\mathrm{d}{\boldsymbol{\mathbf{x}}}_{0:T}\mathrm{d}{\boldsymbol{\mathbf{u}}}.$
(13)

####  4.2.1 elbo Maximization

Next we aim to maximize the elbo functional above with respect to
$q({\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}})$ subject to
the constraint $\int
q({\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}})\
\mathrm{d}{\boldsymbol{\mathbf{x}}}_{0:T-1}\mathrm{d}{\boldsymbol{\mathbf{u}}}=1$.
We can do this by solving the corresponding Euler-Lagrange equation:

$\frac{\partial}{\partial
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})}\Bigg{\\{}-q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\log
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\\\
+q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\log\Big{[}p({\boldsymbol{\mathbf{x}}}_{0})p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})\prod_{t=1}^{T}p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\Big{]}\\\
+q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\mathbb{E}_{q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})}\Big{[}\log\prod_{t=1}^{T}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})\Big{]}\Bigg{\\}}=0.$
(14)

By doing the corresponding derivatives we obtain

$-\log
q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})-1+\log\Big{[}p({\boldsymbol{\mathbf{x}}}_{0})p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})\prod_{t=1}^{T}p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\Big{]}\\\
+\mathbb{E}_{q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})}\Big{[}\log\prod_{t=1}^{T}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})\Big{]}=0.$
(15)

Now we note that the expectation above can be solved in closed form:

$\displaystyle\mathbb{E}_{q({\boldsymbol{\mathbf{f}}}_{1:T}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})}\Big{[}\log\prod_{t=1}^{T}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})\Big{]}$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\mathbb{E}_{p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})}\log
p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\left[\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{u}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\right],$
(16)

where
$p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})$
is the gp predictive distribution over the function values
${\boldsymbol{\mathbf{f}}}_{t}$ at locations ${\boldsymbol{\mathbf{x}}}_{t-1}$
given the inducing variables ${\boldsymbol{\mathbf{u}}}$ at inducing inputs
${\boldsymbol{\mathbf{Z}}}$, i.e.,

$\displaystyle
p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{t};{\boldsymbol{\mathbf{u}}}_{x_{t}},{\boldsymbol{\mathbf{B}}}_{t-1}),$
(17) $\displaystyle{\boldsymbol{\mathbf{u}}}_{x_{t}}$
$\displaystyle:={\boldsymbol{\mathbf{m}}}_{t-1}+{\boldsymbol{\mathbf{A}}}_{t-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z}).$
(18)

Here we have defined

$\displaystyle{\boldsymbol{\mathbf{A}}}_{t-1}$
$\displaystyle:={\boldsymbol{\mathbf{K}}}_{t-1,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1},$
(19) $\displaystyle{\boldsymbol{\mathbf{B}}}_{t-1}$
$\displaystyle:={\boldsymbol{\mathbf{K}}}_{t-1}-{\boldsymbol{\mathbf{K}}}_{t-1,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,t-1},$
(20)

and the cross-covariance term
${\boldsymbol{\mathbf{K}}}_{t-1,Z}:=\kappa_{f}({\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{Z}}};{\boldsymbol{\mathbf{\theta}}})$
and similarly for ${\boldsymbol{\mathbf{K}}}_{Z,t-1}$. Finally,
$\operatorname{Tr}(\cdot)$ is the trace operator and, as defined at the
beginning of § 3, ${\boldsymbol{\mathbf{Q}}}$ is the transition noise
covariance.

####  4.2.2 Optimal Variational Posterior

With this, we obtain the form of the optimal variational distribution
$q^{*}({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{0:T})$ up to a
normalizing constant $Z_{q}$ as:

$\log q^{*}({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{0:T})=\log
p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})+\log
p({\boldsymbol{\mathbf{x}}}_{0})\\\ +\sum_{t=1}^{T}\Big{[}\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})\\\
-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\Big{]}+\log
Z_{q}.$ (21)

Here we note that the function values ${\boldsymbol{\mathbf{f}}}_{1:T}$ have,
effectively, been marginalized variationally. The optimal joint posterior over
inducing variables and state trajectories depends on the prior over the
inducing variables
$p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})$ stemming from the
gp functional prior, the prior over the initial state
$p({\boldsymbol{\mathbf{x}}}_{0})$ and the conditional likelihood terms
$p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})$. It also
depends on the resulting transitions mapping ${\boldsymbol{\mathbf{x}}}_{t-1}$
to ${\boldsymbol{\mathbf{x}}}_{t}$ via the densities
$\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})$,
where ${\boldsymbol{\mathbf{\mu}}}_{x_{t}}$ depends on
${\boldsymbol{\mathbf{x}}}_{t-1}$ in a nonlinear way, as specified by Eq. 19.
The final trace term,
$\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})$,
can be seen as a regularization term acting on state transitions, encouraging
higher transition variances and, therefore, helping prevent overfitting.

####  4.2.3 Alternative Perspective

An alternative way to obtain the optimal joint posterior over state
trajectories and inducing variable is by bounding the true log joint marginal
$\log
p({\boldsymbol{\mathbf{y}}}_{1:T},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})$
using Jensen’s inequality. We give details of such an approach in Appendix E.
This has been used by previous work in standard regression settings (see,
e.g., Rossi et al., 2021, and references therein). However, our setting
considers the more complex case of gpssms. Additionally, our development in
this section shows much more clearly the optimal nature of the variational
posterior, as we have obtained it via calculus of variations.

###  4.3 Posterior Sampling

Having the form of the optimal posterior in Eq. 21, we can then set the latent
variables
${\boldsymbol{\mathbf{\Psi}}}:=\\{{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{0:T}\\}$
and have $\tilde{q}({\boldsymbol{\mathbf{\Psi}}})\propto
q({\boldsymbol{\mathbf{\Psi}}})\approx
p({\boldsymbol{\mathbf{\Psi}}}\,|\,{\boldsymbol{\mathbf{y}}}_{1:T})$. Thus, we
can draw samples from our approximate posterior using stochastic gradient
Hamiltonian Monte Carlo (sghmc; Chen et al., 2014; Havasi et al., 2018) and
the energy function $U({\boldsymbol{\mathbf{\Psi}}})=-\log
p({\boldsymbol{\mathbf{\Psi}}},{\boldsymbol{\mathbf{y}}}_{1:T})=-\log
p({\boldsymbol{\mathbf{\Psi}}}\,|\,{\boldsymbol{\mathbf{y}}}_{1:T})+\log
Z_{q}\approx-\log\tilde{q}({\boldsymbol{\mathbf{\Psi}}})+\log Z_{q}$. Using
this procedure, samples from the target distribution can be obtained even with
noisy gradients (e.g., with mini-batches) without requiring the evaluation of
Metropolis ratios. Importantly, other variables such as gp hyper-parameters
${\boldsymbol{\mathbf{\theta}}}$ and inducing locations
${\boldsymbol{\mathbf{Z}}}$ can be easily included in
${\boldsymbol{\mathbf{\Psi}}}$ using suitable priors and incorporating them in
our objective in Eq. 21.

####  4.3.1 Computational Cost & Prior Whitening

An interesting aspect of gpssm models is that, despite their apparent
Markovian nature, sampling at time $T$ requires conditioning on all the
previous $T-1$ points. This is due to the non-parametric coupled gp prior over
the transition function, making inference in the full model
$\mathcal{O}(T^{3})$ in time. Sparse variational inference approaches, such as
those based on inducing variable approximations, still require expectations
over entire trajectories and, unlike standard supervised i.i.d settings, their
time complexity is inherently dependent on $T$. Evaluation of the stochastic
gradient Hamiltonian Monte Carlo (sghmc) objective in Eq. 21 for sampling in
our free-form variational dynamics (ffvd) algorithm is $\mathcal{O}(M^{2}T)$.
this is the same cost as that attained for elbo evaluation in variationally
coupled dynamic trajectories (vcdt), which like our ffvd, models dependencies
between state trajectories and inducing variables.

As in previous work (see, e.g., Hensman et al., 2015), we have observed that
whitening the prior over the inducing variables improves the performance of
our algorithm. See details of our whitening procedure and the resulting
unnormalized log posterior in Appendices D and E.1.

###  4.4 Smoothing and Predictive Distributions

We are interested in estimating the smoothing distribution
$p({\boldsymbol{\mathbf{x}}}_{0:T}\,|\,{\boldsymbol{\mathbf{y}}}_{1:T})$ and
the predictive distribution
$p({\boldsymbol{\mathbf{y}}}_{T+1:T^{\prime}}\,|\,{\boldsymbol{\mathbf{y}}}_{1:T})$
for $T^{\prime}\geq T$. At the end of our sghmc procedure, we have $S$ samples
from our joint approximate posterior
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})\approx
p({\boldsymbol{\mathbf{x}}}_{0:T}\,|\,{\boldsymbol{\mathbf{y}}}_{1:T})$, i.e.,
$\\{{\boldsymbol{\mathbf{x}}}_{0:T}^{(s)},{\boldsymbol{\mathbf{u}}}^{(s)}\\}_{s=1}^{S}$
and, therefore, the smoothing distribution (i.e., the marginal posterior over
${\boldsymbol{\mathbf{x}}}_{0:T}$) is readily available through this Monte
Carlo approximation.

We can also make one-step-ahead predictions using our posterior samples. In
particular, using Eqs. 17 and 2 we have that:

$p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})=\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{m}}}_{t-1}+{\boldsymbol{\mathbf{A}}}_{t-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z}),\\\
{\boldsymbol{\mathbf{B}}}_{t-1}+{\boldsymbol{\mathbf{Q}}}),$ (22)

$\forall t>T$. Thus, replacing the values of
${\boldsymbol{\mathbf{A}}}_{t-1},{\boldsymbol{\mathbf{B}}}_{t-1}$ using Eqs.
19 and 20 we can make predictions for the next state using samples as

${\boldsymbol{\mathbf{x}}}_{t}^{(s)}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1}^{(s)},{\boldsymbol{\mathbf{u}}}^{(s)}\sim\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{t}^{(s)},{\boldsymbol{\mathbf{\Sigma}}}_{t}^{(s)})\quad\text{with}\\\
{\boldsymbol{\mathbf{\mu}}}_{t}^{(s)}={\boldsymbol{\mathbf{m}}}_{t-1}^{(s)}+{\boldsymbol{\mathbf{K}}}_{t-1,Z}^{(s)}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}^{(s)}-{\boldsymbol{\mathbf{m}}}_{Z}),\\\
{\boldsymbol{\mathbf{\Sigma}}}_{t}^{(s)}={\boldsymbol{\mathbf{K}}}_{t-1}^{(s)}-{\boldsymbol{\mathbf{K}}}_{t-1,Z}^{(s)}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,t-1}^{(s)}+{\boldsymbol{\mathbf{Q}}}.$
(23)

For a general likelihood model, we can sample the (noisy) targets using
${\boldsymbol{\mathbf{y}}}_{t}^{(s)}\,|\,{\boldsymbol{\mathbf{x}}}_{t}^{(s)}\sim
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t}^{(s)},{\boldsymbol{\mathbf{\phi}}}).$
In the case of a Gaussian conditional likelihood, as described in § 3, we can
see that given samples from the latent state, the predictive distribution is a
Gaussian

${\boldsymbol{\mathbf{y}}}_{t}^{(s)}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1}^{(s)},{\boldsymbol{\mathbf{u}}}^{(s)}\sim\mathcal{N}({\boldsymbol{\mathbf{y}}}_{t};{\boldsymbol{\mathbf{C}}}{\boldsymbol{\mathbf{\mu}}}_{t}^{(s)}+{\boldsymbol{\mathbf{d}}},\\\
{\boldsymbol{\mathbf{C}}}{\boldsymbol{\mathbf{\Sigma}}}_{t}^{(s)}{\boldsymbol{\mathbf{C}}}^{\top}+{\boldsymbol{\mathbf{R}}}).$
(24)

##  5 Collapsing Inducing Variables

So far we have described our method to obtain samples from the optimal
variational posterior over the joint distribution of state trajectories and
inducing variables using Eq. 21 and sghmc. In this section we show that we
can, in fact, integrate out the inducing variables ${\boldsymbol{\mathbf{u}}}$
from our joint variational distribution and obtain the optimal marginal
distribution for latent states ${\boldsymbol{\mathbf{x}}}_{0:T}$. We start by
retaking Eq. 21 and isolating the terms that depend on the inducing variables,

$q^{*}({{\boldsymbol{\mathbf{x}}}}_{0:T})=\int_{{{\boldsymbol{\mathbf{u}}}}}q^{*}({{\boldsymbol{\mathbf{u}}}},{{\boldsymbol{\mathbf{x}}}}_{0:T})d{{\boldsymbol{\mathbf{u}}}}\\\
=p({{\boldsymbol{\mathbf{x}}}}_{0})\prod_{t=1}^{T}\left[p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{{\boldsymbol{\mathbf{x}}}}_{t})\exp(-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1}))\right]\\\
\int_{{{\boldsymbol{\mathbf{u}}}}}p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})\prod_{t=1}^{T}\mathcal{N}({{\boldsymbol{\mathbf{x}}}}_{t=1}^{T};{\boldsymbol{\mathbf{\mu}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})\mathrm{d}{\boldsymbol{\mathbf{u}}},$
(25)

where we note that ${\boldsymbol{\mathbf{\mu}}}_{x_{t}}$, as defined as in Eq.
18, depends on the inducing variables ${\boldsymbol{\mathbf{u}}}$. With this,
we can complete the square and identify the terms in the integral as products
of Gaussian distributions, whose normalization constant is a Gaussian.
Therefore, we have

$q^{*}({{\boldsymbol{\mathbf{x}}}}_{0:T})=p({{\boldsymbol{\mathbf{x}}}}_{0})\prod_{t}\Big{[}p({\boldsymbol{\mathbf{y}}}_{t}|{{\boldsymbol{\mathbf{x}}}}_{t})\exp(-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1}))\\\
\mathcal{N}({{\boldsymbol{\mathbf{x}}}}_{t};{\boldsymbol{\mathbf{m}}}_{t-1},{\boldsymbol{\mathbf{Q}}})\Big{]}/\mathcal{N}(\tilde{{\boldsymbol{\mathbf{x}}}};\mathbf{0},\tilde{{\boldsymbol{\mathbf{\Sigma}}}}_{x}),$
(26)

where the Gaussian density in the denominator is determined by

$\displaystyle\tilde{{\boldsymbol{\mathbf{x}}}}$
$\displaystyle=\sum_{t=1}^{T}\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{\top}{\boldsymbol{\mathbf{Q}}}^{-1}({\boldsymbol{\mathbf{x}}}_{t}-{\boldsymbol{\mathbf{m}}}_{t-1}),$
(27) $\displaystyle\tilde{{\boldsymbol{\mathbf{\Sigma}}}}_{x}$
$\displaystyle={\boldsymbol{\mathbf{I}}}+\sum_{t=1}^{T}\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{\top}{\boldsymbol{\mathbf{Q}}}^{-1}\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1},$
(28) $\displaystyle\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{t-1,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{\top})^{-1},$
(29)

where ${\boldsymbol{\mathbf{L}}}_{Z}$ is the Cholesky decomposition of
${\boldsymbol{\mathbf{K}}}_{Z}$, i.e.,
${\boldsymbol{\mathbf{K}}}_{Z}={\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{L}}}_{Z}^{\top}$.
We note here that $\tilde{{\boldsymbol{\mathbf{x}}}}$ is actually a projection
of the cumulative uncorrelated inputs
${\boldsymbol{\mathbf{Q}}}^{-1}({\boldsymbol{\mathbf{x}}}_{t}-{\boldsymbol{\mathbf{m}}}_{t-1})$
via the projection matrix
${\boldsymbol{\mathbf{L}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,t-1}$ and,
therefore, $\tilde{{\boldsymbol{\mathbf{x}}}}\in\mathbb{R}^{M}$.

Furthermore, in this collapsed version, the terms unrelated to
${\boldsymbol{\mathbf{u}}}$ are kept the same as those in the original optimal
joint distribution in Eq. 21. The individual Gaussian transitions for the
latent states are now
$\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{m}}}_{t-1},{\boldsymbol{\mathbf{Q}}})$
since ${\boldsymbol{\mathbf{u}}}$ has been integrated out and we recall that
${\boldsymbol{\mathbf{m}}}_{t-1}=m_{f}({\boldsymbol{\mathbf{x}}}_{t-1})$. The
additional term, which is the inverse of a multivariate Gaussian distribution,
records the cumulative projected trajectory with the corresponding projected
variances. Maximizing $q^{*}({\boldsymbol{\mathbf{x}}}_{0:T})$ would tend to
minimize this density function, which pushes the latent states away from
$\mathbf{0}$ and also decreases the values of the ${\boldsymbol{\mathbf{Q}}}$
regulated variance.

###  5.1 Advantages of Collapsed Algorithm

Collapsing the inducing variables ${\boldsymbol{\mathbf{u}}}$ will generally
tend to improve the convergence of our algorithm, as we are required to sample
from a significantly lower number of latent variables. The computational cost
is similar to that of the uncollapsed algorithm since despite avoiding the
computation of ${\boldsymbol{\mathbf{K}}}_{Z}^{-1}$ in
$p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})$, we require a
similar term, ${\boldsymbol{\mathbf{L}}}_{Z}^{-1}$, in the evaluation of the
projected Gaussian distributions.

It is important to emphasize one particular difference in our approach with
respect to the closely related vcdt algorithm of Ialongo et al. (2019). As
described before, Ialongo et al. (2019) also propose a coupled joint posterior
between state trajectories and inducing variables. Their factorization is
$q_{\lx@glossaries@gls@link{acronym}{VCDT}{\leavevmode\textsc{vcdt}}}=q({\boldsymbol{\mathbf{u}}})q({\boldsymbol{\mathbf{x}}}_{0:T}\,|\,{\boldsymbol{\mathbf{u}}})$,
and they impose additional Gaussian constraints on these densities. Our
implicit assumed factorization is
$q_{\lx@glossaries@gls@link{acronym}{FFVD}{\leavevmode\textsc{ffvd}}}=q({\boldsymbol{\mathbf{x}}}_{0:T})q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$,
which allows us to obtain the optimal variational distribution without
imposing any additional parametric constraints, integrate out the inducing
variables analytically and get a Monte Carlo approximation to the optimal
marginal $q({\boldsymbol{\mathbf{x}}}_{0:T})$ via samples. We also provide an
expression for the conditional
$q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$ in closed-
form. Details can be found in Appendix G.

###  5.2 Particle Markov chain Monte Carlo (pmcmc)

It is clear that the latent states are constructed in a Markovian manner when
the transition function $f(\cdot)$ is given, as the value of the current
latent state is dependent on the previous latent state’s value. Therefore, we
can use pmcmc methods (Andrieu et al., 2010) to infer the posterior
distribution of the Markov structured latent states
${\boldsymbol{\mathbf{x}}}_{0:T}$. This Bayesian treatment might improve
performance over sghmc, as it incorporates the sequential nature of the
problem into the sampling algorithm. The advantages of pmcmc over standard
sequential Monte Carlo approaches have been documented previously, see. e.g.,
Andrieu et al. (2010). Here we note that Frigola et al. (2013) also proposed a
pmcmc treatment for ${\boldsymbol{\mathbf{x}}}_{0:T}$. However, their
algorithm is very different from ours as it is based on the fully-independent
conditional approximation (see, e.g., Quinonero-Candela & Rasmussen, 2005).
More details can be found in Appendix K.

##  6 Related Work

We have already described the main differences between our method and closely-
related approaches throughout the paper, e.g., in §§ 1, 4 and 5. We refer the
reader to Appendix K for more details. Other works have considered the gpssm
in partially observable unstable settings (Curi et al., 2020) or combined
variational inference with the Laplace approximation (Lindinger et al., 2022).
With regards to scalable gps, we note they have been the subject of much
research effort in machine learning, with extensions to more general
frameworks such as compositional models (see, e.g., Wilson & Nickisch, 2015;
Salimbeni & Deisenroth, 2017; Yu et al., 2019; Cutajar et al., 2017; Havasi et
al., 2018; Rossi et al., 2021). Our approach can be seen as a generalization
of the fully-Bayesian supervised learning method proposed by Rossi et al.
(2021) to state-space models, where we have included the non-trivial component
of gp-transition dynamics and have proposed a collapsed optimal variational
distribution for state trajectories.

Other models for gps on sequential data have been proposed, see, for example,
Frigola (2015) for an excellent overview. Of interest here is the state-space
model view of gps that for time-series data with $d_{x}=1$ and Markovian
covariance functions can provide exact inference in linear time
$\mathcal{O}(T)$ (Solin, 2016). This has been extended to non-Gaussian
likelihood models and made more efficient using several computational
primitives (Nickisch et al., 2018). More recently, there has been some work on
using the signature kernel (Toth & Oberhauser, 2020; Salvi et al., 2021)
within gp models. In particular, Lemercier et al. (2021) generalize
variational orthogonal features (Burt et al., 2020; Hensman et al., 2017) to
the sequential case, constructing inducing variables associated with the
signature kernel that yield a variational inference algorithm that does not
require any matrix inversion.

Figure 1: Results on synthetic data. _Left:_ Observations shown as green dots,
Ground truth as a solid black line, and ffvd’s mean fitting as a dashed red
line with one standard deviation error bars. _Right:_ Histograms of p-values
for the hypothesis test that each marginal posterior over states $\\{x_{t}\\}$
(top) and inducing variables $\\{u_{i}\\}$ (bottom) is generated from a Normal
distribution. Figure 2: Training ($t<=150$) and test performance ($t>=150$) on
the Furnace dataset. The black solid line is the underlying ground truth
signal and the dashed red line is ffvd’s mean prediction. The blue solid line
indicates the training/test split. An underlying $d_{x}=4$-dimensional latent
stated was used. Figure 3: Traceplot of the training log-likelihood when using
sghmc (ffvd-m) and when collapsing the inducing variables ffvd-c-m. Collapsing
generally improved convergence.

##  7 Experiments

We evaluate our ffvd method on synthetic data and on six real-world system
identification benchmarks (Ialongo et al., 2019; Doerr et al., 2018),
comparing it with vgpssm (Frigola et al., 2014), prssm (Doerr et al., 2018),
vcdt (Ialongo et al., 2019), and use a LSTM network (Hochreiter & Schmidhuber,
1997) as a baseline non-gp based model. All experiment details can be found in
Appendix I.

Table 2: Mean rmse values $\pm$ one standard deviations on the real-world system identification benchmarks. Our method, ffvd, when using sghmc (ffvd-m); the collapsed version (ffvd-c-m); and when using pmcmc (ffvd-p).  Methods | Actuator | Ballbeam | Drive | Dryer | Flutter | Furnace
---|---|---|---|---|---|---
lstm | $0.586\pm 0.411$ | $0.027\pm 0.023$ | $0.537\pm 0.108$ | $0.115\pm 0.029$ | $0.912\pm 0.562$ | $1.261\pm 0.610$
vgpssm | $0.580\pm 0.274$ | $0.073\pm 0.011$ | $0.722\pm 0.087$ | $0.241\pm 0.023$ | $1.482\pm 0.218$ | $1.115\pm 0.358$
prssm | $0.497\pm 0.381$ | $0.059\pm 0.013$ | $0.813\pm 0.101$ | $\mathbf{0.017}\pm 0.042$ | $1.371\pm 0.156$ | $1.243\pm 0.407$
vcdt | $\mathbf{0.239}\pm 0.040$ | $0.011\pm 0.002$ | $0.585\pm 0.017$ | $0.142\pm 0.003$ | $1.782\pm 0.324$ | $1.166\pm 0.011$
ffvd-m | $0.358\pm 0.242$ | $0.019\pm 0.018$ | $0.673\pm 0.207$ | $0.205\pm 0.313$ | $\mathbf{0.280}\pm 0.193$ | $0.571\pm 0.185$
ffvd-c-m | $\underline{0.259}\pm 0.209$ | $\mathbf{0.009}\pm 0.011$ | $0.775\pm 1.615$ | $\underline{0.065}\pm 0.112$ | $0.663\pm 0.189$ | $\mathbf{0.548}\pm 0.051$
ffvd-p | ${0.388}\pm 0.087$ | ${0.199}\pm 0.045$ | $\mathbf{0.342}\pm 0.057$ | ${0.317}\pm 0.050$ | $0.562\pm 0.088$ | ${0.669}\pm 0.174$

###  7.1 Synthetic Data

We generate data from a sparse gpssm with a squared exponential covariance
function. Our goal here is to investigate the properties of our ffvd
algorithm. Therefore, we fix the value of all parameters to their ground-truth
values except for the latent states ${\boldsymbol{\mathbf{x}}}_{0:T}$ and the
inducing variables ${\boldsymbol{\mathbf{u}}}$. The left panel of Fig. 1
illustrates that ffvd effectively learns the intricate transition function
with two modes within regions containing latent states. As the number of
latent states increases, ffvd more accurately approximates the true function.
A lack of fit in regions without latent states is to be expected. Having a
good fit, we are now interested in knowing whether the true posterior (as
estimated by our algorithm) is close to a Gaussian distribution. The right
panel in Fig. 1 indicates that more than $10\%$ of the latent states
${\boldsymbol{\mathbf{x}}}_{0:T}$ and more than $50\%$ of the inducing
variables ${\boldsymbol{\mathbf{u}}}$ do not provide enough support for the
hypothesis that their marginal distributions follow a Normal distribution (see
details in Appendix I). This brings into question the strong parametric
assumptions over the variational posterior made by previous work (e.g.,
Ialongo et al., 2019; Doerr et al., 2018).

###  7.2 Real-World Data

Here we evaluate the different methods using six system identification
benchmarks with a latent state dimension $d_{x}=4$, as used by Ialongo et al.
(2019); Doerr et al. (2018). See more details in Appendix I. The prediction
performance is shown in Table 2, where we see that ffvd attains the lowest
rmse values in four out of the six benchmarks (in bold). Furthermore, for
_Actuator_ and _Dryer_ , ffvd-c-m ranks second among all the algorithms
(underlined). The performance of vgpssm and prssm is usually worse than
others, which is likely due to their strong mean-field and parametric
assumptions, respectively. lstm obtains good performance on three datasets,
although its deterministic neural network structure is different from our
random function setting. The performance of vcdt is the closest to our ffvd
methods, as it models a coupled
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$.

In addition to this quantitative evaluation, we can see a qualitative
illustration of using our algorithm for predicting the training and test
(future) observations in Fig. 2. This is an example of good generalization,
although it is (of course) not consistent across all problems, given the
limited training data. Finally, we analyze the convergence of our algorithm in
Fig. 3, where we see that ffvd-c-m uses less iterations ($\sim 8\,000$
iterations) than ffvd-m ($\sim 40\,000$ iterations) to achieve similar
performance. This confirms the benefits of collapsing the inducing variables
and sampling only on the lower-dimensional space of state trajectories. We can
see these analyses for all benchmarks in §§ J.1 and J.2.

##  8 Conclusion

We have presented ffvd, a new algorithm for inference in gpssms based on
variational inference under the inducing-variable formalism. Unlike previous
approaches, ffvd does not make any independence or parametric assumptions on
the joint variational posterior over state trajectories and inducing variables
$q({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})$ and, instead,
represents the posterior via samples form the optimal variational
distribution. Furthermore, we have shown that our formulation allows us to
collapse the inducing variables and carry out inference in the lower
dimensional space of state trajectories. Given samples from trajectories, we
can make predictions using the optimal conditional
$q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$, for which
we have derived a closed-form expression. Future work will investigate
scalable multi-output gps in order to deal with very high-dimensional state
representations in a flexible and efficient way.

## References

  * Andrieu et al. (2010) Andrieu, C., Doucet, A., and Holenstein, R.  Particle Markov chain Monte Carlo methods.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 72(3):269–342, 2010.
  * Bauer et al. (2016) Bauer, M., van der Wilk, M., and Rasmussen, C. E.  Understanding probabilistic sparse gaussian process approximations.  _Advances in neural information processing systems_ , 29, 2016.
  * Brown et al. (1998) Brown, E. N., Frank, L. M., Tang, D., Quirk, M. C., and Wilson, M. A.  A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells.  _Journal of Neuroscience_ , 18(18):7411–7425, 1998.
  * Burt et al. (2020) Burt, D. R., Rasmussen, C. E., and van der Wilk, M.  Variational orthogonal features.  _arXiv preprint arXiv:2006.13170_ , 2020.
  * Chen et al. (2014) Chen, T., Fox, E., and Guestrin, C.  Stochastic gradient Hamiltonian Monte Carlo.  In _ICML_ , pp. 1683–1691. PMLR, 2014.
  * Curi et al. (2020) Curi, S., Melchior, S., Berkenkamp, F., and Krause, A.  Structured variational inference in partially observable unstable Gaussian process state space models.  In Bayen, A. M., Jadbabaie, A., Pappas, G., Parrilo, P. A., Recht, B., Tomlin, C., and Zeilinger, M. (eds.), _Proceedings of the 2nd Conference on Learning for Dynamics and Control_ , volume 120 of _Proceedings of Machine Learning Research_ , pp. 147–157. PMLR, 10–11 Jun 2020.
  * Cutajar et al. (2017) Cutajar, K., Bonilla, E. V., Michiardi, P., and Filippone, M.  Random feature expansions for deep Gaussian processes.  In Precup, D. and Teh, Y. W. (eds.), _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pp. 884–893. PMLR, 06–11 Aug 2017.
  * D’agostino & Pearson (1973) D’agostino, R. and Pearson, E. S.  Tests for departure from normality.  _Biometrika_ , 60(3):613–622, 1973.
  * Doerr et al. (2018) Doerr, A., Daniel, C., Schiegg, M., Duy, N.-T., Schaal, S., Toussaint, M., and Sebastian, T.  Probabilistic recurrent state-space models.  In _ICML_ , pp. 1280–1289, 2018.
  * Eleftheriadis et al. (2017) Eleftheriadis, S., Nicholson, T., Deisenroth, M., and Hensman, J.  Identification of Gaussian process state space models.  _Advances in neural information processing systems_ , 30, 2017.
  * Frigola (2015) Frigola, R.  _Bayesian time series learning with Gaussian processes_.  PhD thesis, University of Cambridge, 2015.
  * Frigola et al. (2013) Frigola, R., Lindsten, F., Schön, T. B., and Rasmussen, C. E.  Bayesian inference and learning in Gaussian process state-space models with particle MCMC.  In _NIPS_ , pp. 3156–3164, 2013.
  * Frigola et al. (2014) Frigola, R., Chen, Y., and Rasmussen, C. E.  Variational Gaussian process state-space models.  In _NIPS_ , pp. 3680–3688, 2014.
  * Havasi et al. (2018) Havasi, M., Hernández-Lobato, J. M., and Murillo-Fuentes, J. J.  Inference in deep Gaussian processes using stochastic gradient Hamiltonian Monte Carlo.  _NeurIPS_ , 31, 2018.
  * Hensman et al. (2013) Hensman, J., Fusi, N., and Lawrence, N. D.  Gaussian processes for big data.  In _Uncertainty in Artificial Intelligence_. AUAI Press, 2013.
  * Hensman et al. (2015) Hensman, J., Matthews, A. G., Filippone, M., and Ghahramani, Z.  MCMC for variationally sparse Gaussian processes.  In _NIPS_ , pp. 1648–1656, 2015.
  * Hensman et al. (2017) Hensman, J., Durrande, N., Solin, A., et al.  Variational Fourier features for Gaussian processes.  _J. Mach. Learn. Res._ , 18(1):5537–5588, 2017\.
  * Hernandez et al. (1991) Hernandez, E., Martin, F., and Valero, F.  State-space modeling for atmospheric pollution.  _Journal of Applied Meteorology and Climatology_ , 30(6):793–811, 1991.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Ialongo et al. (2019) Ialongo, A. D., Van Der Wilk, M., Hensman, J., and Rasmussen, C. E.  Overcoming mean-field approximations in recurrent Gaussian process models.  In _ICML_ , pp. 2931–2940, 2019.
  * Kalman (1960) Kalman, R. E.  A New Approach to Linear Filtering and Prediction Problems.  _Journal of Basic Engineering_ , 82(1):35–45, 1960.  doi: 10.1115/1.3662552.
  * Katzfuss et al. (2016) Katzfuss, M., Stroud, J. R., and Wikle, C. K.  Understanding the Ensemble Kalman Filter.  _The American Statistician_ , 70(4):350–357, 2016\.  doi: 10.1080/00031305.2016.1141709.
  * Kitagawa (1987) Kitagawa, G.  Non-Gaussian state—space modeling of nonstationary time series.  _Journal of the American statistical association_ , 82(400):1032–1041, 1987.
  * Lemercier et al. (2021) Lemercier, M., Salvi, C., Cass, T., Bonilla, E. V., Damoulas, T., and Lyons, T. J.  SigGPDE: Scaling sparse Gaussian processes on sequential data.  In _International Conference on Machine Learning_ , pp. 6233–6242. PMLR, 2021.
  * Lindinger et al. (2022) Lindinger, J., Rakitsch, B., and Lippert, C.  Laplace approximated Gaussian process state-space models.  In _The 38th Conference on Uncertainty in Artificial Intelligence_ , 2022.
  * Marmin & Filippone (2022) Marmin, S. and Filippone, M.  Deep gaussian processes for calibration of computer models (with discussion).  _Bayesian Analysis_ , 17(4):1301–1350, 2022.
  * Murphy (2023) Murphy, K. P.  _Probabilistic Machine Learning: Advanced Topics_.  MIT Press, 2023.  URL http://probml.github.io/book2.
  * Nickisch et al. (2018) Nickisch, H., Solin, A., and Grigorevskiy, A.  State space Gaussian processes with non-Gaussian likelihood.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 3789–3798. PMLR, 10–15 Jul 2018.
  * Ogata et al. (2010) Ogata, K. et al.  _Modern control engineering_ , volume 5.  Prentice hall Upper Saddle River, NJ, 2010.
  * Quinonero-Candela & Rasmussen (2005) Quinonero-Candela, J. and Rasmussen, C. E.  A unifying view of sparse approximate Gaussian process regression.  _The Journal of Machine Learning Research_ , 6:1939–1959, 2005.
  * Rossi et al. (2021) Rossi, S., Heinonen, M., Bonilla, E., Shen, Z., and Filippone, M.  Sparse Gaussian processes revisited: Bayesian approaches to inducing-variable approximations.  In _AISTATS_ , pp. 1837–1845, 2021.
  * Salimbeni & Deisenroth (2017) Salimbeni, H. and Deisenroth, M.  Doubly stochastic variational inference for deep Gaussian processes.  In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.
  * Salvi et al. (2021) Salvi, C., Cass, T., Foster, J., Lyons, T., and Yang, W.  The signature kernel is the solution of a goursat PDE.  _SIAM Journal on Mathematics of Data Science_ , 3(3):873–899, 2021.
  * Sauer et al. (2022) Sauer, A., Cooper, A., and Gramacy, R. B.  Vecchia-approximated deep Gaussian processes for computer experiments.  _Journal of Computational and Graphical Statistics_ , 0(0):1–14, 2022.
  * Solin (2016) Solin, A.  _Stochastic differential equation methods for spatio-temporal Gaussian process regression_.  PhD thesis, Aalto University, 2016.
  * Titsias (2009) Titsias, M.  Variational learning of inducing variables in sparse Gaussian processes.  In _AISTATS_ , pp. 567–574, 2009.
  * Toth & Oberhauser (2020) Toth, C. and Oberhauser, H.  Bayesian learning from sequential data using Gaussian processes with signature covariances.  In _International Conference on Machine Learning_ , pp. 9548–9560. PMLR, 2020.
  * Tsay (2005) Tsay, R. S.  _Analysis of financial time series_.  John wiley & sons, 2005.
  * Williams & Rasmussen (2006) Williams, C. K. and Rasmussen, C. E.  _Gaussian processes for machine learning_.  MIT press Cambridge, MA, 2006.
  * Wilson & Nickisch (2015) Wilson, A. and Nickisch, H.  Kernel interpolation for scalable structured Gaussian processes (KISS-GP).  In Bach, F. and Blei, D. (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 1775–1784, Lille, France, 07–09 Jul 2015\. PMLR.
  * Yu et al. (2019) Yu, H., Chen, Y., Low, B. K. H., Jaillet, P., and Dai, Z.  Implicit Posterior Variational Inference for Deep Gaussian Processes.  In Wallach, H., Larochelle, H., Beygelzimer, A., d’Alché Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 14475–14486. Curran Associates, Inc., 2019.

##  Appendix A Basic Results on Gaussian Distributions

###  A.1 Conditional Gaussian Distribution

Assuming the joint Gaussian distribution
$p({\boldsymbol{\mathbf{x}}})=\mathcal{N}({\boldsymbol{\mathbf{x}}};{\boldsymbol{\mathbf{\mu}}},{\boldsymbol{\mathbf{\Sigma}}})$
over the random vector ${\boldsymbol{\mathbf{x}}}$ such that:

${\boldsymbol{\mathbf{x}}}=\begin{bmatrix}{\boldsymbol{\mathbf{x}}}_{a}\\\
{\boldsymbol{\mathbf{x}}}_{b}\end{bmatrix},\quad{\boldsymbol{\mathbf{\mu}}}=\begin{bmatrix}{\boldsymbol{\mathbf{\mu}}}_{a}\\\
{\boldsymbol{\mathbf{\mu}}}_{b}\end{bmatrix},\quad{\boldsymbol{\mathbf{\Sigma}}}=\begin{bmatrix}{\boldsymbol{\mathbf{\Sigma}}}_{aa}&{\boldsymbol{\mathbf{\Sigma}}}_{ab}\\\
{\boldsymbol{\mathbf{\Sigma}}}_{ba}&{\boldsymbol{\mathbf{\Sigma}}}_{bb},\end{bmatrix}$
(30)

where
${\boldsymbol{\mathbf{\Sigma}}}_{ba}={\boldsymbol{\mathbf{\Sigma}}}_{ab}^{T}$,
we have that the conditional distributions are given by:

$p({\boldsymbol{\mathbf{x}}}_{b}\,|\,{\boldsymbol{\mathbf{x}}}_{a})=\mathcal{N}({\boldsymbol{\mathbf{x}}}_{b};{\boldsymbol{\mathbf{\mu}}}_{b}+{\boldsymbol{\mathbf{\Sigma}}}_{ba}{\boldsymbol{\mathbf{\Sigma}}}_{aa}^{-1}({\boldsymbol{\mathbf{x}}}_{a}-{\boldsymbol{\mathbf{\mu}}}_{a}),{\boldsymbol{\mathbf{\Sigma}}}_{bb}-{\boldsymbol{\mathbf{\Sigma}}}_{ba}{\boldsymbol{\mathbf{\Sigma}}}_{aa}^{-1}{\boldsymbol{\mathbf{\Sigma}}}_{ab}).$
(31)

###  A.2 Expectation over log of Normal Distribution

With an approximate marginal posterior $q({\boldsymbol{\mathbf{f}}}_{\star})$
and a Normal distribution
$p({\boldsymbol{\mathbf{y}}}\,|\,{\boldsymbol{\mathbf{f}}}_{\star})$ of the
form:

$\displaystyle q({\boldsymbol{\mathbf{f}}}_{\star})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{\mu}}}_{\star},{\boldsymbol{\mathbf{\Sigma}}}_{\star}),$
(32) $\displaystyle
p({\boldsymbol{\mathbf{y}}}\,|\,{\boldsymbol{\mathbf{f}}}_{\star})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{y}}};{\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{\Sigma}}}_{y}),$
(33)

we can compute

$\displaystyle\mathbb{E}_{q({\boldsymbol{\mathbf{f}}}_{\star})}\log
p({\boldsymbol{\mathbf{y}}}\,|\,{\boldsymbol{\mathbf{f}}}_{\star})$
$\displaystyle=\log\mathcal{N}({\boldsymbol{\mathbf{y}}};{\boldsymbol{\mathbf{\mu}}}_{\star},{\boldsymbol{\mathbf{\Sigma}}}_{y})+\operatorname{Tr}({\boldsymbol{\mathbf{\Sigma}}}_{y}^{-1}{\boldsymbol{\mathbf{\Sigma}}}_{\star}).$
(34)

##  Appendix B Posterior Marginal in Sparse GP Models

In variational sparse GP models we usually have the joint Gaussian model
$p({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{u}}})$ and a
posterior distribution over $q({\boldsymbol{\mathbf{u}}})$ with

$\displaystyle p({\boldsymbol{\mathbf{f}}}_{\star})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{m}}}_{\star},{\boldsymbol{\mathbf{K}}}_{\star}),$
(35) $\displaystyle p({\boldsymbol{\mathbf{u}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{u}}};{\boldsymbol{\mathbf{m}}}_{Z},{\boldsymbol{\mathbf{K}}}_{Z}),$
(36)
$\displaystyle\mathbb{C}\text{ov}({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{u}}})$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{\star,Z},$ (37) $\displaystyle
q({\boldsymbol{\mathbf{u}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{u}}};{\boldsymbol{\mathbf{\mu}}}_{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{\Sigma}}}_{\boldsymbol{\mathbf{u}}}),$
(38) $\displaystyle
q({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{u}}})$
$\displaystyle:=q({\boldsymbol{\mathbf{u}}})p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,{\boldsymbol{\mathbf{u}}}),$
(39)

and we wish to compute:

$q({\boldsymbol{\mathbf{f}}}_{\star})=\int
q({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{u}}})\operatorname{d}\\!{{\boldsymbol{\mathbf{u}}}},$
(40)

where we have omitted conditional dependencies on
${\boldsymbol{\mathbf{x}}}_{\star}$ and ${\boldsymbol{\mathbf{Z}}}$ for
simplicity in the notation. For example this is necessary to compute the
expectations over the conditional log likelihood term during inferencce or to
make predictions on a new test point. Using Eq. 31 we obtain:

$p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,{\boldsymbol{\mathbf{u}}})=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z}),{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}).$
(41)

In order to integrate out over $q({\boldsymbol{\mathbf{u}}})$ we know that the
result is a Gaussian with mean and covariances obtained from the following
linear transformation of ${\boldsymbol{\mathbf{u}}}$:

${\boldsymbol{\mathbf{f}}}_{\star}={\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z})+{\boldsymbol{\mathbf{\epsilon}}},\quad{\boldsymbol{\mathbf{\epsilon}}}\sim\mathcal{N}(\mathbf{0},{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}).$
(42)

Thus, we obtain:

$\displaystyle q({\boldsymbol{\mathbf{f}}}_{\star})$
$\displaystyle=\mathcal{N}\left({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}({\boldsymbol{\mathbf{\mu}}}_{\boldsymbol{\mathbf{\mu}}}-{\boldsymbol{\mathbf{m}}}_{Z}),{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{\Sigma}}}_{\boldsymbol{\mathbf{u}}}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}+{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}\right)$
(43)
$\displaystyle=\mathcal{N}\left({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}({\boldsymbol{\mathbf{\mu}}}_{\boldsymbol{\mathbf{\mu}}}-{\boldsymbol{\mathbf{m}}}_{Z}),{\boldsymbol{\mathbf{K}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}\left({\boldsymbol{\mathbf{\Sigma}}}_{\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{K}}}_{Z}\right){\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}\right).$
(44)

##  Appendix C Graphical Model

Figure 4: Generative process of the gpssm (left panel) and the sparse gpssm
(right panel).

Fig. 4 illustrates gpssm’s graphical model.

##  Appendix D Details of Prior Reparameterization: Whitening Prior over
Inducing Variables

Here we show a more detailed derivation of the prior reparameterization and,
consequently, the new form of the required conditional distribution. We know
that our prior over the inducing variables is
$p({\boldsymbol{\mathbf{u}}})=\mathcal{N}({\boldsymbol{\mathbf{u}}};{\boldsymbol{\mathbf{m}}}_{Z},{\boldsymbol{\mathbf{K}}}_{Z})$
and that we use the new whitened prior
${\boldsymbol{\mathbf{v}}}={\boldsymbol{\mathbf{L}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z})$
with
${\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{L}}}_{Z}^{T}={\boldsymbol{\mathbf{K}}}_{Z}$
and
$p({\boldsymbol{\mathbf{v}}})=\mathcal{N}({\boldsymbol{\mathbf{v}}};\mathbf{0},{\boldsymbol{\mathbf{I}}}_{M})$.
We also have that the marginal prior over ${\boldsymbol{\mathbf{f}}}_{\star}$
is
$p({\boldsymbol{\mathbf{f}}}_{\star})=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star},{\boldsymbol{\mathbf{K}}}_{\star})$.
It follows that $p({\boldsymbol{\mathbf{f}}},{\boldsymbol{\mathbf{v}}})$ is a
Gaussian with cross-covariance:

$\displaystyle\mathbb{C}\text{ov}({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{v}}})$
$\displaystyle=\mathbb{E}[({\boldsymbol{\mathbf{f}}}_{\star}-{\boldsymbol{\mathbf{m}}}_{\star})({\boldsymbol{\mathbf{L}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z}))^{T}]$
(45)
$\displaystyle=\underbrace{\mathbb{E}[({\boldsymbol{\mathbf{f}}}_{\star}-{\boldsymbol{\mathbf{m}}}_{\star})({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z})^{T}]}_{\mathbb{C}\text{ov}({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{u}}})}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}$
(46)
$\displaystyle={\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}.$
(47)

Hence, using Eq. 31, we have that:

$\displaystyle
p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,,{\boldsymbol{\mathbf{x}}}_{\star},{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{Z}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{L}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,t}),$
(48) $\displaystyle
p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,,{\boldsymbol{\mathbf{x}}}_{\star},{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{Z}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}).$
(49)

Alternatively, we can simply obtain this result by replacing
${\boldsymbol{\mathbf{u}}}$ with
${\boldsymbol{\mathbf{u}}}={\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}}$
in the conditional distribution. Here we note that this is _consistent_ with
the gp definition in that if we were to integrate out the variables
${\boldsymbol{\mathbf{v}}}$ from the joint model
$p({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{v}}})$ we would
obtain the exact gp prior $p({\boldsymbol{\mathbf{f}}}_{\star})$.

###  D.1 Posterior Marginal

With the whitened join model
$p({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{v}}})$, Gaussian
approximate posterior $q({\boldsymbol{\mathbf{v}}})$ and a variationally
sparse gp model we have:

$\displaystyle
p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,{\boldsymbol{\mathbf{v}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{K}}}_{\star}-{\boldsymbol{\mathbf{K}}}_{\star,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}),$
(50) $\displaystyle q({\boldsymbol{\mathbf{v}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{v}}};{\boldsymbol{\mathbf{\mu}}}_{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{\Sigma}}}_{\boldsymbol{\mathbf{v}}}),$
(51) $\displaystyle
q({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{v}}})$
$\displaystyle:=q({\boldsymbol{\mathbf{v}}})p({\boldsymbol{\mathbf{f}}}_{\star}\,|\,{\boldsymbol{\mathbf{v}}}).$
(52)

It is easy to show that the posterior marginal
$q({\boldsymbol{\mathbf{f}}}_{\star})=\int
q({\boldsymbol{\mathbf{f}}}_{\star},{\boldsymbol{\mathbf{v}}})\operatorname{d}\\!{{\boldsymbol{\mathbf{v}}}}$
is:

$\displaystyle
q_{v}({\boldsymbol{\mathbf{f}}}_{\star})=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{\star};{\boldsymbol{\mathbf{m}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{\mu}}}_{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{K}}}_{\star}+{\boldsymbol{\mathbf{K}}}_{\star,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}({\boldsymbol{\mathbf{\Sigma}}}_{\boldsymbol{\mathbf{v}}}-{\boldsymbol{\mathbf{I}}}_{M}){\boldsymbol{\mathbf{L}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,\star}).$
(53)

##  Appendix E Free-Form Posterior Estimation By Bounding the Log Marginal
Likelihood Directly

Computing the log marginal:

$\log
p({\boldsymbol{\mathbf{y}}}_{0:T-1},{\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})=\log\left[p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})p({\boldsymbol{\mathbf{x}}}_{0})p({\boldsymbol{\mathbf{y}}}_{0}\,|\,{\boldsymbol{\mathbf{x}}}_{0})\prod_{t=1}^{T-1}p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\right]+\mathcal{L}({\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}}),$
(54)

where

$\mathcal{L}({\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}}):=\left[\log\int_{{\boldsymbol{\mathbf{f}}}_{1:T-1}}\prod_{t=1}^{T-1}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})\
\mathrm{d}{\boldsymbol{\mathbf{f}}}_{1:T-1}\right],$ (55)

which we note is a log of an expectation, which we can bounded using Jensen’s
inequality:

$\displaystyle\mathcal{L}({\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})$
$\displaystyle\geq\mathbb{E}_{\prod_{t=1}^{T-1}p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{1:t-1},{\boldsymbol{\mathbf{x}}}_{0:t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})}\log\prod_{t=1}^{T-1}p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})$
(56)
$\displaystyle=\sum_{t=1}^{T-1}\mathbb{E}_{p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})}\log
p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t}),$ (57)

where
$p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})$
is the gp predictive distribution over ${\boldsymbol{\mathbf{f}}}_{t}$ at
${\boldsymbol{\mathbf{x}}}_{t-1}$ given the (pseudo observations) inducing
variables ${\boldsymbol{\mathbf{u}}}$ at inducing inputs
${\boldsymbol{\mathbf{Z}}}$:

$\displaystyle
p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})$
$\displaystyle=\mathcal{N}({\boldsymbol{\mathbf{f}}}_{t};{\boldsymbol{\mathbf{A}}}_{t-1}{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{B}}}_{t-1})\quad\text{with
}$ (58) $\displaystyle{\boldsymbol{\mathbf{A}}}_{t-1}$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{t-1,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}$
(59) $\displaystyle{\boldsymbol{\mathbf{B}}}_{t-1}$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{t-1}-{\boldsymbol{\mathbf{K}}}_{t-1,Z}{\boldsymbol{\mathbf{K}}}_{Z}^{-1}{\boldsymbol{\mathbf{K}}}_{Z,t-1}.$
(60)

Hence, the expectations in Eq. 57 can be computed in closed-form:

$\mathbb{E}_{p({\boldsymbol{\mathbf{f}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t-1},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}})}\log
p({\boldsymbol{\mathbf{x}}}_{t}\,|\,{\boldsymbol{\mathbf{f}}}_{t})=\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{A}}}_{t-1}{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1}).$
(61)

Thus, we have that:

$\log
p({\boldsymbol{\mathbf{y}}}_{0:T-1},{\boldsymbol{\mathbf{x}}}_{0:T-1},{\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})\geq\log\left[p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})p({\boldsymbol{\mathbf{x}}}_{0})p({\boldsymbol{\mathbf{y}}}_{0}\,|\,{\boldsymbol{\mathbf{x}}}_{0})\right]+\\\
\sum_{t=1}^{T-1}\left[\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{A}}}_{t-1}{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\right].$
(62)

Then setting the latent variables
${\boldsymbol{\mathbf{\Psi}}}:=\\{{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{0:T-1}\\}$,
we can obtain the approximate log unnormalized posterior:

$\log\tilde{q}({\boldsymbol{\mathbf{\Psi}}}):=\log
p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})+\log
p({\boldsymbol{\mathbf{x}}}_{0})+\log
p({\boldsymbol{\mathbf{y}}}_{0}\,|\,{\boldsymbol{\mathbf{x}}}_{0})\ +\\\
\sum_{t=1}^{T-1}\left[\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{A}}}_{t-1}{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\right].$
(63)

where $\tilde{q}({\boldsymbol{\mathbf{\Psi}}})\propto
q({\boldsymbol{\mathbf{\Psi}}})\approx
p({\boldsymbol{\mathbf{\Psi}}}\,|\,{\boldsymbol{\mathbf{y}}}_{0:T-1})$. Thus,
we can draw samples from our approximate posterior using sghmc and the energy
function $U({\boldsymbol{\mathbf{\Psi}}})=-\log
p({\boldsymbol{\mathbf{\Psi}}},{\boldsymbol{\mathbf{y}}}_{0:T-1})=-\log
p({\boldsymbol{\mathbf{\Psi}}}\,|\,{\boldsymbol{\mathbf{y}}}_{0:T-1})+C\approx-\log\tilde{q}({\boldsymbol{\mathbf{\Psi}}})+C$.

###  E.1 Whitened Version

In the whitened version,
${\boldsymbol{\mathbf{v}}}={\boldsymbol{\mathbf{L}}}_{Z}^{-1}({\boldsymbol{\mathbf{u}}}-{\boldsymbol{\mathbf{m}}}_{Z})$
or
${\boldsymbol{\mathbf{u}}}={\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}}$
and the effective prior is

$p({\boldsymbol{\mathbf{v}}})=\mathcal{N}({\boldsymbol{\mathbf{v}}};\mathbf{0},{\boldsymbol{\mathbf{I}}}).$
(64)

Letting
$\tilde{{\boldsymbol{\mathbf{u}}}}=({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}},{\theta}),\tilde{{\boldsymbol{\mathbf{v}}}}=({\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{Z}}},{\theta})$,
the determinant of the Jacobian can be calculated as:

$\displaystyle|\frac{\partial\tilde{{\boldsymbol{\mathbf{u}}}}}{\partial\tilde{{\boldsymbol{\mathbf{v}}}}}|=$
$\displaystyle\left|\begin{array}[]{ccc}\frac{\partial({\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}})}{\partial{\boldsymbol{\mathbf{v}}}}&\frac{\partial({\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}})}{\partial{\boldsymbol{\mathbf{Z}}}}&\frac{\partial({\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}})}{\partial\theta}\\\
\frac{\partial{\boldsymbol{\mathbf{Z}}}}{\partial{\boldsymbol{\mathbf{v}}}}&\frac{\partial{\boldsymbol{\mathbf{Z}}}}{\partial{\boldsymbol{\mathbf{Z}}}}&\frac{\partial{\boldsymbol{\mathbf{Z}}}}{\partial\theta}\\\
\frac{\partial\theta}{\partial{\boldsymbol{\mathbf{v}}}}&\frac{\partial\theta}{\partial{\boldsymbol{\mathbf{Z}}}}&\frac{\partial\theta}{\partial\theta}\end{array}\right|$
(68) $\displaystyle=$
$\displaystyle\left|\begin{array}[]{ccc}{\boldsymbol{\mathbf{L}}}_{Z}&\frac{\partial({\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}})}{\partial{\boldsymbol{\mathbf{Z}}}}&\frac{\partial({\boldsymbol{\mathbf{m}}}_{Z}+{\boldsymbol{\mathbf{L}}}_{Z}{\boldsymbol{\mathbf{v}}})}{\partial\theta}\\\
0&{\boldsymbol{\mathbf{I}}}&0\\\
0&0&{\boldsymbol{\mathbf{I}}}\end{array}\right|=\left|\begin{array}[]{ccc}{\boldsymbol{\mathbf{L}}}_{Z}&0&0\\\
0&{\boldsymbol{\mathbf{I}}}&0\\\
0&0&{\boldsymbol{\mathbf{I}}}\end{array}\right|=|{\boldsymbol{\mathbf{L}}}_{Z}|=|{\boldsymbol{\mathbf{K}}}_{Z}|^{\frac{1}{2}}$
(75)

Thus, we can re-write the approximate log unnormalized posterior as:

$\log\tilde{q}({\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{x}}}_{0:T-1})=\log|{\boldsymbol{\mathbf{K}}}_{Z}|^{\frac{1}{2}}-\log|{\boldsymbol{\mathbf{K}}}_{Z}|^{\frac{1}{2}}-\frac{1}{2}{\boldsymbol{\mathbf{v}}}^{\top}{\boldsymbol{\mathbf{v}}}+\log
p({\boldsymbol{\mathbf{x}}}_{0})+\log
p({\boldsymbol{\mathbf{y}}}_{0}\,|\,{\boldsymbol{\mathbf{x}}}_{0})+\\\
\sum_{t=1}^{T-1}\left[\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{K}}}_{t-1,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\right].\\\
=\log p({\boldsymbol{\mathbf{v}}})+\log p({\boldsymbol{\mathbf{x}}}_{0})+\log
p({\boldsymbol{\mathbf{y}}}_{0}\,|\,{\boldsymbol{\mathbf{x}}}_{0})+\\\
\sum_{t=1}^{T-1}\left[\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{K}}}_{t-1,Z}({\boldsymbol{\mathbf{L}}}_{Z}^{-1})^{T}{\boldsymbol{\mathbf{v}}},{\boldsymbol{\mathbf{Q}}})-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1})\right].$
(76)

##  Appendix F Collapsed Method

We can integrate out ${\boldsymbol{\mathbf{u}}}$ in Eq. 21 to obtain the
marginal distribution of ${\boldsymbol{\mathbf{x}}}_{1:T}$.

$\displaystyle
q^{*}({\boldsymbol{\mathbf{x}}}_{1:T})=\int_{{\boldsymbol{\mathbf{u}}}}q^{*}({\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{x}}}_{1:T})\mathrm{d}{\boldsymbol{\mathbf{u}}}$
$\displaystyle=\int_{{\boldsymbol{\mathbf{u}}}}p({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{Z}}})p({\boldsymbol{\mathbf{x}}}_{0})\prod_{t=1}^{T}\left[p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\exp(-\frac{1}{2}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1}))\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})\right]\mathrm{d}{\boldsymbol{\mathbf{u}}}$
(77)
$\displaystyle=p({\boldsymbol{\mathbf{x}}}_{0})\prod_{t=1}^{T}\left[p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})\exp(-\frac{1}{2}\
\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}^{-1}{\boldsymbol{\mathbf{B}}}_{t-1}))\right]\int_{{\boldsymbol{\mathbf{u}}}}p({\boldsymbol{\mathbf{u}}})\prod_{t=1}^{T}\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{\mu}}}_{x_{t}},{\boldsymbol{\mathbf{Q}}})\mathrm{d}{\boldsymbol{\mathbf{u}}}.$
(78)

In the case of using our prior re-parameterization, i.e., whitened inducing
variables, ${\boldsymbol{\mathbf{v}}}$ and for multidimensional states and,
therefore, multi-output gps, we can write the integral above for dimension $d$
as:

$\int_{{\boldsymbol{\mathbf{v}}}^{(d)}}p({\boldsymbol{\mathbf{v}}}^{(d)})\prod_{t=1}^{T}\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t}^{(d)};{\boldsymbol{\mathbf{m}}}_{t-1}^{(d)}+\tilde{{\boldsymbol{\mathbf{A}}}}^{(d)}_{t-1}{\boldsymbol{\mathbf{v}}}^{(d)},{\boldsymbol{\mathbf{Q}}}_{d})\mathrm{d}{\boldsymbol{\mathbf{v}}}^{(d)}\\\
=\prod_{t=1}^{T}\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t}^{(d)};{\boldsymbol{\mathbf{m}}}_{t-1}^{(d)},{\boldsymbol{\mathbf{Q}}}_{d})\exp\left[\frac{1}{2}({\boldsymbol{\mathbf{g}}}^{(d)})^{\top}({\boldsymbol{\mathbf{H}}}^{(d)})^{-1}{\boldsymbol{\mathbf{g}}}^{(d)}\right]\left(\det({\boldsymbol{\mathbf{H}}}^{(d)})\right)^{\frac{1}{2}}$
(79)

where

$\displaystyle{\boldsymbol{\mathbf{g}}}^{(d)}$
$\displaystyle:=\left({\boldsymbol{\mathbf{I}}}+\sum_{t=1}^{T}(\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)})^{\top}{\boldsymbol{\mathbf{Q}}}_{d}^{-1}(\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)})\right)^{-1}\left(\sum_{t}({\boldsymbol{\mathbf{x}}}^{(d)}_{t}-{\boldsymbol{\mathbf{m}}}_{t-1}^{(d)})^{\top}{\boldsymbol{\mathbf{Q}}}_{d}^{-1}\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)}\right)^{\top},$
(80) $\displaystyle{\boldsymbol{\mathbf{H}}}^{(d)}$
$\displaystyle=\left({\boldsymbol{\mathbf{I}}}+\sum_{t=1}^{T}(\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)})^{\top}{\boldsymbol{\mathbf{Q}}}_{d}^{-1}(\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)})\right)^{-1}$
$\displaystyle\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)}$
$\displaystyle={\boldsymbol{\mathbf{K}}}_{t-1,Z}^{(d)}(({\boldsymbol{\mathbf{L}}}_{Z}^{(d)})^{\top})^{-1},$
(81) $\displaystyle{\boldsymbol{\mathbf{K}}}_{Z}^{(d)}$
$\displaystyle={\boldsymbol{\mathbf{L}}}_{Z}^{(d)}({\boldsymbol{\mathbf{L}}}_{Z}^{(d)})^{\top}.$
(82)

Thus, our optimal log variational posterior marginal over state trajectories
${\boldsymbol{\mathbf{x}}}_{0:T}$ is

$\log q^{*}({\boldsymbol{\mathbf{x}}}_{0:T})={\log
p({\boldsymbol{\mathbf{x}}}_{0})}+\sum_{t=1}^{T}\left[{\log
p({\boldsymbol{\mathbf{y}}}_{t}\,|\,{\boldsymbol{\mathbf{x}}}_{t})}+\log\mathcal{N}({\boldsymbol{\mathbf{x}}}_{t};{\boldsymbol{\mathbf{m}}}_{t-1},{\boldsymbol{\mathbf{Q}}})\right]\\\
-\frac{1}{2}\sum_{d=1}^{d_{x}}\left[\sum_{t=1}^{T}\operatorname{Tr}({\boldsymbol{\mathbf{Q}}}_{d}^{-1}{\boldsymbol{\mathbf{B}}}^{(d)}_{t-1})+\log\det({\boldsymbol{\mathbf{H}}}^{(d)})^{-1}-(\tilde{{\boldsymbol{\mathbf{x}}}}^{(d)})^{\top}{\boldsymbol{\mathbf{H}}}^{(d)}(\tilde{{\boldsymbol{\mathbf{x}}}}^{(d)})\right],$
(83)

where
$\tilde{{\boldsymbol{\mathbf{x}}}}^{(d)}:=\sum_{t=1}^{T}(\tilde{{\boldsymbol{\mathbf{A}}}}_{t-1}^{(d)})^{\top}{\boldsymbol{\mathbf{Q}}}_{d}^{-1}({\boldsymbol{\mathbf{x}}}_{t}^{(d)}-{\boldsymbol{\mathbf{m}}}_{t-1}^{(d)}).$

##  Appendix G Closed-Form Optimal Variational Conditional Distribution over
Inducing Variables

Given a trajectory, the optimal variational conditional over the whitened
inducing variables given a state trajectory is given by:

$q({\boldsymbol{\mathbf{v}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})=\prod_{d=1}^{d_{x}}\mathcal{N}({\boldsymbol{\mathbf{v}}}^{(d)};{\boldsymbol{\mathbf{g}}}^{(d)},{\boldsymbol{\mathbf{H}}}^{(d)}),$
(84)

where ${\boldsymbol{\mathbf{g}}}$ and ${\boldsymbol{\mathbf{H}}}$ are given in
Eqs. 80 and F, respectively. A similar expression can be obtained for the
original (non-whitened) inducing variables ${\boldsymbol{\mathbf{u}}}$. Of
importance here is that our formulation has allowed us to obtain an optimal
conditional variational distribution, given the state trajectories, over the
inducing variables in closed form. We see then that, in effect, our
approximate posterior is
$q_{\lx@glossaries@gls@link{acronym}{FFVD}{\leavevmode\textsc{ffvd}}}({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}})=q({\boldsymbol{\mathbf{x}}}_{0:T})q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$.
In order to make predictions, we run sghmc to obtain samples
$\\{{\boldsymbol{\mathbf{x}}}_{0:T}^{(s)}\\}$ from the marginal
$p({\boldsymbol{\mathbf{x}}}_{0:T})$ and then use the closed-form expression
for the conditional to obtain
$\\{{\boldsymbol{\mathbf{u}}}^{(s)}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T}^{(s)}\\}$.
Running sghmc over the much lower-dimensional space of state trajectories
(instead of the joint space of trajectories and inducing variables) should
generally converge faster. Or experiments with both versions of the algorithm
confirm this.

###  G.1 Optimal Closed-Form Conditional vs Assumed Parametric Factorizations

We contrast our variational distribution with that of the variational Gaussian
process state space model (vgpssm; Frigola et al., 2014) and the variationally
coupled dynamic trajectories (vcdt; Ialongo et al., 2019). While, implicitly,
our variational distribution is
$q_{\lx@glossaries@gls@link{acronym}{FFVD}{\leavevmode\textsc{ffvd}}}({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}):=q({\boldsymbol{\mathbf{x}}}_{0:T})q({\boldsymbol{\mathbf{u}}}\,|\,{\boldsymbol{\mathbf{x}}}_{0:T})$
the vgpssm and vcdt assume
$q_{\lx@glossaries@gls@link{acronym}{VGPSSM}{\leavevmode\textsc{vgpssm}}}({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}):=q({\boldsymbol{\mathbf{x}}}_{0:T})q({\boldsymbol{\mathbf{u}}})$
and
$q_{\lx@glossaries@gls@link{acronym}{VCDT}{\leavevmode\textsc{vcdt}}}({\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}}):=q({\boldsymbol{\mathbf{x}}}_{0:T}\,|\,{\boldsymbol{\mathbf{u}}})q({\boldsymbol{\mathbf{u}}})$,
respectively. Although the former (vgpssm) obtained an “optimal” posterior
(given the factorization assumption), it is a mean-field approximation and
ignores the dependencies between state trajectories and inducing variables.
The latter (vcdt), does not make a mean-field assumption but its factorization
forces parametric constraints over the individual distributions. In fact,
Ialongo et al. (2019) assume Gaussian posteriors for both the marginal
$q({\boldsymbol{\mathbf{u}}})$ and each of the time-dependent conditionals in
$q({\boldsymbol{\mathbf{x}}}_{0:T}\,|\,{\boldsymbol{\mathbf{u}}})$. Our
proposal is the only variational distribution that is theoretically optimal
while yielding a closed-form conditional.

##  Appendix H Particle Markov chain Monte Carlo (mcmc)

We can use the sequential structure of the latent states
${\boldsymbol{\mathbf{x}}}_{0:T}$ for efficient inference. More specifically,
using the transition in Eq. 22. The details of this are given in Algorithm 1.

Algorithm 1 Particle MCMC for inferring latent states
${\boldsymbol{\mathbf{x}}}_{0:T}$

0: Number of particles $S$, observations ${\boldsymbol{\mathbf{y}}}_{1:T}$,
samples of ${\boldsymbol{\mathbf{x}}}_{1:T}^{\prime}$ in the previous
iteration, likelihood parameters
${\boldsymbol{\mathbf{\phi}}}=\\{{\boldsymbol{\mathbf{C}}},{\boldsymbol{\mathbf{d}}},{\boldsymbol{\mathbf{R}}}\\}$,

Initialize state ${\boldsymbol{\mathbf{x}}}_{0}$, weights $W^{(i)}_{t}=1$ for
$i=1,\ldots,S$.

Fix the last particle as setting
${\boldsymbol{\mathbf{x}}}_{1:T}^{(S)}={\boldsymbol{\mathbf{x}}}_{1:T}^{\prime}$

for $t=1,\ldots,T$ do

Generate ${\boldsymbol{\mathbf{x}}}_{t}^{(i)}$ from Eq. 22 for
$i=1,\ldots,S-1$.

Calculate weight
$W_{t}^{i}=p_{{\boldsymbol{\mathbf{\phi}}}}({\boldsymbol{\mathbf{y}}}_{t}|{\boldsymbol{\mathbf{x}}}_{t}^{(i)})$
for $i=1,\ldots,S$

Normalize the weights as
$\overline{W}_{t}^{(i)}=W_{t}^{(i)}/\sum_{i}W_{t}^{(i)}$

if $t<T$ then

For $i^{\prime}=1,\ldots,S-1$, re-sample the index $j_{i^{\prime}}$ from the
categorical distribution, with the event probabilities being
$(\overline{W}_{t}^{(1)},\ldots,\overline{W}_{t}^{(S)})$.

For $i^{\prime}=1,\ldots,S-1$, let
${\boldsymbol{\mathbf{x}}}_{1:t}^{(i^{\prime})}:={\boldsymbol{\mathbf{x}}}_{1:t}^{(j_{i^{\prime}})}$
and set $W_{t}^{i^{\prime}}=1/S$

else

Re-sample the index $j^{*}$ from the categorical distribution, with the event
probabilities being $(\overline{W}_{t}^{(1)},\ldots,\overline{W}_{t}^{(S)})$.

Return
${{\boldsymbol{\mathbf{x}}}}_{0:T}:=={\boldsymbol{\mathbf{x}}}_{1:T}^{(j^{*})}$

end if

end for

##  Appendix I Experiment Details

Regarding the variational Gaussian process state space model (vgpssm),
probabilistic recurrent state space model (prssm), and vcdt, we download the
authors’ implementations from their websites. For a fair comparison, we try to
follow the same treatment for hyper-parameters across all methods. We note,
however, that the Matlab implementation of vgpssm does not provide inference
for gp hyper-parameters ${\boldsymbol{\mathbf{\theta}}}$, conditional
likelihood parameters ${\boldsymbol{\mathbf{\phi}}}$, inducing inputs
${\boldsymbol{\mathbf{Z}}}$, process variance $Q$ or observation variance
${\boldsymbol{\mathbf{R}}}$.

###  I.1 Settings for synthetic data

In generating the synthetic data, we set the kernel’s signal variance
$\sigma=2.0$ and lengthscale $l=0.5$ in the sparse GPSSM. we reduce the impact
of observational error and set the observation variance as small as
$\sigma^{2}=0.01$, i.e., each observation would be generated as
$y_{t}\sim\mathcal{N}(x_{t},0.01)$. We also set the process variance
${\boldsymbol{\mathbf{Q}}}=0.01$. The number of inducing points is set to
$20$, and these $20$ inducing points are evenly spread in the interval
$[-2,2]$. We set the length of the training trajectory as $120$, the number of
iterations as $50\,000$, and the number of posterior samples as $50$.

####  I.1.1 Testing the Marginals for Gaussianity

We conduct a hypothesis test for each individual trajectory state and inducing
variable. Given the $50$ posterior samples for ${\boldsymbol{\mathbf{x}}}_{t}$
(or ${\boldsymbol{\mathbf{u}}}_{m}$), which we denote them as
$\\{{\boldsymbol{\mathbf{x}}}_{t}^{(s)}\\}_{s=1}^{50}$ (or
$\\{{\boldsymbol{\mathbf{u}}}_{m}^{(s)}\\}_{s=1}^{50}$), we use the
implementation of _scipy.stats.normaltest_ from Python’s _scipy_ package to
test whether we have sufficient evidence to reject the hypothesis that
${\boldsymbol{\mathbf{x}}}_{t}$ (or ${\boldsymbol{\mathbf{u}}}_{m}$) follows a
Gaussian distribution. _scipy.stats.normaltest_ is based on the work of
D’agostino & Pearson (1973).

###  I.2 Settings for Real-world Data

We adopt similar settings as in Ialongo et al. (2019) to initialize the hyper-
parameters for all the models, which first uses a factorized nonlinear model
to optimize
${\boldsymbol{\mathbf{\theta}}},{\boldsymbol{\mathbf{x}}}_{0:T},{\boldsymbol{\mathbf{u}}},{\boldsymbol{\mathbf{Z}}},{\boldsymbol{\mathbf{C}}},{\boldsymbol{\mathbf{d}}},{\boldsymbol{\mathbf{R}}},{\boldsymbol{\mathbf{Q}}}$.
We use the identity function as the mean function and the squared exponential
function with automatic relevance determination (ARD) in gpssm. When dealing
with more than $1$ latent dimension, which requires a multi-output gp, we use
a different set of kernel hyper-parameters for each output. We set the number
of inducing points to $M=100$ and the number of dimensions for latent states
${\boldsymbol{\mathbf{x}}}_{0:T}$ to $d_{x}=4$. We set standard diagonal
Gaussian priors for the initial latent state ${\boldsymbol{\mathbf{x}}}_{0}$
and (where applicable) for likelihood parameters
${\boldsymbol{\mathbf{C}}},{\boldsymbol{\mathbf{d}}}$, the logarithm of
observation standard deviation $\log({\boldsymbol{\mathbf{R}}})/2$ and process
variance $\log({\boldsymbol{\mathbf{Q}}})$.

For optimizing these hyper-parameters, we use Adam, with default settings for
the optimizer parameters except for a decayed learning rate. We run our ffvd
algorithm for $50\,000$ iterations and vcdt for $200\,000$. We set the
learning rate as $0.01$ and the decay parameter as $0.05$ during sghmc
sampling. We used $S=10^{5}$ posterior samples for vcdt (as in the results in
the original paper) and $S=100$ for ours. Following the evaluation in Ialongo
et al. (2019), we use the rmse between the models’ predictions and the ground
truth in the future $30$ time steps (not seen in training) as the comparison
criterion for all the methods.

Regarding the benchmarks, we have $1\,024$ observations for _Actuator_ ,
$1\,000$ observations for _Ballbeam_ , $500$ observations for _Drive_ ,
$1\,000$ observations for _Dryer_ , $1\,024$ for _Flutter_ , and $296$
observations for _Furnace_. Training lengths are:
$T_{\text{train}}=500,500,250,500,500,150$, respectively and test lengths are
the rest of the sequences. All these observations are $1$-dimensional but, as
in previous work, we consider $d_{x}=4$-dimensional latent states. Each
benchmark also contains a $1$-dimensional control input vector
${\boldsymbol{\mathbf{a}}}\in\mathbb{R}^{T\times 1}$.

###  I.3 Reproducibility

We will make our code publicly available upon acceptance of the paper.

##  Appendix J Additional Results

###  J.1 Training and Test Performance on All Benchmarks

Fig. 5 shows the training and testing regimes and predictions for the six
system identification benchmarks considered.

Figure 5: Training and test performance on benchmarks. Row 1 (left to right):
_Actuator_ , _Ballbeam_ , _Drive_ ; row 2 (left to right): _Dryer_ , _Flutter_
, _Furnace_.

###  J.2 Traceplot of Log-Likelihood on Training data for ffvd-m and ffvd-c-m

Figure 6 shows the convergence trace plots for all benchmarks.

Figure 6: Traceplot of log-likelihood on training data for FFVD-M and FFVD-
C-M. Row 1 (left to right): _Actuator_ , _Ballbeam_ , _Drive_ ; row 2 (left to
right): _Dryer_ , _Flutter_ , _Furnace_.

##  Appendix K Extended Related Work

Here we give more details about the differences between our method and closely
related approaches for inference in the gpssm. The vgpssm of Frigola et al.
(2014) is the first developing variational approaches to inference in gpssms,
and it uses mean-field assumptions to factorize the variational distribution
of inducing variables ${\boldsymbol{\mathbf{u}}}$ and latent state
trajectories ${\boldsymbol{\mathbf{x}}}_{0:T}$ as
$q_{\lx@glossaries@gls@link{acronym}{VGPSSM}{\leavevmode\textsc{vgpssm}}}({{\boldsymbol{\mathbf{u}}}},{\boldsymbol{\mathbf{x}}}_{0:T})=q({{\boldsymbol{\mathbf{u}}}})q({\boldsymbol{\mathbf{x}}}_{0:T})$.
This assumption clearly overlooks the complex dependencies between
${{\boldsymbol{\mathbf{u}}}}$ and ${\boldsymbol{\mathbf{x}}}_{0:T}$. Doerr et
al. (2018) model these dependencies with their prssm algorithm, but make the
unrealistic assumption of the posterior dynamics over the latent state
trajectories being the same as the prior. Ialongo et al. (2019) show that this
assumption can have critical consequences on the posterior and predictive
distributions and propose a more general algorithm called vcdt. Although a
significant improvement over previous methods, vcdt assumes a Gaussian
distribution over the inducing variables and linear-Gaussian transitions over
the state posterior. Here we lift these assumptions and infer the optimal
joint posterior under the variational framework.

As mentioned in the main paper, Frigola et al. (2013) also proposed a pmcmc
algorithm for inference over state trajectories. However, the algorithm is
based on the fully-independent conditional approximation (fic; Quinonero-
Candela & Rasmussen, 2005). As reported in Frigola et al. (2013), a naïve
implementation of their algorithm has a time complexity of
$\mathcal{O}(M^{2}T^{2})$, which is significantly higher than ours. A much
more efficient implementation with $\mathcal{O}(M^{2}T)$ time complexity that
requires significant tracking of intermediate data structures and
factorizations is hinted at in the original paper. Unfortunately, the
corresponding code has not been made publicly available. Furthermore, model
approximations such as those in fic can have surprising consequences, such as
incorrectly estimating the noise variance to be almost zero or ignoring
additional inducing inputs, see, e.g., Bauer et al. (2016) for a thorough
discussion. Consequently, we see our variational approach as more principled.