The problem defined by Eqs. (REF )-() is often referred to as a latent force model [1]}, [2]}. [1]} 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, [4]} 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.
