Efficient Hessian computation in inverse problems with application to uncertainty quantification

Abstract

This thesis considers the efficient Hessian computation in inverse problems with specific application to the elastography inverse problem. Inverse problems use measurements of observable parameters to infer information about model parameters, and tend to be ill-posed. They are typically formulated and solved as regularized constrained optimization problems, whose solutions best fit the measured data. Approaching the same inverse problem from a probabilistic Bayesian perspective produces the same optimal point called the maximum a posterior (MAP) estimate of the parameter distribution, but also produces a posterior probability distribution of the parameter estimate, from which a measure of the solution's uncertainty may be obtained. This probability distribution is a very high dimensional function with which it can be difficult to work. For example, in a modest application with N = 104 optimization variables, representing this function with just three values in each direction requires 3^10000 U+2248 10^5000 variables, which far exceeds the number of atoms in the universe. The uncertainty of the MAP estimate describes the shape of the probability distribution and to leading order may be parameterized by the covariance. Directly calculating the Hessian and hence the covariance, requires O(N) solutions of the constraint equations. Given the size of the problems of interest (N = O(10^4 - 10^6)), this is impractical. Instead, an accurate approximation of the Hessian can be assembled using a Krylov basis. The ill-posed nature of inverse problems suggests that its Hessian has low rank and therefore can be approximated with relatively few Krylov vectors. This thesis proposes a method to calculate this Krylov basis in the process of determining the MAP estimate of the parameter distribution. Using the Krylov space based conjugate gradient (CG) method, the MAP estimate is computed. Minor modifications to the algorithm permit storage of the Krylov approximation of the Hessian. As the accuracy of the Hessian approximation is directly related to the Krylov basis, long term orthogonality amongst the basis vectors is maintained via full reorthogonalization. Upon reaching the MAP estimate, the method produces a low rank approximation of the Hessian that can be used to compute the covariance.