# Mixed Precision Iterative Refinement with Sparse Approximate Inverse
Preconditioning

Erin Carson  Faculty of Mathematics and Physics, Charles University. Both
authors were supported by Charles University PRIMUS project no.
PRIMUS/19/SCI/11. The first author was additionally supported by Charles
University Research program no. UNCE/SCI/023 and by the Exascale Computing
Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy
Office of Science and the National Nuclear Security Administration.  Noaman
Khan22footnotemark: 2

###### Abstract

With the commercial availability of mixed precision hardware, mixed precision
GMRES-based iterative refinement schemes have emerged as popular approaches
for solving sparse linear systems. Existing analyses of these approaches,
however, are based on using full LU factorizations to construct
preconditioners for use within GMRES in each refinement step. In practical
applications, inexact preconditioning techniques, such as incomplete LU or
sparse approximate inverses, are often used for performance reasons.

In this work, we investigate the use of sparse approximate inverse
preconditioners based on Frobenius norm minimization within GMRES-based
iterative refinement. We analyze the computation of sparse approximate
inverses in finite precision and derive constraints under which user-specified
stopping criteria will be satisfied. We then analyze the behavior of and
convergence constraints for a five-precision GMRES-based iterative refinement
scheme that uses sparse approximate inverse preconditioning, which we call
SPAI-GMRES-IR. Our numerical experiments confirm the theoretical analysis and
illustrate the resulting tradeoffs between preconditioner sparsity and GMRES-
IR convergence rate.

###### keywords:

mixed precision, GMRES, iterative refinement, approximate preconditioners,
sparse linear systems

65F08, 65F10, 65F50, 65G50, 65Y99

##  1 Introduction

We consider the problem of solving linear systems $Ax=b$ with a nonsingular
$n\times n$ matrix $A$. With the recent emergence of commercially available
mixed precision hardware, there has been a renewed interest in the development
of mixed precision algorithms for numerical linear algebra. The benefit of
using low precision is that it is much faster and more energy efficient.
Compared to double precision, for example, using half precision tensor cores
available on modern NVIDIA GPUs leads to a theoretical $16\times$ improvement
in performance. Low precision can also reduce communication and memory
bottlenecks, since we need to move and store fewer bits. The downside is that
with fewer bits, we have less accuracy and a smaller range of representable
numbers. The size in bits, range, and unit roundoff $u$ are given for four
IEEE precisions in Table 1.

Table 1: Parameters for IEEE floating point precisions. The range denotes the order of magnitude of the largest and smallest positive normalized floating point numbers.  Type | Size | Range | Unit Roundoff $u$
---|---|---|---
half | 16 bits | $10^{\pm 5}$ | $2^{-11}\approx 4.9\cdot 10^{-4}$
single | 32 bits | $10^{\pm 38}$ | $2^{-24}\approx 6.0\cdot 10^{-8}$
double | 64 bits | $10^{\pm 308}$ | $2^{-53}\approx 1.1\cdot 10^{-16}$
quad | 128 bits | $10^{\pm 4932}$ | $2^{-113}\approx 9.6\cdot 10^{-35}$

The goal is thus to design mixed precision approaches which use lower
precision in select computations and higher precision in others, in such a way
that both (1) performance is improved and (2) the attainable accuracy remains
sufficient; see
[abdelfattah2021surveyand\cite[cite]{[\@@bibref{}{higham2022mixed}{}{}]}forrecentsurveys].
For the problem of solving linear systems, mixed precision iterative
refinement has been the focus of renewed attention. The general idea behind
iterative refinement is that after computing an initial solution $x_{0}$ to
$Ax=b$, we can refine this solution by repeatedly computing the residual
$r_{i}=b-Ax_{i}$, solving $Ad_{i}=r_{i}$, and then updating the approximate
solution $x_{i+1}=x_{i}+d_{i}$ until some convergence criterion is satisfied.
Typically, an LU factorization of $A$ is used to solve for $x_{0}$, and then
is reused in each iteration to solve for the correction $d_{i}$. We call this
LU-based variant “standard IR” or SIR for short.

There is a long history of using mixed precision within iterative refinement.
What we call “traditional” iterative refinement involves computing the
residuals in double the working precision $u$, which was used already by
Wilkinson in 1948 [41], and was analyzed by Wilkinson [42] and Moler [36].
Fixed precision iterative refinement, in which all computations are performed
in precision $u$, was analyzed by Jankowski and Woźniakowski [31] and Skeel
[40]. Motivated by the faster speed of single precision versus double, in the
early 2000s, there was a flurry of work in using lower precision in the
computation of the LU factorization, which is the most expensive part of the
computation, and the working precision elsewhere; see, for example, [33], [1,
Section 9]. See [12, Table 1.1] for a summary.

The inclusion of half precision in modern GPUs inspired the development of
iterative refinement schemes that use three or more hardware precisions. In
[12], the authors define an iterative refinement scheme which uses three
potentially different precisions: $u_{f}$ for the factorization, $u$ for the
working precision, and $u_{r}$ for the residual computation. To allow for
general solvers for the correction term $d_{i}$, the authors also introduce a
fourth precision $u_{s}$, which denotes the “effective precision” of the
solve. For standard iterative refinement, the effective precision of the solve
(which involves triangular solves with the LU factors computed in precision
$u_{f}$), we have $u_{s}=u_{f}$. Then assuming that $u_{f}\geq u$ and
$u_{r}\leq u^{2}$, the relative forward and backward errors in standard
iterative refinement will converge to level $u$ when $\kappa_{\infty}(A)\leq
u_{f}^{-1}$.

The reason for introducing this effective precision becomes clear when we use
a more accurate solver. The GMRES-based iterative refinement scheme (GMRES-IR)
introduced in [11] uses GMRES left preconditioned by the computed LU factors
in order to solve for $d_{i}$ in each refinement step. Assuming that GMRES is
executed in working precision $u$, with the matrix-vector products with the
preconditioned matrix computed in precision $u^{2}$, GMRES-IR is guaranteed to
give forward and backward errors to the working precision for more ill-
conditioned systems than standard iterative refinement. Again under the
assumption that $u_{f}\geq u$ and $u_{r}\leq u^{2}$, the relative forward and
backward errors in GMRES-IR will converge to level $u$ when
$\kappa_{\infty}(A)\leq u^{-1/2}u_{f}^{-1}$.

The requirement that the preconditioned matrix is applied in double the
working precision within GMRES is unattractive from a performance perspective.
In [2], the authors introduce and analyze a five-precision variant of GMRES-
IR. In addition to the working precision $u$, factorization precision $u_{f}$,
and residual precision $u_{r}$ for the refinement scheme, they also add
precisions $u_{g}$ for the working precision within GMRES and $u_{p}$ for the
precision in which the preconditioned matrix is applied to a vector within
GMRES. Setting $u=u_{g}=u_{p}$ is a variant commonly used in practice. The
cost is that this variant is only guaranteed to converge for a smaller range
of condition numbers; again assuming that $u_{f}\geq u$ and $u_{r}\leq u^{2}$,
one obtains relative forward and backward errors to the level of the working
precision for matrices with $\kappa_{\infty}(A)\leq u^{-1/3}u_{f}^{-2/3}$,
although this restriction is likely overly pessimistic in practice.

Existing analyses of GMRES-IR makes the assumption that an LU factorization is
computed for use as a left preconditioner in GMRES in each refinement step.
This may be undesirable from a performance perspective in the case of sparse
$A$. Even if $A$ is very sparse, its LU factors may have considerable fill-in.
In practice, inexact preconditioners, such as incomplete LU factorization or
sparse approximate inverses are often used. SPAI preconditioners have the
advantage that their construction is highly parallelizable; in theory, each
column of the approximate inverse can be computed independently. Further,
their application to a vector involves only a sparse matrix-vector
multiplication, rather than the two triangular solves required for applying LU
factors.

In this work, we explore the potential for the use of sparse approximate
inverse preconditioners within five-precision GMRES-IR. In Section 2 we give a
brief background on sparse approximate inverse preconditioning and mixed
precision iterative refinement. In Section 3, we analyze the error in
constructing the sparse approximate inverse in finite precision and then
analyze the behavior of five-precision GMRES-based iterative refinement
schemes based on sparse approximate inverses, giving theoretical results on
the attainable accuracy and criteria for convergence of the refinement scheme.
In Section 4 we present numerical experiments and in Section 5 we discuss open
problems and future work.

##  2 Background and related work

###  2.1 Notation

We first introduce notation which will be used in the remainder of the text.
Of particular importance will be various condition numbers. For a given norm
$p$, a matrix $A$, and a vector $x$, we define

$\kappa_{p}(A)=\|A^{-1}\|_{p}\|A\|_{p},\quad\text{cond}_{p}(A)=\||A^{-1}||A|\|_{p},\quad\text{cond}_{p}(A,x)=\frac{\||A^{-1}||A||x|\|_{p}}{\|x\|_{p}},$

where $|A|=(|a_{ij}|)$. In the case that $p$ is not specified, the infinity
norm should be assumed. We will use $u$’s to denote unit roundoffs, where
subscripts on $u$ will distinguish between various precisions. For rounding
error analysis, we will frequently use the notation

$\gamma_{k}=\frac{ku}{1-ku},\quad\tilde{\gamma}_{k}=\frac{cku}{1-cku},$

where $c$ is a small constant independent of problem dimension. A superscript
on $\gamma$ indicates that $u$ has that superscript as a subscript, e.g.,
$\gamma_{k}^{f}=ku_{f}/(1-ku_{f})$. Quantities computed in finite precision
will be denoted with hats. In order to make clear the distinction between
iterative refinement (the outer solver) and GMRES (the inner solver), we will
always use the word “steps” when referring to iterative refinement and
“iterations” when referring to GMRES.

###  2.2 Sparse approximate inverse preconditioners

The idea behind sparse approximate inverse (SPAI) preconditioning is to
explicitly construct a matrix $M\approx A^{-1}$. This has advantages within
Krylov subspace methods since the application of the preconditioner involves
only a matrix-vector product, rather than, e.g., the triangular solves
involved when using LU-based preconditioners. There are many possible
algorithms for computing $M$; see [9, 7] for a survey and references.

For computing a sparse approximate inverse in unfactored form (i.e., a single
matrix $M$), a popular approach is based on Frobenius norm minimization, in
which $M$ is computed as the solution to
$\min_{\mathcal{J}\in\mathcal{S}}\|I-AM\|_{F}$, where
$\mathcal{J}\in\mathbb{B}^{n\times n}$ is a prescribed binary sparsity pattern
in the set of all possible binary sparsity patterns
$\mathcal{S}\in\mathbb{B}^{n\times n}$. The key feature of this minimization
problem is that it can be decoupled as

(1)
$\min_{\mathcal{J}\in\mathcal{S}}\|I-AM\|_{F}^{2}=\sum_{k=1}^{n}\min_{\mathcal{J}_{k}\in\mathcal{S}_{k}}\|e_{k}-Am_{k}\|_{2}^{2},$

where $\mathcal{J}_{k}$, $m_{k}$, and $e_{k}$ are the $k$th columns of
$\mathcal{J}$, $M$, and $I$, respectively. The computation of $M$ thus reduces
to solving a linear least squares problem for each column $m_{k}$ of $M$. A
primary advantage in terms of performance is that these linear least squares
problems are independent, and thus can in principle be solved in parallel.

Early work based on this technique used a fixed sparsity pattern
$\mathcal{J}$. The set $\mathcal{J}_{k}$ is called the extraction set, giving
the column indices of $A$ that are relevant for solving for a column $m_{k}$.
The nonzero rows of the submatrix $A(:,\mathcal{J}_{k})$ are denoted by the
so-called “shadow” of $\mathcal{J}_{k}$,

$\mathcal{I}_{k}=\left\\{i\in\\{1,\ldots,n\\}:\sum_{j\in\mathcal{J}_{k}}|a_{ij}|\neq
0\right\\},$

where $a_{ij}$ is the $(i,j)$ entry of $A$. Thus each term in the summation on
the right in (1) can be reduced to

(2)
$\min_{\mathcal{J}(\bar{m}_{k})=\mathcal{J}_{k}}\|\bar{e}_{k}-\bar{A}_{k}\bar{m}_{k}\|_{2},$

where
$\bar{A}_{k}=A(\mathcal{I}_{k},\mathcal{J}_{k})\in\mathbb{R}^{|\mathcal{I}_{k}|,|\mathcal{J}_{k}|}$,
$\bar{m}_{k}=m_{k}(\mathcal{J}_{k})\in\mathbb{R}^{|\mathcal{J}_{k}|}$,
$\bar{e}_{k}=e_{k}(\mathcal{I}_{k})\in\mathbb{R}^{|\mathcal{I}_{k}|}$, and
$\mathcal{J}(\bar{m}_{k})$ is the binary sparsity pattern of $\bar{m}_{k}$.
Because the resulting least squares problems are small, they can be solved
directly, for example, via QR factorization.

A deficiency of this approach is that it is difficult in general to select a
priori a sparsity pattern which will produce an effective preconditioner.
Common choices include the sparsity pattern of $A$, $A^{T}$, or a power of a
sparsified $A$, although these are not in general guaranteed to produce
effective preconditioners. To overcome this, many authors have developed
iterative approaches, in which one starts with an initial sparsity structure
and adds nonzeros to this pattern until the constraint
$\|e_{k}-Am_{k}\|_{2}\leq\varepsilon$ is satisfied for some threshold
$\varepsilon$ or a maximum number of nonzeros has been reached. For algorithms
of this type, see, e.g., the work of Cosgrove et al. [15], Grote and Huckle
[23], and Gould and Scott [21].

The algorithm of Grote and Huckle [23] is among the most successful and
commonly used approaches for computing an SPAI preconditioner [9], and is the
focus of our present work. Motivated by the difficulty of selecting a sparsity
pattern that results in a good preconditioner a priori, the idea in [23] was
to develop an adaptive technique which captures the most important nonzero
indices dynamically.

A particular variant of Grote and Huckle’s approach, which appears in [39,
Algorithm 4] is shown in Algorithm 1. The algorithm takes as input the matrix
$A$, an initial binary sparsity pattern $\mathcal{J}$, the convergence
tolerance $\varepsilon$, $\alpha$, which gives the maximum number of
iterations for each column, and $\beta$, which gives the maximum number of
nonzeros that are added to the pattern in each iteration.

Given an initial sparsity pattern $\mathcal{J}$, for each column, the
algorithm solves the least squares problem (2) and computes the residual
$\bar{s}_{k}$ (lines 7-12). If the 2-norm of the residual is below the
threshold $\varepsilon$, we are finished for this column. Otherwise, we
proceed with adding entries to $\mathcal{J}$.

The nonzeros entries in $\bar{s}_{k}$ are given by the index set
$\mathcal{L}_{k}$ in line 16. For every entry $\ell$ of $\mathcal{L}_{k}$, we
can define the set $\mathcal{N}_{\ell}$ which contains the column indices of
the nonzero entries in the $\ell$th row of $A$ which are not already in
$\mathcal{J}_{k}$. The union of the sets $\mathcal{N}_{\ell}$ gives the set
$\tilde{\mathcal{J}}_{k}$, which is the set of potential indices to be added
to $\mathcal{J}_{k}$. We want to select only a subset of the “most important”
of these indicies to add.

There are various ways of determining which indices are most important. The
approach of Grote and Huckle is to consider a univariate minimization problem,
through which the quantity $\rho_{jk}$ computed in line 23 gives a measure of
the 2-norm of the new residual if index $j$ is added to $\mathcal{J}_{k}$. A
well-known heuristic (see, e.g., [6]) is to mark indices as “acceptable” if
their $\rho_{jk}$ is less than the arithmetic mean $\bar{\rho}_{k}$ over all
$j$. In each of the $\alpha$ iterations, we add up to $\beta$ of the best
(smallest $\rho_{jk}$) acceptable indices (lines 28-30).

We note that the QR factorization in line 10 does not need to be recomputed
entirely in each step, but can instead be updated using the QR factorization
computed in the previous step and the entries added to $\bar{A}_{k}$; see [23,
Eqns. (14) and (15)]. A benefit of the adaptive approach of Grote and Huckle
is that the algorithm is robust with respect to choice of initial sparsity
pattern $\mathcal{J}$. A common choice, also recommended in [23], is to take
the initial $\mathcal{J}$ to be the identity matrix. Typical choices of the
other parameters are $\varepsilon\in[0.1,0.5]$, $\alpha\in\\{1,\ldots,5\\}$,
and $\beta\in\\{3,\ldots,8\\}$ [39, Section 3.1.3].

Despite that each column can be computed in parallel, constructing an SPAI
preconditioner is often costly, especially for large-scale problems; see,
e.g., [6, 9, 14, 19]. The memory requirements for SPAI scale quadratically and
the computational cost cubically in the number of nonzeros per row [19]. Thus
the use of low precision for storage and computation has the potential to
significantly reduce this cost. The development of efficient sparse
approximate inverse computations for modern hardware, especially GPUs, has
been the subject of much recent work; see, e.g., [20, 35, 17, 24].

We note that there are also techniques based on incomplete biconjugation which
can be used to produce sparse triangular factors, so that the approximate
inverse is a product of two matrices approximating the LU factors of $A$; see,
e.g., [8]. Here we do not consider this approach nor the many other variants
of sparse approximate inverses, including factorized sparse approximate
inverses [29], modified sparse approximate inverses [30], or incomplete sparse
approximate inverses [4]. We note that these could be interesting venues for
extending the present work.

Algorithm 1 Variant of sparse approximate inverse (SPAI) construction ([39,
Algorithm 4])

1:$A\in\mathbb{R}^{n\times n}$, $\mathcal{J}\in\mathbb{B}^{n\times n}$,
$\alpha\geq 0$, $\beta\geq 0$, $\varepsilon>0$

2:Right preconditioner $M\approx A^{-1}$, $M\in\mathbb{R}^{n\times n}$

3:for $k=1$ to $n$ do

4: $e_{k}=I(:,k)$

5: $\mathcal{J}_{k}=\mathcal{J}(:,k)$

6: for step $=0$ to $\alpha$ do

7:
$\mathcal{I}_{k}=\bigg{\\{}i\in\left\\{1,\ldots,n\right\\}\colon\sum_{j\in\mathcal{J}_{k}}\left|a_{ij}\right|\neq
0\bigg{\\}}$

8: $\bar{A}_{k}=A(\mathcal{I}_{k},\mathcal{J}_{k})$

9: $\bar{e}_{k}=e_{k}(\mathcal{I}_{k})$

10: Compute QR factorization $\bar{A}_{k}=\bar{Q}\bar{R}$

11: $\bar{m}_{k}=\bar{R}^{-1}\bar{Q}^{T}\bar{e}_{k}$

12: $\bar{s}_{k}=\bar{A}_{k}\bar{m}_{k}-\bar{e}_{k}$

13: if $\|\bar{s}_{k}\|_{2}$ $\leq$ $\varepsilon$ then

14: break

15: end if

16: $\mathcal{L}_{k}=\mathcal{I}_{k}$ $\cup$ $\left\\{k\right\\}$

17: for $\ell$ $\in$ $\mathcal{L}_{k}$ do

18: $\mathcal{N}_{\ell}=\\{j\colon{a}_{\ell j}\neq 0\\}$

19: end for

20:
$\tilde{\mathcal{J}}_{k}=\bigcup_{\ell\in\mathcal{L}_{k}}\mathcal{N}_{\ell}$

21: $\tilde{\rho}_{k}=0$

22: for $j$ $\in$ $\tilde{\mathcal{J}}_{k}$ do

23:
${\rho}_{jk}=\left(\|\bar{s}_{k}\|^{2}_{2}-\frac{[\bar{r}^{T}_{k}A_{j}(\mathcal{I}_{k})]^{2}}{\|A_{j}(\mathcal{I}_{k})\|^{2}_{2}}\right)^{\frac{1}{2}}$

24: $\tilde{\rho}_{k}=\tilde{\rho}_{k}$ \+ $\tilde{\rho}_{jk}$

25: end for

26: $\tilde{\rho}_{k}=\frac{\tilde{\rho}_{k}}{|\tilde{\mathcal{J}}_{k}|}$

27: for idx $=1$ to $\beta$ do

28: $j={argmin}_{j\in\tilde{\mathcal{J}}_{k}}$ $\rho_{jk}$

29: $\mathcal{J}_{k}={\mathcal{J}}_{k}$ $\cup$
$\left\\{j\colon\rho_{jk}\leq\tilde{\rho}_{k}\right\\}$

30: $\tilde{\mathcal{J}}_{k}=\tilde{\mathcal{J}}_{k}$\$\left\\{j\right\\}$

31: end for

32: end for

33: $m_{k}(\mathcal{J}_{k})=\bar{m}_{k}$

34:end for

###  2.3 Mixed precision iterative refinement

In Algorithm 2 we present a general three-precision iterative refinement
scheme. There are three explicit hardware precisions involved: $u_{f}$, in
which the initial approximate solution $x_{0}$ is computed, $u_{r}$, in which
the residual is computed, and $u$, the working precision in which all other
computations are performed and all quantities are stored.

Algorithm 2 General Iterative Refinement Scheme

1:$n$-by-$n$ nonsingular matrix $A$ and length-$n$ right-hand side vector $b$.

2:Approximate solution $\widehat{x}$ to $Ax=b$.

3:Compute LU factorization of $A$ in precision $u_{f}$.

4:Solve $Ax_{0}=b$ in precision $u_{f}$; store $x_{0}$ in precision $u$.

5:for $i=0,1\ldots$ until convergence do

6: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store $r_{i}$ in precision
$u$.

7: Solve $Ad_{i}=r_{i}$ at effective precision $u_{s}$; store $d_{i}$ in
precision $u$.

8: Update $x_{i+1}=x_{i}+d_{i}$ in precision $u$.

9:end for

The effective precision $u_{s}$ for the solve in line 7 of Algorithm 2 depends
on the precisions and the solver used. In particular, $u_{s}$ is defined such
that the solver satisfies the conditions:

(3) $\displaystyle\widehat{d}_{i}$ $\displaystyle=(1+u_{s}E_{i})d_{i},\quad
u_{s}\|E_{i}\|_{\infty}<1,\quad\text{and}$ (4)
$\displaystyle\|\widehat{r}_{i}-A\widehat{d}_{i}\|_{\infty}$
$\displaystyle\leq
u_{s}(c_{1}\|A\|_{\infty}\|\widehat{d}_{i}\|_{\infty}+c_{2}\|\widehat{r}_{i}\|_{\infty}),$

where $E_{i}$, $c_{1}$, and $c_{2}$ are functions of $n$, $A$,
$\widehat{r}_{i}$, and $u_{s}$. These two conditions are used in analyzing the
relative forward error and normwise backward errors, respectively. For
simplicity, we will only consider normwise forward and backward errors in this
work, but direct the reader to [12, Section 5] for detailed treatment of
componentwise backward error. Given a choice of solver with a particular
$u_{s}$ along with precisions $u_{f}$, $u$, and $u_{r}$, the authors in [12]
prove the following results about the behavior of the resulting iterative
refinement scheme.

######  Theorem 2.1.

[12, Corollary 3.3] Let Algorithm 2 be applied to a linear system $Ax=b$ where
$A$ is $n\times n$ and nonsingular with at most $q$ nonzeros per row and
assume the solver used satisfies (3). Then for refinement step $i\geq 0$, as
long as

(5)
$\phi_{i}=2u_{s}\min({\text{cond}(A)},\kappa_{\infty}(A)\mu_{i})+u_{s}\|E_{i}\|_{\infty}$

is less than 1, the forward error is reduced on the ith step by a factor of
approximately $\phi_{i}$ until an iterate $\widehat{x}$ is obtained for which

$\frac{\|x-\widehat{x}\|_{\infty}}{\|x\|_{\infty}}\lesssim
4qu_{r}{\text{cond}(A,x)}+u.$

The quantity $\mu_{i}$ in the above Theorem is defined such that

$\|A(x-\widehat{x}_{i})\|_{\infty}=\mu_{i}\|A\|_{\infty}\|x-\widehat{x}_{i}\|_{\infty},$

where $\mu_{i}$ satisfies $\kappa_{\infty}(A)^{-1}\leq\mu_{i}\leq 1$. The
insight in [11] is that $\mu_{i}$ will be very small at the beginning of the
refinement process, and will only grow close to its maximum once the process
is close to convergence. Thus the size of the quantity
$u_{s}\|E_{i}\|_{\infty}$ will usually dictate convergence in practice.

######  Theorem 2.2.

[12, Corollary 4.2] Let Algorithm 2 be applied to a linear system $Ax=b$ where
$A$ is $n\times n$ and nonsingular with at most $q$ nonzeros per row and
assume the solver used satisfies (3). Then for refinement step $i\geq 0$, as
long as

(6) $\phi_{i}=(c_{1}\kappa_{\infty}(A)+c_{2})u_{s}$

is less than 1, the residual is reduced in each step by a factor of
approximately $\phi_{i}$ until an iterate $\widehat{x}$ is obtained for which

$\|b-A\widehat{x}\|_{\infty}\lesssim
qu(\|b\|_{\infty}+\|A\|_{\infty}\|\widehat{x}\|_{\infty}),$

which indicates that $\widehat{x}$ is a backward stable solution to the
working precision.

Thus the quantities $u_{s}$, $\|E_{i}\|_{\infty}$, $c_{1}$, and $c_{2}$, which
come from the particular solver, completely determine the convergence
conditions for the forward and backward errors in iterative refinement. The
authors in [12] bound these quantities for two different solvers: 1)
triangular solves using the computed LU factors (SIR), and 2) left-
preconditioned GMRES, where the computed LU factors are used as
preconditioners and the preconditioners and preconditioned coefficient matrix
are applied to vectors in double the working precision (GMRES-IR). The authors
in [2] recently extended this analysis to cover the more general case where
the solver is GMRES left-preconditioned by the computed LU factors run in a
working precision $u_{g}$, where the preconditioners and preconditioned
coefficient matrix are applied to vectors in a precision $u_{p}$. The variant
of GMRES-IR in which $u=u_{g}=u_{p}$ is attractive from a performance
standpoint and is commonly used in practice.

We summarize the size of these quantities in Table 2. Here, $f$ denotes a
generic function, $n$ is the problem dimension, $k$ is the maximum number of
GMRES iterations in any refinement step, $\rho_{n}$ is the growth factor in
LU, and $\tilde{A}=\widehat{U}^{-1}\widehat{L}^{-1}A$ for the computed factors
$\widehat{L}$ and $\widehat{U}$.

Table 2: Summary of sizes of quantities in the assumptions (3)-(4) for various solvers as given in [2]. In the table, $f$ denotes a generic function, $n$ is the problem dimension, $k$ is the maximum number of GMRES iterations in any refinement step, $\rho_{n}$ is the growth factor in LU, and $\tilde{A}=\widehat{U}^{-1}\widehat{L}^{-1}A$ for the computed factors $\widehat{L}$ and $\widehat{U}$. | SIR | GMRES-IR ($u_{g}=u,u_{p}=u^{2}$) | GMRES-IR ($u_{g}=u_{p}=u$)
---|---|---|---
$u_{s}$ | $u_{f}$ | $u$ | $u_{g}+u_{p}\kappa_{\infty}(A)$
$\|E_{i}\|_{\infty}$ | $f(n,\rho_{n})$ | $f(n,k,\rho_{n})\kappa_{\infty}(\tilde{A})$ | $f(n,k,\rho_{n})\kappa_{\infty}(\tilde{A})$
$c_{1}$ | $f(n,\rho_{n})$ | $f(n,k,\rho_{n})\|\tilde{A}\|_{\infty}$ | $f(n,k,\rho_{n})\|\tilde{A}\|_{\infty}$
$c_{2}$ | $f(n,\rho_{n})$ | $f(n,k,\rho_{n})\kappa_{\infty}(A)$ | $f(n,k,\rho_{n})\kappa_{\infty}(A)$

Substituting the values for $u_{s}$, $\|E_{i}\|_{\infty}$, $c_{1}$, and
$c_{2}$ from Table 2 into Theorems 2.1 and 2.2 shows that for SIR, the forward
and backward errors will converge as long as $\kappa_{\infty}(A)\ll
u_{f}^{-1}$, with a limiting value of $qu$ for the backward error and a
limiting value of $qu_{r}\text{cond}(A,x)+u$ for the forward error.

GMRES-IR with $u_{g}=u$ and $u_{p}=u^{2}$ will have the same limiting values,
but looser constraints on condition number. Using the bound
$\kappa_{\infty}(\tilde{A})\ll\kappa_{\infty}(A)^{2}u_{f}^{2}$ (see [12, Eqn.
(8.3)]), for the forward error, we require $\kappa_{\infty}(A)\ll
u^{-1/2}u_{f}^{-1}$. For the backward error, we can use the bound on
$\|\tilde{A}\|_{\infty}$ from [11, Eqn. (3.2)] to give the constraint
$\kappa_{\infty}(A)\ll u^{-1/2}u_{f}^{-1/2}$; see also [2, Table 2.1]. Note
that in the case that $u_{r}=u^{2}$, the constraint for the convergence of the
backward error to the working precision also becomes $\kappa_{\infty}(A)\leq
u^{-1/2}u_{f}^{-1}$ since the backward error is bounded by the forward error.

From the analysis in [2], for GMRES-IR with $u=u_{g}=u_{p}$, using the same
bound on $\kappa_{\infty}(\tilde{A})$, the constraint for convergence of the
forward error to the limiting value is $\kappa_{\infty}(A)\leq
u^{-1/3}u_{f}^{-2/3}$ and that for the backward error is
$\kappa_{\infty}(A)\leq u^{-1/3}u_{f}^{-1/3}$. As above, the constraint for
convergence of the backward error is stricter, and thus we can take the
constraint for the convergence of both forward and backward errors to the
limiting value to be $\kappa_{\infty}(A)\leq u^{-1/3}u_{f}^{-2/3}$ in the case
$u_{r}=u^{2}$. We stress that these constraints are overly strict in practice,
and convergence is often observed even in cases where these conditions are
violated; see, e.g., [38].

In order to simplify the presentation and allow easy extension of results, we
will restrict ourselves to left preconditioning, since this has been used in
previous works on GMRES-based iterative refinement; see, e.g., [2, 12]. In [5,
Remark 3.5], Arioli and Duff comment that right-preconditioned GMRES is only
backward stable if the preconditioner satisfies certain constraints. For the
case of sparse approximate inverses, this condition is
$\|\widehat{M}^{-1}\|_{2}/\|A\|_{2}\approx 1$, where $\widehat{M}$ is the
computed sparse approximate inverse. If we could guarantee this, then it could
be shown that right-preconditioned GMRES is backward stable and we conjecture
that similar bounds could be proven on the quantities in Table 2. We note that
we could also use right-preconditioned FGMRES with an SPAI preconditioner
within GMRES-based iterative refinement, although such a result has not yet
been formally proven. We leave these tasks as future work, but note that it is
likely that right-preconditioning will work in practice. In the context of
SPAI, right-preconditioning may be preferred when $A^{-1}$ has some dense
rows, since this will yield a sparser approximate inverse [23, pp. 843].

###  2.4 Related Work

Along with the resurgence of interest in mixed precision iterative refinement
have come a number of works that build on and expand the work on GMRES-based
iterative refinement, e.g., [11, 12, 1, 2]. In [26], Higham and Mary develop a
new general preconditioning technique based on a low-rank approximation of the
error $E=\widehat{U}^{-1}\widehat{L}^{-1}A-I$, where $\widehat{L}$ and
$\widehat{U}$ are approximate LU factors of $A$. A randomized SVD is used to
produce this low-rank approximation, and the preconditioner is applied in the
context of GMRES-IR. This technique can apply to the case where the
approximation in LU is due to low precision computation, as well as other
sources of error including using an incomplete LU factorization and a block
low-rank LU factorization. See [3, 27] for other theoretical and experimental
results on the use of block low-rank LU factorization preconditioners within
iterative refinement.

In [34], Lindquist, Luszczek, and Dongarra have recently evaluated mixed
precision restarted GMRES (equivalent in some sense to a variant of GMRES-IR)
on a GPU-accelerated node with various preconditioners. They found that the
use of incomplete LU with zero fill-in (ILU(0)) resulted in unfavorable
performance (a slowdown on average). The authors attribute this partially to
the fact that sparse triangular solves have limited parallelism for the GPU to
exploit.

In the recent work [3], Amestoy et al. extend the work on GMRES-based
iterative refinement to LU factorization methods more appropriate for sparse
matrices, including those based on block low-rank approximations and on static
pivoting strategies.

The present work differs from these recent approaches in that we depart from
the restriction to preconditioning based on LU factorization and instead look
at sparse approximate inverses as preconditioners within mixed precision
GMRES-based iterative refinement. In the following section, we will first
analyze the numerical properties of SPAI preconditioners computed in a given
precision with given parameters, and then evaluate the constraints for
convergence of GMRES-IR with SPAI preconditioners.

##  3 Analysis of SPAI-based iterative refinement

Our ultimate goal is to prove conditions under which the SPAI-GMRES-IR
algorithm (Algorithm 3) will converge and to determine how the precision
$u_{f}$ and approximation parameter $\varepsilon$ should be related. We first
give bounds on the quality of the SPAI preconditioner when it is computed in
finite precision, analogous to those given for the case of exact arithmetic in
[23].

Suppose we want to construct a sparse approximate inverse of a nonsingular
matrix $A$ for use as a left preconditioner, which we will denote $P$. That
is, instead of solving $Ax=b$ we will solve $PAx=Pb$. We will construct $P$
using the algorithm of Grote and Huckle [23] (see Algorithm 1), executed in
some precision $u_{f}$. Algorithm 1 as stated produces a right preconditioner
$M$. We can obtain a left preconditioner by executing the algorithm on $A^{T}$
and then transposing the obtained approximate inverse, i.e., $P=M^{T}$. We
will denote the version of $M$ computed in finite precision as $\widehat{M}$.

In contrast to how the method is often used in practice, we will not specify a
maximum number of nonzeros here (i.e., in the extreme case we could allow a
fully dense inverse). More important for our purposes will be the assumption
that for each column $\widehat{m}_{k}$ of $\widehat{M}$, it holds that

(7) $\|e_{k}-A^{T}\widehat{m}_{k}\|_{2}\leq\varepsilon.$

The way the SPAI algorithm works, we keep iterating, computing a QR
factorization of some submatrix $\bar{A}^{T}$ of $A^{T}$ (selected by
specifying a particular set of row and column indices) and solving the least
squares problem
$\min_{\bar{m}_{k}}\|\bar{A}^{T}\bar{m}_{k}-\bar{e}_{k}\|_{2}$. The vector
$\bar{m}_{k}$ then gives the nonzero entries in column $k$ of $M$
(corresponding to the same indices as the rows of $A^{T}$ selected to
construct $\bar{A}^{T}$). Row and column indices are added to the set, and the
iteration is continued until
$\min_{\bar{m}_{k}}\|\bar{A}^{T}\bar{m}_{k}-\bar{e}_{k}\|_{2}\leq\varepsilon$,
which indicates that (7) holds.

We assume that all computations within the SPAI algorithm are performed in a
precision $u_{f}$. We make the simplifying assumption that the QR
factorization of $\bar{A}^{T}$ is computed fully, but note again that in
practice one only needs to compute a QR factorization of a smaller block
corresponding to the newly added indices; see, e.g., [39, Algorithm 5]. We
further assume that the QR factorization is computed using Householder QR (or
a method with similar numerical properties, like TSQR [18]).

Now, assuming we impose no maximum sparsity pattern on $\widehat{M}$, we seek
to determine under what constraint on $u_{f}$ we can guarantee that
$\|\widehat{r}_{k}\|_{2}\leq\varepsilon$ can be satisfied, where
$\widehat{r}=fl_{u_{f}}(e_{k}-A^{T}\widehat{m}_{k})$ for the computed
$\widehat{m}_{k}$. Here we will assume that no overflow or underflow occurs
during the SPAI computation in precision $u_{f}$. We note that to mitigate the
risk of this one could use sophisticated scaling techniques such as those
described in [28] and [13]; we use a similar column scaling technique in our
experiments in Section 4. To account for finite precision error, it suffices
to assume that we are on the final iteration for column $k$, after which we
will quit since the condition $\|\widehat{r}_{k}\|_{2}$ is satisfied. Using
the backward error result [25, Theorem 20.3], which says that the computed
solution $\widehat{m}_{k}$ is the exact least squares solution to

$\min_{\widehat{m}_{k}}\|(e_{k}+\Delta e_{k})-(A^{T}+\Delta
A^{T})\widehat{m}_{k}\|_{2},$

where $|\Delta A^{T}|\leq n\tilde{\gamma}_{n}^{f}ee^{T}|A^{T}|$ and $|\Delta
e_{k}|\leq n\tilde{\gamma}_{n}^{f}ee^{T}|e_{k}|$, we can then bound

(8) $\|e_{k}-A^{T}\widehat{m}_{k}\|_{2}\leq
n^{2}\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}+\left(1+n^{2}\tilde{\gamma}_{n}^{f}\text{cond}_{2}(A)\right)\|r_{k}\|_{2}+O(u_{f}^{2}),$

where $r_{k}=e_{k}-A^{T}m_{k}$. Then using standard rounding error analysis,
we have

$\|\widehat{r}_{k}\|_{2}\leq\|e_{k}-A^{T}\widehat{m}_{k}\|_{2}+\|\Delta\widehat{r}_{k}\|_{2},$

where
$\|\Delta\widehat{r}_{k}\|_{2}\leq\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}$.
Then

(9) $\|\widehat{r}_{k}\|_{2}\leq
n^{2}\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}+\left(1+n^{2}\tilde{\gamma}_{n}^{f}\text{cond}_{2}(A)\right)\|r_{k}\|_{2}+O(u_{f}^{2}),$

where the $\|\Delta\widehat{r}_{k}\|_{2}$ factor has been absorbed into the
$\tilde{\gamma}_{n}^{f}$ in the first term on the right-hand side.

Recall that our goal is to determine under what constraint on $u_{f}$ we can
say that $\|\widehat{r}_{k}\|_{2}\leq\varepsilon$ is achievable, with no
restriction on the number of nonzeros in the columns $\widehat{m}_{k}$. In
theory, we could keep increasing the size of the subproblem in solving for the
column of the inverse until we have the full matrix $A^{T}$, which would mean
that in exact arithmetic, $m_{k}$ is the exact $k$th column of the inverse of
$A^{T}$, and thus $r_{k}=0$. Thus the second term on the right-hand side of
(9) will not limit us from satisfying the stopping criterion
$\|\widehat{r}_{k}\|_{2}\leq\varepsilon$, and it can thus be ignored for our
purposes. We are then left with (ignoring higher order terms)

$\|\widehat{r}_{k}\|_{2}\leq
n^{2}\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}.$

Then in order to guarantee that we will eventually iterate to a solution with
$\|\widehat{r}_{k}\|_{2}\leq\varepsilon$, we should have

(10)
$n^{2}\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}\leq\varepsilon.$

In other words, the problem must not be so ill-conditioned relative to $u_{f}$
such that we can not guarantee that we incur an error less than $\varepsilon$
just computing the residual. We note that the dimensional constants here can
be a large overestimate, since we in practice solve the least squares problem
with the smaller submatrix $\bar{A}^{T}$. We can turn this into the more
descriptive (but less rigorous) bound

$\text{cond}_{2}(A^{T},m_{k})\lesssim\varepsilon u_{f}^{-1},$

as well as an a priori bound of the form

(11) $\text{cond}_{2}(A^{T})\lesssim\varepsilon u_{f}^{-1}.$

We note that this is likely quite pessimistic in practice. Note that
$\text{cond}_{2}(A^{T})$ can be considerably smaller than $\kappa_{2}(A)$.
Perhaps a more useful way to think about (11) is that with a given matrix $A$
and a chosen precision $u_{f}$, in order to guarantee success of the SPAI
algorithm, one must set $\varepsilon$ such that

$u_{f}\text{cond}_{2}(A^{T})\leq\varepsilon.$

This analysis echoes what we may intuitively think: the larger we make
$\varepsilon$ (meaning the more approximate the inverse), the larger the
$u_{f}$ we can tolerate. If we require smaller $\varepsilon$, meaning we want
a better quality preconditioner, then we must have a smaller $u_{f}$ to match.

Given that (10) is satisfied, we know that we can find a computed solution
$\widehat{m}_{k}$ such that $\|\widehat{r}_{k}\|_{2}\leq\varepsilon$. Writing
$\widehat{r}_{k}=e_{k}-A^{T}\widehat{m}_{k}+\Delta\widehat{r}_{k}$, we have

$\displaystyle\|e_{k}-A^{T}\widehat{m}_{k}\|_{2}$
$\displaystyle\leq\|\widehat{r}_{k}\|_{2}+\|\Delta\widehat{r}_{k}\|_{2}$
$\displaystyle\leq\varepsilon+\tilde{\gamma}_{n}^{f}\left\||e_{k}|+|A^{T}||\widehat{m}_{k}|\right\|_{2}$
$\displaystyle\leq 2\varepsilon,$

from which we can write the Frobenius norm bound

$\|I-A^{T}\widehat{M}\|_{F}\leq 2\sqrt{n}\varepsilon.$

Note that assuming exact arithmetic, e.g., using $u_{f}=0$ in (8), the above
bound becomes $\|I-A^{T}M\|_{F}\leq\sqrt{n}\varepsilon$, which is the same
bound given by Grote and Huckle [23, Theorem 3.1]. At this point we can
transpose to get the desired bound for the left preconditioner
$\widehat{P}=\widehat{M}^{T}$, and then convert to the desired infinity norm,
giving

(12) $\|I-\widehat{P}A\|_{\infty}\leq 2n\varepsilon.$

We note that from (12) we then also have an a priori bound on the distance
from the computed $\widehat{P}$ to the true inverse of $A$,

(13)
$\|\widehat{P}-A^{-1}\|_{\infty}=\|(\widehat{P}A-I)A^{-1}\|_{\infty}\leq\|I-\widehat{P}A\|_{\infty}\|A^{-1}\|_{\infty}\leq
2n\varepsilon\|A^{-1}\|_{\infty}.$

An important conclusion that we can draw from this analysis, in particular
from the fact that (12) is guaranteed as long as (11) holds, is that there is
_no advantage_ to using a higher precision in computing the SPAI
preconditioner beyond that dictated by
$u_{f}\approx\varepsilon\text{cond}_{2}^{-1}(A^{T})$. Using higher precision
will have very little effect, since we will end up with a preconditioner of
similar quality. We will demonstrate this observation in Section 4 below. Note
that this is in contrast to the case of using full LU factorization, where the
quality of the factors scales with the precision used.

Intuitively, however, we may think that using higher precision for $u_{f}$
might result in a $\widehat{P}$ with fewer (more accurate) nonzeros than a
$\widehat{P}$ computed in lower precision. This can be the case, but not
necessarily. We demonstrate this in Figure 1 for the matrices saylr1 (left
plot) and steam3 (right plot) from SuiteSparse [16], where we plot the number
of nonzeros in the resulting $\widehat{P}$ computed using the SPAI algorithm
(Algorithm 1) in single precision (blue squares) and double precision (red
circles) for various $\varepsilon$ values. Within Algorithm 1, we use the
identity as the initial sparsity pattern, $\beta=8$, and $\alpha=\lceil
n/\beta\rceil$; this value of $\alpha$ allows the sparse approximate inverse
to fill in as much as needed in order to satisfy (7).

For saylr1, there is very little relative difference between single and double
precision in the resulting size of $\widehat{P}$; if anything, single
precision results in a $\widehat{P}$ with slightly fewer nonzeros. For the
matrix steam3, the relative difference between single and double precision is
more significant; $\widehat{P}$ computed in double precision has, on average,
about 24% fewer nonzeros than the $\widehat{P}$ computed in single precision.
Note that for both matrices, these choices of $\varepsilon$ satisfy (11) for
both single and double precision; see Table 3 for properties of the matrices.
We will explore further examples in Section 4.3.

Figure 1: Number of nonzeros in the computed sparse approximate inverse
preconditioner $\widehat{P}$ when $\widehat{P}$ is computed in single
precision (blue squares) and double precision (red circles), for the matrices
saylr1 (left) and steam3 (right), for various choices of $\varepsilon$.

###  3.1 GMRES-based iterative refinement with SPAI

Our goal is now to bound the relative error in the correction solve when the
solver is GMRES left-preconditioned by the approximate inverse $\widehat{P}$.
For clarity, we present this variant, which we call SPAI-GMRES-IR, in
Algorithm 3. We follow and heavily rely on the analysis of [2], [12], and
[11]. Here we will analyze the general five-precision GMRES-IR variant of [2],
in which in addition to precisions $u_{f}$, $u$, and $u_{r}$, $u_{g}$ denotes
the working precision used within GMRES and $u_{p}$ represents the precision
in which the preconditioner or preconditioned linear system is applied to a
vector within GMRES.

Algorithm 3 GMRES-based Iterative Refinement with SPAI Preconditioning (SPAI-
GMRES-IR)

1:$n\times n$ nonsingular matrix $A$ and length-$n$ right-hand side vector
$b$, maximum number of refinement steps $i_{max}$, GMRES convergenge tolerance
$\tau$, SPAI parameter $\varepsilon$.

2:Approximate solution $x_{i+1}$ to $Ax=b$.

3:Run Algorithm 1 on $A^{T}$ with parameter $\varepsilon$ to obtain $M$; set
$P=M^{T}$.

4:Compute $x_{0}=Pb$ in precision $u_{f}$; store $x_{0}$ in precision $u$.

5:for $i=0:i_{max}-1$ do

6: Compute $r_{i}=b-Ax_{i}$ in precision $u_{r}$; store in precision $u$.

7: Solve $PAd_{i}=Pr_{i}$ via left-preconditioned GMRES with tolerance $\tau$
in working precision $u_{g}$, with matrix-vector products with $P$ and $A$
computed in precision $u_{p}$; store $d_{i}$ in precision $u$.

8: Update $x_{i+1}=x_{i}+d_{i}$ in precision $u$.

9:end for

We first want to bound $\kappa_{\infty}(\tilde{A})$, where
$\tilde{A}=\widehat{P}A$. We stress that the condition number of the
preconditioned matrix will not tell us anything about the convergence rate of
GMRES (see [22]); here we will use it to bound the relative error after we
obtain a backward error result for the preconditioned system. We can write

$\displaystyle\tilde{A}=\widehat{P}A$ $\displaystyle=I-(I-\widehat{P}A),$
$\displaystyle\tilde{A}^{-1}=A^{-1}\widehat{P}^{-1}$
$\displaystyle=(I-(I-\widehat{P}A))^{-1}$ (14) $\displaystyle\approx
I+(I-\widehat{P}A).$

Thus using (12), we have

$\displaystyle\|\tilde{A}\|_{\infty}$ $\displaystyle\leq 1+2n\varepsilon,$
$\displaystyle\|\tilde{A}^{-1}\|_{\infty}$ $\displaystyle\lesssim
1+2n\varepsilon,$

which gives

(15) $\kappa_{\infty}(\tilde{A})\lesssim(1+2n\varepsilon)^{2}.$

We note that this should not be taken as an upper bound but rather as a rough
estimate, since for practical problem sizes, $2n\varepsilon$ will likely be
greater than 1, and thus our dropping of higher order terms in the Taylor
expansion to get (14) may be dubious. Although (15) is not strictly an upper
bound, we note that the growth of $\kappa_{\infty}(\tilde{A})$ does generally
follow this bound asymptotically. In Figure 2, we plot the condition number
$\kappa_{\infty}(\tilde{A})$ versus $\varepsilon$ for two matrices from
SuiteSparse [16], saylr1 (Fig. 2(a)) and steam3 (Fig. 2(b)), where the sparse
approximate inverses are computed in single precision (blue squares) and
single precision (red circles). As before, we use the identity as the initial
sparsity pattern, $\beta=8$, and $\alpha=\lceil n/\beta\rceil$; see Algorithm
1. We additionally plot $(1+2n\varepsilon)^{2}$, the estimate in (15),
demonstrating that (15) well describes the resulting condition number of the
preconditioned coefficient matrix. We again note that there is little
difference between single and double precision, since we expect preconditioner
quality to depend on $\varepsilon$ rather than $u_{f}$ as long as (11) is
satisfied.

(a) saylr1

(b) steam3

Figure 2: $\kappa_{\infty}(\tilde{A})$ versus $\varepsilon$ for matrices
saylr1 (left) and steam3 (right) for SPAI preconditioners computed in single
(blue squares) and double (red circles), compared with the approximation (15).

Now we turn our attention to the application of the preconditioner to the
right-hand side $\widehat{r}_{i}$. Let the (exact) preconditioned right-hand
side be $z_{i}=\widehat{P}\widehat{r}_{i}$. Assuming we apply $\widehat{P}$ to
$\widehat{r}_{i}$ in precision $u_{p}$, we have

$\displaystyle\widehat{z}_{i}$ $\displaystyle=(\widehat{P}+\Delta
P)\widehat{r}_{i},\qquad|\Delta P|\leq\gamma_{n}^{p}|\widehat{P}|,$
$\displaystyle=z_{i}+\Delta P\widehat{P}^{-1}z_{i}.$

We then have the bound

(16)
$\|\widehat{z}_{i}-z_{i}\|_{\infty}\leq\gamma_{n}^{p}\||\widehat{P}||\widehat{P}^{-1}|\|_{\infty}\|z_{i}\|_{\infty}\leq\gamma_{n}^{p}\kappa_{\infty}(\widehat{P})\|z_{i}\|_{\infty}\lesssim\gamma_{n}^{p}\kappa_{\infty}(A)\|z_{i}\|_{\infty},$

under the reasonable assumption that
$\kappa_{\infty}(\widehat{P})\approx\kappa_{\infty}(A)$.

Following [11] and [2], we now seek to bound the error in applying
$\widehat{P}A$ to a vector $\widehat{v}_{j}$ in iteration $j$ of the GMRES
algorithm. We compute this via

$\displaystyle(A+\delta A)\widehat{v}_{j}$
$\displaystyle=\widehat{w}_{j},\qquad|\Delta A|\leq\gamma_{n}^{p}|A|,$
$\displaystyle(\widehat{P}+\Delta\widehat{P})\widehat{w}_{j}$
$\displaystyle=\widehat{y}_{j},\qquad|\Delta\widehat{P}|\leq\gamma_{n}^{p}|\widehat{P}|.$

Then we can write the computed vector $\widehat{y}_{j}$ as

$\widehat{y}_{j}=(\widehat{P}+\Delta\widehat{P})(A+\Delta
A)\widehat{v}_{j}=(\widehat{P}A+\Delta\widehat{P}A+\widehat{P}\Delta
A)\widehat{v}_{j}=(\tilde{A}+\Delta\tilde{A}^{\prime})\widehat{v}_{j},$

where

$\displaystyle\Delta\tilde{A}^{\prime}$
$\displaystyle=\Delta\widehat{P}A+\widehat{P}\Delta A$
$\displaystyle=\Delta\widehat{P}\widehat{P}^{-1}\widehat{P}A+\widehat{P}AA^{-1}\Delta
A$
$\displaystyle=\Delta\widehat{P}\widehat{P}^{-1}\tilde{A}+\tilde{A}A^{-1}\Delta
A,$

for which we can write the bound

$\displaystyle\|\Delta\tilde{A}^{\prime}\|_{F}$
$\displaystyle\leq\gamma_{n}^{p}\big{(}\||\widehat{P}||\widehat{P}^{-1}|\|_{F}+\||A^{-1}||A|\|_{F}\big{)}\|\tilde{A}\|_{F}$
$\displaystyle\leq\gamma_{n}^{p}\big{(}\kappa_{F}(\widehat{P})+\kappa_{F}(A)\big{)}\|\tilde{A}\|_{F}$
$\displaystyle\lesssim\tilde{\gamma}_{n}^{p}\kappa_{F}(A)\|\tilde{A}\|_{F}.$

Thus

$\widehat{y}_{j}=\tilde{A}\widehat{v}_{j}+f_{j},$

with

$\|f_{j}\|_{\infty}\leq\tilde{\gamma}_{n}^{p}\kappa_{F}(A)\|\tilde{A}\|_{F}\|\widehat{v}_{j}\|_{2}.$

Since [2, Eqn. 3.2] is satisfied, we can thus apply Theorem 3.1 from [2],
which gives the result that for some iteration $k\leq n$, the left-
preconditioned MGS-GMRES algorithm will produce an approximate solution
$\widehat{d}_{i}$ for which

$(\tilde{A}+\Delta\tilde{A})\widehat{d}_{i}=\widehat{z}_{i}+\Delta\widehat{z}_{i},$

with

(17)
$\|\Delta\tilde{A}\|_{F}\leq(\tilde{\gamma}_{kn}^{g}+n\tilde{\gamma}_{k^{1/2}n}^{p}\kappa_{\infty}(A))\|\tilde{A}\|_{F},\quad\|\Delta\widehat{z}_{i}\|_{2}\leq\tilde{\gamma}_{kn}^{g}\|\widehat{z}_{i}\|_{2}\lesssim
n^{1/2}\tilde{\gamma}_{kn}^{g}\|z_{i}\|_{\infty}.$

Using this result together with (16), we can bound the normwise relative
backward error by

(18)
$\frac{\|z_{i}-\tilde{A}\widehat{d}_{i}\|_{\infty}}{\|\tilde{A}\|_{\infty}\|\widehat{d}_{i}\|_{\infty}+\|z_{i}\|_{\infty}}\lesssim
f(n,k)\left(u_{g}+u_{p}\kappa_{\infty}(A)\right),$

where $f(n,k)=n^{3}k^{1/2}$. This gives the relative forward error bound

(19) $\frac{\|\widehat{d}_{i}-d_{i}\|_{\infty}}{\|d_{i}\|_{\infty}}\lesssim
f(n,k)\left(u_{g}+u_{p}\kappa_{\infty}(A)\right)\kappa_{\infty}(\tilde{A}).$

Note that from (19), we have that for SPAI-GMRES-IR, the quantity
$u_{s}\|E_{i}\|_{\infty}$ in (3) is on the order
$f(n,k)\left(u_{g}+u_{p}\kappa_{\infty}(A)\right)\kappa_{\infty}(\tilde{A})$,
and by (17), the quantities $c_{1}$ and $c_{2}$ in (4) will be similar to
those stated for five-precision GMRES-IR in [2]. The difference will lie in
the values of $\kappa_{\infty}(\tilde{A})$ and $\|\tilde{A}\|_{F}$, which here
will depend on the $\varepsilon$ used in constructing the SPAI preconditioner.

Using the estimate for the condition number of the preconditioned coefficient
matrix in (15), we can thus give the rough bound

(20) $\frac{\|\widehat{d}_{i}-d_{i}\|_{\infty}}{\|d_{i}\|_{\infty}}\lesssim
f(n,k)\left(u_{g}+u_{p}\kappa_{\infty}(A)\right)(1+2n\varepsilon)^{2}.$

Assume we are in the setting where $u=u_{p}=u_{g}$, which is commonly used in
practice. Combining (20) with (11), we thus must have roughly

$u_{f}\text{cond}_{2}(A^{T})\lesssim\varepsilon\lesssim
u^{-1/2}\kappa_{\infty}(A)^{-1/2}$

in order to guarantee both that the SPAI algorithm will complete and that the
GMRES-based iterative refinement scheme will converge.

We note that the bound (18) on the normwise relative backward error that can
be achieved with GMRES in precisions $u_{g}$ and $u_{p}$ assumes that we do
not stop the iterations early. In practice, we will use some stopping
criterion $\tau$ within GMRES. If $\tau$ is based on the backward error of the
preconditioned system, then (18) becomes

(21)
$\frac{\|z_{i}-\tilde{A}\widehat{d}_{i}\|_{\infty}}{\|\tilde{A}\|_{\infty}\|\widehat{d}_{i}\|_{\infty}+\|z_{i}\|_{\infty}}\lesssim\max\left(f(n,k)\left(u_{g}+u_{p}\kappa_{\infty}(A)\right),\tau\right).$

In theory, we should then set $\tau$ to be on the same order as the other term
in the maximum above if we do not want to affect the convergence conditions.
In practice, however, the convergence conditions are often too strict, and
thus we can use a larger $\tau$ and still maintain convergence of the
iterative refinement process, perhaps with a larger $\phi_{i}$. In practice,
choosing $\tau$ is largely done based on empirical observations, and is a
tradeoff between the number of GMRES iterations per refinement step and the
overall number of refinement steps; see, e.g., the discussion and experiments
in [2, Section 4.2].

Comparing the expected behavior of SPAI-GMRES-IR and GMRES-IR with full LU
factorization, we expect that convergence for GMRES with the SPAI
preconditioner will in general be slower than GMRES using the full LU factors,
at least when we assume common choices of $\varepsilon$. However, the
construction and application of the SPAI preconditioner may be less expensive
and more parallelizable than for the full LU factors. This behavior and the
resulting tradeoffs will be strongly problem dependent.

##  4 Numerical Experiments

We now turn to an experimental evaluation of the numerical behavior of SPAI-
GMRES-IR. Our experiments are performed in MATLAB R2020a. In these experiments
we use four precisions. For half precision, we use the chop library available
at https://github.com/higham/chop and the associated functions available at
https://github.com/SrikaraPranesh/Multi_precision_NLA_kernels. For quadruple
precision, we use the Advanpix Multiprecision Computing Toolbox [37]. For
single and double precision we use the built-in MATLAB datatypes. The code
used to generate plots in this work can be found at
https://github.com/Noaman67khan/SPAI-GMRES-IR.

The matrices tested come from the SuiteSparse Matrix Collection [16]. We list
tested matrices and their relevant properties in Table 3. In all tests, the
right-hand sides are set to the vector with equal components and unit
$2$-norm. For the convergence tolerance in GMRES, we use $\tau=10^{-4}$ when
the working precision is single and $\tau=10^{-8}$ when the working precision
is double, which responds to roughly the square root of the working precision.
These are the default values which have been used in previous works, e.g.,
[2], [12], although as mentioned in Section 3.1, this parameter must be tuned
in practical applications. In Table 3, we split the matrices into two sets,
those used with precisions $(u_{f},u,u_{r})=$ (single, double, quad), and
those used with precisions $(u_{f},u,u_{r})=$ (half, single double), based on
whether GMRES-IR with LU preconditioning converges with the given precisions
and value of $\tau$.

In all cases, we set the initial sparsity pattern for SPAI to be the identity
matrix. As noted by Sedlacek [39, Section 3.1.2], this can result in SPAI
preconditioners with zero columns in the case that $A$ has a zero entry on the
diagonal. We therefore restrict ourselves to problems with nonzero diagonal
entries, but note that this could be remedied by either permuting $A$ or using
the initial sparsity pattern of $A$, which, when SPAI is run on $A^{T}$,
guarantees that we obtain a $P$ with nonzero rows [39, Theorem 3.1].

For SPAI-GMRES-IR, we preprocess the matrices via a column scaling such that
the largest element in absolute value in every column of $A^{T}$ is 1. A one-
sided scaling was proposed in the context of using low-precision QR
factorization to solve least squares problems in [13, Section 3.2], which is
precisely the task we perform in SPAI construction. The goal of this scaling
is to avoid overflow in the low precision computations. To be precise, we run
SPAI on $A^{T}D$ to obtain $M$, where $D$ is the diagonal matrix that performs
the scaling, and then set $P=M^{T}D$. For all tests, we use $\beta=8$, which
is in the range suggested by Sedlacek [39]. Since we want to allow the
approximate inverse to fill in as much as necessary so that (7) is satisfied,
we allow $\alpha$ to be as large as $\lceil n/\beta\rceil$.

As a point of comparison, we also provide results for GMRES-IR with full LU
factorization (which in this section we refer to as LU-GMRES-IR) as well as
GMRES-IR with no preconditioner. To enable a fair comparison, for LU-GMRES-IR,
we first test multiple reorderings of the matrix and then choose the one that
provides the least number of nonzeros in $L$ and $U$. The reorderings tested
are reverse Cuthill-McKee (rcm), approximate minimum degree (amd), nested
dissection (nds), and column count (col). For LU-GMRES-IR, if $u_{f}=$ half
produces LU factors containing NAN or INF, we implement the scaling advised in
[28].

We note that for SPAI-GMRES-IR, we use the natural ordering; as noted in [10,
p. 12], there is no expected benefit to reordering for the case of
nonfactorized sparse approximate inverses. We also note that in all variants
of GMRES-IR tested here, including SPAI-GMRES-IR, we use the variant in which
$u_{g}=u_{p}=u$, which is commonly used in practice.

Table 3: Matrices used in experiments along with their properties. All matrices come from the SuiteSparse collection [16].  Name | $n$ | $nnz$ | $nnz(A^{-1})$ | $\kappa_{\infty}(A)$ | $\text{cond}_{2}(A^{T})$
---|---|---|---|---|---
pores_3 | 532 | 3474 | 213712 | $1.2e+06$ | $1.7e+05$
steam1 | 240 | 2248 | 57599 | $3.1e+07$ | $2.8e+03$
steam3 | 80 | 314 | 6315 | $7.6e+10$ | $5.6e+03$
saylr1 | 238 | 1128 | 56644 | $1.6e+09$ | $5.2e+05$
bfwa782 | 782 | 7514 | 458839 | $6.8e+03$ | $1.3e+03$
cage5 | 37 | 233 | 1369 | $2.9e+01$ | $7.5e+00$
gre_115 | 115 | 421 | 13225 | $1.4e+02$ | $3.7e+01$
orsreg_1 | 2205 | 14133 | 4862025 | $7.0e+03$ | $5.9e+03$
sherman4 | 1104 | 3786 | 298674 | $3.1e+03$ | $1.2e+03$

For each linear system and given combination of precisions, we run SPAI-GMRES-
IR with various values of $\varepsilon$, LU-GMRES-IR, and GMRES-IR with no
preconditioning. We note that it is incredibly difficult to give an accurate
comparison of cost among the different methods; the relative cost of a GMRES
iteration in GMRES-IR without preconditioning and SPAI-GMRES-IR and LU-GMRES-
IR will be different, since the (preconditioned) coefficient matrix has
different structure and density. SPAI-based preconditioners may also be more
inexpensive to apply (and more parallelizable), since their application
involves only a single SpMV rather than two triangular solves. Here we will
give a qualitative comparison of approaches in terms of preconditioner size
and number of GMRES iterations, but we note that high-performance experiments
on large-scale matrices are necessary future work in order to give a
quantitative comparison.

The tables for each setup give the size of the preconditioner in terms of
number of nonzeros, the resulting condition number of the preconditioned
system, as well as information about the number of GMRES-IR refinement steps
and GMRES iterations per refinement step. For LU-GMRES-IR, the size of the
preconditioner is given as $nnz(L+U)$, where $\Pi(QAQ^{T})=LU$ where $Q$ is
the best reordering among those tested and $\Pi$ is a permutation arising from
partial pivoting. The best reorderings are noted in the tables. We stress that
quantitative conclusions about the relative performance of LU-GMRES-IR versus
SPAI-GMRES-IR cannot be drawn from the size of the preconditioner.

In the final column in the tables, the first number gives the total number of
GMRES iterations summed over all refinement steps, and the following tuple
gives in position $i$ the number of GMRES iterations in refinement step $i$.
Note that the number of refinement steps is given by the number of elements in
the tuple.

The best case scenario is that SPAI preconditioning produces a preconditioner
that is sparser than the full LU factorization and still effective in reducing
the total number of GMRES iterations versus the unpreconditioned case. We
again stress, however, that a full picture of expected relative performance
cannot be given by the size of the preconditioner and the number of
iterations. In Section 4.1, we give examples using precisions
$(u_{f},u,u_{r})=$ (single, double, quad) and in Section 4.2, we use
$(u_{f},u,u_{r})=$ (half, single, double). Finally, in Section 4.3, we redo
the experiments in Section 4.2 but now using $(u_{f},u,u_{r})=$ (single,
single, double). These results demonstrate that there is often not a
significant advantage, either in terms of the size of the preconditioner or
the total number of GMRES iterations, to using $u_{f}$ higher than that
dictated by the constraint (11). That is, the resulting SPAI preconditioner
will satisfy (12) regardless of whether half or single precision is used, and
thus we do not expect a significant difference in the behavior of the
iterative refinement process.

###  4.1 Experiments with $(u_{f},u,u_{r})=$ (single, double, quad)

Experiments using the settings $(u_{f},u,u_{r})=$ (single, double, quad) are
shown in Table 4. First, we observe that it is clear that SPAI-GMRES-IR
converges in every case. The potential of SPAI-GMRES-IR is perhaps best
exhibited for the matrix steam1. Here, even with reordering the computed LU
factors contain 14,133 nonzeros, compared to the SPAI preconditioner, which
contains 1,140 nonzeros when $\varepsilon=0.2$ and $1,303$ nonzeros when
$\varepsilon=0.1$ (note that this is fewer than the number of nonzeros in
$A$). In both cases, SPAI-GMRES-IR requires only 14 total GMRES iterations to
converge to double precision accuracy; while this more than required by LU-
GMRES-IR, it is still significantly fewer than the case where no
preconditioner is used, which requires 543 total GMRES iterations.

For the matrices pores_3, steam3, and saylr1, the results are as expected. The
SPAI preconditioner, depending on the $\varepsilon$ parameter used, can have
significantly fewer nonzeros than the LU factors, and the total number of
GMRES iterations required for SPAI-GMRES-IR is somewhere in between that of
LU-GMRES-IR and GMRES-IR with no preconditioning. We again stress that it is
difficult to make quantitative comparisons; although the SPAI preconditioner
may have close to the same number of nonzeros as the LU factors, it may be
less expensive to apply and more parallelizable. We also note that the
tradeoff in choosing $\varepsilon$ is clear: a larger $\varepsilon$ can lead
to a sparser preconditioner, but a greater total number of GMRES iterations
required for the iterative refinement process to converge.

Table 4: Comparison of SPAI-GMRES-IR for different $\varepsilon$ values with LU-GMRES-IR and GMRES-IR with no preconditioner for test matrices using $(u_{f},u,u_{r})=$ (single, double, quad). | Preconditioner | $\kappa_{\infty}(\tilde{A})$ |  Precond. $nnz$ |  | GMRES-IR steps/iterations
---
pores_3 |  SPAI, $\varepsilon=0.5$ | $6.6e+03$ | 3560 | 223(110, 113)
SPAI, $\varepsilon=0.4$ | $3.8e+03$ | 4871 | 174(86, 88)
Full LU (amd) | $1.0e+00$ | 9706 | 4(2, 2)
None | $1.2e+06$ | 0 | 1314(417, 456, 441)
steam1 |  SPAI, $\varepsilon=0.2$ | $1.5e+00$ | 1140 | 14(7, 7)
SPAI, $\varepsilon=0.1$ | $1.5e+00$ | 1303 | 14(7, 7)
Full LU (amd) | $1.9e+00$ | 14133 | 2(2)
None | $3.1e+07$ | 0 | 543(158, 193, 192)
steam3 |  SPAI, $\varepsilon=0.5$ | $3.9e+00$ | 244 | 31(9, 12, 10)
SPAI, $\varepsilon=0.1$ | $1.9e+00$ | 403 | 17(5, 6, 6)
Full LU (amd) | $1.1e+00$ | 483 | 2(2)
None | $7.6e+10$ | 0 | 221(61, 80, 80)
saylr1 |  SPAI, $\varepsilon=0.4$ | $1.9e+04$ | 1932 | 195(64, 66, 65)
SPAI, $\varepsilon=0.3$ | $7.5e+03$ | 3405 | 89(44, 45)
Full LU (amd) | $1.0e+00$ | 3607 | 5(2, 3)
None | $1.6e+09$ | 0 | 658(214, 229, 215)

###  4.2 Experiments with $(u_{f},u,u_{r})=$ (half, single, double)

We now turn to examples using precisions $(u_{f},u,u_{r})=$ (half, single,
double), the results for which are shown in Table 5. We again note that we
have limited ourselves to examples where the use of half precision does not
cause over/underflow in the SPAI construction despite the scaling used.

From Table 5, we can see that SPAI-GMRES-IR indeed converges in every case,
although again we can see that there are clear tradeoffs related to the choice
of $\varepsilon$. When we choose a larger $\varepsilon$, the preconditioner
becomes sparser, but also less effective, requiring a greater total number of
GMRES iterations to converge. For example, for the linear system with the
matrix sherman4, $\varepsilon=0.5$ leads to a very sparse preconditioner
(nearly diagonal), but requires almost as many GMRES iterations across all
refinement steps as the unpreconditioned case.

For the problem with matrix gre_115, we can see that for the choice of the
smaller $\varepsilon$, $\varepsilon=0.3$, the resulting SPAI preconditioner
has a greater number of nonzeros than the LU preconditioner, and is also less
effective in terms of reducing the total number of GMRES iterations.

The matrix orsreg_1 perhaps best exhibits the potential of SPAI-GMRES-IR.
Here, with $\varepsilon=0.3$, the SPAI preconditioner has over $12\times$
fewer nonzeros than the LU factors and still reasonably improves the
convergence rate ($5.8\times$ fewer GMRES iterations than the unpreconditioned
case).

Table 5: Comparison of SPAI-GMRES-IR for different $\varepsilon$ values with LU-GMRES-IR and GMRES-IR with no preconditioner for test matrices using $(u_{f},u,u_{r})=$ (half, single, double). | Preconditioner | $\kappa_{\infty}(\tilde{A})$ |  Precond. $nnz$ |  | GMRES-IR steps/iterations
---
bfwa782 |  SPAI, $\varepsilon=0.5$ | $1.1e+03$ | 6271 | 164(75, 89)
SPAI, $\varepsilon=0.3$ | $5.0e+02$ | 11430 | 114(54, 60)
Full LU (amd) | $2.1e+00$ | 21838 | 7(3, 4)
None | $6.8e+03$ | 0 | 381(172, 209)
cage5 |  SPAI, $\varepsilon=0.5$ | $9.9e+00$ | 101 | 16(8, 8)
SPAI, $\varepsilon=0.3$ | $3.9e+00$ | 213 | 12(6, 6)
Full LU (amd) | $1.0e+00$ | 359 | 2(2)
None | $2.9e+01$ | 0 | 25(13, 12)
gre_115 |  SPAI, $\varepsilon=0.5$ | $5.8e+02$ | 725 | 48(24, 24)
SPAI, $\varepsilon=0.3$ | $1.8e+01$ | 1719 | 21(10, 11)
Full LU (nds) | $1.0e+00$ | 1551 | 2(2)
None | $1.4e+02$ | 0 | 100(49, 51)
orsreg_1 |  SPAI, $\varepsilon=0.5$ | $1.7e+02$ | 9261 | 108(29, 45, 34)
SPAI, $\varepsilon=0.3$ | $1.3e+02$ | 11120 | 61(23, 38)
Full LU (rcm) | $2.2e+00$ | 133634 | 9(4, 5)
None | $7.0e+03$ | 0 | 352(107, 150, 95)
sherman4 |  SPAI, $\varepsilon=0.5$ | $1.6e+03$ | 1386 | 140(67, 73)
SPAI, $\varepsilon=0.3$ | $5.0e+02$ | 8496 | 74(35, 39)
Full LU (amd) | $1.8e+00$ | 14211 | 5(2, 3)
None | $3.1e+03$ | 0 | 178(85, 93)

###  4.3  Experiments with $(u_{f},u,u_{r})=$ (single, single, double)

We now test the same problems used in Section 4.2, now using single precision
instead of half precision for $u_{f}$. The results are presented in Table 6.
The analysis in Section 3 indicates that as long as, say, half precision is
sufficient for producing a preconditioner satisfying (11), there is in general
no expected benefit to using higher precision in the SPAI construction in
terms of preconditioner quality. We do recall, however, that the number of
nonzeros in the preconditioner may be different when higher versus lower
precision is used. For easy comparison, in Table 6, we give in red bracketed
format the corresponding condition number, preconditioner size, and total
number of GMRES iterations from Table 5.

For all test problems, notice that the results for SPAI-GMRES-IR either do not
change at all or change only very slightly. In cases where there is a slight
difference, note that it is not necessarily clear that $u_{f}$ = single or
$u_{f}=$ half is better. For example, for the problem with gre_115 and
$\varepsilon=0.5$, $u_{f}=$ half leads to 48 total GMRES iterations whereas
$u_{f}=$ single leads to 51. Conversely, for the problem with orsreg_1 and
$\varepsilon=0.5$, $u_{f}=$ half leads to 108 total GMRES iterations whereas
$u_{f}=$ single leads to 97. We note that for both these cases, the SPAI
preconditioner has the same number of nonzeros regardless of whether $u_{f}=$
half or $u_{f}=$ single is used.

It is also not necessarily clear that $u_{f}$ = single or $u_{f}=$ half is
better in terms of the size of the preconditioner. For example, for the
problem with bfwa782 and $\varepsilon=0.3$, the SPAI preconditioner has 11,430
nonzeros when $u_{f}=$ half, but 11,470 when $u_{f}=$ single. For the problem
with orsreg_1 and $\varepsilon=0.3$, the SPAI preconditioner has 11,120
nonzeros when $u_{f}=$ half, and 11,025 when $u_{f}=$ single.

In any case, we can see that this data supports our theoretical conclusions:
in general, there is no expected benefit to using a higher precision than that
dictated by (11) for SPAI in terms of preconditioner quality. Note that this
is in sharp contrast to the case of LU preconditioning, where the precision
used for the LU factorization makes a significant difference in preconditioner
quality. This intuitively makes sense since LU factorization is a direct
algorithm rather than an iterative one like SPAI that iterates until the
stopping criterion is met.

One interesting anomaly that we point out is that from Table 6, we can see
that the LU factors can have significantly fewer nonzeros when we use $u_{f}=$
half versus $u_{f}=$ single. This difference is most evident for the problem
with matrix orsreg_1, where $u_{f}=$ half results in preconditioners with
133,634 nonzeros and $u_{f}=$ single results in preconditioners with 330,910
nonzeros. We expect this is due to underflow; the true LU factors likely
contain many small elements which become zeros when half precision is used.

Table 6: Comparison of SPAI-GMRES-IR for different $\varepsilon$ values with LU-GMRES-IR and GMRES-IR with no preconditioner for test matrices using $(u_{f},u,u_{r})=$ (single, single, double). | Preconditioner | $\kappa_{\infty}(\tilde{A})$ |  Precond. $nnz$ |  | GMRES-IR steps/iterations
---
bfwa782 |  SPAI, $\varepsilon=0.5$ |  $1.1e+03$ [$1.1e+03$] |  6261 [6271] |  166(74, 92) [164]
SPAI, $\varepsilon=0.3$ |  $5.0e+02$ [$5.0e+02$] |  11470 [11430] |  114(54, 60) [114]
Full LU (amd) |  $1.0e+00$ [$2.1e+00$] |  21848 [21838] |  1(1) [7]
None | $6.8e+03$ | 0 | 381(172, 209)
cage5 |  SPAI, $\varepsilon=0.5$ |  $9.9e+00$ [$9.9e+00$] |  101 [101] |  16(8, 8) [16]
SPAI, $\varepsilon=0.3$ |  $3.9e+00$ [$3.9e+00$] |  213 [213] |  12(6, 6) [12]
Full LU (amd) |  $1.0e+00$ [$1.0e+00$] |  359 [359] |  1(1) [2]
None | $2.9e+01$ | 0 | 25(13, 12)
gre_115 |  SPAI, $\varepsilon=0.5$ |  $6.2e+02$ [$5.8e+02$] |  725 [725] |  51(24, 27) [48]
SPAI, $\varepsilon=0.3$ |  $1.7e+01$ [$1.8e+01$] |  1739 [1719] |  20(10, 10) [21]
Full LU (nds) |  $1.0e+00$ [$1.0e+00$] |  1556 [1551] |  1(1) [2]
None | $1.4e+02$ | 0 | 100(49, 51)
orsreg_1 |  SPAI, $\varepsilon=0.5$ |  $1.4e+02$ [$1.7e+02$] |  9261 [9261] |  97(25, 40, 32) [108]
SPAI, $\varepsilon=0.3$ |  $1.1e+02$ [$1.3e+02$] |  11025 [11120] |  60(22, 38) [61]
Full LU (rcm) |  $1.0e+00$ [$2.2e+00$] |  330910 [133634] |  1(1) [9]
None | $7.0e+03$ | 0 | 352(107, 150, 95)
sherman4 |  SPAI, $\varepsilon=0.5$ |  $1.6e+03$ [$1.6e+03$] |  1385 [1386] |  140(67, 73) [140]
SPAI, $\varepsilon=0.3$ |  $5.0e+02$ [$5.0e+02$] |  8499 [8496] |  74(35, 39)[74]
Full LU (amd) |  $1.0e+00$ [$1.8e+00$] |  14211 [14211] |  1(1) [5]
None | $3.1e+03$ | 0 | 178(85, 93)

##  5 Conclusions and future work

In this work we explored the use of sparse approximate inverse preconditioners
within mixed precision GMRES-based iterative refinement as a method of solving
sparse linear systems.

We first analyzed the computation of an SPAI preconditioner in some finite
precision $u_{f}$. The stopping criterion used in constructing each column of
the sparse approximate inverse $M$ is based on the size of the measured
residual norm $e_{k}-A^{T}m_{k}$ being less than some user-specified parameter
$\varepsilon$. Thus for a given $\varepsilon$, our analysis shows that we must
use $u_{f}\lesssim\varepsilon\text{cond}_{2}^{-1}(A^{T})$ in order to
guarantee that this stopping criterion can be achieved. As long as this
constraint on $u_{f}$ holds, then the computed $\widehat{M}$ satisfies the
bounds proved for the exact arithmetic case in [23] up to a constant factor.
An interesting point is that there is no significant benefit in terms of
resulting preconditioner quality to using a precision higher than dictated by
this level.

We then prove conditions under which GMRES left preconditioned by the SPAI
preconditioner will be backward stable with respect to the preconditioned
system. As long as the constraint
$u_{f}\lesssim\varepsilon\text{cond}_{2}^{-1}(A^{T})$ is satisfied, meaning
that an approximate preconditioner with parameter $\varepsilon$ is computable
in precision $u_{f}$, then the constraints for convergence of five-precision
SPAI-GMRES-IR are essentially the same as those for the five-precision GMRES-
IR given in [2]. Our numerical experiments confirm the expected behavior of
SPAI-GMRES-IR.

While we have started with one variant of sparse approximate inverses here, it
is clear that other approximate preconditioners, such as incomplete LU
factorization, factorized sparse approximate inverses, or sparse approximate
inverses based on a fixed sparsity pattern may be suitable to use within mixed
precision Krylov-subspace-based iterative refinement schemes. In particular,
we believe that extension of the analysis here to modified SPAI (see, e.g.,
[30]) as well as factorized SPAI (see, e.g., [32], [43]) should be
straightforward.

## References

  * [1] A. Abdelfattah, H. Anzt, J. Dongarra, M. Gates, A. Haidar, J. Kurzak, P. Luszczek, S. Tomov, I. Yamazaki, and A. YarKhan, Linear algebra software for large-scale accelerated multicore computing, Acta Numerica, 25 (2016), pp. 1–160.
  * [2] P. Amestoy, A. Buttari, N. J. Higham, J.-Y. L’Excellent, T. Mary, and B. Vieublé, Five-precision GMRES-based iterative refinement, MIMS EPrint 2021.5, Manchester Institute for Mathematical Sciences, The University of Manchester, Manchester, UK, apr 2021.
  * [3] P. Amestoy, A. Buttari, N. J. Higham, J.-Y. L’Excellent, T. Mary, and B. Vieublé, Combining sparse approximate factorizations with mixed precision iterative refinement.  working paper or preprint, Jan. 2022, https://hal.archives-ouvertes.fr/hal-03536031.
  * [4] H. Anzt, T. K. Huckle, J. Bräckle, and J. Dongarra, Incomplete sparse approximate inverses for parallel preconditioning, Par. Comput., 71 (2018), pp. 1–22.
  * [5] M. Arioli and I. S. Duff, Using FGMRES to obtain backward stability in mixed precision, Elec. Trans. Numer. Anal., 33 (2009), pp. 31–44.
  * [6] S. T. Barnard, L. M. Bernardo, and H. D. Simon, An MPI implementation of the SPAI preconditioner on the T3E, Int. J. High Perf. Comput. Appl., 13 (1999), pp. 107–123.
  * [7] M. Benzi, Preconditioning techniques for large linear systems: a survey, J. Comput. Phys., 182 (2002), pp. 418–477.
  * [8] M. Benzi and M. Tůma, A sparse approximate inverse preconditioner for nonsymmetric linear systems, SIAM J. Sci. Comput., 19 (1998), pp. 968–994.
  * [9] M. Benzi and M. Tůma, A comparative study of sparse approximate inverse preconditioners, Appl. Numer. Math., 30 (1999), pp. 305–340.
  * [10] M. Benzi and M. Tůma, Orderings for factorized sparse approximate inverse preconditioners, SIAM J. Sci. Comput., 21 (2000), pp. 1851–1868.
  * [11] E. Carson and N. J. Higham, A new analysis of iterative refinement and its application to accurate solution of ill-conditioned sparse linear systems, SIAM J. Sci. Comput., 39 (2017), pp. A2834–A2856.
  * [12] E. Carson and N. J. Higham, Accelerating the solution of linear systems by iterative refinement in three precisions, SIAM J. Sci. Comput., 40 (2018), pp. A817–A847.
  * [13] E. Carson, N. J. Higham, and S. Pranesh, Three-precision GMRES-based iterative refinement for least squares problems, SIAM J. Sci. Comput., 42 (2020), pp. A4063–A4083.
  * [14] E. Chow, Parallel implementation and performance characteristics of sparse approximate inverse preconditioners with a priori sparsity patterns, Int. J. High Perform. Comput. Appl., 15 (2001), pp. 10–1177.
  * [15] J. Cosgrove, J. Díaz, and A. Griewank, Approximate inverse preconditionings for sparse linear systems, Int. J. Comput. Math., 44 (1992), pp. 91–110.
  * [16] T. A. Davis and Y. Hu, The university of Florida sparse matrix collection, ACM Trans. Math. Soft., 38 (2011), pp. 1–25.
  * [17] M. M. Dehnavi, D. M. Fernandez, J.-L. Gaudiot, and D. D. Giannacopoulos, Parallel sparse approximate inverse preconditioning on graphic processing units, IEEE Trans. Par. Dist. Sys., 24 (2012), pp. 1852–1862.
  * [18] J. Demmel, L. Grigori, M. Hoemmen, and J. Langou, Communication-optimal parallel and sequential Qr and LU factorizations, SIAM J. Sci. Comput., 34 (2012), pp. A206–A239.
  * [19] J. Gao, Q. Chen, and G. He, A thread-adaptive sparse approximate inverse preconditioning algorithm on multi-GPUs, Par. Comput., 101 (2021), p. 102724.
  * [20] J. Gao, K. Wu, Y. Wang, P. Qi, and G. He, GPU-accelerated preconditioned GMRES method for two-dimensional Maxwell’s equations, Int. J. Comput. Math., 94 (2017), pp. 2122–2144.
  * [21] N. I. Gould and J. A. Scott, Sparse approximate-inverse preconditioners using norm-minimization techniques, SIAM J. Sci. Comput., 19 (1998), pp. 605–625.
  * [22] A. Greenbaum, V. Pták, and Z. Strakoš, Any nonincreasing convergence curve is possible for gmres, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 465–469.
  * [23] M. J. Grote and T. Huckle, Parallel preconditioning with sparse approximate inverses, SIAM J. Sci. Comput., 18 (1997), pp. 838–853.
  * [24] G. He, R. Yin, and J. Gao, An efficient sparse approximate inverse preconditioning algorithm on GPU, Concurrency and Computation: Practice and Experience, 32 (2020), p. e5598.
  * [25] N. J. Higham, Accuracy and Stability of Numerical Algorithms, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second ed., 2002.
  * [26] N. J. Higham and T. Mary, A new preconditioner that exploits low-rank approximations to factorization error, SIAM J. Sci. Comput., 41 (2019), pp. A59–A82.
  * [27] N. J. Higham and T. Mary, Solving block low-rank linear systems by LU factorization is numerically stable, IMA J. Numer. Anal., (2020), pp. 1–30.
  * [28] N. J. Higham, S. Pranesh, and M. Zounon, Squeezing a matrix into half precision, with an application to solving linear systems, SIAM J. Sci. Comput., 41 (2019), pp. A2536–A2551.
  * [29] T. Huckle, Factorized sparse approximate inverses for preconditioning, J. Supercomput., 25 (2003), pp. 109–117.
  * [30] T. Huckle and A. Kallischko, Frobenius norm minimization and probing for preconditioning, Int. J. Comput. Math., 84 (2007), pp. 1225–1248.
  * [31] M. Jankowski and H. Woźniakowski, Iterative refinement implies numerical stability, BIT Numer. Math., 17 (1977), pp. 303–311.
  * [32] L. Y. Kolotilina and A. Y. Yeremin, Factorized sparse approximate inverse preconditionings I. Theory, SIAM J. Matrix Anal. Appl., 14 (1993), pp. 45–58.
  * [33] J. Langou, J. Langou, P. Luszczek, J. Kurzak, A. Buttari, and J. Dongarra, Exploiting the performance of 32 bit floating point arithmetic in obtaining 64 bit accuracy (revisiting iterative refinement for linear systems), in Proc. 2006 ACM/IEEE Conf. Supercomput., 2006.
  * [34] N. Lindquist, P. Luszczek, and J. Dongarra, Accelerating restarted GMRES with mixed precision arithmetic, IEEE Trans. Par. Dist. Sys., 33 (2022), pp. 1027–1037.
  * [35] M. Lukash, K. Rupp, and S. Selberherr, Sparse approximate inverse preconditioners for iterative solvers on GPUs, in Proc. 2012 Symp. High Perf. Comput., Society for Computer Simulation San Diego, CA, USA, 2012, p. 13.
  * [36] C. B. Moler, Iterative refinement in floating point, J. ACM, 14 (1967), pp. 316–321.
  * [37] Multiprecision Computing Toolbox.  Advanpix, Tokyo.  http://www.advanpix.com.
  * [38] E. Oktay and E. Carson, Multistage mixed precision iterative refinement, Numerical Linear Algebra with Applications, (2022), p. e2434.
  * [39] M. Sedlacek, Sparse approximate inverses for preconditioning, smoothing, and regularization, PhD thesis, Technische Universität München, 2012.
  * [40] R. D. Skeel, Iterative refinement implies numerical stability for Gaussian elimination, Math. Comput., 35 (1980), pp. 817–832.
  * [41] J. H. Wilkinson, Progress report on the Automatic Computing Engine, Report MA/17/1024, Mathematics Division, Department of Scientific and Industrial Research, National Physical Laboratory, Teddington, UK, Apr. 1948\.
  * [42] J. H. Wilkinson, Rounding Errors in Algebraic Processes, Notes on Applied Science No. 32, Her Majesty’s Stationery Office, London, 1963.  Also published by Prentice-Hall, Englewood Cliffs, NJ, USA. Reprinted by Dover, New York, 1994.
  * [43] A. Y. Yeremin and A. Nikishin, Factorized-sparse-approximate-inverse preconditionings of linear systems with unsymmetric matrices, J. Math. Sci., 121 (2004), pp. 2448–2457.