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 the numerical error. In this paper we propose a novel
statistical model for this numerical 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