Accelerated Gibbs sampling of normal distributions using matrix splittings and polynomials

The Cholesky factorization is the conventional way to produce samples from a moderately sized multivariate normal distribution (Rue, 2001; Rue and Held, 2005), and is also the preferred method for solving moderately sized linear systems. For large linear systems, iterative solvers are the methods of choice due to their inexpensive cost per iteration, and small computer memory requirements.

Gibbs samplers applied to normal distributions are essentially identical to stationary iterative methods from numerical linear algebra. This connection was exploited by Adler (1981), and independently by Barone and Frigessi (1990), who noted that the component-wise Gibbs sampler is a stochastic version of the Gauss-Seidel linear solver, and accelerated the Gibbs sampler by introducing a relaxation parameter to implement the stochastic version of the successive over-relaxation (SOR) of Gauss-Seidel. This pairing was further analyzed by Goodman and Sokal (1989).

This equivalence is depicted in panels A and B of Figure 1. Panel B shows the contours of a normal density $\pi(\mathbf{x})$, and a sequence of coordinate-wise conditional samples taken by the Gibbs sampler applied to $\pi$. Panel A shows the contours of the quadratic minus $\log\left(\pi(\mathbf{x})\right)$ and the Gauss-Seidel sequence of coordinate optimizations¹¹1Gauss-Seidel optimization was rediscovered by Besag (1986) as iterated conditional modes., or, equivalently, solves of the normal equations $\nabla\log\pi(\mathbf{x})=0$. Note how in Gauss-Seidel the step sizes decrease towards convergence, which is a tell-tale sign that convergence (in value) is geometric. In Section 4 we will show that the iteration operator is identical to that of the Gibbs sampler in panel B, and hence the Gibbs sampler also converges geometrically (in distribution). Slow convergence of these algorithms is usually understood in terms of the same intuition; high correlations correspond to long narrow contours, and lead to small steps in coordinate directions and many iterations being required to move appreciably along the long axis of the target function.

Roberts and Sahu (1997) considered forward then backward sweeps of coordinate-wise Gibbs sampling, with relaxation parameter, to give a sampler they termed the REGS sampler. This is a stochastic version of the symmetric-SOR (SSOR) iteration, which comprises forward then backward sweeps of SOR.

The equality of iteration operators and error polynomials, for these pairs of fixed-scan Gibbs samplers and iterative solvers, allows existing convergence results in numerical analysis texts (for example Axelsson, 1996; Golub and Loan, 1989; Nevanlinna, 1993; Saad, 2003; Young, 1971) to be used to establish convergence results for the corresponding Gibbs sampler. Existing results for rates of distributional convergence by fixed-sweep Gibbs samplers (Adler, 1981; Barone and Frigessi, 1990; Liu et al., 1995; Roberts and Sahu, 1997) may be established this way.

The methods of Gauss-Seidel, SOR, and SSOR, give stationary linear iterations that were used as linear solvers in the 1950’s, and are now considered very slow. The corresponding fixed-scan Gibbs samplers are slow for precisely the same reason. The last fifty years has seen an explosion of theoretical results and algorithmic development that have made linear solvers faster and more efficient, so that for large problems, stationary methods are used as preconditioners at best, while the method of preconditioned conjugate gradients, GMRES, multigrid, or fast-multipole methods are the current state-of-the-art for solving linear systems in a finite number of steps (Saad and van der Vorst, 2000).

Linear iterations derived from a symmetric splitting may be sped up by polynomial acceleration, particularly Chebyshev acceleration that results in optimal error reduction amongst methods that have a fixed non-stationary iteration structure (Fox and Parker, 1968; Axelsson, 1996). The Chebyshev accelerated SSOR solver and corresponding Chebyshev accelerated SSOR sampler (Fox and Parker, 2014) are depicted in panels C and D of Figure 1. Both the solver and sampler take steps that are more aligned with the long axis of the target, compared to the coordinate-wise algorithms, and hence achieve faster convergence. However, the step size of Chebyshev-SSOR solving still decreases towards convergence, and hence convergence for both solver and sampler is still asymptotically geometric, albeit with much improved rate.

Fox and Parker (2014) considered point-wise convergence of the mean and variance of a Gibbs SSOR sampler accelerated by Chebyshev polynomials. In this paper we prove convergence in distribution for Gibbs samplers corresponding to any matrix splitting and accelerated by any polynomial that is independent of the Gibbs iterations. We then apply a polynomial accelerated sampler to solve a massive Bayesian linear inverse problem that is infeasible to solve using conventional techniques.

Chebyshev acceleration requires estimates of the extreme eigenvalues of the error operator, which we obtain via a conjugate-gradient (CG) algorithm at no significant computational cost (Meurant, 2006). The CG algorithm itself may be adapted to sample from normal distributions; the CG solver and corresponding sampler, depicted in panels E and F of Figure 1, were analysed by Parker and Fox (2012) and is discussed in the supplementary material.

In Section 2 we review efficient methods for sampling from normal distributions, highlighting Gibbs sampling in various algorithmic forms. Standard results for stationary iterative solvers are presented in Section 3. Theorems in Section 4 establish equivalence of convergence and convergence factors for iterative solvers and Gibbs samplers. Application of polynomial acceleration methods to linear solvers and Gibbs sampling is given in Section 5, including a proof of convergence of the first and second moments of a polynomial accelerated sampler. Numerical verification of convergence results is presented in Section 6.

Iterative samplers, such as Gibbs, are an attractive option when drawing samples from high dimensional multivariate normal distributions due to their inexpensive cost per iteration and small computer memory requirements (only vectors of size $n$ need be stored). If the precision matrix is sparse with $\mathcal{O}(n)$ non-zero elements, then, regardless of the bandwidth, iterative methods cost only about $2n$ flops per iteration, which is comparable with sparse Cholesky factorizations. However, when the bandwidth is $\mathcal{O}(n)$, the cost of the Cholesky factorization is high at $\mathcal{O}(n^{3})$ flops, while iterative methods maintain their inexpensive cost per iteration. Iterative methods are then preferable when requiring significantly fewer than $\mathcal{O}(n^{2})$ iterations for adequate convergence. In the examples presented in section 6 we find that $\mathcal{O}(n)$ iterations give convergence to machine precision, so the iterative methods are preferable for large problems.

A common class of methods for solving (2) are the linear iterative methods based on a splitting of $\mathbf{A}$ into $\mathbf{A}=\mathbf{M}-\mathbf{N}$. The matrix splitting is the standard way of expressing and analyzing linear iterative algorithms, following its introduction by Varga (1962). The system (2) is then transformed to $\mathbf{M}\mathbf{x}=\mathbf{N}\mathbf{x}+\mathbf{b}$ or, if $\mathbf{M}$ is nonsingular, $\mathbf{x}=\mathbf{M}^{-1}\mathbf{Nx}+\mathbf{M}^{-1}\mathbf{b}.$ The iterative methods use this equation to compute successively better approximations $\mathbf{x}^{\left(k\right)}$ to the solution using the iteration step

We follow the standard terminology used for these methods (see e.g. Axelsson, 1996; Golub and Loan, 1989; Saad, 2003; Young, 1971). Such methods are termed linear stationary iterative methods (of the first degree); they are stationary²²2This use of stationary corresponds to the term homogeneous when referring to a Markov chain. It is not to be confused with a stationary distribution that is invariant under the iteration. Later we will develop non-stationary iterations, inducing a non-homogeneous Markov chain that will, however, preserve the target distribution at each iterate. because the iteration matrix $\mathbf{G}=\mathbf{M}^{-1}\mathbf{N}$ and the vector $\mathbf{g}=\mathbf{M}^{-1}\mathbf{b}$ do not depend on $k$. The splitting is symmetric when both $\mathbf{M}$ and $\mathbf{N}$ are symmetric matrices. The iteration, and splitting, is convergent if $\mathbf{x}^{\left(k\right)}$ tends to a limit independent of $\mathbf{x}^{(0)}$, the limit being $\mathbf{A}^{-1}\mathbf{b}$ (see, e.g. (Young, 1971, Theorem 5.2)).

The iteration (3) is often written in the residual form so that convergence may be monitored in terms of the norm of the residual vector, and emphasizes that $\mathbf{M}^{-1}$ is acting as an approximation to $\mathbf{A}^{-1}$, as in Algorithm 4.

In computational algorithms it is important to note that the symbol $\mathbf{M}^{-1}\mathbf{r}$ is interpreted as “solve the system $\mathbf{Mu}=\mathbf{r}$ for $\mathbf{u}$” rather than “form the inverse of $\mathbf{M}$ and multiply $\mathbf{r}$ by $\mathbf{M}^{-1}$” since the latter is much more computationally expensive (about $2n^{3}$ flops (Watkins, 2002)). Thus, the practicality of a splitting depends on the ease with which one can solve $\mathbf{Mu}=\mathbf{r}$ for any vector $\mathbf{r}$.

A fundamental theorem of linear stationary iterative methods states that the splitting $\mathbf{A}=\mathbf{M}-\mathbf{N}$, where $\mathbf{M}$ is nonsingular, is convergent (i.e., $\mathbf{x}^{(k)}\to\mathbf{A}^{-1}\mathbf{b}$ for any $\mathbf{x}^{(0)}$) if and only if $\varrho\left(\mathbf{M}^{-1}\mathbf{N}\right)<1$, where $\varrho\left(\cdot\right)$ denotes the spectral radius of a matrix (Young, 1971, Theorem 3.5.1). This characterization is often used as a definition (Axelsson, 1996; Golub and Loan, 1989; Saad, 2003).

The error at step $k$ is $\mathbf{e}^{(k+1)}=\mathbf{x}^{(k+1)}-\mathbf{x}^{\ast}$, where $\mathbf{x}^{\ast}=\mathbf{A}^{-1}\mathbf{b}$. It follows that

and hence the asymptotic average reduction in error per iteration is the multiplicative factor

(Axelsson, 1996, p. 166). In numerical analysis this is called the (asymptotic average) convergence factor (Axelsson, 1996; Saad, 2003). Later, we will show that this is exactly the same as the quantity called the geometric convergence rate in the statistics literature (see e.g. Robert and Casella, 1999), for the equivalent Gibbs sampler. We will use the term ‘convergence factor’ throughout this paper to avoid a clash of terminology, since in numerical analysis the rate of convergence is minus the log of the convergence factor (see e.g. Axelsson, 1996, p. 166).

We now present the matrix splittings corresponding to some common stationary linear iterative solvers, with details for the case where $\mathbf{A}$ is symmetric, as holds for precision or covariance matrices.

We have seen that the Gauss-Seidel iteration uses the splitting $\mathbf{M}_{\text{GS}}=\mathbf{L}+\mathbf{D}$ and $N_{\text{GS}}=-\mathbf{L}^{T}$. Gauss-Seidel is one of the simplest splittings and solvers, but is also quite slow. Other splittings have been developed, though the speed of each method is often problem specific. Some common splittings are shown in Table 1, listed with, roughly, greater speed downwards. Speed of convergence in a numerical example is presented later in Section 6.

The method of successive over-relaxation (SOR) uses the splitting

in which $\omega$ is a relaxation parameter chosen with $0<\omega<2$. SOR with $\omega=1$ is Gauss-Seidel. For optimal values of $\omega$ such that $\varrho(\mathbf{M}^{-1}_{\text{SOR}}\mathbf{N}_{\text{SOR}})<\varrho(\mathbf{M}^{-1}_{\text{GS}}\mathbf{N}_{\text{GS}})$, SOR is an accelerated Gauss-Seidel iteration. Unfortunately, there is no closed form for the optimal value of $\omega$ for an arbitrary matrix $\mathbf{A}$, and the interval of values of $\omega$ which admits accelerated convergence can be quite narrow (Young, 1971; Golub and Loan, 1989; Saad, 2003).

The symmetric-SOR method (SSOR) incorporates both a forward and backward sweep of SOR so that if $\mathbf{A}$ is symmetric then the splitting is symmetric (Golub and Loan, 1989; Saad, 2003),

We will make use of symmetric splittings in conjunction with polynomial acceleration in Section 5.

When the matrix $\mathbf{A}$ is dense, Gauss-Seidel and SOR cost about $3n^{2}$ flops per iteration, with $2n^{2}$ due to multiplication by the matrix $\mathbf{A}$ (in order to calculate the residual) and another $n^{2}$ for the forward substitution to solve $\mathbf{Mu=r}$. Richardson incurs no cost to solve $\mathbf{Mu=r}$, while a solve with the diagonal Jacobi matrix incurs $n$ flops. Iterative methods are particularly attractive when the matrix $\mathbf{A}$ is sparse, since then the cost per iteration is only $O(n)$ flops.

Many theorems establish convergence of splittings by utilizing properties of $\mathbf{A}$ in specific applications. Some general conditions that guarantee convergence when $\mathbf{A}$ is SPD are given in the right column of Table 1 (Golub and Loan, 1989; Saad, 2003; Young, 1971).

The striking similarity between the Gibbs sampler (1) and the Gauss-Seidel iteration (4) is no coincidence. It is an example of a general equivalence between the stationary linear solver derived from a splitting and the associated stochastic iteration used as a sampler. We will make explicit the equivalence in the following theorems and corollary. In the first theorem we show that a splitting is convergent (in the sense of stationary iterative solvers) if and only if the associated stochastic iteration is convergent in distribution.

Convergence of the stochastic iteration (10) could also be established via the more general theory of Diaconis and Freedman (1999) that allows the iteration operator $\mathbf{G}=\mathbf{M}^{-1}\mathbf{N}$ to be random, with convergence in distribution guaranteed when $\mathbf{G}$ is contracting on average; see Diaconis and Freedman (1999) for details.

The following theorem shows how to design the noise distribution $\pi$ so that the limit distribution $\Pi$ has a desired mean $\mu$ and covariance ${\boldsymbol{\Sigma}}=\mathbf{A}^{-1}$.

Additionally, the mean and covariance converge geometrically, with convergence factors given by the convergence factors for the linear iterative method, as established in the following corollary.

Note that the matrix splitting has allowed an explicit construction of the noise covariance to give a desired precision matrix of the target distribution. We see from Theorem 2 that the stochastic iteration may be designed to converge to a distribution with non-zero target mean, essentially by adding the deterministic iteration (9) to the stochastic iteration (10). This is particularly useful when the mean is defined implicitly via solving a matrix equation. In cases where the mean is known explicitly, the mean may be added after convergence of the stochastic iteration with zero mean, giving an algorithm with faster convergence since the covariance matrix converges with factor $\varrho(\mathbf{M}^{-1}\mathbf{N})^{2}<\varrho(\mathbf{M}^{-1}\mathbf{N})$ (this was also noted by Barone et al., 2002). Convergence in variance for non-normal targets was considered in Fox and Parker (2014).

Using Theorems 1 and 2, and Corollary 3 we can draw on the vast literature in numerical linear algebra on stationary linear iterative methods to find random iterations that are computationally efficient and provably convergent in distribution with desired mean and covariance. In particular, results in Amit and Grenander (1991); Barone and Frigessi (1990); Galli and Gao (2001), Roberts and Sahu (1997), and Liu et al. (1995) are special cases of the general theory of matrix splittings presented here.

We now restrict attention to the case of normal target distributions.

Using Corollary 4, we found normal sampling algorithms corresponding to some common stationary linear solvers. The results are given in Table 2.

A sampler corresponding to a convergent splitting is implemented in Algorithm 5.

The practicality of a sampler derived from a convergent splitting depends on the ease with which one can solve $\mathbf{M}\mathbf{y}=\mathbf{r}$ for any $\mathbf{r}$ (as for the stationary linear solver) and also the ease of drawing iid noise vectors from $\operatorname{N}(\mathbf{0},\mathbf{M}^{T}+\mathbf{N})$. Sampling the noise vector is simple when a matrix square root, such as the Cholesky factorization, of $\mathbf{M}^{T}+\mathbf{N}$ is cheaply available. Thus, a sampler is at least as expensive as the corresponding linear solver since, in addition to operations in each iteration, the sampler must factor the $n\times n$ matrix $\mathbf{Var(c}^{(k)})=\mathbf{M}^{T}+\mathbf{N}$. For the samplers listed in Table 2 it is interesting that the simpler the splitting, the more complicated is the variance of the noise. Neither Richardson nor Jacobi splittings give useful sampling algorithms since the difficulty of sampling the noise vector is no less than the original task of sampling from $\operatorname{N}(\mathbf{0},\mathbf{A}^{-1})$. The Gauss-Seidel splitting, giving the usual Gibbs sampler, is at a kind of sweet spot, where solving equations in $\mathbf{M}$ is simple while the required noise variance is diagonal, so posing a simple sampling problem.

The SOR stationary sampler uses the SOR splitting $\mathbf{M}_{\text{SOR}}$ and $\mathbf{N}_{\text{SOR}}$ in (7) for $0<\omega<2$ and the noise vector $\mathbf{c}^{(k)}\sim\operatorname{N}(\mathbf{0},\mathbf{M}_{SOR}^{T}+\mathbf{N}_{SOR}=\frac{2-\omega}{\omega}\mathbf{D})$ (Table 2). This sampler was introduced by Adler (1981), rediscovered by Barone and Frigessi (1990), and has been studied extensively (Barone et al., 2002; Galli and Gao, 2001; Liu et al., 1995; Neal, 1998; Roberts and Sahu, 1997). For $\omega=1$, the SOR sampler is a Gibbs (Gauss-Seidel) sampler. For values of $\omega$ such that $\varrho(\mathbf{M}_{\text{SOR}}^{-1}\mathbf{N}_{\text{SOR}})<\varrho(\mathbf{M}^{-1}_{GS}\mathbf{N}_{GS})$), the SOR-sampler is an accelerated Gibbs sampler. As for the linear solver, implementation of the Gibbs and SOR samplers by Algorithm 5 requires multiplication by the upper triangular $\mathbf{N}$ and forward substitution with respect to $\mathbf{M}$ at a cost of $2n^{2}$ flops. In addition, these samplers must take the square root of the diagonal matrix $\frac{2-\omega}{\omega}\mathbf{D}$ at a mere cost of $O(n)$ flops.

Implementation of an SSOR sampler instead of a Gibbs or SOR sampler is advantageous since the Markov chain $\{\mathbf{y}^{(k)}\}$ is reversible (Roberts and Sahu, 1997). SSOR sampling uses the symmetric-SOR splitting $\mathbf{M}_{\text{SSOR}}$ and $\mathbf{N}_{\text{SSOR}}$ in (8). The SSOR stationary sampler is most easily implemented by forward and backward SOR sampling sweeps as in Algorithm 6, so the matrices $\mathbf{M}_{\text{SSOR}}$ and $\mathbf{N}_{\text{SSOR}}$ need never be explicitly formed.

We first encountered restricted versions of Corollary 4 for normal distributions in Amit and Grenander (1991) and in Barone and Frigessi (1990) where geometric convergence of the covariance matrices was established for the Gauss-Seidel and SOR splittings. These and the SSOR splitting were investigated in Roberts and Sahu (1997) (who labelled the sampler REGS).

Corollary 4 and Table 1 show that the Gibbs, SOR and SSOR samplers converge for any SPD precision matrix $\mathbf{A}$. This summarizes results in Barone and Frigessi (1990); Galli and Gao (2001) and the deterministic sweeps investigated in Amit and Grenander (1991); Roberts and Sahu (1997); Liu et al. (1995). Corollary 4 generalizes these results for any matrix splitting $\mathbf{A=M-N}$ by guaranteeing convergence of the random iterates (10) to $\operatorname{N}(\mathbf{0},\mathbf{A}^{-1})$ with convergence factor $\varrho(\mathbf{M}^{-1}\mathbf{N})$ (or $\varrho(\mathbf{M}^{-1}\mathbf{N})^{2}$ if $\mu=\mathbf{0}$).

A common scheme in numerical linear algebra for accelerating a stationary method when $\mathbf{M}$ and $\mathbf{N}$ are symmetric is through the use of polynomial preconditioners (Axelsson, 1996; Golub and Loan, 1989; Saad, 2003). Equation (5) shows that after $k$ steps the error in the stationary method is a $k^{th}$ order polynomial of the matrix $\mathbf{I}-\mathbf{G}=\mathbf{M}^{-1}\mathbf{A}$. The idea behind polynomial acceleration is to implicitly implement a different $k^{th}$ order polynomial $P_{k}(\mathbf{M}^{-1}\mathbf{A})$ such that $\varrho(P_{k}(\mathbf{M}^{-1}\mathbf{A}))<\varrho(\left(\mathbf{I}-\mathbf{M}^{-1}\mathbf{A}\right)^{k})$. The coefficients of $P_{k}(\mathbf{M}^{-1}\mathbf{A})$ are functions of a set of acceleration parameters $\{\{\alpha_{k}\},\{\tau_{k}\}\}$, introduced by the second order iteration

At the first step, $\alpha_{0}=1$ and $\mathbf{x}^{(1)}=\mathbf{x}^{(0)}+\tau_{0}\mathbf{M}^{-1}(\mathbf{b}-\mathbf{Ax}^{(0)})$. Setting $\alpha_{k}=\tau_{k}=1$ for every $k$ yields a basic un-accelerated stationary method. The accelerated iteration in (11) is implemented at a negligible increase in cost of $\mathcal{O}(n)$ flops per iteration (due to scalar-vector multiplication and vector addition) over the corresponding stationary solver (3).

It can be shown that (e.g., Axelsson (1996)) that the $(k+1)^{st}$ order polynomial $P_{k+1}$ generated recursively by the second order non-stationary linear solver (11) is

This polynomial acts on the error $\mathbf{e}^{(k)}=\mathbf{x}^{(k)}-\mathbf{A}^{-1}\mathbf{b}$ by $\mathbf{e}^{(k+1)}=P_{k}(\mathbf{M}^{-1}\mathbf{A})\mathbf{e}^{(0)}$, which can be compared directly to (5).

When estimates of the extreme eigenvalues $\lambda_{\min}$ and $\lambda_{\max}$ of $\mathbf{I}-\mathbf{G}=\mathbf{M}^{-1}\mathbf{A}$ are available ($\lambda_{\min}$ and $\lambda_{\max}$ are real when $\mathbf{M}$ and $\mathbf{N}$ are symmetric), then the coefficients $\{\tau_{k},\alpha_{k}\}$ can be chosen to generate the scaled Chebyshev polynomials $\{Q_{k}\}$, which give optimal error reduction at every step. The Chebyshev acceleration parameters are

where $\alpha_{0}=1$ and $\beta_{0}=\tau_{0}$ (Axelsson, 1996). Note that these parameters are independent of the iterates $\{\mathbf{x}^{(k)}\}$. Since $\mathbf{M}$ is required to be symmetric, applying Chebyshev acceleration to SSOR is a common pairing; its effectiveness as a linear solver is shown later in Table 3.

Whereas the stationary methods converge with asymptotic average convergence factor $\varrho(\mathbf{M}^{-1}\mathbf{N})$, the convergence factor for the Chebyshev method depends on $\operatorname{cond}(\mathbf{M}^{-1}\mathbf{A})=\lambda_{\max}/\lambda_{\min}$. Specifically the scaled Chebyshev polynomial $Q_{k}(\mathbf{\lambda})$ minimizes $\max_{\lambda\in[\lambda_{\min},\lambda_{\max}]}P_{k}(\lambda)$ over all $k^{th}$ order polynomials $P_{k}$, with

Since the error at the $k^{th}$ step of a Chebyshev accelerated linear solver is $\mathbf{e}^{(k+1)}=Q_{k}(\lambda)\mathbf{e}^{(0)}$, then the asymptotic convergence factor is bounded above by

(p.181 Axelsson, 1996). Since $\sigma\in[0,1)$, the polynomial accelerated scheme is guaranteed to converge even if the original splitting was not convergent. Further, the convergence factor of the stationary iterative solver is bounded below by $\rho=\frac{1-{\lambda_{\min}/\lambda_{\max}}}{1+{\lambda_{\min}/\lambda_{\max}}}$ (see e.g. Axelsson, 1996, Thm 5.9). Since $\sigma<\rho$ (except when $\lambda_{\min}=\lambda_{\max}$ in which case $\sigma=0$), polynomial acceleration always reduces the convergence factor, so justifies the term acceleration. The Chebyshev accelerated iteration (11) is amenable to preconditioning that reduces the condition number, and hence reduces $\sigma$, such as incomplete Cholesky factorization or graphical methods (Axelsson, 1996; Saad, 2003). Axelsson also shows that after

iterations of the Chebyshev solver, the error reduction is $||\mathbf{e}^{(k^{*})}||_{\mathbf{A}^{\nu}}/||\mathbf{e}^{(0)}||_{\mathbf{A}^{\nu}}\leq\varepsilon$ for some real number $\nu$ and any $0<\varepsilon<1$ (Axelsson, 1996, eqn 5.32).

Any acceleration scheme devised for a stationary linear solver is a candidate for accelerating convergence of a Gibbs sampler. For example, consider the second order stochastic iteration

analogous to the linear solver in (11) but now the vector $\mathbf{b}$ has been replaced by a random vector $\mathbf{c}^{(k)}$. The equivalence between polynomial accelerated linear solvers and polynomial accelerated samplers is made clear in the next three theorems.

Given a second order linear solver (11) that converges, Theorem 5 makes clear how to construct a second order sampler (17) that is guaranteed to converge. The next Corollary shows that the polynomial $P_{k}$ that acts on the linear solver error $\mathbf{x}^{(k)}-\mathbf{A}^{-1}\mathbf{b}$ is the same polynomial that acts on the errors in the first and second moments of the sampler, $\operatorname{E}(\mathbf{y}^{(k)})-\mathbf{A}^{-1}{\boldsymbol{\nu}}$ and $\operatorname{Var}(\mathbf{y}^{(k)})-\mathbf{A}^{-1}$ respectively. In other words, the convergence factors for a polynomial accelerated solver and sampler are the same.

Corollary 6 allows a direct comparison of the convergence factor for a polynomial accelerated sampler ($\sigma$, or $\sigma^{2}$ if $\mu=\mathbf{0}$) to the convergence factor given previously for the corresponding un-accelerated stationary sampler ($\varrho(\mathbf{M}^{-1}\mathbf{N})$, or $\varrho(\mathbf{M}^{-1}\mathbf{N})^{2}$ if $\mu=\mathbf{0}$). In particular, given a second order linear solver with accelerated convergence compared to the corresponding stationary iteration, the corollary guarantees that the second order Gibbs sampler (17) will converge faster than the stationary Gibbs sampler (10).

Just as Chebyshev polynomials are guaranteed to accelerate linear solvers, Corollary 6 assures that Chebyshev polynomials can also accelerate a Gibbs sampler. Using Theorem 5, we derived the Chebyshev accelerated SSOR sampler (Fox and Parker, 2014) by iteratively updating parameters via (13) and then generating a sampler via (17). Explicit implementation details of the Chebyshev accelerated sampler are provided in the supplementary materials. The polynomial accelerated sampler is implemented at a negligible increase in cost of ${\mathcal{O}}(n)$ flops per iteration over the cost ($4n^{2}$ flops) of the SSOR sampler (Algorithm 6). The asymptotic convergence factor is given by the next Corollary, which follows from Corollary 6 and equation (15).

Corollary 7 and (14) show that a Chebyshev accelerated normal sampler is guaranteed to converge faster than any other acceleration scheme that has the parameters $\{\{\tau_{k},\alpha_{k}\}\}$ independent of the iterates $\{\mathbf{y}^{(k)}\}$. This result also shows that the preconditioning ideas presented in section 5.1 to reduce ${\rm cond}(\mathbf{M}^{-1}\mathbf{A})=\lambda_{\max}/\lambda_{\min}$ can also be used to speed up Chebyshev accelerated samplers. We do not investigate such preconditioning here.

Corollary 7 and equation (16) suggest that, for any $\varepsilon>0$, after $k^{*}$ iterations the Chebyshev error reduction for the mean is smaller than $\varepsilon$. But even sooner, after $k^{**}=k^{*}/2$ iterations, the Chebyshev error reduction for the variance is predicted to be smaller than $\varepsilon$ (Fox and Parker, 2014).

A first order locally linear sparse precision matrix $\mathbf{A}$, considered by Higdon (2006); Rue and Held (2005), is

The discrete points $\{s_{i}\}$ are on a regular $10\times 10$ lattice ($n=100$) over the two dimensional domain ${\cal S}=[1,10]\times[1,10]$. Thus $\mathbf{A}$ is $100\times 100$, $||\mathbf{A}||_{2}=7.8$ and $||{\boldsymbol{\Sigma}}=\mathbf{A}^{-1}||_{2}=10^{4}$. The sparsity of $\mathbf{A}$ is shown in the left panel of Figure 2. The scalar $n_{i}$ is the number of points neighbouring $s_{i}$, i.e., with distance $1$ from $s_{i}$. Although $n^{2}=10^{4}$, the number of non-zero elements of $\mathbf{A}$ is $\mathcal{O}(n)$ ($460$ in this example). Since the bandwidth of $\mathbf{A}$ is $\mathcal{O}(n^{1/2})$, a Cholesky factorization costs $\mathcal{O}(n^{2})$ flops (Rue, 2001) and each iteration of an iterative method costs $\mathcal{O}(n)$ flops.

To provide a comparison between linear solvers and samplers, we solved the system $\mathbf{Ax}=\mathbf{b}$ using linear solvers with different matrix splittings (Table 1), where $\mathbf{b}$ is fixed and non-zero, all initialized with $\mathbf{x}^{(0)}=\mathbf{0}$. The results are given in Table 3. The Richardson method does not converge (DNC) since the spectral radius of the iteration operator is greater than 1. The SOR iteration was run at the optimal relaxation parameter value of $\omega=1.9852$. SSOR was run at its optimal value of $\omega=1.6641$. Chebyshev accelerated SSOR (Cheby-SSOR), CG accelerated SSOR (CG-SSOR) (both run with $\omega=1.6641$) and CG utilize a different implicit operator for each iteration, and so the spectral radius given in these cases is the geometric mean spectral radius of these operators (estimated using (5)). Even for this small example, Chebyshev acceleration reduces the computational effort required for convergence by about two orders of magnitude, while CG acceleration reduces work by nearly two more orders of magnitude.

We investigated the following Gibbs samplers: SOR, SSOR, and the Chebyshev accelerated SSOR. These samplers are guaranteed to converge since the corresponding solver converges (Theorem 1). Since the convergence factor for a sampler is equal to the convergence factor for the corresponding solver (Corollaries 3 and 6) then Gibbs samplers implemented with any of the matrix splittings in Table 1 exhibit the same convergence behavior as shown for the linear solvers in Table 3. Convergence of the sample covariance $\mathbf{S}^{(k)}_{y}\approx\mathbf{Var(y}^{(k)})\to\mathbf{A}^{-1}$, calculated using $10^{4}$ samples, is shown in the right panel of Figure 2 that displays the relative error $||\mathbf{A}^{-1}-\mathbf{S}^{(k)}_{y}||_{2}/||\mathbf{A}^{-1}||_{2}$ as a function of the flop count. Each sampler iteration costs about $2.24\times 10^{3}$ flops. This performance is compared to samples constructed by a Cholesky factorization, which cost $1.34\times 10^{4}$ flops (depicted as the green vertical line in the right panel of Figure 2). Since the sample means were uniformly close to zero, error in the mean is not shown.

The benchmark for evaluation of the of the convergence of the iterative samplers in finite precision is the Cholesky factorization, its relative error is depicted as the green horizontal line in the right panel of Figure 2. For this example, the iterative samplers produce better samples than a Cholesky sampler since the iterative sample covariances become more precise with more computing time.