A Bayesian Conjugate Gradient Method (with Discussion)

Abstract

A fundamental task in numerical computation is the solution of large linear systems. The conjugate gradient method is an iterative method which offers rapid convergence to the solution, particularly when an effective preconditioner is employed. However, for more challenging systems a substantial error can be present even after many iterations have been performed. The estimates obtained in this case are of little value unless further information can be provided about, for example, the magnitude of the error. In this paper we propose a novel statistical model for this error, set in a Bayesian framework. Our approach is a strict generalisation of the conjugate gradient method, which is recovered as the posterior mean for a particular choice of prior. The estimates obtained are analysed with Krylov subspace methods and a contraction result for the posterior is presented. The method is then analysed in a simulation study as well as being applied to a challenging problem in medical imaging.

Probabilistic linear solvers: a unifying view

“…This corollary suggests a natural way to select a prior covariance still linked to the linear system, though this choice is still not computationally convenient. Furthermore, in the case that A is symmetric positive-definite, this recovers the prior which replicates CG described in Cockayne et al [2018]. Note that each of H and P can be stated explicitly as H = (A A) 1 2 and P = A(A A) − 1 2 .…”

“…Conjugate gradients has been studied from a probabilistic point of view before by Hennig [2015] and Cockayne et al [2018]. This section generalizes the results of Hennig [2015] and leverages Proposition 6 for new insights on BayesCG.…”

“…where A ∈ R d×d is an invertible matrix, b ∈ R d is a given vector and x * ∈ R d is an unknown to be determined. Recent work [Hennig, 2015, Cockayne et al, 2018 has constructed iterative solvers for this problem which output probability measures, constructed to quantify uncertainty due to terminating the algorithm before the solution has been identified completely. On the surface the approaches in these two works appear different:…”

“…In the matrix-based inference (MBI) approach of Hennig [2015], a posterior is constructed on the matrix A −1 , while in the solution-based inference (SBI) method of Cockayne et al [2018] a posterior is constructed on the solution vector x * . These algorithms are instances of probabilistic numerical methods (PNM) in the sense of and Cockayne et al [2017].…”

“…The BayesCG algorithm proposed by Cockayne et al [2018] encompasses conjugate gradients as a special case. BayesCG uses left-multiplied observations and was derived in the solution-based perspective.…”