The Adaptive $s$-step Conjugate Gradient Method

TL;DR

This paper introduces an adaptive $s$-step conjugate gradient method that automatically chooses the block size to balance convergence and efficiency, reducing synchronization costs on parallel computers without sacrificing accuracy.

Contribution

The paper develops a variable $s$-step CG method that adaptively determines block size based on residual growth bounds, improving reliability and usability over fixed $s$-step methods.

Findings

Achieves similar accuracy to classical CG

Reduces total global synchronizations

Automatically adapts block size for optimal performance

Abstract

On modern large-scale parallel computers, the performance of Krylov subspace iterative methods is limited by global synchronization. This has inspired the development of $s$-step Krylov subspace method variants, in which iterations are computed in blocks of $s$, which can reduce the number of global synchronizations per iteration by a factor of $O(s)$. Although the $s$-step variants are mathematically equivalent to their classical counterparts, they can behave quite differently in finite precision depending on the parameter $s$. If $s$ is chosen too large, the $s$-step method can suffer a convergence delay and a decrease in attainable accuracy relative to the classical method. This makes it difficult for a potential user of such methods - the $s$ value that minimizes the time per iteration may not be the best $s$ for minimizing the overall time-to-solution, and further may cause an unacceptable decrease in accuracy. Towards improving the reliability and usability of $s$-step Krylov subspace methods, in this work we derive the \emph{adaptive $s$-step CG method}, a variable $s$-step CG method where in block $k$, the parameter $s_k$ is determined automatically such that a user-specified accuracy is attainable. The method for determining $s_k$ is based on a bound on growth of the residual gap within block $k$, from which we derive a constraint on the condition numbers of the computed $O(s_k)$-dimensional Krylov subspace bases. The computations required for determining the block size $s_k$ can be performed without increasing the number of global synchronizations per block. Our numerical experiments demonstrate that the adaptive $s$-step CG method is able to attain up to the same accuracy as classical CG while still significantly reducing the total number of global synchronizations.

Tables

Table 1. Table 1: Test Matrix Properties.

Matrix	$n$	nnz	$\Vert A\Vert$	$\kappa \left(A\right)$
gr_30_30	900	7744	1.5	1.9e2
mesh3e1	289	1377	1.8	8.6e0
nos6	675	3255	2.0	3.5e6
bcsstk09	1083	18437	2.0	1.0e4
ex5	27	279	3.8	5.7e7

Equations

$${\mathcal{K}}_i(A,r_1)=\text{span}(r_1,A{r}_1,\dots ,A^{i-1}{r}_1)$$

$$p_{sk+j+1},r_{sk+j+1}$$

$$x_{sk+j+1}-x_{sk+1}$$

$${\mathcal{P}}_{k,s}$$

$${\mathcal{R}}_{k,s}$$

$${\rho }_0(z)$$

$${\rho }_j(z)$$

$$A{\underline{\mathcal{Y}}}_{k,s}={\mathcal{Y}}_{k,s}{\mathcal{B}}_{k,s},$$

$$p_{sk+j+1}$$

$$r_{sk+j+1}$$

$$x_{sk+j+1}-x_{sk+1}$$

$$x_{k,j+1}^{\prime }=x_{k,j}^{\prime }+{\alpha }_{sk+j}$$

$$\text{and}{p}_{k,j+1}^{\prime }$$

$${\alpha }_{sk+j}$$

$${\beta }_{sk+j}$$

$$\Vert b-A{\hat{x}}_i\Vert \le \Vert b-A{\hat{x}}_i-{\hat{r}}_i\Vert +\Vert {\hat{r}}_i\Vert .$$

$$A{\underline{\hat{\mathcal{Y}}}}_{k,s_k}$$

$$\text{where}\Vert E_k\Vert \le \epsilon \left((3+N_A)\Vert A\Vert \Vert {\underline{\hat{\mathcal{Y}}}}_{k,s_k}\Vert +4\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\mathcal{B}}_{k,s_k}\Vert \right),$$

$${\hat{x}}_{k,j+1}^{\prime }$$

$$\text{where}\Vert {\xi }_{k,j+1}\Vert \le \epsilon \Vert {\hat{x}}_{k,j+1}^{\prime }\Vert ,$$

$${\hat{r}}_{k,j+1}^{\prime }$$

$$\text{where}\Vert {\eta }_{k,j+1}\Vert \le \epsilon (\Vert {\hat{r}}_{k,j+1}^{\prime }\Vert +3\Vert {\mathcal{B}}_{k,s_k}\Vert \Vert {\hat{q}}_{k,j}^{\prime }\Vert ),$$

$${\hat{x}}_{m+j+1}$$

$$\text{where}\Vert {\varphi }_{m+j+1}\Vert \le \epsilon (\Vert {\hat{x}}_{m+j+1}\Vert +3j\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{x}}_{k,j+1}^{\prime }\Vert ),\text{ and}$$

$${\hat{r}}_{m+j+1}$$

$$\text{where}\Vert {\psi }_{m+j+1}\Vert \le 3j\epsilon \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{r}}_{k,j+1}^{\prime }\Vert ,$$

$${\delta }_{m+j+1}-{\delta }_{m+1}=-\sum _{i=1}^j\left(A{\underline{\hat{\mathcal{Y}}}}_{k,s_k}{\xi }_{k,i+1}+{\hat{\mathcal{Y}}}_{k,s_k}{\eta }_{k,i+1}-E_k{\hat{q}}_{k,i}^{\prime }\right)-A{\varphi }_{m+j+1}-{\psi }_{m+j+1}.$$

$$\Vert {\delta }_{m+j+1}-{\delta }_{m+1}\Vert$$

$$+\epsilon \left((3+N_A)\Vert A\Vert \Vert {\underline{\hat{\mathcal{Y}}}}_{k,s_k}\Vert +7\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\mathcal{B}}_{k,s_k}\Vert \right)\sum _{i=1}^j\Vert {\hat{q}}_{k,i}^{\prime }\Vert$$

$$+\epsilon \Vert A\Vert \Vert {\hat{x}}_{m+j+1}\Vert +3j\epsilon \Vert A\Vert \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{x}}_{k,j+1}^{\prime }\Vert +3j\epsilon \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{r}}_{k,j+1}^{\prime }\Vert$$

$$\le \epsilon \left((7+2N_A)\Vert A\Vert \Vert {\underline{\hat{\mathcal{Y}}}}_{k,s_k}\Vert +14\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\mathcal{B}}_{k,s_k}\Vert \right)\sum _{i=1}^j\Vert {\hat{x}}_{k,i+1}^{\prime }\Vert$$

$$+\epsilon \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \sum _{i=1}^j\Vert {\hat{r}}_{k,i+1}^{\prime }\Vert$$

$$+\epsilon \Vert A\Vert \Vert {\hat{x}}_{m+j+1}\Vert +3j\epsilon \Vert A\Vert \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{x}}_{k,j+1}^{\prime }\Vert +3j\epsilon \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{r}}_{k,j+1}^{\prime }\Vert ,$$

$$\Vert {\delta }_{m+j+1}-{\delta }_{m+1}\Vert$$

$$+\epsilon \Vert A\Vert \Vert {\hat{x}}_{m+j+1}\Vert +3j\epsilon \Vert A\Vert \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{x}}_{k,j+1}^{\prime }\Vert +3j\epsilon \Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{r}}_{k,j+1}^{\prime }\Vert .$$

$$\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{r}}_{k,i+1}^{\prime }\Vert$$

$$\le \kappa ({\hat{\mathcal{Y}}}_{k,s_k})\Vert {\hat{\mathcal{Y}}}_{k,s_k}{\hat{r}}_{m+i+1}-{\psi }_{m+i+1}\Vert$$

$$=\kappa ({\hat{\mathcal{Y}}}_{k,s_k})\Vert {\hat{r}}_{m+i+1}\Vert +O(\epsilon ),\text{and},$$

$$\Vert {\hat{\mathcal{Y}}}_{k,s_k}\Vert \Vert {\hat{x}}_{k,i+1}^{\prime }\Vert$$

$$\le \kappa ({\hat{\mathcal{Y}}}_{k,s_k})\Vert {\hat{\mathcal{Y}}}_{k,s_k}{\hat{x}}_{k,i+1}^{\prime }\Vert .$$

$${\hat{\mathcal{Y}}}_{k,s_k}{\hat{x}}_{k,i+1}^{\prime }$$

$$=(x-{\hat{x}}_{m+1})-(x-{\hat{x}}_{m+i+1})-{\varphi }_{m+i+1}$$

$$=A^{-1}(b-A{\hat{x}}_{m+1})-A^{-1}(b-A{\hat{x}}_{m+i+1})-{\varphi }_{m+i+1}$$

$$=A^{-1}({\hat{r}}_{m+1}+{\delta }_{m+1})-A^{-1}({\hat{r}}_{m+i+1}+{\delta }_{m+i+1})-{\varphi }_{m+i+1}$$

$$=A^{-1}({\hat{r}}_{m+1}-{\hat{r}}_{m+i+1})+A^{-1}({\delta }_{m+1}-{\delta }_{m+i+1})-{\varphi }_{m+i+1},$$

Full text

Although the $s$-step variants are mathematically equivalent to their classical counterparts, they can behave quite differently in finite precision depending on the parameter $s$. If $s$ is chosen too large, the $s$-step method can suffer a convergence delay and a decrease in attainable accuracy relative to the classical method. This makes it difficult for a potential user of such methods - the $s$ value that minimizes the time per iteration may not be the best $s$ for minimizing the overall time-to-solution, and further may cause an unacceptable decrease in accuracy.

Towards improving the reliability and usability of $s$-step Krylov subspace methods, in this work we derive the adaptive $s$-step CG method, a variable $s$-step CG method where in block $k$, the parameter $s_{k}$ is determined automatically such that a user-specified accuracy is attainable. The method for determining $s_{k}$ is based on a bound on growth of the residual gap within block $k$, from which we derive a constraint on the condition numbers of the computed $O(s_{k})$-dimensional Krylov subspace bases. The computations required for determining the block size $s_{k}$ can be performed without increasing the number of global synchronizations per block. Our numerical experiments demonstrate that the adaptive $s$-step CG method is able to attain up to the same accuracy as classical CG while still significantly reducing the total number of global synchronizations.

1 Introduction

In this paper, we consider the use of Krylov subspace methods for solving large, sparse linear systems $Ax=b$, where $A\in\mathbb{R}^{n\times n}$. We will focus on the conjugate gradient (CG) method [25], which is used when $A$ is symmetric positive definite. Given an initial approximate solution $x_{1}$ and corresponding residual $r_{1}=b-Ax_{1}$, the CG method iteratively updates the approximate solution using a sequence of nested Krylov subspaces $\mathcal{K}_{1}(A,r_{1})\subset\mathcal{K}_{2}(A,r_{1})\subset\cdots$, where

[TABLE]

denotes the $i$-dimension Krylov subspace with matrix $A$ and starting vector $v$. In iteration $i$, the updated approximate solution $x_{i+1}\in x_{1}+\mathcal{K}_{i}(A,r_{1})$ is chosen by imposing the Galerkin condition $r_{i+1}=b-Ax_{i+1}\perp\mathcal{K}_{i}$.

Thus each iteration of CG requires a matrix-vector product with $A$ in order to expand the dimension of the Krylov subspace and a number of inner products to perform the orthogonalization. On modern computer architectures, the speed with which the matrix-vector products and inner products can be computed is limited by communication (i.e., the movement of data). This limits the potential speed of individual iterations attainable by an implementation of CG. To perform a sparse matrix-vector product in parallel, each processor must communicate entries of the source vector and/or the destination vector that it owns to neighboring processors. Inner products require a global synchronization, i.e., the computation can not proceed until all processors have finished their local computation and communicated the result to other processors. For large-scale sparse problems on large-scale machines, the cost of synchronization between parallel processors can dominate the run-time (see, e.g., the exascale computing report [19, pp. 28]).

Research efforts toward removing the performance bottlenecks caused by communication in CG and other Krylov subspace methods have produced various approaches. One such approach are the $s$-step Krylov subspace methods (also called “communication-avoiding” Krylov subspace methods); for a thorough treatment of background, related work, and performance experiments, see, e.g., the theses [5, 27]. In $s$-step Krylov subspace methods, instead of performing one iteration at a time, the iterations are performed in blocks of $s$; i.e., in each iteration, the Krylov subspace is expanded by $O(s)$ dimensions by computing $O(s)$ new basis vectors and then all inner products between the new basis vectors needed for the next $s$ iterations are computed in one block operation. In this way, computing the inner products for $s$ iterations only requires a single global synchronization, decreasing the synchronization cost per iteration by a factor of $O(s)$. This approach has been shown to lead to significant speedups for a number of problems and real-world applications (see, e.g., [35, 43]). In the remainder of the paper, we will refer to the matrices whose columns consist of the $O(s)$ basis vectors computed in each block as $s$-step basis matrices. Further details of the $s$-step CG method are discussed in section 2. We emphasize that our use of the overloaded term “$s$-step methods” here differs from other works, e.g., [14, 30] and [21, §9.2.7], in which ‘$s$-step method’ refers to a type of restarted Lanczos procedure.

In exact arithmetic the $s$-step CG method produces the exact same iterates as the classical CG method, but their behavior can differ significantly in finite precision. In both $s$-step and classical Krylov subspace methods, rounding errors due to finite precision arithmetic have two basic effects: a decrease in attainable accuracy and a delay of convergence. It has long been known that for $s$-step Krylov subspace methods, as $s$ is increased (and so the condition numbers of the $s$-step basis matrices increase), the attainable accuracy decreases and the convergence delay increases relative to the classical CG method (see, e.g., [10]). At the extreme, if the parameter $s$ is chosen to be too large, the $O(s)$-dimensional bases computed for each block can be numerically rank deficient and the $s$-step method can fail to converge.

This sensitive numerical behavior poses a practical obstacle to optimizing the performance of $s$-step methods, and diminishes their usability and reliability. In a setting where the performance of CG is communication-bound, we expect that up to some point, increasing $s$ will decrease the time per iteration. If we pick $s$ only based on minimizing the time per iteration, however, we can run into problems. First, the finite precision error may cause a large convergence delay, negating any potential performance gain with respect to the overall runtime. Since the number of iterations required for convergence for a given $s$ value is not known a priori, choosing the $s$ value that results in the fastest time-to-solution is a difficult problem. Second, the chosen $s$ parameter may cause $s$-step CG to fail to converge to the user-specified accuracy. In this case, the particular problem is unsolvable by the $s$-step CG method.

Requiring the user to choose the parameter $s$ thus diminishes the practical usability of $s$-step Krylov subspace methods. It is therefore imperative that we develop a better understanding of the convergence rate and accuracy in finite precision $s$-step CG and other $s$-step Krylov subspace methods. Our hope is that by studying the theoretical properties of methods designed for large-scale computations in finite precision, we can develop methods and techniques that are efficient, capable of meeting application-dependent accuracy constraints, and which do not require that the user have extensive knowledge of numerical linear algebra.

Toward this goal, in this paper we develop a variable $s$-step CG method in which the $s$ values can vary between blocks and are chosen automatically by the algorithm such that a user-specified accuracy is attained. We call this the adaptive $s$-step CG method. The condition for determining $s$ in each block is derived using a bound on the size of the gap between the true residuals $b-Ax_{i}$ and the updated residuals $r_{i}$ computed in finite precision. We show that in order to prevent a loss of attainable accuracy, in a given block of iterations, the condition number of the $s$-step basis matrix should be no larger than a quantity related to the inverse of the 2-norms of the residuals computed within that block. In other words, as the method converges, more ill-conditioned bases (i.e., larger $s$ values) can be used without a loss of attainable accuracy. Furthermore, we show that the adaptive $s$-step CG method can be implemented in a way that does not increase the synchronization cost per block versus (fixed) $s$-step CG.

Whereas the fixed $s$-step CG method may not attain the same accuracy as classical CG, our numerical experiments on a set of example matrices demonstrate that our adaptive $s$-step CG method can achieve the same accuracy as classical CG, while still reducing the number of synchronizations required by up to over a factor of $4$.

In section 2 we give a brief background on $s$-step Krylov subspace methods and detail the $s$-step CG method, and in section 3 we discuss related work. Section 4 presents our main contribution; in section 4.1, we derive a theoretical bound on how large the $s$-step basis matrix condition number can be in terms of the norms of the updated residuals without affecting the attainable accuracy, and in section 4.2, we detail the adaptive $s$-step CG algorithm that makes use of this bound in determining the number of iterations to perform in each block. Numerical experiments for a small set of test matrices are presented in section 5. We discuss possible extensions and future work and briefly conclude in section 6.

2 Background: $s$-step Krylov subspace methods

The basic idea behind $s$-step Krylov subspace methods is to grow the underlying Krylov subspace $O(s)$ dimensions at a time and perform the necessary orthogonalization with a single global synchronization. Within a block of $s$ iterations, the vector updates are performed implicitly by updating the coordinates of the iteration vectors in the $O(s)$-dimensional basis spanned by the columns of the $s$-step basis matrix. The algorithmic details vary depending on the particular Krylov subspace method considered, but the general approach is the same. There is a wealth of literature related to $s$-step Krylov subspace methods. We give a brief summary of early related work below; a thorough bibliography can be found in [27, Table 1.1].

The $s$-step approach was first used in the context of Krylov subspace methods by Van Rosendale [41], who developed a variant of the parallel conjugate gradient method which minimizes inner product data dependencies with the goal of exposing more parallelism. The term “$s$-step Krylov subspace methods” was first used by Chronopoulos and Gear, who developed an $s$-step CG method [9, 10]. This work was followed in subsequent years by the development of other $s$-step variants of Orthomin and GMRES [11], Arnoldi and Symmetric Lanczos [31, 32], MINRES, GCR, and Orthomin [8, 13], Nonsymmetric Lanczos [33], and Orthodir [12]. Walker similarly used $s$-step bases, but with the goal of improving stability in GMRES by replacing the modified Gram-Schmidt orthogonalization process with Householder QR [42].

Many of the earliest $s$-step Krylov subspace methods used monomial bases (i.e., $[v,Av,A^{2}v,\ldots]$) for the computed $s$-step basis matrices, but it was found that convergence often could not be guaranteed for $s>5$ (see, e.g., [10]). This motivated research into the use of other more well-conditioned bases for the Krylov subspaces, including scaled monomial bases [26], Chebyshev bases [29, 16, 17], and Newton bases [1, 20]. The growing cost of communication in large-scale sparse problems has created a recent resurgence of interest in the implementation, optimization, and development of $s$-step Krylov subspace methods; see, e.g., the recent works [18, 27, 35, 46, 23, 34, 44, 45, 43].

2.1 The $s$-step conjugate gradient method

In this work, we focus on the $s$-step conjugate gradient method (Algorithm 2), which is equivalent to the $s$-step CG method (Algorithm 2) in exact arithmetic. In the remainder of this section, we give a brief explanation of the $s$-step CG method; further details on the derivation, implementation, and performance of the $s$-step CG method can be found in, e.g., [5, 27].

The $s$-step CG method consists of an outer loop, indexed by $k$, which iterates over the blocks of iterations, and an inner loop, which iterates over iterations $j\in\{1,\ldots,s\}$ within a block. For clarity, we globally index iterations by $i\equiv sk+j$. It follows from the properties of CG that for $j\in\{0,\ldots,s\}$,

[TABLE]

Then the CG vectors for the next $s$ iterations lie in the sum of the column spaces of the matrices

[TABLE]

where $\rho_{j}(z)$ is a polynomial of degree $j$ satisfying the three-term recurrence

[TABLE]

We define the $s$-step basis matrix $\mathcal{Y}_{k,s}=[\mathcal{P}_{k,s},\mathcal{R}_{k,s}]$ and define $\underline{\mathcal{Y}}_{k,s}$ to be the same as $\mathcal{Y}_{k,s}$ except with all zeros in columns $s+1$ and $2s+1$. Note that the $s$’s in the subscripts are unnecessary here, but will be useful in describing the variable $s$-step CG method later on. Under certain assumptions on the nonzero structure of $A$, $\mathcal{Y}_{k,s}$ can be computed in parallel in each outer loop for the same asymptotic latency cost as a single matrix-vector product using the so-called “matrix powers kernel” (see [18] for details). Since the columns of $\mathcal{Y}_{k,s}$ satisfy (3), we can write

[TABLE]

where $\mathcal{B}_{k,s}$ is a $(2s+1)\times(2s+1)$ block tridiagonal matrix defined by the 3-term recurrence coefficients in (3).

By (1), there exist vectors $p^{\prime}_{k,j+1}$, $r^{\prime}_{k,j+1}$, and $x^{\prime}_{k,j+1}$ that represent the coordinates of the CG iterates $p_{sk+j+1}$, $r_{sk+j+1}$, and $x_{sk+j+1}-x_{sk+1}$, respectively, in $\mathcal{Y}_{k,s}$ for $j\in\{0,\ldots,s\}$. That is,

[TABLE]

Therefore, in the inner loop of $s$-step CG, rather than update the CG vectors explicitly, we instead update their coordinates in $\mathcal{Y}_{k,s}$, i.e., for $j\in\{1,\ldots,s\}$,

[TABLE]

Thus the matrix-vector product with $A$ in each iteration becomes a matrix-vector product with the much smaller matrix $\mathcal{B}_{k,s}$. This along with the length-$(2s+1)$ vector updates can be performed locally by each processor in each inner loop iteration.

We can also reduce communication in computing the inner products. Letting $G_{k,s}=\mathcal{Y}_{k,s}^{T}\mathcal{Y}_{k,s}$, and using (5) and (4), $\alpha_{sk+j}$ and $\beta_{sk+j}$ can be computed by

[TABLE]

The matrix $G_{k,s}$ can be computed with one global synchronization per outer loop iteration, which, in terms of asymptotic parallel latency, costs the same as a single inner product. As $\mathcal{B}_{k,s}$ and $G_{k,s}$ are dimension $(2s+1)\times(2s+1)$, $\alpha_{sk+j}$ and $\beta_{sk+j}$ can be computed locally by each processor within the inner loop.

3 Related work

In this section, we discussed related work in the areas of variable $s$-step Krylov subspace methods, the analysis of maximum attainable accuracy in finite precision CG, and inexact Krylov subspace methods.

3.1 Variable $s$-step Krylov subspace methods

Williams et al. [43] use a variable $s$-step BICGSTAB as the coarse grid solver in a multigrid method. The technique used here, termed “telescoping”, is motivated by the observation than some coarse grid solves converge after only a few iterations, whereas other solves take significantly longer. By starting with a small $s$ value and gradually increasing it, they ensure that $s$-step BICGSTAB with larger $s$ is only used when the solve takes enough iterations to amortize the additional costs associated with $s$-step CG.

Imberti and Erhel [28] have recently developed an $s$-step GMRES method that allows variable $s$ sizes. They recommend choosing $s_{k}$ values according to the Fibonacci sequence up to some maximum value $s_{\text{max}}$, i.e., starting with $s_{0}=s_{1}=1$ and for $k\geq 2$, $s_{k}=\min{(s_{\text{max}},s_{k-1}+s_{k-2})}$. In this way, the sequence $\{s_{i}\}$ used by Imberti and Erhel is predetermined.

In contrast, in our approach, the sequence $\{s_{i}\}$ is dynamically chosen based on the 2-norm of the updated residual (which is not necessarily monotonically decreasing in CG). In addition, our method is designed such that a user-specified accuracy $\varepsilon^{*}$ can be attained when $\varepsilon^{*}\geq\varepsilon_{\text{CG}}$, where $\varepsilon_{\text{CG}}$ is the accuracy attained by classical CG for the same problem.

3.2 Maximum attainable accuracy in Krylov subspace methods

In both $s$-step and classical variants of Krylov subspace methods, finite precision roundoff error in updates to the approximate solution $x_{i}$ and the residual $r_{i}$ in each iteration can cause the updated residual $r_{i}$ and the true residual $b-Ax_{i}$ to grow further and further apart as the iterations proceed. If this deviation grows large, it can limit the maximum attainable accuracy, i.e., the accuracy with which we can solve $Ax=b$ on a computer with unit round-off $\varepsilon$. Analyses of maximum attainable accuracy in CG and other classical KSMs are given by Greenbaum [22], van der Vorst and Ye [40], Sleijpen, van der Vorst, and Fokkema [37], Sleijpen, van der Vorst, and Modersitzki [38], Björck, Elfving, and Strakoš [2], and Gutknecht and Strakoš [24]. One important result of these analyses is the insight that loss of accuracy can be caused at a very early stage of the computation, which can not be corrected in later iterations. Analyses of the maximum attainable accuracy in $s$-step CG and the $s$-step biconjugate gradient method (BICG) can be found in [6, 5].

3.3 Inexact Krylov subspace methods

The term “inexact Krylov subspace methods” refers to Krylov subspace methods in which the matrix-vector products $Av$ are computed inexactly as $(A+E)v$ where $E$ represents some error matrix (due to either finite precision computation or intentional approximation). It is assumed that all other computations (orthogonalization and vector updates) are performed exactly. Using an analysis of the deviation of the true and updated residuals, it is shown that under these assumptions the size of the perturbation term $E$ can be allowed to grow inversely proportional to the norm of the updated residual without adversely affecting the attainable accuracy. This theory has potential application in mixed-precision Krylov subspace methods, as well as in situations where the operator $A$ is expensive to apply and can be approximated in order to gain performance.

Much work has been dedicated to inexact Krylov subspace methods. In [36], Simoncini and Szyld provide a general theory for inexact variants applicable to both symmetric and nonsymmetric linear systems and eigenvalue problems and use this to derive criteria for bounding the size of $E$ in such a way that does not diminish the attainable accuracy of the method. A similar analysis was given by van de Eschof and Sleijpen  [39]. These analyses confirm and improve upon the earlier empirical results for the inexact GMRES method of Bouras and Fraysseé [3] and the inexact conjugate gradient method of Bouras, Frayssé, and Giraud [4].

Our work is very closely related to the theory of inexact Krylov subspace methods in that our analysis results in a somewhat analogous “relaxation strategy”, where instead of relaxing the accuracy of the matrix-vector products we are instead relaxing a constraint on the condition numbers of the computed $s$-step basis matrices. In addition, our analysis assumes that all parts of the computation are performed in a fixed precision $\varepsilon$.

4 Variable $s$-step CG

As discussed, the conditioning of the $s$-step basis matrices in $s$-step Krylov subspace methods affects the rate of convergence and attainable accuracy. In the $s$-step CG method (Algorithm 2), the Krylov subspace basis is computed in blocks of size $s$ in each outer loop. This method can be generalized to allow the blocks to be of varying size; in other words, a different $s$ value can be used in each outer loop iteration. We denote the block size in outer iteration $k$ by $s_{k}$, where $1\leq s_{k}\leq s_{\text{max}}$ for some maximum $s$ value $s_{\text{max}}$. We call this a variable $s$-step Krylov subspace method. A general variable $s$-step CG method is shown in Algorithm 3.

In this section, we derive a variable $s$-step CG method which automatically determines $s_{k}$ in each outer loop such that a user-specified accuracy can be attained, which we call the adaptive $s$-step CG method. Our analysis is based on the maximum attainable accuracy analysis for $s$-step CG in [6], which shows that the attainable accuracy of $s$-step CG depends on the condition numbers of the $s$-step basis matrices $\mathcal{Y}_{k,s}$ computed in each outer loop. In summary, we show that if in outer loop $k$, $s_{k}$ is selected such that the condition number of the basis matrix $\mathcal{Y}_{k,s_{k}}$ is inversely proportional to the maximum 2-norm of the residual vectors computed within outer loop $k$, then the approximate solution produced by the variable $s$-step CG method can be as accurate as the approximate solution produced by the classical CG method. Further, our method allows for the user to specify the desired accuracy, so in the case that a less accurate solution is required, our method automatically adjusts to allow for higher $s_{k}$ values. This effectively exploits the tradeoff between accuracy and performance in $s$-step Krylov subspace methods.

In this section we derive our adaptive $s$-step approach, which stems from a bound on the growth in outer loop $k$ of the gap between the true and updated residuals in finite precision.

4.1 A constraint on basis matrix condition number

In the remainder of the paper, hats denote quantities computed in finite precision, $\varepsilon$ denotes the unit round-off, $\|\cdot\|$ denotes the 2-norm, and $\kappa(A)$ denotes the 2-norm condition number $\|A^{-1}\|\|A\|$. To simplify the analysis, we drop all $O(\varepsilon^{2})$ terms (and higher order terms in $\varepsilon$). Writing the true residual $b-A\widehat{x}_{i}=b-A\widehat{x}_{i}-\widehat{r}_{i}+\widehat{r}_{i}$, we can write the bound

[TABLE]

Then as $\|\widehat{r}_{i}\|\rightarrow 0$ (i.e., as the method converges), we have $\|b-A\widehat{x}_{i}\|\approx\|b-A\widehat{x}_{i}-\widehat{r}_{i}\|$. The size of the residual gap, i.e., $\|b-A\widehat{x}_{i}-\widehat{r}_{i}\|$, therefore determines the attainable accuracy of the method (see the related references in section 3.2).

We begin by considering the rounding errors made in the variable $s$-step CG method (Algorithm 3). In outer loop $k$ and inner loops $j\in\{1,\ldots,s_{k}\}$ of finite precision variable $s$-step CG, using standard models of floating point error (e.g., [21, §2.4]) we have

[TABLE]

where $N_{A}$ denotes the maximum number of nonzeros per row in $A$ and $m=\sum_{i=0}^{k-1}s_{i}$. Let $\delta_{i}\equiv b-A\widehat{x}_{i}-\widehat{r}_{i}$ denote the residual gap in iteration $i$. Similar to the analysis in [6], for variable $s$-step CG we can write the growth of the residual gap in iteration $m+j+1$ in outer loop $k$ as

[TABLE]

Our goal is now to manipulate a bound on this quantity in order to derive a method for determining the largest $s_{k}$ value we can use in outer loop $k$ such that the desired accuracy is attained. Using standard techniques, we have the norm-wise bound

[TABLE]

where we have used $\sum_{i=1}^{j}\|\widehat{q}^{\prime}_{k,i}\|\leq(2+\varepsilon)\sum_{i=1}^{j}\|\widehat{x}^{\prime}_{k,i+1}\|$. To simplify the notation, we let $\tau_{k}=\|\mathcal{B}_{k,s_{k}}\|/\|A\|$ and $N_{k}=7+2N_{A}+14\tau_{k}$. We note that $\|\mathcal{B}_{k,s_{k}}\|$ depends on the polynomial basis used and is not expected to be too large; e.g., for the monomial basis, $\|\mathcal{B}_{k,s_{k}}\|=1$. The above bound can then be written

[TABLE]

Note that assuming $\widehat{\mathcal{Y}}_{k,s_{k}}$ is full rank, we have

[TABLE]

To bound $\|\widehat{\mathcal{Y}}_{k,s_{k}}\widehat{x}^{\prime}_{k,i+1}\|$ in terms of the size of the residuals, notice that by (9),

[TABLE]

so

[TABLE]

Thus together with (13) we can write

[TABLE]

This allows us to write the bound (11) in the form

[TABLE]

Using $x-\widehat{x}_{m+j+1}=A^{-1}\widehat{r}_{m+j+1}+\delta_{m+j+1}$, we have

[TABLE]

and writing $\widehat{x}_{m+j+1}=\widehat{x}_{m+j+1}-x+x$, we have

[TABLE]

Then

[TABLE]

and letting $c_{k,j}=6N_{k}j\kappa(A)$, this bound can be written

[TABLE]

This gives a norm-wise bound on the growth of the residual gap due to finite precision errors made in outer iteration $k$.

Letting $m_{\ell}=\sum_{i=0}^{\ell-1}s_{i}$, notice that we have

[TABLE]

Suppose that at the end of the iterations we want the size of the residual gap to be on the order $\varepsilon^{*}$ (where $\varepsilon^{*}\gtrsim O(\varepsilon)\|A\|\|x\|$, as this is the best we could expect). Then assuming the number of outer iterations $k$ is not too high, this can be accomplished by requiring that in each outer loop, for $j\in\{1,\ldots,s_{k}\}$,

[TABLE]

From (15), this means that we must have, for all inner iterations $j\in\{1,\ldots,s_{k}\}$,

[TABLE]

The bound (16) tells us that as the 2-norm of the residual decreases, we can tolerate a more ill-conditioned basis matrix $\widehat{\mathcal{Y}}_{k,s}$ without detriment to the attainable accuracy. Since $\kappa(\widehat{\mathcal{Y}}_{k,s})$ grows with increasing $s_{k}$, this suggests that $s_{k}$ can be allowed to increase as the method converges at a rate proportional to the inverse of the residual 2-norm. This naturally gives a relaxation strategy that is somewhat analogous to those derived for inexact Krylov subspace methods (see section 3.3). We discuss the implementation of this strategy in the following subsection.

4.2 The adaptive $s$-step CG method

In this section, we propose a variable $s$-step CG method that makes use of the constraint (16) in determining $s_{k}$ in each outer loop $k$, which we call the adaptive $s$-step CG method. Our approach uses the largest $s_{k}$ value possible in each outer iteration $k$ such that (16) is satisfied up to some $s_{\text{max}}$ set by the user (e.g., $s_{\text{max}}$ could be selected by auto-tuning to find the $s$ value that minimizes the time per iteration in $s$-step CG). Since a variable $s$-step CG method produces the same iterates as classical CG when $s_{k}=1$ for all outer loops $k$, we expect that our adaptive $s$-step CG method can attain accuracy $\varepsilon^{*}$ when $\varepsilon^{*}\geq\varepsilon_{\text{CG}}$, where $\varepsilon_{\text{CG}}$ is the attainable accuracy of classical CG. When the right-hand-side of (16) is less than one (due to, e.g., $\varepsilon^{*}$ being set too small), the inequality (16) cannot be satisfied. In this case, we default to using $s_{k}=1$ (i.e., we perform an iteration of classical CG).

In CG, it is the $A$-norm of the error rather than the $2$-norm of the residual that is minimized, which means that in some inner loop iteration $j$, we can have $\|\widehat{r}_{m+j+1}\|>\|\widehat{r}_{m+1}\|$. Since $\|\widehat{r}_{m+i+1}\|$ for $i>1$ are not known at the start of the outer loop iteration, we must determine $s_{k}$ online, i.e., during the iterations. We therefore first construct a basis matrix that satisfies (16) under the assumption that $\|\widehat{r}_{m+i+1}\|\leq\|\widehat{r}_{m+1}\|$. Within each inner loop, we can perform an inexpensive check to ensure that (16) is still satisfied after computing the next residual; if not, we break from the inner loop, discard the current basis matrix, and begin a new outer loop iteration. We elaborate on this approach in the remainder of this section. An optimized high-performance implementation of our variable $s$-step approach is outside the scope of the present work, although we will briefly discuss how one might implement this strategy efficiently such that each outer loop still requires only a single global synchronization and the number of wasted matrix-vector products is kept to a minimum.

The bound (16) indicates that when $\widehat{\mathcal{Y}}_{k,s_{k}}$ is constructed at the beginning of outer loop $k$, we should enforce

[TABLE]

Therefore we first construct an $\bar{s}_{k}$-step basis matrix $\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}$ using some $\bar{s}_{k}\leq s_{\text{max}}$, compute the Gram matrix $\widehat{G}_{k,\bar{s}_{k}}=\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}^{T}\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}$, and then use the condition numbers of the leading principle submatrices of the blocks of $\widehat{G}_{k,\bar{s}_{k}}$ to compute the condition numbers of the leading columns of the blocks of $\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}$. To be concrete, for $i\in\{1,2,\ldots,\bar{s}_{k}\}$, we have (using MATLAB notation)

[TABLE]

and

[TABLE]

where $i_{s}=\bar{s}_{k}+2$. Then since $\sqrt{\kappa(\widehat{G}_{k,i})}=\kappa(\widehat{\mathcal{Y}}_{k,i})$, we can easily compute the condition number of the $i$-step basis matrix $\widehat{\mathcal{Y}}_{k,i}$, i.e., the basis matrix needed to perform $i$ inner loop iterations, for $i\in\{1,\ldots,\bar{s}_{k}\}$. We therefore let $\tilde{s}_{k}\in\{1,\ldots,\bar{s}_{k}\}$ be the largest value possible such that

[TABLE]

holds.

To handle the case where $\|\widehat{r}_{m+j+1}\|>\|\widehat{r}_{m+1}\|$ for some inner iteration $j$, we can add a check within the inner loop iterations: after computing $\widehat{r}^{\prime}_{k,j+1}$, we check if

[TABLE]

is still satisfied. Note that since $\|\widehat{r}_{m+j+1}\|=\left(\widehat{r}_{k,j+1}^{{}^{\prime}T}\widehat{G}_{k,\tilde{s}_{k}}\widehat{r}^{\prime}_{k,j+1}\right)^{1/2}+O(\varepsilon)$, performing this check does not require any additional communication. If (19) is satisfied for all $j\in\{1,\ldots,\tilde{s}_{k}\}$, we will perform all $s_{k}=\tilde{s}_{k}$ inner loop iterations. If this is not satisfied for some $j$, we break from the current inner loop iteration and begin the next outer loop; i.e., we only perform $s_{k}=j$ inner loop iterations in outer loop $k$.

Although performing the checks for determining when to break from the current outer loop do not incur any additional communication cost in terms of the number of messages sent between processors, there is, however, potential wasted effort in terms of computation and the number of words moved when $\bar{s}_{k}>\tilde{s}_{k}$ and/or $\tilde{s}_{k}>s_{k}$, since in this case we have computed more basis vectors than we need to perform $s_{k}$ inner loop updates. One way to mitigate this is to only allow $\bar{s}_{k}$ to grow by at most $f$ vectors in each outer loop, i.e., $\bar{s}_{k}\leq s_{k-1}+f\leq s_{\text{max}}$, where $f$ is some small positive integer (e.g., $f=1$ or $f=2$).

The approach outlined above is summarized in Algorithm 4. Note that in Algorithm 4, we have specified that the function $c_{k,j}$ can be input by the user. Because our bounds are not tight, in many cases using $c_{k,j}=6N_{k}j\kappa(A)$ (from the bound (16)) overestimates the error, which results in the use of smaller $s_{k}$ values than necessary.

Note that although there is no expectation that $\|\widehat{r}_{m+j+1}\|$ decreases between iterations in CG, as long as the method is converging, we still expect the residual norm to trend downward. Therefore one could make the simplifying assumption that $\|\widehat{r}_{m+1}\|\geq\|\widehat{r}_{m+j+1}\|$ for $j\in\{1,\ldots,\tilde{s}_{k}\}$. In this case, the value $s_{k}=\tilde{s}_{k}$ can be determined at the beginning of outer loop $k$ before performing any inner iterations, and so the checks (19) are unnecessary. This simplification comes at the potential cost of reduced reliability of the adaptive $s$-step method, although in our experimenting we have observed that taking $s_{k}=\tilde{s}_{k}$ usually works well even when the residual norm is nonmonotonic.

Because our strategy is based on a loose upper bound on the deviation of residuals, we do not require particularly accurate estimates of $\kappa(\widehat{\mathcal{Y}}_{k,i})$. Therefore alternative methods for inexpensively approximating this quantity could also be used in practice.

5 Numerical experiments

In this section, we test the adaptive $s$-step CG method (Algorithm 4) on a number of small examples from the University of Florida Sparse Matrix Collection [15] using MATLAB. The test examples were chosen to exemplify various convergence trajectories of CG, which allows us to demonstrate the validity of our relaxation condition in which $s_{k}$ depends on the rate of convergence. The test problems used here are small, so we would not use (fixed or variable) $s$-step methods here in practice; nevertheless, they still serve to exhibit the numerical behavior of the methods.

All tests were run in double precision, i.e., $\varepsilon\approx 2^{-53}$, and use a right-hand side $b$ with entries $1/\sqrt{n}$ where $n$ is the dimension of $A$. We use matrix equilibration, i.e., $A\leftarrow D^{-1/2}AD^{-1/2}$ where $D$ is a diagonal matrix of the largest entries in each row of $A$. Properties of the test matrices (post-equilibration) are shown in Table 1.

Our results are plotted in Figures 2-6. For each test matrix, we ran classical CG (Algorithm 1), as well as both (fixed) $s$-step CG (Algorithm 2) and adaptive $s$-step CG (Algorithm 4) for $s_{\text{max}}$ values $4$, $8$, and $10$ (for fixed $s$-step CG, $s=s_{\text{max}}$). For each matrix and $s_{\text{max}}$ value, we tested two different $\varepsilon^{*}$ values: $\varepsilon^{*}=\varepsilon_{\text{CG}}$, the maximum attainable accuracy of classical CG, and $\varepsilon^{*}=1e-6$. The monomial basis was used for both variable and fixed $s$-step methods in all tests (note that there is no technical restriction to the monomial basis; the analysis used in deriving the adaptive $s$-step CG method still applies regardless of the polynomials used in constructing $\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}$). In all tests for adaptive $s$-step CG, we set $\bar{s}_{0}=f=s_{\text{max}}$ (see Algorithm 4), as this gives the optimal choice of $s_{k}$ in each iteration. (We also tried using smaller values of $f$, e.g., $f=1$ and $f=2$; this did not significantly change the total number of outer loop iterations performed). We used $c_{k,j}=1$ in our condition for determining the number of inner loop iterations to perform in all experiments except for tests with matrices bcsstk09 (Figure 5) and ex5 (Figure 6); for these problems, convergence is more irregular, and thus the bound (16) is tighter.

Each figure has a corresponding table (Tables 2-6) which lists the number of outer loop iterations (a proxy for the latency cost or number of global synchronizations) required for convergence of the true residual $2$-norm to level $\varepsilon^{*}$. For classical CG, the number reported gives the total number of iterations. Dashes in the tables indicate that the true residual diverged or stagnated before reaching level $\varepsilon^{*}$.

5.1 Observations

In all experiments, the adaptive $s$-step CG method was able to converge to the desired tolerance $\varepsilon^{*}$, whereas for $\varepsilon^{*}=\varepsilon_{\text{CG}}$ and higher $s$ values ($s=8$ and/or $s=10$), fixed $s$-step CG may not achieve convergence to this level. We point out that even when the matrix is very well-conditioned, choosing $s$ too large can very negatively affect the convergence rate in fixed $s$-step CG; see, e.g., the case with $s=10$ and $\varepsilon^{*}=1e-6$ in Figure 3, where $s$-step CG requires more than 8 times the number of synchronizations than classical CG. In tests where both fixed $s$-step CG and adaptive $s$-step CG converge, adaptive $s$-step CG takes at most $2$ additional outer iterations versus fixed $s$-step CG (although in most cases it takes as many or fewer outer iterations). We note that in all cases, variable $s$-step CG uses fewer outer loop iterations (synchronization points) than classical CG, so we still expect a performance advantage on latency-bound problems.

The experiments using $\varepsilon^{*}=1e{-}6$ demonstrate the ability of the adaptive $s$-step CG method to automatically adjust to use the largest $s_{k}$ values possible (according to our bound) that will eventually attain the desired accuracy. For example, in Figure 2 for gr_30_30, for all tests using $\varepsilon^{*}=1e{-}6$, adaptive $s$-step CG uses $s_{k}=s_{\text{max}}$ in every outer iteration, defaulting to the fixed $s$-step method. In contrast, in plots on the left using $\varepsilon^{*}=\varepsilon_{\text{CG}}$, adaptive $s$-step CG uses smaller $s_{k}$ values at the beginning and then increases up to $s_{k}=s_{\text{max}}$ and is able to attain accuracy $\varepsilon_{\text{CG}}$ whereas fixed $s$-step CG fails when for $s=8$ and $s=10$.

Below we list the sequence of $s_{k}$ values chosen in adaptive $s$-step CG for tests with $\varepsilon^{*}=\varepsilon_{\text{CG}}$ and $s_{\text{max}}=10$ (plots in the lower left in each figure):

• •

gr_30_30: 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6, 8, 10, 10

• •

mesh3e1: 1, 1, 2, 4, 6, 9, 10

• •

nos6: 6, 1, 2, 3, 4, 5, 4, 4, 1, 4, 4, 5, 6, 7, 8, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10

• •

bcsstk09: 2, 2, 3, 3, 3, 3, 3, 1, 2, 3, 3, 3, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 7, 1, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10

• •

ex5: 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 2, 1, 1, 1, 3, 3, 2, 1, 1, 1, 2, 1, 1, 2, 4, 6, 2, 4, 5, 7, 2, 6, 1, 7, 6

There are a few things to notice. First, as expected, the sequence of $s_{k}$ values is monotonically increasing when the residual 2-norm decreases monotonically or close to it (matrices gr_30_30 and mesh3e1), but when the residual 2-norm oscillates, we no longer have monotonicity of the $s_{k}$ values (matrices nos6, bcsstk09, and ex5). Second, it may be that $s_{\text{max}}$ is never reached, as is the case for matrix ex5. It is also worth pointing out that the sequence of $s_{k}$ values need not start with $s_{0}=1$; for nos6, $s_{0}=6$ and for bcsstk09, $s_{0}=2$.

It may be that on a particular machine for a particular matrix structure, it is not beneficial in terms of minimizing time per iteration to go above a small value of $s_{\text{max}}$, say $s_{\text{max}}=4$. In this case, as long as the matrix is reasonably well-conditioned and the accuracy required is not too high, the fixed $s$-step CG method can probably be used without too much trouble. However, even in this case the adaptive $s$-step CG method can provide reliability while still maintaining the reduction in synchronizations of the fixed $s$-step method, at the cost of potentially wasting a few matrix-vector multiplications in each outer iteration.

Looking at the number of outer loop iterations in Tables 2-6, we see that in our tests for adaptive $s$-step CG it is often the case that increasing the $s_{\text{max}}$ value past a certain point gives diminishing returns in terms of reducing synchronization. This may represent a fundamental limit of the possible performance using the monomial basis for these problems. In these cases it may be that better performance (higher $s_{k}$ values) can be achieved only by using a better-conditioned polynomial basis for constructing $\widehat{\mathcal{Y}}_{k,\bar{s}_{k}}$ (e.g., Newton or Chebyshev bases).

6 Conclusions and future work

In this work, we developed the adaptive $s$-step CG method, which is a variable $s$-step CG method where the parameter $s_{k}$ is determined automatically in each outer loop such that a user-specified accuracy is attainable. The method for determining $s_{k}$ is based on a bound on the residual gap within outer loop $k$, from which we derive a constraint on the condition number of the $s$-step basis matrix generated in outer loop $k$. The computations required for determining $s_{k}$ can be performed without introducing any additional synchronizations in each outer loop. Our numerical experiments demonstrate that the adaptive $s$-step CG method is able to attain up to the same accuracy as classical CG while still significantly reducing the number of outer loop iterations (i.e., synchronizations) performed. In contrast, the fixed $s$-step CG method can fail to converge to this accuracy when $s>4$.

Additionally, the adaptive $s$-step CG method adjusts to use larger $s_{k}$ when an application requires a less accurate solution, automatically optimizing the tradeoff between speed and accuracy. In this way, a user can use the matrix structure and machine parameters to estimate the value $s_{\text{max}}$ that will minimize the time per iteration (e.g., by guessing, or offline auto-tuning), input the accuracy they require, and the adaptive $s$-step CG method will choose $s_{k}$ as large as it can (i.e., reducing as much synchronization as it can) while still converging to the desired level. This takes the burden of parameter selection off the user. Our adaptive $s$-step CG method thus makes significantly advances in improving reliability and usability in $s$-step Krylov subspace methods.

Our future work includes a high-performance implementation and performance experiments in order to compare classical CG, fixed $s$-step CG, and adaptive $s$-step CG on large-scale problems. Another potential improvement is to obtain a tighter bound to replace (15), which could avoid the requirement of adjusting the function $c_{k,j}$ (which represents a multiplicative factor in the bounds (15) and (16)) in Algorithm 4.

The adaptive $s$-step CG method is designed to maximize $s_{k}$ subject to a constraint on accuracy, but obviously it would also be useful to optimize $s_{k}$ to achieve the overall lowest latency/synchronization cost. This would require a better understanding of the finite precision convergence behavior of both $s$-step and classical Krylov subspace methods, which is an active area of research.

We believe that this adaptive approach based on a bound on the residual gap can be extended to other $s$-step Krylov subspace methods as well. The same approach outlined here can be used in $s$-step BICG (see, e.g., [7]), which uses the same recurrences for the solution and updated residual. We conjecture that a similar approach applies to $s$-step variants of other recursively computed residual methods (see [22]), like $s$-step BICGSTAB [7], as well as methods like $s$-step GMRES. In the case of $s$-step GMRES, it may be possible to simplify the adaptive algorithm: since the residual norm monotonically decreases, the sequence of $s_{k}$ values should monotonically increase. Such a sequence of $s_{k}$ values was used in the variable $s$-step GMRES method of Imberti and Erhel [28].

As a final comment, we note that in our analysis in section 4, we used equation  (15) to determine a bound on $\kappa(\widehat{\mathcal{Y}}_{k,s_{k}})$ that guarantees accuracy on the level $\varepsilon^{*}$. We could have just as easily used equation  (15) to develop a mixed precision $s$-step CG method where $s$ is fixed, but the precision used in computing in outer loop $k$ varies. In other words, rearranging the inequality (15) we see that, given that we’ve computed an $s$-step basis with condition number $\kappa(\widehat{\mathcal{Y}}_{k,s})$, if we want to achieve accuracy $\varepsilon^{*}$, we must have that

[TABLE]

Investigating the potential for mixed precision $s$-step Krylov subspace methods based on this analysis remains future work.