Adjoint-aided inference of Gaussian process driven differential equations
Authors: Paterne Gahungu
Christopher W Lanyon
Mauricio A Álvarez
Engineer Bainomugisha
Michael Thomas Smith
Richard David Wilkinson
Venue: NeurIPS 2022 Accept (2022)
Keywords: Uncertainty quantification
latent force
Gaussian process
differential equation
inverse problem

Abstract: Linear systems occur throughout engineering and the sciences, most notably as differential equations. In many cases the forcing function for the system is unknown, and interest lies in using noisy observations of the system to infer the forcing, as well as other unknown parameters. In differential equations, the forcing function is an unknown function of the independent variables (typically time and space), and can be modelled as a Gaussian process (GP). In this paper we show how the adjoint of a linear system can be used to efficiently infer forcing functions modelled as GPs, after using a truncated basis expansion of the GP kernel. We show how exact conjugate Bayesian inference for the truncated GP can be achieved, in many cases with substantially lower computation than would be required using MCMC methods. We demonstrate the approach on systems of both ordinary and partial differential equations, and show that the basis expansion approach approximates well the true forcing  with a modest number of basis vectors. Finally, we  show how to infer point estimates for the non-linear model parameters, such as the kernel length-scales, using Bayesian optimisation.

Adjoint-aided inference of Gaussian process driven differential equations
Paterne Gahungu"
Christopher W. Lanyon"
Quantum Leap Africa
Department of Computer Science
African Institute for Mathematical Sciences, Rwanda paterne gahungu@aims.. aw.
University of Sheffield, UK
c.w. w.w.lanyon@sheffielle ac uk
Mauricio A. Álvarez Department of Computer Science University of Manchester, UK
Engineer Bainomugisha
Department of Computer Science Makerere University, Uganda
Michael T. Smith Department of Computer Science University of Sheffield, UK
Richard D. Wilkinson School of Mathematical Sciences University of Nottingham, UK
Abstract
Linear systems occur throughout engineering and the sciences, most notably differential equations. In many cases the forcing function for the system is unknown, and interest lies in using noisy observations of the system to infer the forcing, well as other unknown parameters. In differential equations, the forcing function is unknown function of the independent variables (typically time and space), and can be modelled as Gaussian process (GP). In this paper we show how the adjoint of a linear system can be used to efficiently infer forcing functions modelled GPs, using truncated basis expansion of the GP kernel. We show how exact conjugate Bayesian inference for the truncated GP can be achieved, in many cases with substantially lower computation than would be required using MCMC methods. We demonstrate the approach equations, and show that the basis expansion approach approximates well the true foreing with a modest number of basis vectors. Finally, we show how infer point estimates for the non- non-linear model parameters, such as the kernel length-scales, using Bayesian optimisation.
systems of both ordinary and partial differential
1
Introduction
Linear systems are used as models throughout the sciences and engineering, encompassing wide range of both ordinary and partial differential equations (including the heat, wave, Schrödinger's, Maxwell's equations etc), as well as systems of linear algebraic equations (such as eigenvalue problems). To fix notation, let L 1 V be linear operator between Banach spaces U/ and V (i.c., complete normed spaces, Kreyszig 1991h) prototypical linear system is then the form
Lu=
(1)
where u is the quantity of interest being modelled, and f E is the forcing function of the system. Given fully specified operator L and forcing function ff (and possibly initial and boundary conditions), solving the system for u is referred to as the forward problem. Typically, this is
"Equal Contribution
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

computationally intensive task. For example, consider modelling air pollution as it moves through the atmosphere. In this case, C may be a partial differential operator describing the advection, diffusion, and reaction of the pollution, and ʃ will be function describing the source of the pollution at each location and time. The forward problem refers to computing the concentration, u, given the emission sources and will usually require the use of numerical integration methods.
In many applications, both the linear operator L and the foreing f may not be fully specified, and we may face the statistical task of learning L and f from noisy observations of u:
= hb(u) +
Here, E IR are the observations, h the observation operator, ni the number of observations, and E IR a zero-mean observation error. This is often referred to as the inverse problem in applied maths and statistics (Stuart, 2010), or sometimes as latent force model in machine learning (Alvarez et al. 2009). In the air pollution example this would equate to finding the distribution of pollution sources, given a set of concentration measurements 2.
We focus on the situation where
1. f is modelled as a Gaussian process (GP) when V is an infinite dimensional Banach space, with Gaussian distribution in the finite dimensional case. E.g., if L is ordinary differential operator with independent variable I, then ʃ will be an unknown function of t.
2. The observation operator h U + I" is affine. In the finite dimensional case, this implies hu()) Hu for some constant matrix H and vector c. but in the infinite dimensional case (when LI is a space of functions) includes pointwise evaluation of u, i.e. u(x, t) for some values of æ and t, and integral and derivative observations of u, c.g. f u drdf, de etc. 3. The observation error 6 has a Gaussian distribution. This assumption can be relaxed for maximum likelihood (ML) estimation.
de
The full specification of the statistical model is then
Lu f. ((t)))
== h(u) +e
(3)
N,,,00
where m and & are the prior mean and covariance (kernel) functions of the GP. Note that the linear system in Eqs. (3) may include initial and boundary conditions for differential operators. Our aim is to infer f (and possibly u, L, and k) given & either via
* ML estimation, by solving the constrained optimization problem
min /
(2 h()) (z h(u)) + all/lb
s.t. Lu
Bayesian inference, by computing the posterior distribution
(f, u|zz x (()lu)((l/f++(()
(4)
The prior distribution for (f),, and the regularization parameter, o, play similar role in the two approaches. To solve either problem numerically is likely to require many solves of the forward problem (Eq. [T. For example, for ML we may seek solution using numerical optimization (e.g. Arellano et al. 2007 Tröltzsch 2010). whereas with Bayes, we might use an MCMC scheme (c.g. Kopacz et al.) 2009 Cotter et al. 2013 Sengupta et al., 2016 Albani et al. (021) or variational
approach (Chappell et al. 2008). ). All of these approaches require multiple solves of the forward
problem. Here we develop an approach to minimize the computational cost of inference.
1.1 Contribution
In this paper we show that implementing an adjoint of the linear system can result in much faster statistical inference. Instead of using numerical approaches to solve either the ML or Bayesian inference problem, we can do inference for f at the cost forward model solves, where the number of data points. In many (but not all) cases, this will incur substantially lower computational cost than competing methods, such as MCMC. More specifically, we show that
1. if f depends linearly on parameters q, without resorting to numerical integration methods;
estimate q or its distribution analytically, i.e.
2. if we model / as a Gaussian process, then by using truncated basis expansion we can efficiently infer the posterior distribution for f.

---

The paper is structured as follows. In the next section we discuss related work before introducing adjoints in Section 3.1 We derive the main results in Section (..2 and in Section (.33 we show how linearizing Gaussian processes via basis expansion reduces inference for GPs in linear systems to simple linear algebra. Finally, in Section 4 we demonstrate the approach on two linear systems: ordinary (ODE) and partial differential equations (PDE).
2 Related work
The problem defined by Eqs. 2013). Alvarez al. (2009) showed how the posterior distribution 4 can be computed by using the G( f) where Gi-) is the Green's 's function associated with the differential operator L. Due to the linearity of the integral transform, placing a GP prior over f leads to joint GP over f and u. From this joint GP, the posterior distributions T (u (f) and T Tuu) can be computed in closed form However, many situations, particularly for non-trivial differential equation models, the expressions for the covariances are cumbersome and lead to the use of error funetions with complex arguments or functions like the Faddeeva function that can be numerically unstable to compute. Guarnizo & Alvarez (2018) proposed representing using random Fourier features (RFFs) to reduce the number of integrations necessary to be solved analytically. Here, rather than using Green's 's functions, we instead use adjoints to write the problem as a linear model and then use a reduced-rande Gaussian process formulation, leading to numerically stable and fast approximations to the posterior distribution.
B is often referred to as latent force model (Alvarez et al. 2009
integral formulation of Lu f. i.. u(x)=
More specifically, estimation of forcing functions differential equation models has been extensively studied, for example, in the field of modelling atmospheric advection-diffusion (e.g. Yee| 2008
Boryssewicz et al., 2012, Singh Rani [0014] Rajaanaa 2016 Yeo et al.] 2019) with some authors solving the inverse problem using an adjoint approach in combination with MCMC to compute the posterior distribution (c.g. Yee, 2008; Hwang et al., 2019; Luhar et al. 2020 Albani et al. (02I). Of particular relevance is Hwang et al. (2019) who use point source model of pollution, and use the adjoint to write the backward finite difference approximation, noting that this can be written as a linear model, where the features are conjugate fields associated with each sensor. MCMC sampling is still a limiting factor, restricting the extension of the approach to more complex situations such as time-varying pollution sources.
Other related work includes the stochastic PDE approach of Lindgren et al. (001)) in which the underlying function is modelled as a Gauss Markov random field (GMRF), which can then be formulated as stochastic partial differential equation. This allows finite element methods to be used to efficiently compute the posterior distribution for u given z. Similarly, Hartikainen & Särkkä (2010) exploit the link between GMRFs and dynamical systems to convert inference for u to form in which Kalman filtering methods can be used, which scale linearly with Sigrist et al. (2015) focus exclusively on the advection diffusion PDE considered in Section 4.2 and use white noise for the forcing funetion, f. to create stochastic PDE model for They solving the PDE in the Fourier domain, to develop an efficient algorithm for statistical inference. Whilst attractive, it is difficult to generalize this approach from white noise models of f to correlated Gaussian process models. Jidling et al. (2017) ) show how to infer statistical inference for linearly constrained systems, i.e., where Lu 0. These approaches all model u, whereas our focus is on inferring the forcing function f.
spectral approach,
3 Methods
We first recap the definition of adjoints, before deriving our main result illustrating how they can be used to accelerate inference. We then show how using a truncated basis approximation to Gaussian processes allows us to find the GP posterior distribution without resorting to MCMC methods.
3.1 Adjoints
Recall that L 24 - V is linear mapping between Banach spaces 24 and V. Let lt* and V* denote the dual spaces of 2/ and V (Kreyszig) 1991). We can construct the adjoint to a continuous linear
For more detailed explanation of the Green's function method, please the supplementary material.
3

---

operator L as follows. Let v* E V* and define F - R by F + v Then F is a bounded linear functional on U, i... u² for some u² E ut. Thus for all o E V we' 've associated unique u² E U
:
= C.
(5)
We call L the adjoint of C, and L is itself a bounded linear operator Estep 2004). By construction, we have that for all u U and v* EV*
(L"v")u)) ( (Cu).
(6)
a result known as the bilinear identity. In the case where U1 and V are also real Hilbert spaces with respect to inner products (:-) and (--)v,, then we can identify the dual spaces with their underlying space: by the Riesz representation theorem if o" E V* then there exists v E V such that v" (+) (: v.v. In this case, the bilinear identity reduces its more familiar form
(Cu, 4/v v² (Lu) (C* ("))))) (u, u)u
where we now consider L V - U. We only consider real vector spaces here, resulting in symmetric inner product.
Generally, the adjoint L- will be the same type of operator as L (e.g. if L is a differential operator then L will be too), so solving an adjoint system of the form L will have similar computational complexity as solving Lu = See Hestep (0044) for an introduction to adjoints, and Section 4 for examples of adjoint systems.
3.2 Benefits of adjoints
How does the development of an adjoint to a linear system help us perform statistical inference for that system? Consider the situation where uncertainty about the unknown forcing function, I in Eq. (1, can be characterized by linear dependence upon unknown parameters That is, we can write
M
f())
Σ fm0m(()) (.). 1==1
(8)
In the infinite dimensional case where 24 and V are spaces of functions some set X' the om will also be functions on X. In the finite-dimensional case, the om will vectors of fixed length.
In the situation where the observation operator 22 linear, then we can write the ith observation as (u) (h,, u) plus noise, for some h2 E U. Consider the different adjoint systems
L h; for i 1...., n.
Then using the bilinear identity 7 we have that
h.(w) = (h,, u) (L",,, u) = (ws, Lu) (06 f),
ie., the jth observation is the inner product between the unknown forcing function and the solution of the jth adjoint system. At first, the introduction of the adjoint doesn't appear to have helped. To evaluate the likelihood (or sum of squares) we have gone from needing a single solve of the forward problem, to requiring the solution to n adjoint systems. The benefit arises if we now use the assumption of a linear dependence upon the parameters, and linearity of inner products:
M
h₂(u)= (v.,
Im 2m))
Σ
Om (T. pm)
The complete observation vector z can then be written
(U1, Φ1)
(V1, (M)
91
2
+ bbg
(9)
(wa Φ1)
(Uns 2x))
gMs
where y E M the parameter vector, and + E RRXM the matrix of inner products between the adjoint solutions and M basis vectors.

---

This can be recognized as a linear model. Thus, standard results can be used to compute the least squares/ML and Bayesian estimators. For ML, the minimum of S( (2 h(u)) (2 h(u)) subject Lu = f. is obtained at ( T ) - with Var(q) = o² + Φ) in the case where the observation errors are uncorrelated and homoscedastic with variance o². Standard results can be used from regularized least squares if we need to regularize g.
In a Bayesian setting, we assume a priori that q ~ Nu (Ho, Eg), then the posterior for q given 2
(and other parameters) is
q 2 ~ (Ns En)
(10)
where
1 o2
En - (1200 5
Pn En,
$ 2
Do Ho)
(11)
See, e.g.. e.g. O'Hagan & Forster (2004). Note that the solutions of the adjoint equation, vi, are present in the posterior distribution of the matrix 4.
3.3 Inference of Gaussian forcing functions
We now consider the case where the forcing function given a Gaussian process prior distribution:
((-) ~ GP(m())
(12)
where m(.) and -): the prior mean and covariance functions respectively (Rasmussen & Williams 2006). Our approach is to use a reduced-rank representation of the GP, as in Eq. derived by truncating an expansion for k of the form
k(I,2') Σ pm (z)))m2(2)) (2').
(13)
There are many possible choices for the basis vectors (-). including the Karhunen-Loee (KL) (Deheuvels& & Martynov 2008), Laplacian (Solin & Särkkäl 2020 Coveney et al., 2020 Borovitskiy et al. 2020), and random Fourier feature (RFF) (Rahimi et al. 2007) expansions. The Karhunen
Loeve basis is formed by finding the spectral expansion of the integral operator defined in Mercer's theorem:
Txf(x)) - J/ (a,z'))(z')da
In the case where dim(x) is small (such as the ODE example below), the eigenfunctions of Ti are easy to compute either analytically (Rasmussen & Williams 2006, Section 4.3.2, p1 p6 or numerically (Greengard & O' O'Neil 2021). Using the eigenfunctions of Tk gives the L L'-opii approximation (Kosambil 2016), but the eigenfunctions can be difficult to compute even in three dimensions. A simpler approach that extends easily to higher dimensional problems, is to use RFFs (Rahimi et al. 2007). This relies upon Bochner's 's theorem, which can be used to express stationary kernels k( x k(x - x') as the Fourier transform of a positive measure p
k(x - x') = exp(-i - d')))
If an isotropic exponentiated quadratic (EQ) kernel k, then the measure p corresponds to multivariate Gaussian, Rahimi et al.] (2007) give expressions for the Fourier feature basis vectors that can be used to approximate the Gaussian process:
1
(m) = - V2r
2T² M
1 cos( w + bm) (14)
k(# - x') = 72 expepppppppp
2x (x - (x-xx===)
where Wm ~ N (0, and bm ~ (0. Substituting this om into equation alllow us to write the forcing function ff as truncated Gaussian process with an EQ kernel (as such the true Gaussian process is never calculated). The extension to anisotropic kernels is straight-forward.
Although the RFF expansion will require more basis vectors (a larger M) than the KL expansion to achieve the same accuracy, the computational complexity of our adjoint approach is dominated by the number of adjoint solves, not the number of features (which only affects the cost of computing the relatively low-cost Eq. [3]. and so including more terms has a minor effect on overall computational cost.
5

---

3.4 Time complexity of Algorithm 1
Let G be the number of grid elements, n the number of observations, and M the number of features. There are five stages in our algorithm:
1. Solving the adjoint systems, which requires (Gn) operations.
2. Computing each basis vector, Φ, over the grid for each feature requires (GGM) operations.
3. Computing the matrix 2, requires (G)M)) operations.
4. Finding & & requires O(nM²) operations.
5. Finally solving the matrix inverse will require (MM)) operations.
This results an algorithm that scales linearly in the number of data points, n. Empirically, for the problems we' ve experimented with, we found computing over the large grid was the most time consuming step. The overall time complexity is OGGMn) + Onn2²) + O(M)), but note that the constants associated with each term are important in most problems.
4 Experiments
4.1
Ordinary differential equation (ODE)
Consider the non-homogeneous linear ODE:
duu di
du
P2
P1 dt + pou = f(t)
(15)
on the domain [0, T] with initial conditions 2()) n1(0) 0. The right hand side, f(), is the unknown forcing function that we wish to estimate, and pn. P1, p2 are parameters in the linear operator. We model f as a zero-mean GP with an EQ kernel (Eq. 14). We assume observations are obtained as noisy averages over short time windows:
A+ 1
hu(u)=
u(t)dt = (u,
(16)
At
where hū is the indicator function h₂(() [(td (4) (t). This includes direct measurements of u(t;) we set (t) = 8(t - t;), the Dirac delta function. To generate synthetic data, we simulate a realization from the GP model, solve Eqq. (15) for u(I), and then simulate observations from Eq. (16) with T P2 0.5, P1 1, Po
f; iT adding zero-mean Gaussian noise with standard deviation o = 0.1. For simplicity, we solve the ODE with a simple forward Euler approximation, but higher order schemes and should be used in real applications. We approximate the GP using Eq. (8), using 200 RFFs generated using Eq. 14 with À v0.6 and 7² 4. The linear operator in this case is
At
d dt
d po
Lu
P2
+
+
u.
dt
To derive the adjoint operator we use the bilinear identity (Eq. 7]. and integrate by parts twice:
(Cu, v) = S' (Cu)udt /'' (pzü + pru u + pou)uct = 1/ (pzū - prú + vjudt (u, L*v)
when b(7) = ((7) = 0. So the adjoint of L is
d² P2
C'=
dtz
Pi dt
+ Po
Note that rather than an initial condition, the adjoint has a final condition: to solve the system we have integrate backwards in time from T to
See Algorithm
Fig. shows the effects of the number of training points n and the number of features M on the posterior distribution of the forcing function. As expected, more data results in more confident

---

Algorithm Computing the posterior distribution of 4
for = 11...M n do
Solve adjoint system L² h4 with appropriate final and boundary conditions. end for
for m 1. M do
Sample an RFF basis vector 6m using Eq. 14
for i ..... n do
Compute {+) (Vi, pm)-
end for end for Compute the posterior distribution for q using Eqs. (00) and
1/
18
15
16
85
16
1/
15
166
18
46
8
160
125
166
18
Figure 1: Posterior distribution for the unknown forcing function (with 95% credibility interval). True forcing in red. The number of training points and features are, (left) = 10 and M = 100; (middle) n 100 and 10: (right) n 100 and M 100. The overconfident and wrong posterior when M = 10 is a consequence of the model being heavily misspecified.
Note though the danger of using too few features: with M 10 the approximation has limited expressive power and cannot capture the true form of f. i... the model is heavily misspecified. This can result in the uncertainty collapsing upon the most likely, but wrong, value. This can be diffieult spot, users should check the sensitivity of the posterior with respect to M (which can be done with no additional forward solves).
posterior.
In the supplementary material provide comparison between the adjoint method, the Green's function methods as in Alvarez et al. (2099) and Guarnizo Alvarez (2018), and a vanilla Gaussian process when applied to this ODE problem.
4.2 Partial differential equation (PDE)
We now demonstrate the approach on PDE in which there are two spatial variables, x E C R and time, E [0, 7 c IR, so that the solution u = u(x, () is function of three independent variables. We consider the advection-diffusion equation
Ou Øt
+ Vu V (22 in X [0,
(17)
with initial and boundary conditions u(x; 0) for x X' and V, the unknown forcing function f f(x t) is function of space and time, and we model it as zero-mean Gaussian process, f(z t) ~ GP (x, (x',t') (')) with EQ kernel & (Eq. We use an RFF approximation to A with random weights Wm No I) and bm ~ (0, 2T)).
for E OX. Here,
Observations are assumed arise from sensors which take an average of u(x, 1) over small spatial and temporal window R4 x T C X x [0,T]:
if E R, and e T
44 = (u, h4) +
with hg = [...I.
otherwise.
The adjoint of the linear operator Lu
Ou dt
+ P1. Vu V (p2 Vu) is
dr
L'
P1 Vu V Vv).
7

---

10
10
0
2
8 10
0
2
8
10
0
2
8 10
Ground truth
Four sensors
16 Sensors
Figure 2: Spatial maps of the forcing function single slice: f(r 5). Left shows ground truth,
middle shows the posterior mean with 4 sensors, right shows the posterior mean using 16 sensors.
Table 1: Posterior mean and standard deviation of the q parameters estimated by MCMC and the adjoint method.
91
q3
95
MCMC
-3.62 (0.02)
-0.64 (0.01)
1.68 (0.01)
-0.09 (0.01)
4.31 (0.02)
Adjoint -3.62 (0.02)
-0.64 (0.01)
1.68 (0.01)
-0.09 (0.01)
4.31 (0.02)
Our adjoint-aided approach then requires the solution of
L 0/ h; in X' [0, 7] for i n. The final and boundary conditions on the adjoint system are 0s (x, T) = 0 for x E X, and (r,;(z,t) + p2 Ve( = 0 for x E an and f [0,
For details of the adjoint derivation see the supplementary material and Estep (004)).
Inference then proceeds as before. After solving the n adjoint systems, we compute the inner product of these solutions with the RFF basis vectors form the matrix & (Eq. 1). We can then use Eq. [I] to compute the posterior with minimal additional computational cost.
Data was simulated the spatial domain X' [0, 10 for t [0, 10] by first randomly generating a forcing function fz, t) (generated from a GP using an EQ kernel with À = 2, 7² = 2), and then solving the forward problem (Eq. 17 find u(x,t) t) using PDE parameters P1 (0.4, 0.4) and P2 = 0.01. We generate n observations using sensors that record averages over short time windows equally spaced across the domain [0, 10] at the locations shown in Fig. Zero-mean Gaussian distributed noise is added to the true sensor readings with standard deviation o = 0.05 (note that this relatively small compared the signal, which often create problems for sampling methods). We then use Algorithm| to calculate the posterior distribution for y, hence giving the posterior for f. By sampling forcing functions from this posterior and simulating forward, posterior predictive accuracy of the model.
evaluate the
To validate the posterior estimates from the adjoint method, we also used MCMC to compute the posterior distribution for the PDE model using just M 10 RFFs. Table T shows the posterior mean and variance of the first parameters determined using both methods. which can be seen to be in close agreement. Fig. in the supplementary material shows the trace plots. We used a batch-update random-walk Metropolis- Hastiing (MH) sampler, which failed to converge (after hour of computation) when using a larger number of features. Even with M = 10, we still required 10,000s of forward model evaluations to reach convergence (whereas in comparison, the adjoint method required just 75 forward model solves in this case).
To illustrate the effects of sensor density, observations were generated for five time windows using arrays of and 16 sensors at the locations shown in Fig. 3 i.c., a total either n = 20 or n 80 observations. The ground truth forcing and inference results (with M 200) are shown in Fig. As expected, more sensors results improved estimates. Fig. 3 shows the posterior standard deviation for ff in the 4 and 16 sensor cases. Here advection occurs parallel to the (i.e., as there is wind blowing from the south west). We can see that standard deviation is smallest upwind of the sensors, with more uncertainty downwind, expected.
line p1 (0.4. 0.4)
To investigate the effects of varying feature and sensor numbers we performed a posterior predictive check using held-out data and used Monte Carlo estimation calculate the posterior predictive mean

---

10
10
x
8
x
1.0
6
X
X
X
x
x
x
0.5
x
x
x
X
0
10
0
2
4
10
Sensors
16
Sensors
Posterior standard deviation for the examples shown in Fig. The left image is for the case where 4 sensors are used, and the right for 16 sensors. Sensor locations are shown as white crosses.
Figure 3:
Table 2: The median MSE function of number of sensors and RFFs. The ground truth was generated using a GP with À = 2, 7² 2 and PDE parameters P1 = (0.01, 0.01), 0.01. The MH algorithm did not converge after 20,000 iterations for 50 more RFFs. The numbers brackets are the 95% confidence intervals computed from 10 repeated experiments (for the adjoint approach) and (for MH).
Sensors
Features
10
200
300
(2.82,4.19)
3.27 (3. (3.138888
3.24 (3.10,3.37)
3.27 (3.17,3.44)
3.24 (3.12.336)))
3.27 (3.17.3.35)
7.12 (1.7.888))))
(2.06.262)
(2.13,2.6)
2.45 (2.32.2.57)
(2.41.2.9)
(2.32,2.60)
2.38 (1,44,4440
(1.4...))))
1.70 (1.49.207))
1.48 (1.40.1.72)
1.47 (1.32,1.61)
1.45 (1.40.1.50)
1.73 (1.23...33
3.99 (2.32.10.90)
2.18 (1.72,3.54)
(1.02.))))))
1.12 (0.98.1.37)
1.12 (1.02, 1.21))
25
1.35 (1.1,,0)))))
8.93 (4.92.99888
4.36 (2.53.8.20)
1.86 (1433.2775
(1.07...)
1.05 (0.89.1.45)
25 (MH)
3.27
(1.73.6.12)
squared error (MSE). Table 2 gives the MSE as the number of sensors and RFFs vary. As both increase, so does the accuracy of our estimates. In general, the accuracy depends upon a variety of factors, including the PDE parameters (ratio of diffusion advection), kernel parameters (decreasing lengthscale makes the problem more challenging). and sensor locations. The speed and efficiency of the proposed adjoint-aided approach allows us investigate these effects in way that would not be possible if we were using MCMC (as each estimate of the posterior requires tens of thousands of simulator evaluations, rather than just the evaluations required for the adjoint- -ided approach).
Although not our focus here, we note that we can infer the remaining parameters p (PDE) and (r² (GP hyperparameters) in variety of ways. For example, using Bayesian optimization the maximum likelihood estimates of these parameters with relatively minimal computational cost. See the supplementary material for details and examples.
find
In the supplementary material we include some examples of the effect of varying the Gaussian process lengthscale on inference quality, example in which apply method to advection diffusion problem using the Roundhill II dataset (Cramer & Record 1957) and apply our method to shift operators non-differential linear operator).
4.3 Real time cost analysis
We compared the run time of the adjoint method to basic MH MCMC algorithm (recorded on a laptop with 16GB RAM and an Intel 17-1065G7 CPU @ 1.50 GHz). Exact inference of the posterior of ff in the ODE model using the adjoint method observations using 100 RFFs took 569 + 72.5 ms. This is approximately the time took to perform single iteration of the MH algorithm. Inference of the posterior of f in the PDE model using the adjoint method on a 50 x 30 x 30 spatiotemporal grid with 100 observations using 100 RFFs took 3.32 + 0.18 In this time were able to 1 iterations of the MH algorithm.
100 element temporal grid with 100
The Metropolis Hastings algorithm comparatively inefficient due to the matrix multiplication required to compute the forcing function from a given using the basis vectors (see equation 8). In the adjoint method this only needs to be calculated once, but is required at each step of the MH algorithm.

---

Whilst more sophisticated MCMC implementations will achieve faster convergence than random walk samplers, any MCMC algorithm will require 10,000s of forward model evaluations to approximate the same posterior the adjoint method can compute exactly B in fraction of the time.
5
Discussion
This work was motivated by the problem of inferring the distribution of spatially and temporally varying sources of pollution across city from noisy observations using model of the pollution's atmospheric transport. Estimating the pollution concentration and its sources can help local authorities reduce the population's 's exposure, and motivate policy interventions. Linear systems such as this, are typical of the type of challenge faced throughout the sciences and engineering, and still pose computational challenges that make principled statistical inference almost intractable, if they tractable, the computational cost (e.g. of using MCMC) is such that only a limited range of models and situations can be analysed. In this paper developed an approach that results in conjugate Bayesian inference of an unknown GP forcing function, with a computational cost that scales linearly with the number of observations, n. As n increases, the approach may eventually require more computation than competitor methods such as MCMC (which in theory has cost independent of n)), but given that MCMC typically requires 10 10th iterations even for low-dimensional problems" there is still a large range of problems for which an adjoint may be beneficial. In our PDE example, the adjoint-aide
approach required orders of magnitude fewer PDE solves than MCMC. We should note that a key limitation of our approach is that it only applies to linear systems and cannot be used to determine forcing functions for non-linear systems such as the Navier-Stokes equations.
We only briefly touched upon the problem of estimating the operator's non-linear parameters (e.g., advection and diffusion rates) and the GP hyperparameters. Adjoint sensitivity methods (see, e.g. BBradley, A. .M 20100: Margossian 2019 Yashchuk (020) can be used to estimate gradients of the log-likelihood with respect to the non-linear parameters in a numerically stable alternative to weeommati-difffeenniie
frameworks (such as TensorFlow), which can be unstable when back propagating through long iterative loops. Access gradient information allows inference of the additional parameters to be performed efficiently within the preferred statistical paradigm (e.g. maximum likelihood using Bayesian optimization see the supplementary information; Hamiltonian Monte Carlo (Neal) 2011); or a variational autoencoder (Kingma & Welling, 2014) framework etc.).
Finally, there are many ways in which this approach may be accelerated, for example, by the use of intelligent numerical solvers that reuse solution trajectories and adaptive step size solvers, multi- fidelity methods that varying grid sizes, and stochastic approaches which use only subset of the data at each stage. The core computational tasks in our approach (solving the adjoint systems) are embarrassingly parallelisable, enabling easy deployment on HPC facilities if required.
Software:
The algorithm has been implemented (for both the ODE and PDE problems) in a Python module available notebooks that are used to produce the figures and tables in this paper.
The repository also contains Jupyter
Acknowledgements
This work was directly funded by EPSRC projects EP/T00343X/1 and EP/T00343X/2. In addi- tion, EB was supported by funding from Google.org and RW by EPSRC projects EP//010741/1, EP/W00009 1/1, and EP/X012603/1.
'Here "exactly" refers to our ability to write down and compute the posterior for given basis truncation of the GP, rather than approximating the posterior with MCMC.
'In the best case scenario of independent Gaussian posteriors, the number of iterations required for random
walk Metropolis Hasting scales as ( (Gelman et al. 1997) and (M i) for Hamiltonian Monte Carlo, [Beskos et al. (013)) where M is the dimension of the parameter, with each iteration requiring solve of the forward problem. For more complex problems, scaling rates will be worse.
10

---

References
Albani, R. A., Albani, V. V., Migon, H. S., and Neto, A. J. S. Uncertainty quantification and atmospheric source estimation with discrepancy- and state-dependent adaptative MCMC. Environmental Pollution, 290:1 18039, 2021.
Alvarez, M., Luengo, D., and Lawrence, N. D. Latent force models. In Artificial Intelligence and Statistics, pp. 9-16. PMLR, 2009.
Alvarez, M. A., Luengo, D., and Lawrence, N. D. Linear latent force models using Gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(11):2693-2777 nov 2013, ISSN 1939-3539. doi: 10. 10.1109/PAM..2013.88
Arellano Jr, A., Raeder, K., Anderson, J.. Hess, P., Emmons, L., Edwards, D., Pfister, G., Campos, T. and Sachse, G. Evaluating model performance of an ensemble-based chemical data assimilation system during INTEX-B field mission. Atmospheric Chemistry and Physics, 7(21):5695-5711 2007.
Beskos, A., Pillai, N., Roberts, G., Sanz- J.-M., and Stuart, A. Optimal tuning of the hybrid monte carlo algorithm. Bernoulli, 1965A):501-1534, 2013.
Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. Matérn Gaussian processes on Riemannian manifolds. arXiv preprint arXiv:2006. 10160, 2020.
Borysiewicz, M., Wawrzynczak, A., and Kopka, P. Bayesian-based methods for the estimation of the unknown model's parameters in the case of the localization of the atmospheric contamination source. Foundations of Computing and Decision Sciences, 37(4):253, 2012.
Bradley, A. M. PDE-constrainee optimization and the adjoint method. Online tutorial, https //c8. stwnsordd.d...o.. attaaford.edd//"abraaaaaaad////d aatprr///adoointtttttiiaall pd. 2010. Accessed: 2022-01-27.
Chappell, M. A., Groves, A. R., Whitcher, B., and Woolrich, M. W. Variational Bayesian inference for a nonlinear forward model. IEEE Transactions on Signal Processing, 57(11)223-22) 236, 2008
Cotter, S. I.., Roberts, G. O., Stuart, A. M., and White, D. Meme methods for functions: modifying old algorithms make them faster. Statistical Science, 28(3):424 446, 2013.
Coveney, S., Corrado, C., Roney, C. H., O'Hare, D., Williams, S. E., O'Neill, M. D., Niederer, S. A., Clayton, R. H., Oakley, J. E., and Wilkinson, R. D. Gaussian process manifold interpolation for probabilistic atrial activation maps and uncertain conduction velocity. Philosophical Transactions of the Royal Society A, 378(2173):20190345 2020.
Cramer, H. E. and Record, F. A. Field studies of atmospheric diffusion and the structure of turbulence. American Industrial Hygiene Association Quarterly, 18(2):126-131) 126-131, 1957.
Deheuvels, P. and Martynov, G. V. A Karhunen- decomposition of Gaussian process generated by independent pairs of exponential random variables. Journal of Functional Analysis. 255(9): 2363 2363-2394, 2008
Estep, D. A short course duality, adjoint operators, Green's functions, and posteriori error analysis. Lecture Notes, 2004.
Gelman, A., Gilks, W. R., and Roberts, G. O. Weak convergence and optimal scaling of random walk metropolis algorithms. The annals of applied probability. 7(1:1110-120 1997.
Greengard, P. and O'Neil, M. Efficient reduced-ram methods for Gaussian processes with eigenfunc- tion expansions, arXiv preprint Stt.:::200000922222 2021.
Guarnizo, C. and Alvarez, M. A. Fast kernel approximations for latent force models and convolved multiple-output Gaussian processes. In Uncertainty in Artificial Intelligence, pp. 835-844, 2018
Hartikainen, J. and Särkkä, S. Kalman filtering and smoothing solutions to temporal Gaussian process regression models. 20/0 IEEE international workshop on machine learning for signal processing, pp. 379-384, IEEE, 2010.

---

Hwang, Y., Kim, H. J., Chang, W., Yeo. K., and Kim, Y. Bayesian pollution source identification via an inverse physies model. Computational Statistics & Data Analysis, 134:76-92 2019.
Jidling, C., Wahlström, N., Wills, A., and Schön, T. B. Linearly constrained Gaussian processes. Advances in Neural Information Processing Systems (NIPS), 2017.
Kingma, D. and Welling, Stochastic gradient VB and the variational auto-encoder. In Second International Conference on Learning Representations, ICLR, volume 19, pp. 121, 2014.
Kopacz, M., Jacob, D. J., Henze, D. K., Heald, C. L., Streets, D. G., and Zhang, Q. Comparison of adjoint and analytical Bayesian inversion methods for constraining Asian sources of carbon monoxide using satellite (MOPITT) measurements of CO columns. Journal of Geophysical Atmospheres, 114(D4), 2009.
Research:
Kosambi, D. Statistics in function space. In DD Kosambi, pp. 115-123. Springer, 2016.
Kreyszig, E. Introductory functional analysis with applications, volume 17. John Wiley & Sons, 1991.
Lindgren, F. Rue, H., and Lindström, J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodoloyy), 73(4):423-44) 498, 2011.
Luhar, A. K., Etheridge, D. M., Loh, Z. M., Noonan, J.. Spencer, D., Smith, L.,, and Ong, C. Quantifying methane emissions from Queensland's 's coal seam gas producing Surat Basin using inventory data and regional Bayesian inversion. Atmospheric Chemistry and Physics, 20(23): 15487-1551 1, 2020.
Margossian, C. C. A review of automatic differentiation and its efficient implementation. Wiley interdisciplinary reviews: data mining and knowledge discovery, 9(4):e1305 2019.
Neal, R. M. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 2(11): 2, 2011.
O' O'Hagan, A. and Forster, J. Kendall's advanced theory of statistics, volume 2B: Bayesian inference, volume 2. Amoll, 2004.
Rahimi, A., Recht, B., al. Random features for large-scale kernel machines. In NIPS, volume 3, pp. 5. Citeseer, 2007.
Rajaona, H. Inférence Bayésienne adaptative pour la reconstruction de source en dispersion atmo- sphérique. PhD thesis, Lille 1, 2016.
Rasmussen, C. E. and Williams, C. K. I. Gaussian processes for machine learning. MIT press Cambridge, MA, 2006.
Sengupta, B., Friston, K. and Penny, W. D. Gradient-based MCMC samplers for dynamic causal 2016.
modelling, Neurolmage,
Sigrist, F., Künsch, H. R., and Stahel, W. A. Stochastic partial differential equation based modelling of large space-time data sets. Journal of the Royal Statistical Society: Series Statistical Methodology. pp. 3-33, 2015.
Singh, S. K. and Rani, R. A least-squares inversion technique for identification of point release: Application to fusion field trials 2007. Atmospheric environment, 92:104-117, 2014.
Solin, A. and Särkkä, S. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing. 30(2):419446, 2020.
Stuart, A. M. Inverse problems: A bayesian perspective. Acta Numerica, 19:451-559, 2010. doi: 10.101771096242210000060
Tröltzsch, Optimal control of partial differential equations: theory, methods, and applications, volume 112. American Mathematical Soc., 2010.
12

---

Yashchuk. I. Bringing PDEs to JAX with forward and reverse modes automatic differentiation. In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020.
Yee, E. Theory for reconstruction an unknown number of contaminant sources using probabilistic inference. Boundary-layer meteorology, 127(3)355--94,, 2008
Yeo, K., Hwang. Y., Liu, X., and Kalagnanam, J. Development of hp-inverse model by using generalized polynomial chaos. Computer Methods in Applied Mechanics and Engineering, 347: 1-20, 2019.
13

---

Checklist
For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [es]] The two main claims in the abstract and introduction are (a) The application of the adjoint approach, which is described in (b) The benefits over, for example MCMC: Comparisons are made with MCMC. For example, Table
and the relative times the two approaches took are discussed in Section
(b) Did you describe the limitations of your work? [Yes] In the introduction (Section we emphasise the approach is only for linear systems. In the discussion (Section we raised the main issues around computational complexity.
(c) Did you discuss any potential negative societal impacts of your work? [Yes] At the start of Section [ we briefly consider that, because this approach could improve source attribution it might be used by city authority to motivate policy interventions that discriminatory.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] The three point list in the introduction (Section [ states the assumptions made.
(b) Did you include complete proofs of all theoretical results? [Yes] We step through the derivation of the inference approach in Section 3
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URLY? [Yes] We will link, in the discussion, to the github repository from where the python module can be installed and where jupyter notebooks (that made the figures etc) can be downloaded. The link has been temporarily redacted to avoid breaking the double-blin blind review process.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] The datasets were synthetic, but the hyperparameters, kernels, etc that were used to both generate them and model them were recorded in the paper. For example Section 4.2) records the noise scale (a = 0.05), the kernel hyperparameters
(A = 72 2), the number of RFFs (M 200) and the advection and diffusion parameters.
(c) Did you report error bars (e.g., with respect to the random seed after running ex- periments multiple times)? [Yes] Although the actual inference is deterministic, a stochastic aspect is caused by the sampling of the RFFs. Table report the 95% confidence intervals over the multiple runs.
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, or cloud provider)? [Yes] Section includes detailed runtime analysis of the adjoint approach and MCMC iterations.
4. If you are using existing assets (c.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] The algorithm as implemented does not require dittoddifferentation. and the MCMC implementation does not use a library. The only tool we use is GPyOpt, which we cite.
(b) Did you mention the license of the assets? [Yes] We included this information with the reference the bibliography.
(c) Did you include any new assets either in the supplemental material or as URL? [Yes] We link to our code in the discussion.
(d) Did you discuss whether and how consent was obtained from people whose data you're 're using/curating? [N/A] No data was used from people.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] Only synethetic data was used in this paper.
5. If you used crowdsourcing conducted research with human subjects...
14

---

(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] We didn't use crowdsourced data human subjects.
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] We didn't use crowdsourced data or human subjects.
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A] We didn't 't use crowdsourced data human subjects.
15


Reviews:
[{"review_id": "usTgsRpmbb", "reviewer": "Reviewer_mA14", "summary": "The authors developed a method to learn unknown forcing functions $f$ of linear systems $Lu=f$ from (affine) observations of $u$, also for PDEs. For finite dimensional learning, this is done symbolically by means of the adjoint. For infinite dimensional learning, a Gaussian process prior is made \"finite dimensional\" by means of a truncated basis expansion.", "questions": "**Questions**\n- I might have failed to understand the explanation for the linear runtime when using a GP and would like the authors to elaborate on this. Is the GP never calculated, but rather $z$ is estimated and using the adjoint, the resulting $f$ can be estimated?\n- How general is the suggested method really? Did you try shift operators?\n- Is your method applicable to the case when u is not given via observations but as a GP itself? I.e., can you solve linear operator equations for GPs?\n\n**Suggestions**\n- A quick reminder in the beginning of 3.1 that $L:\\mathcal{U} \\to \\mathcal{V}$, might improve the flow of reading\n- Inconsistencies in how the authors write equations\n\t- line 172 Eq. 13) and line 179 (Eq. 14) \n- The typesetting of the small brackets in the left part of equation (11) is awkward.", "limitations": "The authors have not brought up limitations of their work. I do not see a necessity, apart from the above weaknesses.", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\nThe authors leverage non-trivial ideas from functional analysis to construct better models in a relevant and significant area of research.\nThe technique is not original in mathematics, but it is a non-trivial extension to machine learning.\nComparison to the state of the art is clearly given (even though, of course it could be longer, but the number of pages is limited).\nThe language is clear and readable.\n\n**Weaknesses**\n1) After reading the paper, there remains a big gap/error in the functional analysis part of the paper. (This might just be a gap/error in the theoretical justification of the method. I do not claim that the method is wrong.)\n   Under one viewpoint, three kinds of function spaces appear:\n     (i) Hilbert space. They are used a lot, e.g. to compute adjoints. And the realizations of a GP with EQ kernel lie in such spaces (non-trivially, see arXiv:1807.02582).\n\t (ii) Banach spaces. I do not know why these spaces appear in the paper.\n\t (iii) Frechet spaces, which is the class of spaces your GP-realizations with the EQ kernel have a measure on (see arXiv:2205.03185). (These realizations lie in the space of smooth functions, which is not Banach w.r.t. its usual topology.) \n   Of course, there are Banach spaces suitable for GPs (probably for the Matern kernel, though I do not have a reference).\n   Furthermore, the computation of the adjoint as elements in a Hilbert space (using the Riesz representation theorem) depends on the choice of the scalar product. In the example in Section 4.1 such a scalar product is chosen. Hence, the result depends on this choice. This choice (and hence the adjoint, and hence the entire computation) are arbitrary. Even as we have learned in arXiv:1807.02582, there seems to be no canonical RKHS, which contains the GP realizations, and hence no canonical choice of a scalar product. I would assume that the choice of the topology to define the adjoint is clearly justified.\n   \n   Perhaps I am missing a trivial point here, but I did not find this point in the paper. A summary of the paper according to my feeling is: \"Talk about Hilbert and Banach Spaces. Some math magic happens, which we do not explain. This justifies our method.\" Again, I am not claiming that the method is wrong. But please either explain this math magic or do not try to sound mathmatical.\n   If the authors can clarify this point convincingly in a response, I am willing to drastically improve my scores, as otherwise (apart from the following points (2) and (3), which are not too damning) the paper is rather good.\n  \n2) In my opinion, the authors did not clearly show the benefit of their method experimentally. Neither a comparison with a standard GP was made, nor did the authors compare themselves to related work that might have tackled the same task. Adding this to the experiments would greatly improve the papers value.\n3) Experiments on real data would have helped. When dealing with differential equations, synthetic data is usually to \"nice\".", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "axp4izgBIf", "reviewer": "Reviewer_bJxr", "summary": "The paper discusses a method that infers the forcing function $f$ of the differential equation $L u = f$ from linear observations of $u$, $z = H u + \\epsilon$, $\\epsilon \\sim N(0, \\sigma^2)$.\nIt assumes a Gaussian process prior over $f$, $f \\sim GP(0, k)$, which is approximated with basis functions, $f(x) = \\sum_k q_k \\phi_k(x)$, for instance random Fourier features.\nThe essence of the paper is that in the above setup, the observation model can be written as a linear model $z = \\Phi q + \\epsilon$, where $\\Phi = [\\langle v_j, \\phi_k \\rangle]_{i,j}$ is the inner product matrix of the basis elements $(\\phi_k)$ and the solutions $(v_j)_j$ of the adjoint equation, $L^* v_j = h_j$, where $h_j$ is the observation functional of the $j$th data point.\nThis linear model $z = \\Phi q + \\epsilon$ can be used to infer $q$ from $z$ with standard Gaussian inference tools.\n", "questions": "See \"weaknesses\" above. ", "limitations": "-", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n* The topic is very relevant \n* The methodology is sound and could be useful\n* The experiments suggest that the method works\n* The related work is more or less covered adequately (with some open questions, see below)\n\n\n## Weaknesses\nWhile the methodology itself is sound, \nI am struggling with the premise of this work. \nMore specifically, I have questions about related work:\n\n* Line 58f states: \"All of these approaches require multiple forward solves of the forward problem\", but the proposed algorithm also requires multiple forward solves ($n$ solves of the adjoint equation for $n$ data points). For very few data points, this might be cheaper than an MCMC approach, but competing algorithms (see next point) don't seem to require multiple simulations.\n\n* It would be great if the paper discussed why one cannot/should not replace $\\Phi q$ with $H L^{-1} f$ for some, say, finite element discretisation of $L$. This seems to be a common approach in the literature on statistical finite element methods; refer to, e.g., Eqs. (21) and (22) in [1].  It is similar to SPDE approaches for Gaussian processes. I think positioning itself against this class of ideas would improve the manuscript. Or have I missed something?\n\n* Section 2 (\"Related work\") explains how an approach with Green's functions is often not feasible, because Green's functions are difficult to find and potentially unstable to apply. While I don't disagree with this statement, I would appreciate it if the experiments would include a case study in which the Green's function approach cannot be used: the ODE example seems to fit the setup from Alvarez et al.; the advection-diffusion equation appears to also admit a Green's function (see e.g., Eq. 4-3-13 in [2]; however, its representation is quite complicated)\n\n\n\n\n**References**\n\n[1] Girolami, Mark, et al. \"The statistical finite element method (statFEM) for coherent synthesis of observation data and model predictions.\" Computer Methods in Applied Mechanics and Engineering 375 (2021): 113533.\n\n\n[2] Xu, Zhanjie, Wolfgang Breitung, and J. R. Travis. Green's function method and its application to verification of diffusion models of GASFLOW Code. Karlsruhe, Germany: Forschungszentrum Karlsruhe, 2007.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "QgQgF6qAehb", "reviewer": "Reviewer_pYjt", "summary": "This paper concerns statistical inference in linear systems. A typical approach uses an integral formulation of the problem via Green's function of the differential operator. The approach taken in this paper is to instead use adjoints turning the problem into a linear model and then use a reduced-rank Gaussian process for posterior approximation. In more detail, once the adjoint of the linear operator is found, n (number of observations) equations composing a linear system are set up and the unknown parameters are solved for using the standard least squares approach.\n", "questions": "In line 256, you say \"compute the posterior exactly\" but doesn't the adjoint method rely on a truncated basis expansion? \n\nThe adjoints you have used in your models are easy to derive analytically, are there any cases where it is trickier doing this? Are there any models where you can't use the adjoint method and have to use MCMC? I don't have a good understanding of the scope of this new method compared to the old ones.\n\nWhy are adjoint solves so fast compared to the forward problem, or at least it seems this way? The use of a single MH step requires two (or one, if the previous one is saved) likelihood evaluations. At the same time in section 4.3 you mention that you did posterior inference for 100 observations, i.e. requiring 100 adjoint solves, in the time it took you to do a single MH step. Why is your MH so inefficient? I understand that the cost of the adjoint solves, once paid, leads to a nice linear system setup and posterior approximation and this may well win over MCMC when the MCMC has to be run for a long time - which is pretty much always. However, I don't think the comparison should be this optimistic.\n\nFinally, you say that M = 10 leads to a poor approximation, but you run your Metropolis sampler with it?\n", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I think this is a very nice paper, and seems like a clean approach to inference in the class of problems the paper addresses. The solutions produced enjoy stability and are conceptually easier to work with and require no sophisticated MCMC methods. In this sense, I believe it is a worthwhile addition to the literature. I am generally happy with the paper, with the only note is that the MCMC comparisons are not done particularly carefully, even so, there is enough evidence in favour of the presented method at least on the demo problems. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "L8Kz4w7AqV", "reviewer": "Reviewer_jfqV", "summary": "This paper is concerned with the inverse problem of estimating the forcing function of f of a linear system Lu=f, with affine linear observations, where the prior on f is a Gaussian process. The authors introduce a trick. They propose to use the adjoint linear system L* v = h in order to rewrite the observation vector as a linear transformation of the parameters of f (plus noise). This of course assumes that f can be decomposed into such an expansion. If f is modeled as a GP, then a truncated basis expansion can play this role - such that the posterior can be inferred.  \nThis trick has the advantage of turning the model into a standard linear model. Now, standard results can be used to compute least squares (ML) or Bayesian estimators. This is advantageous because these estimators are now only simple linear algebra. The authors derive the trick for Hilbert spaces, but claim that it can be executed for more general Banach spaces.\nIn the experiments, this approach is spelled out both for ODEs and PDEs. A real time cost analysis is given, comparing with MCMC. Which approach is best, will depend on n. The authors discuss that, if n is not huge, their adjoint approach could outperform MCMC.", "questions": "1. Table 1: It is curious that there is no difference between MCMC and Adjoint in this data. Why is that?\n2. The linear-parameter assumption of Eq. (8) is common for ODE inverse problems. See e.g. Assumption 1 in https://arxiv.org/pdf/2002.09301.pdf and Eq. (10) in https://proceedings.neurips.cc/paper/2017/file/e71e5cd119bbc5797164fb0cd7fd94a4-Paper.pdf . Maybe comment on that? In any case, it would be good for the reader to know how much of a restriction that assumption entails. What else does the existing literature propose to do with such an assumption?\n3. Lines 141-151: Why is the solution to the adjoint equation v_i not part of the ML and Bayesian estimators.", "limitations": "There are no limitations and potential negative social consequences.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength: \n1. The paper is easy to read and straightforward. The authors do a good job of highlighting the core of their paper -- how the adjoint trick turns the system into a standard linear model, solvable by linear algebra. The derivation of the adjoint equation reads quite nicely.\n2. The literature review is comprehensive.\n3. The experiments are instructive and interesting. \n4. The discussion is honest and transparent over the trade-offs in the variable n, compared to MCMC.\n5. The method is novel, as far as I know. \n6. The work is significant because (Bayesian) inverse problems are of great interest to the (probabilistic) ML community.\n\nWeaknesses: \n1. A few more experiments could have been given. Right now there is just one ODE and one PDE. More experiments would seem appropriate.\n2. The authors compare with none of the methods cited in the Related Work section. It would be great to compare with some of the most recent methods.\n3. Comments on the Computational Complexity are spread out over the paper. It would be good to write one section where all of these issues are discussed.\n4. In line 152, the authors claim that their trick extends to \"more general Banach spaces\". I don't understand how. More justification for this claim should be given; perhaps in an Appendix.\n5. Section 3.3. is a bit difficult to contextualize. I assume it shows how Bayesian inference of f works, when it is a GP. The crucial bit seems to be that the GP prior can be expressed by a truncated random series, similar to Eq. (8). If so, this connection should be more explicitly made - to improve the clarity of the writing.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept