Physics-Informed Variational State-Space Gaussian Processes
Authors: Oliver Hamelijnck
Arno Solin
Theodoros Damoulas
Venue: NeurIPS 2024 poster (2024)
Keywords: gaussian processes
variational approximations
state space gaussian processes
physics informed gaussian processes
Primary area: probabilistic_methods

Abstract: Differential equations are important mechanistic models that are integral to many scientific and engineering applications. With the abundance of available data there has been a growing interest in data-driven physics-informed models. Gaussian processes (GPs) are particularly suited to this task as they can model complex, non-linear phenomena whilst incorporating prior knowledge and quantifying uncertainty. Current approaches have found some success but are limited as they either achieve poor computational scalings or focus only on the temporal setting. This work addresses these issues by introducing a variational spatio-temporal state-space GP that handles linear and non-linear physical constraints while achieving efficient linear-in-time computation costs. We demonstrate our methods in a range of synthetic and real-world settings and outperform the current state-of-the-art in both predictive and computational performance.

Physics-Informed Variational State-Space Gaussian Processes
Oliver Hamelijnck
University of Warwick
oliver hhw.hhhhewwwwkckkwww.. ac. uk
Arno Solin
Theodoros Damoulas University of Warwick
Aalto University arno wono.solin@aalt...
t. ww..oulas@warriickk.
Abstract
Differential equations are important mechanistic models that are integral to many scientific and engineering applications. With the abundance of available data there has been a growing interest in data-driven physics-informed models. Gaussian processes (GPs) are particularly suited this task as they model complex, non- linear phenomena whilst incorporating prior knowledge and quantifying uncertainty. Current approaches have found some success but are limited as they either achieve poor computational scalings or focus only on the temporal setting. This work addresses these issues by introducing a variational spatio-temporal state-space GP that handles linear and non-linear physical constraints while achieving efficient linear-in-time computation costs. We demonstrate our methods in range of synthetic and real-world settings and outperform the current state-of-thed of-the-art in both predictive and computational performance.
Introduction
Naive vs. Plysiic-Ifformed State-Space GP (PHYSS-GP)
Physical modelling is integral modern science and engineering with applications from climate modelling [62] to options pricing [6]. Here, the key formalism to inject mechanistic physical knowl- edge differential equations (DEs), which given initial and/or boundary values, are typically solved numerically [8]. In contrast machine learning is data-driven, and aims to learn latent functions from observations. However the increasing availability of data has spurred interest combining these traditional mechanistic mod- els with data-driven methods through physics-informe machine learning. These hybrids approaches aim improve predictive accuracy, computational efficiency by leveraging both physical inductive biases with observations [30, 44].
Br
No.
=
Ingut,
Time,
Figure 1: The state-space formal- ism allows for linear-time infer- through Gaussian processes (GPS). GPS stochastic processes ence in the temporal dimension. and are data-centric approach that facilitates the quantification of uncertainty. Recently AUTOIP was proposed in order to integrate non-linear physies into [40], where solutions to ordinary and partial differential equations (ODEs, PDEs) are observed at finite set of collocation points. This is an extension of the probabilistic meshless method (PMM, [12]) to the variational setting such that non linear equations can be incorporated. Similarly, [4] introduced HELMHOLTZ-GP. that constructs GP priors that adhere curl and divergence-fre constraints. Such properties are required for the
principled way to incorporate prior physical knowledge
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

successful modelling of electromagnetic fields [61] and ocean currents through the Helmholtz decom- position [4]. These approaches enable the incorporation of physics but incur cubic computational complexity from needlessly computing full covariance matrices, as illustrated in Fig. 1. For ODEs (time-series setting), extended Kalman smoothers incorporate non-linear physics (FKS) [65, 34] and recover popular ODE solvers whilst achieving linear-in-time complexity through state-space GPS [58, 25].
In this work we propose unified physics informed state-space GP (PIYSSSGP)) that is a probabilistic models where mechanistic/physies knowledge is incorporated an inductive bias. We can handle both linear and non-linear PDEs and ODEs whilst maintaining linear-in meear-in-time in-time computational efficiency. We additionally derive state-space variational inference algorithm that further reduces the computational cost in the spatial dimension. We recover EKS, PMM. and HELMHOLTZ-GP special cases, and outperform AUTOIP in terms of computational efficiency and predictive performance. summary:
1. We derive state-space GP that handle spatio-temporal derivatives with computational complexity that is linear in the temporal dimension.
2. With this we derive a unifying state-space variational inference framework that allows the incorporation of both linear and non-linear PDEs whilst achieving linear-in-tine complexity and recovering state-of-the-at methods such as EKS, PMM and HELMHOLTZ-GH
3. We further explore three approximations, namely a structured variational posterior, spatial sparsity, and spatial minibatching, that reduce the cubic spatial computational costs to linear.
4. We showcase our methods on a variety of synthetic and real-world experiments and out- perform the current state-of-the-art methods AUTOIP and HELMHOLTZ-GP both in terms computational and predictive performance.
Code to reproduce experiments available at https: /ttt::://githbb... wtattthuu.co//haaaaaaeeee com/ http///hhmmlijjc/ppppaaaaae
Background on Gaussian Processes
Gaussian processes A GP is distribution infinite collection of random variables such that any finite subset is jointly Gaussian [50]. Given observations X E Nx and y E IRN then
ply, f10) II, p(un O( p(f(()
(I)
joint model where p is zero mean GP prior with kernel K(:-), f(X) p([() 0, K(X,X) and are (hyper) parameters. We are primarily concerned with the spatio- temporal setting where we observe N, temporal and N₂ spatial observations Tz,s E IR Vt.s E I on spatio-temporal grid. Under Gaussian likelihood, all quantities for inference and training are available analytically and, naïvely, carry a dominant computational cost of O(N, N₂ ))) For time series data, an efficient way construct GP over f (and its time derivatives) is through the state-space representation of GPs. Given a Markov kernel, the temporal GP prior can be written as the solution of discretised linear innee-ivvariatt stochastic differential equation (LTI-SDE), which at time l is
f A Qk and yk fr P(U Hf)))
where A is transition matrix, 4s is Gaussian noise, H is an observation matrix, and is a d- afc-) z²5() . With appropriately de- dimensional vector of temporal derivatives f [f(-), 8z 8
IT. signed states, matrices and densities, SDEs of this form represent large class of GP models, and Kalman smoothing enables inference in ON, d3), see [56]. In the spatio-temporal setting, when the kernel matrix decomposes as Kronecker product K K, @ K,, then with Markov time kernel, a state space form is admitted. This takes a particularly convenient form where the state f. ((((()))
f((X.) No. t)] and inference requires (N, d))), see [60].
Derivative Gaussian processes One main appeal of GPS is that they are closed under linear [1] D.D. be linear functional that computes D = d, d, space-time derivatives Dt 3 ... and D, [-]
operators. Let D
with D, [-]
8- 82
0 d..... then at finite set of index points, the
joint prior between and its time and spatial derivatives is
p(f(X)) N (ff O, DK((XX)D'')
(3)

---

where f(X) D f() and D is the adjoint of D, meaning it operates on the second argument of the kernel [54]. When joindy modelling a single time and space derivative (d₂ dg 1) the latent functions are f [f, 8f of 8²r 1 and the kernel
05
K
T
K
K
= D=DK(XXX)D*=
K
0s
K
Us
a
8 K
T
DE
Us
d: 8²
K
8
K
K
a
St
318
This is a multi-output prior whose samples are paths of f with its corresponding derivatives. This prior is commonly known as derivative GP and has found applications in monotonic GPs [51], input-dependent noise [41, 67] and explicitly modelling derivatives [59, 17, 43]. State-space GPs can be employed in the temporal setting since the underlying state computes (x) with its corresponding time derivatives. In Sec. 3.1, we extend this to the spatio-temporal setting.
Physics-Informee State-Space Gaussian Processes (PHYSS-GP)
We now propose flexible generative model for incorporating information from both data observations and (non-linear) physical mechanics. We consider general non-linear evolution equations of the form
0f
g(No f)
No f =0
Ot
with appropriate boundary conditions, where f RF - IR is the latent quantity of interest and No is non-linear differential operator [49]. We assume that g R.P.D + is measurable, and is well-defined such that there are sensible solutions to the differential equation [25]. We wish to place a GP prior over f and update our beliefs after 'observing' that should follow the solution of the differential equation. In general this is intractable and can only be handled approximately. By viewing Eqn. (4) as a loss function that measures the residual between af and the operator Ni then the right hand side (0) are virtual observations. The PDE can now be observed at a finite set of locations known collocation points. This is soft constraint (i.e. f is not guaranteed to follow the differential equation), but can handle non-linear and linear mechanisms. However, there are special cases, namely curl and divergence-free constraints, that can be solved exactly. This follows from properties of vectors fields, where f defines a potential function where linear combinations of its partial derivatives define vector fields that enforce these properties. To handle both of these situations we propose the following generative model
F..
W. [X,,
J ~ CP Ky),
(5)
Linear Mixing
Independent GP Priors
Ho F, + Eo.
(C)
FF)))
H8 Fr, + EB.
(6)
Data
Collocation Points
Boundary Values
where fg are derivative GPS (see Eqn. (3)) that are linearly mixed by W ((PD)X(QQ)) and Y()), ,00) RNP are observations and collocation points over the P outputs and ( (RNx D) are boundary values over the derivatives of each output. The observation matrices Ho. H8 simply select the relevant parts of F.. For further details on notation see App. A. In many case we want to observe the solution of the differential equation exactly, however in some cases it may be required to add observation noise Eç to the collocation points, whether for numerical reasons or to model inexact mechanics. This is flexible generative model where different assumptions and approximations will lead to various physics informed methods such as AUTOIP, EKS, PMM, and HELMHOLTZ-CE GP that will develop state space algorithms for. Additionally it is possible to learn missing physics by parameterising unknown terms in Eqn. (4) through the GP priors in Eqn. (6) (see App. B.2).
Example 3.1 (EKS Prior and PMM). We recover EKS style generative models (see Hennig et al. (25)) when the mixing weight is identity W I, and 6C, 6g _ 0, and the non-linear trans- form g is linearised. Let the prior be Markov p(f) I p(f (1) with marginals p(f.)
N ffx m P ). By taking a first-order Taylor linearisation g(f.) ~ g(m ) + dg(m; 8fg with dm,
3

---

ōff. ~ N (0, P
the joint is
p(
f gR
=
fg gi
m
0g(m,
P
0g 0m,
99k
(m
0m
This is now form that can directly be implemented into extended Kalman smoothing algorithm [63]. When Q the state f constructed by stacking the individual states of each latent [56]. With linear ODEs EKS coincides with PMM.
Vector Fields in 2D). Let v = [U2, Vas Vey! denote a 3D-vector field, then curl indicates the tendency of vector field rotate and divergence at specific point indicates the tendency of the field to spread out. Curl and divergence-free fields follow
Example 3.2 (HELMHOLTZ-GP and Curl and Divergence-
XXV=00 0 (curl free), V =0 (div. free)
(8)
where V
De, Dag Two basic properties of vector fields state that the divergence of a curl field and the curl of a derivative field are zero [3]. Let [fi, f2] be scalar potential functions then
Dt
Vurl V11 (curl free), Vdiv V V f2 (div. free)
define curl and divergence- free fields. 2D this simplifies using the grad and rot operators over [Vey Us2] (see [4]). Placing GP priors over Sq we incorporate this into Eqn. (6) by defining
0
0 1
05 0s Os,
W grad
H.
W.. =
H where H selects
(10)
0 1
-1 0
HELMHOLTZ-GP is defined as the sum of GP priors over 2D curl and divergence-free fields [4].
3.1 A Pratio-Temporal State-Space Prior
The generative model in Eqn. (6) contains complications: i) it includes potential non-lineararities, and ii) the independent priors are defined over latent functions with their partial derivatives which substantially increases the computational complexity. We wish to tackle both issues through state- space algorithms that are linear-in-time. We begin by deriving a state-space model that observes derivatives across space and time (see App. A.3 for the simpler time-series setting). In Sec. 3.2 we further derive state-space variational lower bound that will enable computational speeds up in the spatial dimension.
First, we show how Kronecker structure in the kernel allows us to rewrite the model as the solution to an LTI-SDE. From the definition of D. the separable covariance matrix has repetitive structure that can be represented through a Kronecker product. The gram matrix is
K((,, D' =K. (Xt, xx)
X K
(Xa, X²)
(11)
K.
K. D. D. and D. (D. excludes the underlying latent function.
where K.P
D.K.
To find a Kronecker form of the gram matrix over X, we will exploit the fact that X is on spatio- temporal grid and that the kernel is separable. Due the separable structure derivative over either the spatio (or temporal) dimension only affects the corresponding kernel, and so when considering X, the gram matrix is still Kronecker structured:
a Os
K K(x,x) x) (x(,x) x)
(x,
K(X, KX X,)
K, (X,, X,). (12)
Os
as
Os
The full prior over (a permuted) X is now given as
p(f()))= (o. (P K. (X,,
This is the form of a spatio-tempone Gaussian process with derivative kernels that can be immediately cast into a spaee form as Eqn. (2) where H I, as we want to observe the whole state, not just f. The marginal likelihood and the GP posterior can now be computed using standard Kalman filtering and smoothing algorithms with computational time of O(N, (N₂ d, d)}s). Inference in PHYSS-GP now follows Ex. 3.1 by recognising that the filtering state consists of the spatial points with there spatio-temporal derivatives. The EKS prior Ex. can now be simply extended the PDE setting by placing colocation points on spatio-temporal grid [35].

---

3.2 A State-Space Variational Lower Bound (PHYSS-VGP and PHYSS-EKS)
We now derive a variational lower bound for PHYSS-GP that maintains the computational benefits of state-space GPS. This acts as an alternative way of handling the non-linearity of in Eqn. (6), and will also enable the reduction of the cubic spatial computation complexity Sec. 4. We start by focusing on the single latent function setting (Q= = 1) and collect all terms that relate to observations in Eqn. (6) with Y f]- II, p(y,,) Ho F..) p(O (Fx)) )) p(y! [Hg F,). VI frames inference as the minimisation of the Kullback- Leibler divergence between the true posterior and an approximate posterior, which leads the optimisation of the ELBO [28]:
p( ( ())
arg max L C= E ((()) log a( <)
(13)
q(f)
where we define the approximate posterior gf {) N( m, S) free-form Gaussian with & = (m, S) and m E RDDNx11 S E DDNND, The aim is to represent the approximate posterior as state-space GP posterior, which will enable efficient computation of the whole evidence lower bound (ELBO). We will achieve this through the use of natural gradients. The natural gradient preconditions the standard gradient with the inverse Fisher matrix, meaning the information geometry of the parameter space is taken into account, leading to faster convergence and superior performance [2, 31, 27J. For Gaussian approximate posteriors the natural gradient has simple form [26]
DC 0pks
ÖELL 0pp
A AR- + 3
(1 3) A +
+ n - -
(14)
where À -S- m, 1/28-1) and μ = (m, m m T + S) are the natural and expectation parameteri- sations. This known as conjugate variational inference (CVT) as X represent the natural parameters for the conjugate prior 17 [31, 10, 20, 72]. For now, we will assume that the likelihood is conjugate to ensure that [Au]e is p.s.d, this will be relaxed in Sec. 5. The derivative of the ELL is
ÖELL NNN.N. N. a Ju
aluj E, [log (o ((.) ].
(15)
where the expectation is under (((.s))) a D dimensional Gaussian over the spatio-tempore deriva- tives at location Xt.s- Within the sum, the only elements of [ul] whose gradient will propagate through the expectation are the D x D elements corresponding to these locations. These points are JELL has some (permutated) plock-diagonal structure, hence Eqn. (14) can be written unique and so pp as
(f)
x II NO [N(Ý, f. V.)] p(F)
(16)
where Y, is D-dimensional. The natural gradient update, i.e. (fc)) in moment parameterisation, can now be computed using Kalman smoothing in ON,. (N-dd.-d3)) Collecting Y = vec( (Y.)). V blkdiag (Va)). then the ELBO can also be computed efficiently by substituting this form of ((f.) in
No.Vo
No Σ
C Σ E ((R [logr( Y (c.s) [(.ss)]] 6,8
E alf))
log N(Y, F.. (.F..VV) V.) + logp(Y V)
(17)
where the first two terms only depend on (DD..) and the final term is simply a by-product of running the Kalman filter, leading to dominant computational complexity of O(N (N₂ - de d)3).
This cost is linear in the datapoints (N) because the expected log likelihood above decomposes across all opatio-temporal locations. In summary have shown that natural gradient is equivalent updating a block -iiggoaa likelihood that decomposes across time; hence the approximate posterior is computable via Kalman smoothing algorithms. Extending multiple latent functions (Q > 1) we define a full Gaussian approximate posterior that captures all correlations between the latent functions (f...... .fQ) & N (ft..... .. m, S) where m E (Nx))x I ,S E R ([Nxx))(NNX)))
All the observation models in Eqn. (6) decompose across data points, hence Eqn. (16) is still block-diagonal and decomposes across time, except now each component is of dimension Q x N, as it encodes the correlations of spatial points and their spatio-temporal derivatives across the latent functions. We denote this model PHYSS-VGP and PHYSS-EKS when using EKS prior (see Ex. 3.1.
5

---

Theorem 3.1. Let the approximate posterior be (full) Gaussian "ff,,... .t) N (f...... .f. m, 8) where R(NXQ)X1, E [(NxQ)((XX)) When g is linear single natural gradient step with 3 recovers the optimal solution p(f,,..., fo Y).
We prove this in App. A.5.4. This result not only demonstrates the optimality of our proposed inference scheme in the linear Gaussian setting, but confirms that we recover batch models like PMM and HELMHOLTZ-GP, as well as EKS (see Ex. 3.
4 Reducing the Spatial Computational Complexity
We now propose three approaches that reduce the cubic computational complexity in the number of spatial derivatives and locations. The first augments the process with inducing points that alleviate cubic costs associated with Ng. The second is a structured variational approximation that defines the approximate posterior only over the temporal prior and alleviates cubic costs associated with d,. Finally, we introduce spatial mini-batching that alleviates linear N, costs. When used in conjunction, the dominant computation cost is O (N, d2- (M, d These approximations are not only useful for the state-space setting and can readily be applied to reduce the computational complexity for batch variational models (such as AUTOIP). See App. B.1 for more details.
this first approximation, denoted by PHYSS- SVGP, augment the full prior p(f) with inducing points. By defining these inducing points spatio-temporal grid, we will show that we can still exploit Markov conjugate operations through natural gradients. Let a Du E R RMX be inducing points at locations Z E RMXF, From the standard SVGP formulation [27], the ELBO is
Spatio-Temporal Inducing Points (PHYSS-SVGP)
P(Y f (()) q(u)
L = E q(T.a)
log
(18)
where q(f, ú) A p(f ū) q(ū) By defining the inducing points on spatio-temporal grid at temporal locations X, E TEE and spatial Zg (M then the marginal p(f ū) is Gaussian with mean
PF P I (X,,,,,,2 Z.) (K? ( (Zg, Z,))
(19)
and variance given in Eqn. (41). This Kronecker structure allows us to again 'decouple" space and time, leading to natural gradient updates with block size M, D, reducing the computational complexity to O(N (M, d. ))²). For full details, see App. A.5.I.
Structured Variational Inference (PHYSS-SVGP.) This second approximation, denoted as PHYSSS- SVGP,, defines the inducing points only over the temporal derivatives. This is a useful approximation as it can drastically reduce the size of the filter state, making it more computationally and memory efficient. We begin by defining the joint prior as
(FF D, f) = p(F D. D.f)) P(D,
where p(F
|D f) is Gaussian conditional with mean
K. KX Z,) De (X.,,
E F D.f
I ② K. Z,, Z.) Z.)-11 D. with k,,X,,
We then define a structured variational posterior
(f, D. A ((F)D,))
Substituting this into the ELBO we see that all the terms with the prior spatial derivatives cancel
P P(Y
E q
log
Again, the marginal g(D, f) maintains Kronecker structure, enabling Markov conjugate operations, leading to a computational cost of O(N d, (N₂ d)3), see App. A.5.2. These variational approxi- mations can simply be applied to non-state-space variational approximation, App. B.

---

Spatial Mini-Batchiin A standard approach for handling big data is through mini- mini-batching where the ELL is approximated using only data subsample [27]. Directly appling mini-batching would be of little computation benefit because computation of the EIBBO requires running a Kalman smoother that iterates through all time points. Instead, we mini-batch by subsampling B3 spatial points
N B
Ba Σ Bg
ELL 2
[lgg(YY Lts [f,))
(20)
where i is uniformly sampled. We used in conjunction with PHYSS-SVGP and PHYSS-SVGP: this results in dominant costs O(N, (M, d, d)²) and O (N, d, .M, d) when B, < N..
5 Handling the PSD Constraint
As discussed in Sec. 3.2 when the differential equation is non-linear, the model is no longer conjugate and the resulting natural gradients are not guaranteed to result in p.s.d updates. This issue has received some attention in the literature [53, 64, 39), but these approaches do not maintain an efficient conjugate representation. One distinction [72], which uses the Gauss-Newton approximation to maintain conjugate operations. We now extend this to support spatial inducing points and non-linear transformations. Due to space we focus on PHYSS-SVGP. but see App. A.5.3 for further details. The troublesome term for the natural gradient update in Eqn. (14) is the Jacobian of the ELL w.r.f. to the second expectation parameter; which is not guaranteed to be p.s.d unless the ELL is log convex [39]. Focusing at single location (t, s):
dELLa 0| as..
E p(t. (a.) log |logp(Y (Yn [fm)
we apply the Bonnet's and Price's 's theorem [38] bring the differential inside the expectation and make a Gauss-Newton [19] approximation ensuring that the Jacobian is p.s.d
ŌELL
N Σ n.P
8gn (Jaa) Oú,
dolog pl (Y, a gn)
E (()))))
Jn.p Hn.p Jn where Ja
Hn,
d? 9n
(21)
and In.p g (,) (Eqn. (4)) and Mn is the mean of p(f, ū.) (Eqn. (19)). When using spatial mini-batching Eqn. (21) also subsampled.
6
Related Work
From the optimality of natural gradients, in the conjugate setting, we exactly recover batch GP based models such as [68, 29, 4]. Our inference scheme also applies models that do not require derivative information i.e. in d, = d, = 1. As a special case, we recover [20], but we have extended the inference scheme support spatial mini-batching, allowing big spatial datasets be used. The linear weighting matrix can be used to define a linear model of coregionalisation and its variants [7, 77, 42, 66] and through appropriately designed functionals also non-linear variants [73].
In Álvarez et al. [76] GP priors over the solution of differential equations obtained through stochastic forcing term but they only consider situations where the Greens function available. In [22, 23, 57,33), efficient statesppace algorithms are derived but are limited to the temporal setting only. Similarly, Heinonen et al. [24], Iearn a 'free-form ODE'. In the spatio-temporal setting Krämer et al. [35] and Duffin et al. [14] (which builds [18]) derive extended Kalman filter algorithms. Additionally there are approaches to constraining GPS by linear differential equations [37, 1, SJ. More generally than [4] in [21] GP priors over the solutions to linear PDEs with constant coefficients are derived.
Beyond GP based models, physics informed neural networks (PINNs) incorporate physics by con- structing a loss function between the network and the differential equation at a finite set of collocation points [48]. This amounts to highly complex optimisation problem [36] bringing difficulties for training 770, 71] and uncertainty quantification (UQ) [16]. Current approaches to quantifying uncer- tainty in PINNS are based on dropout [75] and conformal predictions [47]. recent years UQ and deep learning has received much attention however is limited by its computational cost [45].
7

---

Table Test performance on the simulated damped pendulum. Time is the total wall clock time in seconds.
Table 2: Test performance on the magnetic field strength experiment. Results are computed w.r.f. to the first output. Time the average epoch time in seconds.
MODEL
WHITEN
C
TIME
RMSE
NLPD
112.2
0..2
TIME
SQUARED
1000
144.29
0.79
SPATIAL SIZE
10
20
5
20
10
153.55 52
0.35
0.41
100
195.25
0.2
-0.08
HEI MHOLFZ-GP
0.21
0.46
2.37
0.23
0.97
0.97
500
1011.88
0.33
0.32
MIYSS-GP
0.43
0.60
1.16
0.21
0.97
0.97
AUTOIP
1000
5134.13
0.36
0.41
PHYSS-SVGP
0.41
0.41
0.31
0.23
0.96
0.97
10
164.58
0.16
-0.30
HGYS
0.29
0.40
0.35
0.65
0.86
0.93
100
208.8 81
0.05
-0.41
0.29
0.28
0.16
0.65
0.61
0.74
500
1088.31
0.05
-0.75
1000
5656.62
0.05
-1.39
Scalar Potential (f)
Magnetic Field Strengh (H (H-field: fV
0.02
2
Spacce
0.00
0
-0.02
-2
2
3
'time'
'time"
Figure 2: Curl free synthetic example. The left panel displays the learnt scalar potential functions by PHYSS-GP with N₂ 20, and the right panel illustrates the associated vector field.
7 Experiments
We now examine the performance of our PHYSS-GP methods on multiple synthetic and real-world datasets. We compare against batch GP (no physical knowledge) and current state-of-art methods AUTOIP and HELMHOLTZ-GPP We provide more details on all experiments App. B.
Non-linear Damped Pendulum In this first synthetic example, we consider learning the non-linear dynamics of damped swinging pendulum. This is described by a second-order differential equation
d²0 dt2
d0 sin (0) + dt
(22)
where b > The first term is non-linear foreing term, and the third the damping term. With b 0.2 we simulate solution using Fuler's method [9] and generate 20 points in L. E [0, 6] for training and 200 E [6, 30) for testing, with additive Gaussian noise of variance 0.01.
We are interested in i) the effect of the number of collocation points, ii) the effect of the optimisation algorithm. To answer these questions, we compare against AUTOIP on [10, 100, 500, 1000] collocation points, with and without whitening, Results are tabulated in Table 1. As expected, the predictive RMSE of all models decreases as the number of collocation points increases. Due to the cubic complexity of AUTOIP, the total time significantly increases the number of collocation points increases. For example, when using 1000 collocation points, AUTOIP is z 39 times slower than PHYSS-GP. Interestingly, the un-whine on-teiteeed case performs poorly, possibly due to the nonlinearity of the differential equation making optimisation difficult. This indicates that either whitening or natural gradients are required to handle the non-linearity arising due to the differential equation.
Curl-free Magnetic Field Strength In this experiment, we consider modelling the magnetic field strength of dipole H(r) = - -vrr), where (r) m-r/t=/" is scalar potential function
[11]. Labelling the input dimensions as 'time", 'space' and 'z", we let m = [0, 1,0) and generate observations from spatio-temporal grid with N, = 50, and N. [5, 10, 20), at 1. H(r) is curl-free field and so we compare curl free part of HHOL against PHYSS-GP and its variants. HELMHOLTZ-GP and PHYSS-GP equivalent models (as this is the conjugate setting,

---

Table 3: Test performance on the diffusion-reaction system. Time the total wall clock time in seconds. PHYSS-EKS significantly outperforms all models, and due to the EKS prior only requires 1 epoch for inference. PHYSS-SVGP, achieves the same performance as AUTOIP over twice as fast.
MODEL
RMSE
NLPD
CRPS
TIME
EPOCHS
PHYSS-EKS BKS
0.09
-1.26
0.038
1.1 x 10
PHYSS-SVOP
0.19
0.093
1.4x 10
20000
PHYSSS-SVGF
0.17
0.077
10a
20000
AUTOIP
0.17
-0.29
0.063
1.1 X 10
20000
Theorem 3.1), and recover the same posterior and predictive distribution (up numerical precision). However, due to the cubic-in-time complexity HELMHOLTZ-GP, at larger spatial sizes, is over 2-times slower. The hicrarchical approximation substantially faster than PHYSS-GP and performs similarly. As expected when introducing sparsity both PHYSS-SVGP and PHYSS-SVGP, are even faster; however, this is compensated by slight drop in predictive performance. See Fig. and Table
PHYSS-EKS
AUTOIP
PHYSS-SVGPb
Variunce
0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Figure 3: Results the diffusion reaction system. The top row denotes the predictive mean, and the bottom the 95% confidence intervals. The white line denotes where the training data ends. Only PHYSS-EKS captures the sharp boundaries, due the IWP kernel. PHYSS-SVGPP recovers similar solution to AUTOIP but at half the computational cost.
Iiffusion-Reaction System
Consider diffusion-reaction system given by an Allen-Cahn equation d²u 0.00001
Ou Ot
dx + + 5u=
(23)
de (1, 1). Following [40]. we use the solution provided by [49] and sample 256 training examples from t [0,00.2 00.0.8 We compare PHYSS-EKS (where g is linerized the EKS prior). PHYSS-SVGP and PHYSS-SVGP against AUTOIP. Following [40].
where x E [-1, / [0, 1]. u(0, .) x² cos( x)), u(t, 1) = u(l,, 1) and 0u (1, -1)
learning rate of 0.001 for Adam. For AUTOIP, we place 100 collocation points across the whole input domain on a regular grid. For both PHYSS-SVGP, and PHYSS-SVGP,
require more collocation points in the temporal dimension and place them on a regular grid of size 20 x 10. For PHYSS-EKS we use an integrated Wiener kernel (IWP) time [56] and place 100 x 40 collocation points. We unable to place more collocation for AUTOIP due to computational limits. Results are presented in Fig. 3 and Table PHYSS-EKS requires only single epoch and better handle the sharp boundaries. Our method PHYSS-SVGP, is over twice as fast as AUTOIP whilst achieving similar predictive RMSE.
Ocean Currents We now model oceanic currents in the Gulf Mexico. We follow [4] and use the dataset provided by D' Asaro et al. [15] that has information from over 1, 000 buoys. We focus on region long. [-90,- 84.5), lat. [26, 30] on 2016-02-25. by computing hourly averages. This results in N = 42, 243 observations, and we construct a test-train split on 0.1 per cent of the data. infeasible to run HHL.TOTTTZ- due data size Berlinghieri al. [4], observations from only 19 buoys are used with N 55). However,
PHYSS-SVGP, with 50 spatial inducing

---

X10
2016-02-25 00:00:00
2016-02-25 20:00:00
3.4
Northings
3.2
3.0
-9.7
-9.4
Eastings
x 10
Fastings
x10
Figure 4: Predicted ocean currents by PHYSS-SVGP; True observations are in grey, and predictions in green. The thickness of the line represents uncertainty and computed by the L2 norm of the standard deviations across both outputs.
points and a spatial mini-batch size of 10, and plot results Fig. 4. Our predictions are in excellent agreement with the test data, achieving RMSE of 0.14, NLPD of -0.52, CRPS of 0.078, and an average run-time 1.86(s) per epoch.
8 Conclusion
introduced a physics-informed state-space GP that integrates observational data with physical knowledge. Within the variational inference framework, we derived computationally efficient algorithm that uses Kalman smoothing to achieve linear-in-time costs. To gain further computational speed-ups, proposed three approximations with inducing points, spatial mini batching and structured variational posteriors. When used in conjunction, they allow us to handle large-scale spatiotemporal problems. The bottleneck is always the state size, where nearest neighbours GPS [13, 74] could be explored. For highly non-linear problems, future directions could explore deep approaches [52] more flexible kernel families [69]. One limitation is the use of the collocation method which is only enforcing the differential equation point wise, whilst future work could look at the more general methods of weighted residuals [46].
Acknowledgments and Disclosure of Funding
OH acknowledges funding from The Alan Turing Institute PhD fellowship programme and the UKRI Turing Al Fellowship (EP/VO2678X/I), AS acknowledges support from the Research Council of Finland (339730). TD acknowledges support from UKRI Turing AI Acceleration Fellowship (EP/V02678X/1) and Turing Impact Award from the Alan Turing Institute. The authors acknowl- edges the University of Warwick Research Technology Platform (aquifer) for assistance in the research described in this paper. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) license any Author Accepted Manuscript version arising from this submission.
References
[1] C. G. Albert. Gaussian processes for data fulfilling linear differential equations. In Proceedings, volume 33, page 5. MDPI, 2019.
[2] S.-i. Amari. Natural gradient works efficiently learning. Neural Computation, 10(2):251-276 1998.
[3] G. B. Arfken, H. J. Weber, and F. E. Harris. Mathematical Methods for Physicists: Compre- hensive Guide. Academic Press, 2011.
[4] R. Berlinghieri, B. L. Trippe, D. R. Burt, R. Giordano, K. Srinivasan, T. Özgökmen, J. Xia, and Broderick, Gaussian processes at the Helm (holtz): more fluid model for ocean currents. arXiv preprint arXiv:2302.10366 2023,
10

---

[5] A. Besginow and M. Lange-Horr Hegermann. Constraining Gaussian processes to systems of linear ordinary differential equations. In Advances in Neural Information Processing Systems 35 (NeurIPS), pages
Curran Associates, Inc., 2022.
[6] Black and M. Scholes. The pricing of options and corporate liabilities. Journal of Political Economy, 81 (3):6337 1973.
[7] E. Bonilla, K. Chai, and C. Williams. Multi-task Gaussian process prediction. In Advances in Neural Information Processing Systems 20 (NIPS). Curran Associates, Inc., 2008
[8] D. Borthwick. Introduction to Partial Differential Equations. Universitext. Springer Interna- tional Publishing, 2017.
[9] J. Butcher. Numerical Methods for Ordinary Differential Equations. Wiley, 2016.
[10] P. E. Chang. J. J. Wilkinson, M. E. Khan, and A. Solin. Fast variational learning in state-space Gaussian process models. In 30th IEEE International Workshop on Machine Learning for Signal Processing (MLSP), pages 1 1-6. IEEE, 2020.
[11] T. Chow. Introduction to Electromagnetic Theory: A Modern Perspective. Jones and Bartlett Publishers, 2006.
[12] J. Cockayne, C. Oates, T. Sullivan, and M. Girolami. Probabilistic numerical methods for PDE-constrained Bayesian inverse problems. AIP Conference Proceedings. 1853(I), 06 2017.
[13] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand. Hierarchical nearest-neighhor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514):800-811 2016.
[14] C. Duffin, E. Cripps, Stemler, and M. Girolami. Low-rank statistical finite elements for scalable model-data synthesis. Journal of Computational Physics, 463:11 (88::11266) 1261, 2022.
[15] E. D' Asaro, C. Guigand, A. Haza, Huntley, G. Novelli, T. Özgökmen, and Ryan. La- grangian submesoscale experiment (LASER) surface drifters, interpolated to 15-minute intervals, 2017. URL https: / /ataa
org/data/R4. wtt.aaataaai.....aaaa/a///t. x265. 237: :0011.
[16] C. Edwards. Neural networks learn to speed up simulations. Communications of the ACM, 65: 27-29, 04 2022.
[17] D. Eriksson, K. Dong, Lee, D. Bindel, and A. G. Wilson. Scaling Gaussian process regression with derivatives. In Advances in Neural Information Processing Systems 31 (NeurIPS). Curran Associates, Inc., 2018
[18] M. Girolami, E. Febrianto, G. Yin, and F. Cirak. The statistical finite element method (statFEM) for coherent synthesis of observation data and model predictions. Computer Methods in Applied Mechanics and Engineering, 375:1 13533, 2021.
[19] G. H. Golub and V. Pereyra. The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. SIAM Journal on Numerical Analysis, 10(2):413-43) 1973.
[20] O. Hamelijnck, W. J. Wilkinson, N. A. Loppi, A. Solin, and T. Damoulas. Spatio-temporal vari- ational Gaussian processes. In Advances Neural Information Processing Systems (NeurIPS). Curran Associates, Inc., 2021.
[21] M. Harkonen, M. linear partial differential equations with constant coefficients. In International Conference on Machine Learning, pages 12587-12615. PMIR, 2023.
and B. Raita. Gaussian process priors for systems of
[22] J. Hartikainen and S. Särkkä. Sequential inference for latent force models. In Proceedings of the Twenty-Sevnnth Conference on Uncertainty Artificial Intelligence, UAI'11, page 311-318, Arlington, Virginia, USA, 201 L. AUAI Press.

---

[23] J. Hartikainen, M. Seppänen, and S. Särkkä. State-space inference for non-linear latent force models with application to satellite orbit prediction. In Proceedings of the 29th International Coference on International Conference on Machine Learning, ICML' page 723-730. Omni- press, 2012.
[24] M. Heinonen, Çagatay Yildiz, H. Mannerström, J. Intosalmi, and H. Lähdesmäki. Learning unknown ODE models with Gaussian processes. In International Conference on Machine Learning (ICML), 2018.
[25] P. Hennig. M. Osborne, and H. Kersting. Probabilistic Numerics: Computation Machine Learning Cambridge University Press, 2022.
[26] J. Hensman, M. Rattray, and N. D. Lawrence. Fast variational inference in the conjugate exponential family. In Advances in Neural Information Processing Systems (NIPS), pages 2888 2896. Curran Associates Inc., 2012.
[27] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence (UAI), pages 282-290. AUAI Press, 2013.
[28] M. D. Hoffman, D. M. Blei, C. Wang, and Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14, 2013.
[29] C. Jidling, N. Wahlström, A. Wills, and T. B. Schön. Linearly constrained Gaussian processes. In Advances in Neural Information Processing Systems 30 (NeurIPS). Curran Associates, Inc., 2017.
[30] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics- informed machine learning. Nature Reviews Physics, 3(6):422 440, 2021.
[31] M. E. Khan and W. Lin. Conjugate-computatee variational inference: Converting variational inference in non conjconju models to inferences in conjugate models. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), 2017.
[32] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXv preprint arXiv:1412.6980 2014.
[33] N. Krämer and P. Hennig. Linear-time probabilistic solution of boundary value problems. In Advances in Neural Information Processing Systems 34 (NeurIPS), pages 11160-11171. Curran Associates, Inc., 2021.
[34] N. Krämer and P. Hennig. Stable implementation of probabilistic ODE solvers. Journal of Machine Learning Research, 25(111):1 2024.
[35] N. Krämer, J. Schmidt, and P. Hennig. Probabilistic numerical method of lines for time- dependent partial differential equations. In Proceedings of the 25th International Conference on Artificial Intelligence and Statisties, volume 151 of Proceedings of Machine Learning Research, pages 625-639, PMLR, 2022.
[36] A. S. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, and M. W. Mahoney. Characterizing possible failure modes in physics-ifformed neural networks. Advances in Neural Information Processing Systems 34 (NeurIPS), 2021.
[37] M. Lange- Hegermmeg Algorithmic linearly constrained Gaussian processes. In Advances in Neural Information Processing Systems 31 (NeurIPS), pages 2137-2148. Curran Associates, Inc., 2018
[38] W. Lin, M. E. Khan, and M. Schmidt. Stein's lemma for the reparameterizztion trick with exponential family mixtures. arXiv preprint arXiv: 13398, 2019.
[39] W. Lin, M. Schmidt, and M. E. Khan. Handling the positive-definite constraint in the Bayesian learning rule. In Proceedings of the 37th International Conference on Machine Learning, Proceedings of Machine Learning Research. PMLR, 2020.
12

---

[40] D. Long, Z. Wang, A. Krishnapriyan, R. Kirby, S. Zhe, and M. Mahoney. AutolP: A united framework to integrate physics into Gaussian processes. Proceedings of the 39th International Conference Machine Learning, volume 162 of Proceedings of Machine Learning Research. PMI.R, 17-23 Jul 2022.
[41] A. Mchutchon and C. Rasmussen. Gaussian process training with input noise. In Advances in Neural Information Processing Systems 25 (NeurIPS). Curran Associates, Inc., 2011.
[42] P. Moreno-Muñoz, A. Artés, and M. Álvarez. Heterogeneous multi-output Gaussian process prediction. In Advances in Neural Information Processing Systems 31 (NeurIPS). Curran Associates, Inc., 2018.
[43] M. Padidar, X. Zhu, L. Huang, J. R. Gardner, and D. S. Bindel. Scaling Gaussian processes with derivative information using variational inference. Advances Neural Information Processing Systems (NeurIPS). Curran Associates, Inc., 2021.
[44] I. Pan, L. R. Mason, and O. K. Matar. Data-centric engineering: integrating simulation, machine learning and statistics. challenges and opportunities. Chemical Engineering Science, 249: 1 117271, 2022.
[45] T. Papamarkou, M. Skoularidou, K. Palla, L. Aitchison, J. Arbel, D. Dunson, M. Filippone, V. Fortuin, P. Hennig, J. M. Hemández-Lobato, A. Hubin, A. Immer, T. Karaletsos, M. E. Khan, A. Kristiadi, Y. Li, Mandt, Nemeth, M. A. Osborne, T. G. J. Rudner, D. Rügamer, Y. W. Teh, M. Welling, A. G. Wilson, and R. Zhang. Position: Bayesian deep learning is needed in the age of lcage-scale AI. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 39556-39586 PMLR, 2024.
[46] M. Pförtner, I. Steinwart, P. Hemnig, and J. Wenger. Physics-ifformed Gaussian process regression generalizes linear PDE solvers. arXiv preprint arXiv:2212.12474, 2022.
[47] L. Podina, M. T. Rad, and M. Kohandel. Conformalized physies-informed neural networks. In ICLR 2024 Workshop A14Differential Equations in Science, 2024. URL https:// openreview net/Eorum/id--ooWW8112
[48] M. Raissi and G. E. Karniadakis. Hidden physics models: Machine learning of nonlinear partial differential equations. Journal of Computational Physics, 2017.
[49] M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 178:686-707, 2019.
[50] C. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA, USA, 2006.
[51] J. Rilhimäki and A. Vehtari. Gaussian processes with monotonicity information. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS), volume Proceedings of Machine Learning Research. PMLR, 2010.
[52] H. Salimbeni and M. Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. In Advances in Neural Information Processing Systems 30 (NeurIPS), pages 4588- 4599, Curran Associates, Inc., 2017.
[53] H. Salimbeni, S. Fleftheriadis, and J. Hensman. Natural gradients in practice: Non-conjugate variational inference in gaussian process models. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTATS), 2018.
[54] S. Särkkä. Linear operators and stochastic partial differential equations in gaussian process regression. In Proceedings of the 21st International Conference on Artificial Neural Networks Volume Part II, ICANN' 11, page 151 1-1815 Springer-Verlas 2011.
[55] Särkkä and F. Careía-Fernández Temporal parallelization of bayesian smoothers, IEEE Transactions Automatic Control, 1021.
13

---

[56] S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Cambridge University Press, 2019.
[57] J. Schmidt, N. Krämer, and P. Hennig. A probabilistic state space model for joint inference from differential equations and data. In Advances in Neural Information Processing Systems 34 (NeurIPS), pages 12374-12385. Curran Associates, Inc., 2021.
[58] M. Schober, D. K. Duvenaud, and P. Hennig. Probabilistic ODE solvers with Runge- Kutta means. In Advances in Neural Information Processing Systems 27 (NeurIPS). Curran Associates, Inc., 2014.
[59] E. Solak, R. Murthy-smith, W. Leithead, D. Leith, and C. Rasmussen. Derivative observations in Gaussian process models of dynamic systems. In Advances in Neural Information Processing Systems 15 (NIPS). MIT Press, 2002.
[60] A. Solin. Stochastic Differential Equation Methods for Spatio-Temporal Gaussian Process Regression. Doctoral thesis, Aalto University, 2016.
[61] A. Solin, M. Kok, N. Wahlström, T. B. Schön, and S. Särkkä. Modeling and interpolation of the ambient magnetic field by Gaussian processes. Transactions on Robotics, 34(4):1122-117,, 2018.
[62] T. Stocker. Introduction to Climate Modelling. Springer Science & Business Media, 2011.
[63] S. Särkkä. Bayesian Filtering and Smoothing. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2013.
[64] M.-N. Tran, D. H. Nguyen, and D. Nguyen. Variational Bayes on manifolds. Statistics and Computing, 31:1-17, 2021.
[65] F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: new perspective. Statistics and Computing. 29: 1297-1315, 2018.
[66] J. Vanhatalo, M. Hartmann, and L. Veneranta. Additive Multivariate Gaussian Processes for Joint Species Distribution Modeling with Heterogeneous Data. Bayesian Analysis, 15:415- 447, 2020.
[67] C. Villacampa-Cale Calvo, B. Zaldívar, E. C. Garrido-Merchán, and D. Hernández-Lobato. Multi- class Gaussian process classification with noisy inputs. Journal of Machine Learning Research, 22(1), 2021.
[68] N. Wahlström, M. Kok, T. Schön, and F. Gustafsson. Modeling magnetic fields using Gaussian processes. In IEEE International Conference Acoustics, Speech and Signal Processing, pages 3522-3526, 2013.
[69] K. Wang, O. Hamelijnck, T. Damoulas, and M. Steel. Non-separable non-stationary random fields. In Proceedings of the 37th International Conference on Machine Learning (ICML), volume 119 of Proceedings of Muchine Learning Research, pages 9887-9897. PMLR, 2020.
[70] S. Wang, X. Yu, and P. Perdikaris. When and why PINNs fail to train: A neural tangent kernel perspective. Journal of Computational Physics, 449: 10768, 2022.
[71] S. Wang. S. Sankaran, and P. Perdikaris. Respecting causality for training physies-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 421:1168133 2024.
[72] W. J. Wilkinson, S. Särkkä, and A. Solin. Bayes-Newtor methods for approximate Bayesian inference with PSD guarantees. Journal of Machine Learning Research, 24(83):1-50, 2023.
[73] A. G. Wilson, D. A. Knowles, and Z. Ghahramani. Gaussian process regression networks. In Proceedings of the 29th International Coference on International Conference on Machine Learning, page 1139-1144 1199-116,111 2012.
14

---

[74] L. Wu, G. Pleiss, and J. P. Cunningham. Variational nearest neighbor Gaussian process. In Proceedings of the 39th International Conference on Machine Learning (ICML), volume 162 of Proceedings of Machine Learning Research, pages 24114-24130. PMLR, 2022.
[75] D. Zhang, L. Lu, L. Guo, and G. E. Karniadakis. Quantifying total uncertainty in physics- informed neural networks for solving forward and inverse stochastic problems. Journal of Computational Physics, 397:108850, 2019.
[76] M. Álvarez, D. Luengo, and N. D. Lawrence. Latent force models. In Proceedings of the Twelth International Conference Artificial Intelligence and Statistics (AISTATS), volume 5 of Proceedings of Machine Learning Research, pages 9-16. PMLR, 2009.
[77] M. A. Álvarez, L. Rosasco, and N. D. Lawrence. Kernels for vector-valued functions: A review. Foundations and Trends in Machine Learning, 4(3):195-260 2012.
15

---

Appendices
A Variational Approximation Derivation
A.1 Overview of Notation
Table 4: Table of Notation
Symbol
Size
Description
N
Number of observations.
O F N.
Number of latent functions.
Number of latent oulputs.
Number of input features.
Number of spatial points.
N
Number of temporal points.
d,
Number of spatial derivatives. Number of temporal derivatives.
de D de
de
Total number of temporall derivatives. State dimension.
d
B
Spatial batch size.
M.
Number of spatial inducing points.
X
N X F
Input data matrix.
X
N. X F
Spatial Locations of training data.
X
N
Temporal locations of training data.
x, X
X1.
F
Single training input.
X:
Temporal axis of single training input location x. Spatial axes of a single training input location x.. Output data matrix. Single training output. Filtering state. Mixing matrix between Q latent GPs Random vector of the D derivatives at location X,, Output linearly mixed GPS.
X
F-1 NxP P
Y
Y Y
f
N. d
W
(P D) X (QQD)
(X..)
D (P
F
R P.D.
Differential equation defined using D spatio-temporal derivatives and P outputs/states.
Ze
((F-F1
Spatial Inducing Points.
KX, X.)
N. X N.
Spatial Kernel.
X.)
N X N
Temporal Kemel.
K)X,X)
N
Temporal Kernel.
K DK(X,X) X)
(N D) (N D)
Spatio-tempora keree over all N locations and D derivatives.
(X,,, (.),
(N, de) (N. x de)
Gram matrix over temporal derivatives.
(N, x d. (N. d.)
Gram matrix over spatial derivatives.
A.2 Layout Vectors and Matrices
We use numerator layout for derivatives. Let Q denote the number of independent latent functions and the number of derivatives computed, and let fa.a denote the latent GP for d'th derivative the th latent function. We will need to keep track of the permutation of our data w.r.f. to space, time, and latent functions. Inspired by 'row-major' and 'column-major" layouts, we will use the following terminology that describes the ordering of the data across latent functions and time and space:
latent-data: F Fia [ .... (X)) with F, (X) [C,, (X). (X)]
which is ordered by stacking each of the latent functions on top of each other.
data-latent: Fai F (X.))... (X(X))) ) which is ordered by stacking the latent functions evaluated at each data point across all data points.
time-space: [f (x) which is ordered by stacking each the input points at each time
point on top of each other. This is only applicable when there is single latent function.
latent-times pient-time-sppce: F (x())),
......F...
(x (
Fo((X((1))
(st) N.
16

---

time-latent-space: [F (x.
The default order will be latent-data (and latent-time-space for spatio-temporal problems). Since all of these are just simple permutations of each other, there exists permutation matrix that permutes between any of the layouts above. We the function denote function that performs this permutation such that:
F Ald di (Fw)
(24)
Fia -dl-id (Fa)
A.3 Timeseries Setting- Single Latent Function
Let X RNXD, E Nx be input-output observations across P outputs, where N = Nv. For now, we only consider the case where Q = I. We assume that f has a state-space representation, and
denote its state with its D time derivatives F(X) [f(X), a e ... in latent-data format. The vector F(X) is of dimension (N x D). We also use the notation F,, = F(X,), which D-dimensional vector of the derivatives at location X,. The joint model is
a(x)) 8*f(X)
N Fr,
| DE)
(FF).
(25)
At this point, place no particular restriction the form of the likelihood, aside from decomposing across N. The prior p(F) is a multivariate GP of dimension N x D
FFF)
N 0, D D x)
(26)
Let (F) be a free-form multivariate Gaussian of the same dimension as FFF) then the corresponding ELBO is:
P
F, (D)
L = E (F)
log
(F)
E (F) [log p(
F, DE)] - KL [F(() [[F(
(27)
N
E ((F.) [logp(Y, Y, Fn, DE)] - (lp(F),
and the marginal ((F,)) w) D-dimensional Gaussian corresponding to the n n'th observation. The natural gradients are
ØELL 0j
A- (1 33) 8
Surrogate likelihood update
(28)
Surrogate model update
(29)
where Ā = [ [7]2] and [7]1 is an (N x D) vector and [7]2 an (N X D) X (N x D) matrix.
Eqn. (29) is a sum natural parameters, and so is the conjugate Bayesian update. Naively computing this would yield no computation speed up the computation cost would be cubic O(N.). However, the natural parameters of the likelihood (A) are guaranteed to be block diagonal, one block per data point (if X is initialised as so). This immediately implies that Eqn. (29) can be computed using efficient Kalman filter and smoothing algorithms. The structure of depends on the gradient of the expected log-likelihood BCLL Expanding this out By.
ÖELL 0H422
N
DE ((F.) [log (Y, F,, DE)]
N Σ 70
Σ
(30)
0[p4|e
where each component f., is (N x D) x (N x D) matrix that only has D x D non-zero entries; these the only elements that directly affect F.. Collecting all these submatrices into block diagonal matrix, we have matrix in data-latent format, however, ØELL is in latent-data, and so all we dp
need to do is permute by P:
ØELL
mal-id (blkdiag 1o F 1)]
(31)
0μ
17

---

Converting from natural to moment paramerisation the surrogate update is:
((F)
x N Y F, p(F)
N
(32)
= II
Y
Fn, V.. P(F)
where Y, is a D-dimensional vector, and V, is D) x D matrix, and efficient Kalman filtering and smoothing algorithms can be used to compute the surrogate model update. Substituting ((F) back into the ELBO it further simplifies:
p(Y Y F. DE) F P Y V))
log
((F)
N
Y
F.
pFF
N
E (F,,)) [log (Y,, F.. DE)]
Σ Eg v(F
log N YYn Fr, Vn)] log log V)
(33)
each term be computed efficiently as the by-product of the Kalman filtering and smoothing algorithm used to compute (F)).
A.4 Timeseries Setting- Multiple Latent Functions
We now generalise the previous section to handle multiple independent latent functions, i.e. Q > 0. The model prior now has the form
(F) I[
(34)
where F(F,) is a prior over f and its D partial deriatives. We consider two approaches: mean-field approximate posterior and full Gaussian.
The first approach defined mean-field approximate posterior (F) 11, q(F,) where each (F is
free-form Gaussian of dimension (N D). The natural gradient updates are now simply applied each component (F,) q) separately, and we essentially follow the update set out in App. A.3.
The second approach is full-Gaussian approximate posterior where ((F) is (Q D x N)- dimensional free-form Gaussian. In this case the ELI. is
N ELL Σ Eq(F.) Mogg((Y, F,, DE)
(35)
where F(F,)) is of dimension (Q x D). This implies that the gradient of the ELL BELL Blule Σ pn
where in now has (Q x D) x (Q D) non-zero entries. Switching to moment parameterisation
(F) x II Y F., V..VV..
| p(F)
(36)
where Y, is of dimension (Q x D) and V, is (Q D) (Q D). We can still use space algorithms by simply stacking the states corresponding each F, [56].
A.5
Spatio-temporal Data Single Latent Function
We now turn to the spatio-temporal setting where X, are spatio-temporal observations on a spatio- temporal grid ordered in time-space format. We now derive the conjugate variational algorithm for PHYSS-SVGP and PHYSS-SVGPy- The algorithms for PHYSS-GP and PHYSS-VGP, are recovered as special cases when Z = X.-
18

---

A.5.1 Spatial Derivative Inducing Points
We follow the standard sparse variational GP procedure and augment that prior with inducing points U D u at locations Z. We require that the inducing points are defined on spatial-temport grid at the same temporal points as the data X, such that Z = [Z.] N This is required to ensure Kronecker structure between the inducing points and the data. The joint model is
p(Y F) P(F ( UU)
(37)
where
P(U) N (U O, KP (X, X,) (2,
(38)
P [ N MFjU, EFIU)
and the conditional mean and covariance are given by
/F U K X, X.) 8 K. , 2,) K (X,, K2 Z,) -
(39)
I Z.) (K (Z,, Z.)) !|
and
EF = (X,, x,) X x
(40)
(X,, X,) @ (Z 2.)]
[ [ T
where K X = KP (X,, (X,, X,)
(x., Z+) which simplifies to
EFIU= TPjU=K?(XX) (, X,) @
(Z. K (X X.) K2 (x,, Z.) (Z,, (2 X.)
(41)
Due to the Kronecker structure, the marginal at time only depends on the inducing points in that time slice so we can still get CVI-style update that can be computed using a state-space model. To see why we again look at the Jacobian of the ELL: SELL
Σn fl, where ft, now has (D x M) x (D x M) non-zero entries, which corresponding to needed all M spatial inducing points with there derivatives predict at single time point. This is similar to the time series setting, except we have now predicted in space to compute marginals of (F). To be complete, we write that the marginal (U)) is
(UU) & (V. F., V.
P(F)
(42)
where Y. and F, are vectors of dimension (D x M), and V, is a matrix of dimension (D x M) x (D x M). The marginals (U,,, and the corresponding marginal likelihood p(Y V) can be computed by running a Kalman filter and smoother O(N, (M, . dg d)8). The marginal ((F) = N(F /F: Ep) where
18 (X,, Z,) (K? Z,,
(Z,, Z,)) m
(43)
and
Σp EFjU+ KP 4,) (K? (Zs, Z,)) 2] S
[ T K. (X 2,) (Z, 2,) "] (44)
A.5.2
Structured Approximate Posterior With Spatial Inducing Points
We now derive the algorithm for the case of the structured approximate posterior with spatial inducing points. The key is to define the free-form approximate posterior over the inducing points and their temporal derivatives and then use the model conditional compute the spatial derivatives. The model is
((Y|F) p),) u).
(45)
Each term is
p(D, u) = N O, D. Z)) F(F D. u) = N(F / Uo EFIU.)
(46)
where
/FIU. K.P XX, X, X,, ] K. 2 (X,, x.) (2, (2,)
D,
(47)
[
X X.) K, ZZ, - D. u
19

---

and
EFIU. = (X,, X1) & K. (X,, X.)]
(48)
-[Kx,z x K (X,, X,) @ K, (Z,, Z.)
[x 1T
where K
X X.
-KP x,) KP¹ (Xg: 2ex) and
K2 Z,, De Do 2o)
K P (X,, X,)
(49)
The approximate posterior is defined as
(F, D. u) FFF Di u) q(D. u)
(50)
where (D, u is a free-form Gaussian of dimension (Nd x N: x M,). The rest of the derivation simply follows App. 1. by simpling substituting K2 (X.. Zez) into the corresponding conditionals.
The final result is that the approximate posterior decomposes as
N. (U) x II N ( (I,(X.1FL, Fr, V. ) D(F)
(51)
where Y, and F, are vectors of dimension (d, x M,), and V. is a matrix of dimension (de x M,) x (dg x M,). The marginals q( (U)), and the corresponding marginal likelihood p( V) can be computed by running Kalman filter and smoother in O(N, (M. - dd)),, which compared to App. A.5.1 is not cubic in the number of spatial derivatives.
A.5.3 Gauss-Newton Natural Gradient Approximation
We now provide the full derivation of the Gauss-Newton approximation of the natural gradient used to ensure p.s.d updates. We will make use of the following identities, known the Bonnet and Price theorems (see, [38]).
0
EF (f E)
[E(f)] E 19tfl
f(f) Of
(52)
0μ
02
E g (( [f({)]
2 IE o(F (1
or orT E(T)
(53)
which describes how to bring derivatives inside expectations. To ease notations, we work with more general description of the model presented in the main paper, where have multiple independent latent functions and use Tp to denote -speeiii---pccfiiccc functions which, for example, can be used to represent DE or as the identity of standard Gaussian likelihoods. The model is
p(u,) = N ( [u 0, Kg)
p(f, mg) = N (f, Mr, uu; Dr, |u,)
(54)
Y n.g p( pYY,,,
(Y,.. n.g (fn. fn.
where the shapes are Ug E RM fg N Tp Rº -R.. R...Y. E NXP, and Ynp E R. The variational approximation is
(U) = NU m, S)
(55)
where U ju,,... uo]. m E QM and S E RQMXQM Let F (fi....,fo) fo). The expected
log-likelihood of the variational approximation is
ELL = E ((U)
E EMF UU)
Σ u.P
(Y ap (T
Σ E ((U) [EPE, U) [logr( Y n. 1T n,
(56)
E ((U)) [EMF_ [logp(Y p YY,,) n. T3 []
n.p
20

---

where & = t(n) is the time period associated with data X,, Here U3 are the spatial inducing points time t(n) and hence U RQM, We need to compute
ÖELL Øs
Σ
E ((2))
[ p(F. UL) [logp(Yn,)
aS
Σ n.P
Pk
E. (UU))
[PF, UW) [logp(Y (p((F,)) (p,]]) Pk
OSk
0²
Σ n.p
Pk E ((()))
E PFF, U2) loggp(y m.p [T (P(F),)))]
Pk
aU OU
0²
Pk E ((2))
logp(Y Yn (P(FF
Pk
aU 0U
n.P
Gauss-Newton
2
Σ
Pk ((+))
OT, (Tp n.p)
8'log p a.p
Tp)
TTT (F n.p)
P
OU,
OT, OT,
OU
R.P
Σ Pk E {(.)
J H&J
P
niP
(57)
where the shapes are
E QM, x 1 H₃ exx1
(58)
(59)
and Pxc is permutation matrix that permutes from data-latent to latent-data format. In implementa-
BELL but write it here tion, we do not need to perform this permutation as we only require the blocks 8S. for completeness.
A.5.4 Optimality of Natural Gradients In Linear Models
Theorem A.1. Consider a linear multi-task model of the form
fa(-) ~ GP(O, K,)
f(x) = [f)(
(60)
Y, = wf(X.) + where 1/ ~ N (0, blkdligg (1o2)p=))
then under a full Gaussian variational approximate posterior
( 2 N (f [f], mm
(61)
where m E R(NXQ)x1 S E (NX then the natural gradient update with a learning rate 1f recovers the optimal solution pf YY.
Proof. To prove this first derive the natural parameters of the posterior p(f Y). We then derive the closed form expression of the natural parameter update and show that they recover that of the posterior. Let
F(X) WF(X) ~ N [( W WT)
(62)
where Kx.x bikdiag
[K, (k1 and W w then the posterior (X Y) is Gaussian
p(F(X) Y) = N /FIY: EFJY)
(63)
with
MP = Wxxx W
WKxx WT + Y
(64)
Ep Y
WXxx W WT]
-
21

---

The covariance matrix can be simplified by invoking Woodbury's 's identity twice
EF [ [xx W
WKxx W
+ +
W Wxxx W
WI]
W xxxx W
=W Kxxx - Kx.x W
Wx,x W
WKxx
W
W W NoWwT
W
+ K
Kx
W
x.X
(65)
and the mean can be expressed
MF EFlY 'Y
(66)
W W
W
+ Kx X,X
W Φ YY
Now can immediately read off the posterior p(f transformed version
Y) (F() Y) = =p( Y) is simply
W Kx )1)
p(f
Y)
N
W
W
Xx
W W Φ b
W
(67)
N f "i Y
whose natural parameters are
W W
T-
W + W
(68)
2
X.X
We now derive the closed form expression of the natural gradient update with a learning rate of 1, and show that it recovers Aj Aiir. The expected log likelihood (ELL)
log N
Y Wf. $ })
ELL = E 9(t)
(69)
= log N
Y
Wf. +
Tr WSW
2
The required derivatives
ÖELL
2
(Y
W m)
- (
W m)
dm
3m
(70)
W Φ (Y
W m)
where the last follows because symmetric and
dELL
1 ô
Tr ②
wsw wSw
as
2 as
(71)
W
Φ
W
2
The natural gradient is now given
ØELL ÖUFLY
ØELL 8m
ØELL es
W $
Y
Wm) m) W -- W m
ØELL
W + W
8s
(72)
W - Y
W -w WT
The natural gradient update with learning rate of is
ÖELL 0pF
X ((F)
(73)
(
Y
where X X
-
the natural parameters of the prior pf)), hence after the update the
natural parameters are
W - Y W
(74)
9
W
1K-1
which recover those of pf Y), and hence recover the optimal posterior.
22

---

B Further Experimental Details and Results
EKS methods were run on CPUs. State-space methods running on GPU used the parallel form of the Kalman smoother (see [55, 20]).
B.1 An extension of AUTOIP
If one drops the requirement for state-space representations then the approximations proposed in Sec. 4 directly define approximations to the variational GP defined by Eqn. (13), and hence directly extend AUTOIP. For example on the non-linear damped pendulum in Sec. 7 run this extension of AUTOIP with whitening and 50 inducing points for C' = 1000 and achieve an RMSE of 0.06 + 0.001 and running time of 158.16 + 0.34, clearly improving the running time against AUTOIP. However the benefit of our methods is that PHYSS-GE GP remains linear in temporal dimensions which is vital for applications that are highly structured in time [20].
B.2 Modelling Unknown Physies
Modelling of missing physics can be handled by parameterising unknown terms with GPs. For example take a simple non-linear pendulum
d²0
sin 0.
(75)
dta
Now consider that the the sin(0) is unknown and we would like to learn it. we define the our differential equation in Eqn. (4)
d g dt
+ (2
(76)
where both fi(), f2(-) are latent GPS that wish to learn. We now construct 300 observations for training from the solution of Eqn. (75) across the range [0, 30) and 1000 for testing. We run PHYSS-GP and compare the similarity of the learnt latent GP f2(-) to the true function at the test locations and achieve an RMSE of 0.068 indicating we have recovered the latent force/unknown physics well.
B.3 Monotonic Timeseries
This first example showcases the effectiveness of PHYSS-GP GP in learning monotonic functions. Mono- tonicity information expressed by regularising the first derivative to be positive at set collocation points [51]:
Of 01
Of
p( = N f. o3) P(0)
4
where ((-) is a Gaussian cumulative distribution function, and v 1e - 1 a tuning parameter that controls the steepness of the step function. We plot predictive distributions of (batch) GP and PHYSS-GP in Fig. 5. The GP fits data and does not learn a monotonic function. However, using 300 collocation points, PHYSS-GP able to include the additional information and learn monotonic function whilst running 1.5 times faster.
B.4 Non-linear Damped Pendulum
All models were run using Nvidia Titan RTX GPU and an Intel Core i5 CPU. All were optimised for 1000 epochs using Adam [32] with learning rate of 0.01. Both the GP and AUTOIP had an RBF kemel (following Long al. [40]) and PHYSS-GP used a Matérn- 72; all with lengthscale of 1.0. The observation noise was initialised to 0.01 and the collocation 0.001. Both were fixed for the first 40% of training and then released. Predictive distribution of PHYSS-GP and AUTOIP plotted in Fig. 6.
B.5 Curl-free Magnetic Field Strength
All models were run using an Nvidia Titan RTX GPU and an Intel Core CPU. All models are run for 5000 epochs using Adam with learning rate of 0.01, and use Matérn-32 kernel on time,
23

---

100
75
50
25
0.0
0.2
0.1
0.6
0.8
1.0
Figure 5: Predictive distributions of and PIYSS-GG on the monotonic function in App. B. The cannot incorporate monotonicity information and fits the data.
PHYSS-GP- C-
2
00000000
-1
PHYSS-GP- C= 1000
AUTOIP- C= 1000
2-
10
25
30
Figure 6: Predictive distributions on the Damped Pendulum.
with ARD RBF kernels on the spatial dimensions, with a lengthscale 0.1 across all. The Gaussian likelihood is initialised with variance of 0.01 and held for 40% of training. All our methods used natural gradient learning rate of 1.0 as this is the conjugate setting.
B.6 Iiffusion-Reactionn System
We use data provided by [49] under an MIT license. All models were run using an Nvidia Titan RTX GPU and an Intel Core i5 CPU. Our method PHYSS-SVGP and PHYSS-SVGH use a Matérn 72 kernel on time and an RBF of space, both initialised with a lengthscale
0.1. We place the collocation points on regular grid of size 20 x 10 and use M, 20 spatial inducing points. We pretrain for 100 iterations using natural gradient learning rate of 0.01 and after use a learning rate 0.1 for remaining 19000 iterations. AUTOIP uses RBF kernel on both time and space with a lengthscale of 0.1. We place the collocation points on regular grid of size 10 learning rate 0.001 and train for a total of 20000 iterations.
All models use Adam with

---

B.7 Ocean Currents
Our method PHYSS-SVGP, was run using an Nvidia Titan RTX GPU and an Intel Core i5 CPU. We ran for 10000 iterations, using Adam with learning rate of 0.01. For natural gradients with used learning rate of 0.1. We used Matérn-3 kernel on time and RBF kernels on both spatial dimensions with lengthscales [24.0, 1.0, 1.0]. We used 100 spatial inducing points and a spatial mini-batch size
of 10.
25

---

NeurIPS Paper Checklist
The checklist is designed encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], , [No] or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the
relevant information is Not Available.
Please provide short (1-2 sentence) justificatin right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation.
While "[Yes] " is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] "provided proper justification is given (c.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used").
general, answering "[No] " or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS paper checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only the provided macros for your answers.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We have both theoretically (Theorem 1), and empirically demonstrated our claims (Sec. 7). We provide computational complexities for all proposed models showing
linear the temporal dimension performance.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
26

---

Justification: A discussion of two main limitations (collocation method, and spatial scaling) and future work provided in the conclusion.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "I "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: The assumptions of Theorem 3.1 are provided in Sec. 3 with full proof (correct to the best of our knowledge) provided in the appendix App. A.5.4.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: Experimental details, sufficient to reproduce the results, are provided both in the main paper Sec. and the appendix App. B. Code reproducing all methods and results provided in the supplementary material.
27

---

Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g.. a large language model), then there should either way to access this model for reproducing the results or way to reproduce the model (e.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Code for data downloading and processing for train test splits, as well as code reproducing all methods and results is provided in the supplementary material.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / /nipp.cco wwtpppbllcgggpiiieeeeeeeeeooeeee www.////e/e//eee///ee///e/eeee (oodeSummiisssss wwtpiiimssssooooooooioooooi for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
28

---

At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer:
[Yes]
Justification:
All required details are provided in App. B.
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core paper to a level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: We followed established published results for experimental setups, which did not include error bars.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: Specific details on the CPUs and GPUs used for experiments are provided in App. B.
Guidelines:
29

---

The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: tttps://neuripe: http://s.ecccpbllechcccccceeeeee
Answer: [Yes]
Justification: We have reviewed the NeurIPS Code of Ethics and we conform with this in every respect.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We discuss the benefits of physics informed machine learning in our introduc- tion. We do not envisage any negative societal implications.
Guidelines:
The answer NA means that there no socictal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
30

---

Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging. and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All existing methods and datasets used are open and cited throughout. Any licenses are explicitly stated in the appendix.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website). the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. conddttaeess has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors are encouraged reach out to the asset's 's creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We release our code under CC-BY 4.0.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
31

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
32


Reviews:
[{"review_id": "WOi6q0dNoM", "reviewer": "Reviewer_X3gM", "summary": "The submission discusses conditioning spatiotemporal Gaussian processes using differential equation constraints and observational data.\nSpecifically, the temporal component is handled by a Markovian prior (to achieve linear complexity), and the spatial component is dealt with by variational methods.\nAs such, the proposed algorithm extends the work by Hamelijnck et al. [18] to differential-equation constraints (and to handle spatial mini-batching); or, from the opposite perspective, it implements a version of a probabilistic numerical solver for spatiotemporal PDEs with variational methods (in space).\nThe resulting algorithm is evaluated on a damped pendulum ODE, a curl-free magnetic strength field, a one-dimensional reaction-diffusion PDE, and an ocean-current problem.", "questions": "None", "limitations": "Discussed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "All components are technically sound, and the submission fits into recent literature on dynamical systems, Gaussian processes, and physics-informed machine learning.\nOverall, I think this is a nice paper, and despite a few minor weaknesses (see \"Weaknesses\"), I recommend acceptance.", "weaknesses": "I identify two weaknesses:\n\n1. The proposed algorithm's scientific novelty as a combination of known techniques is limited. Hamelijnck et al. essentially cover spatiotemporal variational inference with a Markovian prior in the temporal dimension. Conditioning Gaussian processes on differential equation constraints via collocation has become standard practice in recent years. Berlingheri et al. developed the curl- and divergence-free spatial priors. The combination of these techniques is novel, but the increment to the existing literature is relatively small.\n2. The clarity of the submission's relation to existing approaches would improve by discussing more articles related to probabilistic numerical methods. Currently, the manuscript cites the PMM by Cockayne et al., the ODE (initial value problem) solvers by Schober et al., Krämer et al., and Tronarp et al., and the GP-PDE paper by Pförtner et al. (all references are in the submission). Beyond those papers, Krämer et al. (below) describe solving time-dependent, nonlinear PDEs with collocation and Markovian priors; Schmidt et al. (below) combine nonlinear collocation constraints with observational data (for ODEs); and Krämer and Hennig (below) solve boundary value problems with collocation. There are other related papers, but I believe the former articles should be discussed in the manuscript. Further, the literature on statistical finite element methods is related and should be acknowledged: see Duffin et al. (below). \n\n\nThe first weakness is more significant than the second one. I expect that the second weakness will be relatively straightforward to resolve. \nIn any case, I believe the strengths outweigh the weaknesses and lean towards recommending acceptance.\n\n\n**References:**\n\nSchmidt, Jonathan, Nicholas Krämer, and Philipp Hennig. \"A probabilistic state space model for joint inference from differential equations and data.\" Advances in Neural Information Processing Systems 34 (2021): 12374-12385.\n\nKrämer, Nicholas, Jonathan Schmidt, and Philipp Hennig. \"Probabilistic numerical method of lines for time-dependent partial differential equations.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2022.\n\nKrämer, Nicholas, and Philipp Hennig. \"Linear-time probabilistic solution of boundary value problems.\" Advances in Neural Information Processing Systems 34 (2021): 11160-11171.\n\nDuffin, Connor, et al. \"Statistical finite elements for misspecified models.\" Proceedings of the National Academy of Sciences 118.2 (2021): e2015006118.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive review. We address your two comments below.\n\n**W1: Novelty of the proposed approach.**  \nWhilst we agree that some special cases of our framework are known, the 'beauty and originality of the proposed framework is its unifying property' (reviewer **z1FQ**). Not only do we provide a variational state-space framework unifying AUTOIP, HELMHOLTZ-GP, PMM, and EKS (Thm. 1), but we additionally extend this framework in three ways\n    - Spatial-temporal derivative inducing points which are necessary for data sets that do not align to a grid and for reducing the cubic cost associated with the filtering state\n    - A structured VI approximation in which the state-space prior is only defined over the space-time points and temporal derivatives, which significantly reduces the size of the state\n    - Spatial minibatching which is necessary for large spatial data sets\n\nEach of these points is novel, and when used in conjunction alleviates many of the limitations of previous methods. For example, on our ocean currents experiment, running any competitor (AUTOIP, HELMHOLTZ, PMM) is infeasible due to the sheer size of the data set ($N=42243$). \n\t\n**W2: Relation to existing works.**  \nThank you for the additional literature, they are indeed related works, and we will include them in our literature review. We briefly reiterate their connection here:\n\n- Both Krämer and Hennig (2021) and Schmidt et al. (2021) are interesting works that are limited to the ODE setting. \n\n- Krämer et al. (2022): This work derives a spatio-temporal PDE solver that can be encapsulated in our framework by seeing it as an extension of the EKF prior described in Ex. 3.1. As such we do not see it as an alternative to our work but a method that can be combined with ours. This could lead to interesting directions where finite differences can be used in space whilst used in conjunction with our variational approximations enabling the application to large spatio-temporal datasets. We will add this discussion to our related work and future work section. \n\n- Duffin et al. (2021) build on the StatFEM work of Girolami et al. (2021) and Conrad et al. (2017) by deriving an extended Kalman filter where space is discretised through finite-elements. Like Krämer et al., this would be an interesting approach that could be used in conjunction with our work.", "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": "", "extra_scores": {}}, {"review_id": "R8kjDEfQhk", "reviewer": "Reviewer_z1FQ", "summary": "In this paper, the authors present a physics-informed Gaussian process based approach to learn the solution of ODE and PDE systems. In particular, they address the challenge of the cubic computational complexity with respect to the number of spatial observations. It is shown that multiple state-of-the-art approaches can be recovered as special cases of the presented approach, i.e., it is a unifying framework. With additional approximation techniques, the cubic spatial computational cost is reduced to linear. In multiple simulations, the proposed method achieves similar or better results compared to SOTA methods but with a significant reduced computational time.", "questions": "- I assume that there exist PINN based approaches with some form of uncertainty quantification. It would be interesting to compare the GP-based method to NN-based approaches (not necessarily for this paper but in general). Do you have any thoughts on that?\n- Is it possible to use the same approximations in AUTOIP to reduce the cubic computational cost? If so, is there still a benefit of the proposed method?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The proposed method demonstrates a significant improvement of the computational complexity with respect to the number of spatial observations compared to SOTA methods. \n- The beauty and originality of the proposed model is its unifying property, including existing approaches such as HelmholzGP and AUTOIP.\n- The techniques and ideas in the paper are clearly presented. I really enjoyed reading the paper.", "weaknesses": "- There are some typos such as line 103 “develop sate space algorithms for”, line 167 “complexity of O(N  (N_s” (missing t), and eq (17) Nt instead of N_t. Furthermore, I think that line 257  “RMSE of all models increases as the number of collocation points increases” make no sense as the RMSE is decreasing. Please check the paper carefully.\n- In the experimental section, some phenomena are rarely discussed. For instance, the RMSE of PHYSS-GP is increasing with 1000 collocation points.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and positive view of our work. We address the points you raised below.\n\n**W1: Typos.**  \nThank you for highlighting these typos, we have addressed them all. However, line 167 is *not a typo*; the computational complexity is $O(N (N_s \\, d_s \\, d)^3)$ because the expected log-likelihood decomposes across all spatio-temporal datapoints (as shown by the first term on Eqn 17), and hence is linear in $N$. This motivates the spatial-minibatching approximation where this is reduced to be $O(N_t \\, (N_s \\, d_s \\,d)^3)$. We will add this explanation to line 167.\n\t\n**W2: Increasing RMSE.**  \nThis is an artefact of the stochastic nature of the optimisation algorithms used for inference. We have rerun our Non-linear Damped Pendulum experiment across 5 datasets constructed with different random seeds. We report the results below.\n\n| Model    | C | Time |RMSE | NLPD |\n| -------- | ------- | ------- | ------- | ------- |\n| PHYSS-GP  | $10$ | \t$59.27 \\pm 38.524$ |\t$0.20 \\pm 0.007$ |\t$-0.28 \\pm 0.102$ |\n| | $100$ |\t$112.82 \\pm 0.400$ |\t$0.05 \\pm 0.001$ |\t$-0.44 \\pm 0.134$ |\n| | $500$ |\t$139.86 \\pm 0.541$ |\t$0.05 \\pm 0.003$ |\t$-0.79 \\pm 0.414$ |\n| | $1000$ |\t$122.94 \\pm 40.184$ |\t$0.05 \\pm 0.003$ |\t$-0.87 \\pm 0.327$ |\n|AUTOIP  | $10$ |\t$117.16 \\pm 27.378$ |\t$0.16 \\pm 0.001$ |\t$-0.41 \\pm 0.070$ |\n||\t$100$ \t|$154.10 \\pm 31.590$ |\t$0.05 \\pm 0.001$ |\t$-0.51 \\pm 0.129$ |\n|| \t$500$ |\t$1058.08 \\pm 17.497$ |\t$0.05 \\pm 0.001$ |\t$-0.88 \\pm 0.087$ |\n|| \t$1000$ |\t$5600.90 \\pm 36.231$ |\t$0.05 \\pm 0.001$  |\t$-1.39 \\pm 0.091$ |\n\n\n**Q1: Uncertainty in PINNs.**  \nPINNs have become a popular method for solving differential equations and amount to a highly complicated optimisation problem that requires specialised training regimes (see Wang et al., 2024). Current approaches to quantifying uncertainty (UQ) in PINNs are based on dropout (see Papamarkou et al.) and conformal predictions (see Podina et al., 2024). In recent years UQ for deep learning has received much attention however is limited by its computational cost (see Papamarkou et al., 2024). Exciting avenues of work could be to explore combinations of state-space algorithms and PINNs, to achieve linear time complexities but with the flexibility of PINNs. We will add this discussion to the related work section of the main paper.\n\n- Podina et al. (2024). Conformalized physics-informed neural networks. In *ICLR 2024 Workshop on AI4Differential Equations In Science*.\n- Papamarkou et al. (2018). Position: Bayesian deep learning is needed in the age of large-scale AI. arXiv preprint.\n- Zhang et al. (2018). Quantifying total uncertainty in physics-informed neural networks for solving forward and inverse stochastic problems. *Journal of Computational Physics*.\n- Wang et al. (2024). Respecting causality for training physics-informed neural networks. *Computer Methods in Applied Mechanics and Engineering*.\n\t\n**Q2: Application to AUTOIP.**  \nThe derivation of the state-space algorithms hinges on the approximate posterior being represented as a posterior with likelihoods/sites that decompose across time, which is guaranteed within the natural gradient framework (as shown in Eq. 16). AUTOIP has no such guarantees since it optimises the natural parameters in Euclidean space with optimizers like Adam, and uses a whitened representation. However, if one drops the requirements of deriving a state-space algorithm, all of the approaches for reducing the spatial computational complexity in Sec. 4 can equally be applied to AUTOIP. Indeed, this is simple to do within our codebase (that will be released on publication). For example, on the Non-linear damped pendulum, we run this extension of AUTOIP with whitening and 50 inducing points for $C=1000$ and achieve an RMSE of $0.06 \\pm 0.001$ and running time of $158.16 \\pm 0.34$, clearly improving the running time against AUTOIP. This is still slower than our methods as it is cubic in the number of inducing derivatives/points. We will add this example to the appendix and this discussion into the main paper. The benefit of our proposed approach is that we remain linear in the temporal dimension which is vital for applications that are highly structured in time.", "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": "", "extra_scores": {}}, {"review_id": "d3jBckyBYa", "reviewer": "Reviewer_4ngK", "summary": "The paper introduces a novel approach for solving partial differential equations (PDEs) using a Gaussian process prior. It leverages established methods for approximate inference and provides a cohesive framework that unifies existing results in probabilistic differential equation solving.", "questions": "- In the background section please define N and F. And the relations between al the other variables N_t, N_s x_{t,s}, y_{t,s}, d, etc. I find this very confusing – is there also a clash in notation?\n-\tIn Eq (3) can you elaborate if $\\mathbf f$ is multidimensional and how does it relate to $f$. This notation is again used in line 78.\n-\tIn line 66 is this definition correct? Somehow I have a hard time coming up with the correct dimensions.\n-\tIn line 70 $\\bar{\\mathbf f} $ is a time depdent function. But in line 78 it is a space time depdent function. Could you elaborate on this?\n-\tIn Line 87 is the non-linear differential operator only dependent on space derivatives, or can the also be higher order time derivatives?\n-\tEq (5): $\\bar f_q(\\mathbf X_n)$ should be a vector of dimensionality D, but in line 98 W is of dimensionality $PD \\times QD$. Can you elaborate on this.\n-\tIn eq(7) can you elaborate what $g(\\mathbf F_n)$ is? Is it an element-wise operation for each element in $\\mathbf F_n$?\n-\tIn line 265-266 can you please elaborate what $z$ is?", "limitations": "-", "rating": 4, "confidence": 2, "soundness": 1, "presentation": 1, "contribution": 2, "strengths": "The paper introduces a novel approach for solving partial differential equations (PDEs) using a Gaussian process prior. It leverages established methods for approximate inference and provides a cohesive framework that unifies existing results in probabilistic differential equation solving.", "weaknesses": "The paper’s presentation is challenging to follow. Despite my familiarity with continuous-time filtering, smoothing, Gaussian process (GP) literature, and continuous-time stochastic optimal control (where some of these partial differential equations (PDEs) appear), I struggle to determine the paper’s solidity. Numerous mathematical inaccuracies and inconsistencies hinder comprehension. For instance, the abrupt switch between notation for space variables (e.g., $s$ and $t$) and the combined space-time variable ($x$) lacks clear definition. While I recognize that machine learning sometimes necessitates complex notation, in this case, it impedes understanding.\n\nThe main model section exacerbates the issue. Equations 5 and 6 remain elusive; the dimensions of $Q$ and $P$ are unspecified. The definition of $F_n$ confuses me, as it prevents a sensible matrix product computation. Additionally, the collocation point definition remains unclear. Even in Example 3.1, the meaning of $g_k$ remains unknown to me. Consequently, assessing the paper’s solidity proves challenging.\nIn my view, this paper requires substantial revision before acceptance. Although I appreciate the effort invested, from a reviewer’s standpoint, the current presentation poses significant hurdles.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank you for your review and for raising questions that have helped improve the presentation of the paper. \n\nWe understand your concern related to the complicated notation. To aid understanding we will use the additional page in the camera-ready stage to include a notation section that will clarify the dimensions of all quantities. Additionally, we will provide a table in the appendix reporting the same. However, we would like to emphasize that throughout the paper we have attempted to use standard notation.\n\n**W1: Definition of Q/P.**  \nQ/P are scalars denoting the number of latent functions and outputs. Q as the number of latent functions is established notation in multi-output GPs (see, e.g., [67, 37]). Additionally in the App. A.3 we have already provided an expanded form which provides the dimensions of all the relevant quantities. We will explicitly state these dimensions in Sec. 3 of the main paper.\n\t\n**W2: Abrupt switch in s/t.**  \nWe are unsure what you mean by this. Throughout we have been clear that we are operating in the spatio-temporal setting (explicitly stated on line 58/59) and throughout 't' refers to time and 's' to space.\n\t\n**W3: Collocation point definition.**  \nThe collocation points are defined in Eq. (6). Intuitively we would want the function that we learn (which we have placed a GP prior over) to coincide exactly with the solution of the differential equation at hand, i.e. that the residual between them is zero. In practice, we can only ever enforce this at a finite set of locations, and these are called the collocation points. \n\t\n**W4: Definition of $g$.**  \nThe function $g$ is computing the residual described in [Collocation point definition] above. It measures the (pointwise) error between the current function and the solution to the differential equation. Equivalently if $f$ followed the differential equation exactly, then, by definition, its first time derivate would equal $N_\\theta(f)$ and $g(f)$ would be zero.\n\t\n\n**Q1: Definitions.**  \nAll quantities $N_s, N_t, x_{t,s}, y_{t, s}$ are defined on line 59. Again this is in accordance with the state-space GP literature (see [18]).\n\t\n**Q2: Multi-output.**  \nAs explained on line 79, Eq. (3) is a multi-output prior over a latent function f and its space/time derivates. Here we use established notation (see [47]) where a sample from $\\bar{f}$ is of dimension $N \\cdot D$ and corresponds to the multiple-outputs being stacked. We will clarify this further in the main paper. \n\t\n**Q3: Dimensions.**  \nYes, the equation is correct. As described on line 65/66 this is a $d$-dimensional vector (where $d$ is the number of time-derivates) and $\\bar{f}$ is the corresponding state vector. See for example [48].\n\t\n**Q4: Notation.**  \nThroughout we use $\\bar{f}$ to denote $f$ together with its spatial and/or temporal derivates. On lines 66-71 we highlight that standard state-space GPs construct a state over $f$ and its time derivates, and hence use the notation $\\bar{f}$. Specifically on line 66 we discuss timeseries models, and hence the state only depends on time. On lines 70/71 we discuss how spatio-temporal GPs can be represented as a state-space model, where now the state is modelling the temporal dynamics at each spatial point and hence the state is constructed over the spatial points (see [48, Sec. 12.5]). On line 78 we are constructing a GP prior over $\\bar{f}$ at all the input locations $X$, as defined in Eq. (3). We will further clarify this in Sec. 2 of the main paper.\n\t\n**Q5: Higher-order derivatives.**  \nThe non-linear operator can be dependent on any order of spatial or temporal derivatives. The only requirement is that the latent multi-output GP is defined over those derivatives as well (which will require sufficiently smooth temporal and spatial kernels). We will add this discussion to Sec. 3 to clarify.\n\t\n**Q6: Dimensionality.**  \nAs discussed in (1-1) above this follows standard multi-output GP notation. Here $f_q$ are (stacked) multi-ouput GPs. At a single input location a sample will have dimension $D$ (as defined on line 73). $F_n$ is defined by stacking all outputs from the $f_q$ GPs, which for a single input location will have dimension $Q \\times D$. This is then linearly mixed to create $P \\times D$ outputs. To make this clearer we will write the stacking explicitly so that $F_n (X_n) = W \\, \\text{vec}[f_1(X_n), \\cdots, f_Q(X_n)]^\\top$ and state all these dimensions in Sec. 3 of the main paper.\n\t\n**Q7: What is $g$ doing?**  \nThe function $g: \\mathbb{R}^{PD} \\to \\mathbb{R}$ computes the residual between the latent process and the differential equation (see W4. in **4ngK**). This follows the notation set out in [22, 42, 35]. \n\t\n**Q8: What is z?**  \nThe full magnetic field is defined in a 3-dimensional space on which we place our spatio-temporal model PHYSS-GP. Here 'z' is simply a label to denote the third dimension and can be thought of as 'depth'. We will be more explicit about this.", "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": "", "extra_scores": {}}, {"review_id": "da7pfejSi9", "reviewer": "Reviewer_zAxE", "summary": "The paper introduces a formulation for state-space GP which aims to capture the behaviour of PDE-based systems. The paper suggests multiple techniques to speed up the GP inference, including through decomposed kernel, variational methods on natural parameterisation of distribution, use of inducing points", "questions": "- Are there any intuition as to why the decomposition of the kernel is done between the spatial and the temporal variables? Is it for computational reasons or can it be linked to properties of physical phenomena?\n- Does (or can) the method perform in the case that the full form of the PDE is unknown? Or is the paper mostly focused on cases of just generating the PDE solution? \n- In the formulation given by Eq 6, the PDE and BC residuals are assumed to arise from a Gaussian distribution. How can they be related to the physical realisation of the PDE, or how can they be interpreted properly as a probability distribution? How realistic is this assumption? \n- Are there any additional assumptions in the form of g(f) should take for this method to work?", "limitations": "Limitations are adequately mentioned.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides good motivation and use cases for physics-informed GPs, and the need for scaling the inference process. The paper also provides adequate results and demonstrations of the GP for prediction and uncertainty quantifications, which well motivates the benefits of GPs over other solver types.", "weaknesses": "- The paper may be a bit hard to follow if one is not familiar with GPs already. For example, it may be difficult to immediately see how the PDE or the BV are incorporated into the probability distribution during the inference stage.\n- It may also be arguable that some parts of the paper can be seen as just application of standard sparse GP techniques and variational methods into the context of physics-informed GPs.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and positive view of our work. We address the points you raised below.\n\n**W1: Clarity.**  \nTo aid understanding we will provide an additional notation section and nomenclature table in the appendix as well as the proposed additions mentioned in the reply to reviewer 4ngK. To clarify, the boundary values are incorporated as observations (as defined in Eq. 6). At the inference stage this is incorporated through an additional likelihood. This means for  PHYSS-VGP the likelihood term in Eq. (13) can be written as $p(Y_n | \\bar{f}) = p(Y^\\mathcal{O} | H_\\mathcal{O} \\, F_n) \\, p(O | g(F_n)) \\, p(Y^\\mathcal{B} | H_\\mathcal{B} \\, F_n)$ where the differential equation is used to define $g$. We will add this discussion into Sec. 3.\n\n\n\n\t\n**W2: Novelty.**  \nPlease, see the response to reviewer X3gM W1.\n\n**Q1: Decomposition of kernels.**  \nThe primary motivation for the decomposition of the kernel into space and time is computation. This is a standard assumption for GPs and in this paper, we show that we can exploit this separability to derive state-space inference algorithms that are linear in the temporal dimension. One would expect that for linear and autonomous differential equations (time-invariant in this case) this will be a reasonable assumption, however, this is an interesting avenue for future research.\n\t\n**Q2: Unknown physics.**  \nThank you for this excellent question. This paper has focused on modelling solutions to differential equations; however, we can simply extend our framework (as in [35]) to model missing physics through GPs. This is similar to the latent force model of [66] except we now jointly model the latent force and the solution as GPs. To demonstrate this we construct 300 observations from a non-linear pendulum $\\frac{d^2 \\theta}{dt^2} + \\sin(\\theta) = 0$. We will consider the $\\sin(\\theta)$ term as missing, which we would like to learn. We now define $g=\\frac{d^2 \\, f_1}{dt^2} + f_2(t) = 0$. We achieve an RMSE of $0.068$ indicating we have recovered the latent force/unknown physics well. Going beyond this simple example is an exciting avenue for future work. We will add this discussion to the main paper with this example in the appendix.\n\t\n**Q3: PDE residual.**  \nThis is a common modelling assumption for collocation based methods. Ideally, we would model the residual as a noise-free observation (arising from a delta likelihood) and enforce it exactly. Treating the residual as arising from a Gaussian likelihood is a relaxation of this. We will add this discussion to the main paper. \n\n**Q4: Assumptions on $g$.**  \nIn general, for the probabilistic model to be well defined we require that g is measurable. To ensure that we are modelling something 'useful' we require that there be a unique well-defined solution to the differential equation. For ODEs, one can appeal to the Picard-Lindelöf theorem (see Thm. 36.4 in [22]) which gives general conditions for unique solutions to exist and requires that the non-linear operator be locally Lipsitz continuous. However, we are unaware of any results for general PDEs. Extending to problems that have multiple solutions would be an interesting avenue for future work. We will add this discussion to the main paper.", "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": "", "extra_scores": {}}, {"review_id": "pb7aA6oYMH", "reviewer": "Reviewer_BaEt", "summary": "The paper introduces a method for training spatio-temporal Gaussian processes which incorporate physics constraints including satisfying the governing equations at a number of collocation points and satisfying curl / divergence free constraints. Their approach scales both linearly in time and space by leveraging a number of approximations.", "questions": "- As I understand it, you can always choose to satisfy the differential equations at enough collocation points such that you overwhelm the data. For example, say I have very limited data so I choose the number of collocation points to be much greater than the number of data points. Clearly, if I use too many collocation points, my model will ignore the data and just follow the differential equations. How do you think about this perspective in the context of uncertainty quantification?\n- As a follow up to this previous point, what are the implications of assuming that $0_n^{(C)} = g(F_n) + \\epsilon_C$? Is it correct to think about this as implicitly assuming that the PDE is stochastic? \n- Are there situations where you believe the Gaussian process assumption could be limiting?", "limitations": "I think the authors have done a good job of identifying some potential limitations. I think the paper would greatly benefit from a more in depth discussion on using collocation points to enforce governing equation constraints and the limitations this perspective brings.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "- While many of the individual components (such as incorporating curl / divergence free constraints and satisfying governing equations at a number of collocation points) have been done before, I believe the combination of ideas here is novel. \n- Showing how to achieve linear complexity in the space and time dimension is valuable.\n- Showing how your approach recovers some prior works as a special case from a more general perspective is also valuable.\n- The bevy of test cases demonstrating the efficacy of your approach on both synthetic and real-world problems convincingly demonstrate the advantages offered by your approach.", "weaknesses": "- The problem statement is poorly motivated in my perspective. It would be helpful to be more specific about what you hope to achieve with hybrid modeling. For example, is your goal to improve the predictive accuracy of mechanistic models by incorporating data? What computational efficiency do you hope to improve? i.e. do you want to reduce the amount of data needed to train physics surrogates? Do you argue your approach is more efficient than solving the mechanistic model using traditional methods? etc. While I understand space is limited, being more specific about the problem statement will also help guide more targeted numerical studies in future works.\n- I think you need slightly more of a discussion on virtual observations (see questions below). When virtual observations are competing with actual data it would be helpful to discuss in detail how you think about $\\epsilon_C$ and choosing the number of collocation points.\n- Since you state that one of the goals of your approach is to be useful in uncertainty quantification, the numerical studies would have been strengthened by comparing more than just RMSE. For example, comparing the continuous ranked probability score (CRPS) would have given an indication of how well your approach is estimating uncertainty.\n- Minor:\n  - L94 space cases -> special cases?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review. We address the points you raised below.\n\n**W1: Problem statement.**  \nOur work contributes to the growing field of data-driven physics-informed models. These are hybrid models that aim to exploit physical inductive biases and data observations. Recently there have been principled ways to approach this problem through GPs, where the current state of the art includes AUTOIP and HELMHOLTZ-GP.  However, these approaches are cubic in the number of data points limiting them to small-scale studies. In this work we provide a unifying view of such GP approaches and derive state-space algorithms that are linear in the temporal dimension, bringing computational complexity closer to standard ODE solvers, and enable the application to large scale spatio-temporal problems. \n\t\n**W2: Virtual observations.**  \nPlease, see Q1 below.\n\t\n**W3: Uncertainty quantification.**  \nWe agree that providing uncertainty metrics is valuable. In Table 1, we have already provided negative log predictive distributions (NLPDs) however we will report them and CRPS for all experiments. For the Diffusion-Reaction experiment they are \n\n| Model    | RMSE | NLPD | CRPS |\n| -------- | ------- | ------- | ------- |\n| PHYSS-EKS  |  $0.06$  | $-1.26$ | 0.038 |\n| PHYSS-SVGP$_{H}$ | $0.17$ | $1.69$ | 0.077 |\n| AUTOIP    | $0.17$    | $-0.29$ | 0.065 |\n\n and for the Ocean Current experiment we achieve an NLPD of $-0.52$ and CRPS of $0.078$. \n\t\n**W4: Typo.**  \nThank you for pointing out this typo, we have addressed this now in the manuscript.\n\n**Q1: Uncertainty.**  \nIn this paper, we assume that the data observations and mechanistic descriptions align (although we can simply extend to handle missing physics, see review zAxE Q2). In this setting, we can view the uncertainty as representing our beliefs on 'how well we have solved the system'. Since we can only ever place a finite number of colocation points there will always be some error in our solution, and uncertainty quantification is vital in representing that. Moving forward having capable (both computationally and predictively) uncertainty-equipped models will be vital for handling misspecification (perhaps in prior or the likelihood) and for parameter learning.\n\n**Q2: Collocation points?**  \nPlease, see the discussion in reply to reviewer **zAxE** Q3. \n\t\n**Q3: GP as a limitation.**  \nUnder non-linear differential equations the true resulting posterior (of the generative model in Eqs. 5-6) will be non-Gaussian, and with highly non-linear equations (and or chaotic systems) approximating the true posterior with a Gaussian could lead to underestimation of the uncertainty (Turner et al., 2011). However, for many problems, the Gaussian assumption is effective as a Gaussian system is much simpler to solve [54,22].\n\n- Richard Turner and Maneesh Sahani (2011). Two problems with variational expectation maximisation for time-series models. *Bayesian Time Series Models.*", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank all five reviewers for their time and constructive reviews. This work *introduces a novel approach for solving partial differential equations (PDEs) using a Gaussian process prior* (**4ngK**) that *fits into recent literature on dynamical systems, Gaussian processes, and physics-informed machine learning* (**X3gM**), whose *beauty and originality of the proposed model is its unifying property* (**z1FQ**). We showcase our methods under a *bevy of test cases* (**BaEt**) that *convincingly demonstrate the advantages offered* (**BaEt**).

Based on the reviews, we have:
* Conducted additional experiments to showcase extensions of our proposed framework (unknown physics (**zAxE**, Q2) and extensions to AUTOIP (**z1FQ**, Q2)) and to provide additional uncertainty metrics (**BaEt**, W3).
* Improved the presentation by clarifying our problem statement (**BaEt**, W1) and our notation (see **4ngK**).
* Fixed the typos pointed out by the reviewers.

Meta-review: The paper presents a novel variational spatio-temporal Gaussian process (GP) approach for solving partial differential equations (PDEs) that efficiently handles both spatial and temporal dimensions with linear-in-time and reduced spatial complexity. This method, which incorporates physics constraints and variational techniques, outperforms existing state-of-the-art methods in both predictive accuracy and computational efficiency across various synthetic and real-world applications. Reviewers have generally praised the paper's strengths, and most concerns were addressed in the rebuttal phase. One remaining issue is readability; however, I found the paper understandable despite its technical nature. I concur with the four reviewers who recommend acceptance.

For the camera-ready version, address these two points: (1) clarify that Equation (4) represents a general form of only an evolution equation (and not a genreal PDE) and (2) note that the linear case can be handled without the need for potentials, as detailed in https://arxiv.org/abs/2212.14319.
Decision: Accept (poster)