Constraining Gaussian Processes to Systems of Linear Ordinary Differential Equations
Authors: Andreas Besginow
Markus Lange-Hegermann
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian Process
Ordinary Differential Equations
Machine Learning
Probabilistic Model

Abstract: Data in many applications follows systems of Ordinary Differential Equations (ODEs).This paper presents a novel algorithmic and symbolic construction for covariance functions of Gaussian Processes (GPs) with realizations strictly following a system of linear homogeneous ODEs with constant coefficients, which we call LODE-GPs. Introducing this strong inductive bias into a GP improves modelling of such data. Using smith normal form algorithms, a symbolic technique, we overcome two current restrictions in the state of the art: (1) the need for certain uniqueness conditions in the set of solutions, typically assumed in classical ODE solvers and their probabilistic counterparts, and (2) the restriction to controllable systems, typically assumed when encoding differential equations in covariance functions. We show the effectiveness of LODE-GPs in a number of experiments, for example learning physically interpretable parameters by maximizing the likelihood.

Constraining Gaussian Processes to Systems of Linear Ordinary Differential Equations
Andreas Besginow
Department of Electrical Engineering and Computer Science
OWL University of Applied Sciences and Arts
Lemgo, Germany
Institute industrial IT Lemgo, Germany
andreas bes.inowtth-ooll.coo
Markus Hannee-Hgerrann
Department of Electrical Engineering and Computer Science
OWL University of Applied Sciences and Arts
Lemgo, Germany
Institute industrial IT Lemgo, Germany
markus.. hangg--hggemmmn@tthoorr
Abstract
Data in many applications follows systems of Ordinary Differential Equations (ODEs). This paper presents novel algorithmic and symbolic construction for covariance functions Gaussian Processes (GPs) with realizations strictly fol- lowing system of linear homogeneous ODEs with constant coefficients, which we call LÖDE-GPs. Introducing this strong inductive bias into a GP improves modelling of such data. Using smith normal form algorithms, symbolic tech- nique, we overcome two current restrictions in the state of the art: (1) the need for certain uniqueness conditions in the set of solutions, typically assumed in clas- sical ODE solvers and their probabilistic counterparts, and (2) the restriction to controllable systems, typically assumed when encoding differential equations in covariance functions. We show the effectiveness of LODE-GP a number experiments, for example learning physically interpretable parameters by maxi- mizing the likelihood.
Introduction
Many real world tasks have underlying dynamic behavior, for example chemical reactions [2S]]
systems in bioprocess engineering [RZ]], population dynamics [Z0]] One currently widely debated
application are compartmental models in epidemiology [20, BS 52 Many such systems are linear or can be decently linearized, such in control theory [dE],, biology [IG]], process engineering [I], $9], or engines [7]] Including prior knowledge in the form of differential equations benefits the model fit and enhances interpretability for model. Hence, modelling differential equations in Machine Learning (ML) has therefore been the focus of much research in both Gaussian Processes (GPs) (e.g.
B88.
30
and Deep Learning (DL) (e.g.
The class of probabilistic Ordinary Differential Equation (ODE) solvers allow to estimate the solu- tion of ODE initial value problems through approximations, often based e.g. Runge- Kutta meth- ods or Kalman filters [5 5z 53, [T), and thereby do not strictly guarantee to yield solutions of the ODE. This class of algorithms is commonly used to solve non-linear ODEs, but typically require that systems are well well-posed, in particular the solution needs to be unique once a finite number of
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

initial conditions is known. While this second limitation is typically irrelevant in physical or biolog- ical systems, it strongly relevant in systems in engineering, such as in control systems with their freely choosable inputs.
With early works like [zZ]], the introduction of derivatives as linear operators in GPs became promi- nent for partial differential equations. This inspired several works that encode differential equations in the covariance structure of GPs [10
X Mil to introduce a strong inductive bias into GPs. This motivated [ل to make these approaches algorithmic, build a mathematical foundation, and show the hidden assumption in previous works of being limited to controllable systems.
This paper overcomes the necessity of approximations, the restriction to systems where initial con- ditions lead to unique solutions, and the restriction to controllable systems. For that purpose, algorithmically construct LODE-GPs (Theorem M1, a novel class of GPs whose realizations strictly follow given system of homogenuous linear ODEs with constant coefficients.
Our construction of LODE-GPs starts by describing our system of ODEs via so-called operator matrix A. Calculating its Smith Normal Form (SNF) gives us a decoupled representation of the system via a diagonal operator matrix D, at the cost of going through two invertible base change matrices U, V with the relationship U A V = D. We can easily construct multi-output GP g (0, k) (cf. Lemma W with as many outputs as the system channels and whose covariance function encodes the decoupled system. Then, applying the base change matrix V to g yields the pushforward GP Vig = GP(0, VkV') whose covariance function encodes the original system, i.e. the LODE-GP. We refer to Theorem @ and its algorithmic proof for details.
Our paper makes the following contributions:
1. It develops and proves an algorithmic and symbolic construction of LODE-GPs, GPs with the strong inductive bias that their realizations strictly satisfy any given system of homoge- neous linear ODEs with constant coefficients (see Theorem M).
2. It demonstrates that the constructed GPs numerically stable and robust (see Section B
3. It automatically includes ODE system parameters as GP parameters, learns them during the training process of the GP, and allows interpretation of the data (see Subsection 5)2).
4. It provides a free and public implementation based on the GP library GPyTorch [22]]
We successfully test our approach on controllable and non-controllable systems of differential equa- tions, with and without system parameters. The LODE-GP strictly satisfies the ODEs and outper- forms GPs by several magnitudes in its precision, additionally it correctly reconstructs the system parameters used generate the data with small relative error, even despite noise. We also discuss how this strict behavior influences results for metrics like the RMSE with e.g. Figure !!
2 Related Work
Using physical information in MI. systems became central research topic over recent years in both DL and GPs. With the GP research in two important categories: (1) Works that introduce derivative information as inductive bias into the GP
A |il and (2) works that train GPs without inductive bias on physical data, with specially selected standard kernels, as an approximate H. 122
[RO.
BR
Ry
5
Of all the works GPs, [B8, R9] are the closest to us, as they also provide general algorithms to introduce inductive bias into GPs, based on Gröbner bases. These approaches only applicable controllable systems. By basing our algorithm on the SNF instead of Gröbner bases, we overcome this limitation to controllable systems, but are restricted to ODEs.
model for additional processing steps [6 R.
[T].
Probabilistic ODE solvers like [52 of ODEs but sometimes require a number of comparably complex calculations and some approx- imations. We avoid any approximation through calculation the SNF and the application of the pushforward with linear differential operator. Hence, work is limited to linear ODEs, whereas probabilistic ODE solvers are able to work with non-linear ODEs.
are able to include inductive bias
For GP priors for decoupled ODE systems with constant coefficients and right hand side functions see [R]] which consideres the right hand side functions as latent, hence these models are called Latent

---

Force Model (LFM). A GP prior on these latent forces and is pushed forward through differential operators and Green's 's operator. We empirically compare LODE-GPs LFMs in the appendix.
Of the other works that introduce inductive bias, [ and [-TI]] are the earliest works and already discuss general formulation of differential equations by viewing derivatives linear operators that can be applied to GPs. Subsequent works target specific ODEs and create covariance functions with inductive bias T [ They can be considered special cases of our approach with specific manually constructed pushforward. Finally, others constructed covariance functions for partial differential equations from standard covariance functions with e.g. curl-free behaviour [69]] 1
to model physical behaviour 669
Y
X or use the GP as a surrogate model to learn physical or a variety of applications [+].
problems like the inverse pole problem
Most DL. models explore physically informed models through additional loss terms that punish solutions that deviate from the system [9
SW. [K], in the case of [[С]] show the equivalence of
the loss to their assumed physical behavior. This is often combined with different modifications to
the networks like additional features [TM]] specific network structures [4%
3 Background
3.1 Gaussian Processes
A GP GP(u, k) defines probability distribution over the space of functions R2 R such that the outputs g(z;) at any set of E R are jointly Gaussian [72]] Such (multi-output) GP is defined by its mean function (often set to zero)
p R
R x B(g(z))
and its (multi-outputt) positive semi-definite covariance function (also called kernel)
12d x Rd
R (I,x') (((((())~(=))(00))) - -u(x))²²)
- m(x))(((xx))
popular kernel is the Squared Exponential (SE) kernel ka x') 0² exp
(- [ 22
with its
signal variance and lengthscale f. It models smooth data very well and is thus usable for many datasets [|Z p. 83), as its realizations are dense in the space of smooth, i.e. infinitely differentiable, functions th (TR, RR) 20, Prop. 1].
Manipulating existing GPs allows introduce indutive biases in various ways
One important example is the application of matrices of linear operators to GP. Formally, this is the pushforward operation of an operator matrix on the GP g as B.g = GP (p(x), Bk(r, (B')) [I. where B' denotes the operation of B on the second argument of k(r, x²) [B0 Lemma 2.2]. These operators can, for example, be differential operators. The matrix B induces the strong bias such that all realizations lie in the image of B. This pushforward is typical for applications in differential equations [B0,, 3g. pR. bll or geometry [[8].
Example GPs can be constrained realizations satisfying system of linear equations given by matrix A, e.g. A = [2 -3]
solz(A) = (1A(w)) fa(x)) e FeX1 A (f((()) f2(x)) =0-
33
with J = 0000 (R,, R.) the space of smooth input functions. The matrix B
is maximal in that
2
it solves A B- = 0 and hence solz(A) F {B [(x) [xxx e FJ. Taking a GP prior
g = GP(, k) for (xx E F and applying the pushforward then yields a new, constrained GP prior
B.g GD BK(B)) gP ((0)
9 ksF KSF
6 ksF kSF
for B f(x) (A), whose realizations guaranteed lie in sol
This example is similar to our use case of systems of differential equations for LODE-GPs, where we just replace the matrix B of numbers with suitable matrix of differential operators.
3

---

(1)
1
82
f1(1)
(t
Figure 1: A visualization of the bipendulum with its components (rope lengths 4,,2) and states
(angles (1),, [2((), rod acceleration x" (())).
3.2 Smith normal form
The SNF
is a normal form for a matrix A E R[x])xx over the polynomial ring R(x), such that U A V D. Here, D E [[x]mx is (not necessarily square) diagonal matrix of same size as and the base change matrices U E [x and E [[æ]"*" are invertible square matrices, i.e. det(U), det(V)) det(()) E IR (0) EE IE Algorithms to construct the SNF are implemented
in many computer algebra systems such as Matlab [H3], Maple [[]]], or SageMath [50] (based on
PARI [69]] which is a free and open source python library for computer algebra. Intuitively, the SNF can be compared to the eigendecomposition a IR x matrix, which if exists producess a square matrix of eigenvectors W and the diagonal matrix of eigenvalues A. This analogy is lacking since e.g. the SNF produces two independent base change operator matrices V resp. U instead of single eigenvector matrix W resp. its inverse W-1 Still, the matrix D decouples and thereby simplifies system given by the input matrix A, as does the diagonal eigenvalue matrix A [S0,, [e]]
Since the SNF exists for matrices over polynomial rings over any field, can also compute it over a polynomial ring over the function field ;.,,,..... (x). Hence, we can model differential equations The heating example in Subsection b demonstrates how such
containing parameters di,..., dg. parameters are being used in the SNF and subsequently learned from data.
Example 2. bipendulum acts as our running example. It consists of rod with pendulum attached to each end of the rod, see Figure E The linearized equations of an idealized bipendulum:
ft)
I²(t) + (.)(()) + 9 x²(1) (2f(() +
0
81
fe(t)
(1)
0
82
u(tt)
with the operator matrix A, L 0 the length of the i-th pendulum, g 9.81 the gravitational constant, u(t) -x" (t) the rod's acceleration, and fi(t). = 1,2 the angle between the i-th pendulum and the vertical axis. The SNF of reveals its properties, specifically when & = lg and L * Eg. An algorithm calculates the SNF D of A and the base change matrices U, V as follows:
case E = 1 + lo=
4
20 tg
1
0
0
+ 9
0
1
0
0
0
0
1
2
48
+
4g
(0)
(
g)
g
case
1/2 1:
0
07 0
0
1
0 0
0
1
1
-
07
+
06
+
0
07
In case of unequal rope length, the controllability intuitive when imagining that, by accelerating the rod at the right time, the two ropes can reach every combination of valid positions. Formally,

---

this follows that D (ht(t) h2(t) ha(t)]² = 0 if and only if 1u (4) he(() = 0 and h((()) is an arbitrarily choosable output. In case of equal rope length, the system not controllable, indicated by the diagonal entry 0g + 9 of D PE p. 154].
4 Constructing a GP for differential equations
We introduce a polynomial time algorithm to construct LODE-GPs, a class of GPs with realizations dense in the space of solutions of linear ODEs with constant coefficients. This ensures that the LODE-GP is able to produce all possible solutions to the ODEs and nothing but solutions for the ODEs. All this guaranteed be strictly accurate, since we make no approximations.
Consider system of linear homogenous ordinary differential equations with constant coefficients
A f(t) =
(2)
with operator matrix A E R[O,]mx determining the relationship between the smooth functions fi(t) E C (I, R) of E(t) = (fi(t)
fn((t)) For such systems our main result holds.
Theorem 1. (LODE-GPS For every system as in Equation (D) there exists a GP such that the set of realizations of g dense in the of solutions of A f(t) Af((t) = O.
The following lemma play crucial role in the proof of this theorem by constructing base covariance functions for system as in Equation (l) that decoupled (via the SNF) scalar equations. The proof in Appendix D makes use of the solution space being decomposed, in particular it deals with finite dimensional vector spaces like Bayesian linear regression.
Lemma 1. Covariance functions for solutions of the scalar linear differential equations d f with constant coefficients, i.e. d E [[ are given by Table for d is primary, i.e. a power of an irreducible real polynomial. In the case of non-primary d, each primary factor d, of d = II d, first translated to its respective covariance function kú separately before they are added up to give the full covariance function k(l1, = -'' ki (1(1112))
Table 1: Primary operators d and their corresponding covariance function k(t1, t2).
d
k(t1,tz)
(0: aa ((0: a- + 2b))
E fits expla (t1 + t2))
E 1- tits) exp(a (t1 tz)) cos(b (t1 t2)) t2)
For the case of single (i.e. j 1) zeroes, the sum in the covariance function simplifies a factor 1. Example 3. The ODE d 0 has solutions of the form exp(t) for E R. This one- dimensional solution space is described by the covariance function k((1,, 12) = exp(tz) exp(t2) = exp(t1 + tz).
Proof of Theorem M
We begin with neutral multiplication of A = with V V-1 and left multiplication by to decouple the system using the SNF A V (cf. Equation (D)):
A. V - f
D. f
D D.p=0
(3)
min(n.m)
D P 0 ^
0 0. =
minn (1.m)++
for decoupled latent vector p V- of functions.
5

---

We assume the GP- GP-prior h ~ GP(O, k) for this vector p via GP, where & is a multi- diagonal covariance function such that each diagonal entry of given by Lemma By this Lemma, the set of realizations of this prior h is dense in the set of solutions of D p 0.
The pushforward of h with V yields GP g that can learn the original system of ODEs.
g V.h GP V k. V")
(4)
with V = VT the operation applied on the second entry of the kernel (i.e. t2) using [3 Lemma 2.2]. As V is invertible, the set of realizations of g is dense in the set of solutions of A- =
This proof shows that diagonalizing the system to D p = 0 decouples the components of p to expose the intrinsic behavior of the system. We can also easily interpret the controllability of a system in this decoupled form. A zero column (or 0- P, 0) stands for a freely choosable function in p. i.e. an output of the system. A one in a column (or p. = 0) stands for something we have no choice in, e.g. an input or state that is fixed once the desired output was chosen. Any other entry leads to manusiidal-expoeettia behavior, which is uncontrollable in the sense that it cannot be influenced. This behaviour is classic in systems of differential equations, as most of them can be split into a controllable and uncontrollable solution space, which can be clearly seen in the entries of D. Summing up, system is controllable iff all diagonal entries of D are either zero one.
This paper improves upon [B8] for two reasons: first, it decouples the uncontrollable part of a system via the SNF as a direct summand and, second. it constructs a covariance function for this uncontrol- lable summand. This generally impossible for partial differential equations due to [B. H6]
The runtime of the construction of LODE-GPs in Theorem M is dominated by computing the SNF. This has deterministic polynomial complexity in the size of the (usually small) operator matrices [ with quicker probabilistic Las Vegas algorithms [[8], and it can be parallelized [66, 58). For our examples, the SNF terminates instantaniously. The application of the base change matrix V is cubic in the (usually small) size of matrices, assuming constant time to apply operators to the covariance functions from Lemma M1 Once the covariance function constructed, the complexity of O for GPs applies, where n is the (potentially big) number of data points.
A [55, 57, 5G6
Example 4. We continue Example and consider the SNF of A for 2 # l2- Since the three diagonal entries of the matrix D are 1, 1, 0 conclude that latent kernel & has diagonal entries O, O, kSE with Lemma M. Now the pushforward creates the LODE-GP, which requires the following application of two operator matrices; we simplify by removing irrelevant zeroes of the latent covariance.
202 +9
V&V'
80 tg
[KSE]
28 2
+
(0P, + (00, + +9]
(0), + (OP + g)
Applying these two operators to the kernel KsE results in the covariance function for the bipendulum.
Experimental evaluation
We demonstrate the effectiveness of LODE-GPs in three examples, where constrain LODE-GP using the system description and train it with datapoints, similar to solving an initial value problem. The Gröbner basis approach of [[8] yields precisely the same results for the two controllable ex- amples and is not applicable to the uncontrollable example. The only other method that could deal with our class of differential equations is [ for which comparison with the three tank system is discussed in the appendix, in the following we compare our model mainly to classic GPs.
Our comparison includes the error in satisfying the ODEs, specified by the median error the GPs posterior mean function has in satisfying the ODEs at evenly spread points, where we calculate derivatives through finite differences. The LODE-GPs hyperparameters (lengthscale and signal vari- ance) are randomly initialized from a uniform distribution on [-3, 3] and were trained using Adam [ The 25 training datapoints are created uniformly in the intervals [1, 6] and [-5, 5] for the bipendulum and three tank system, resp. the heating system. Training and evaluation is repeated
2 Most probabilistic ODE solvers are not applicable to systems with free functions in their solution set, with the exception of [5 which can only estimate free functions of parameters.

---

Table 2: The median results and standard dev. for training and evaluation RMSE, loss value (negative marginal log likelihood per datapoint) and the mean ODE satisfaction error, over 20 experiments for the GP and LODE-GP on noisy training data. The mean ODE enror is the average of each ODE error for function in a system. The RMSEs were calculated with noiseless data in the training and evaluation interval, the ODE error only with the latter, for each respective experiment. Smaller is better.
training RMSE
evaluation RMSE
loss
mean ODE error
Bipendulum
LODE-GP
0.060 + 0.060
0.106 + 0.258
-1.638 + 0.567
5.135e-07 + 1.607e-06
GP
0.008 0.030
0.044 + 0.0175
-2.072
0.215
0.064 0.023
Heating
LODE-GP
0.006 + 0.001
0.170 + 0.068
-2.089 + 0.100
0.008 + 0.007
GP
0.010 0.001
0.333 + 0.055
-1.628 + 0.051
0.218 0.070
Three tank
LODE-GP
0.023 0.004
0.078
0.046
-0.949
0.063
1.360e-05
4.919e-06
GP
0.028 + 0.005
0.058 + 0.031
-0.974 + 0.071
0.040 + 0.020
A
Noisy
Noisy data Posiziar maso
Noisy
022
Posteriox mean
0.2
Posictior mcaan
Noise tree data. confidence
Noise free data confidence
Nouse fiee data
confilence
0.1
-0.
-1.1
-0.2
-02
00.2
0.4
6
2
Figure 2: A posterior LODE-GP trained on noisy bipendulum data (black stars), its confidence (blue transparent), posterior mean (blue dashed line) and the original noise-free function (red line). The posterior mean fits the original data very well and shows that the behavior is learned despite noise.
ten times in each experiment using a GPyTorch |2A]|] implementation of our LODE- GP construction with SageMath [60] to symbolically calculate the SNF. The training loss is the negative marginal log likelihood fc.f. re Eq. (2.30))) which, in GPyTorch, is additionally divided by the number of train- ing data points. Details of the training, data generation, verification process, additional experiments, e.g. without noise, and further system details can be found the appendix.
5.1 Bipendulum Controllable
We continue with the controllable variant of the bipendulum from Example with the rope lengths 4 and &g = 2. We create 25 points of training data from a solution of the ODEs as shown in Figure We add white noise the data with standard deviation of 2% of the maximal signal.
Figure compares the error satisfying the ODE LODE-GP and a GP. The LODE-GP is pro- ducing samples that satisfy the ODE with median error of 2c-6, for both differential equations. In comparison, the symbolically computed solution used to generate the training data, has an error of 1.75e-7, only one order of magnitude better than the LODE-GP, Hence, we conclude that, up to numerical precision. the LODE-GP produces samples that strictly satisfy the ODE, whereas re- gression models like GPs do not satisfy the differential equations. The LODE-GP produces roughly 30000 times more precise samples (0.06 to 2e-6) and deviates only slightly from this error.
Further investigate the RMSE of the models by comparing their posterior mean noiseless data in the training and evaluation intervals. The LODE-GP RMSE mostly performs similar compared to the GP RMSE (see Table D) For the training interval, the LODE-GP achieves basically the
'The LODE CPS trains into a local minimum at its value 10 variance. This models constant behaviour which satisties the ODEs, but does not fit the training data.
10 with high lengthscale and miniscule signal
7

---

GP ODE2
GP ODE1
LODE-GP ODE2
LODE-GP ODEI
10
10
10
10
10
10
10
10
10
10-3
10-2
10
100
Figure 3: The enror in satisfying the bipendulum ODEs in Equation (M) for GP (top) and LODE- GP (bottom), for noisy training data. Smaller better. The green dashed line shows the error of solution in symbolic form. The low error shows that the LODE-GP strictly satisfies the ODEs, up to numerical precision.
fi(t)
fi(t)
2
f2(1)
f2(1)
m(1)
1
x(1)
0
10
15
20
0
5
10
15
20
Figure 4: The posterior mean functions of LODE-GP (left) and GP (right) trained on the bipen- dulum data as described in the text.
same performance as the GP. For the evaluation interval the performance worsens due to SE GP extrapolating to its (constant zero) prior mean. The GP can easily reach this mean zero, by smoothly modifying all channels without regard for the system equations. Since the ODE solution also approaches zero, this leads system equations to move the mean function back to mean zero. Hence, the system produces large spike the input channel (t), before the system can reach mean zero, as seen Figure
small error extrapolation. The LODE-GP has to satisfy the
We also analyse the training runtime and see that the LODE-GP is as fast as the GP during the (nn model inference and slower by factor smaller than two in the ((n²) covariance matrix calculation.
5.2 Heating system Controllable with parameters
In this experiment we use a parameterized heating system to show that LODE-GP can learn physi- cally interpretable system parameters during GP hyperparameter training. The controllable heating system is depicted in Figure B It uses an input u(t) to control a heating element fi()) which ex- changes heat with an object fa(t). Two physical parameters a and determine how quickly the heat moves between the objects according to the ODEs:
(1) = a a (1i(() f2(()) + u(t)
(2(1) = -b. ((2(1) - (1)
We generate training data using a solution the differential equation (see Figure B with parameter values a = 3, b- 1. The LODE-GP can be constructed with parameters in the ODEs.
Learning the physical parameters and reconstructs the original values successfully with max- imal relative етrоr of less than 2.8% in ten training runs on data without noise. After adding white noise the data with standard deviation of 1% of the maximal signal, the parameters were success- fully reconstructed with maximal relative error of 5.3% in ten training runs. Again, the LODE-GP satisfies the original ODEs (with parameters on the noisy data, we refer the reader to Appendix B for visualization similar to Figure A This is bigger than in the previous example, as the LODE-GP only uses approximate parameter values for
and b 1) with an median error of 1e-2, trained
'The third data channel is particularly large due to its additional factor of = 9.81, which has a different scale by order of magnitude as the otherwise small data.

---

2
fi(t)
f2(1)
S2(1)
m(1)
f1(t)
0
-1
u(t)
0
2
6
8
10
12
14
Figure 5: (left) A visualization of the heating system with input u(() and states (1)))(((() (right) solution of the differential equations, which used to create the data.
fi(t) ....... f2(()
2z2 (t)
fs(1) (1)
0.5
22(
Ui
fi(1)
(2(1)
Sa(1)
8
10
Figure 6: (left) A sketch of the three tank system and (right) solution of the system.
and b, but still smaller than the еттоr of GP of 0.20. Our interpretation of these results is that the LODE-GP can learn the physically interpretable parameters.
5.3 Three tank system Non-controllable
We conclude the experiments with non-controllable fluid system where the water level in three tanks is changed by two pipes. The system non-controllable due to pipes' overlap over the center tank, whose changes directly affect the other tanks. This system requires multiple non-zero covari- ance functions in the latent GP to describe both the non-controllable subsystem and the two degrees of freedom, which corresponds the columns the SNFs matrix D via Lemma
1
0
07
D
1
0
0
0
0
We use solution the system of ODEs (see Figure to generate 25 datapoints, which we add white noise with standard deviation of 10% of the maximal signal.
The system ODEs again strictly satisfied by the LODE-GP with median error of 1e-5 in each ODE. We refer the reader to Appendix E for a visualization similar to Figure Hence, LODE-GP handle larger, non-controllable systems of ODEs, despite high noise.
Conclusion
In this paper we have introduced an algorithmical approach to automatically and symbolically create LODE-GP, class of covariance functions for GPs such that their realizations strictly follow a given system of linear homogeneous ODEs. We have proven that this approach is mathematically sound and verified its effectiveness in experiments. We have demonstrated that the learned posterior mean strictly satisfies the given system of ODEs up to numerical precision and discussed the consequences of this strict behavior e.g. in their extrapolation, leading to drastic behavior when both following the

---

ODEs and going back to the prior mean. Additionally, we automatically trained the physically inter- pretable parameters of the system and were able reconstruct their original values with low relative emor. The runtime of the LODE-GP is still asymptotically dominated by the On+)) Cholesky de- composition of the covariance matrix, despite bigger constants the ((n²) part of constructing the covariance matrix and a (usually very small) constant time precomputation of the covariance func- tion. Of course, approximate GPs are applicable to our covariance functions. Beyond applications where the dynamic behavior is known through corresponding ODEs, the combination with kernel search methods [HII] T9A 29, E. DJ has the potential to detect unknown dynamical behavior in data by systematic exploration and use these findings to generate new interpretable knowledge about the data through the corresponding system.
Acknowledgments and Disclosure of Funding
This research was supported by the research training group "Dataninja" (Trustworthy for Seam- less Problem Solving: Next Generation Intelligence Joins Robust Data Analysis) funded by the German federal state of North Rhine- Shine-Westphalia..
We want to thank the anonymous reviewers for their valuable feedback.
References
[1] William A Adkins and Mark G Davidson. Linear systems of differential equations. In Ordinary Differential Equations, pages 629-721. 2012.
[2] Mauricio Alvarez, David Luengo, and Neil D Lawrence. Latent force models. In Artificial Intelligence and Statistics, pages 9-16. PMLR, 2009.
[3] Mohamed Barakat. Purity Filtration and the Fine Structure of Autonomy. In Proceedings of the 19th International Symposium Mathematical Theory of Networks and Systems, 2010.
[4] Alain Berlinet and Christine Thomas-Agnan. Reproducing kernel Hilbert spaces in probability and statistics. Springer Science & Business Media, 2011.
[5] Fabian Berns and Christian Beecks. Automatic Gaussian process model retrieval for big data. In Proceedings of the 29th ACM International Conference on Information & Knowledge Man- agement, 2020.
[6] Fabian Berns, Kjeld Schmidt, Ingolf Bracht, and Christian Beecks. 3es algorithm for efficient Gaussian process model retrieval. 2020 25th International Conference on Pattern Recogni- tion (ICPR), 2021.
[7] Yves Bertin, Etienne Videcoq, Sophie Thieblin, and Daniel Petit. Thermal behavior of an electrical motor through a reduced model. IEEE Transactions energy conversion, 15(2): 129-134, 2000.
[8] LT Biegler, JJ Damiano, and GE Blau. Nonlinear parameter estimation: case study compari- son. AIChE Journal, 32(1):29-45, 1986.
[9] Ilias Bilionis, Nicholas Zabaras, Bledar Konomi, and Guang Lin. Multi-output separable Gaussian process: Towards an efficient, fully Bayesian paradigm for uncertainty quantification. Journal of Computational Physics, 241:212-239, 2013.
[10] Nathanael Bosch, Philipp Hennig, and Filip Tronarp. Calibrated adaptive probabilistic ODE solvers. In AISTATS, 2021.
[11] Roberto Calandra, Jan Peters, Carl Edward Rasmussen, and Mare Peter Deisenroth. Manifold Gaussian processes for regression. In International Joint Conference on Neural Nenworks (IJCNN) 2016.
[12] Ben Calderhead, Mark Girolami, and Neil Lawrence. Accelerating Bayesian inference over nonlinear differential equations with Gaussian processes. NeurIPS, 2009.
10

---

[13] Kian Chai, Christopher Williams, Stefan Klanke, and Vijayakumar Sethu. Multi-tast multi-task Gaussian process learning of robot inverse dynamics. NeurIPS, 2008
[14] Thomas Cluzeau, Victorita Dolean, Frédéric Nataf, and Alban Quadrat. Symbolic precondi- tioning techniques for linear systems of partial differential equations. 2011. URL https:
77hal. -es
fr/ha1-00644092
[15] Elizabeth J Cross and Timothy J Rogers. Physics-derived covariance functions for machine learning in structural dynamics. PrAC-PPpprsOnniie,, 54(7):168-773, 2021.
[16] Michiel JL. De Hoon, Seiya Imoto, Kazuo Kobayashi, Naotake Ogasawara, and Satoru Miyano. Inferring gene regulatory networks from time-ordered gene expression data of bacillus subtilis using differential equations. In Biocomputing 2003, pages 17-28.
[17] Claudia Drygala, Benjamin Winhart, Francesca Mare, and Hanno Gottschalk, Generative modeling of turbulence. Physics of Fluids, 34(3):035114, 2022.
[18] David Duvenaud. Automatic model construction with Gaussian processes. Phl) thesis, Univer- sity of Cambridge, 2014.
[19] David Duvenaud, James Lloyd, Roger Grosse, Joshua Tenenbaum, and Ghahramani Zoubin. Structure discovery in nonparametric regression through compositional kernel search. ICML, 2013.
[20] Duccio Fanelli and Francesco Piazza. Analysis and forecast of covid-19 spreading in china, italy and france. Chaos, Solitons & Fractals, 134:109761, 2020.
[21] Paterne Gahungu, Christopher W Lanyon, Mauricio A Alvarez, Engineer Bainomugisha, Michael Smith, and Richard D Wilkinson. Adjoint-aided inference of gaussian process driven differential equations. Advanves in Neural Information Processing Systems, 2022.
[22] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wil- son. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Advances in Neural Information Processing Systems, 2018.
[23] Alexandra Goeke, Christian Schilli, Sehastian Walcher, and Eva Zerz. Computing quasi-steady state reductions. Journal of Mathematical Chemistry, 50(6): 2012.
[24] Laurenz Göllmann. Lineare Algebra. Springer, 2017.
[25] ftalo Gomes Gonçalves, Felipe Guadagnin, Sissa Kumaira, and Saulo Lopes Da Silva. A machine learning model for structural trend fields. Computers & Geosciences, 149:104715, 2021.
[26] Thore Graepel. Solving noisy linear operator equations by Gaussian processes: Application to ordinary and partial differential equations. In ICML, 2003.
[27] Tanja Hernández Rodríguez, Anton Sekulic, Markus Lange-Hegermen and Björn Frahm. Designing robust biotechnological processes regarding variabilities using multi-objective opti- mization applied to biopharmaceutical seed train design. Processes, 10(5):883, 2022.
[28] Michael Hutchinson, Alexander Terenin, Viacheslav Borovitskiy, So Takao, Yee Teh, and Mare Deisenroth. Vector-valued Gaussian processes on Riemannian manifolds via gauge indepen- dent projected kernels. NeurlPS 2021, 34, 2021.
[29] Jan David Hüwel, Fabian Berns, and Christian Beecks. Automated kernel search for Gaussian processes data streams. In 2021 IEEE International Conference Big Data (Big Data), 2021.
[30] Carl Jidling, Niklas Wahlström, Adrian Wills, and Thomas B Schön. Linearly constrained Gaussian processes. In NeurIPS 2017, 2017.
[31] Hyunjik Kim and Yee Whye Teh. Scaling up the automatic statistician: Scalable structure discovery using Gaussian processes. In AISTATS, 2018

---

[32] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. ICLR, 2015.
[33] Edgar D Klenske, Melanie N Zeilinger, Bernhard Schölkopf, and Philipp Hennig. Gaussian process-based predictive control for periodic error correction. IEEE Transactions on Control Systems Technology, 24(1):1 10-121, 2015.
[34] Nicholas Krämer and Philipp Hennig. Linear-time probabilistic solution of boundary value problems. NeurlPS, 2021.
[35] Naresh Kumar and Seba Susan. COVID-19 pandemic prediction using time series forecasting models. In International Conference on Computing, Communication and Networking Tech- nologies (ICCCNT), 2020.
[36] Salah Labhalla, Henri Lombardi, and Roger Marlin. Algorithmes de calcul de la réduction de hermite d'une matrice à coefficients polynomiaux. Theoretical Computer Science, 161 69-92, 1996.
[37] Isaac E Lagaris, Aristidis C Likas, and Dimitris Papageorgiou. netroll methods for boundary value problems with irregular boundaries. IEEE Transactions on Neural Networks, 1165)104111049, 2000.
[38] Markus
Algorithmic linearly constrained Gaussian processes. In NeurlPS 2018, pages 2141-2152, 2018.
[39] Markus Lange-Hegeraann Linearly constrained Gaussian processes with boundary condi- tions. In AISTATS, 2021.
[40] Markus
and Daniel Robertz. On boundary conditions parametrized by an- alytic functions. In François Boulier, Matthew England, Timur M. Sadykov, and Evgenii V. Vorozhtsov, editors, Computer Algebra in Scientific Computing, pages 225-222222 245, Cham, 2022. Springer International Publishing.
[41] Maplesoft, a division of Waterloo Maple Inc.. Maple.
[42] Alonso Marco, Philipp Hennig, Jeannette Bohg, Stefan Schaal, and Sebastian Trimpe. Au- tomatic LQR tuning based on Gaussian process optimization: Early experimental results. In Second Machine Learning in Planning and Control of Robot Motion Workshop at International Conference on Intelligent Robots and Systems, 2015.
[43] MATLAB. The MathWorks Inc., Natick, Massachusetts.
[44] Morris Newman. The Smith normal form. Linear algebra and its applications, 254(1-3): 367-381, 1997.
[45] Ulrich Oberst. Multidimensional constant linear systems. Acta Applicandae Mathematica, 20 (1):1-175, 1990.
[46] Alban Quadrat. Grade Filtration of Linear Functional Systems. Acta Appl. Math., 127:27-86. 2013.
[47] Maziar Raissi, Paris Perdikaris, and George Em Karniadakis. Machine learning of linear differ- ential equations using Gaussian processes. Journal of Computational Physics, 348:683-693, 2017.
[48] Maziar Raissi, Paris Perdikaris, and George E. Karniadakis. Physics-informed neural networks: deep learning framework for solving forward and inverse problems involving nonlinear par- tial differential equations. Journal of Computational physics, 378:686-707, 2019.
[49] Sascha Ranftl. A connection between probability, physics and neural networks. In Physical Sciences Forum, volume 5, page 11. MDPI, 2022.
[50] Daniel Robertz. Formal computational methods for control theory. PhD thesis, RWTH Aachen University, 2006.
12

---

[51] Simo Särkkä. Linear operators and stochastic partial differential equations in Gaussian process regression. In International Conference on Artificial Neural Networks, pages 151-158, 2011.
[52] Jonathan Schmidt, Nicholas Krämer, and Philipp Hennig. A probabillistic state space model for joint inference from differential equations and data. NeurIPS 2021, 34, 2021.
[53] Michael Schober, David Duvenaud, and Philipp Hennig. Probabilistic ODE solvers with Runge-Kutte Kutta means. NeurIPS 2014, 2014.
[54] Michael Schober, Simo Särkkä, and Philipp Hennig. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing, 29(1):99-122, 2019.
[55] Henry John Stephen Smith. I. on systems of linear indeterminate equations and congruences. Proceedings of the Royal Society of London, (11):86-89, 1862.
[56] Edward Snelson, Zoubin Ghahramani, and Carl Rasmussen. Warped Gaussian processes. NeurlPS, 2003.
[57] Arno Solin, Manon Kok, Niklas Wahlström, Thomas B Schön, and Simo Särkkä. Modeling and interpolation of the ambient magnetic field by Gaussian processes. IEEE Transactions on robotics, 34(4):1112-113) 2018.
[58] Arne Storjohann and George Labahn. A fast Las Vegas algorithm for computing the Smith normal form of polynomial matrix. Linear Algebra and its Applications, 253((-3):155-173 1997.
[59] PARI/GP version 21.1 2. The PARI Group, Univ. Bordeaux, 2019. available from http:
/ waari.math.u-aa www.m.u-bordeaux..e....
[60] The Sage Developers. SageMath, the Sage Mathematics Software System (Version 9.4), 2021. available from https://wwis wttp:///wwwwwgwmtthhhacaaama www.math.eeaa.ee...eemmm........
org.
[61] Silja Thewes, Markus Lange-Hegermann, Christoph Reuber, and Ralf Beck. Advanced Gaus- sian process modeling techniques. Design of Experiments (DoE) in Engine Development. Ex- pert Verlag, 2015.
[62] Filip Tronarp. Simo Särkkä, and Philipp Hennig. Bayesian ode solvers: The maximum a posteriori estimate. Statistics and Computing, 31(3);1-18, 2021.
[63] Selvakumar Ulaganathan, Ivo Couckuyt, Tom Dhaene, Joris Degroote, and Eric Laermans. Performance study of gradient-enhanced Kriging. Engineering with computers, 32(1):15-34, 2016.
[64] Ph van Milligen, V Tribaldos, and JA Jiménez. Neural network differential equation and plasma equilibrium solver. Physical review letters, 75(20):3594, 1995.
[65] Gilles Villard. Computation of the Smith normal form of polynomial matrices. In International symposium on Symbolic and algebraic computation, pages 209-217, 1993.
[66] Gilles Villard. Fast parallel computation of the Smith normal form of polynomial matrices. In International symposium on Symbolic and algebraic computation, pages 312-317, 1994.
[67] Gilles Villard. Generalized subresultants for computing the Smith normal form of polynomial matrices. Journal of Symbolic Computation, 20(3):269-286, 1995.
[68] Gilles Villard. Fast parallel algorithms for matrix reduction to normal forms. Applicable Algebra in Engineering, Communication and Computing, 866)5511-337, 1997.
[69] Niklas Wahlström, Manon Kok, Thomas B Schön, and Fredrik Gustafsson. Modeling magnetic fields using Gaussian processes. In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 3522-3526, 2013.
[70] Peter Wangersky. Lotka-volterra population models. Annual Review of Ecology and System- atics. 91):119-2188 1978.
13

---

[71] Jon Wilkening and Jia Yu. A local construction of the Smith normal form of a matrix polyno- mial. Journal of Symbolic Computation, 46(1):1-22, 2011.
[72] Christopher Williams and Carl Edward Rasmussen. Gaussian processes for machine learn- ing. MIT press Cambridge, MA, 2006.
[73] Eva Zerz. Topics in Multidimensional Linear Systems Theory, volume 256 of Lecture Notes in Control and Information Sciences. London, 2000.
14


Reviews:
[{"review_id": "ZZ9mIr2gAKU", "reviewer": "Reviewer_s3hy", "summary": "The paper presents a way to algorithmically construct Gaussian Process\ncovariance functions in such a way, that realizations strictly follow a given\nlinear ordinary differential equation. The proposed construction via smith\nnormal form algorithms overcomes existing limitations on controllability of the\nsystem, which exists in other approaches. The resulting LODE-GPs are compared to\nclassic GPs in three experiments. \n", "questions": "My most relevant question / suggestion is about the relation to parameter inference problems, stated above under \"weaknesses\".\n\nSome questions / comments on probabilistic numerics (PN):\n- l. 23 writes that probabilistic ODE solvers allow to \"learn\" non-linear ODEs,\n  but it is not quite clear to me in which sense this holds (i.e. inferring the\n  solution or inferring parameters). To my understanding, those methods \"solve\"\n  or \"estimate the solution of\" ODE initial value problems (IVPs)\n- Solving IVPs with non-linear ODEs _requires_ approximations, since they can, in\n  general, not be solved in closed form. It might thus be misleading to state\n  approximations as limitations of the solver, rather than a necessity of the\n  general problem. In contrast, IVPs with linear ODEs can be solved in closed\n  form, and therefore PN does not apply since no numerical error needs to be\n  estimated.\n- l. 70: What is meant with \"complex calculations\"? It could be helpful to\n  clarify.\n- Page 6, footnote: What is meant with \"free functions in their solution set\"?\n  Depending on the meant object, [47] might possibly be in disagreement since it\n  describes inference of time-varying system parameters (= unknown functions).\n\nI want to highlight that these comments are only meant as helpful suggestions and not as necessary changes.\nThe resolution about the relation to parameter inference, commented on above, as\nwell as possibly the extension of the baseline for the experimental evaluation\nare of higher relevance.\n\nMinor comments:\n- l. 18: \"such _as_ in\"\n- l. 88 contains a lengthscale l, but it is not mentioned.\n- Around Lemma 1 it would be helpful to point the reader to the proof in the\n  appendix.\n", "limitations": "The authors adequately discussed the limitations of the proposed method.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\nThe presented method appears to be novel and I believe that it is a significant\ncontribution to research on Gaussian processes and generally machine learning in\nthe context of dynamical systems. The technical statements in the paper appear\nsound. Related work is adequately cited (some comments below). Additionally, the\npaper is very well-written and the method is presented in a clear way; the\nbipendulum example used throughout the paper is very helpful, and the\nexperimental setup is discussed in much detail in the appendix.\n\n# Weaknesses\nThe main criticism I have for the paper is about the comparison only to a\nclassic Gaussian process in the experiments. While I believe that this _might_ be\nappropriate, it would be helpful to comment more extensively to why this is the\nonly reasonable comparison, and to define the problem setting more formally /\nprecisely.\n\nFor example, the experiment settings in section 5 can also be interpreted as\nparameter / initial value inference problems, which could be approached with\ne.g. least-squares regression with classic ODE solvers (e.g. [Biegler]) or\nprobabilistic numerics [Tronarp], or (probabilistic) gradient matching (e.g.\n[Wenk]). Note that the former two approaches would further come with the\nadvantage of having linear O(n) scaling (due to not modeling covariances or to\nusing Gauss--Markov priors, respectively). Clarfying the problem setting further\ncould be helpful for readers to understand why these do not apply; otherwise\nsome of these approaches might be relevant baselines.\n\n[Biegler] Nonlinear parameter estimation: A case study comparison, Biegler et al, 1986  \n[Wenk] ODIN: ODE-informed regression for parameter and state inference in time-continuous dynamical systems, Wenk et al, 2020  \n[Tronarp] Fenrir: Physics-Enhanced Regression for Initial Value Problems, Tronarp et al, 2022\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "CMwNP3HrfM2", "reviewer": "Reviewer_JWZD", "summary": "The authors introduce a new method which tackles the important and much discussed problem of learning the solutions to systems of ordinary differential equations from noisy data. The authors construct a Gaussian process prior for the solutions to the equation, which has been a common approach in previous work in the area. In particular the authors attempt to solve the problem in the case of linear differential operators for systems which are underspecified and uncontrollable. The proposed method relies on the symbolic calculation of the Smith normal form of the differential operator representing the system, which the authors use to construct a covariance over both the system outputs and inputs, ensuring that samples obey the differential equation in question. The proposed method is validated on 3 synthetically constructed physical systems, providing convincing evidence that the method works.\n", "questions": "Below are my questions/suggestions for the authors, which are mostly a summary of what is discussed in the weaknesses section. I've numbered the points separately to allow the authors to respond more easily. \n\n1. Please can the authors comment on the significance of the ability to model noncontrollable systems, and provide some real world motivation? In doing so could the authors expand the discussion of how this work represents a significant contribution over [34]?\n2. How does the proposed method compare to existing models of this type from problems that are fully specified?  \n3. Can the authors comment on the relationship between this model and LFMs? If they are indeed equivalent, the authors should discuss this in the text, and add a comparison for system 3 if possible?\n4. I suggest that the authors add standard errors/deviations to table 2, and indicate if the standard errors are narrow enough to indicate one model being significantly better than the other. \n   \nIf the above points are addressed I would be more than happy to raise my score. \n\n\nBelow are some very minor typos/formatting/wording issues:\n\n* Line 23 \"allow to\"\n* Line 193 not sure what \"All experiments ran numerically stable and robust using Adam [28].\" means, please rephrase.\n* Line 233 \"The LODE-GP can be constructed using parameters.\" this is not clear, please rephrase.\n* Many of the figures don't have titles/axis labels, so are difficult to read, although mostly the relevant information is in the caption/legend. For figure 2 particularly it is not at all clear which plot is which output, presumably it is f1, f2, u from left to right, but please clarify via titles.", "limitations": "I do not believe there are any ethical concerns with this work or problems with potential negative societal impact. \n\nAside from what I have already mentioned, I believe the authors reasonably fairly represent the limitations of their work, although I do think the systems shown in the evaluation are somewhat similar, and it would be good to get some idea of how these models would perform on real world data, not just synthetic systems. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths \n---------\n\nOverall I enjoyed reading the paper, and think it has a number of strengths. The paper for the most part is very clearly written, and introduces the problem, and the pre existing solutions, well. I particularly like the use of explicit examples at key points in the exposition which help the reader to understand how the proposed method is applied in practice. The figures are mostly clear and useful. \n\nThe idea presented in the paper is quite an elegant one in my opinion, and seems to be a sensible way to tackle the problem. The fact that the method can be used to construct priors which consist of exact solutions to the given system is compelling. \n\nThe method seems to work on a number of different problems, and in some cases shows good performance relative to standard GPs, although I do have some concerns about the evaluation, as listed below. \n\nI am not expert in the theoretical tools such as the Smith normal form used in the paper, and so I can not vouch for the correctness of those parts with 100% certainty. That being said, the results given in the paper are correct to the best of my knowledge.\n\nWeaknesses\n----------\n\nThere are a few weaknesses with the paper, mostly related to the evaluation of, and motivation for, the proposed methods. \n\n**Motivation for uncontrollability**\nThe authors say that their methods is an improvement on [34] because it can deal with uncontrollable systems.  As the authors note, in 2 of the experiments which are controllable, the proposed method generates exactly the same results as the previous work of [34]. To me, there is not sufficient discussion of the reasons that the uncontrollable case is an important problem, and thus why the proposed method represents an important contribution. I am not an expert in control theory, so it may be that I am unaware, but I think it would be good to motivate this problem more clearly in the text, and explain why it is important to be able to model such systems. Additionally for the experiment where the proposed model is novel (three tanks), the performance seems to be worse than a vanilla GP in terms of it's representation of the data. \n\n**Comparison with related work**\nPresumably the proposed method can work for problems that are well specified, as well as those discussed in the text that are underspecified. Currently there is no indication of how this method performs in that case, and how well it performs relative to other prior work. I think it would improve the paper considerably to include some comparisons. \n\nThe authors claim that \"No other state of the art in probabilistic ML exists for our class of differential equations\" I am not entirely sure that this is correct. I believe that is is possible to solve some of the problems using the framework of latent force models (LFMs) [Alvarez et al. 2009] which represent the input functions as GPs, and use the Green function of the differential operator to map to exact outputs. Most work on LFMs deals ODE systems which are not coupled, i.e. the part of $A$ relevant to the outputs is diagonal, so it may be additional effort to get systems 1 & 2 running in the LFM framework, but I believe system 3 can easily be recast in terms of LFMs. System 3 is equivalent to the model in equation 3 of [Alvarez et al. 2009], when $D_q = 0$, $B_q=0$ and $S_{rq}=\\\\begin{pmatrix}\n    -1 & 0\\\\\\\\\n    -1 & -1 \\\\\\\\\n    0 & -1\n\\end{pmatrix}$, the expressions given in later equations can then be used to compute the covariance of the outputs. It would be great to have a comparison between the proposed model and LFMs for system 3, at least i think a comment on this equivalence is necessary. Again I may be mistaken in this point, but I would appreciate a comment from the authors if so. \n\n**Evaluation**\nA minor point, but I'm not sure \"ODE error\" shown table is interesting to compare with the GP. The metric is useful as a check that the model is working, but it doesn't seem fair to compare with the GP since the information about the differential equation is not available to the GP, and the error should be small for the proposed model almost by definition. \n\nTable 2 shows results for 10 repeats, with bold text indicating best performance, but does not include the standard deviation or error on these results. This means it is not possible for the reader to assess the significance of the results. Please add sandard errors or deviations to this table\n\n##### References\n- [Alvarez et al. 2009] Alvarez, Mauricio, David Luengo, and Neil D. Lawrence. \"Latent force models.\" Artificial Intelligence and Statistics. PMLR, 2009. ", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "AqtyiIZeZbD", "reviewer": "Reviewer_6STF", "summary": "This paper proposes a method, based on Smith normal form, to construct Gaussian processes indexed by R such that they satisfy a specified differential equation. \n\n", "questions": "\n\nQ1 - I guess \\ell_1 should really be \\ell_2 in entry (3,2) of A in eq. (1) ? ", "limitations": "Its fine. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\nS1 - The paper is clearly written. \n\nW1 - While the proof of Lemma 1 is indeed trivial enough to defer to appendix, it may be useful to point out in the main text that (a) the proof is indeed given in appendix, and (b) the covariances are obtained by expanding f in a basis of the null-space of d with i.i.d standard normal coefficients. \n\nW2 - The entire noise free data is not shown in fig 2. It woudl certainly be interesting to see how well the method extrapolates. \n\n\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "14xNJrzYk-", "reviewer": "Reviewer_Wp1C", "summary": "The paper describes how to constrain realisations of a Gaussian process $f \\sim GP(m, k)$ to satisfy $A f = 0$, where $A$ is a linear, matrix-valued differential operator with constant coefficients.\nThis is done by computing the Smith normal form of $A$, which is a (symbolic) decomposition of $A$ into\n$U D V = A$, or equivalently $D = U^{-1} A V^{-1}$, and identifying the null space of $D$.\n\nThis is most closely related to the work by Hegermann (2018) and Hegermann (2021), who tackle a similar problem using a Gröbner basis.\nThe difference between the submitted work and the aforementioned papers is the use of a Smith normal form, which allows for non-controllable systems, and can be implemented with standard software. The Smith normal form also restricts the class of problems that can be tackled to linear, homogenous, ordinary differential equations with constant coefficients.", "questions": "\n* Abstract: what does \"novel algorithmic and symbolic construction\" mean? It is written frequently, but I don't understand \"algorithmic and symbolic construction\".\n* Line 87: Please explain the following notation in the paper $\\mathbb{R}_\\preceq$. I suppose it refers to some sort of positive-definiteness?\n* Line 223f: \"We additionally analyse the training runtime and see that...\": are those results in some figure/table? I can't seem to find the results\n", "limitations": "It would be useful if the limitations of the method would be discussed more centrally. What do I mean by this? Line 68 states \"but are restricted to ODEs\". Lines 72f. describe the limitations as \"restricted to linear ODEs\". Line 140 states: \"Consider a system of linear homogenous ODEs with constant coefficients\". Line 172 says: \"This is generally impossible for PDEs\". Line 255 states: \"follow a given system of linear homogeneous ODEs\".\nIt feels like a reader has to search the entire document to learn what the limitations of the method are. \n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\nOverall, the methodology seems to be sound, the related work seems to be discussed adequately, and the experiments support the claims.\nI would expect that the method described in the paper will be useful to some people. The paper is at the interface of differential equation models and Gaussian processes, which seems to be a relevant subject for Neurips.\n\n\n## Weaknesses\n* One criticism of the method is that the algorithm is restricted to linear, homogenous ODEs with constant coefficients and that there may not be much hope to generalise to nonlinear/nonhomogenous problems or even partial differential equations.\n* From a didactical viewpoint, I would appreciate it if the paper could provide more intuition for Lemma 1. Maybe a few equations would help.\n* This is minor, but I would like to see samples from the constrained Gaussian process, perhaps in one of Figures 2, 4, 5, or 6?\n* Why does only Figure 2 show the credibility intervals? The remaining figures only discuss the mean. Is there a specific reason?", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept