
\(f\)  is modelled as a zero-mean Gaussian process, with kernel function \(k_\psi (\cdot , \cdot )\) , in the case where \( is an infinite dimensional Banach space, or with a Gaussian distribution in the finite dimensional case. E.g., if \) L\( is an ordinary differential operator with independent variable \) t\(, then \) f\( will be an unknown function of \) t\(.\item The linear operator \) L\( depends (possibly non-linearly) on unknown parameters \) p\(. We write \) Lp\( to emphasise this dependence.\item The observation operator \) h:URn\( is linear (or affine). In the finite dimensional case, \) h(u)=Hu\( for some matrix \) H\(, but in the infinite dimensional case (i.e., when \) U\( is a space of functions) this includes pointwise evaluation, as well as integral and derivative observations.\item The observation error \)\( has a Gaussian distribution. This assumption can be relaxed for maximum likelihood estimation.\)
The full specification of the statistical model is then
\({\mathcal {L}}_p u&=f\quad \mbox{(possibly including ICs/BCs)} \\z&=h(u)+\epsilon \\f&\sim GP(0, k_\psi (\cdot , \cdot ))\\\epsilon &\sim \mathcal {N}_n(\mu , \sigma ^2 I_n).\) 
Our aim is to infer \(f\)  (and possibly \(p\) , \(u\) , and \(\psi \) ) given \(z\)  either via

maximum likelihood (ML) estimation, by solving the constrained optimization problem
\(\min _{p, f} \quad &(z-h(u))^\top (z-h(u))+||f||^2_{eqn:sos}\\\mbox{ s.t. }\;\; &{\mathcal {L}}_p u=f\nonumber \) 

Bayesian inference, where we find
\(\pi (p, f, u | z) \propto \pi (z|u) \pi (u | f,p) \pi (f) \pi (p).\) 
Here, \(\pi (p)\)  and \(\pi (f)\)  denote the prior distributions for \(p\)  and \(f\) .

To solve either problem numerically is likely to require many solves of the forward problem (). For example, for ML we may seek a solution using numerical optimization [1]}, whereas with Bayes, we might use an MCMC scheme [2]}[3]}[4]} or variational approach [5]}. All of these approaches require multiple solves of the forward problem. Our aim is to reduce the number of forward solves required during inference, resulting in lower computational costs.
Contribution
In this paper we show that implementing an adjoint of the linear system can result in much faster statistical inference for this problem. Instead of using numerical approaches to solve either the ML or Bayesian inference problem, we can do inference for \(f\)  at the cost of \(n\)  forward model solves, where \(n\)  is the number of data points. In many (but not all) cases, this will incur a substantially lower computational cost than competing methods, such as MCMC. More specifically, we show that

if \(f\)  depends linearly on parameters \(q\) , we can estimate \(q\)  or its distribution analytically, i.e., without resorting to numerical integration methods

if we model \(f\)  as a Gaussian process, then by using a 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 REF . We derive the main results in Section REF , and in Section REF  we show how linearizing Gaussian processes via a basis expansion reduces inference for GPs in linear systems to simple linear algebra. Finally, in Section  we demonstrate the approach on two linear systems: ordinary (ODE) and partial differential equations (PDE).

Related work
Source term estimation in differential equation models has been extensively studied, for example, in the field of modelling atmospheric advection-diffusion [6]}, [7]}, [8]}, [9]}, [10]}. with many papers solving the inverse problem using an adjoint approach combined with MCMC to compute the posterior distribution [6]}, [12]}, [13]}, [4]}. Of particular relevance is [12]} who use a 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 still is still a limiting factor, restricting the extension of the approach to more complex situations such as time-varying pollution sources.
The problem defined by Eqs. (REF )-() is often referred to as a latent force model [16]}, [17]}. [16]} showed how the posterior for this model can be computed by using the integral formulation of \({\mathcal {L}}_p u = f\) , \(u(x)=\int G_p(x-v)f(v)dv\) , where \(G(\cdot )\)  is the Green's function associated with the differential operator \({\mathcal {L}}_p\) . Due to the linearity of the integral transform, placing a GP prior over \(f\)  leads to a joint GP over \(f\)  and \(u\) . From this joint GP, the posterior distributions \(\pi (u \mid f, p)\)  and \(\pi (u \mid p)\)  can be computed in closed form.
Unfortunately, in many situations, particularly for non-trivial differential equation models, the expressions for the covariances are cumbersome and lead to the use of error functions with complex arguments or functions like the Faddeeva function which can be numerically unstable to compute. Recently, [19]} proposed representing \(f\)  using random Fourier features to reduce the number of integrations necessary to be solved analytically. Here, rather than using Green's functions, we instead use adjoints to write the problem as a linear model and then combining this with a reduced-rank Gaussian process formulation, leading to numerically stable and fast approximations to the posterior distribution.
Other related work includes the stochastic PDE approach of [20]}, in which the underlying function \(u\)  is modelled as a Gauss Markov random field (GMRF), which can then be formulated as a stochastic partial differential equation. This allows finite element methods to be used to efficiently compute the posterior distribution for \(u\)  given \(z\) . Similarly, [21]} exploit the link between GMRFs and dynamical systems to convert inference for \(u\)  to a form in which Kalman filtering methods can be used, which scale linearly with \(n\) . [22]} focus exclusively on the advection diffusion PDE considered in Section REF , and use white noise for the forcing function, \(f\) , to create a stochastic PDE model for \(u\) . They use a spectral approach, 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.
[23]} show how to infer statistical inference for linearly constrained systems, i.e., where \({\mathcal {L}}u =0\) .
These approaches all model \(u\) , whereas our focus is on inferring the forcing function \(f\)  to the system. There is likely to be benefit in combining these approaches with ours to further accelerate inference in these models.

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.
Adjoints
Let \(\mathcal {U}^*\)  and \(*\)  denote the dual spaces of \(\mathcal {U}\)  and \(. We can construct the adjoint to a continuous linear operator \) Lp\( as follows. Let \) v* *\( and define \) F:UR\( by\begin{equation}F: u \mapsto v^*({\mathcal {L}}_p u).\end{equation}Then \) F\( is a bounded linear functional on \) U\(, i.e., \) F=u*\( for some \) u*U*\(. Thus for all \) v**\( we^{\prime }ve associated a unique \) u*U*\(.\begin{equation}{\mathcal {L}^*_p}: v^* \mapsto u^*=v^*\circ {\mathcal {L}}_p\end{equation}We call \) L*p\( the {\it adjoint} of \) Lp\(, and \) L*p\( is itself a bounded linear operator \cite {estep2004short}. By construction, we have that for all \) uU\( and \) v* *\(\begin{equation}({\mathcal {L}^*_p}v^* )(u) = v^* ({\mathcal {L}}_p(u)),\end{equation}a result known as the {\it bilinear identity}.\) In the case where \(\mathcal {U}\)  and \( are also real Hilbert spaces, with respect to inner products \) , U\( and \) , , then we can identify the dual spaces with their underlying space: by the Riesz representation theorem if \(v^*\in *\)  then there exists \(v\in  such that \) v*()=, v. In this case, the bilinear identity reduces to its more familiar form
\(\langle {\mathcal {L}}_p u, v\rangle _= v^*({\mathcal {L}}_p u)=({\mathcal {L}}_p^*v^*)(u)\\&= \langle u, {\mathcal {L}}_p^* v\rangle _\mathcal {U}\) 

where we now consider \({\mathcal {L}^*_p}: \mathcal {U}\) . We only consider real vector spaces here, resulting in a symmetric inner product.
Generally, the adjoint \({\mathcal {L}}_p^*\)  will be the same type of operator as \({\mathcal {L}}_p\)  (e.g. if \({\mathcal {L}}_p\)  is a differential operator then \({\mathcal {L}}_p^*\)  will be too), so solving an adjoint system of the form \({\mathcal {L}}^*_p v=g\)  will have similar computational complexity as solving \({\mathcal {L}}_p u=f\) . See [24]} for an introduction to adjoints, and the experiments section for examples of adjoints.

Benefits of adjoints
How does the development of an adjoint to a linear system help us perform statistical inference for that system? There are two main advantages. Consider first the situation where uncertainty about the unknown forcing function, \(f\)  in Eq. (), can be characterized by a linear dependence upon unknown parameters \(q\) . That is, we can write
\(f(\cdot )= \sum _{m=1}^M q_m \phi _m(\cdot ).\) 
In the infinite dimensional case where \(\mathcal {U}\)  and \( are spaces of functions on some set \) X\(, the\) m\( will also be functions on \) X\(. In the finite-dimensional case, the \) m\( will be vectors of length \) n\(.\) In the situation where the observation operator (REF ) is linear, then we can write the \(i^{th}\)  observation as \(h_i(u)=\langle h_i, u \rangle \)  plus noise, for some \(h_i \in \mathcal {U}\) . Consider the \(n\)  different adjoint systems
\({\mathcal {L}^*_p}v_i = h_i \mbox{ for }i=1, \ldots , n.\) 
Then using the bilinear identity (REF ) we get
\(h_i(u)=\langle h_i, u \rangle &=\langle {\mathcal {L}^*_p}v_i, u\rangle \\&= \langle v_i, {\mathcal {L}}_p u\rangle \\&= \langle v_i, f\rangle ,\) 

i.e., the \(i^{th}\)  observation is the inner product between the unknown forcing function \(f\)  and the solution of the \(i^{th}\)  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: an \(n\) -fold increase in computational cost. The benefit arises if we now use the assumption of a linear dependence upon the parameters, and linearity of inner products:
\(h_i(u) = \langle v_i, \sum _{m=1}^M q_m \phi _m \rangle =\sum _{m=1}^M q_m \langle v_i, \phi _m\rangle .\) 
The complete observation vector \(z\)  can then be written as
\(z &= \begin{pmatrix} \langle v_1, \phi _1 \rangle & \ldots & \langle v_1, \phi _M \rangle \\\vdots &&\vdots \\\langle v_n, \phi _1 \rangle & \ldots & \langle v_n, \phi _M \rangle \end{pmatrix}\begin{pmatrix} q_1\\\\q_M\end{pmatrix}+\epsilon \\&=\Phi q+\epsilon \nonumber \) 

where \(q\in \mathbb {R}^M\)  is the parameter vector, and
\(\Phi \in \mathbb {R}^{n\times M}\)  is the matrix of inner products between the \(n\)  adjoint solutions and \(M\)  basis vectors.
This can be recognized as a standard linear model. Thus, standard results can be used to compute the least squares/maximum likelihood and Bayesian estimators. For maximum likelihood, the minimum of \(S(q) = (z-h(u))^\top (z-h(u))\)  subject to \({\mathcal {L}}_p u=f\) ,
is obtained at
\(\hat{q} = (\Phi ^\top \Phi )^{-1} \Phi ^\top z\) 
with
\(\mathbb {V}\mbox{ar}(\hat{q}) = \sigma ^2(\Phi ^\top \Phi )^{-1}\)  in the case where the observation errors \(\epsilon _i\)  are uncorrelated and homoscedastic with variance \(\sigma ^2\) . Standard results can be used from regularized least squares if we need to regularize \(q\) .
In a Bayesian setting, if we assume a priori that \(q\sim \mathcal {N}_M(\mu _0, \Sigma _0)\) , then the posterior for \(q\)  given \(z\)  (and other parameters) is
\(q \mid z \sim \mathcal {N}_M(\mu _n, \Sigma _n)\) 
where
\(\mu _n = \Sigma _n (\frac{1}{\sigma ^2}\Phi ^\top z +\Sigma _0^{-1} \mu _0), \;\; \Sigma _n=\left(\frac{1}{\sigma ^2}\Phi ^\top \Phi + \Sigma _0^{-1}\right)^{-1}.\) 
See, e.g., [25]}.
Note that for readability, we have just presented the approach for Hilbert spaces here, but the idea works in more general Banach spaces.
A second advantage of using adjoints is that the gradient of a cost function depending on the solution \(u\)  can be computed with respect to parameters \(p\) , using the adjoint sensitivity approach [26]}, [27]}, [28]}.
For example, consider the unregularized sum of squares/log-likelihood \(S(p,q) = (z-h(u))^\top (z-h(u))\)  from Eq. (), where \(u\)  satisfies \({\mathcal {L}}_p u =f_q\) . We've shown how we can estimate \(\hat{q}(p)=\arg \max S(p,q)\)  or find the posterior \(\pi (q\mid z, p)\) , but using the adjoint sensitivity approach also gives us \(\frac{\partial }{\partial p}S(p,q)\) . This can be used either by a gradient-based optimizer to learn \(\hat{p}=\arg \max S(p, \hat{q}(p))\) , or in a variational inference scheme that requires gradients to estimate the variational parameters such as the variational autoencoder [29]}, [30]}, or in a Hamiltonian Monte Carlo algorithm [31]} to efficiently target the posterior (Eq. REF ). We do not explore these approaches in detail here.

Inference of Gaussian forcing functions
We now consider the situation where the forcing function is given a Gaussian process prior distribution:
\(f(\cdot )\sim GP(m(\cdot ),k(\cdot ,\cdot )),\) 
where \(m(\cdot )\)  and \(k(\cdot ,\cdot )\)  are the prior mean and covariance functions respectively [32]}. Our approach is to use a reduced-rank representation of the GP, as in Eq. (REF ), derived by truncating an expansion for \(k\)  of the form
\(k(x,x^{\prime })=\sum _{m=1}^{\infty }\phi _m(x)\phi _m(x^{\prime }).\) 
There are many possible choices for the basis vectors \(\phi _i(\cdot )\) , including the Karhunen-Loève (KL) [33]}, Laplacian [34]}, [35]}, [36]}, and random Fourier feature [37]} expansions.
The Karhunen Loève basis is formed by finding the spectral expansion of the integral operator defined in Mercer's theorem:
\(T_k f(x) = \int _\mathcal {X} k(x,x^{\prime }) f(x^{\prime }) \mathrm {d}x^{\prime }.\) 
In the case where \(\dim (x)\)  is small (such as the ODE example below), the eigenfunctions of \(T_k\)  are easy to compute either analytically [32]} or numerically [39]}. Using the eigenfunctions of \(T_k\)  gives the \(L^2\) -optimal approximation [40]}, 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 random Fourier features (RFFs) [37]}. This relies upon Bochner's theorem, which can be used to express stationary kernels \(k(x,x^{\prime }):=k(x-x^{\prime })\)  as the Fourier transform of a positive measure \(p\)
\(k(x-x^{\prime })=\int e^{-iw(x-x^{\prime })}\mathrm {d}p(w).\) 
If we use an isotropic exponentiated quadratic kernel,
\(k(x-x^{\prime })= \tau ^2 \exp (-\frac{1}{2\lambda ^2}(x-x^{\prime })^\top (x-x^{\prime }))\) 
then the measure corresponds to a multivariate Gaussian. Using results from [37]}, we can approximate the Gaussian process with Fourier feature vectors:
\(\phi _i(x) = \sqrt{\frac{2\tau ^2}{M}} \cos ( \frac{1}{\lambda }w_i^\top x+b_i)\) 
where
\(w_i\sim N(0,1)\)  and
\(b_i \sim U(0, 2\pi )\) . The extension anisotropic kernels is straight-forward.
<FIGURE>Although the RFF expansion will requires 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. REF ), and so including more terms has a minor effect on overall computational cost.
<FIGURE>
Experiments
Ordinary differential equation (ODE)
Consider the non-homogeneous linear ODE:
\(p_2\frac{d^2u}{dt^2}+p_1\frac{du}{dt}+p_0 u=f(t) \) 
on the domain \([0, T]\)  with initial conditions \(u(0)=u^{\prime }(0)=0\) . The right hand side, \(f(t)\) , is the unknown forcing function that we wish to estimate, and \(p=(p_0, p_1, p_2)^\top \)  are parameters in the linear operator. We model \(f\)  as a zero-mean GP with the exponentiated quadratic covariance function:
\(f(\cdot )\sim GP(0,k(t,t^{\prime }))\) .We assume observations are obtained as noisy time averages over short time windows:
\(h_i(u) = \int _{t_i}^{t_i+\Delta t} \frac{1}{\Delta t} u(t) \mathrm {d}t=\langle u, \tilde{h}_i\rangle \) 
where \(\tilde{h}_i\)  is the indicator function \(\tilde{h}_i(t) = \mathbb {I}_{[t_i, t_i+\Delta t]}(t)\) . We could include direct measurements of \(u(t_i)\)  by setting \(\tilde{h}_i(t) = \delta (t-t_i)\) . To generate synthetic data, we first simulate a realization \(f\)  from the GP model, solve Eq. (REF ) for \(u(t)\) , and then simulate \(n\)  observations from Eq. (REF ) with \(T=1\) , \(p_2=0.5, p_1=1, p_0=5\) , \(\Delta t = \frac{T}{n}\) , \(t_i = \frac{iT}{n}\) , adding zero-mean Gaussian noise with standard deviation \(\sigma =0.1\) . For simplicity, we solve the ODE with a simple forward Euler approximation, but higher order schemes can and should be used in real applications (although see the discussion for comments about numerical difficulties). We approximate the GP using Eq. (REF ), using 200 RFFs generated using Eq. (REF ) with \(\lambda =\sqrt{0.6}\)  and \( \tau ^2=4\) . The linear operator in this case is
\({\mathcal {L}}_p u = \left(p_2 \frac{\mathrm {d}^2}{\mathrm {d}t^2} + p_1 \frac{\mathrm {d}}{\mathrm {d}t} + p_0\right)u.\) 
To derive the adjoint operator we use the bilinear identity (Eq. REF ), and integrate by parts twice:
\(\langle \mathcal {L}u, v \rangle &= \int _0^T \mathcal {L}u(t) v(t) \mathrm {d}t= \int _0^T (p_2 \ddot{u} + p_1 \dot{u} + p_0 u)v \mathrm {d}t\\&=\left[p_2\dot{u} v+p_1 u v\right]_0^T+\int _0^T -p_2\dot{u}\dot{v} - p_1 u\dot{v} + p_0 u v \mathrm {d}t\\&= \left[p_2(\dot{u} v+u\dot{v})+p_1 u v\right]_0^T+\int _0^T p_2 u\ddot{v} - p_1 u\dot{v} + uv \mathrm {d}t\\&= \int _0^T (p_2\ddot{v} - p_1\dot{v}+v)u \mathrm {d}t =\langle u, \mathcal {L}^*v \rangle \) 

when \(v(T)=\dot{v}(T)=0\) . So the adjoint of \({\mathcal {L}}_p\)  is
\({\mathcal {L}^*_p}v = \left(p_2 \frac{\mathrm {d}^2}{\mathrm {d}t^2} - p_1 \frac{\mathrm {d}}{\mathrm {d}t} + p_0\right)v.\) 
Note that rather than an initial condition, the adjoint has a final condition, and so to solve the system we have to integrate it backwards in time from \(t=T\)  to \(t=0\) .
A summary of the approach is given in Algorithm REF .
[ht]
Adjoint-aided inference for the ODE model

[0]
For \(i=1,\ldots , n\) :

Solve the adjoint system \({\mathcal {L}^*_p}v_i = \tilde{h}_i\)  with final condition \(v_i(T)=\dot{v}_i(0)=0.\)

For \(m=1, \ldots , M\) :

Sample an RFF basis vector \(\phi _m(t)\)  using Eq. (REF )

Compute \(\langle v_i, \phi _m\rangle \)  for \(1\le i\le n\)

Form the matrix \(\Phi \) , and compute the posterior distribution for \(q\)  using Eqs. (REF ) and (REF ).
Fig. REF  shows the posterior distribution of the forcing function when using \(n=50\)  and \(n=500\)  data points. As expected, more data results in a more confident posterior.
Fig. REF  shows the effect of increasing the number of features in the GP approximation. Note the danger of using too few features: with \(M=10\)  the approximation to \(f\)  has limited expressive power and cannot capture the true form of \(f\) , i.e., the model is heavily misspecified. This can result in the uncertainty collapsing upon the most likely, but wrong, value, which can be difficult to spot in real applications.

Partial differential equation (PDE)
We now demonstrate the approach on a PDE in which there are two spatial variables, \(x\in \mathcal {X}\subset \mathbb {R}^2\) , and time, \(t\in [0,T]\subset \mathbb {R}\) , so that the solution \(u\equiv u(x,t)\)  is a function of three independent variables. We consider the advection-diffusion equation:
\(\frac{\partial u}{\partial t} + p_1\nabla u -\nabla .(p_2 \nabla u)=f \mbox{ in } \mathcal {X}\times [0,T]\) 
with the following initial and boundary conditions:
\(u(x,0)&=0 \mbox{ for } x\in \mathcal {X}\\\nabla _{n}u&=0 \mbox{ for }x\in \partial \mathcal {X}.\) 

Here, the unknown forcing function \(f\equiv f(x,t)\)  is a function of space and time, and we model it as a zero-mean GP
\(f(x,t) \sim GP(0,k((x,t),(x^{\prime },t^{\prime }))\) 
with exponentiated quadratic kernel \(k\)  (Eq. REF ). We use a random Fourier feature approximation to \(k\)  (Eq. REF ), where in this case the random weights are \(w_m\sim \mathcal {N}_3(0, I)\)  and \(b_m \sim U[0,2\pi ]\) .
Adjoint-aided inference for the PDE model

[0]
For \(i=1,\ldots , n\) :

Solve the adjoint system \({\mathcal {L}^*_p}v_i = \tilde{h}_i\)  with final and boundary conditions:
\(v(x,T)&=0 \mbox{ for }x\in \mathcal {X}\\p_1v+p_2\nabla v&=0 \mbox{ for } x\in \partial \Omega \mbox{ and } t\in [0,T].\) 

For \(m=1, \ldots , M\) :

Sample an RFF basis vector \(\phi _m(t)\)  using Eq. (REF )

Compute \(\langle v_i, \phi _m\rangle \)  for \(1\le i\le n\) .

Form the matrix \(\Phi \) , and compute the posterior distribution for \(q\)  using Eqs. (REF ) and (REF ).
Observations are assumed to arise from sensors which take an average of \(u(x,t)\)  over a small spatial and temporal window \(\mathcal {R}_i \times \mathcal {T}_i\subset \mathcal {X}\times [0,T]\) :
\(z_i &= \langle u, \tilde{h}_i\rangle + e_i\\\mbox{with}\quad \tilde{h}_i &= {\left\lbrace \begin{array}{ll}\frac{1}{|\mathcal {R}_i|.|\tau _i|} & \mbox{ if } x\in \mathcal {R}_i \mbox{ and } t\in \mathcal {T}_i \\0 &\mbox{ otherwise}.\end{array}\right.}\) 

The adjoint of the linear operator \({\mathcal {L}}_p u = (\frac{\partial u}{\partial t} + p_1\nabla u -\nabla .(p_2 \nabla u)\) 
is
\({\mathcal {L}^*_p}v=-\frac{\partial v}{\partial t}-p_1\nabla v-\nabla .(p_2\nabla v).\) 
This can be derived by integrating by parts as for the ODE, but for partial differential operators, we rely upon Green's theorem to do so; see [43]} for details. Our adjoint-aided approach then requires the solution of
\({\mathcal {L}^*_p}v_i = \tilde{h}_i \mbox{ in } \mathcal {X} \times [0,T]\) 
with the following initial and boundary conditions:
\(v(x,T)&=0 \mbox{ for }x\in \mathcal {X}\\p_1v+p_2\nabla v&=0 \mbox{ for } x\in \partial \Omega \mbox{ and } t\in [0,T]\) 

for \(i = 1, \ldots , n\) . Inference then proceeds as before. After solving the \(n\)  adjoint systems, we compute the inner product of these solutions with the random Fourier feature basis vectors to form the matrix \(\Phi \)  (Eq. REF ). We can then use Eqs. (REF ) and (REF ) to compute the posterior with minimal additional computational cost. The approach is summarized in Algorithm REF .
Data was simulated on the spatial domain \(\mathcal {X}=[0,10]^2\)  for \(t\in [0,10]\)  by first randomly generating a forcing function \(f(x,t)\) , and then solving the forward problem (Eqs. REF -REF ) to find \(u(x,t)\) . We generate \(n\)  observations using sensors that record averages over short time windows equally spaced across the domain \([0,10]\) , and at the locations shown in Fig. REF . Zero-mean Gaussian distributed noise is added to the true sensor readings with standard deviation \(\sigma =0.05\)  (note that this is relatively small compared to the signal, which can often create problems for sampling methods).
We then use Algorithm REF  to calculate the posterior distribution for \(q\) , hence giving the posterior for \(f\) . By sampling forcing functions from this posterior and simulating forward, we can test the posterior predictive accuracy of the model.
Fig. REF  shows the ground-truth forcing function \(f\)  (generated from an isotropic GP with an exponentiated quadratic kernel with \(\lambda =2\) , \(\tau ^2=2\) ) as a spatial map at three different time points. Observations were generated for four time windows using both arrays of 4 and 16 sensors at the locations shown in Fig. REF , i.e., a total of either \(n=16\)  or \(n=64\)  observations.
We used \(M=200\)  RFFs to infer \(f\)  using the adjoint-aided approach. The posterior mean for \(f\)  is shown as spatial maps in Figs. REF  and REF  at 3 different time slices. As expected, more sensors results in better estimates.
<FIGURE>Figure REF  shows the posterior standard deviation for \(f\)  in the 4 and 16 sensor cases. Here we used PDE parameters \(p_1=(0.4,0.4)^\top \)  (and \(p_2=0.01\) ), so that advection occurs parallel to the \(y=x\)  line. We can see that standard deviation is greatest downwind of the sensors, with a reduction in the variance upwind, as expected.
<FIGURE>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 to calculate the posterior predictive mean squared error (MSE). Figure REF  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 to 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 to investigate these effects in a 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 \(n\)  evaluation required for the adjoint-aided approach).
<FIGURE>Finally, we can infer the remaining parameters (PDE and GP hyperparameters) in a variety of ways. By way of illustration,
Figure REF  shows the results from using Bayesian optimization to infer lengthscale in a problem where data are generated with true lengthscale \(\lambda =2\) . We used GPyOpt to maximise the negative log-likelihood using the expected improvement acquisition function. In this case, we are easily able to determine the correct lengthscale with relatively minimal computational cost.
<FIGURE>
Discussion
This work was motivated by the problem of inferring (with uncertainty) the distribution of spatially and temporally varying sources of pollution (e.g. across a city), given a small set of (noisy) observations, and a model of the pollution's atmospheric transport.
This is typical of the type of challenge faced throughout the sciences and engineering. Linear systems such as this, still pose computational challenges that make any form of principled statistical inference almost intractable, or if they are tractable, the computational cost (e.g. of using MCMC) is such that only a limited range of models and situations can be analysed. This paper combines the adjoint approach, which allows us to write the problem as a standard linear model, with a linear-sum formulation of a Gaussian process for the forcing term. The result is conjugate Bayesian inference for an unknown GP forcing function, which can vary in both time and space. The computational cost scales linearlySpecifically, with \(G\)  grid elements and \(M\)  features, the algorithm has computational complexity \(O(GMn + M^2 n + M^3)\) . in time 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^5-10^6\)  iterations even for low-dimensional problems, there is still a large range of problems for which an adjoint may be beneficial, more so if it used to also compute gradients. In our PDE example with relatively few observations, the adjoint-aided approach requires orders of magnitude fewer PDE solves than MCMC. It is also worth noting that the computational tasks in our approach are embarrassingly parallelisable.
One issue that is only briefly touched upon in this paper, is the problem of estimating the non-linear operator parameters \(p\) , such as advection and diffusion rates, as well as the GP kernel hyperparameters. Using adjoint sensitivity to estimate gradients with respect to these parameters is a numerically stable alternative to using auto-diff frameworks (such as TensorFlow), which can be unstable when back propagating through long iterative loops. Knowledge of gradients then allows inference for the additional parameters to be performed efficiently within the preferred statistical paradigm. A cautionary note is that when the observation operator \(\tilde{h}\)  has small support, such as for pointwise evaluation (i.e., \(h(u) = u(x,t)\) ), then care is needed when using numerical methods to solve the adjoint systems, as adaptive step size methods can easily miss small non-zero regions in \(\tilde{h}\) .
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 use varying grid sizes, and stochastic approaches which use only a subset of the data at each stage.

Acknowledgements
This work is funded by EPSRC/GCRF grant EP/T00343X/1, with additional support provided by GOOGLE.ORG. RW is supported by EPSRC projects EP/P010741/1 and EP/S014985/1.

Supplementary material

PDE Inference Examples
Here we show the effect of the GP lengthscale, \(\lambda \) , on the quality of inference in the PDE problem. Figures REF  to REF  show the ground truth forcing function, \(f(x,t)\) , simulated concentration, \(u\) , and the forcing function inferred with 200 random Fourier features, in the cases where \(\lambda =0.2\) , 1 and 3, respectively. It should be noted that at \(t=0\) , \(u\)  is 0 everywhere, due to the established initial conditions. In the shorter lengthscale cases, the quality of inference is worse than that at \(\lambda =3\)  (and at \(\lambda =2\) , as shown in main paper). However, the quality of inference in shorter lengthscale cases can be improved by increasing the number of random Fourier features used to perform the adjoint method.
<FIGURE><FIGURE><FIGURE>
Time complexity
Let \(G\)  be the number of grid elements, \(n\)  the number of observations and \(M\)  the number of features. There are five specific operations our algorithm requires. (1) Solving the \(n\)  adjoint systems takes \(O(Gn)\) . (2) Computing each \(\phi \)  vector over the grid for each feature takes \(O(GM)\) . (3) Computing the \(\Phi \)  matrix, takes \(O(GMn)\)  time (having to dot product each feature computed over the grid with each \(q\)  associated with each observation). (4) Finding \(\Phi ^\top \Phi \)  takes \(O(nM^2)\) . Finally solving the matrix inverse will take \(O(M^3)\) . The upshot is an algorithm linear in \(n\) . Empirically, for the problems we've experimented with, we found computing \(\phi \)  over the large grid was the most time consuming step. The overall time complexity is \(O(GMn + M^2 n + M^3)\) , but the constants associated with each term are important to consider.
