GMRES Polynomials

• GMRES polynomials are residual-minimizing polynomials normalized at the origin, capturing the action of A⁻¹ and linking convergence bounds with harmonic Ritz values.
• They are derived from the Arnoldi process and expressed in factored forms to enhance numerical stability and approximate inversion for nonsingular systems.
• They underpin practical methods like polynomial preconditioning, communication-avoiding Krylov solvers, and deflation strategies to improve efficiency in large-scale linear algebra problems.

GMRES polynomials are the residual-minimizing polynomials implicitly generated by the generalized minimal residual method for a nonsingular linear system $Ax=b$. At iteration $k$, GMRES selects a polynomial of degree at most $k$, normalized by its value at the origin, so that the residual $r_k$ is of the form $r_k=p_k(A)r_0$ and has minimal norm. In this form, the polynomial viewpoint unifies residual minimization, harmonic Ritz values, worst-case and ideal approximation problems, polynomial preconditioning, inverse approximation, convergence bounds on spectral sets, and communication-avoiding Krylov implementations (Faber et al., 2013, Embree et al., 2018, Loe et al., 2019, Embree, 2022).

1. Residual-polynomial formulation

Let $A\in\mathbb C^{n\times n}$ be nonsingular, let $x_0$ be an initial guess, and let $r_0=b-Ax_0$. GMRES seeks $x_k\in x_0+\mathcal K_k(A,r_0)$ by minimizing the residual norm. Equivalently, there exists a polynomial $q_k$ of degree at most $k$0 with $k$1 such that

$k$2

A standard reparametrization writes

$k$3

so that the same GMRES iterate can be expressed through a degree-$k$4 polynomial $k$5 applied to $k$6, or equivalently through a degree-$k$7 polynomial acting as an approximate inverse on $k$8 (Embree, 2022, Loe et al., 2019).

This formulation extends directly to the best possible residual bound that is independent of the initial vector. For real nonsingular $k$9, Faber, Liesen, and Tichý define

$k$0

where

$k$1

For fixed $k$2, the minimizing polynomial is the usual GMRES polynomial; $k$3 is the sharp, $k$4-independent upper bound on all GMRES residuals at step $k$5 (Faber et al., 2013).

The same polynomial viewpoint also underlies “ideal GMRES,” which minimizes the operator norm rather than the residual for one starting vector: $k$6 The distinction between $k$7 and $k$8 is fundamental in the theory of GMRES polynomials, because the two quantities can coincide, differ mildly, or differ sharply depending on the matrix and iteration index (Faber et al., 2013).

2. Arnoldi, harmonic Ritz values, and factored representations

The Arnoldi process provides the finite-dimensional representation from which the GMRES polynomial is extracted. After $k$9 steps,

$r_k$0

where $r_k$1 spans the Krylov subspace and $r_k$2 is upper Hessenberg. If $r_k$3, GMRES solves the least-squares problem

$r_k$4

and the residual satisfies

$r_k$5

This is the standard bridge between the infinite-dimensional polynomial problem and the small Arnoldi least-squares problem (Embree et al., 2018, Loe et al., 2019).

A central structural fact is that the roots of the GMRES residual polynomial coincide with the harmonic Ritz values. If $r_k$6 are these roots, then

$r_k$7

The factored form is computationally important: it avoids explicit monomial expansion, can be applied by successive factors $r_k$8, and admits modified Leja reordering for numerical stability. For real matrices, complex roots occur in conjugate pairs and may be combined in real arithmetic through

$r_k$9

These representations are used both in polynomial-preconditioned GMRES and in polynomial-preconditioned Arnoldi (Embree et al., 2018, Loe et al., 2019).

The harmonic-Ritz factorization also clarifies the approximate-inverse interpretation. If

$r_k=p_k(A)r_0$0

then

$r_k=p_k(A)r_0$1

Thus the same polynomial that minimizes the residual induces a polynomial approximation to $r_k=p_k(A)r_0$2 on the Krylov-generated subspace. This identity is the basis for later uses of GMRES polynomials as reusable approximate inverses and as filters for eigenvalue computations (Embree et al., 25 Feb 2025, Henson et al., 16 Oct 2025).

3. Worst-case, ideal, and singular-vector characterizations

The worst-case GMRES problem asks for the largest attainable $r_k=p_k(A)r_0$3-step minimal residual over all unit right-hand sides. Faber, Liesen, and Tichý show that the worst-case behavior of GMRES for $r_k=p_k(A)r_0$4 and $r_k=p_k(A)r_0$5 is the same and that initial vectors attaining the worst case satisfy a “cross equality” (Faber et al., 2013). If $r_k=p_k(A)r_0$6 is a unit vector with $r_k=p_k(A)r_0$7, $r_k=p_k(A)r_0$8, and if

$r_k=p_k(A)r_0$9

then

$A\in\mathbb C^{n\times n}$0

This relation identifies worst-case initial vectors as eigenvectors of a composite residual operator.

The stationarity analysis goes further. One can choose $A\in\mathbb C^{n\times n}$1, so that

$A\in\mathbb C^{n\times n}$2

Hence $A\in\mathbb C^{n\times n}$3 is a right singular vector of the residual matrix $A\in\mathbb C^{n\times n}$4 associated with its largest singular value $A\in\mathbb C^{n\times n}$5. If

$A\in\mathbb C^{n\times n}$6

then

$A\in\mathbb C^{n\times n}$7

This recasts worst-case GMRES in the language of singular-value approximation problems (Faber et al., 2013).

The worst-case polynomial need not be unique. A necessary and sufficient condition for uniqueness is the simplicity of the top singular value of $A\in\mathbb C^{n\times n}$8 at the maximizer $A\in\mathbb C^{n\times n}$9. If that singular value has multiplicity greater than one, distinct right singular vectors can generate different polynomials that achieve the same max–min value. In Toh’s $x_0$0 example,

$x_0$1

one has $x_0$2, and both $x_0$3 and $x_0$4 solve the worst-case problem, with $x_0$5 numerically (Faber et al., 2013).

For comparison with ideal GMRES,

$x_0$6

Equality holds in particular if $x_0$7 is normal for any $x_0$8, or if $x_0$9 for any nonsingular $r_0=b-Ax_0$0. A more general characterization is that equality holds exactly when the corresponding max–min and min–max problems coincide, equivalently when

$r_0=b-Ax_0$1

has a saddle-point in $r_0=b-Ax_0$2. The same paper also establishes that, for real matrices,

$r_0=b-Ax_0$3

with both inequalities potentially strict, whereas $r_0=b-Ax_0$4 is invariant under all four choices of real or complex polynomials and vectors (Faber et al., 2013).

4. Polynomial preconditioning and approximate inverses

Polynomial preconditioning uses the GMRES minimal-residual polynomial itself as a preconditioner. In one formulation, GMRES constructs

$r_0=b-Ax_0$5

and applies $r_0=b-Ax_0$6 as a left preconditioner: $r_0=b-Ax_0$7 Because $r_0=b-Ax_0$8, the preconditioner remains nonsingular for typical spectra, and eigenvalues $r_0=b-Ax_0$9 are mapped to $x_k\in x_0+\mathcal K_k(A,r_0)$0. The intended effect is that small or ill-conditioned portions of the spectrum are flattened or moved away from the origin, while the bulk of the spectrum is clustered near $x_k\in x_0+\mathcal K_k(A,r_0)$1 (Loe et al., 2019).

The same mechanism can be expressed through the approximate-inverse polynomial. If

$x_k\in x_0+\mathcal K_k(A,r_0)$2

then a small GMRES residual implies that $x_k\in x_0+\mathcal K_k(A,r_0)$3 is a good inverse on the subspace spanned by $x_k\in x_0+\mathcal K_k(A,r_0)$4. Embree and coauthors state the bounds

$x_k\in x_0+\mathcal K_k(A,r_0)$5

and, when $x_k\in x_0+\mathcal K_k(A,r_0)$6 is diagonalizable with $x_k\in x_0+\mathcal K_k(A,r_0)$7,

$x_k\in x_0+\mathcal K_k(A,r_0)$8

They also note that once $x_k\in x_0+\mathcal K_k(A,r_0)$9 has been built, each new right-hand side can be solved by one polynomial matrix-vector evaluation, and they identify variance reduction in multilevel Monte Carlo methods as another application (Embree et al., 25 Feb 2025).

High-degree and composite constructions are common. In polynomial-preconditioned Arnoldi, one defines

$q_k$0

but in practice applies the factored residual polynomial through blocks $q_k$1 rather than inverting $q_k$2. A two-level “double polynomial preconditioning” strategy composes two GMRES-derived filters: if $q_k$3 is built first and $q_k$4, then a second polynomial $q_k$5 on $q_k$6 gives the composite polynomial

$q_k$7

of degree $q_k$8, without forming the monolithic high-degree polynomial explicitly (Embree et al., 2018).

For indefinite matrices, GMRES-derived polynomial preconditioners require additional control. One approach is to enforce a balanced polynomial through

$q_k$9

so that $k$00 for $k$01 near zero. The same work also describes range-restricted GMRES polynomials, composite inner–outer constructions, and convergence estimates for real indefinite spectra, with tests on interior eigenvalue computations such as diagonal problems, Af23560, and a $k$02-Hermitian QCD matrix (Henson et al., 16 Oct 2025).

5. Stability, high degree, and communication-avoiding polynomial bases

The main numerical difficulty with GMRES polynomials is stability at high degree. Monomial expansions are unstable, and even factored products can become ill-conditioned if one root lies very close to an eigenvalue. A practical monitor used in polynomial-preconditioned Arnoldi is

$k$03

Large values warn that the polynomial may exhibit severe amplification before cancellation. The remedy proposed in several works is root-adding: duplicate a troublesome root so that the derivative of the polynomial vanishes there, flattening the factorization locally. Practical thresholds include adding one extra copy when $k$04, and adding further copies when values exceed $k$05, with analogous rules expressed in terms of $k$06 in related formulations (Embree et al., 2018, Loe et al., 2019, Embree et al., 25 Feb 2025).

Additional stability controls are problem dependent. One may damp the starting vector by building the polynomial from $k$07 or $k$08 instead of $k$09, average two starting vectors through a block-diagonal system $k$10, or, in the indefinite case, add copies of roots only on the long side of the spectrum, deflate approximate eigenvectors associated with large $k$11, or apply a few GMRES iterations to cleanse the residual (Embree et al., 2018, Henson et al., 16 Oct 2025).

Communication-avoiding GMRES changes not only how the polynomial is used but also which polynomial basis is employed. In $k$12-step GMRES, the matrix-powers kernel generates

$k$13

where $k$14 can be the monomial basis

$k$15

the Newton basis

$k$16

or the scaled Newton basis

$k$17

The scaled-Newton basis is designed so that the block norms remain $k$18, slowing the growth of $k$19. Reported practice is that scaled-Newton allows $k$20 in the hundreds or even thousands, whereas monomials or unscaled Newton are limited to $k$21 (Xu et al., 2023).

Automatic stability control in adaptive $k$22-step GMRES combines BCGS2 with partial Cholesky–QR and an incremental condition estimator. The factorization is truncated as soon as the estimate exceeds a user-set bound $k$23, guaranteeing $k$24 loss of orthogonality under the stated conditioning requirement. An initial step-size estimator based on

$k$25

selects

$k$26

typically with $k$27 (Xu et al., 2023).

6. Polynomial convergence bounds and deflation

The classical convergence theory of GMRES is itself a polynomial approximation theory. One always has

$k$28

and the remaining question is how to bound $k$29 from information about spectral sets. Three standard paradigms place the min–max problem on the spectrum, the field of values, or the pseudospectra (Embree, 2022).

For diagonalizable $k$30, the spectrum-based estimate yields

$k$31

up to the usual condition-number factor, and for normal matrices it is exact. The field-of-values estimate uses Crouzeix’s theorem,

$k$32

leading to a field-of-values bound. The pseudospectral estimate uses the Dunford integral and an infimum over $k$33, producing the form

$k$34

Embree’s analysis shows, through six examples, that none of these three classes of bounds is uniformly descriptive: spectrum, field-of-values, and pseudospectral bounds each succeed in some cases and fail in others (Embree, 2022).

A related framework replaces the actual spectral set by a rectangle. For preconditioned, weighted, or deflated GMRES, one introduces

$k$35

where $k$36 contains the effective field of values. Then

$k$37

In the SPD case this reduces to the classical Chebyshev rate

$k$38

For rectangles, explicit bounds include Elman-type, disk-segment, ellipse-based, conformal-mapping, and Faber-polynomial estimates (Spillane et al., 8 Apr 2025).

Deflation modifies the polynomial problem by shrinking the relevant set. With the standard deflation projector

$k$39

GMRES is applied to $k$40, and the residual again has the form

$k$41

If the deflation space captures the modes responsible for the large skew-Hermitian part, the containing rectangle narrows from $k$42 to

$k$43

and the corresponding $k$44 is markedly smaller (Spillane et al., 8 Apr 2025).

7. Generalized settings and application domains

GMRES polynomials also appear in settings that are not the standard complex-linear Krylov method. In the R-linear GMRES for systems of the form

$k$45

the convergence bound is controlled by a polynomial approximation problem over a finite subset of the spectrum, now involving polynomials in $k$46. In the CSYM case, one obtains a discrete inner product

$k$47

and a family of orthogonal polynomials satisfying the three-term recurrence

$k$48

with the same coefficients forming a complex-symmetric Jacobi matrix. This creates a mathematical framework analogous to the Hermitian Lanczos setting, but for the complex-symmetric case (Huhtanen et al., 2011).

In multigrid for Boltzmann transport, low-order GMRES polynomials are used as explicit approximations to $k$49. For fixed order $k$50, one constructs

$k$51

and then uses $k$52 to define

$k$53

as well as F-point relaxation updates. With fixed sparsity, the truncated polynomial approximation preserves the sparsity pattern of $k$54, limiting memory and communication (Dargaville et al., 2023).

Large-scale parallel solvers exploit the same polynomial ideas for communication reduction. Polynomial-preconditioned GMRES moves work from global reductions into sparse matrix-vector products. Reported examples include the matrix “e20r0100,” where unpreconditioned GMRES(100) needs $k$55 million dot products to converge or stalls, while PP-GMRES with $k$56 reduces this to $k$57 thousand dot products, and convection–diffusion tests where PP-GMRES with auto-selected degree attains substantial speedups while shrinking the fraction of runtime spent in dot-product reductions (Loe et al., 2019). This suggests that GMRES polynomials are not only an abstract convergence device but also a concrete mechanism for reducing synchronization costs in large Krylov solvers.

Across these variants, the common object is unchanged: a polynomial normalized at the origin and adapted to the operator through Arnoldi data, harmonic Ritz values, or residual minimization. What changes is the role assigned to that polynomial—residual certificate, worst-case extremizer, approximate inverse, spectral filter, communication-avoiding basis, deflation-aware bound, or multigrid surrogate for a block inverse.