Efficient algorithms for Bayesian Inverse Problems with Whittle--Mat\'ern Priors

TL;DR

This paper develops efficient computational methods for Bayesian inverse problems using Whittle--Matérn Gaussian priors with noninteger exponents, enabling scalable sampling, MAP estimation, and variance approximation.

Contribution

It introduces a novel finite element and Krylov subspace-based approach for handling all admissible noninteger exponents in Whittle--Matérn priors, improving computational efficiency.

Findings

Efficient discretization of covariance operators for noninteger exponents.

Scalable algorithms for sampling and MAP estimation.

Successful application to tomography and heat equation inverse problems.

Abstract

This paper tackles efficient methods for Bayesian inverse problems with priors based on Whittle--Mat\'ern Gaussian random fields. The Whittle--Mat\'ern prior is characterized by a mean function and a covariance operator that is taken as a negative power of an elliptic differential operator. This approach is flexible in that it can incorporate a wide range of prior information including non-stationary effects, but it is currently computationally advantageous only for integer values of the exponent. In this paper, we derive an efficient method for handling all admissible noninteger values of the exponent. The method first discretizes the covariance operator using finite elements and quadrature, and uses preconditioned Krylov subspace solvers for shifted linear systems to efficiently apply the resulting covariance matrix to a vector. This approach can be used for generating samples from the distribution in two different ways: by solving a stochastic partial differential equation, and by using a truncated Karhunen-Lo\`eve expansion. We show how to incorporate this prior representation into the infinite-dimensional Bayesian formulation, and show how to efficiently compute the maximum a posteriori estimate, and approximate the posterior variance. Although the focus of this paper is on Bayesian inverse problems, the techniques developed here are applicable to solving systems with fractional Laplacians and Gaussian random fields. Numerical experiments demonstrate the performance and scalability of the solvers and their applicability to model and real-data inverse problems in tomography and a time-dependent heat equation.