Probabilistic Linear Solvers: A Unifying View

Several recent works have developed a new, probabilistic interpretation for numerical algorithms solving linear systems in which the solution is inferred in a Bayesian framework, either directly or by inferring the unknown action of the matrix inverse. These approaches have typically focused on replicating the behavior of the conjugate gradient method as a prototypical iterative method. In this work surprisingly general conditions for equivalence of these disparate methods are presented. We also describe connections between probabilistic linear solvers and projection methods for linear systems, providing a probabilistic interpretation of a far more general class of iterative methods. In particular, this provides such an interpretation of the generalised minimum residual method. A probabilistic view of preconditioning is also introduced. These developments unify the literature on probabilistic linear solvers, and provide foundational connections to the literature on iterative solvers for linear systems.


Consider the linear system
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]. PNM are numerical methods which output posterior distributions that quantify uncertainty due to discretisation error. An interesting property of PNM is that they often result in a posterior distributions whose mean element coincides with the solution given by a classical numerical method for the problem at hand. The relationship between PNM arXiv:1810.03398v2 [stat.CO] 17 Oct 2018 and classical solvers has been explored for integration [e.g. Karvonen and Sarkka, 2017], ODE-solvers [Schober et al., 2014, Kersting et al., 2018 and PDE solvers [Cockayne et al., 2016] in some generality. For linear solvers, attention has thus far been restricted to the conjugate gradient (CG) method. Since CG is but a single member of a larger class of iterative solvers, and applicable only if the matrix A is symmetric and positivedefinite, extending the probabilistic interpretation is still an interesting endeavour. Probabilistic interpretations provide an alternative perspective on numerical algorithms, and can also provide extensions such as the ability to exploit noisy or corrupted observations. The probabilistic view has also been used to the develop new numerical methods [Xi et al., 2018], and Bayesian PNM can be incorporated rigorously into pipelines of computation [Cockayne et al., 2017].
Preconditioning-mapping Eq. (1) to a better conditioned system whith the same solution-is key to the fast convergence of iterative linear solvers, particularly those based upon Krylov methods [Liesen and Strakos, 2012]. The design of preconditioners has been referred to as "a combination of art and science" [Saad, 2003, p. 283]. In this work we also provide a new, probabilistic interpretation of preconditioning as a form of prior information.

Contribution
This text contributes three primary insights: 1. It is shown that, for particular choices of the generative model, matrix-based inference (MBI) and solution-based inference (SBI) can be equivalent (Section 2). 2. A general probabilistic interpretation of projection methods [Saad, 2003] is described (Section 3.1), leading to a probabilistic interpretation of the generalised minimum residual method (GMRES; Saad and Schultz [1986], Section 6). The connection to CG is expanded and made more concise in Section 5. 3. A probabilistic interpretation of preconditioning is presented in Section 4.
Most of the proofs are presented inline; lengthier proofs are deferred to Appendix B. While an important consideration, the predominantly theoretical contributions of this paper will not consider the impact of finite numerical precision.

Notation
For a symmetric positive-definite matrix M ∈ R d×d and two vectors v, w ∈ R d , we write v, w M = v M w for the inner product induced by M , and v 2 M = v, v M for the corresponding norm.
A set of vectors s 1 , . . . , s m is called M -orthogonal or M -conjugate if s i , s j M = 0 for i = j, and Morthonormal if, in addition, s i M = 1 for 1 ≤ i ≤ m.
For a square matrix A = a 1 . . . a d ∈ R d×d , the vectorisation operator vec : R d×d → R d 2 stacks the rows 1 of A into one long vector: The Krylov space of order m generated by the matrix A ∈ R d×d and the vector b ∈ R d is We will slightly abuse notation to describe shifted and scaled subspaces of R d : Let S be an m-dimensional linear subspace of R d with basis {s 1 , . . . , s m }. Then for a vector v ∈ R d and a matrix M ∈ R d×d , let v + M S = span(v + M s 1 , . . . , v + M s m ).

Probabilistic Linear Solvers
Several probabilistic framework describing the solution of Eq. (1) have been constructed in recent years. They primarily differ in the subject of inference: SBI approaches such as Cockayne et al. [2018], of which BayesCG is an example, place a prior distribution on the solution x * of Eq. (1). Conversely, the MBI approach of Hennig [2015] and Bartels and Hennig [2016] places a prior on A −1 , treating the action of the inverse operator as an unknown to be inferred 2 . This section reviews each approach and adds some new insights. In particular, SBI can be viewed as strict special case of MBI (Section 2.4). Throughout this section, we will assume that the search 1 Stacking the columns is equivalently possible and common. It is associated with a permutation in the definition of the Kronecker product, but the resulting inferences are equivalent.
2 Hennig [2015] also discusses inference over A. This model class will not be discussed further in the present work. It has the disadvantage that the associated marginal on x * is nonanalytic, but more easily lends itself to situations with noisy or otherwise perturbed matrix-vector products as observations. directions S m in S m Ax = S m b are given a-priori ; Section 5 examines algorithms which iteratively generate search directions adapted to the problem at hand.

Background on Gaussian conditioning
The propositions in this section follow from the following two classic properties of Gaussian distributions.
Lemma 2 Let x ∈ R d be distributed as in Lemma 1, and let observations y ∈ R n be generated from the conditional density with M ∈ R n×d , z ∈ R n , and Λ ∈ R n×n again positivesemidefinite. Then the associated conditional distribution on x after observing y is again Gaussian, with This formula also applies if Λ = 0, i.e. observations are made without noise, with the caveat that if M ΣM is singular, the inverse should be interpreted as a pseudoinverse.

Solution-Based Inference
To phrase the solution of Eq. (1) as a form of probabilistic inference, Cockayne et al. [2018] consider a Gaussian prior over the solution x * , and condition on observations provided by a set of search directions s 1 , . . . , s m , m < d. Let S m ∈ R d×m be given by S m = [s 1 , . . . , s m ], and let information be given by y m := S m Ax * = S m b. Since the information is clearly a linear projection of x * , the posterior distribution is a Gaussian distribution on x * : Lemma 3 (Cockayne et al. [2018]) Assume that the columns of S m are linearly independent. Consider the prior The posterior from SBI is then given by The following proposition establishes an optimality property of the posterior mean x m . This is a relatively wellknown property of Gaussian inference, but has not appeared before in the literature on these methods and will prove useful in subsequent sections.
Proposition 4 If S m = range(S m ), then the posterior mean in Lemma 3 satisfies the optimality property Proof With the abbreviations X = Σ 0 A S m and y = x * − x 0 the mean in Lemma 3 can be written as is the solution of the weighted least squares problem [Golub and Van Loan, 2013, Section 6.1] This is equivalent to the desired statement.

Matrix-Based Inference
In contrast to SBI, the MBI approach of Hennig [2015] treats the matrix inverse A −1 as the unknown in the inference procedure. As in the previous section, search directions S m yield matrix-vector products Y m ∈ R d×m . In Hennig [2015] these arise from right-multiplying 3 A with S m , i.e. Y m = AS m . Note that Thus S m is a linear transformation of A −1 and Lemma 2 can again be applied: Lemma 5 (Lemma 2.1 in Hennig [2015] 4 ) Consider the prior Then the posterior given the observations For linear solvers, the object of interest is − − → A −1 , and again using Lemma 1, we see that the associated marginal is also Gaussian, and given by In the Kronecker product specification for the prior covariance on A −1 , the first matrix, here Σ 0 , describes the dependence between the columns of A −1 . The second matrix, W 0 , captures the dependency between the rows of A −1 . Note that in Theorem 5, the posterior covariance has the form Σ 0 ⊗ W m . When compared to the prior covariance, Σ 0 ⊗ W 0 , it is clear that the observations have conveyed no new information to the first term of the Kronecker product covariance.

Equivalence of MBI and SBI
In practise Hennig [2015] notes that inference on A −1 should be performed only implicitly, avoiding the d 2 storage cost and the mathematical complexity of the operations involved in Lemma 5. This raises the question of when MBI is equivalent to SBI. Although, based on Lemma 1, one might suspect SBI and MBI to be equivalent, in fact the posterior from Lemma 5 is structurally different to the posterior in Lemma 3: After projecting into solution space, the posterior covariance in Lemma 5 is a scalar multiple of the matrix Σ 0 , which is not the case in general in Lemma 3. However, the implied posterior over the solution vector can be made to coincide with the posterior from SBI if one considers observations in MBI as That is, as left-multiplications of A. We will refer to the observation model of Eq.
Proposition 6 Consider a Gaussian MBI prior conditioned on the left-multiplied information of Eq. (5).
The associated marginal on x is identical to the posterior on x arising in Lemma 3 from p(x) = N (x; x 0 , Σ 0 ) under the conditions Proof See Appendix B.
The first of the two conditions requires that the prior mean on the matrix inverse be consistent with the prior mean on the solution, which is natural. The second condition demands that, after projection into solution space, the relationship between the rows of A −1 modelled by W 0 does not inflate the covariance Σ 0 . Note that this condition is trivial to enforce for an arbitrary covariancē

Remarks
The result in Proposition 6 shows that any result proven for SBI applies immediately to MBI with left-multiplied observations. Though MBI has more model parameters than SBI, there are situations in which this point of view is more appropriate. Unlike in SBI, the information obtained in MBI need not be specific to a particular solution vector x * and thus can be propagated and recycled over several linear problems, similar to the notion of subspace recycling [Soodhalter et al., 2014]. Secondly, MBI is able to utilise both left-and rightmultiplied information, while SBI is restricted to leftmultiplied information. This additional generality may prove useful in some applications.

Projection Methods as Inference
This section discusses a connection between probabilistic numerical methods for linear systems and the classic framework of projection methods for the iterative solution of linear problems. Section 3.1 reviews this established class of solvers, while Section 3.2 presents the novel results.

Background
Many iterative methods for linear systems, including CG and GMRES, belong to the class of projection methods [Saad, 2003, p. 130f.]. Saad describes a projection method as an iterative scheme in which, at each iteration, a solution vector x m is constructed by projecting More formally, each iteration of a projection method is defined by two matrices X m , U m ∈ R d×m , and by a starting point x 0 . The matrices X m and U m each encode the solution and constraint spaces as X m = range(X m ) and U m = range(U m ). The projection method then constructs x m as From this perspective CG and GMRES perform only a single step with the number of iterations m fixed and determined in advance. For CG the spaces are

Probabilistic Perspectives
In this section we first show, in Proposition 7, that the conditional mean from SBI after m steps corresponds to some projection method. Then, in Proposition 8 we prove the converse: that each projection method is also the posterior mean of a probabilistic method, for some prior covariance and choice of information.
Proposition 7 Let the columns of S m be linearly independent. Consider SBI under the prior and with observations y m = S m b. Then the posterior mean x m in Lemma 3 is identical to the iterate from a projection method defined by the matrices U m = S m and X m = Σ 0 A S m , and the starting vector x 0 .
Proof Substituting U m = S m and X m = Σ 0 A S m into Lemma 3 gives Eq. (7), as required.
The converse to this also holds: Proposition 8 Consider a projection method defined by the matrices X m , U m ∈ R d×m , each with linearly independent columns, and the starting vector x 0 ∈ R d . Then the iterate x m in Eq. (7) is identical to the SBI posterior mean in Lemma 3 under the prior when search directions S m = U m are used.
Proof Abbreviate Z = X m A U m and write the projection method iterate from Eq. (7) as Multiply the middle matrix by the identity, and insert this into the expression for x 0 , Setting U m = S m gives the mean in Lemma 3.
Including a basis of the solution space in the prior may seem problematic. A direct way to enforce the posterior occupying the solution space is by placing a prior on the coefficients α in x = x 0 + X m α. Under a unit Gaussian prior α ∼ N (0, I), the implied prior on x naturally has the form of Eq. (8). However, this prior is nevertheless unsatisfying both since it requires the solution space to be specified a-priori, precluding adaptivity in the algorithm, and, perhaps more worryingly, because the posterior uncertainty over the solution is a matrix of zeros even though the solution is not fully identified. Again taking Z = X m A U m : [ Hennig, 2015] and [Bartels and Hennig, 2016] each proposed to address this issue by adding additional uncertainty in the null space of X m . This empirical uncertainty calibration step has not yet been analysed in detail. Such analysis is left for future work. Nevertheless, the proposition provides a probabilistic view for arbitrary projection methods and does not require knowledge of A −1 , unlike some of the results presented in [Hennig, 2015, Cockayne et al., 2017 and in the following propositions.
This prior is not unique. The next proposition establishes more restrictive conditions under which a projection method may have a probabilistic interpretation and still result in a nonzero posterior uncertainty.
Proposition 9 Consider a projection method defined by X m , U m ∈ R d×m and the starting vector x 0 . Further suppose that U m = RX m for some invertible R ∈ R d×d , and that A R is symmetric positive-definite. Then under the prior and the search directions S m = U m = RX m , the iterate in the projection method is identical to the posterior mean in Lemma 3.
Proof First substitute X m = R −1 U m into Eq. (7) to obtain A corollary which provides further insight arises when one considers the polar decomposition of A. Recall that an invertible matrix A has a unique polar decomposition A = P H, where P ∈ R d×d is orthogonal and H ∈ R d×d is symmetric positive-definite.
Corollary 10 Consider a projection method defined by X m , U m ∈ R d×m and the starting vector x 0 , and suppose that U m = P X m , where P arises from the polar decomposition A = P H. Then under the prior and the search directions S m = U m = P X m , the iterate in the projection method is identical to the posterior mean in Lemma 3.
Proof This follows from Proposition 9. Setting R = P aligns the search directions in Corollary 10 with those in Proposition 9. Since P is orthogonal, P −1 = P , and since H is symmetric positive-definite, A P = P A = H by definition of the polar decomposition, which gives the prior covariance required for Proposition 9. This is an intuitive analogue of similar results in Hennig [2015] and Cockayne et al. [2017] which show that CG is recovered under certain conditions involving a prior Σ 0 = A −1 . When A is not symmetric and positive definite it cannot be used as a prior covariance. 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 . Thus in the case of symmetric positive-definite A we have that H = A and P = I, so that the prior covariance Σ 0 = A −1 arises naturally from this interpretation.

Preconditioning
This section discusses probabilistic views on preconditioning. Preconditioning is a widely-used technique accelerating the convergence of iterative methods [Saad, 2003, Sections 9 and 10]. A preconditioner P is a nonsingular matrix satisfying two requirements: 1. Linear systems P z = c can be solved at low computational cost (i.e. "analytically") 2. P is "close" to A in some sense.
In this sense, solving systems based upon a preconditioner can be viewed as approximately inverting A, and indeed many preconditioners are constructed based upon this intuition. One distinguishes between right preconditioners P r and left preconditioners P l , depending on whether they act on A from the left or the right. Two-sided preconditioning with nonsingular matrices P l and P r transforms implicitly Eq. (1) into a new linear problem The preconditioned system can then be solved using arbitrary projection methods as described in Section 3.1, from the starting point z 0 defined by x 0 = P r z 0 . The probabilistic view can be used to create a nuanced description of preconditioning as a form of prior information. In the SBI framework, Proposition 11 below shows that solving a right-preconditioned system is equivalent to modifying the prior, while in Proposition 12 shows that left-preconditioning is equivalent to making a different choice of observations.
Proposition 11 (Right preconditioning) Consider the right-preconditioned system SBI on Eq. (10) under the prior is equivalent to solving Eq. (1) under the prior x ∼ N (x; P r z 0 , P r Σ 0 P r ).
Proof Let p(x) = N (x; x 0 , Σ r ). Lemma 3 implies that after observing information from search directions S m , the posterior mean equals where r 0 = b − Ax 0 . Setting x 0 = P r z 0 and letting Σ r = P r Σ 0 P r gives shows that this is equivalent to Thus z m is the posterior mean of the system Bz * = b with prior Eq. (11) after observing search directions S m .
Proposition 12 (Left preconditioning) Consider the left-preconditioned system And the SBI prior Then the posterior from SBI on Eq. (12) under search directions S m is equivalent to the posterior from SBI applied to the system Eq.
(1) under search directions P l S m .
Proof Lemma 3 implies that after observing search directions T m , the posterior mean over the solution of Eq.
(1) equals where B := P l A andr 0 = P l b − P l Ax 0 . Thus, x m is the posterior mean of the system Bx * = P l b after observing search directions S m .
If a probabilistic linear solver has a posterior mean which coincides with a projection method (as discussed in Section 3.1), the Propositions 11 and 12 show how to obtain a probabilistic interpretation of the preconditioned version of that algorithm. Furthermore, the equivalence demonstrated in Section 2.4 shows that the reasoning from Propositions 11 and 12 carries over to MBI based on left-multiplied observations: right-preconditioning corresponds to a change in prior belief, while left-preconditioning corresponds to a change in observations. We do not claim that this probabilistic interpretation of preconditioning is unique. For example, when using MBI with right-multiplied observations, the same line of reasoning can be used to show the converse: right-preconditioning corresponds to a change in the observations and left-preconditioning to a change in the prior.

Conjugate Gradients
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. For this Section (but not thereafter) assume that A is a symmetric and positive definite matrix.

Left-multiplied view
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.
The posterior in Lemma 3 does not immediately result in a practical algorithm as it involves the solution of a linear system based on the matrix S m AΣ 0 A S m ∈ R m×m , which requires O(m 3 ) arithmetic operations. BayesCG avoids this cost by constructing search directions that are AΣ 0 A -orthonormal, as shown below, see [Cockayne et al., 2018, Proposition 7].
With these search directions constructed, BayesCG becomes an iterative method: Proposition 14 (Proposition 6 of Cockayne et al. [2018]) Using the search directions from Proposition 13, the posterior from Lemma 3 reduces to: In Proposition 4 of Cockayne et al. [2018] it was shown that the BayesCG posterior mean corresponds to the CG solution estimate when the prior covariance is taken to be Σ 0 = A −1 , though this is not a practical choice of prior covariance as it requires access to the unavailable A −1 . Furthermore, in Proposition 9 it was shown that when using the search directions from Proposition 13, the posterior mean from BCG has the following optimality property: Note that this is now a trivial special case of Proposition 4. The following proposition leverages these results along with Proposition 6 to show that there exists an MBI method which, under a particular choice of prior and with a particular methodology for the generation of search directions, is consistent with CG.
Proposition 15 Consider the MBI prior where W 0 ∈ R d is symmetric positive-definite and so that b W 0 b = 1. Suppose left-multiplied information is used, and that the search directions are generated sequentially according to: Then it holds that the implied posterior mean on solution space, given by A −1 m b, corresponds to the CG solution estimate after m iterations, with starting point Proof First note that, by Proposition 6, since left-multiplied observations are used and since b W 0 b = 1, the implied posterior distribution on solution space from MBI is identical to the posterior distribution from SBI under the prior It thus remains to show that the sequence of search directions generated is identical to those in Proposition 13 for this prior. Fors 1 : as required. Fors j : where the second line uses that A −1 j−1 b = x j−1 . Thus, the search directions coincide with those in Proposition 13. It therefore holds that the implied posterior mean on solution space, A −1 m b, coincides with the solution estimate produced by CG.

Right-multiplied view
Interpretations of CG (and general projection methods) that use right-multiplied observations seems to require more care than those based on left-multiplied observations. Nevertheless, Hennig [2015] provided an interpretation for CG in this framework, essentially showing 5 that Algorithm 1 reproduces both the search directions and solution estimates from CG under the prior where α ∈ R \ {0}, β ∈ R + and ⊗ denotes the symmetric Kronecker product (see Section A.1). The posterior under such a prior is described in Lemma 2.2 of Hennig [2015] (see Lemma 21), though we note that the sense in which the solution estimate x m output by this algorithm is related to the posterior over A −1 differs from that in the previous section, in the sense that as the CG estimate is corrected by the step size computed in line 6. Fixing this rank-1 discrepancy would complicate the exposition of Algorithm 1 and yield a more cumbersome algorithm). The following proposition generalizes this result.
For all choices α ∈ R \ {0} and β, γ ∈ R +,0 with β + γ > 0, Algorithm 1 is equivalent to CG, in the sense that it produces the exact same sequence of estimates x i and scaled search directions s i .
Proof The proof is extensive and has been moved to Appendix B.
Algorithm 1 The algorithm referred to by Proposition 16, which reproduces the search directions and solution estimates from CG.
12 end for 13 return x m 5 Algorithm 1 is not included in this form in the op.cit.
Note that, unlike previous propositions, Proposition 16 proposes a prior that does not involve A −1 for the case when γ = 0.

GMRES
The Generalised Minimal Residual Method [Saad, 2003, Section 6.5] applies to general nonsingular matrices A. At iteration m, GMRES minimises the residual over the affine space , this corresponds to minimizing the error in the A A norm. We present a brief development of GMRES, starting with Arnoldi's method (Section 6.1) and the GMRES algorithm (Section 6.2), before presenting our Bayesian interpretation (Section 6.3).

Arnoldi's Method
GMRES uses Arnoldi's method [Saad, 2003, Section 6.3] to construct orthonormal bases for Krylov spaces of general, nonsingular matrices A. Starting with q 1 = r 0 / r 0 2 , Arnoldi's method recursively computes the orthonormal basis Q m = q 1 . . . q m ∈ R d×m for K m (A, r 0 ). The basis vectors satisfy the relations and

GMRES
GMRES computes the iterate based on the optimality condition in Eq. (13), which can equivalently be expressed as Thus confirming that GMRES is a projection method with X m = Q m and U m = AQ m . GMRES solves the least squares problem in Eq. (15). efficiently by projecting it to a lower dimensional space via Arnoldi's method. To this end, express the starting vector in the Krylov basis, r 0 = r 0 2 q 1 = r 0 2 Q m+1 e 1 , and exploit the Arnoldi recursion from Eq. (14), The computations are summarized in Algorithm 2.

Bayesian Interpretation of GMRES
We now present probabilistic linear solvers with posterior means that coincide with the solution estimate from GMRES.

Left-multiplied view
Proposition 17 Under the SBI prior and the search directions U m = AQ m , the posterior mean is identical to the GMRES iterate x m in Eq. (16).
Proof Substitute R = A and U m = AQ m into Proposition 9.
Proposition 17 is intuitive in the context of Proposition 4: Setting Σ 0 = (A A) −1 ensures that the norm being minimised coincides with that of GMRES, as does the solution space X m = AQ m . This interpretation exhibits an interesting duality with CG for which Σ 0 = A −1 . Another probabilistic interpretation follows from Proposition 8.
Corollary 18 Under the prior and with observations y m = Q m b, the posterior mean from SBI is identical to the GMRES iterate x m in Eq. (16).
Note that Proposition 17 has a posterior covariance which is not practical, as it involves A −1 . [Cockayne et al., 2017] proposed replacing A −1 in the prior covariance with a preconditioner to address this, which does yield a practically computable posterior, but this extension was not explored here. Furthermore, that approach yields poorly calibrated posterior uncertainty, as described in that work. Corollary 18 does not have this drawback, but the posterior covariance is a matrix of zeroes.

Right-multiplied view
As for CG in Section 5.2, finding interpretations of GMRES that use right-multiplied observations appears to be more difficult.
Proposition 19 Under the prior and given Y m = AQ m , the implied posterior mean on the solution space given by A −1 m b is equivalent to the GMRES solution. This correspondence breaks when x 0 = 0.
Proof Under this prior, b applied to the posterior mean is which is the GMRES projection step if x 0 = 0.

Simulation Study
In this section the simulation study of Cockayne et al. [2018] will be replicated to demonstrate that the uncertainty produced from GMRES in Proposition 17 is similarly poorly calibrated, owing to the dependence of Q m on x * by way of its dependence on b. Throughout the size of the test problems is set to d = 100. The eigenvalues of A were drawn from an exponential distribution with parameter γ = 10, and eigenvectors uniformly from the Haar-measure over rotation-matrices (see Diaconis and Shahshahani [1987]). In contrast to Cockayne et al. [2018] the entries of b are drawn from a standard Gaussian distribution, rather than x * . By Lemma 1, the prior is then perfectly calibrated for this scenario, providing justification for the expectation that the posterior should be equally well-calibrated for m ≥ 1. Figure 6.4 shows on the left the convergence of GM-RES and on the right the convergence rate of the trace of the posterior covariance. Figure 6.4 repeats the uncertainty quantification study of Cockayne et al. [2018]. Cockayne et al. [2018] argue that if the uncertainty is well-calibrated then x * can be considered as a draw from the posterior. Under this assumption, i.e. Σ − 1 /2 m (x * − x m ) ∼ N (0, I) they derive the test statistic: It can be seen that the same poor uncertainty quantification occurs in BayesGMRES; even after just 10 iterations, the empirical distribution of the test statistic exhibits a profound left-shift, indicating an overly conservative posterior distribution. Producing well-calibrated posteriors remains an open issue in the field of probabilistic linear solvers.

Discussion
We have established many new connections between probabilistic linear solvers and a broad class of iterative methods. Matrix-based and solution-based inference were shown to be equivalent in a particular regime, showing that results from SBI transfer to MBI with left-multiplied observations. Since SBI is a special case of MBI, future research will establish what additional benefits the increased generality of MBI can provide.
We also established a connection between the wide class of projection methods and probabilistic linear solvers. The common practise of preconditioning has an intuitive probabilistic interpretation, and all probabilistic linear solvers can be interpreted as projection methods. While the converse was shown to hold, the conditions under which generic projection methods can be reproduced are somewhat restrictive; however, GM-RES and CG, which are among the most commonly used projection methods, have a well-defined probabilistic interpretation. Probabilistic interpretations of other widely used iterative methods can, we anticipate, be established from the results presented in this work.
Posterior uncertainty remains a challenge for probabilistic linear solvers. Direct probabilistic interpretations of CG and GMRES yield posterior covariance matrices which are not always computable, and even when the posterior can be computed the uncertainty remains poorly calibrated. This is owed to the dependence of the search directions in Krylov methods on Ax * = b, resulting in an algorithm which is not strictly Bayesian. Mitigating this issue without sacrificing the fast rate of convergence provided by Krylov methods remains an important focus for future work.
Definition 20 (symmetric Kronecker-product) The symmetric Kronecker-product for two square matrices A, B ∈ R N ×N of equal size is defined as Proposition 21 (Theorem 2.3 in Hennig [2015]) Let W ∈ R d×d be symmetric and positive definite. Assume a Gaussian prior of symmetric mean A −1 0 and covariance W ⊗ W on the elements of a symmetric matrix A −1 . After m linearly independent noise-free observations of the form S = A −1 Y , Y ∈ R d×m , rk(Y ) = m, the posterior belief over A −1 is a Gaussian with mean and posterior covariance where G := (Y W Y ) −1 .
Remark 22 Since A −1 0 is symmetric and the symmetric prior places mass only on symmetric matrices, the posterior mean A −1 m is also symmetric.

B.1 Proposition 6
Proof (Proof of Proposition 6) Let H = A −1 and let A −1 0 = H 0 . First note that by right-multiplying the information in Eq. (5) by H: Now note that where the second line uses Eq. (K2) and the third uses Eq. (K4). Thus where the second line is again using Eq. (K2) and Eq. (K4), while the third line uses Eq. (K3). We conclude that From these expressions it is straightforward to simplify the expressions for − − → H m : where the last line follows from K1. For Ω m : where the last line is from application of K5. It remains to project the posterior into R d by performing the matrix-vector product Hb.
where in the last line we have used that where in the second line we have used K2 and the fact that b W 0 b is a scalar, while in the third line we have used that b W 0 b = 1 and that Y m = A S m .
Note that x m =x m and Σ m =Σ m , as defined in Cockayne et al. [2018]. Thus, the proof is complete.

B.2 Theorem 16
Proof (Proof of Theorem 16.) Denote by x CG i the conjugate gradient estimate in iteration i and with p i the search direction in that iteration. From one iteration to the next, the update to the solution can be written as [Nocedal and Wright, 1999, p. 108] Comparing this update to lines 7 to 10 in Algorithm 1 it is sufficient to show that d i ∝ p i which follows from Lemma 23.
Lemma 23 Assume that CG does not terminate before d iterations. Using the prior of Theorem 16 in Algorithm 1, the directions d i are scaled conjugate gradients search directions, i.e.
Proof The proof proceeds by induction. Throughout we will suppress the superscript CG on the CG search directions, i.e. p CG i = p i . For i = 1, A −1 i−1 = αI by assumption and therefore d i = αr 0 which is the first CG search direction scaled by γ 1 = α = 0.
For the inductive step, suppose that the search directions s 1 , ..., s i−1 are scaled CG directions and that the vectors x 1 , . . . , x i−1 are the same as the first i − 1 solution estimates produced by CG. We will prove that s i is the i th CG search direction, and that x i is the i th solution estimate from CG. Lemma 25 states that d i can be written as where ν j ∈ R, j = 1, . . . , i. Under the prior, the posterior mean A −1 i is always symmetric as stated in Remark 22. This allows application of Lemma 24, so that {s 1 , . . . , s i−1 , d i } is an A-conjugate set. Thus we have, for < i: 0 = s Ad i = ν s As + ν i s Ar i−1 = ν s As + ν i y r i−1 .
This follows from Line 10 of Algorithm 1, from which it is clear that y = r −r −1 . Recall that the CG residuals r j are orthogonal [Nocedal and Wright, 1999, p. 109], and that from the inductive assumption, Algorithm 1 is equivalent to CG up to iteration i − 1). Thus, for < i − 1 we have that y r i−1 = 0 =⇒ s Ad i = ν s As = 0 ∀ < i − 1 where the second line is from application of the first line in Eq. (24). However, A is positive definite and by assumption the algorithm has not converged, so d = 0. Furthermore clearly s As = 0. Hence we must have that ν = 0 ∀ j < i − 1.
Equation (23) where the second line again applies the inductive assumption, that d i−1 and s i−1 are proportional to the CG search direction p i−1 , noting that the proportionality constants on numerator and denominator cancel. The term inside the brackets is precisely the i th CG search direction. This completes the result.
Lemma 24 If the belief over A −1 m is symmetric for all m = 0, . . . , d and A is symmetric and positive definite, then Algorithm 1 produces A-conjugate directions.
Proof The proof is by induction. Note that the case i = 1 is irrelevant since a set consisting of one element is trivially A-conjugate. On many occasions the proof relies on the consistency of the MBI belief, i.e. A −1 i z k = d k for k ≤ i and by symmetry z k A −1 i = d k . Thus, for the base case i = 2 we have: where the second line is by Line 10 of Algorithm 1. Now recall that α 1 = − d 1 r0 /d 1 Ad1 to give: = 0.
Here, the symmetry of the estimator A −1 i is used in Eq. (27). For the inductive step, assume {d 0 , . . . , d i−1 } are pairwise A-conjugate. For any k < i we have: where the second line follows from the fact that r i = r i−1 + y i . Thus, we have: Now, applying the conjugacy from the inductive assumption: where the second line rearranges line 6 of the algorithm to obtain α i d i z i = −d i r i−1 . The third line again uses that r i = r i−1 + y i , while the fourth line is from the assumed conjugacy.
Lemma 25 Under the prior in Theorem 16 and given scaled CG search directions p 1 , ..., p i , it holds that A −1 i r i ∈ span{p 1 , ..., p i , r i }.