Compositional Modeling of Nonlinear Dynamical Systems with ODE-based Random Features
Year: 2021

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.

heading: 1 Introduction
text: 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 of 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 a 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 to a model which is able to readily extrapolate beyond the training input space, thereby retaining one of 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., 2020], 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 as a 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] are 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 [LeCun et al., 2015]. Whilst a number of recent works have incorporated random Fourier features into deep models [Cutajar et al., 2017, Mehrkanoon and Suykens, 2018], the only work we are aware of which does so in the context of a physics-informed DGP is that of Lorenzi and Filippone [2018]. The authors impose physical structure on a DGP with
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
random features by constraining the function values within the model as a means of performing ODE parameter inference. Rather than following the approach of Lorenzi and Filippone [2018] and placing constraints on function values whilst using an exponentiated quadratic (EQ) kernel, we instead integrate a physics-informed prior into the structure of the DGP by utilising a kernel based on an ODE.
Our main contribution in this paper is the introduction of a novel approach to incorporating physical structure into a deep probabilistic model, whilst providing a 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 function associated with a first order ODE. These features are then incorporated into each layer of a DGP, as shown in Figure 1. To ensure the scalability of our model to large datasets, stochastic variational inference is employed as a method for approximate Bayesian 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.

heading: 2 Background
text: 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 to compute those covariances using random Fourier features.

heading: 2.1 Deep Gaussian Processes with Random Feature Expansions
text: 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 creation being the ability of deep networks to reuse features and allow for higher levels of abstraction in said features [Bengio et al., 2013], which results in such models having more representational power and flexibility than shallow models such as GPs. DGPs are effectively a composite function, where the input to the model is transformed to the output by being passed through a series of latent mappings (i.e. multivariate GPs). If we consider a supervised learning problem with inputs denoted by X = {xn}Nn=1 and targets denoted by y = {yn}Nn=1, we can write the analytic form of the marginal likelihood for a DGP with Nh hidden layers as p(y|X, θ) = ∫ p(y|F(Nh))p(F(Nh)|F(Nh−1), θ(Nh−1)) ... p(F(1)|X, θ(0))dF(Nh) ... dF(1), where
F(`) and θ(`) represent the latent values and covariance parameters respectively at the `-th layer, where ` = 0, ..., Nh. However, due to the need to propagate densities through non-linear GP covariance functions within the model, this integral is intractable [Damianou, 2015]. 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. [2017] outline an alternative approach to tackling this problem, which involves replacing the GPs present at each layer of the network with their two layer weight-space approximation, forming a Bayesian neural network which acts as an approximation to a DGP. If the `-th layer of such a DGP (consisting of zero-mean GPs with exponentiated quadratic kernels) receives an input F(`) (where F(0) = X), the random features for this layer are denoted by Φ(`) ∈ RN×2N(`)RF (where N (`)RF denotes the number of random features used), and are given by,
Φ(`) =
√ (σ2)(`)
N (`) RF
[ cos(F(`)Ω(`)), sin(F(`)Ω(`)) ] , (1)
where (σ2)(`) is the marginal variance kernel hyperparameter, N (`)RF is the number of random features used and Ω(`) ∈ RDF (`)×N (`) RF is the matrix of spectral frequencies used to determine the random features. This matrix is assigned a prior p(Ω(`)d ) = N (Ω (`) d | 0, ( l(`) )−2
) where l(`) is the lengthscale kernel hyperparameter, DF (`) is the number of GPs within the layer, and d = 1, ..., DF (`) . The random features then undergo the linear transformation, F(`+1) = Φ(`)W(`), where W(`) ∈ R2N (`) RF×DF (`+1) is a weight matrix with each column assigned a standard normal prior. Training is achieved via stochastic variational inference, which involves establishing a 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.

heading: 2.2 Random Feature Approximations for Latent Force Models
text: In this paper, we take a similar approach to Cutajar et al. [2017], but we utilise physically-inspired random features. Physics-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 Álvarez [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-inspired 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 byQ latent forces. Given an input data-point t and a set ofD output variables {fd(t)}Dd=1, a LFM expresses each output as fd(t) = ∑Q q=1 Sd,q ∫ t 0 Gd(t− τ)uq(τ)dτ , where Gd(·) represents the Green’s function for a certain form of linear differential equation, uq(t) ∼ GP(0, kq(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 kfd,fd′ (t, t ′) = ∑Q q=1 Sd,qSd′,q ∫ t 0 Gd(t− τ) ∫ t′ 0 Gd′(t
′ − τ ′)kq(τ, τ ′)dτ ′dτ . Typically, an EQ form is assumed for the kernel governing the latent forces, kq(·). Due to the linear nature of the convolution used to compute fd(t), the outputs can also be described by a GP. Exact inference in LFMs is tractable, however as with GPs, it scales with O(N3) complexity [Rasmussen and Williams, 2006]. Furthermore, computing the double integration for kfd,fd′ (t, t
′) leads to terms that include the numerical approximation of complex functions which burden the computation of the kernel function. By providing a random Fourier feature representation for the EQ kernel kq(t, t′), Guarnizo and Álvarez [2018] were able to reduce this cubic dependence on the number of data-points to a linear dependence. This representation arises from Bochner’s theorem, which states,
kq(τ, τ ′) = e
−(τ−τ′)2
`2q = ∫ p(ω)ej(τ−τ
′)ωdω ≈ 1 NRF
NRF∑
s=1
ejωsτe−jωsτ ′ , (2)
where `q is the lengthscale of the EQ kernel, NRF is the number of random features used in the approximation and ωs ∼ p(ω) = N (ω|0, 2/`2q). Substituting this form of the EQ kernel into kfd,fd′ (t, t ′) leads to a fast approximation to the LFM covariance kfd,fd′ (t, t ′) ≈ ∑Q q=1 Sd,qSd′,q NRF [∑NRF s=1 φd(t, θd, ωs)φ ∗ d′(t ′, θd′ , ωs) ] , where
φd(t, θd, ω) =
∫ t
0
Gd(t− τ)ejωτdτ, (3)
with θd representing the Green’s function parameters. When the Green’s function is a real function, φ∗d′(t ′, θd′ , ω) = φd′(t ′, θd′ ,−ω). Guarnizo and Álvarez [2018] refer to φd(t, θd, ω) as random Fourier response features (RFRFs).

heading: 3 Deep Latent Force Models
text: 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.

heading: 3.1 Model Formulation
text: 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 [2018] derived a 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)dt + γf(t) = ∑Q q=1 Squq(t), where γ is a decay parameter associated with the ODE and Sq is a sensitivity parameter which weights the effect of each latent force. For simplicity of exposition, we assume D = 1. The Green’s function associated to this ODE has the form G(x) = e−γx. By using equation (3), the RFRFs associated with this ODE follow as
φ(t, γ, ωs) = ejωst − e−γt γ + jωs , (4)
where compared to the expression (3), the parameter θ of the Green’s function corresponds to the decay parameter γ. From here on, we redefine the spectral frequencies as ωq,s to emphasise the fact that the values sampled are dependent on the lengthscale of the latent force by way of the prior, ωq,s ∼ N (ω|0, 2/`2q). We can collect all of the random features corresponding to the q-th latent force into a single vector, φcq(t, γ,ωq) = √ S2q/NRF [φ(t, γ, ωq,1), · · · , φ(t, γ, ωq,NRF )]> ∈ CNRF×1, where
ωq = {ωq,s}NRFs=1 , C refers to the complex plane and the super index c inφcq(·) makes explicit that this vector contains complex-valued numbers. By including the random features corresponding to all Q latent forces within the model, we obtain φc(t, γ,ω) = [(φc1(t, γ,ω1) > , · · · , ( φcQ(t, γ,ωQ) )> ]> ∈ CQNRF×1, with ω = {ωq}Qq=1. These random features will be denoted as φcLFM (t, γ,ω) to differentiate them from the features φ(·) computed from a generic EQ kernel.
Higher-dimensional inputs Although the expression for φcLFM (t, γ,ω) 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-dimensional vector x = [x1, x2, · · · , xp]> ∈ Rp×1. As will be noticed later, such an extension is also necessary if we attempt to use such features at intermediate layers of the composition. Essentially, we compute a vector φcLFM (xm, γm,ωm) for each input dimension xm leading to a set of vectors {φcLFM (x1, γ1,ω1), · · · ,φcLFM (xp, γp,ωp)}. Notice that the samples ωm can also be different per input dimension, xm. Although there are 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 as φcLFM (x,γ,Ω) = ∑p m=1 φ c LFM (xm, γm,ωm), where Ω = {ωm}pm=1 and γ = {γm}pm=1. An alternative to explore for future work involves expressing φcLFM (x,γ,Ω) as φ c LFM (x,γ,Ω) = ∑p m=1 αmφ c LFM (xm, γm,ωm), with αm ∈ R a parameter that weights the contribution of each input feature differently. Although we allow each input dimension to have a different decay parameter γm in the experiments in Section 5, for ease of notation we will assume that γ1 = γ2 = · · · = γp = γ. For simplicity, we write φcLFM (x, γ,Ω). Therefore, φcLFM (x, γ,Ω) is a vector-valued function that maps from Rp×1 to CQNRF×1.
Real version of the RFRFs Rather than working with the complex-value random features φcLFM (x, γ,Ω), we can work with their real-valued counterpart by using φLFM (x, γ,Ω) =
[(Re{φcLFM (x, γ,Ω)}) > , (Im{φcLFM (x, γ,Ω)}) > ]> ∈ R2QNRF×1 Guarnizo and Álvarez [2018], where Re(a) and Im(a) take the real component and imaginary component of a, respectively. For an input matrix X = [x1, · · · ,xN ]> ∈ RN×p, the matrix ΦLFM (X, γ,Ω) = [φLFM (x1, γ,Ω), · · · ,φLFM (xN , γ,Ω)]> ∈ RN×2QNRF .
Composition of RFRFs We now use ΦLFM (X, γ,Ω) as a building block of a layered architecture of RFRFs. We write Φ(`)LFM (F
(`), γ(`),Ω(`)) ∈ RN×2Q(`)N(`)RF and follow a similar construction to the one described in Section 2.1, where F(`+1) = Φ(`)LFM (F
(`), γ(`),Ω(`))W(`), and W(`) ∈ R2Q (`)N (`) RF×DF (`+1) . As before, F(0) = X. Figure 2 is an example of this compositional architecture of RFRFs, which we refer to as 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 to computing F(L)d = Φ (L−1) LFM (F (`), γ (L−1) d ,Ω (L−1) d )W (L−1) d , where d = 1, ..., D.
Model interpretation In a 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 our model can be interpreted as a filter, which selectively filters out some of 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.

heading: 3.2 Variational Inference
text: As previously mentioned, exact Bayesian inference is intractable for models such as ours, therefore in order to train the DLFM 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 of the model. We seek to optimise the variational distributions over Ω and W whilst also optimising Θ, 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 Ω(`)ij such that Ω (`) ij = s (`) ij (`) rij +m (`) ij , where m (`) ij and (s 2) (`) ij represent the means and variances associated with
the variational distribution over Ω(`)ij , and ensuring that the standard normal samples (`) rij are fixed throughout training rather than being resampled at each iteration. If we denote Ψ = {W,Ω} and consider training inputs X ∈ RN×Din and outputs y ∈ RN×Dout , we can derive a tractable lower bound on the marginal likelihood using Jensen’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 Cutajar et al. [2017], takes the form,
log[p(y|X,Θ] = Eq(Ψ) log[p(y|X,Ψ,Θ)]− DKL[q(Ψ)||p(Ψ|Θ)] (5)
≈ [ N
M
∑
k∈IM
1
NMC
NMC∑
r=1
log[p(yk|xk, Ψ̃r,Θ)] ] − DKL[q(Ψ)||p(Ψ|Θ)], (6)
where DKL denotes the KL divergence, Ψ̃r ∼ q(Ψ), the minibatch input space is denoted by IM and NMC is the number of Monte Carlo samples used to estimate Eq(Ψ). q(Ψ) and p(Ψ|Θ) 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 overW 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.

heading: 4 Related Work
text: As mentioned previously, the work of 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 a physics-informed 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 Section 5, and thus we do not attempt to perform parameter inference; our primary aim is to construct a robust, physics-inspired model with extrapolation capabilities and the ability to quantify uncertainty in its predictions.
The work of Mehrkanoon and Suykens [2018] is another example of a deep model in which each layer consists of a feature mapping computed using random Fourier features, followed by a linear transformation. Whilst this bears a similarity to the architecture of our model and to that of Lorenzi and Filippone [2018], the authors only consider an EQ kernel for their feature mapping, with no consideration given to physically-inspired features or constraining function dynamics. Wang et al. [2019] also propose a scalable deep probabilistic model capable of modeling non-linear functions (they specifically consider multivariate time series), but their approach relies on a deep neural network to model global behaviour, whilst relying on GPs at each layer to capture random effects and local variations specific to individual time series. Zammit-Mangion et al. [2021] is another example of recent work which aims to model nonstationary processes using a deep probabilistic architecture, with their primary focus being spatial warping. In Duncker and Sahani [2018], 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 to compute fd(t) in Section 2.2, our model differs from the convolutional DGP construction outlined by Dunlop et al. [2018], f (`+1)(t) = ∫ f (`)(t − τ)u(`+1)(τ)dτ , which the authors argue
results in trivial behaviour with increasing depth. Instead, our model takes the form, f (`+1)(f (`)) = ∫ f(`) 0
G(f (`) − τ (`+1))u(`+1)(τ (`+1))dτ (`+1), which is more akin to the compositional construction outlined by the Dunlop et al. [2018], just with a kernel which happens to involve a convolution.

heading: 5 Experiments
text: All of the experimental results in this section were obtained using a single node of a 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 a learning rate of 0.01 and a batch size of 1000.1 The DLFMs and DGPs with random feature expansions [Cutajar et al., 2017] tested all utilised a single hidden layer of dimensionality DF (`) = 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.

heading: 5.1 Compositional System Toy Problem
text: Firstly, to verify that our model is capable of modeling the compositional systems which its architecture resembles, we train the model on the vector-valued inputs t and noise-corrupted outputs y2(t) = f2(t) + of a system characterised by f1(t) = ∫ t 0 G1(t − τ)u(τ)dτ and
f2(f1) = ∫ f1 0 G2(f1 − τ ′)u(τ ′)dτ ′, where G1(·) and G2(·) represent the Green’s functions corresponding to the first order ODEs which these two integrals represent the solutions to, and ∼ N (0, 0.04). We evaluate an exact GP, a DLFM (with a single latent force and 100 random features) and a DGP with random features, with a batch size equal to the size of the training set. The latent function is sinusoidal in nature, resulting in a signal consisting of pulses with varying amplitude. From the model predictions shown in Figure 3, we see that all three models capture the behaviour 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 a 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.

heading: 5.2 PhysioNet Multivariate Time Series
text: To assess the ability of our model to capture highly nonlinear behaviour, we evaluate its performance on a subset of 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 at 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, one with two latent forces and 50 random features 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
1Our code is publicly available in the repository: https://github.com/tomcdonald/Deep-LFM. The code was also included within the supplemental material at the time of review.
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 a 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 heterogeneous multi-output GP trained using stochastic variational inference (VAR-GP) [Moreno-Muñoz 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 training-input 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 a test set.
The results in Table 1 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 a 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 4, we can see that the DLFM is more able to accurately extrapolate the nonlinear dynamics present in the system than the DGP, especially the variations in ABP.

heading: 5.3 UCI Regression Benchmarks
text: Finally, we also evaluated the performance of the model on two regression datasets from the UCI Machine Learning Repository [Dua and Graff, 2017], ’Powerplant’ and ’Protein’, which both consist of a multivariate input and a univariate target. For both of these experiments, we evaluate the same DLFM, DGP and DNN architectures as for the PhysioNet experiment, alongside a multi-task variational GP (VAR-GP) [Hensman et al., 2015] implemented using GPyTorch [Gardner et al., 2018]. Achieving high performance on benchmark regression datasets was not the primary motivation behind this work, however from the results shown in Figure 5 and Table 2, 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.

heading: 6 Conclusion
text: We have presented a novel approach to modeling highly nonlinear dynamics with a sound quantification of uncertainty, using compositions of random features derived from a first order ordinary differential equation (ODE). Our results show that our model is able to effectively capture nonlinear 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 a kernel 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 more expressive kernels based on higher order ODEs [Guarnizo and Álvarez, 2018] and even partial differential equations. Integrating such kernels
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 is that whilst we giveW and Ω 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 to avoid situations which could cause possible harm to patients. For example, clinical decisions should be made with input from a trained professional, rather than purely based upon information from an ML 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.

Review summaries:
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.