--- title: GMRES Polynomials url: https://www.emergentmind.com/topics/gmres-polynomials type: topic --- # GMRES Polynomials 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 [1302.5535][1806.08020][1907.00072][2209.01231]. ## 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\) with \(q_k(0)=1\) such that \[ r_k = q_k(A)\,r_0, \qquad \|r_k\|_2 = \min_{p\in\mathbb P_k,\;p(0)=1}\|p(A)r_0\|_2. \] A standard reparametrization writes \[ q_k(z)=1-z\,\pi_{k-1}(z), \] so that the same GMRES iterate can be expressed through a degree-\((k-1)\) polynomial \(\pi_{k-1}\) applied to \(r_0\), or equivalently through a degree-\((k-1)\) polynomial acting as an approximate inverse on \(b\) [2209.01231][1911.07065]. This formulation extends directly to the best possible residual bound that is independent of the initial vector. For real nonsingular \(A\in\mathbb R^{n\times n}\), Faber, Liesen, and Tichý define \[ \Psi_k(A)=\max_{\|b\|=1}\min_{p\in\pi_k}\|p(A)b\|, \] where \[ \pi_k=\{\,p(z)=1-c_1z-\cdots-c_k z^k\mid c\in\mathbb R^k\,\}. \] For fixed \(b\), the minimizing polynomial is the usual GMRES polynomial; \(\Psi_k(A)\) is the sharp, \(b\)-independent upper bound on all GMRES residuals at step \(k\) [1302.5535]. The same polynomial viewpoint also underlies “ideal GMRES,” which minimizes the operator norm rather than the residual for one starting vector: \[ \varphi_k(A)=\min_{p\in\pi_k}\|p(A)\| =\min_{p\in\pi_k}\max_{\|v\|=1}\|p(A)v\|. \] The distinction between \(\Psi_k(A)\) and \(\varphi_k(A)\) 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 [1302.5535]. ## 2. Arnoldi, harmonic Ritz values, and factored representations The Arnoldi process provides the finite-dimensional representation from which the GMRES polynomial is extracted. After \(m\) steps, \[ A\,V_m = V_{m+1}\,H_{m+1,m}, \qquad V_{m+1}^*V_{m+1}=I, \] where \(V_m\) spans the Krylov subspace and \(H_{m+1,m}\) is upper Hessenberg. If \(\beta=\|r_0\|\), GMRES solves the least-squares problem \[ y_m=\arg\min_{y\in\mathbb C^m}\|\beta e_1-H_{m+1,m}y\|_2, \] and the residual satisfies \[ r_m(A)\,r_0 =\beta\,V_{m+1}\Bigl(e_1-\tfrac{1}{\beta}H_{m+1,m}y_m\Bigr). \] This is the standard bridge between the infinite-dimensional polynomial problem and the small Arnoldi least-squares problem [1806.08020][1907.00072]. A central structural fact is that the roots of the GMRES residual polynomial coincide with the harmonic Ritz values. If \(\{\theta_i\}_{i=1}^m\) are these roots, then \[ r_m(z)=\prod_{i=1}^m\Bigl(1-\frac{z}{\theta_i}\Bigr), \qquad p_{m-1}(z)=\frac{1-r_m(z)}{z}. \] The factored form is computationally important: it avoids explicit monomial expansion, can be applied by successive factors \((I-A/\theta_i)\), 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 \[ \Bigl(1-\frac{z}{a+bi}\Bigr)\Bigl(1-\frac{z}{a-bi}\Bigr) =1+\frac{z^2-2az}{a^2+b^2}. \] These representations are used both in polynomial-preconditioned GMRES and in polynomial-preconditioned Arnoldi [1806.08020][1911.07065]. The harmonic-Ritz factorization also clarifies the approximate-inverse interpretation. If \[ r_m = \pi_m(A)b, \qquad \pi_m(z)=1-z\,p_{m-1}(z), \] then \[ x_m=p_{m-1}(A)b. \] Thus the same polynomial that minimizes the residual induces a polynomial approximation to \(A^{-1}\) 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 [2502.18317][2510.14816]. ## 3. Worst-case, ideal, and singular-vector characterizations The worst-case GMRES problem asks for the largest attainable \(k\)-step minimal residual over all unit right-hand sides. Faber, Liesen, and Tichý show that the worst-case behavior of GMRES for \(A\) and \(A^T\) is the same and that initial vectors attaining the worst case satisfy a “cross equality” [1302.5535]. If \(b_0\) is a unit vector with \(\|r_k\|=\Psi_k(A)\), \(r_k=p_k(A)b_0\), and if \[ c_0=\frac{r_k}{\|r_k\|}, \qquad s_k=q_k(A^T)c_0, \qquad \|s_k\|=\Psi_k(A), \] then \[ q_k(A^T)\,p_k(A)\,b_0=\Psi_k(A)^2\,b_0. \] This relation identifies worst-case initial vectors as eigenvectors of a composite residual operator. The stationarity analysis goes further. One can choose \(q_k=p_k\), so that \[ p_k(A)^T\,p_k(A)\,b_0=\Psi_k(A)^2\,b_0. \] Hence \(b_0\) is a right singular vector of the residual matrix \(R_k=p_k(A)\) associated with its largest singular value \(\sigma_1=\Psi_k(A)\). If \[ R_k=U\Sigma V^T, \] then \[ Ve_1=b_0,\qquad \Sigma_{11}=\Psi_k(A),\qquad Ue_1=\frac{r_k}{\|r_k\|}. \] This recasts worst-case GMRES in the language of singular-value approximation problems [1302.5535]. The worst-case polynomial need not be unique. A necessary and sufficient condition for uniqueness is the simplicity of the top singular value of \(p_k(A)\) at the maximizer \(b_0\). 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 \(4\times4\) example, \[ A(\omega,\varepsilon)= \begin{pmatrix} 1 & \varepsilon & & \\ & -1 & \omega/\varepsilon & \\ & & 1 & \varepsilon \\ & & & -1 \end{pmatrix}, \] one has \(\Psi_3(A)<\varphi_3(A)\), and both \(p_3(z)\) and \(p_3(-z)\) solve the worst-case problem, with \(p_3(z)\neq p_3(-z)\) numerically [1302.5535]. For comparison with ideal GMRES, \[ \Psi_k(A)\le \varphi_k(A). \] Equality holds in particular if \(A\) is normal for any \(k\), or if \(k=1\) for any nonsingular \(A\). A more general characterization is that equality holds exactly when the corresponding max–min and min–max problems coincide, equivalently when \[ f(c,v)=\|(1-c_1A-\cdots-c_kA^k)v\|^2 \] has a saddle-point in \((c,v)\in\mathbb R^k\times S^{n-1}\). The same paper also establishes that, for real matrices, \[ \Psi_{k,\mathbb R,\mathbb R}(A)=\Psi_{k,\mathbb C,\mathbb R}(A) \le \Psi_{k,\mathbb C,\mathbb C}(A) \le \Psi_{k,\mathbb R,\mathbb C}(A), \] with both inequalities potentially strict, whereas \(\varphi_k\) is invariant under all four choices of real or complex polynomials and vectors [1302.5535]. ## 4. Polynomial preconditioning and approximate inverses Polynomial preconditioning uses the GMRES minimal-residual polynomial itself as a preconditioner. In one formulation, GMRES constructs \[ \phi_m(x)=\arg\min_{p\in\mathcal P_m,\;p(0)=1}\|p(A)b\|_2, \] and applies \(\phi_m(A)\) as a left preconditioner: \[ \phi_m(A)\,A\,x=\phi_m(A)\,b. \] Because \(\phi_m(0)=1\), the preconditioner remains nonsingular for typical spectra, and eigenvalues \(\lambda\) are mapped to \(\lambda\,\phi_m(\lambda)\). 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 \(1\) [1907.00072]. The same mechanism can be expressed through the approximate-inverse polynomial. If \[ r_m=b-A\,p_{m-1}(A)b=(I-A\,p_{m-1}(A))b, \] then a small GMRES residual implies that \(p_{m-1}(A)\) is a good inverse on the subspace spanned by \(b\). Embree and coauthors state the bounds \[ \frac{\|r\|}{\kappa(A)} \le \frac{\|A^{-1}-p(A)\|}{\|A^{-1}\|}, \] and, when \(A=Z\Lambda Z^{-1}\) is diagonalizable with \(b=\sum \widehat\beta_i z_i\), \[ \frac{\|A^{-1}-p(A)\|}{\|A^{-1}\|} \le \kappa(Z)\,\frac{\|r\|}{\min_i|\beta_i|}. \] They also note that once \(p(A)\approx A^{-1}\) 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 [2502.18317]. High-degree and composite constructions are common. In polynomial-preconditioned Arnoldi, one defines \[ p_{m-1}(A)=A^{-1}(I-r_m(A)), \] but in practice applies the factored residual polynomial through blocks \((I-A/\theta_i)\) rather than inverting \(A\). A two-level “double polynomial preconditioning” strategy composes two GMRES-derived filters: if \(\pi_1\) is built first and \(\tau(A)=I-\pi_1(A)\), then a second polynomial \(\pi_2\) on \(\tau(A)\) gives the composite polynomial \[ \pi_{\rm comp}(z)=\pi_2\bigl(1-\pi_1(z)\bigr), \] of degree \(d_1d_2\), without forming the monolithic high-degree polynomial explicitly [1806.08020]. For indefinite matrices, GMRES-derived polynomial preconditioners require additional control. One approach is to enforce a balanced polynomial through \[ \varphi_k'(0)=\sum_{i=1}^k\frac{1}{\theta_i}=0, \] so that \(\varphi_k(\lambda)\approx1\) for \(\lambda\) 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 \(\gamma_5\)-Hermitian QCD matrix [2510.14816]. ## 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 \[ P_{\neg j}=\prod_{i\ne j}\bigl|1-\theta_j/\theta_i\bigr|. \] 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 \(P_{\neg j}>10^4\), and adding further copies when values exceed \(10^{18}\), with analogous rules expressed in terms of \(\mathit{pof}(k)\) in related formulations [1806.08020][1911.07065][2502.18317]. Additional stability controls are problem dependent. One may damp the starting vector by building the polynomial from \(A b\) or \(A^k b\) instead of \(b\), average two starting vectors through a block-diagonal system \(\operatorname{diag}(A,A)\), or, in the indefinite case, add copies of roots only on the long side of the spectrum, deflate approximate eigenvectors associated with large \(pof(\theta_j)\), or apply a few GMRES iterations to cleanse the residual [1806.08020][2510.14816]. Communication-avoiding GMRES changes not only how the polynomial is used but also which polynomial basis is employed. In \(s\)-step GMRES, the matrix-powers kernel generates \[ V(Q_i,s)=\bigl[p_0(A)Q_i,\;p_1(A)Q_i,\;\dots,\;p_s(A)Q_i\bigr], \] where \(p_j\) can be the monomial basis \[ p_j(A)=A^j, \] the Newton basis \[ p_j(A)=\prod_{k=1}^j(A-\theta_k I), \] or the scaled Newton basis \[ p_j(A)=\prod_{k=1}^j\frac{A-\theta_k I}{\gamma_k}, \qquad \gamma_k=\bigl|\bar\theta-\theta_k\bigr|, \qquad \bar\theta=\frac1s\sum_{k=1}^s\theta_k. \] The scaled-Newton basis is designed so that the block norms remain \(O(1)\), slowing the growth of \(\kappa(V(Q_i,s))\). Reported practice is that scaled-Newton allows \(s\) in the hundreds or even thousands, whereas monomials or unscaled Newton are limited to \(s\approx5\text{–}20\) [2303.08953]. Automatic stability control in adaptive \(s\)-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 \(\Omega\sim10^7\), guaranteeing \(O(\epsilon)\) loss of orthogonality under the stated conditioning requirement. An initial step-size estimator based on \[ E_{i,j} = \prod_{k=1}^{j-1} \frac{|\theta_i-\theta_k|}{|\bar\theta-\theta_k|} \times \epsilon^{\,\delta_{i0\), producing the form \[ \|r_k\|_2 \le \|r_0\|_2\, \inf_{\epsilon>0} \Bigl[ \min_{p\in\mathbb P_k,\;p(0)=1}\max_{z\in\Lambda_\epsilon(A)}|p(z)| +\epsilon \Bigr]. \] 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 [2209.01231]. A related framework replaces the actual spectral set by a rectangle. For preconditioned, weighted, or deflated GMRES, one introduces \[ K_k(R)=\min_{p\in\mathcal P_k,\;p(0)=1}\max_{z\in R}|p(z)|, \] where \(R=[\alpha,\beta]+i[-\gamma,\gamma]\) contains the effective field of values. Then \[ \frac{\|r_k\|}{\|r_0\|} \le (1+\sqrt2)\,K_k(R). \] In the SPD case this reduces to the classical Chebyshev rate \[ K_k([\lambda_{\min},\lambda_{\max}]) = \Bigl(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\Bigr)^k, \qquad \kappa=\frac{\lambda_{\max}}{\lambda_{\min}}. \] For rectangles, explicit bounds include Elman-type, disk-segment, ellipse-based, conformal-mapping, and Faber-polynomial estimates [2504.05723]. Deflation modifies the polynomial problem by shrinking the relevant set. With the standard deflation projector \[ P_D=I-AU(U^*AU)^{-1}U^*, \] GMRES is applied to \(P_DA\), and the residual again has the form \[ r_k=p_k(P_DA)r_0, \qquad p_k(0)=1. \] If the deflation space captures the modes responsible for the large skew-Hermitian part, the containing rectangle narrows from \(\Omega_1=[\lambda_{\min},\lambda_{\max}]+i[-\rho,\rho]\) to \[ \Omega_1^\tau=[\lambda_{\min},\lambda_{\max}]+i[-\tau,\tau], \qquad \tau\ll\rho, \] and the corresponding \(K_k(\Omega_1^\tau)\) is markedly smaller [2504.05723]. ## 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 \[ \kappa z + M\overline z=b, \] the convergence bound is controlled by a polynomial approximation problem over a finite subset of the spectrum, now involving polynomials in \((\lambda,\overline\lambda)\). In the CSYM case, one obtains a discrete inner product \[ \langle p,q\rangle = \sum_{j=1}^n p(\lambda_j)\,\overline{q(\lambda_j)}\,r_j^2 \] and a family of orthogonal polynomials satisfying the three-term recurrence \[ \beta_{k+1}p_{k+1}(z) = z\,\overline{p_k(z)} - \alpha_{k+1}p_k(z) - \beta_k p_{k-1}(z), \] 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 [1111.5167]. In multigrid for Boltzmann transport, low-order GMRES polynomials are used as explicit approximations to \(A_{\rm ff}^{-1}\). For fixed order \(m\), one constructs \[ x^m=q_{m-1}(A_{\rm ff})\,b, \qquad p_m(\lambda)=1-\lambda q_{m-1}(\lambda), \] and then uses \(q_{m-1}(A_{\rm ff})\) to define \[ P=\begin{pmatrix} -\,q_{m-1}(A_{\rm ff})A_{\rm fc}\\ I \end{pmatrix}, \qquad R=\bigl[-A_{\rm cf}q_{m-1}(A_{\rm ff})\;\;I\bigr], \] as well as F-point relaxation updates. With fixed sparsity, the truncated polynomial approximation preserves the sparsity pattern of \(A_{\rm ff}\), limiting memory and communication [2301.05521]. 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 \(\sim1.3\) million dot products to converge or stalls, while PP-GMRES with \(\deg=10\) reduces this to \(\sim10.8\) 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 [1907.00072]. 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. Source: https://www.emergentmind.com/topics/gmres-polynomials