Unbiased Estimation of the Gradient of the Log-Likelihood in Inverse Problems

We consider the problem of estimating a parameter associated to a Bayesian inverse problem. Treating the unknown initial condition as a nuisance parameter, typically one must resort to a numerical approximation of gradient of the log-likelihood and also adopt a discretization of the problem in space and/or time. We develop a new methodology to unbiasedly estimate the gradient of the log-likelihood with respect to the unknown parameter, i.e. the expectation of the estimate has no discretization bias. Such a property is not only useful for estimation in terms of the original stochastic model of interest, but can be used in stochastic gradient algorithms which benefit from unbiased estimates. Under appropriate assumptions, we prove that our estimator is not only unbiased but of finite variance. In addition, when implemented on a single processor, we show that the cost to achieve a given level of error is comparable to multilevel Monte Carlo methods, both practically and theoretically. However, the new algorithm provides the possibility for parallel computation on arbitrarily many processors without any loss of efficiency, asymptotically. In practice, this means any precision can be achieved in a fixed, finite constant time, provided that enough processors are available.


Introduction
The problem of inferring unknown parameters associated to the solution of (partial) differential equations (PDEs) is referred to as an inverse problem. In such a context, when the forward problem is well-posed, the inverse problem is often ill-posed and challenging to solve, even numerically. The area has a long history and a large literature (see e.g. [7,14]) yet the intersection with statistics is still comparatively small, particularly considering the significant intersection, in terms of both methods and algorithms as well as objectives. If one adopts a Bayesian approach to solution of the inverse problem then the object of interest is a posterior distribution, and in particular expectations with respect to this distribution [8,12]. While this provides an elegant solution and quantified uncertainty via well-defined target distribution, it is more challenging to solve than its deterministic counterpart, requiring at least a Hessian in addition to a maximum a posteriori estimator for a Laplace approximation, if not more expensive Monte Carlo methods. Here we assume solution of the Bayesian inverse problem (BIP) requires computationally intensive Monte Carlo methods for accurate estimation. We furthermore assume that the statistical model can only be defined up to some unknown parameters.
Consider a BIP with unknown u ∈ X and data y ∈ Y, related through a PDE, and assume that the statistical model is known only up to some parameter θ ∈ Θ ⊆ R d θ (assumed finite dimensional). In other words, the posterior distribution takes the form p(du, θ|y) ∝ p(y|u, θ)p(du|θ)p(θ) .
Due to sensitivity with respect to the parameter θ and strong correlation with the unknown u, such posterior distribution can be highly complex and very challenging to sample from, even using quite advanced Markov chain Monte Carlo (MCMC) algorithms. In this article, the unknown u is treated as a nuisance parameter and the goal is to maximize the marginal likelihood of the parameters p(y|θ) = X p(y|u, θ)p(du|θ) .
In such a scenario one is left with a finite-dimensional optimization problem, albeit with an objective function that is not available analytically. This intractability arises from two sources: • first, for a given (u, θ) only a discretization of the likelihood p(y|u, θ) can be evaluated; • second, the discretized marginal likelihood is a high-dimensional integral which itself must be approximated.
Moreover, the associated gradient of the log-likelihood is not available, which may be of interest in optimization algorithms. In the following we will suppress the notation for fixed observation y and present the method generally.
In particular, we use the notation γ θ (u) = p(y|u, θ)p(u|θ), where du represents the finite measure of an infinitesimal volume element, which may or may not be Lebesgue measure, and p(u|θ) = (p(du|θ)/du)(u). We will also denote its integral Z θ = p(y|θ), and the posterior by η θ (du).
In this article, we present a new scheme to provide finite variance estimates of the gradient of the log-likelihood that are unbiased. To be precise, let E θ = ∇ θ log(Z θ ) denote the gradient of the log-likelihood with no discretization bias. The proposed method provides an estimatorÊ θ such that E[Ê θ ] = E θ , where E is the expectation with respect to the randomization induced by our numerical approach. Moreover, the estimatorÊ θ is constructed so that one only needs access to finite resolution (discretized) approximations of the BIP. This scheme is of interest for several reasons: 1. Unbiased estimates of gradients help to facilitate stochastic gradient algorithms; 2. The method is easy to parallelize; 3. The method helps to provide a benchmark for other computations.
In terms of the first point, it is often simpler to verify the validity of stochastic gradient algorithms when the estimate of the noisy functional is unbiased. Whilst this is not always needed (see [13] for a special case, which does not apply in our context), it at least provides the user a peace-of-mind when implementing optimization schemes. The second point is of interest, in terms of efficiency of application, especially relative to competing methods. The third point simply states that one can check the precision of biased methodology. We now explain the approach in a little more detail.
The method that we use is based upon a technique developed in [9]. In that article the authors consider the filtering of a class of diffusion processes, which have to be discretized. The authors develop a method which allows one to approximate the filtering distribution, unbiasedly and without any discretization error. The methodology that is used in [9] is a double randomization scheme based upon the approaches in [10,11]. The work in [10,11] provides a methodology to turn a sequence of convergent estimators into an unbiased estimator, using judicious randomization across the level of discretization. It is determined for the problem of interest in [9] that an additional randomization is required in order to derive efficient estimators, that is, estimators that are competitive with the existing state-of-the-art methods in the literature. In this article we follow the basic approach that is used in [9], except that one cannot use the same estimation methodology for the current problem. An approach is introduced in [2] which enables application of the related deterministic multilevel Monte Carlo identity [15] to a sequential Monte Carlo (SMC) sampler [5,6] for inference in the present context. In this article, we consider such a strategy to allow the application of the approach in [9] to unbiasedly estimate the gradient of the log-likelihood for BIPs. The method of [2] is one of the most efficient techniques that could be used for estimation of the gradient of the log-likelihood for BIPs. However, this method is subject to discretization bias. In other words, suppose E l θ is the gradient of the log-likelihood with a choice of discretization bias level, e.g. 2 −l . The original method would produce an estimatê E l θ for which E[Ê l θ ] = E l θ . On the other hand, under assumptions, it is proven that the new method introduced here can produce an estimate E θ with finite variance and without bias, i.e. E[Ê θ ] = E ∞ θ . We also show that the cost to achieve a given variance is very similar to the multilevel SMC (MLSMC) approach of [2], with high probability. This is confirmed in numerical simulations. We furthermore numerically investigate the utility of our new estimator in the context of stochastic gradient algorithms, where it is shown that a huge improvement in efficiency is possible. Our approach is one of the first which can in general provide unbiased and finite variance estimators of the gradient of the log-likelihood for BIPs. A possible alternative would be the approach of [1], however, the methodology in that article is not as general as is presented here and may be more challenging to implement.
This article is structured as follows. In Section 2 we explain the generic problem to which our approach is applicable. In particular, a concrete example in the context of Bayesian inverse problems is described. In Section 3 we present our methodology and the proposed estimator. In Section 4 we show that our proposed estimator is unbiased and of finite variance and we consider the cost to obtain the estimate. In Section 5 several numerical examples are presented to investigate performance of the estimator in practice, including the efficiency of the estimator when used in in the relevant context of a stochastic gradient algorithm for parameter estimation. In appendix A the proofs of some of our theoretical results can be found.
Let P P (p), p ∈ N 0 , be a positive probability mass function with P P (p) = ∞ q=p P P (q). Now if and P ∼ P P (·), then will allow (Ξ l θ ) l∈N0 to satisfy (6)- (7), where expectations are understood to be with respect to P P yet P is suppressed in the notation. Moreover (Ξ l θ ) l∈N0 will have finite variances. This result follows as we are simply using the coupled sum estimator as in [11] and using [15,Theorem 5], for instance, to verify the conditions required.

Details of the Approach
We will now describe how to obtain the sequence (Ξ l,p θ ) p∈N0 for l ∈ N 0 fixed.

MLSMC Method of [2]
To introduce our approach, we first consider the MLSMC method in [2] which will form the basis for our estimation procedure. Define for l ∈ N 0 Define for µ ∈ P(X) (the collection of probability measures on (X, X )), l ∈ N Noting that equations (14) and (15) lead to the recursion Consider N ∈ N, and slightly modify the MLSMC algorithm used in [2] to keep the number of samples across layers fixed, up to some given level l ∈ N. Details are given in Algorithm 1. This algorithm yields samples distributed according to the following joint law where One can compute an estimate of η 0 θ (ϕ 0 θ ) as Algorithm 1 A Multilevel Sequential Monte Carlo Sampler with a fixed number of samples N ∈ N and a given level l ∈ N 0 . 1. Initialization: For i ∈ {1, . . . , N } sample U i 0 from η 0 θ . If l = 0 stop; otherwise set s = 1 and go-to step 2.

Resampling and Sampling: For
). This consists of sampling a i s ∈ {1, . . . , N } with probability mass function , and then sampling U i s from M s θ (u θ ) p∈N0 In order to calculate our approximation, we will consider the following approach, which was also used in [9]. Given any (l, P ) ∈ N 2 0 , we will run Algorithm 2 in order to obtain (Ξ l,p θ ) p∈{0,1,...,P } .
The joint probability law of the samples simulated according to Algorithm 2 is where N −1 = 0 and P Np−Np−1 θ is as defined in (17). For (l, P ) ∈ N 2 0 given, consider running Algorithm 2. Then for any s ∈ {0, 1, . . . , (l − 1) ∨ 0} and any p ∈ {0, . . . , P } we can construct the following empirical probability measure on (X, X ) Note the recursion η s,N0:p θ where E θ is used to denote expectation associated to the probability P θ in (18).

Method
The new method is now presented in Algorithm 3.
For i = 1, . . . , M : 1. Generate L i ∼ P L and P i ∼ P P .
2. Run Algorithm 2 with l = L i and P = P i .

Compute:
Return the estimate: The estimate of η θ (ϕ θ ) is given by (23). In Section 4 we will prove that it is both unbiased and of finite variance, as well as to investigate the cost of computing the estimate.
There are several points of practical interest to be made at this stage (the first two were noted already in [9]). First, the loop over the number of independent samples i in Algorithm 3 can be easily parallelized. Second, one does not need to make L and P independent; this is only assumed for simplicity of presentation, but is not required. Third, the current method uses only the level l − 1 marginal of (18). All the samples for s = 0, . . . , l − 2 and associated empirical measures (19) are discarded and only the level l − 1 empirical measure is utilized. This differs from [2] where all the lower level empirical measures are used. It is possible these samples could be utilized to improve the accuracy of the method, but it is not necessary and so is not investigated further here. The potential efficiency of the double randomization scheme, as well as a discussion of the overall efficiency of the approach, is given in [9, Section 2.5].
To continue our discussion, to complete our proof, we must know something about the quantities in terms of a possible decay as a function of l. To that end, we shall assume that these terms are O(h β l ) for some β > 0. This assumption can be verified for the example in Section 2.2. Recall from Section 2.2 that h l = 2 −l . Proposition 4.2. Assume (A1). Then there is C > 0, depending on f and u * , and β > 0 such that for all u ∈ X p l (·; u) , p(·; u) < C , where the norm is L 2 (D). Given a function F : N × X → R n , suppose that there is a C > 0 which does not depend on (l, u) such that where the first norm is understood as the n−dimensional Euclidean norm, while the second norm is L 2 (D), and F (∞, ·) := lim l→∞ F (l, ·). Then there is another C > 0 which does not depend on (l, u) such that by Proposition 4.1 that for (l, s) ∈ N 2 Also, by using Cauchy-Schwarz, (l, p, q) and E θ [ Ξ l,s θ Ξ l,q θ ] are O(1) so one can find a C such that the following upper-bound holds) . Now if one chooses, for instance N p = 2 p and P P (p) ∝ 2 −p (p + 1) log 2 (p + 2) 2 and P L (l) ∝ h αβ l for any α ∈ (0, 1) then (8) is satisfied and hence the proof is completed.
In most cases of practical interest, it is not possible to choose P L , P P and (N p ) p∈N0 so that (23) is an unbiased and finite variance estimator, as well as having finite expected cost. Suppose, in the case of Section 2.2, the cost to evaluate G l θ is O(h −1 l ) and β = 1. Then, just as in [9], if we choose P L (l) ∝ 2 −l (l + 1) log 2 (l + 2) 2 , N p = 2 p , and P P (p) ∝ 2 −p (p+1) log 2 (p+2) 2 , then to achieve a variance of O( 2 ) (for > 0 arbitrary) the cost is O( −2 | log( )| 2+δ ) for any δ > 0, with high probability. For the MLSMC method in [2], the cost to obtain a mean square error of O( 2 ) is O( −2 log( ) 2 ), which is a mild reduction in cost. However, we note that this discussion is constrained to the case of a single processor. The unbiased method is straightforwardly parallelized.

Numerical Results
First we will consider a toy example where we can analytically calculate the marginal likelihood and investigate the performance of the resulting estimator in comparison to the estimator obtained using the original MLSMC method of [2] (not presented here). Subsequently we will consider the example from Section 2.2. Finally, for both examples we will explore the potential applicability of our estimators within the context of parameter optimization using the stochastic gradient descent (SGD) method.
The forward model is the same for both problems, hence the anticipated rate of convergence is the same, and is estimated as β = 4, just as in [2]. The cost would be O(h −γ l ) in general for the problem in Section 2.2, and for the particular setup described in 2.2.1 γ = 1. We choose h l = 2 −l and P L (l) ∝ 2 −2.5l . This is far into the so-called canonical regime (β > γ), and therefore we allow unbounded L, i.e. L max = ∞ in the terminology of [9]. The reason for this is basically that the sum (8) and the corresponding cost series both easily converge, if the cost is deterministic and O(h −γ l ) as a function of h l . However, in this case the cost depends upon the randomized estimator of the series in p. Since the rate of convergence is borderline in the p direction, β p = 1 = γ p , as in [9] we impose a maximum P max on P . This is necessary to prevent the possibility of the algorithm getting stuck with an extremely expensive sample. It is discussed further in that work. In particular, we choose N p = 2 p+3 and The piecewise definition of P P ensures that it has the correct asymptotic behaviour but is also monotonically decreasing. Note that in this regime, i.e. strongly canonical convergence in L, or large β > γ, the MLSMC method easily achieves the optimal complexity O( −2 ). However, since the convergence rate in P is necessarily subcanonical, our method therefore suffers from a logarithmic penalty, i.e. O( −2 log( ) 2+δ ), for any δ > 0. This cannot be observed in the simulations though. Empirically we observe that we can set P max rather small, which is perhaps afforded by the very fast convergence in the L direction. This may also be why we cannot see the theoretically predicted log penalty in the simulations.

Toy example
We first consider an example where the marginal likelihood is analytically calculable. Consider the following PDE on D: . Define the observation operator as Suppose the observation takes the form y = G(u) + ξ, ξ ∼ N (0, θ −1 · I M ), ξ ⊥ u, and θ follows a log-normal distribution with mean 0 and variance σ 2 . Then the unnormalized likelihood is given by and the marginal likelihood is So the log-likelihood is then given by and the derivative of the log-likelihood is First, the performance of the unbiased algorithm for a single gradient estimation The data is generated with M=50, and observation operator G(u) = [p(x 1 ; u), p(x 2 ; u), . . . , p(x M ; u)] with x i = i/(M + 1). θ * is set to be 2, and the true value of the derivative of the log-likelihood at this point is calculated using the above analytical solution. For each L, the MLSMC estimator is realized 50 times and the MSE is reported. Similarly, the MSE of unbiased algorithm is calculated based on 50 realizations. The results are presented in Figure  1. The cost reported in the plot is proportional to the sum of the cost per forward solve at level l (tridiagonal linear system), h −1 l , multiplied by the total number of samples at level l.

Example of Sec. 2.2
Now that we understand the behaviour of the estimator we consider the example from Sec. 2.2. Here we do not have an analytical solution, so the true value of the target was first estimated with the MLSMC algorithm with L=12. This sampler was realized 50 times and the average of the estimator is taken as the ground truth. Now for each L, the MLSMC estimator is realized 50 times and the MSE is reported. Similarly, the MSE of unbiased algorithm is calculated based on 50 realizations. The results are presented in Figure 2. The cost in the plot is proportional to the sum of the cost per forward solve at level l (tridiagonal linear system), h −1 l , multiplied by the total number of samples at level l.

Stochastic gradient descent method
In this section, we investigate the potential to use our unbiased estimators within the SGD method. Recall we want to find the MLE of θ by minimizing − log p(y|θ) = − log( γ θ (u)du). Our estimator given in equation (23) provides an unbiased estimator η θ (ϕ θ ) of ∇ θ log p(y|θ), for any choice of M ≥ 0. In other words E η θ (ϕ θ ) = ∇ θ log p(y|θ).
We will see that it is most efficient to choose M = 1. To ensure the output of the SGD algorithm satisfies θ > 0, we let θ = exp(ξ) and optimize ξ. The details are given in Algorithm 4.
Algorithm 4 SGD using new unbiased estimator. 1. Initialize ξ 1 and choose a sequence {α k } ∞ k=1 and a value M ∈ N.
As above, it makes sense to first explore the toy model with analytical solution. As before, the MSE is calculated based on 50 realizations, and the cost in the plot is again proportional to the sum of the cost per forward solve at level l (tridiagonal linear system), h −1 l , multiplied by the total number of samples at level l. In Figure 3 we explore the performance of the unbiased estimator with different choices of α k = α 1 /k, α 1 ∈ {0.1, 0.025}, and different choices of the number of samples M used to construct η θ (ϕ e ξ k ) using (23) in step 2 of Algorithm 4. The two takeaways from this experiment are that (1) it is more efficient to take fewer samples M (in particular M = 1), and (2) it is more efficient to choose a larger constant α 1 = 0.1. In particular, the dynamics of the algorithm experiences a phase transition as one varies the constant α 1 . A large enough value appears to provide gradient descent type exponential convergence O(e −cost ), while a value which is too small yields Monte Carlo (MC) type O(1/cost). It is notable that the exponential convergence eventually gives way to MC type convergence, and that the point where this occurs increases proportionally to the additional constant in cost incurred with larger sample size M , so that the error curves for different values of M eventually intersect.
Natural questions are then whether there is a limit to how large one can choose α 1 and at which value precisely does the phase transition occur. These questions are partially answered by the experiments presented in Figure 4 (a), where we see that α 1 should not be chosen larger than 0.2 and the phase transition happens in between 0.025 and 0.05. Figure 4 (b) illustrates the benefits and drawbacks of using a constant α. In particular, the algorithm may converge more quickly at first, but plateaus when it reaches the induced bias.
In Figure 5 we explore various choices of P max , for M = 1 fixed and α 1 = 0.1. It is apparent that it is preferable to choose a smaller value of P max . We note however that there will be an induced bias, which will be larger for smaller P max . However, for this particular problem we do not even observe that bias over the range of MSE and cost considered.
As a last experiment with the toy example, we compare the convergence of SGD using our unbiased algorithm with P max = 0, M = 1, and α 1 = 0.1 to the analogous algorithm where an MLSMC estimator with various L (single gradient estimator MSE ∝ 2 −L ) replaces the unbiased estimator in step 2 of Algorithm 4. Similar behaviour was observed for MLSMC relative to different choices of α k as compared to the unbiased estimator. The results are shown in Figure 6. Here it is clear that over a wide range of MSE the unbiased estimator provides a significantly more efficient alternative to the MLSMC estimator.
Finally, we consider the same last experiment except with the example of Sec. 2.2. Again, over a wide range of MSE the unbiased estimator provides a significantly more efficient alternative to the MLSMC estimator. Here one can already observe the induced P max = 1 bias for the unbiased estimator around 2 −20 . Adjusting the various tuning parameters resulted in similar behaviour as was observed in the earlier experiments with the toy example. These results are not presented.

A Proofs
We begin with some technical results that have been proved in previous works (e.g. [2,6]). Recall that E N θ,m is an expectation w.r.t. the probability with finite dimensional law (17), associated to the simulation of Algorithm 1.