Physics-Informed Variational State-Space Gaussian Processes

Authors: oliver  hamelijnck
arno  solin
theodoros  damoulas
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 23:50:09

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.

Main paper: 
 background on gaussian processes: 
 Gaussian processes A GP is a distribution on an infinite collection of random variables such that any finite subset is jointly Gaussian [50] . Given observations X ∈ R N ×F and y ∈ R N thenEQUATIONis a joint model where p(f | θ) is a zero mean GP prior with kernel K(•, •), f (X) ∼ p(f (X) | 0, K(X, X)), and θ are (hyper) parameters. We are primarily concerned with the spatiotemporal setting where we observe N t temporal and N s spatial observations x t,s ∈ R, y t,s ∈ R on a spatio-temporal grid. Under a Gaussian likelihood, all quantities for inference and training are available analytically and, naïvely, carry a dominant computational cost of O((N t N s ) 3 ). For time series data, an efficient way to construct a 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 a discretised linear time-invariant stochastic differential equation (LTI-SDE), which at time k isEQUATIONwhere A is a transition matrix, q k is Gaussian noise, H is an observation matrix, and f is a d-dimensional vector of temporal derivatives f = [f (•), ∂f (•) ∂x , ∂ 2 f (•) ∂x 2 , • • • ] ⊤ .With appropriately designed states, matrices and densities, SDEs of this form represent a large class of GP models, and Kalman smoothing enables inference in O(N t d 3 ), see [55] . In the spatio-temporal setting, when the kernel matrix decomposes as a Kronecker product K = K t ⊗ K s , then with a Markov time kernel, a state space form is admitted. This takes a particularly convenient form where the state is ft = [ f ((X s ) 1 , t), • • • , f ((X s ) N s , t)] ⊤ , and inference requires O(N t (N s d) 3 ), see [59] .Derivative Gaussian processes One main appeal of GPs is that they are closed under linear operators. Let D[•] = D t D s [•] be linear functional that computes D = d t d s space-time derivatives with D t [•] = •, ∂• ∂t , ∂ 2 • ∂t 2 , • • • and D s [•] = •, ∂• ∂s , ∂ 2 • ∂s 2 , • •• , then at a finite set of index points, the joint prior between f and its time and spatial derivatives isEQUATIONwhere f (X) = D f (X) and D * is the adjoint of D , meaning it operates on the second argument of the kernel [54] . When jointly modelling a single time and space derivative (d t = d s = 1) the latent functions are f = [f , ∂f ∂s , ∂f ∂t , ∂ 2 f ∂t ∂s ] ⊤ and the kernel isK = D K(X, X) D * =      K - - - ∂ ∂s K ∂ ∂s K ∂ ∂s ⊤ - - ∂ ∂t K ∂ ∂t K ∂ ∂s ⊤ ∂ ∂t K ∂ ∂t ⊤ - ∂ 2 ∂t ∂s K ∂ 2 ∂t ∂s K ∂ ∂s ⊤ ∂ 2 ∂t ∂s K ∂ ∂t ⊤ ∂ 2 ∂t ∂s K ∂ 2 ∂t ∂s ⊤     . This is a multi-output prior whose samples are paths of f with its corresponding derivatives. This prior is commonly known as a derivative GP and has found applications in monotonic GPs [51] , input-dependent noise [41, 66] and explicitly modelling derivatives [58, 17, 43] . State-space GPs can be employed in the temporal setting since the underlying state computes f (x) with its corresponding time derivatives. In Sec. 3.1, we extend this to the spatio-temporal setting.3 Physics-Informed State-Space Gaussian Processes (PHYSS-GP)We now propose a flexible generative model for incorporating information from both data observations and (non-linear) physical mechanics. We consider general non-linear evolution equations of the formEQUATIONwith appropriate boundary conditions, where f : R F → R is the latent quantity of interest and N θ is a non-linear differential operator [49] . We assume that g : R P •D → R 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 it 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 ∂f ∂t and the operator N θ f then the right hand side (0) are virtual observations. The PDE can now be observed at a finite set of locations known as collocation points. This is a soft constraint (i.e. f is not guaranteed to follow the differential equation), but it 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 modelF n = W • fq (X n ) ⊤ Linear Mixing, fq ∼ GP(0, Kq )EQUATIONEQUATIONwhere fq are derivative GPs (see Eqn. ( 3)) that are linearly mixed by W ∈ R (P D)×(Q D) , and We recover EKS style generative models (see Hennig et al. [25] ) when the mixing weight is identity W = I, and ϵ C , ϵ B → 0, and the non-linear transform g is linearised. Let the prior be MarkovY (O) , 0 (C)EQUATIONThis is now a form that can directly be implemented into an extended Kalman smoothing algorithm [62] . When Q > 1 the state f is constructed by stacking the individual states of each latent [55] . With linear ODEs EKS coincides with PMM. Example 3.2 (HELMHOLTZ-GP and Curl and Divergence-Free Vector Fields in 2D).Let v = [v t , v s1 , v s2] denote a 3D-vector field, then curl indicates the tendency of a vector field to rotate and divergence at a specific point indicates the tendency of the field to spread out. Curl and divergence-free fields followEQUATIONwhere∇ = [ ∂ ∂t , ∂ ∂s1 , ∂ ∂s2 ]. 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 [f 1 , f 2 ] be scalar potential functions thenEQUATIONdefine curl and divergence-free fields. In 2D this simplifies to using the grad and rot operators over v = [v s1 , v s2 ] (see [4] ). Placing GP priors over f q we incorporate this into Eqn. ( 6) by definingEQUATION)HELMHOLTZ-GP is defined as the sum of GP priors over 2D curl and divergence-free fields [4] . 
a spatio-temporal state-space prior: 
 The generative model in Eqn. (6) contains two 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 statespace 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 a 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 a repetitive structure that can be represented through a Kronecker product. The gram matrix isEQUATIONwhereK D• • = K • K • D * • D• K • D• K • D * • and D• [•] = (D • [•]) 1:excludes the underlying latent function.To find a Kronecker form of the gram matrix over X, we will exploit the fact that X is on a spatiotemporal grid and that the kernel is separable. Due to the separable structure a 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:EQUATION)The full prior over (a permuted) X is now given asp( f (X)) = N 0, K D t (X t , X t ) ⊗ K D s (X s , X s ). This is the form of a spatio-temporal Gaussian process with derivative kernels that can be immediately cast into a state-space form as in 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 a computational time of O(N t • (N s • d s • d) 3 ). 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 in Ex. 3.1 can now be simply extended to the PDE setting by placing colocation points on a spatio-temporal grid [35] .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 g in Eqn. (6) , and will also enable the reduction of the cubic spatial computation complexity in 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) withp(Y | f ) = N n p(y (O) n |H O F n ) p(0 (C) n |g(F n )) p(y (B)n |H B F n ). 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] :EQUATIONwhere we define the approximate posterior q(f | ξ) ≜ = N ( f | m, S ) as a free-form Gaussian with ξ = (m, S) and m ∈ R D N ×1 , S ∈ R D N ×D N . The aim is to represent the approximate posterior as a 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, 27] . For Gaussian approximate posteriors the natural gradient has a simple form [26] EQUATIONwhere λ = (S -1 m, 1 /2 S -1 ) and µ = (m, m m ⊤ + S) are the natural and expectation parameterisations. This is known as conjugate variational inference (CVI) as λ represent the natural parameters for the conjugate prior η [31, 10, 20, 71] . For now, we will assume that the likelihood is conjugate to ensure that [λ k ] 2 is p.s.d, this will be relaxed in Sec. 5. The derivative of the ELL isEQUATIONwhere the expectation is under q( f(t,s) ), a D dimensional Gaussian over the spatio-temporal derivatives at location x t,s . Within the sum, the only elements of [µ] 2 whose gradient will propagate through the expectation are the D × D elements corresponding to these locations. These points are unique and so ∂ELL ∂[µ]2 has some (permutated) block-diagonal structure, hence Eqn. ( 14) can be written as q( f )EQUATIONwhere Y t is D-dimensional. The natural gradient update, i.e. q( ft ) in moment parameterisation, can now be computed using Kalman smoothing inO(N t •(N s •d s •d) 3 ). Collecting Y = vec( [ Y t ] ), V = blkdiag [ V t ], then the ELBO can also be computed efficiently by substituting this form of q( ft ) inEQUATIONwhere the first two terms only depend on q(D f t ) and the final term is simply a by-product of running the Kalman filter, leading to a dominant computational complexity ofO(N • (N s • d s • d) 3 ).This cost is linear in the datapoints (N ) because the expected log likelihood above decomposes across all spatio-temporal locations. In summary we have shown that natural gradient is equivalent updating a block-diagonal likelihood that decomposes across time; hence the approximate posterior is computable via Kalman smoothing algorithms. Extending to multiple latent functions (Q > 1) we define a full Gaussian approximate posterior that captures all correlations between the latent functions q( f1 , ×Q) . 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 × N t as it encodes the correlations of spatial points and their spatio-temporal derivatives across the latent functions. We denote this model as PHYSS-VGP and PHYSS-EKS when using a EKS prior (see Ex. 3.1). Theorem 3.1. Let the approximate posterior be (full) Gaussian q( f1 , • • • , fQ )• • • , fQ ) ≜ = N f1 , • • • , fQ | m, S where m ∈ R (N ×Q)×1 , S ∈ R (N ×Q)×(N≜ = N f1 , • • • , fQ | m, S where m ∈ R (N ×Q)×1 , S ∈ R (N ×Q)×(N ×Q) . When g is linear a single natural gradient step with β = 1 recovers the optimal solution p( f1 , • • • , fQ | 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.1). 
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 N s . The second is a structured variational approximation that defines the approximate posterior only over the temporal prior and alleviates cubic costs associated with d s . Finally, we introduce spatial mini-batching that alleviates linear N s costs. When used in conjunction, the dominant computation cost is 3 . 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.O N t • d s • (M s • d t )In this first approximation, denoted by PHYSS-SVGP, we augment the full prior p( f ) with inducing points. By defining these inducing points on a spatio-temporal grid, we will show that we can still exploit Markov conjugate operations through natural gradients. Let ū = D u ∈ R M ×D be inducing points at locations Z ∈ R M ×F . From the standard SVGP formulation [27] , the ELBO isEQUATIONwhere q( f , ū) = p( f | ū) q(ū). By defining the inducing points on a spatio-temporal grid at temporal locationsEQUATIONand 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 s × D, reducing the computationalcomplexity to O(N (M s • d s • d) 3). For full details, see App. A.5.1.This second approximation, denoted as PHYSS-SVGP H , 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 asp(F, D t f ) = p(F | D t f ) p(D t f ) where p(F | D t f ) is a Gaussian conditional with mean E [F | D t f ] = I ⊗ K D s (X s , X s ) K s (Z s , Z s ) -1 D t f , with K D s (X s , X s ) = K s (X s , Z s ) D s K s (X s , Z s ) .We then define a structured variational posteriorq( f , D t f ) ≜ = p(F | D t f ) q(D t f ).Substituting this into the ELBO we see that all the terms with the prior spatial derivatives cancelE q log p(Y | f ) $ $ $ $ $ p( f | D t f ) p(D t f ) $ $ $ $ $ p(F | D t f ) q(D t f )Again, the marginal q(D t f ) maintains Kronecker structure, enabling Markov conjugate operations, leading to a computational cost ofO(N • d s • (N s • d) 3), see App. A.5.2. These variational approximations can simply be applied to non-state-space variational approximation, see App. B.1.Spatial Mini-Batching A standard approach for handling big data is through mini-batching where the ELL is approximated using only a data subsample [27] . Directly appling mini-batching would be of little computation benefit because computation of the ELBO requires running a Kalman smoother that iterates through all time points. Instead, we mini-batch by subsampling B s spatial pointsEQUATIONwhere i is uniformly sampled. We used in conjunction with PHYSS-SVGP and PHYSS-SVGP H , this results in dominant costs ofO(N t (M s • d s • d) 3 ) and O N t • d s • (M s • d) 3 when B s ≪ N s . 
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, 63, 39] , but these approaches do not maintain an efficient conjugate representation. One distinction is [71] , 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.t. to the second expectation parameter; which is not guaranteed to be p.s.d unless the ELL is log convex [39] .Focusing at a single location n = (t, s):∂ELL n ∂[µ k ] 2 = ∂ ∂S u E q(ūt) E p( fn| ūt) log p(Y n | fn )we apply the Bonnet's and Price's theorem [38] to bring the differential inside the expectation and make a Gauss-Newton [19] approximation ensuring that the Jacobian is p.s.dEQUATIONand g n,p = g(ū n ) (Eqn. ( 4)) and µ n is the mean of p( fn | ūt ) (Eqn. (19) ). When using spatial mini-batching Eqn. ( 21) is also subsampled. 
related work: 
 From the optimality of natural gradients, in the conjugate setting, we exactly recover batch GP based models such as [67, 29, 4] . Our inference scheme also applies to models that do not require derivative information i.e. in d t = d s = 1. As a special case, we recover [20] , but we have extended the inference scheme to support spatial mini-batching, allowing big spatial datasets to be used. The linear weighting matrix can be used to define a linear model of coregionalisation and its variants [7, 76, 42, 65] and through appropriately designed functionals also non-linear variants [72] .In Álvarez et al. [75] GP priors over the solution of differential equations are obtained through a stochastic forcing term but they only consider situations where the Greens function is available. In [22, 23, 56, 33] , efficient state-space algorithms are derived but are limited to the temporal setting only. Similarly, Heinonen et al. [24] , learn 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, 5] . 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 constructing a loss function between the network and the differential equation at a finite set of collocation points [48] . This amounts to a highly complex optimisation problem [36] bringing difficulties for training [69, 70] and uncertainty quantification (UQ) [16] . Current approaches to quantifying uncertainty in PINNs are based on dropout [74] and conformal predictions [47] . In recent years UQ and deep learning has received much attention however is limited by its computational cost [45] . Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, 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 not include experiments.• The experimental setting should be presented in the core of the paper to a level of 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, or as supplemental material. 
experiments: 
 We now examine the performance of our PHYSS-GP methods on multiple synthetic and real-world datasets. We compare against a batch GP (no physical knowledge) and current state-of-art methods AUTOIP and HELMHOLTZ-GP. We provide more details on all experiments in App. B.Non-linear Damped Pendulum In this first synthetic example, we consider learning the non-linear dynamics of a damped swinging pendulum. This is described by a second-order differential equationEQUATIONwhere b > 0. The first term is a non-linear forcing term, and the third is the damping term. With b = 0.2 we simulate a solution using Euler's method [9] and generate 20 points in t ∈ [0, 6] for training and 200 in t ∈ [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 as the number of collocation points increases. For example, when using 1000 collocation points, AUTOIP is ≈ 39 times slower than PHYSS-GP. Interestingly, the un-whitened 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.In this experiment, we consider modelling the magnetic field strength of a dipole H(r) = -∇ψ(r), where ψ(r) = m•r /|r| 3 is a scalar potential function [11] . Labelling the input dimensions as 'time', 'space' and 'z', we let m = [0, 1, 0] and generate observations from a spatio-temporal grid with N t = 50, and N s = [5, 10, 20] , at z = 1. H(r) is a curl-free field and so we compare the curl free part of HELMHOLTZ-GP against PHYSS-GP and its variants. HELMHOLTZ-GP and PHYSS-GP are equivalent models (as this is the conjugate setting, Theorem 3.1), and recover the same posterior and predictive distribution (up to numerical precision). However, due to the cubic-in-time complexity HELMHOLTZ-GP, at larger spatial sizes, is over 2-times slower. The hierarchical approximation is substantially faster than PHYSS-GP and performs similarly.As expected when introducing sparsity both PHYSS-SVGP and PHYSS-SVGP H are even faster; however, this is compensated by a slight drop in predictive performance. See Fig. 2 and Table 2 . Diffusion-Reaction System Consider a diffusion-reaction system given by an Allen-Cahn equationEQUATIONwhere x ∈ [-1, 1], t ∈ [0, 1], u(0, x) = x 2 cos(π x), u(t, -1) = u(t, 1) and ∂u ∂x (t, -1) = ∂u ∂x (t, 1). Following [40] , we use the solution provided by [49] and sample 256 training examples from t ∈ [0, 0.28]. We compare PHYSS-EKS (where g is linerized in the EKS prior), PHYSS-SVGP and PHYSS-SVGP H against AUTOIP. Following [40] , we use a 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 H we require more collocation points in the temporal dimension and place them on a regular grid of size 20 × 10. For PHYSS-EKS we use an integrated Wiener kernel (IWP) on time [55] and place 100 × 40 collocation points. We are unable to place more collocation for AUTOIP due to computational limits. Results are presented in Fig. 3 and Table 3 . PHYSS-EKS requires only a single epoch and can better handle the sharp boundaries. Our method PHYSS-SVGP H is over twice as fast as AUTOIP whilst achieving similar predictive RMSE.Ocean Currents We now model oceanic currents in the Gulf of 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 the region in 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. It is infeasible to run HELMHOLTZ-GP due to data size (in Berlinghieri et al. [4] , observations from only 19 buoys are used with N = 55). However, we run PHYSS-SVGP H with 50 spatial inducing points and a spatial mini-batch size of 10, and plot results in Fig. 4 . Our predictions are in excellent agreement with the test data, achieving an RMSE of 0.14, NLPD of -0.52, CRPS of 0.078, and an average run-time of 1.86(s) per epoch.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?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, confidence 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). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration 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 introduction. We do not envisage any negative societal implications. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or 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 a 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 is being used as 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 to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
conclusion: 
 We introduced a physics-informed state-space GP that integrates observational data with physical knowledge. Within the variational inference framework, we derived a computationally efficient algorithm that uses Kalman smoothing to achieve linear-in-time costs. To gain further computational speed-ups, we 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, 73] could be explored. For highly non-linear problems, future directions could explore deep approaches [52] or more flexible kernel families [68] . 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] . 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?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 nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main 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. 
: 
 Physical modelling is integral in modern science and engineering with applications from climate modelling [61] to options pricing [6] . Here, the key formalism to inject mechanistic physical knowledge are 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 in combining these traditional mechanistic models with data-driven methods through physics-informed machine learning. These hybrids approaches aim to improve predictive accuracy, computational efficiency by leveraging both physical inductive biases with observations [30, 44] .A principled way to incorporate prior physical knowledge is through Gaussian processes (GPs). GPs are stochastic processes and are a data-centric approach that facilitates the quantification of uncertainty. Recently AUTOIP was proposed in order to integrate non-linear physics into GPs [40] , where solutions to ordinary and partial differential equations (ODEs, PDEs) are observed at a 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 to curl and divergence-free constraints. Such properties are required for the successful modelling of electromagnetic fields [60] and ocean currents through the Helmholtz decomposition [4] . These approaches enable the incorporation of physics but incur a 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 (EKS) [64, 34] and recover popular ODE solvers whilst achieving linear-in-time complexity through state-space GPs [57, 25] .In this work we propose a unified physics informed state-space GP (PHYSS-GP) that is a probabilistic models where mechanistic/physics knowledge is incorporated as an inductive bias. We can handle both linear and non-linear PDEs and ODEs whilst maintaining linear-in-time computational efficiency. We additionally derive a state-space variational inference algorithm that further reduces the computational cost in the spatial dimension. We recover EKS, PMM, and HELMHOLTZ-GP as special cases, and outperform AUTOIP in terms of computational efficiency and predictive performance. In summary:1. We derive a state-space GP that can handle spatio-temporal derivatives with a 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 a linear-in-time complexity and recovering state-of-the-art methods such as EKS, PMM and HELMHOLTZ-GP. 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 outperform the current state-of-the-art methods AUTOIP and HELMHOLTZ-GP both in terms computational and predictive performance.Code to reproduce experiments is available at https://github.com/ohamelijnck/physs_gp. 
 
 Appendix: A.1 Overview of Notation Nt × Nt Temporal Kernel. K(X, X) N × N Spatio-Temporal Kernel. K = D K(X, X) D * (N • D) × (N • D) Spatio-We use a numerator layout for derivatives. Let Q denote the number of independent latent functions and D the number of derivatives computed, and let f q,d denote the latent GP for d'th derivative of the q'th latent function. We will need to keep track of the permutation of our data w.r.t. 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 = F ld = [F 1 (X), • • • , F Q (X)] with F q (X) = [f q,1 (X), • • • , f q,D (X)]which is ordered by stacking each of the latent functions on top of each other.• data-latent:F dl = [F 1 (X n ), • • • , F Q (X n )]N n which is ordered by stacking the latent functions evaluated at each data point across all data points.• time-space: [f (X (st) t )] Nt t which is ordered by stacking each of the input points at each time point on top of each other. This is only applicable when there is a single latent function.• latent-time-space:[F 1 (X (st) 1 ), • • • , F 1 (X (st) Nt ), • • • , F Q (X (st) 1 ), • • • , F Q (X (st) Nt )] • time-latent-space: [F 1 (X (st) t ), • • • , F Q (X (st) t )] Nt tThe 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 a permutation matrix that permutes between any two of the layouts above. We use the function π to denote a function that performs this permutation such that:EQUATIONA.3 Timeseries Setting -Single Latent Function Let X ∈ R N ×D , Y ∈ R N ×P be input-output observations across P outputs, where N = N t . For now, we only consider the case where Q = 1. We assume that f has a state-space representation, and we denote its state with its D time derivatives asF(X) = [f (X), ∂f (X) ∂X , ∂ 2 f (X) ∂X 2 , • • • ] in latent-data format. The vector F(X) is of dimension (N × D).We also use the notation F n = F(X n ), which is a D-dimensional vector of the derivatives at location X n . The joint model isEQUATIONAt this point, we place no particular restriction on the form of the likelihood, aside from decomposing across N . The prior p(F) is a multivariate GP of dimensionEQUATION) Let q(F) be a free-form multivariate Gaussian of the same dimension as p(F) then the corresponding ELBO is:EQUATIONand the marginal q(F n ) is a D-dimensional Gaussian corresponding to the n'th observation. The natural gradients areEQUATIONλ ← λ + η Surrogate model update (29) where 29) is a sum of natural parameters, and so is the conjugate Bayesian update. Naively computing this would yield no computation speed up as the computation cost would be cubic O(N 3 ). However, the natural parameters of the likelihood ( λ) are guaranteed to be block diagonal, one block per data point (if λ 0 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 ∂ELL ∂µ . Expanding this outλ = [ λ] 1 , [ λ] 2 ⊤ and [ λ] 1 is an (N × D) vector and [ λ] 2 an (N × D) × (N × D) matrix. Eqn. (EQUATIONwhere each component µ n is a (N × D) × (N × D) matrix that only has D × D non-zero entries; as these are the only elements that directly affect F n . Collecting all these submatrices into a block diagonal matrix, we have a matrix in data-latent format, however, ∂ELL ∂µ is in latent-data, and so all we need to do is permute by P:EQUATIONConverting from natural to moment paramerisation the surrogate update is:EQUATIONwhere Y n is a D-dimensional vector, and V n is a D × D matrix, and efficient Kalman filtering and smoothing algorithms can be used to compute the surrogate model update. Substituting q(F) back into the ELBO it further simplifies:EQUATION) each term can be computed efficiently as the by-product of the Kalman filtering and smoothing algorithm used to compute q(F).We now generalise the previous section to handle multiple independent latent functions, i.e. Q > 0. The model prior now has the formEQUATIONwhere p(F q ) is a prior over f q,1 and its D partial deriatives. We consider two approaches: a mean-field approximate posterior and a full Gaussian.The first approach defined mean-field approximate posterior q(F) ≜ = Q q q(F q ) where each q(F q ) is a free-form Gaussian of dimension (N × D). The natural gradient updates are now simply applied to each component q(F q ) separately, and we essentially follow the update set out in App. A.3.The second approach is a full-Gaussian approximate posterior where q(F) is a (Q × D × N )dimensional free-form Gaussian. In this case the ELL isEQUATIONwhereEQUATION)whereY n is of dimension (Q × D) and V n is (Q × D) × (Q × D).We can still use state-space algorithms by simply stacking the states corresponding to each F q [55] .We now turn to the spatio-temporal setting where X, Y are spatio-temporal observations on a spatiotemporal grid ordered in time-space format. We now derive the conjugate variational algorithm for PHYSS-SVGP and PHYSS-SVGP H . The algorithms for PHYSS-GP and PHYSS-VGP H are recovered as special cases when Z = X s .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 a spatial-temporal grid at the same temporal points as the data X, such that Z = [Z t ] Nt t . This is required to ensure Kronecker structure between the inducing points and the data. The joint model isEQUATIONwhereEQUATION)and the conditional mean and covariance are given byEQUATIONandEQUATION)whereEQUATIONDue to the Kronecker structure, the marginal at time t only depends on the inducing points in that time slice so we can still get a CVI-style update that can be computed using a state-space model. To see why we again look at the Jacobian of the ELL:∂ELL ∂[µ]2 =N n µ n where µ n now has (D × M ) × (D × M ) non-zero entries, which corresponding to needed all M spatial inducing points with there derivatives to predict at a single time point. This is similar to the time series setting, except we have now predicted in space to compute marginals of q(F). To be complete, we write that the marginal q(U) isEQUATIONwhere Y t and F t are vectors of dimension (D × M ), and V t is a matrix of dimension (D × M ) × (D × M ). The marginals q(U t ), and the corresponding marginal likelihood p( Y | V) can be computed by running a Kalman filter and smoother inEQUATIONandΣ F = Σ F | U + I ⊗ K D s (X s , Z s ) (K D s (Z s , Z s )) -1 S I ⊗ K D s (X s , Z s ) (K D s (Z s , Z s )) -1 ⊤ (44)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 to compute the spatial derivatives. The model isEQUATIONEQUATION)whereEQUATIONandEQUATION)The approximate posterior is defined asEQUATIONwhere q(D t u is a free-form Gaussian of dimension (N d × N t × M s ). The rest of the derivation simply follows App. A.5.1 by simpling substituting K D s (X s , Z ss ) into the corresponding conditionals. The final result is that the approximate posterior decomposes asEQUATIONwhere Y t and F t are vectors of dimension (d t × M s ), and V t is a matrix of dimension(d t × M s ) × (d t × M s ).The marginals q(U t ), and the corresponding marginal likelihood p( Y | V) can be computed by running a Kalman filter and smoother in O(N t • (M s • d) 3 ), which compared to App. A.5.1 is not cubic in the number of spatial derivatives.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 as the Bonnet and Price theorems (see, [38] ),EQUATIONEQUATIONwhich describes how to bring derivatives inside expectations. To ease notations, we work with a more general description of the model presented in the main paper, where we have multiple independent latent functions and use T p to denote likelihood-specific functions which, for example, can be used to represent DE or as the identity of standard Gaussian likelihoods. The model isEQUATIONwhere the shapes areEQUATIONwhereU = [u 1 , • • • , u Q ], m ∈ R QM ×1 and S ∈ R QM ×QM . Let F = [f 1 , . . . , f Q ].The expected log-likelihood of the variational approximation isEQUATIONwhere k = t(n) is the time period associated with data X n . Here U k are the spatial inducing points at time t(n) and hence U k ∈ R QMs . We need to computeEQUATION)where the shapes areEQUATIONand P k is a permutation matrix that permutes from data-latent to latent-data format. In implementation, we do not need to perform this permutation as we only require the blocks ∂ELL ∂S k but write it here for completeness.Theorem A.1. Consider a linear multi-task model of the formEQUATIONthen under a full Gaussian variational approximate posteriorEQUATIONwhere m ∈ R (N ×Q)×1 , S ∈ R (N ×Q)×(N ×Q) then the natural gradient update with a learning rate of 1 recovers the optimal solution p( f | Y).Proof. To prove this we 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. LetEQUATIONwhereEQUATIONwithEQUATIONThe covariance matrix can be simplified by invoking Woodbury's identity twiceΣ F | Y = W K X,X W ⊤ -W K X,X W ⊤ Φ + W K X,X W ⊤ -1 W K X,X W ⊤ = W K X,X -K X,X W ⊤ Φ + W K X,X W ⊤ -1 W K X,X W ⊤ = W W Φ -1 W ⊤ + K -1 X,X -1W ⊤ (65) and the mean can be expressed asEQUATIONNow we can immediately read off the posteriorEQUATIONwhose natural parameters areEQUATIONWe now derive the closed form expression of the natural gradient update with a learning rate of 1, and show that it recovers λ f | Y . The expected log likelihood (ELL) isEQUATIONThe required derivatives areEQUATIONwhere the last follows because Φ is symmetric andEQUATIONThe natural gradient is now given asEQUATIONThe natural gradient update with a learning rate of 1 isEQUATIONwhere λ p( f ) = 0, -1 2 K -1 ⊤ are the natural parameters of the prior p( f ), hence after the update the natural parameters areEQUATIONwhich recover those of p( f | Y), and hence we recover the optimal posterior.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 we 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-GP remains linear in temporal dimensions which is vital for applications that are highly structured in time [20] .Modelling of missing physics can be handled by parameterising unknown terms with GPs. For example take a simple non-linear pendulumEQUATION)Now consider that the the sin(θ) is unknown and we would like to learn it. If we define the our differential equation in Eqn. (4) asEQUATIONwhere both f 1 (•), f 2 (•) are latent GPs that we 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 f 2 (•) 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.This first example showcases the effectiveness of PHYSS-GP in learning monotonic functions. Monotonicity information is expressed by regularising the first derivative to be positive at a set of collocation points [51] :p(Y | f ) = N Y | f , σ 2 y , p(0 | ∂f ∂t ) = Φ( ∂f ∂t • 1 v )where Φ(•) is a Gaussian cumulative distribution function, and v = 1e -1 is 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 is able to include the additional information and learn a monotonic function whilst running 1.5 times faster. All models were run using an Nvidia Titan RTX GPU and an Intel Core i5 CPU. All were optimised for 1000 epochs using Adam [32] with a learning rate of 0.01. Both the GP and AUTOIP had an RBF kernel (following Long et al. [40] ) and PHYSS-GP used a Matérn-7 /2; all with a 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 are plotted in Fig. 6 .All models were run using an Nvidia Titan RTX GPU and an Intel Core i5 CPU. All models are run for 5000 epochs using Adam with a learning rate of 0.01, and use a Matérn-3 /2 kernel on time, with ARD RBF kernels on the spatial dimensions, with a lengthscale of 0.1 across all. The Gaussian likelihood is initialised with a variance of 0.01 and held for 40% of training. All our methods used a natural gradient learning rate of 1.0 as this is the conjugate setting.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-SVGP H use a Matérn 72 kernel on time and an RBF of space, both initialised with a lengthscale of 0.1. We place the collocation points on a regular grid of size 20 × 10 and use M s = 20 spatial inducing points. We pretrain for 100 iterations using a natural gradient learning rate of 0.01 and after use a learning rate 0.1 for the remaining 19000 iterations. AUTOIP uses a RBF kernel on both time and space with a lengthscale of 0.1. We place the collocation points on a regular grid of size 10 × 10. All models use Adam with a learning rate 0.001 and train for a total of 20000 iterations.Our method PHYSS-SVGP H was run using an Nvidia Titan RTX GPU and an Intel Core i5 CPU. We ran for 10000 iterations, using Adam with a learning rate of 0.01. For natural gradients with used a learning rate of 0.1. We used a Matérn-3 /2 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.• 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 a 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 a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to 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 to 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 a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an 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.In the case of closed-source models, it 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 to have some path to reproducing or verifying the results. • 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) is recommended, but including URLs to data and code is permitted. 


Decision: 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.

Reviews: review: paper_summary: paper_summary: The submission discusses conditioning spatiotemporal Gaussian processes using differential equation constraints and observational data.
Specifically, the temporal component is handled by a Markovian prior (to achieve linear complexity), and the spatial component is dealt with by variational methods.
As 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).
The 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.
strength_weakness: strength_weakness: I identify two weaknesses:

1. 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.
2. 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). 


The first weakness is more significant than the second one. I expect that the second weakness will be relatively straightforward to resolve. 
In any case, I believe the strengths outweigh the weaknesses and lean towards recommending acceptance.


**References:**

Schmidt, 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.

Krä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.

Krämer, Nicholas, and Philipp Hennig. "Linear-time probabilistic solution of boundary value problems." Advances in Neural Information Processing Systems 34 (2021): 11160-11171.

Duffin, Connor, et al. "Statistical finite elements for misspecified models." Proceedings of the National Academy of Sciences 118.2 (2021): e2015006118.
questions: questions: None
limitations: limitations: Discussed.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- In the experimental section, some phenomena are rarely discussed. For instance, the RMSE of PHYSS-GP is increasing with 1000 collocation points.
questions: 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?
- 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: limitations: Yes
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.

The 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.
In 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.
questions: 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?
-	In 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.
-	In line 66 is this definition correct? Somehow I have a hard time coming up with the correct dimensions.
-	In 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?
-	In Line 87 is the non-linear differential operator only dependent on space derivatives, or can the also be higher order time derivatives?
-	Eq (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.
-	In eq(7) can you elaborate what $g(\mathbf F_n)$ is? Is it an element-wise operation for each element in $\mathbf F_n$?
-	In line 265-266 can you please elaborate what $z$ is?
limitations: limitations: -
paper_summary: paper_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
strength_weakness: strength_weakness: - 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.
- 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.
questions: 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?
- 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? 
- 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? 
- Are there any additional assumptions in the form of g(f) should take for this method to work?
limitations: limitations: Limitations are adequately mentioned.
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.
- 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.
- 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.
- Minor:
  - L94 space cases -> special cases?
questions: 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?
- 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? 
- Are there situations where you believe the Gaussian process assumption could be limiting?
limitations: 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.

Comments: comment: Thank you for the clarification. My recommendation is to accept the paper.
Dear Reviewers,

Thank you for the reviews and for taking the time to respond to our rebuttal. If you have any additional comments or requests we would be happy to address them in the remaining reviewer-author discussion period, and if you believe we have already addressed your corrections/questions we would appreciate it if this could be reflected in your updated evaluation.

Many Thanks,
Authors
We thank you for your comment and your score increase.
Thank you for your comment and for clarifying W1.

Our proposed approach fits into building probabilistic models where mechanistic/physics knowledge is incorporated as an inductive bias to improve performance (such as the latent quantity of interest should follow a differential equation, or the vector field we are modelling should be curl-free). We agree that this allows one to use less data because the physics has been explicitly incorporated. However our approach can also be extended to the setting where there is partial physical knowledge, see discussion Q2 in reply to zAxE.  We do not see this work as an alternative to classical PDE solvers. We will add this discussion to our introduction and agree that it will help better place our work within the wider field of physics-informed ML.
I have read the rebutal and thank the authors for their answers. I do not have any further questions. In view of the other reviews and the promised improvement in clarification of the paper, I will increase my score by one point, though I am still of the opinion that the is work hard to follow.
Thank you for the reply!

For W1, I agree that the combination of the tools is new, even though each component may be known. My review assessed that the novelty of combining these tools may be somewhat limited, not that it doesn't exist. However, this perspective is subjective, which is why I recommend acceptance regardless.

For W2, thank you for the clarification. I agree with your perspective on those four papers.

Again, thank you for iterating. Since my assessed strengths and weaknesses remain, I will keep my already positive score.
Thank you for your thoughtful response.

Regarding the W1, let me try to rephrase. I agree that your unifying view of GP 
approaches for hybrid modeling is itself a valuable contribution. In my understanding 
of physics informed learning, it is often useful to understand the downstream tasks 
which you hope to achieve using a specific approach as this will guide algorithm design. 
For example, some approaches may seek to build predictive models 
which can learn from less data because they 
explicitly incorporate physics knowledge. Other approaches may seek to improve the 
accuracy of flawed mechanistic models using data. Other approaches may seek to 
use data to approximately solve differential equations more quickly 
with error estimates. Other approaches may seek to build cheap, approximate models 
of differential equations so that they can be used in optimization and control.
My point with W1 was that I think it would be very helpful to discuss where you 
think your approach fits into the broader context of physics informed ML.

Thank you for providing CRPS scores. I think this will be helpful for those 
hoping to extend/compare to your work in the future. 

Regarding Q1 (and this relates back to W1), thank you for clarifying. I think making it 
clear in the main text that the Gaussian likelihood is just a relaxation of a "delta" 
likelihood is helpful in understanding the motivations of your approach.