---
title: 'GMRES-IR: Mixed-Precision Iterative Refinement'
url: https://www.emergentmind.com/topics/gmres-based-iterative-refinement-gmres-ir
type: topic
---

# GMRES-IR: Mixed-Precision Iterative Refinement

GMRES-based iterative refinement (GMRES-IR) is a numerical solver paradigm for large sparse linear systems and related problems, leveraging the Generalized Minimal Residual (GMRES) algorithm as an inner solver within an iterative refinement loop. It is specifically designed to exploit mixed-precision or low-precision arithmetic (e.g., fp16/fp32), prevalent in modern accelerator hardware, to achieve double-precision (fp64) or high-accuracy solutions at a reduced computational and energy cost. The method’s central rationale is that only a subset of operations (notably residual computation and final solution updates) must be executed in high precision to ensure backward and forward stability, while the bulk of computationally intensive linear algebra kernels can be performed in lower precision, substantially accelerating the overall solve without compromising the final solution accuracy when proper algorithmic safeguards are provided [2109.01232].

## 1. Algorithmic Framework and Mathematical Foundation

GMRES-IR is structured as an outer iterative-refinement loop, with each iteration comprising a residual computation, an (approximate) solution of a linear system using GMRES (possibly preconditioned), and a high-precision update:

Let $A\in\mathbb{R}^{n\times n}$, $b\in\mathbb{R}^n$.

1. High-precision computation of the initial residual:
   $$
   r_0 = b - A x_0 \quad \text{(fp64)}
   $$
2. For $k=0,1,2,...$ until convergence:
   - **Inner solve in lower precision (e.g., fp32):**
     $$
     \text{Find } u_k \approx A^{-1} r_k \; \text{(via GMRES$(m)$ in fp32)}
     $$
   - **Correction update in high precision:**
     $$
     x_{k+1} = x_k + \text{cast}_{32\to 64}(u_k) \quad \text{(fp64)}
     $$
   - **New residual in high precision:**
     $$
     r_{k+1} = b - A x_{k+1} \quad \text{(fp64)}
     $$
   - **Stopping criterion:**
     $$
     \|r_{k+1}\|_2 / \|b\|_2 \leq \text{tol}
     $$
Here, all GMRES subspace operations, SpMV (sparse matrix-vector multiplication), inner products, and (explicit) preconditioning are performed in fp32, while the residual evaluation and the update are computed in fp64 [2109.01232], [2105.07544].

## 2. Mixed-Precision Strategy and Theoretical Rationale

Performance gains on modern architectures are realized as fp32 operations typically yield 2–3× throughput compared to fp64. However, scientific computing often mandates fp64 accuracy. The iterative refinement theory (Higham, Moler; Carson & Higham) establishes that, provided the system condition satisfies $\kappa(A)\cdot\epsilon_\text{fp32} < 1$ (for single-precision inner solves), the outer iteration converges rapidly to double precision accuracy. Contemporary refinements show that only the outer residual computation and update require fp64; the inner linear system correction can be safely performed in fp32 without loss of final accuracy [2109.01232].

For large, moderately to well-conditioned $A$, the accumulation of roundoff in the fp32 inner solve is corrected at each outer step by a recomputed fp64 residual; in practice, only a few outer iterations are needed to achieve double precision accuracy.

## 3. Implementation, Preconditioning, and Parallelism

GMRES-IR is architected for performance portability and high concurrency, with all computationally heavy kernels (SpMV, orthogonalization, preconditioning) executed in low precision. Preconditioners must be highly parallelizable: classical incomplete LU (ILU) may introduce scalability bottlenecks due to triangular solves; thus, block–Jacobi and polynomial preconditioning are preferred [2109.01232], [2105.07544].

- **Block–Jacobi**: Decompose $A$ into block diagonal $s\times s$ blocks, invert each in fp32, and apply in parallel to enable full GPU occupancy.
- **Polynomial preconditioning**: Employs a low-degree polynomial $p(A)$ (often a truncated GMRES polynomial), all in fp32, to approximate $A^{-1}$, yielding a sequence of fully parallel SpMVs.

All preconditioner applications are realized in fp32: casting fp64 vectors to fp32, applying the preconditioner, and casting back as necessary. GPU implementations (e.g., via Kokkos Kernels) further streamline memory transfers and enable fused operations.

## 4. Parameter Selection and Convergence

Key parameters influencing GMRES-IR performance include:

- **Restart length ($m$)** for GMRES: Typical values are $25 \leq m \leq 100$. Larger $m$ reduces outer iterations but increases orthogonalization overhead; smaller $m$ produces more refinement iterations [2109.01232].
- **Preconditioner strength**: Block size $s$ for block–Jacobi or polynomial degree $d$ are chosen to minimize the total number of inner iterations, targeting values that keep iteration counts moderate ($500$–$2000$) without bloating the cost of applying the preconditioner.
- **Refinement tolerance**: Set to the target fp64 accuracy (e.g., $10^{-10}$).

Convergence is typically observed in 1–3 outer iterations, with the total fp32 iteration count closely tracking that of a full-fp64 GMRES(m) with only a minor overhead.

## 5. Performance Evaluation and Practical Outcomes

Extensive benchmarking on NVIDIA V100 GPUs with a Power9 host has demonstrated that GMRES-IR can deliver systemic performance gains [2109.01232], [2105.07544]:

- **Kernel-level speedups**:
  - SpMV: $2.4$–$2.6\times$ over fp64 for matrices with $\leq 15$ nonzeros per row.
  - Orthogonalization (GEMV): $1.3$–$1.6\times$ speedup.
- **Full-solver speedups**:
  - Unpreconditioned GMRES(50): $1.32$–$1.44\times$ faster.
  - Polynomial preconditioning: up to $1.58\times$ speedup.
  - Block–Jacobi: $1.5\times$ faster for CFD-type matrices (block size 42).
- **Accuracy**: In all experiments, the double precision residual norm $\|b-Ax\|_2$ obtained via GMRES-IR matched full-fp64 GMRES down to round-off (e.g., $10^{-10}$–$10^{-12}$), even for systems with $\kappa(A)\sim 10^6$–$10^8$. Provided $\kappa(A)\cdot\epsilon_\text{fp32} < 10^{-1}$, convergence is robust.
- **Scalability**: When the number of nonzeros per row increases ($>15$), SpMV and thus overall speedup declines modestly but remains significant.

## 6. Algorithmic Variants and Extensions

### 6.1. SPAI-GMRES-IR and Adaptive-Precision Variants

Variants using adaptive or sparse approximate inverse (SPAI) preconditioning further tune precision at the preconditioner level. BSPAI-GMRES-IR stores SPAI components in variable precision “buckets,” enabling a potentially significant reduction in storage and computational cost, with only a modest increase in the number of GMRES iterations. Provided the lowest bucket precision matches the working precision, convergence and accuracy are essentially unchanged [2307.03914], [2202.10204].

### 6.2. Integer Arithmetic GMRES-IR

An integer-arithmetic-based variant substitutes all inner GMRES computations with fixed-point arithmetic, relying on precise scaling and logical operand shifts to prevent overflow. The outer (floating-point) refinement ensures that accuracy remains competitive with pure floating-point solvers, as any drift in the low-precision GMRES is corrected at each iteration [2009.07495].

### 6.3. Krylov Subspace Recycling

For problems where multiple refinement steps are needed, recycling previous Krylov spaces (e.g., via GCRO-DR) across refinement iterations can reduce the computational effort per outer step, particularly for ill-conditioned matrices, by reusing spectral information [2201.09827].

## 7. Applications and Demonstrated Scenarios

GMRES-IR and its mixed-precision extensions are immediately applicable in large-scale computational science and engineering workflows that demand both high accuracy and hardware efficiency. Key contexts include:

- Solution of large, sparse, nonsymmetric systems arising from PDE discretizations, model reduction, or CFD.
- Accelerated iterative refinement for dense and sparse overdetermined least-squares, including augmented system formulations and weighted LS problems, with mixed or low-precision QR preconditioners [2410.06319], [2401.03755], [2406.16499].
- Matrices encountered in practical settings, such as convection-diffusion or CFD, even for challenging condition numbers up to $10^{8}$, provided proper parameter selection is observed.

## Summary Table: Core GMRES-IR Workflow

| Step                        | Precision   | Main Operations                      |
|-----------------------------|-------------|--------------------------------------|
| Residual computation        | High (fp64) | $r_k = b - A x_k$                    |
| Inner GMRES solve           | Low (fp32)  | Iterative correction, SpMV, Precon   |
| Preconditioning             | Low (fp32)  | Block–Jacobi, polynomial, or SPAI    |
| Solution update             | High (fp64) | $x_{k+1} = x_k + u_k$                |
| Convergence check           | High (fp64) | Residual norm test                   |

## References

- "A Study of Mixed Precision Strategies for GMRES on GPUs" [2109.01232]
- "Experimental Evaluation of Multiprecision Strategies for GMRES on GPUs" [2105.07544]
- "Mixed Precision Iterative Refinement with Sparse Approximate Inverse Preconditioning" [2202.10204]
- "Mixed Precision Iterative Refinement with Adaptive Precision Sparse Approximate Inverse Preconditioning" [2307.03914]
- "An Integer Arithmetic-Based Sparse Linear Solver Using a GMRES Method and Iterative Refinement" [2009.07495]
- "Mixed Precision GMRES-based Iterative Refinement with Recycling" [2201.09827]
- "Mixed precision sketching for least-squares problems and its application in GMRES-based iterative refinement" [2410.06319]
- "Mixed Precision FGMRES-Based Iterative Refinement for Weighted Least Squares" [2401.03755]
- "Mixed precision iterative refinement for least squares with linear equality constraints and generalized least squares problems" [2406.16499]

Source: https://www.emergentmind.com/topics/gmres-based-iterative-refinement-gmres-ir