  This paper presents a probabilistic perspective on iterative methods for
approximating the solution $\mathbf{x}_* \in \mathbb{R}^d$ of a nonsingular
linear system $\mathbf{A} \mathbf{x}_* = \mathbf{b}$. In the approach a
standard iterative method on $\mathbb{R}^d$ is lifted to act on the space of
probability distributions $\mathcal{P}(\mathbb{R}^d)$. Classically, an
iterative method produces a sequence $\mathbf{x}_m$ of approximations that
converge to $\mathbf{x}_*$. The output of the iterative methods proposed in
this paper is, instead, a sequence of probability distributions $\mu_m \in
\mathcal{P}(\mathbb{R}^d)$. The distributional output both provides a "best
guess" for $\mathbf{x}_*$, for example as the mean of $\mu_m$, and also
probabilistic uncertainty quantification for the value of $\mathbf{x}_*$ when
it has not been exactly determined. Theoretical analysis is provided in the
prototypical case of a stationary linear iterative method. In this setting we
characterise both the rate of contraction of $\mu_m$ to an atomic measure on
$\mathbf{x}_*$ and the nature of the uncertainty quantification being provided.
We conclude with an empirical illustration that highlights the insight into
solution uncertainty that can be provided by probabilistic iterative methods.
