An Adaptive $s$-step Conjugate Gradient Algorithm with Dynamic Basis Updating

TL;DR

This paper introduces an improved adaptive s-step conjugate gradient algorithm that dynamically updates eigenvalue estimates to enhance convergence and reduce synchronization costs in solving large sparse linear systems.

Contribution

It develops a method to iteratively estimate Ritz values for better basis conditioning and automatic parameter tuning, advancing the efficiency of s-step CG algorithms.

Findings

Improved convergence behavior in numerical experiments.

Reduced synchronization costs compared to previous methods.

Enhanced basis conditioning through dynamic parameter updates.

Abstract

The adaptive $s$-step CG algorithm is a solver for sparse, symmetric positive definite linear systems designed to reduce the synchronization cost per iteration while still achieving a user-specified accuracy requirement. In this work, we improve the adaptive $s$-step conjugate gradient algorithm by use of iteratively updated estimates of the largest and smallest Ritz values, which give approximations of the largest and smallest eigenvalues of $A$, using a technique due to Meurant and Tich{\' y} [G. Meurant and P. Tich{\' y}, Numer. Algs. (2018), pp.~1--32]. The Ritz value estimates are used to dynamically update parameters for constructing Newton or Chebyshev polynomials so that the conditioning of the $s$-step bases can be continuously improved throughout the iterations. These estimates are also used to automatically set a variable related to the ratio of the sizes of the error and residual, which was previously treated as an input parameter. We show through numerical experiments that in many cases the new algorithm improves upon the previous adaptive $s$-step approach both in terms of numerical behavior and reduction in number of synchronizations.

Tables

Table 1. Table 1: Test matrix properties

Matrix	$N$	nnz	$\Vert A\Vert$	$\kappa \left(A\right)$
494bus	494	1666	2.00	$7.90\cdot {10}^4$
bcsstk09	1083	18437	1.98	$1.04\cdot {10}^4$
gr_30_30	900	7744	1.49	$1.95\cdot {10}^2$
nos6	675	3255	2.00	$3.49\cdot {10}^6$
mhdb416	416	2312	1.52	$3.02\cdot {10}^3$
nos1	237	1017	2.00	$3.94\cdot {10}^6$

Equations

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

$$T_i=\left[\begin{array}{ccccc}\frac{1}{{\alpha }_0}& \frac{\sqrt{{\beta }_0}}{{\alpha }_0}& & & \\ \frac{\sqrt{{\beta }_0}}{{\alpha }_0}& \frac{1}{{\alpha }_1}+\frac{{\beta }_0}{{\alpha }_0}& \frac{\sqrt{{\beta }_1}}{{\alpha }_1}& & \\ & \frac{\sqrt{{\beta }_1}}{{\alpha }_1}& \frac{1}{{\alpha }_2}+\frac{{\beta }_1}{{\alpha }_1}& \ddots & \\ & & \ddots & \ddots & \frac{\sqrt{{\beta }_{i-2}}}{{\alpha }_{i-2}}\\ & & & \frac{\sqrt{{\beta }_{i-2}}}{{\alpha }_{i-2}}& \frac{1}{{\alpha }_{i-1}}+\frac{{\beta }_{i-2}}{{\alpha }_{i-2}}\end{array}\right];$$

$$p_{sk+\ell },r_{sk+\ell }$$

$$x_{sk+\ell }-x_{sk}$$

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

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

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

$${\rho }_{\ell }^{(k)}(z)$$

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

$${\mathcal{B}}_{k,s}=\left[\begin{array}{cc}{B}_{s+1}^{(k)}& 0\\ 0& B_s^{(k)}\end{array}\right],\text{with}{B}_i^{(k)}\equiv \left[\begin{array}{ccccc}{\theta }_0^{(k)}& {\mu }_0^{(k)}& & & \\ {\gamma }_0^{(k)}& {\theta }_1^{(k)}& \ddots & & \\ & {\gamma }_1^{(k)}& \ddots & {\mu }_{i-3}^{(k)}& \\ & & \ddots & {\theta }_{i-2}^{(k)}& 0\\ & & & {\gamma }_{i-2}^{(k)}& 0\end{array}\right].$$

$$[x_{sk+\ell }-x_{sk},r_{sk+\ell },p_{sk+\ell }]={\mathcal{Y}}_{k,s}[x_{k,\ell }^{\prime },r_{k,\ell }^{\prime },p_{k,\ell }^{\prime }],$$

$$x_{k,j+1}^{\prime }$$

$$r_{k,j+1}^{\prime }$$

$$p_{k,j+1}^{\prime }$$

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

$$=\frac{r_{k,j}^{\prime T}{G}_{k,s}{r}_{k,j}^{\prime }}{p_{k,j}^{\prime T}{G}_{k,s}{\mathcal{B}}_{k,s}{p}_{k,j}^{\prime }},$$

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

$$=\frac{r_{k,j+1}^{\prime T}{G}_{k,s}{r}_{k,j+1}^{\prime }}{r_{k,j}^{\prime T}{G}_{k,s}{r}_{k,j}^{\prime }},$$

$$G_{k,s}={\mathcal{Y}}_{k,s}^T{\mathcal{Y}}_{k,s}$$

$${\theta }_0$$

$${\theta }_{\ell }$$

$${\mu }_{\ell }$$

$${\gamma }_0$$

$$\Vert {\delta }_{m+j+1}-{\delta }_m\Vert \le c\epsilon \left(\underset{0\le \ell \le j+1}{\max }\Vert {\hat{r}}_{m+\ell }\Vert \right),$$

$$\Vert {\delta }_{m+j+1}-{\delta }_m\Vert \le c_k\epsilon \kappa ({\hat{\mathcal{Y}}}_{k,s})\left(\underset{0\le \ell \le j+1}{\max }\Vert {\hat{r}}_{m+\ell }\Vert \right)+\epsilon \Vert A\Vert \Vert x\Vert ,$$

$$c_{k}=2s\big{(}2(3+N_{A})\nu t+(6+8t)\tau_{k}+2t^{3}+3\big{)}\kappa(A),$$

$$\kappa ({\hat{\mathcal{Y}}}_{k,s})\le \frac{{\epsilon }^{\ast }}{c_k\epsilon \Vert {\hat{r}}_{m+\ell }\Vert },\text{for}0\le \ell \le s.$$

$$\kappa ({\hat{\mathcal{Y}}}_{k,{\tilde{s}}_k})\le \frac{{\epsilon }^{\ast }}{c_k\epsilon \Vert {\hat{r}}_m\Vert }$$

$$\Vert r_{m+j+1}\Vert =\sqrt{r_{k,j+1}^{\prime T}{G}_{k,{\tilde{s}}_k}{r}_{k,j+1}^{\prime }}.$$

$$L_i^T\equiv \left[\begin{array}{cccc}{\zeta }_0& {\eta }_0& & \\ & \ddots & \ddots & \\ & & \ddots & {\eta }_{i-2}\\ & & & {\zeta }_{i-1}\end{array}\right]=\left[\begin{array}{cccc}\frac{1}{\sqrt{{\alpha }_0}}& \sqrt{\frac{{\beta }_0}{{\alpha }_0}}& & \\ & \ddots & \ddots & \\ & & \ddots & \sqrt{\frac{{\beta }_{i-2}}{{\alpha }_{i-2}}}\\ & & & \frac{1}{\sqrt{{\alpha }_{i-1}}}\end{array}\right],$$

$${\tilde{\lambda }}_{\text{max}}={\lambda }_{\text{max}}(T_i)=\Vert L_i{\Vert }^2,{\tilde{\lambda }}_{\text{min}}={\lambda }_{\text{min}}(T_i)=\Vert L_i^{-1}{\Vert }^{-2}.$$

$$\Vert A\Vert \Vert {\hat{x}}_{m+j+1}-x\Vert$$

$$\le \kappa (A)\Vert A({\hat{x}}_{m+j+1}-x)\Vert =\kappa (A)\Vert {\hat{r}}_{m+j+1}\Vert +O(\epsilon ),$$

$${\xi }_{m+j+1}\equiv \frac{\Vert A\Vert \Vert {\hat{x}}_{m+j+1}-x\Vert }{\Vert {\hat{r}}_{m+j+1}\Vert },$$

$$\Vert {\hat{x}}_{m+j+1}-x\Vert \le \frac{\Vert {\hat{r}}_{m+j+1}\Vert }{{\mu }^{1/2}}{\psi }_{m+j+1}^{1/2},$$

Full text

Keywords: conjugate gradient, iterative methods, high-performance computing

MSC 2010: 65F10, 65F50, 65Y05, 65Y20

1 Introduction

In this work we focus on the problem of solving linear systems $Ax=b$, where $A\in\mathbb{R}^{N\times N}$ is symmetric positive definite (SPD). When $A$ is large and sparse, the iterative conjugate gradient method (CG), which is a Krylov subspace method, is commonly-used as a solver. Given an initial approximate solution $x_{0}$ with initial residual $r_{0}=b-Ax_{0}$, Krylov subspace methods construct a sequence of nested Krylov subspaces $\mathcal{K}_{1}(A,r_{0})\subset\mathcal{K}_{2}(A,r_{0})\subset\cdots\subset\mathcal{K}_{i}(A,r_{0})$ where

[TABLE]

In CG, the approximate solution $x_{i}\in x_{0}+\mathcal{K}_{i}(A,r_{0})$ is constructed according to the orthogonality constraint $r_{i}=b-Ax_{i}\perp\mathcal{K}_{i}(A,r_{0})$, which is equivalent to selecting the vector $x_{i}\in x_{0}+\mathcal{K}_{i}(A,r_{0})$ that minimizes the $A$-norm (energy norm) of the error, i.e., $\|x-x_{i}\|_{A}=\sqrt{(x-x_{i})^{T}A(x-x_{i})}$.

Perhaps the most well-known algorithm for CG is due to Hestenes and Stiefel [25], which uses three coupled 2-term recurrences for recursively updating the approximate solution $x_{i}$, the residual $r_{i}$, and $A$-conjugate search direction vectors $p_{i}$. We refer to this particular algorithm, displayed in Algorithm 1, as “HSCG” in this work.

The CG method is closely related to the Lanczos tridiagonalization method, which iteratively constructs a symmetric tridiagonal matrix $T_{i}\in\mathbb{R}^{i\times i}$ such that $T_{i}=V_{i}^{T}AV_{i}$ where $V_{i}$ is an orthonormal basis for the Krylov subspace $\mathcal{K}_{i}(A,r_{0})$. The $i$ eigenvalues of $T_{i}$, called Ritz values, provide estimates of the eigenvalues of $A$. The matrix $T_{i}$ can be written in terms of the coefficients $\alpha_{i}$ and $\beta_{i}$ computed in lines 5 and 9 of Algorithm 1 via

[TABLE]

see, e.g., [36, Sec. 6.7.3].

In the setting of large-scale problems on parallel machines, the performance of HSCG is limited by communication, i.e., data movement, due to the sparse matrix vector product (SpMV) and inner products in each iteration, both of which have low computational intensity; see, e.g., [16]. This has led to the development of a number of algorithmic variants of CG which aim to reduce the communication and/or synchronization cost over a fixed number of steps. One such variant is called $s$-step CG (also called communication-avoiding CG; see, e.g., [27], [8] and the historical references therein). The $s$-step CG algorithm works by computing $O(s)$ new basis vectors for the Krylov subspace at a time and then computing a block inner product between these computed basis vectors. The former can be accomplished with $O(1)$ messages between parallel processors (usually assumed to involve communication only between neighboring processors) under some constraints on the sparsity structure and parallel partition of the matrix $A$; see [14] for details. The latter can be accomplished with a single global synchronization point. Over $s$ iterations, this approach can reduce the number of synchronizations from $O(s)$ to $O(1)$. We elaborate on the mathematics behind this approach in Section 2.

In practice, however, we are not concerned only with the number of synchronizations over a fixed number of steps, but rather the total number of synchronizations required to achieved the prescribed convergence criterion. It is well-known that $s$-step CG (and $s$-step variants of other Krylov subspace methods) can exacerbate the delay of convergence and decrease in attainable accuracy that are characteristic of finite precision HSCG. This behavior generally grows worse for increasing $s$. In the extreme case, this can lead to a situation where the prescribed accuracy is no longer attainable, making the $s$-step approach inapplicable.

In $s$-step CG, it has been shown that the loss of attainable accuracy can be bounded in terms of the condition numbers of the computed $O(s)$-dimensional bases for the Krylov subspaces generated at the beginning of each block of $s$ iterations [6]. This rounding error analysis led to the insight that if we wish to achieve a certain accuracy, then the condition numbers of the bases must be controlled to be less than a certain quantity inversely proportional to the largest residual norm within the current block of $s$ iterations. This naturally suggests a variable $s$-step approach; when the size of the residual is still large (at the beginning of the iterations), $s$ should be small so as to keep the bases well-conditioned, but as the method converges and the size of the residual decreases, $s$ can be gradually increased without detriment to the maximum attainable accuracy.

This inspired the adaptive $s$-step CG algorithm [9], in which the value of $s$ is automatically adjusted according to the basis condition number and the user-specified accuracy requirement. It was shown in [9] that this adaptive approach can provide improved reliability in terms of numerical behavior while still reducing the overall number of synchronizations. Given the importance of the condition of the bases, one particular gap (which applies to most $s$-step CG algorithms in the literature) is the reliance on user-supplied parameters for generating the polynomial bases for the $O(s)$-dimensional Krylov subspaces. If this information is not known a priori, either basis parameters must be computed before execution of the $s$-step algorithm (by precomputing information about the spectrum of $A$ or running a number of iterations of HSCG). A simple monomial basis can be used as a default, it is known that the condition number of the monomial basis grows exponentially with $s$ [18]. The previous method also required the user to heuristically set a certain parameter used in determining how large $s$ can be in each outer loop iteration.

In this work, we develop an improved version of the adaptive $s$-step CG algorithm which makes use of estimates of the largest and smallest eigenvalues of $A$ (the largest and smallest Ritz values) obtained automatically and inexpensively as a byproduct of the iterations using the technique of Meurant and Tichý [33]. The eigenvalue estimates are used for two purposes: (1) to dynamically update coefficients for Newton or Chebyshev polynomials so that the conditioning of the $s$-step bases can be continuously improved throughout the iterations and (2) to automatically set the previously heuristically-chosen parameter based on information obtained during the iterations.

In Section 2, we briefly review $s$-step CG. We then outline the idea of the adaptive $s$-step CG algorithm in Section 3. In Section 4, we review the work of Meurant and Tichý [33] and use this to develop an improved adaptive $s$-step CG algorithm. Section 5 presents numerical experiments for a variety of small test problems, which demonstrate the benefits of the improved approach. We conclude and discuss further challenges in Section 6.

2 The $s$-step CG algorithm

The idea of $s$-step Krylov subspace algorithms is not new. The first known appearance of $s$-step CG in the literature was due to van Rosendale [40], although the name “$s$-step CG” was first coined later by Chronopoulos and Gear [10]. There have since been many efforts towards developing $s$-step formulations of CG and other Krylov subspace methods; for a thorough overview of related works, see [27, Table 1.1]. Much of the early work in this area was motivated by reducing the amount of I/O and/or by increasing the potential for parallelism in CG. These algorithms (and other variants of Krylov subspace methods designed to reduce communication overhead, such as pipelined algorithms [19], [20]) have recently come back into vogue for their potential to reduce data movement, both between levels of the memory hierarchy on a single processor and between processors in the parallel setting, which becomes increasingly important in efforts to scale to larger problem sizes and larger machines [15].

The $s$-step approach can be thought of as blocking the iterations into sets of $s$. For one block of $s$ iterations, one first expands the underlying Krylov subspace by $O(s)$ dimensions and subsequently performs a block orthogonalization using only a single global synchronization point. The vector updates for the block of $s$ iterations can then be performed by updating the $O(s)$ coordinates of the vectors in the generated $O(s)$-dimensional Krylov subspace. The particular details of the algorithm depend on the particular underlying Krylov subspace method, but the general concept is the same. In order to establish notation, we give a brief overview of the $s$-step CG algorithm.

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

[TABLE]

Then the CG vectors for the next $s$ iterations to be computed within this block lie in the union of the column spaces of the matrices

[TABLE]

where $\rho^{(k)}_{\ell}(z)$ is a polynomial of degree $\ell$ satisfying the three-term recurrence

[TABLE]

The need for the superscripts $(k)$ above will become clear later when we introduce the improved adaptive $s$-step CG algorithm, in which the coefficients in the recurrence (2.3) are dynamically updated between outer loop iterations. Under certain constraints on the sparsity structure and the partition of the data, the matrices (2.2) can be computed with $O(1)$ messages per processor in a parallel algorithm; see [14] for details.

We define the $s$-step basis matrix $\mathcal{Y}_{k,s}=[\mathcal{P}_{k,s},\mathcal{R}_{k,s}]$ and we define $\underline{\mathcal{Y}}_{k,s}$ to be the same as $\mathcal{Y}_{k,s}$ except with columns $s+1$ and $2s+1$ set to zero. We can then write the recurrence relation

[TABLE]

where

[TABLE]

For $\ell\in\{0,\ldots,s\}$ we can then represent the vectors $x_{sk+\ell}-x_{sk}$ , $r_{sk+\ell}$, and $p_{sk+\ell}$ by their $2s+1$ coordinates in the basis spanned by the columns of $\mathcal{Y}_{k,s}$, i.e.,

[TABLE]

and the updates to these coordinate vectors in the inner loop from $j\in\{0,\ldots,s-1\}$ become

[TABLE]

In practical applications we expect $s\ll N$, so the updates to the length-$(2s+1)$ coordinate vectors can be accomplished locally on each processor without any further communication. The coefficients $\alpha_{sk+j}$ and $\beta_{sk+j}$ can also be computed locally on each processor without communication; using (2.4) and (2.6), we have

[TABLE]

where

[TABLE]

is the $(2s+1)\times(2s+1)$ Gram matrix which is computed only once per outer loop iteration (requiring a single global synchronization) and stored locally on each processor.

In (2.5), the choice $\theta^{(k)}_{\ell}=0$, $\gamma^{(k)}_{\ell}=1$, and $\mu^{(k)}_{\ell}=0$ for all $\ell$ corresponds to the monomial basis. Without further information about the spectrum, this is the simplest choice we have, although it is known that the monomial basis quickly becomes ill-conditioned with $s$ [18], and this thus limits the $s$ we can choose. Already in early works on $s$-step CG it was observed that the conditioning of the $s$-step basis matrices plays a large role in the resulting finite precision behavior, which led many to experiment with more well-conditioned polynomial bases such as Newton or Chebyshev bases; see, e.g, [26, 12, 29, 1, 17, 13].

Estimates of the maximum and minimum eigenvalues of $A$, $\tilde{\lambda}_{\text{max}}$ and $\tilde{\lambda}_{\text{min}}$, can be used to construct either Newton or Chebyshev polynomials, which will in general result in a better-conditioned basis and thus allow the use of larger $s$ values without loss of accuracy. We note that the idea of adaptively improving basis conditioning using Ritz values is not new. A method for generating parameters for Newton and Chebyshev polynomial bases for Krylov subspaces based on Ritz values is described in [34]; see also the works [31] and [4]. The modern software package Trilinos [24] currently includes the ability to automatically generate Newton or Chebyshev polynomials in their implementations of (fixed) $s$-step Krylov subspace algorithms. The approach they use is to first perform a fixed number of iterations of the classical algorithm (with $s=1$) and then use the Ritz value estimates (eigenvalues of $T_{i}$) to generate basis parameters.

We note that, in contrast with the use of estimates of $\tilde{\lambda}_{\text{max}}$ and $\tilde{\lambda}_{\text{min}}$ within the Chebyshev semi-iterative method, for the purposes of generating a Krylov subspace basis, $\tilde{\lambda}_{\text{max}}$ and $\tilde{\lambda}_{\text{min}}$ do not need to be particularly accurate in order to provide an improvement over the monomial basis; see, e.g., [34, pp. 12], where it is stated (referring to the Chebyshev basis) that “it is not important that the ellipse be determined to high accuracy.”

Given $\tilde{\lambda}_{\text{max}}$ and $\tilde{\lambda}_{\text{min}}$, the Newton basis parameters can be taken as (for now dropping the superscripts) $\mu_{\ell}=0$, $\gamma_{\ell}=1$, and

[TABLE]

which corresponds to a Leja ordering of the points on the real line between $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$. The Leja ordering is known to improve the accuracy of operations on polynomials; see, e.g., [1], [5], and [35]. Parameters for a simplified version of the Chebyshev basis can be chosen as (see [29, Section 4.4])

[TABLE]

2.1 The $s$-step CG algorithm in finite precision

It is well-known that finite precision roundoff errors can cause a delay of convergence and a decrease of attainable accuracy in CG algorithms; see, e.g., [32], [30] and references therein. In this work, we focus on the maximum attainable accuracy, as achieving a certain prescribed accuracy is the goal of the adaptive $s$-step approaches.

The size of the true residual is often used as a computable measure of accuracy. The mechanism by which accuracy is lost is the deviation of the recursively updated residual $\widehat{r}_{i}$ and the true residual $b-A\widehat{x}_{i}$, where $\widehat{r}_{i}$ and $\widehat{x}_{i}$ denote the quantities computed in finite precision (in general, we will now use hats to denote quantities computed in finite precision). Writing $b-A\widehat{x}_{i}=(b-A\widehat{x}_{i}-\widehat{r}_{i})+\widehat{r}_{i}$, it is clear that as the size of the recursively updated residual becomes very small, the upper bound on the size of $b-A\widehat{x}_{i}$ depends on the size of the residual gap $\delta_{i}\equiv b-A\widehat{x}_{i}-\widehat{r}_{i}$.

There is a large literature on analyses of maximum attainable accuracy in HSCG and related CG algorithms, including the works of Greenbaum [22], van der Vorst and Ye [39], Sleijpen and van der Vorst [38], and Gutknecht and Strakoš [23]. Modifying slightly the bound on the residual gap in HSCG derived by Sleijpen and van der Vorst [38, Eqn. (6)], we can bound the growth of the residual gap starting from some iteration $m\equiv sk$ to iteration $m+j+1$, $j\in\{0,\ldots,s-1\}$, in HSCG by

[TABLE]

where $c=2sN_{A}\nu\kappa(A)$, $\kappa(A)=\|A^{-1}\|\|A\|$ denotes the condition number, $N_{A}$ is the maximum number of nonzeros per row in $A$, $\nu=\||A|\|/\|A\|$, $\|\cdot\|$ denotes the 2-norm, and $\varepsilon$ is the machine unit roundoff.

In [9, Eqn. (16)] it is shown that the growth of the residual gap within one outer loop iteration $k$ (which begins at global iteration $m$) of $s$-step CG can be bounded by

[TABLE]

for $j\in\{0,\ldots,s-1\}$. The constant $c_{k}$ comes from the rounding error analysis and can be written as

[TABLE]

where $t=\sqrt{2s+1}$ and $\tau_{k}=\||\mathcal{B}_{k,s}|\|/\|A\|$. The most notable difference between (2.9) and (2.10) is the appearance of the term $\kappa(\widehat{\mathcal{Y}}_{k,s})=\|\widehat{\mathcal{Y}}_{k,s}^{+}\|\|\widehat{\mathcal{Y}}_{k,s}\|$, where $\widehat{\mathcal{Y}}_{k,s}^{+}$ denotes the Moore-Penrose pseudoinverse. In other words, the local roundoff errors made in $s$-step CG are amplified by the condition numbers of the computed $s$-step basis matrices. This theoretically confirms observations regarding the effect of the condition numbers of the $s$-step bases on the numerical behavior of $s$-step CG compared to HSCG; an ill-conditioned $s$-step basis can cause an increase in the residual gap and thus can decrease the attainable accuracy.

3 The adaptive $s$-step CG algorithm

Rearranging (2.10), it can be shown that if the application requires a relative residual norm of $\varepsilon^{*}$, the condition number of the basis matrix $\widehat{\mathcal{Y}}_{k,s}$ must satisfy

[TABLE]

This naturally suggests that $s$ should be allowed to vary in each outer loop $k$; when the residuals are large, $\kappa(\widehat{\mathcal{Y}}_{k,s})$ and thus $s$ must be small, but as the residual is reduced, the condition number of the basis and thus $s$ can be larger without detriment to the attainable accuracy111 We note that this insight is similar to that behind the development of the so-called “inexact Krylov subspace methods”; see, e.g., [37] as well as the technical report [2] which was later published as [3]. . We therefore introduce the subscripted quantity $s_{k}$ to denote the number of inner iterations in outer loop $k$. Quantities which depend on $s_{k}$ will now also have a subscript $k$, e.g., $\widehat{\mathcal{Y}}_{k,s_{k}}$, $\mathcal{B}_{k,s_{k}}$, $t_{k}=\sqrt{2s_{k}+1}$, and the expression for $c_{k}$ now containing the relevant quantities dependent on the value of $s_{k}$ in outer loop $k$.

This idea led to the adaptive $s$-step CG algorithm published in [9], displayed in Algorithm 2. It is shown in [9] that, assuming that the algorithm converges, as long as (3.1) is satisfied, then adaptive $s$-step CG can attain a solution to the same level of accuracy as HSCG. We now give a brief description of the algorithm.

In CG there is no guarantee that the residual norms are monotonically decreasing (in exact arithmetic, the method rather minimizes the $A$-norm of the error). Therefore, for an outer loop beginning at iteration $m=\sum_{\ell=0}^{k-1}s_{\ell}$, we use $\bar{s}_{k}$ to denote our initial guess for $s_{k}$ and construct the basis $\widehat{\mathcal{Y}}$ only based on the current residual $\widehat{r}_{m}$. We then 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}}$. From this, we find the largest value $\tilde{s}_{k}\leq\bar{s}_{k}$ such that condition (3.1) holds, i.e., $\tilde{s}_{k}$ is the maximum value in $\{1,\ldots,\bar{s}_{k}\}$ such that

[TABLE]

holds. The values $\kappa(\widehat{\mathcal{Y}}_{k,i})$ can be estimated using the square roots of the leading principle submatrices of the constructed $\widehat{G}_{k,\bar{s}_{k}}$, since since $\kappa(\widehat{\mathcal{Y}}_{k,\ell})\approx\sqrt{\kappa(\widehat{G}_{k,\ell})}$. The quantities $\sqrt{\kappa(\widehat{G}_{k,\ell})}$ are inexpensive to compute; this involves $O(s^{3})$ floating point operations to compute the eigenvalues of $\widehat{G}_{k,\bar{s}_{k}}$ and no additional data movement as $\widehat{G}_{k,\bar{s}_{k}}$ is stored locally on each processor. We then check if the condition (3.2) is violated within each inner loop, which can occur if we encounter a large intermediate residual norm. We use $s_{k}\leq\tilde{s}_{k}$ to denote the actual number of inner loop iterations which occurred. We note that the residual norms can be estimated cheaply (without communication) within the inner loop since (in exact arithmetic)

[TABLE]

The algorithm also requires the user to input some $\bar{s}_{0}$ as an initial value, a value $\sigma$ which is the maximum value for $\bar{s}_{k}$ (which could be determined by offline auto-tuning and should be based on the matrix nonzero structure, machine parameters, and matrix partition), and a value $f$, which is the maximum $\bar{s}_{k}$ can be allowed to grow in each iteration (e.g., $f$ could be made small to reduce the amount of wasted flops if $s_{k}$ is much smaller than $\bar{s}_{k}$ in each outer loop, or $f$ could be made larger to maximize the potential $s_{k}$ in each outer loop).

In the experiments in [9], it was found that the value of $c_{k}$ in (2.11) is often a large overestimate, resulting in smaller $s_{k}$ values than necessary to achieve the desired accuracy which results in more outer loop iterations than necessary. It was found that in most cases, taking $c_{k}=1$ in (3.1) worked well, although there was no theoretical justification to support this. In the following section, we describe a way to adaptively and automatically set this parameter based on existing quantities obtained during the iterations.

We briefly comment on related work in the area of using a variable $s$ value in $s$-step Krylov subspace algorithms. Also motivated by improving numerical behavior, Imberti and Erhel used a variable $s$ value in their $s$-step GMRES algorithm [28], although their approach requires the user to prescribe a priori the sequence of $s_{k}$ values. A variable $s$ value was also used within a $s$-step BICGSTAB algorithm used as the coarse grid solve routine within a geometric multigrid method [41]. In [41] this approach was termed a “telescoping $s$”, in which the value of $s$ starts small and is allowed to grow as the outer loops proceed. This was done for performance reasons rather than numerical ones; when the coarse grid problem is easy (converges in a few iterations), we do not waste effort computing a larger $s$-step basis than need be.

4 An improved adaptive $s$-step CG algorithm

As stated, the adaptive $s$-step CG algorithm described in [9] requires the user to supply parameters for the polynomial recurrence (2.3) used in constructing the $s$-step basis matrices as well as the function $c_{k}$ used in determining when and if to break from the inner loop. In this section we present an improved adaptive $s$-step CG algorithm which uses Ritz value estimates computed via the results of Meurant and Tichý [33], which we now briefly summarize. The $\alpha$ and $\beta$ coefficients computed during the CG iterations can be composed to form the Cholesky factor $L_{i}^{T}$ of the Lanczos tridiagonal matrix (see (1.1)) $T_{i}=L_{i}L_{i}^{T}$, where

[TABLE]

with $i$ now denoting the global iteration index (the total number of inner loop iterations). It is expected that the eigenvalues of $T_{i}$, i.e., the Ritz values, give decent approximations for the extremal eigenvalues of $A$. The extremal eigenvalues of $A$ can be estimated via the relations

[TABLE]

As described in [33], the norms $\|L_{i}\|^{2}$ and $\|L_{i}^{-1}\|^{-2}$ can computed using incremental norm estimation of the matrices $L_{i}^{T}$ and $L_{i}^{-T}$ without the need to explicitly construct these matrices; this requires insignificant extra work and no communication. The algorithms for estimating $\|L_{i}\|^{2}$ and $\|L_{i}^{-1}\|^{-2}$ are displayed as Algorithms 3 and 4, respectively, which appear as Algorithms 4 and 5 in [33, Sec. 5] with slightly modified notation and indexing.

Thus in each inner loop iteration, we can update estimates of the extremal eigenvalues of $A$, $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$, using Algorithms 3 and 4. We note that in each inner loop iteration, updating $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ requires executing only a single additional for-loop in Algorithms 3 and 4 and storing only a small number of scalar quantities from the previous iteration.

The estimates $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ are then used in two ways to improve the adaptive $s$-step CG algorithm. The first improvement comes from improving the quality of the computed Krylov subspace bases. The computation begins by using a monomial basis to construct the $\bar{s}_{k}$-step polynomial bases. Once at least two iterations have finished (i.e., $\tilde{\lambda}_{\text{min}}\neq\tilde{\lambda}_{\text{max}}$), the estimates $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ can be used in the subsequent outer loop $k$ to estimate the basis parameters $\theta_{i}^{(k)}$, $\gamma_{i}^{(k)}$, and $\mu_{i}^{(k)}$ for either the Newton basis (2.7) or the Chebyshev basis (2.8). The process continues, with the latest updated estimates $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ being used to generate the basis parameters after every outer loop, improving the quality of the bases as the iterations proceed.

The second improvement comes from using the estimates $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ to eliminate the heuristic choice of the parameter $c_{k}$. As mentioned, the adaptive $s$-step CG algorithm in [9] required that the quantity $c_{k}$ in the condition (3.1) be set by the user, as using the full value from (2.11) was in most cases too restrictive. This is because the $\kappa(A)$ term arises in the quantity $c_{k}$ as defined in (2.11) as a result of the effort to bound the residual gap solely in terms of the size of the residuals, for example,

[TABLE]

which can give a rather pessimistic bound in practice. We note that a similar observation was made by Sleijpen and van der Vorst [38], who argued that the $\kappa(A)$ term will only appear in ‘unusual’ cases. The numerical experiments in [9] found that using $c_{k}=1$ (i.e., ignoring the $\kappa(A)$ term) worked well in most (but not all) cases.

We make further use of the results of Meurant and Tichý [33] and develop a way to adaptively set $c_{k}$ that removes the burden of setting this parameter heuristically. Instead of the potentially loose upper bound (4.3), we will define $\xi_{m+j+1}$ to be the exact quantity defined by $\|A\|\|\widehat{x}_{m+j+1}-x\|=\xi_{m+j+1}\|\widehat{r}_{m+j+1}\|$, i.e.,

[TABLE]

where we will have $1\leq\xi_{m+j+1}\leq\kappa(A)$. It is shown in [33, Section 3] that the $A$-norm of the error can be upper bounded as

[TABLE]

where $0<\mu\leq\lambda_{\text{min}}$ and the quantity $\psi_{m+j+1}\equiv\|\widehat{r}_{m+j+1}\|^{2}/\|\widehat{p}_{m+j+1}\|^{2}$ can be incrementally updated in each iteration by

[TABLE]

Thus in each iteration we take $\mu=\tilde{\lambda}_{\text{min}}$ as computed in (4.2) and approximately bound $\xi_{m+j+1}$ in (4.4) via

[TABLE]

Thus in each iteration we set $c_{m+j+1}=\tilde{\xi}_{m+j+1}$ (note the subscript notation change which indicates that this quantity can now change in each inner iteration). Updating $\tilde{\xi}_{m+j+1}$ in (4.6) requires only the additional scalar operation in (4.5) to update $\psi_{m+j+1}$. Our experiments in Section 5 confirm that $\tilde{\xi}_{m+j+1}$ can be far below $\kappa(A)$.

We make a further small improvement to the adaptive $s$-step CG algorithm in [9]. If we are at the end of some inner iteration $j<\tilde{s}_{k}-1$, notice that by (2.1), the iterates that will be updated in the next iteration $j+1$ depend on only a subset of the basis vectors. The previous adaptive $s$-step CG algorithm (Algorithm 2) breaks from the inner loop if in some inner iteration $j$,

[TABLE]

in other words, if the residual norm $\|r_{m+j+1}\|$ is so large that (3.1) will not hold in some future inner iteration (not only just the next one). This is overly pessimistic and can cause the algorithm to quit the inner loop iterations unnecessarily early. We can easily and inexpensively modify this approach as follows.

In each inner loop iteration, we keep track of the maximum residual norm that we have encountered so far in outer loop $k$, and store this as a variable $\phi$. At the end of inner loop iteration $j<\tilde{s}_{k}-1$, we check whether

[TABLE]

and if so, we break from the current inner loop. If not, we continue with the next inner loop iteration $j+1$. As previously described, it is easy and inexpensive (involving no communication) to estimate all the condition numbers $\kappa(\widehat{\mathcal{Y}}_{k,\ell+1})$ for $\ell\in\{1,\ldots,\tilde{s}_{k}-1\}$ by computing the square roots of the condition numbers of the appropriate leading principal submatrices of $\widehat{G}_{k,\tilde{s}_{k}}$. The resulting improved adaptive $s$-step CG algorithm is displayed as Algorithm 5. The changes compared with Algorithm 2 are highlighted in red.

5 Numerical experiments

In this section we present experiments run in MATLAB (version R2017a) to compare the numerical behavior of HSCG, fixed $s$-step CG, adaptive $s$-step CG with a monomial basis, and improved adaptive $s$-step CG with dynamically updated Newton and Chebyshev bases for small SPD matrices from the SuiteSparse collection [11]. For each matrix, we test the $s$ values $5,10,$ and $15$ (which are the maximum allowable $s_{k}$ values $\sigma$ in the adaptive algorithms). We test two different values of $\varepsilon^{*}$; the first being the relative true residual $2$-norm attainable by HSCG as determined experimentally, and the second being $\varepsilon^{*}=10^{-6}$. The adaptive algorithms all use $f=\sigma$ to allow for the largest possible $s_{k}$ values. For all experiments in this section, we use two-sided diagonal preconditioning, where the resulting preconditioned matrix is $D^{-1/2}AD^{-1/2}$ where $D$ is a diagonal matrix of the largest entries in each row of $A$. The (unpreconditioned) right hand side $b$ is set with entries $1/\sqrt{N}$ and the initial guess is $x_{0}=0$. We use double precision in all tests ($\varepsilon\approx 2^{-53}$). We discuss the relative performance of the (fixed and adaptive) $s$-step algorithms in terms of the number of outer loop iterations, which can be seen as a proxy for the number of global synchronizations. We stress, however, that the experiments here are on very small matrices, and thus the aim is solely to demonstrate the numerical behavior rather than evaluate potential parallel performance improvements.

5.1 The benefit of adaptively setting $c_{m+j+1}$

We first demonstrate the benefits of our approach for dynamically setting the constant $c_{m+j+1}$. In Figure 1, we plot the convergence (in terms of relative residual $2$-norm) of HSCG, fixed $s$-step CG, and the improved adaptive $s$-step CG algorithm (Algorithm 5) for the matrix nos1 from SuiteSparse [11] using various values of $c_{m+j+1}$, with the specified tolerance $\varepsilon^{*}=10^{-6}$ and $s=10$ ($\sigma=10$ in the adaptive algorithm). Properties of the nos1 matrix can be found in Table 1. For this problem, HSCG takes 510 iterations to converge. Fixed $s$-step CG suffers from significantly delayed convergence, requiring a total of 7134 iterations (corresponding to 714 outer loop iterations) to converge. We note that the old adaptive $s$-step CG (Algorithm 2) in this case chooses $s_{k}=\sigma$ in each outer loop iteration and thus also requires $714$ outer iterations to converge. Thus for these algorithms we would expect performance to be worse than HSCG.

For the improved adaptive $s$-step CG algorithm, we show convergence for 3 different choices of $c_{m+j+1}$. The plots on the left show convergence of the relative true residual and the plots on the right show the $s_{k}$ values used throughout the iterations. The top row of plots uses $c_{m+j+1}=\tilde{\lambda}_{\text{max}}/\tilde{\lambda}_{\text{min}}$ where $\tilde{\lambda}_{\text{max}}$ and $\tilde{\lambda}_{\text{min}}$ are updated in every iteration, the middle row uses a constant $c_{m+j+1}=1$, and the bottom row uses the new adaptive approach which automatically sets $c_{m+j+1}=\tilde{\xi}_{m+j+1}$; see (4.6).

In the top row of plots using $c_{m+j+1}=\tilde{\lambda}_{\text{max}}/\tilde{\lambda}_{\text{min}}$, the convergence of the improved adaptive $s$-step CG algorithm is relatively close to that of the HSCG algorithm for both Newton and Chebyshev bases. However, this comes at the cost of the use of small $s_{k}$ values, with $s_{k}=1$ in many outer loop iterations (see the top right plot). Here the improved adaptive $s$-step CG algorithm requires 291 and 255 outer iterations with the Newton and Chebyshev bases, respectively.

It is clear from the middle plots that blindly using $c_{m+j+1}=1$ as previously suggested in [9] is not always sufficient in practice (note the stretched $x$-axis in the left-hand plot for this experiment). Using $c_{m+j+1}=1$ results in larger $s_{k}$ values (with the maximum $s_{k}=10$ occurring in most outer loop iterations), but this results in a significant convergence delay due to the use of ill-conditioned bases. This results in a greater total number of outer loop iterations required for convergence, with the implementations with Newton and Chebyshev bases requiring 646 and 357 outer loop iterations, respectively. Note that for the Newton basis case, this is more “global synchronizations” than in HSCG!

Finally, in the bottom row, we see the advantage of our approach for automatically setting $c_{m+j+1}$ based on information gained from the iterations. Here the total number of iterations required for convergence is greater than when we use $c_{m+j+1}=\tilde{\lambda}_{\text{max}}/\tilde{\lambda}_{\text{min}}$, but still much smaller than the naive choice $c_{m+j+1}=1$. But notice that the $s_{k}$ values used are much closer to the $c_{m+j+1}=1$ case, with the maximum value $s_{k}=10$ occurring in the majority of iterations. This results in fewer overall outer loop iterations; here the implementations with the Newton and Chebyshev bases required only 134 and 187 outer loop iterations, respectively. For reference, in Figure 2, we show the value of $c_{m+j+1}=\tilde{\xi}_{m+j+1}$ computed in each iteration. The dashed horizontal lines mark 1 (the minimum value of $\xi_{m+j+1}$) and $\kappa(A)$ (the maximum value of $\xi_{m+j+1}$). We see that the computed $\tilde{\xi}_{m+j+1}$ is actually somewhere between $1$ and $\kappa(A)$; by using this tighter, computable bound, we effectively balance the goals of minimizing the number of iterations required for convergence and maximizing the $s_{k}$ values used.

5.2 Test problems from SuiteSparse

We now present a few select test problems which demonstrate the behavior of the improved adaptive $s$-step CG algorithm. For all problems, the experimental setup is as described at the beginning of Section 5. The matrices used in the experiments are shown in Table 1, where the norm and condition number are those of the diagonally preconditioned system matrices. For space purposes, we only include plots for two of the test matrices, 494bus and bcsstk09. However, for each test matrix, we include a table which shows the number of iterations for HSCG and the number of outer loop iterations (global synchronizations) for fixed $s$-step CG, adaptive $s$-step CG, and improved adaptive $s$-step CG with both Newton and Chebyshev bases, for the $s$ (or $\sigma$) values 5, 10, and 15, and for two different choices of $\varepsilon^{*}$. For the $s$-step variants, the first number gives the number of outer loop iterations and the number inside the parentheses gives the total number of (inner loop) iterations. A dash alone in the table indicates that the algorithm diverged. A dash along with a number in square brackets indicates that the algorithm did eventually converge (i.e., the residual norm stagnated), but to a relative residual norm less than $\varepsilon^{*}$; the number in square brackets gives the final attained relative residual $2$-norm. In the plots, markers are used to denote the outer iterations in the (fixed and adaptive) $s$-step algorithms, which are the points at which global synchronization would occur.

We show results for the matrix 494bus in Figure LABEL:fig.494busn_0 and the corresponding Table LABEL:tab.494busn0. This problem represents the ideal case for the improved adaptive $s$-step algorithm and nicely highlights the benefits of the improved approach. When $\varepsilon^{*}=2.2\cdot 10^{-10}$, the improved algorithm outperforms the fixed $s$-step algorithm and the old adaptive $s$-step algorithm in all cases, even for small $\sigma$. When $\varepsilon^{*}=10^{-6}$, the improved approach takes about the same number of inner and outer loop iterations as $s$-step CG and old adaptive $s$-step CG for $s,\sigma=5$, but clearly outperforms both algorithms for higher $s,\sigma$ values.

Using the improved bases generated using $\tilde{\lambda}_{\text{min}}$ and $\tilde{\lambda}_{\text{max}}$ has a clear benefit; in all cases, the total number of iterations required by the improved adaptive approach is about the same as for HSCG. In contrast, for $\varepsilon^{*}=2.2\cdot 10^{-10}$, the fixed $s$-step and old adaptive approaches require more than 200 more iterations to converge compared to HSCG even for $s,\sigma=5$. For larger $s,\sigma$ values, fixed $s$-step CG no longer converges to the specified tolerance. Old adaptive $s$-step CG eventually converges for $\sigma=10$, but requires more than twice the total number of iterations as HSCG, limiting the potential benefit of any $s$-step approach. For $s,\sigma=15$, both fixed $s$-step CG and old adaptive $s$-step CG diverge. In the case $\varepsilon^{*}=10^{-6}$, the number of iterations required for convergence in old adaptive $s$-step CG grows drastically with $\sigma$, more than doubling from $\sigma=5$ to $\sigma=10$ and more than quadrupling from $\sigma=10$ to $\sigma=15$! For $\sigma=15$, this results in a greater number of outer loops iterations (i.e., more global synchronizations) than HSCG! In contrast, the total number of iterations required for the improved approach stays constant with $\sigma$ also in the $\varepsilon^{*}=10^{-6}$ case, and the number of outer loop iterations required also continues to decrease with increasing $\sigma$. For $\sigma=15$, the improved approach with both Newton and Chebyshev bases exhibits a decrease in outer loop iterations of more than $12\times$ versus HSCG.

We next show results for the text problem with matrix bcsstk09 (Figure LABEL:fig.bcsstk09n_0 and Table LABEL:tab.bcsstk09n0). For both tested $\varepsilon^{*}$ values, when $s,\sigma=5$ the improved adaptive approach requires more outer iterations to converge than both fixed $s$-step CG and old adaptive $s$-step CG. However, again the benefit of the more well-condition bases is clear for larger $s,\sigma$ values; as $s,\sigma$ increases, the total number of outer iterations actually increases rather than decreases for both fixed $s$-step CG and old adaptive $s$-step CG. Neither of those approaches converges to the prescribed level $\varepsilon^{*}=1.9\cdot 10^{-12}$ when $s,\sigma=15$. The old adaptive approach eventually converges to the prescribed level when $\varepsilon^{*}=10^{-6}$ and $\sigma=15$, although it requires almost $10\times$ the number of total iterations as HSCG; in fact, looking at how the total number of outer loop iteration required grows as $\sigma$ is increased, we see around a $2\times$ increase in the data movement cost for each increase in $\sigma$. In contrast, the number of total iterations required for the improved approach is constant with $\sigma$ using both Newton and Chebyshev bases, and thus the number of outer loop iterations decreases with increasing $\sigma$ as desired.

For gr_30_30, nos6, and mhdb416 (Tables LABEL:tab.gr3030n0,  LABEL:tab.nos6n0, and LABEL:tab.mhdb416n0, respectively), the story is similar; fixed $s$-step CG and old adaptive $s$-step CG fail to converge to the requested accuracy in many cases, especially for larger $s,\sigma$ values. In cases where these algorithms do converge for higher $s,\sigma$ values, the number of total iterations required (and also the number of outer loop iterations required) generally grow with increasing $s,\sigma$, sometimes drastically. For example, for the tests with mhdb416 for $\varepsilon^{*}=10^{-6}$ and $\sigma=15$, the number of total iterations required for convergence grows by almost $40\times$ versus HSCG, and thus requires more than $3\times$ the number of global synchronizations versus HSCG. In all cases the improved adaptive approach converges to the prescribed level and in most cases the number of iterations required for convergence remains relatively constant, resulting in a decrease in the number of outer loop iterations with increasing $\sigma$. We note that in some cases, the total number of outer loop iterations stagnates with increasing $\sigma$, or at least the rate of decrease is decreasing; see, for example, the experiments with gr_30_30 with $\sigma=10,15$ for both $\varepsilon^{*}$ values. Here, there is no decrease in total number of outer loop iterations going from $\sigma=10$ to $\sigma=15$ for either $\varepsilon^{*}$. This behavior can occur due to some fundamental limit on how large $s$ can be before the bases become too ill-conditioned. Regardless, this demonstrates another benefit of the adaptive approach; even if the user inputs (or an autotuner chooses) a $\sigma$ value that is too large, the adaptive approach can still be successful.

The experiments with the matrix nos1 (Table LABEL:tab.nos1n0, see also Section 5.1) are a notable exception. This is a particularly difficult problem even for HSCG, which requires more than $2N$ iterations to converge. Even for $\sigma=5$, the improved adaptive algorithm requires more than $2\times$ the number of iterations to converge than HSCG. We are not sure why this behavior occurs. Regardless, the benefits of the improved approach are still clearly observed. Neither fixed $s$-step CG nor old adaptive $s$-step CG converges to the prescribed accuracy when $s,\sigma=15$ is used, and even when convergence does occur, the number of outer loop iterations required grows quickly with $s,\sigma$. In all cases, the improved adaptive approach still converges to the prescribed level and still reduces the number of outer iterations required by at least a factor of $2$ versus HSCG. We also see that the number of outer loop iterations required in the improved approach still reliably decreases with increasing $\sigma$.

5.3 Discussion

The experimental results clearly show the benefit of the improved adaptive approach over both fixed $s$-step CG and the old adaptive $s$-step CG. First, the improved algorithm is more reliable. Whereas a poor choice of $s$ or $\sigma$ can cause failure to converge to the prescribed accuracy in the other algorithms, the improved algorithm converges to the desired tolerance in all tested cases. The number of outer loop iterations (global synchronizations) decreases with increased $\sigma$ in the improved approach, and in all cases, this cost is reduced relative to HSCG. Further, the improved algorithm does not require a heuristic choice of $c_{m+j+1}$, but instead automatically sets this parameter based on an inexpensive computation using quantities already available. Another benefit of the improved approach is that unlike the old adaptive approach, the total number of iterations required for convergence is independent of $\sigma$.

In the improved approach, we use information generated during the iterations to construct more well-conditioned polynomial bases, which allows for larger $s_{k}$ values to be used while still achieving the desired accuracy. We note that the adaptive approaches (both old and new) are developed only with the goal of attaining a prescribed accuracy rather than improving the convergence behavior. However, a side effect of using more well-conditioned bases is that the convergence behavior improves and becomes closer to that of HSCG. This is to be expected based on the theoretical analysis of the $s$-step Lanczos algorithm in [7], which is currently being investigated in order to extend Greenbaum’s classic results on the convergence of HSCG in finite precision [21] to $s$-step variants.

We also note that it is not clear whether the improved algorithm performs better using Newton or Chebyshev bases. In many cases, we observe that the Chebyshev basis provides the fastest convergence (even faster than HSCG in some cases), however, this is not universal; see, for example, the experiments for nos1 with $\sigma=10$ and $\varepsilon^{*}=10^{-6}$ (the final row of plots in Figure 1).

We note that in most experiments, the number of outer loop iterations required by fixed $s$-step CG, old adaptive $s$-step CG, and improved adaptive $s$-step CG is about the same for $s,\sigma=5$, with the fixed $s$-step and old adaptive approaches outperforming the improved approach in a few cases. This suggests that for small values like $s,\sigma=5$, it is unlikely to be beneficial to use the improved adaptive approach. In this case, a simple monomial basis is likely good enough to ensure that the desired accuracy is attained. Depending on the structure and size of the problem and the particular hardware on which the problem is solved, it may be that a small $s$ value is the best option for minimizing the time per iteration anyway. Providing a concrete answer regarding which approach to use when is nearly impossible due to the highly problem-dependent nature of the performance of iterative methods. Here we have only presented a few experiments on small matrices in order to demonstrate the numerical behavior of the algorithms. In order to provide a clearer picture, we stress that these algorithms should be implemented and compared on a high-performance parallel machine for problems from a wide variety of domains.

6 Conclusion and future work

In this work, we presented an improved adaptive $s$-step CG algorithm for solving SPD linear systems. The primary improvement over the previous adaptive $s$-step CG algorithm is due to the use of the method of Meurant and Tichý [33] for incremental estimation of the largest and smallest Ritz values, which are used to dynamically improve the conditioning of the $s$-step basis matrices (which has a large influence on the numerical behavior of the $s$-step methods) as the iterations proceed and also to automatically set a parameter that represents the ratio between the size of the error and the size of the residual. The improved algorithm can provide convergence behavior closer to that of HSCG while also allowing for the use of larger $s_{k}$ values without unacceptable loss of accuracy. We made an additional small improvement to the criterion for determining $s_{k}$ based on the observation that only a subset of the basis vectors are needed to compute the iterate updates in each subsequent step. Our numerical experiments verify that the improved algorithm provides increased reliability and, in many cases, a reduction in the total number of global synchronizations required to converge to the prescribed tolerance versus the old adaptive $s$-step CG algorithm.

The adaptive approach in $s$-step CG algorithms is an example of how finite precision analysis can provide insight leading to more numerically stable algorithms. It remains to show, however, that the adaptive $s$-step approaches provide a benefit to performance over classical Krylov subspace methods like HSCG. In practical implementations, there remain a few parameters in the algorithm which must be tuned to maximize performance, including the maximum $s_{k}$ value $\sigma$ and the maximum basis growth factor $f$. A high-performance parallel implementation and a thorough exploration of the design space remains critical future work.