Compositional Modeling of Nonlinear Dynamical Systems with ODE-based Random Features
Authors: Thomas M. McDonald
Mauricio A Álvarez
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian processes

Abstract: Effectively modeling phenomena present in highly nonlinear dynamical systems whilst also accurately quantifying uncertainty is a challenging task, which often requires problem-specific techniques. We present a novel, domain-agnostic approach to tackling this problem, using compositions of physics-informed random features, derived from ordinary differential equations. The architecture of our model leverages recent advances in approximate inference for deep Gaussian processes, such as layer-wise weight-space approximations which allow us to incorporate random Fourier features, and stochastic variational inference for approximate Bayesian inference. We provide evidence that our model is capable of capturing highly nonlinear behaviour in real-world multivariate time series data. In addition, we find that our approach achieves comparable performance to a number of other probabilistic models on benchmark regression tasks.

Compositional Modeling of Nonlinear Dynamical Systems with ODE-based Random Features
Thomas M. McDonald Department of Computer Science University of Sheffield tmmcdonald1 @sheff ac. uk
Mauricio A. Álvarez Department of Computer Science University of Sheffield
mauricio wlwaareesseee
Abstract
Effectively modeling phenomena present in highly nonlinear dynamical systems whilst also accurately quantifying uncertainty challenging task, which often re- quires problem-specific techniques. We present novel, domain-agnostic approach tackling this problem, using compositions of physics-informed random features, derived from ordinary differential equations. The architecture of our model lever- ages recent advances in approximate inference for deep Gaussian processes, such as layer-wise wise weight-space approximations which allow us to incorporate random Fourier features, and stochastic variational inference for approximate Bayesian inference. We provide evidence that our model is capable of capturing highly nonlinear behaviour real-world multivariate time series data. In addition, we find that our approach achieves comparable performance to number of other probabilistic models on benchmark regression tasks.
Introduction
Dynamical systems are ubiquitous across the natural sciences, with many physical and biological processes being driven on a fundamental level by differential equations. Inferring ordinary differential equation (ODE) parameters using observational data from such systems is an active area research [Meeds et al. 2019 Ghosh et al. 2021], however, in particularly complex systems it is often infeasible to characterise all of the individual processes present and the interactions between them. Rather than attempt to fully describe complex system, latent force models (LFMs) [Alvarez et al. 2009] specify a simplified mechanistic model of the system which captures salient features of the dynamics present. This leads model which is able readily extrapolate beyond the training input space, thereby retaining one the key advantages of mechanistic modeling over purely data-driven techniques.
Modeling nonlinear dynamical systems presents an additional challenge, and whilst nonlinear differential equations have been incorporated into LFMs [Ward et al. 22200,, shallow models such as LFMs are generally less capable than deep models of modeling the non-stationarities often present in nonlinear systems. Deep model architectures have greater representational power than shallow models result of their hierarchical structure, which typically consists of compositions of functions [LeCun et al. 2015) Deep probabilistic models such as deep Gaussian processes (DGPs) [Damianou and Lawrence 2013 able to leverage this representational power, with the additional advantage of being able to accurately quantify uncertainty.
In this paper, we aim to model nonlinear dynamics by constructing a DGP from compositions of physics-informed random Fourier features, with the motivation behind this approach being that many real-world systems are compositional hierarchies [Leuun et al. 2015). Whilst number of recent works have incorporated random Fourier features into deep models [Cutajar et al. 2017 Mehrkanoon and Suykens 2018), the only work are aware of which does so in the context of a physics-informed DGP is that Lorenzi and Filippone 2018). The authors impose physical structure on a DGP with
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

2 f1
f2
fL
(a) Deep Gaussian process (DGP)
2' fi
f2
ft
fr.
(e,
0 U1.
2 (..
(r.
(b) Deep latent force model (DLFM)
Figure 1: A conceptual explanation of how our proposed DLFM differs from DGP. At each layer, we perform the operation ( {xr, G} Sod G r)u(T)dT,, where G' is the Green's function corresponding to an ODE, and u(-) represents an exponentiated quadratic GP prior. For example, the second operation in the model shown above would take the form, ( (f1, G} = S6' - (f))()))
random features by constraining the function values within the model means of performing ODE parameter inference. Rather than following the approach of Lorenzi and Filippone (0118) and placing constraints on function values whilst using exponentiated quadratic (EQ) kernel, we instead integrate a -nnfored prior into the structure the DGP by utilising a kernel based on an ODE.
Our main contribution in this paper the introduction of novel approach to incorporating physical structure into a deep probabilistic model, whilst providing sound quantification of uncertainty. We achieve this through derivation of physics-informed random Fourier features via the convolution of an EQ GP prior with the Green's 's function associated with first order ODE. These features are then incorporated into each layer of a DGP, as shown in Figure [ To ensure the scalability of our model to large datasets, stochastic variational inference is employed inference. We provide experimental evidence that our modeling framework is capable of capturing highly nonlinear dynamics effectively in both toy examples and real world data, whilst also being applicable to more general regression problems.
method for approximate Bayesian
2 Background
In this section, we review the theory behind using random Fourier features for deep Gaussian processes. We also review the idea of building covariance functions for GPs that encode the dynamics of ODEs and how compute those covariances using random Fourier features.
2.1 Deep Gaussian Processes with Random Feature Expansions
Gaussian processes (GPs) are non-parametric probabilistic models which offer a great degree of flexibility in terms of the data they can represent, however this flexibility is not unbounded. The hierarchical networks of GPs now widely known as deep Gaussian processes (DGPs) were first formalised by Damianou and Lawrence 2013). with the motivating factor behind their ere- ation being the ability of deep networks to reuse features and allow for higher levels of ab- straction in said features (Bengio et al. 2013). which results in such models having more rep- resentational power and flexibility than shallow models such as GPs. DGPs are effectively a composite function, where the input the model transformed to the output by being passed through a series of latent mappings (i.e. multivariate GPs). If we consider a supervised learn- ing problem with inputs denoted by X {xn) 1 and targets denoted by y {ur we can write the analytic form of the marginal likelihood for a DGP with N₂, hidden layers as _[ p(F(1) ((0) (dF(Ns) (..) where p(y X, 0) = f[y(( (p))) ( (
( and g( represent the latent values and covariance parameters respectively at the L-th layer, where .N,, However, due the need propagate densities through non-linear GP covari- ance functions within the model, this integral is intractable [Damianou] 2013]. As this precludes us from employing exact Bayesian inference in these models, various techniques for approximate inference have been applied to DGPs in recent years, with most approaches broadly based upon either variational inference (Salimbeni and Deisenroth 2017 Salimbeni et al. 2019 Yu et al. 2019 or
Monte Carlo methods [Havasi et al. 2018

---

Cutajar et al. [ outline an alternative approach to tackling this problem, which involves replacing the GPs present at each layer the network with their two layer weight-space approximation, forming a Bayesian neural network which acts as an approximation to a DGP. If the E-th layer of such DGP (consisting of zero-mean GPs with exponentiated quadratic kernels) receives an input F(4) (where F(O) X), the random features for this layer are denoted by ((()) RR Nx2N (where JRr denotes the number of random features used), and are given by,
(Fs(F(()) (F(() (2)))) sin(F())
NO
where (a²)(0) is the marginal variance kernel hyperparameter, N10)
"RF is the number of random features
used and (2) E P N is the matrix of spectral frequencies used to determine the random features. This matrix assigned a prior ) =V 0, (1(9) =3 where ((0) is the lengthscale kernel hyperparameter, Deco is the number of GPs within the layer, and df = 1,.,,, Dyco. The random features then undergo the linear transformation, F((+1) (((),,,, where W(() E
R2N xD p((11) is a weight matrix with each column assigned a standard normal prior. Training is achieved via stochastic variational inference, which involves establishing tractable lower bound for the marginal likelihood and optimising said bound with respect to the mean and variance of the variational distributions over the weights and spectral frequencies across all layers of the network. The bound is also optimised with respect to the kernel hyperparameters across all layers.
2.2 Random Feature Approximations for Latent Force Models
In this paper, we take a similar approach to Cutajar et al. (017). but we utilise physically-inspired random features. Physies-inspired covariance functions for GPs have been proposed by [Alvarez et al. 2009] under the name of latent force models, and more recently, the authors of Guarnizo and Alvarez 2018 studied the construction of random features related to such covariance functions.
Latent force models (LFMs) [Alvarez et al. 2009 are GPs which incorporate a physically-inspinn kernel function, which typically encodes the behaviour described by a specific form of differential equation. Instead of taking a fully mechanistic approach and specifying all the interactions within a physical system, the kernel describes a simplified form of the system in which the behaviour is determined by Q latent forces. Given an input data-point I and a set of D output variables (fa())
a LFM expresses each output as fa(t) = Σの 1 Sa.a Ie (T(t-))u((rdd where Ga(+) represents the Green's 's function for a certain form of linear differential equation. (t(() GP(O, kg (t, t')) represents the GP prior over the the q-th latent force, and Sd.q is a sensitivity parameter weighting the influence of the q-th latent force on the d-th output. The general expression for the covariance of
a LFM is given by ffa.far (t, t') = Ed=1 Sa.g .g fo Galt - So Gz - (Tg +') dr.
Typically, an EQ form is assumed for the kernel governing the latent forces, (.). Due to the linear nature of the convolution used to compute f(((), the outputs can also be described by GP. Exact inference in LFMs is tractable, however as with GPs, scales with ((N3) complexity [Rasmussen and Williams 2006). Furthermore, computing the double integration for ffa.... (t, t') leads to terms that include the numerical approximation of complex functions which burden the computation of the kernel function. By providing random Fourier feature representation for the EQ kernel (t, t"). Guarnizo and Álvarez (018)) were able reduce this cubic dependence on the number of data-points to a linear dependence. This representation arises from Bochner's 's theorem, which states,
kg(T, kg(T,T')=e =e -(r-r'12
-f p(w)eilt- de
NRF
Σ jUST e iius"
(2)
NRF
where Ey is the lengthscale of the EQ kernel, NR is the number of random features used in the approximation and We ~ p(w) = N(w 2/e). Substituting this form of the EQ kernel into kfarfa (t,)') leads to fast approximation to the LFM covariance kfa (t, se
Sa NEF salt, Oa, wa)d (t', Ov w..). where EDAAr
NRF
dalt, Oa, = 1/' Galt r) dr,
(3)
3

---

with 0g representing the Green's function parameters. When the Green's 's function is a real function,
(t'(', (r.w) = or ,t',, 0x, -w). Guarnizo and Alvarez (2018) refer to @a(t, Od w) as random Fourier response features (RFRFs).
3 Deep Latent Force Models
In this section, we present the deep latent force model (DLFM), a novel approach to incorporating physically-inspired prior beliefs into a deep Gaussian process. Firstly, we will outline the model architecture, before discussing our approach to training the model via stochastic variational inference.
3.1 Model Formulation
Rather than deriving the random features ((() within the DGP from an EQ kernel (Cutajar et al. 2017]. we instead populate this matrix with features derived from an LFM kernel. The exact form of the features derived is dependent on the Green's function used, which in turn depends on the form of differential equation whose characteristics we wish to encode within the model. Guarnizo and Álvarez (0188) derived number of different forms corresponding to various differential equations, with the simplest case being that of a first order ordinary differential equation (ODE) of the form, df(t) + +() = Σg- (tu,()), where 7 is decay parameter associated with the ODE and Sy is sensitivity parameter which weights the effect of each latent force. For simplicity of exposition, we assume D = 1. The Green's 's function associated to this ODE has the form G(x) = e e_ By
using equation 33. the RFRFs associated with this ODE follow as
juHt 72
ott, w,,)
(4)
+ +jj
where compared to the expression (3). the parameter 0 of the Green's function corresponds to the decay parameter y. From here on, we redefine the spectral frequencies as wq.s to emphasise the fact that the values sampled are dependent on the lengthscale of the latent force by way of the prior, 4ig.s ~ (w 2/ 2/()). We can collect all of the random features corresponding to the q-th latent force into a single vector, (t 7: wyy = / / / [neloltt ", Wg,1,, ... ,tt, y, (g.N 11 E Narx1, where w/ (w ) C refers to the complex plane and the super index in dG(()) makes explicit that this vector contains complex-vaie valued numbers. By including the random features corresponding to all Q latent forces within the model, obtain & (t, 7, w) [( 7, w1) ..... (tg 7, wQ)) '1T CQNa=x1 x1, with w [ These random features will be denoted as (i. 7, w) to differentiate them from the features (6.) computed from a generic EQ kernel.
nigher-dimensional inputs Although the expression for $i.F (t, 7, w) was obtained in the context of an ODE where the input is the time variable, in this paper we exploit this formalism to propose the use of these features even in the context of a generic supervised learning problem where the input is a potentially high high-ddmeesioonl vector later, such an extension is also necessary we attempt to use such features at intermediate layers of the composition. Essentially, leading to a set of vectors (Φr,Far (,,,7, 7,, W1), &i. (Tp, Tp: Wp)}. Notice that the samples wm also be different per input dimension, I'. Although there different ways in which these feature vectors can be combined, in this paper, we assume that the final random feature vector computed over the whole input vector x is given ΦLFM (x
[x1, 22, Fpp
R As will be noticed
compute vector $LFM (Im, Yma Wm) for each input dimension Vi
Σm $LFM (Ims Tm wm). where Ω = {wm 1m==1 and 7 = {7m P m= An alternative to explore for future work involves expressing of f (x: 52 as df (x, 7, S2) -E Omm &f Ff. (Xmy Tm, Wm), with Onn E IR parameter that weights the contribution of each imput feature differently. Although we allow each input dimension have a different decay parameter 7' the experiments in Section [5] for ease of notation we will assume that 71 Therefore, & (x, 52) is a vector-valued function that maps from IRxX1 to C
= = 7p = y. For simplicity, we write ΦLFM (x, 7, Ω).
Real version of the RFRFs Rather than working with the complex-value random features ΦirM (x, 7, 22), we can work with their real-valued counterpart by using (x,

---

(
(a) Deep Gaussian process (DGP)
C
Lue (ar)
f(
em
P (fra)
f
Ouseleent well))
Purnlle 2 mgol
Muselff (2 4211
Dreudes 25 253
(b) Deep latent force model (DLFM)
Figure 2: An illustration of how our proposed model differs from a DGP with random feature expansions, with this example containing dimension, NE random features of the form shown in (4) are computed for each the Q latent
layers. At each layer of the DLFM, for each input
(e) forces. The random feature vector LLFM is then formed by taking the sum of these features across the input dimensions. This summation is shown in the Figure by the block containing E.
[Rif (x., (2))) (Jm (mm(@LF (L 7x E 2 Guarnizo and Álvarez 2018 where Re(a) and Jm(a) take the real component and imaginary component of a, respec- tively. For an input matrix X [x..... _NN] RNxp, the matrix LLFM S2) (ØLFM(XI, S2,.)... ,LL. (XN, 2))T RX2QNRR
Composition of RFRFs We now use + LFA (X; 7, 2) as a building block of a layered architecture of RFRFs. We write & (e)
LEM (F(t) (0) (2)) N x N and follow similar construction
to the one described in Section 2.1 where F((+1) PI (F(*), (0), (2() (((),()) and W e
22Q") Nie p(e+))) As before, F00) = - Figure 2 is an example of this compositional architecture of RFRFs, which we refer to a deep latent force model (DLFM). When considering multiple-output problems, we allow extra flexibility to the decay parameters and lengthscales at the final (L-th) layer such that they vary not only by input dimension, but also by output. Mathematically, this corresponds (L-1) s2 (2(1--1) ( (L-1) where d = l.... D.
to computing pLL) 4 (L-1) (F((), 7
Model interpretation In DGP, the layers of GPs successively warp the input space, propagating uncertainty through the model. The DLFM also performs a similar warping of the inputs, however, in addition, the Green's function involved in the convolution which takes place within model can be interpreted as filter, which selectively filters out some the frequencies of the GPs at each layer. Therefore, the DLFM is not only warping the intermediate inputs, but also filtering components in the latent GPs, altering their degree of smoothness.
3.2 Variational Inference
As previously mentioned, exact Bayesian inference is intractable for models such ours, therefore in order train the DIFM we employ stochastic variational inference [Hoffman et al.] 2013). For notational simplicity, let W, Ω and @ represent the collections of weight matrices, spectral frequencies and kernel hyperparameters respectively, across all layers the model. We seek to optimise the variational distributions over Q and W whilst also optimising O, however we do not place a variational distribution over these hyperparameters. Our approach resembles the VAR-FIXED training strategy described by Cutajar et al. 2017 which involves reparameterizing 2 such that (E) (E)
(2) 2 Sij Erij + m. where mi and (s (2)
represent the means and variances associated with

---

the variational distribution over 2.3 and ensuring that the standard normal samples Frij are fixed throughout training rather than being resampled at each iteration. If we denote v = (WW S2) and consider training inputs X E RN Die and outputs y E IR Nxour, we can derive a tractable lower bound on the marginal likelihood using Jensen's 's inequality, which allows for minibatch training using a subset of M observations from the training set of N total observations. This bound, derived by
(e)
Cutajar et al. 2017], takes the form,
log(p(y|X, O E log[p(y(X, V. O) (KLlg())lll(()))
(5)
Nuc
Σ
loglp(y>|xe V., 6) DKLla [ (o)ll
(6)
M
NMC
KEZMr
where DKL. denotes the KL divergence, V q( v) the minibatch input space is denoted by Za and NMC is the number of Monte Carlo samples used to estimate Eg(D). o( v()) and ( () denote the approximate variational distribution and the prior distribution over the parameters respectively, both of which are assumed to be Gaussian in nature. A full derivation of this bound and the expression for the KL divergence between two normal distributions are included in the supplemental material.
We mirror the approach of Cutajar et al. 2017] and Kingma and Welling 2014 by reparameterizing the weights and spectral frequencies, which allows for stochastic optimisation of the means and variances of our distributions over W and Ω via gradient descent techniques. Specifically, we use the AdamW optimizer (Loshchilov and Hutter, 2018), implemented in PyTorch, as empirically it led to superior performance compared to other alternatives such as conventional stochastic gradient descent.
4 Related Work
As mentioned previously, the work Lorenzi and Filippone 2018 also aims to incorporate physical structure into a deep probabilistic model, however the authors achieve this by constraining the dynamics of their DGP with random feature expansions, rather than specifying physics-nformmed kernel as in our approach. Additionally, the model developed by Lorenzi and Filippone| (2018) is primarily designed for tackling the problem of ODE parameter inference in situations where the form of the underlying ODEs driving the behaviour of a system are known. In contrast, our DLFM does not assume any knowledge of the differential equations governing the complex dynamical systems we study in 5) and thus we do not attempt to perform parameter inference; our primary aim is to construct robust, physics-inspired model with extrapolation capabilities and the ability to quantify uncertainty in its predictions.
The work of Mchrkanoon and Suykens 2018 is another example of deep model in which each layer consists of a feature mapping computed using random Fourier features, followed by a linear transformation. Whilst this bears similarity to the architecture of our model and to that of Lorenzi and Filippone| (0188), the authors only consider an EQ kernel for their feature mapping, with no consideration given to physically-inspire
features constraining function dynamics. /Wang et al. [219] also propose a scalable deep probabilistic model capable of modeling non-linear functions (they specifically consider multivariate time series), but their approach relies to model global behaviour, whilst relying on GPs at each layer to capture random effects and local variations specific to individual time series.
deep neural network
et al. [2021] is another example of recent work which aims to model nonstationary processes using deep probabilistic architecture, with their primary focus being spatial warping. In Duncker and Sahani (018)) the authors propose a nested GP model in which the inner GP warps the temporal input it is provided. If we consider a univariate time input, our model performs a similar process, with the weight-space GP layer warping the temporal input and passing the output into a second GP layer. However, unlike in Duncker and Sahani 2018), our framework allows for an arbitrary number of warping transformations.
We wish to clarify that although our model is based on LFMs defined by convolutions of the form used compute fa(t) in Section 2.2 model differs from the convolutional DGP construction outlined by Dunlop et al. [2018). ())) J f(0 (t- - (f+) (7)d7, which the authors argue results in trivial behaviour with increasing depth. Instead, our model takes the form, f(++1) (f(1) =
" G(F(t) ((++)))((++1) (e+1) ((+11) which more akin to the compositional construction outlined by the Dunlop et al. 20188) just with a kernel which happens to involve convolution.

---

-6
-0
[
(a) GP
(b) DGP
(c) DLFM
Figure 3: Three model fits data from toy dynamical system with the latent function u(t) cO.(0.))) + 6 6sin (3t). Grey, orange and red data-points represent training data, interpolation test data and extrapolation test data respectively. The purple curves represent the predictive means, and the shaded areas represent +20.
Experiments
All of the experimental results in this section were obtained using a single node of cluster, consisting of a 40 core Intel Xeon Gold 6138 CPU and a NVIDIA Tesla V100 SXM2 GPU with 32GB of RAM. Unless otherwise specified, all models in this section were implemented in pure PyTorch, trained using the AdamW optimizer with learning rate of 0.01 and batch size of 1000 The DLFMs and DGPs with random feature expansions [Cutajar et al. 2017 tested all utilised single hidden layer of dimensionality Dvv 3, 100 Monte Carlo samples and 100 random features per layer. Further details regarding experimental setups are provided in the supplemental material, alongside additional experimental results.
5.1 Compositional System Toy Problem
Firstly, to verify that our model is capable of modeling the compositional systems which its architecture resembles,
train the model on the vector-valued inputs t and noise-corrupted outputs (() f2(() + e of system characterised by f1(1) So Gt r)u(r)dr and fe(fi) Jo G2(f1 ='''('')d''' where G+)+) and Gz(-) represent the Green's functions cor- responding to the first order ODEs which these two integrals represent the solutions to, and ~~, 0.04). We evaluate exact GP, DLFM (with single latent force and 100 random features) and DGP with random features, with a batch size equal to the size of the training set. The latent function sinusoidal in nature, resulting in signal consisting of pulses with varying amplitude. From the model predictions shown in Figure 3 we see that all three models capture the hehaviour of the system well in the training regime, but only the deep models are able to extrapolate. The DLFM and DGP have similar normalised mean squared errors (NMSE) on the extrapolation test set of 1.5 and 1.4 respectively, however, the DLFM appears to provide more realistic quantification of predictive uncertainty, which is supported by the fact that the mean negative log likelihood (MNLL) evaluated on the extrapolation test set for the DLFM was 2.4, whereas for the DGP it was 3.7.
5.2 PhysioNet Multivariate Time Series
To assess the ability of our model to capture highly nonlinear behaviour, evaluate its performance on a subset the CHARIS dataset (ODC-BY 1.0 License) [Kim et al. 2016). which can be found on the PhysioNet data repository [Goldberger et al. 2000]. The data available for each patient consists of an electrocardiogram (ECG), alongside arterial blood pressure (ABP) and intracranial pressure (ICP) measurements; all three of these signals are sampled regular intervals and are nonlinear in nature. We use a subset of this data consisting of the first 1000 time steps for a single patient. We test two variations of the DLFM in this section, per force (Q = 2), and one with a single latent force and 100 random features (Q = 1). To clarify, both of these models contain 100 random features per layer, but in the Q 2 case, 50 features are
with latent forces and 50 random features
'Our code is publicly available in the repository: https: //gi thub. conm//ooccdooaa//DDee---e The code was also included within the supplemental material at the time of review.
7

---

PL r0-
DGP
ABP
1
05
00
25
EGG
60
10
8
Time
Time
DLFM Predictiva Mean
Training Data
Test
DGP Predictive Mean
Figure 4: Predictions generated for each of the three outputs within the CHARIS dataset from DLFM with one latent force (in the left-hand column), and a DGP with an exponentiated quadratic kernel (in the right-hand column). The shaded grey areas in each plot represent 22..
Table 1: Extrapolation test set results for each output within the CHARIS dataset, with standard error in brackets.
ABP
ECG
ICP
NMSE
MNLL
NMSE
MNLL
NMSE
MNLL
DLFM (Q= 1)
0.39 (0.05)
0.75 (0.07)
0.45 (0.04)
0.99 (0.05)
0.36 (0.04)
0.90 (0.06)
DLFM (Q- 2)
0.42 (0.07)
0.77 (0.08)
(0.98)
1.16 (0.16)
0.32 (0.06)
0.88 (0.11)
DGP-EQ
0.73 (0.002)
1.21 (0.001)
(0.01)
1.22 (0.01)
0.44 (0.01)
1.00 (0.01)
DGP-ARC
1.01 (0.001)
2.60 (0.02)
1.03 (0.001)
2.73 (0.04)
1.36 (0.003)
3.25 (0.02)
VAR-GP
1.01 (0.0002)
26.90 (9.02)
1.02 (0.001)
4.98 (0.10)
1.33 (0.003)
15.50 (6.48)
DNN
1.04 (0.01)
N/A
1.08 (0.02)
N/A
2.41 (0.09)
N/A
LFM-RFF
1.00 (0.003)
N/A
1.11 (0.02)
N/A
1.02 (0.01)
N/A
derived from two distinct latent forces. We compare these models to two DGPs with random features, employing EQ (DGP-EQ) and arccosine (DGP-ARC) kernels, as well as deep neural network (DNN) with a single hidden layer of dimensionality 300, ReLU activations, and a dropout rate of 0.5. In addition to deep models, we also evaluate the performance of a heterogencous multi-output GP trained using stochastic variational inference (VAR-GP) [Moreno-Muñoe et al. 2019] and a shallow
LFM with random Fourier features derived from the same first order ODE as our model (LFM-RFF)
Guarnizo and Álvarez 2018), consisting of two latent forces and 20 random features. Experimental results comparing the ability of each of these models to impute values within the saining-inpur space are presented in the supplemental material, however we focus here on the more challenging task of extrapolating beyond the training input-space by training the aforementioned models on the first 700 observations and withholding the remaining 300 as test set.
The results in Table [ show that our model outperforms the other techniques tested, with the single latent force variant in particular converging to a significantly lower NMSE and MNLL than the other approaches across all three outputs. From these results we can conclude that the improved performance we see is result of the ODE-based random features in the DLFM allowing the model to extrapolate accurately beyond the training input-space, rather than purely due to the additional flexibility afforded by allowing an increased number of lengthscales per layer, as the Q = 1 model contains the same number of lengthscales per layer as the DGP. We also find that the deep models tested consistently outperform the shallow models. From Figure able to accurately extrapolate the nonlinear dynamics present the system than the DGP, especially the variations in ABP.
4 we can see that the DLFM is more
8

---

NMSE
MNLI
Powerplant
101
007
n= n=9568, d=4
-0.2
003
-0.3
0.80
1.4
Protein 4=#3530, d=9
0.70
1.2
0.60
Time
Time
DLFM
DGP-EQ
DGP-ARC
VAR-GP
Figure 5: Progression of validation set metrics the UCT benchmarks, averaged over three folds.
Table 2: Test set results on the UCI benchmarks, with standard error in brackets.
Powerplant
Protein
NMSE
MNLL
NMSE
MNLL
DLFM
0.0557 (0.0005)
0.017 (0.005)
0.698 (0.004)
0.754 (0.006)
DGP-EQ
0.0581 (0.0002)
0.016 (0.003)
0.553 (0.005)
1.101 (0.002)
DGP-ARC
0,0619 (0.0002)
0.039 (0.002)
0.644 (0.027)
1.203 (0.008)
VAR-GP
0,0576 (0.0001)
0.050 (0.001)
0.520 (0.005)
1.174 (0.002)
DNN
0.0754 (0.0011)
N/A
0.664 (0.010)
N/A
5.3 UCI Regression Benchmarks
Finally, we also evaluated the performance of the model on two regression datasets from the UCI Machine Learning Repository (Dua and Graff] [017], 'Powerplant' and 'Protein', which both consist of a multivariate input and univariate target. For both of these experiments, we evaluate the same DLFM, DGP and DNN architectures as for the PhysioNet experiment, alongside multi-task variational GP (VAR-GP) [Hensman et al., 2015) implemented using GPyTorch [Gardner et al.) (018)) Achieving high performance benchmark regression datasets not the primary motivation behind this work, however from the results shown in Figure 5 and Table we find that the DLFM exhibits comparable performance to the models described above. Most notably, our model consistently converges to a superior mean negative log likelihood (MNLL), both on the validation set during training and also on the held-out test set.
Conclusion
We have presented a novel approach to modeling highly nonlinear dynamics with a sound quan- tification of uncertainty, using compositions of random features derived from first order ordinary differential equation (ODE). Our results show that our model is able to effectively capture non- linear dynamics in both single and multiple output settings, with the added benefit of competitive performance on benchmark regression datasets.
Whilst we did find that utilising kemnel based on a first order ODE yielded superior results compared to competing models, the simplicity of the ODE on which our model is based limits its flexibility. However, random features can also be derived from expressive kernels based on higher order
ODEs Guarnizo and Álvarez 2018 and even partial differential equations. Integrating such kemelss

---

into our modeling framework, and perhaps even varying the kernel choice between layers, would be interesting avenues for future research. Another limitation of our approach that whilst we give W and Q a variational treatment, we do not extend this to the kernel hyperparameters θ this could also be considered in future work.
In this work, we have demonstrated the applicability of our model to healthcare data in Section 5.2 Whilst our results demonstrate the capability of our model to learn complex dynamics within a highly nonlinear system, it is important to note that relying heavily on machine learning (ML) software within the domain of healthcare presents the possibility for negative societal impact, specifically ethical and practical issues. Care must be taken to ensure that such software is thoroughly tested and that appropriate safeguards are in place avoid situations which could cause possible harm to patients. For example, clinical decisions should be made with input from trained professional, rather than purely based upon information from an MI. system.
Acknowledgments and Disclosure of Funding
Thomas M. McDonald thanks the Department of Computer Science at the University of Sheffield for financial support. Mauricio A. Álvarez has been financed by the EPSRC Research Projects EP/R034303/1, EP/T00343X/2 and EP/V029045/1.
References
Ted Meeds, Geoffrey Roeder, Paul Grant, Andrew Phillips, and Neil Dalchau. Efficient amortised Bayesian inference for hierarchical and nonlinear dynamical systems. In International Conference on Machine Learning, pages 4445 4455. PMLR, 2019.
Sanmitra Ghosh, Paul Birrell, and Daniela De Angelis. Variational inference for nonlinear ordinary differential equations. In International Conference on Artificial Intelligence and Statistics, pages 2719-2727. PMLR, 2021.
Mauricio Alvarez, David Luengo, and Neil D Lawrence. Latent force models. In Artificial Intelligence and Statistics, pages 16. PMLR, 2009.
Wil Ward, Tom Ryder, Dennis Prangle, and Mauricio Alvarez. Black-box inference for non- non-linear latent force models. In International Conference on Artificial Intelligence and Statistics, pages 3088-3098. PMLR, 2020.
Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. nature, 521(7553):436-444 2015.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In Artificial intelligence and statistics, pages 207-215. PMLR, 2013.
Kurt Cutajar, Edwin Bonilla, Pietro Michiardi, and Maurizio Filippone. Random feature expansions for deep Gaussian processes. In International Conference on Machine Learning, pages 883. PMLR, 2017.
Siamak Mehrkanoon and Johan AK Suykens. Deep hybrid -ernnll networks using random Fourier features. Neurocomputing, 298:46-54, 2018.
Marco Lorenzi and Maurizio Filippone. Constraining the dynamics of deep probabilistic models. In International Conference on Machine Learning, pages 3227-3236. PMLR, 2018
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Representation learning: A review and new perspectives. IEEE Transactions on Pattern Analysis and Machine Intelligence, 3598):1798-1888, 2013.
Andreas Damianou. Deep Gaussian processes and variational propagation of uncertainty. PhD thesis, University of Sheffield, 2015.
Salimbeni and MP Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. In NIPS, volume 31, pages 4591-4602. Neural Information Processing Systems (NIPS), 2017.
10

---

Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Mare Deisenroth. Deep Gaussian processes with seportance-weigheee variational inference. In International Conference on Machine Learning. pages 5589-5598. PMLR, 2019.
Haibin Yu, Yizhou Chen, Zhongxiang Dai, Kian Hsiang Low, and Patrick Jaillet. Implicit posterior variational inference for deep Gaussian processes. arXiv preprint arXiv: 1 1998, 2019.
M Havasi, JM Hemmández-Lobat,, and JJ Murillo-Fuentes. Inference in deep Gaussian processes using stochastic gradient Hamiltonian Monte Carlo. Advances in Neural Information Processing Systems, 2018:7506-2 -516, 2018.
Cristian Guarnizo and Mauricio A Álvarez. Fast kernel approximations for latent force models and convolved multiple-output Gaussian processes. In Uncertainty in Artificial Intellience: Proceedings of the Thirty-Fourth Conference (2018). AUAI Press, 2018.
Carl E. Rasmussen and Christopher Williams. Gaussian Processes for Machine Learning. MIT Press, orggg/////rgggg///////m/m/
2006. Uttpttpp:/www htt.aappp//ww..gguuuiappooooooo waaaussiapprocsss.ooo.
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference.
Journal of Machine Learning Research, 14(5), 2013.
Diederik P. Kingma and Max Welling. Auto-encoding variational Bayes. In 2nd International Conference on Learning Representations, ICLR 2014, Banff. AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014. Uttp:t/:://rrat org/abs/ 1112.1112222221 6114
Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. In International Confer- ence on Learning Reprexentations, 2018.
Yuyang Wang, Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski.
Deep factors for forecasting. In International Conference on Machine Learning, pages 6607-6617. PMLR, 2019.
Andrew Zammit-Mangion, Tin Lok James Ng. Quan Vu, and Maurizio Filippone. Deep compositional spatial models. Journal of the American Statistical Association, pages 1-47, 2021.
Lea Duncker and Mancesh Sahani. Temporal aligmment and latent Gaussian process factor inference in population spike trains. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 10466-10476, 2018
Matthew M Dunlop, Mark A Girolami, Andrew M Stuart, and Aretha L Teckentrup. How deep are deep Gaussian processes? Journal of Machine Learning Research, 19(54): 2018.
Nam Kim, Alex Krasner, Colin Kosinski, Michael Wininger, Maria Qadri, Zachary Kappus, Shabbar Danish, and William Craelius. Trending autoregulatory indices during treatment for traumatic brain injury. Journal of clinical monitoring and computing, 30(6):821-831) 2016.
Ary L Goldberger, Luis AN Amaral, Leon Glass, Jeffrey M Hausdorff, Plamen Ch Ivanov, Roger G Mark, Joseph E Mietus, George B Moody, Chung Kbung-Kang Peng, and H Eugene Stanley. PhysioBank, PhysioToolkit, and PhysioNet: components new research resource for complex physiologic signals. circulation, 101(23)=e15-<222, 2000.
Moreno-Muñoz, A Artés Srdés-Rodrigue, and MA Álvarez. Heterogeneous multi-output Gaussian process prediction. In Advances in Neural Information Processing Systems, volume 31, pages 6711-6720. Neural Information Processing Systems Foundation, 2019.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL httht:://arcciie.ic.c. uci edu/ml
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pages 351-360. PMLR, 2015.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. Advances Neural Information Processing Systems, 31:7576-7586, 2018.


Reviews:
[{"review_id": "oyI7cMQ_hIb", "reviewer": "Reviewer_YdjF", "summary": "This paper introduces a deep GP for nonlinear dynamical systems. Crucially, each of the GPs uses a physics-informed kernel a la the latent force model.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "# Strengths\nLearning good models of nonlinear dynamical systems is an important topic that spans many fields such as RL and state-space models to name a few. I think this is a novel combination of deep GPs + random features + latent force models that allows for fast and efficient learning and inference.\n\n# Weaknesses\nI think the biggest weakness is the structure of the paper. Firstly, I don't think the idea of using physics-informed random features was well motivated. Specifically, for a simple latent force model the functional form of the dynamics is assumed--providing a very strong inductive bias--and a GP is used to learn the force components. Or for deep CNN, the pro of using multiple convolutional layers is that they extract a hierarchy of visual features. It would have been nice if the authors would have demonstrated some empirical evidence that these physics-informed random features give the model an inductive bias towards certain types of dynamics.\n\nNext, I think there was unnecessary notation that made the paper harder to follow. Specifically, in section 2.2 the authors introduce random feature expansions for latent force models with a generic Green's functions but in section 3.1 they end up using the Green's function for a linear ODE; this led to an increase in, my opinion, unnecessary notation that just made the approach harder to follow. \n\nLastly, I am puzzled by the choice of the UCI regression benchmarks task in the experiments section. The goal of this paper was to introduce a method for learning dynamical systems. While I am impressed that the proposed approach works well here, these types of datasets isn't what the approach was motivated for (which ties back to my point about lack of motivation). For instance, there are a number of synthetic dynamical systems that are commonly used to benchmark sequential VAEs (these can be adapted to this instance by having the latent dynamics be fully observed): FitzHugh-Nagumo, Lorenz attractor, Cart-pole, different robots from MuJoCo, etc. I would have been more impressed seeing these types of results in a paper about learning nonlinear dynamical systems.\n\n# Conclusion\nWhile I do think this is a cool idea, I think more time needs to be spent on the writing of the paper and applying the proposed approach on common nonlinear dynamical systems.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}, {"review_id": "fwW17Uoiyh1", "reviewer": "Reviewer_aNaL", "summary": "This paper introduces a method for incorporating physical dynamics into Deep GPs for the purpose of modeling nonlinear dynamical systems. The dynamics are incorporated through a kernel built via the Green's function for an associated dynamical system. ", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, the paper was well motivated and well written.  My main concerns are have to do with related work:\n\nFirst, I would have appreciated additional summary of how this method compares to previous work. The introduction and related work section discuss a number of methods and techniques that this work builds off of or is similar to (LFMs, DGPs, random Fourier features, etc) as well as different axes along which the related works differ (shallow vs deep, if or how they incorporate physical dynamics, etc). Perhaps include a table that enumerates these related works and which combination(s) of the techniques they studied?\n\nSecond, the experimental section could make additional comparisons to other methods for modeling nonlinear dynamical systems. The paper mentions the work of Lorenzi and Filippone as being quite similar, but having a different approach to incorporating dynamical constraints. Shouldn't their method also be included in the experimental evaluation?\n\nFinally, the paper could provide more guidance on how to choose the number of latent forces and the number of random features, in addition to other hyperparameters. How do these choices affect the results?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "eBMvFevRRrs", "reviewer": "Reviewer_Q3vp", "summary": "This manuscript describes Deep Latent Force Models (DLFM), a nonlinear dynamical system with nicely calibrated uncertainty quantification. Similar to deep Gaussian processes (DGPs), the presented framework stacks layers of GP layers with physics-inspired covariance functions following ordinary differential equations (ODEs). The inference is based on the weight space approximation of GPs that maintains variational distributions over the weights and spectral frequencies. The model outperforms (deep) GPs on time series prediction tasks and performs slightly better in regression tasks.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This submission proposes a very natural extension to LFMs by \"making it deeper\". Although the idea and motivation are straightforward, the authors have done a great job at describing, formulating, and evaluating the method. The experiments are not very exhaustive but adequate to show clear improvements over both DGPs and LFMs. Also, the manuscript is written very clearly; I had no difficulty following the background and the methodology. All in all, the paper does not lack to address any major concerns of mine and I recommend an accept. \n\nBelow are more detailed comments and a few questions. I would be happy if these are addressed in the rebuttal (and possibly in the paper).\n\n- The authors did not specify a recipe for choosing the form of the differential function. Since it directly impacts the covariance, gaining some intuition would be valuable. Also, is the same differential equation (line 132) used in all experiments?\n-  Without any latent layer, the input to the function $f_d(t)$ was time and the ODE helped including the physics into covariance function, both of which make sense for learning dynamical systems. Now in this deep model, how can we interpret the differential equation in the latent space? \n- Why is the VAR-FIXED trick needed? Also, are the same Monte Carlo samples used both for training and testing? If so, does not it imply potential overfitting to these samples?\n- I wonder if the model has passed a stress test. How would deeper architectures, wider hidden layer(s), more/fewer random features, and so on affect the results, training speed/ease, etc?\n- Why was LFM-RFF trained with 20 random features instead of 100 as in the deeper variant?\n- Looking at Figure3b-3c and Figure5, the model seems to achieve similar MSEs as the competing methods but certainly better at uncertainty quantification. Would the authors agree with that statement? If so, stating this explicitly and a discussion on the reason(s) would help the reader better comprehend the methodology.\n- In the top-left panel of Figure4, the uncertainty drops at test time (particularly when ABP<-2). Can the authors clarify why might this be happening?\n- A note on the computational complexity would be nice.\n- This is more of a comment than a question: Instead of this (potentially overly) complicated model, what if we had a simpler two-layer architecture that first transforms the time input non-linearly and then the final output comes from an LFM? Why should the proposed formulation be preferred over this simpler architecture (or something similar)? \n\n**Post-rebuttal comment:** Thanks the authors for the rebuttal. I didn't have any major concerns, so I keep my rating as it is. Nice piece of work!", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Both are addressed in the text.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "3ypz_D86_eT", "reviewer": "Reviewer_L4Ej", "summary": "This paper introduces a model aimed to learn time series stemming from Dynamical Systems while incorporating uncertainty. The model is a Deep Gaussian Process where the random features expansion is done using a physics informed approximation derived from the Green Kernels of certain families of ODEs.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The model presented here is novel and its derivation seems principled enough in my opinion. The results show that DLFM compares favourably to DGPs both in terms of accuracy and uncertainty calibration.\n\nHowever, I am quite skeptical about the baselines used here: how is a NN with a single hidden layer deep? There are many other architectures (deeper, residual, recurrent, Neural ODE-like models,...) which can be compared in terms of accuracy and sometimes uncertainty calibration to the proposed model. More specifically, I think that the model presented in https://openreview.net/forum?id=Skey4eBYPS is particularly relevant to this work and should be presented as a baseline.\n\nThus, while I think the ideas proposed here are promising indeed, I don't think that the baselines compared to are strong enough to fully validate them.\n\nEdit after rebuttal: After reading other reviews and the answer of the authors, I am happy to raise my score: This paper brings a novel method to GP-based estimation of non-linear DS which should be of interest to many in the community.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "I don't think there are any other than what is already acknowledged in the paper. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: There was initially some spread in the reviewer scores and perception of this paper, but the reviewer consensus did lean towards acceptance in the end. The additions proposed in the rebuttal are quite extensive but written out in sufficient detail that the reviewers also believe the paper can be accepted. Please, carefully pay attention to the reviewer comments and make sure to include the promised changes in the camera-ready version. In terms of additional experiments, including the results you presented during the rebuttal/discussion phase is sufficient.

Decision: Accept (Poster)