Comments on the article "A Bayesian conjugate gradient method"

TL;DR

This paper comments on a Bayesian conjugate gradient method that provides probabilistic estimates and uncertainty quantification for solving linear systems, suggesting future research directions.

Contribution

It offers critical insights, questions, and potential research directions for the Bayesian conjugate gradient approach introduced in the original article.

Findings

Provides a critique and analysis of the Bayesian conjugate gradient method.

Raises questions about the method's assumptions and applicability.

Suggests avenues for further theoretical and practical research.

Abstract

The recent article "A Bayesian conjugate gradient method" by Cockayne, Oates, Ipsen, and Girolami proposes an approximately Bayesian iterative procedure for the solution of a system of linear equations, based on the conjugate gradient method, that gives a sequence of Gaussian/normal estimates for the exact solution. The purpose of the probabilistic enrichment is that the covariance structure is intended to provide a posterior measure of uncertainty or confidence in the solution mean. This note gives some comments on the article, poses some questions, and suggests directions for further research.

Equations

$$A^{†}\mathbf{b}=\mathrm{argmin}\left\{\Vert \mathbf{x}{\Vert }_{{\mathbb{R}}^d}|\mathbf{x}\in \mathrm{argmin}\left\{\Vert A{\mathbf{x}}^{\prime }-\mathbf{b}{\Vert }_{{\mathbb{R}}^c}|{\mathbf{x}}^{\prime }\in {\mathbb{R}}^d\right\}\right\}$$

$${\mu }_m:=\mathcal{N}(P_m{\mathbf{x}}^{\ast },0)\otimes \mathcal{N}(P_m^{\perp }{\mathbf{x}}_0,P_m^{\perp }{\Sigma }_0{P}_m^{\perp })$$

Full text

Abstract:

The recent article “A Bayesian conjugate gradient method” by Cockayne, Oates, Ipsen, and Girolami proposes an approximately Bayesian iterative procedure for the solution of a system of linear equations, based on the conjugate gradient method, that gives a sequence of Gaussian/normal estimates for the exact solution. The purpose of the probabilistic enrichment is that the covariance structure is intended to provide a posterior measure of uncertainty or confidence in the solution mean. This note gives some comments on the article, poses some questions, and suggests directions for further research.

Keywords:

probabilistic numerics, linear systems, Krylov subspaces, uncertainty quantification.

2010 Mathematics Subject Classification:

62C10, 62F15, 65F10.

1 Overview

In “A Bayesian conjugate gradient method”, Cockayne, Oates, Ipsen, and Girolami add to the recent body of work that provides probabilistic/inferential perspectives on deterministic numerical tasks and algorithms. In the present work, the authors consider a conjugate gradient (CG) method for the solution of a finite-dimensional linear system $A\boldsymbol{x}^{\ast}=\boldsymbol{b}$ for $\boldsymbol{x}^{\ast}\in\mathbb{R}^{d}$, given an assuredly invertible symmetric matrix $A\in\mathbb{R}^{d\times d}$ and $\boldsymbol{b}\in\mathbb{R}^{d}$, with $d\in\mathbb{N}$.

The authors derive and test an algorithm, BayesCG, that returns a sequence of normal distributions $\mathcal{N}(\boldsymbol{x}_{m},\Sigma_{m})$ for $m=0,\dots,d$, starting from a prior distribution $\mathcal{N}(\boldsymbol{x}_{0},\Sigma_{0})$. This sequence of normal distributions is defined using a recursive relationship similar to that defining the classical CG method, and indeed the BayesCG mean $\boldsymbol{x}_{m}$ coincides with the output of CG upon choosing $\Sigma_{0}\coloneqq A^{-1}$ — this choice is closely related to what the authors call the “natural prior covariance”, $\Sigma_{0}\coloneqq(A^{\top}A)^{-1}$ . The distribution $\mathcal{N}(\boldsymbol{x}_{m},\Sigma_{m})$ is intended to be an expression of posterior belief about the true solution $\boldsymbol{x}^{\ast}$ to the linear system under the prior belief $\mathcal{N}(\boldsymbol{x}_{0},\Sigma_{0})$ given the first $m$ BayesCG search directions $\boldsymbol{s}_{1},\dots,\boldsymbol{s}_{m}$. Like CG, BayesCG terminates in at most $d$ steps, at which point its mean $\boldsymbol{x}_{d}$ is the exact solution $\boldsymbol{x}^{\ast}$ and it expresses complete confidence in this belief by having $\Sigma_{d}=0$. The convergence and frequentist coverage properties of the algorithm are investigated in a series of numerical experiments.

The field of probabilistic perspectives on numerical tasks has been enjoying a resurgence of interest in the last few years; see e.g. Oates and Sullivan (2019) for a recent survey of both historical and newer work. The authors’ contribution is a welcome addition to the canon, showing as it does how classical methods (in this case CG; cf. the treatment of Runge–Kutta methods for ODEs by Schober et al. (2014)) can be seen as point estimators of particular instances of inferential procedures. it is particularly encouraging to see contributions coming from authors with both statistical and numerical-analytical expertise, and the possibilities for generalisation and further work are most interesting.

2 Questions and directions for generalisation

The article raises a number of natural questions and directions for generalisation and further investigation, which Cockayne et al. might use their rejoinder to address.

Symmetry and generalisations.

It is interesting that, for the most part, BayesCG does not require that $A$ be symmetric, and works instead with $A\Sigma_{0}A^{\top}$. This prompts a question for the authors to consider in their rejoinder: How does BayesCG behave in the case that $A\in\mathbb{R}^{d\times d}$ is square but not invertible? Could one even show that the BayesCG posterior concentrates, within at most $d$ iterations, to a distribution centred at the minimum-norm solution (in some norm on $\mathbb{R}^{d}$) of the linear system? One motivation for this question is that, if such a result could be established, then BayesCG could likely be applied to rectangular $A\in\mathbb{R}^{c\times d}$ and produce a sequence of normally-distributed approximate solutions to the minimum-norm least-squares problem, i.e. a probabilistically-motivated theory of Moore–Penrose pseudo-inverses. (Recall that the Moore–Penrose pseudo-inverse $A^{\dagger}\in\mathbb{R}^{d\times c}$ can be characterised as the solution operator for the minimum-norm least squares problem, i.e.,

[TABLE]

for the usual norms on $\mathbb{R}^{c}$ and $\mathbb{R}^{d}$.)

The authors could also comment on whether they expect BayesCG to generalise easily to infinite-dimensional Hilbert spaces, analogously to CG (Fortuna, 1977, 1979; Málek and Strakoš, 2015), since experience has shown that analysis of infinite-dimensional statistical algorithms can yield powerful dimension-independent algorithms for the finite-dimensional setting (Cotter et al., 2013; Chen et al., 2018). A more esoteric direction for generalisation would be to consider fields other than $\mathbb{R}$. The generalisation of BayesCG to $\mathbb{C}$ would appear to be straightforward via a complex normal prior, but how about fields of finite characteristic?

Conditioning relations.

Could the authors use their rejoinder to provide some additional clarity about the relationship of BayesCG to exact conditioning of the prior normal distribution $\mathcal{N}(\boldsymbol{x}_{0},\Sigma_{0})$, and in particular whether BayesCG is exactly replicating the ideal conditioning step, approximating it, or doing something else entirely?

To make this question more precise, fix a weighted inner product on $\mathbb{R}^{d}$, e.g. the Euclidean, $A$-weighted, $\Sigma_{0}$-weighted, or $A\Sigma_{0}A^{\top}$-weighted inner product. With respect to this inner product, let $P_{m}$ be the (self-adjoint) orthogonal projection onto the Krylov space $\mathcal{K}_{m}$ and $P_{m}^{\perp}\coloneqq I-P_{m}$ the (self-adjoint) orthogonal projection onto its orthogonal complement. The product measure

[TABLE]

on $\mathcal{K}_{m}\oplus\mathcal{K}_{m}^{\perp}$ is also a normal distribution on $\mathbb{R}^{d}$ that expresses complete confidence about the true solution to the linear system in the directions of the Krylov subspace and reverts to the prior in the complementary directions; $\mu_{0}$ is the prior, and $\mu_{d}$ is a Dirac on the truth. The obvious question is, for some choice of weighted inner product, does BayesCG basically output $\mu_{m}$, but in a clever way? Or is the output of BayesCG something else?

Uncertainty quantification: cost, quality, and terms of reference.

Cockayne et al. point out that the computational cost of BayesCG is a small multiple (a factor of three) of the cost of CG, and this would indeed be a moderate price to pay for high-quality uncertainty quantification. However, based on the results presented in the paper, the UQ appears to be quite poorly calibrated. One could certainly try to overcome this shortcoming by improving the UQ. However, an alternative approach would be to choose a prior covariance structure $\Sigma_{0}$ that aggressively sparsifies and hence accelerates the linear-algebraic operations that BayesCG performs (particularly lines 7, 9, and 11 of Algorithm 1), while preserving the relatively poor UQ. Does this appear practical, in their view?

In fact, the whole question of “the UQ being well calibrated” is somewhat ill defined. Some might argue that, since $\boldsymbol{x}^{\ast}$ is deterministic, there is no scope for uncertainty in this setting. However, It certainly makes sense to ask — as the authors do in Section 6.1 – whether, when the problem setup is randomised, the frequentist coverage of BayesCG lines up with that implied by the randomisation of the problem. The statistic $Z$ that Cockayne et al. introduce is intersting, and already captures several scenarios in which the UQ is well calibrated and several in which it is not; I strongly encourage Cockayne et al. to address in their rejoinder the Bayesian accuracy of BayesCG, e.g., to exhaustively sample the true posterior on $\boldsymbol{x}^{\ast}$ given $\boldsymbol{y}_{1},\dots,\boldsymbol{y}_{m}$ and see whether this empirical distribution is well approximated by the BayesCG distribution $\mathcal{N}(\boldsymbol{x}_{m},\Sigma_{m})$.

As a related minor question, can BayesCG be seen as an (approximate Gaussian) filtering scheme for the given prior $\mathcal{N}(\boldsymbol{x}_{0},\Sigma_{0})$ and the filtration associated to the data stream $\boldsymbol{y}_{1},\boldsymbol{y}_{2},\dots$?

Precision formulation.

As formulated, BayesCG expresses a Gaussian belief about the solution of the linear system in terms of mean and covariance; Gaussian measures can also be expressed in terms of their mean and precision. Do the authors have a sense of whether the BayesCG algorithm be formulated as a sequence of precision updates, or does the singularity of the covariance in the Krylov directions essentially forbid this?

One motivation for seeking a precision formulation would be to render the “natural prior” $\Sigma_{0}\coloneqq(A^{\top}A)^{-1}$ of Section 4.1 more tractable, since this prior has an easily-accessible precision while accessing its covariance involves solving the original linear problem.

As the authors note, their “natural prior” is closely related to one introduced by Owhadi (2015) and applied in Cockayne et al. (2017). It seems that working with images of Gaussian white noise, such as this “natural prior”, is presently producing considerable analytical and computational strides forward (Owhadi, 2017; Chen et al., 2018), and so this seems to be a topic worth further attention in the statistical community as a whole.

3 Minor comments

Rank and trace estimates.

Proposition 3, which states that $\operatorname{tr}(\Sigma_{m}\Sigma_{0}^{-1})=d-m$, seems to miss the point slightly. It would be good to have a companion results to the effect that $\operatorname{rank}\Sigma_{m}=d-m$, and a more quantitative result for the trace such as $\operatorname{tr}\Sigma_{m}\leq C_{\Sigma_{0}}(d-m)$ for some constant $C_{\Sigma_{0}}\geq 0$ depending only on $\Sigma_{0}$.

Posterior and Krylov spaces.

It seems natural to ask whether the posterior nullspace $\ker\Sigma_{m}$ coincides with the Krylov space $\mathcal{K}_{m}$. Put another way, is the posterior column space the same as the orthogonal complement of the Krylov space, in the Euclidean or $A$-weighted inner product?

Square roots.

Is the square root $M^{1/2}$ introduced just before Proposition 5 required to be unique? Does it even matter whether a (unique) symmetric positive semidefinite square root or a Cholesky factor is chosen? This is related to the above discussion of symmetry and generalisations.

Interpretation of termination criteria.

In Section 5, the authors refer to “probabilistic termination criteria”. As a matter of semantics, the termination criterion that they propose is in fact deterministic, albeit based on a probabilistic interpretation of the algorithmic quantities.

Acknowledgements

TJS is supported by the Freie Universität Berlin within the Excellence Strategy of the German Research Foundation (DFG), including the ECMath/MATH+ transition project CH-15 and project TrU-2 of the Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID 390685689).