[leburgel]

Note

Work in progress. This draft has not been reviewed in full yet and may still change. It is opened early to run CI.

This adds deflated restarting, GMRES-DR(m, k) (R. B. Morgan, SIAM J. Sci. Comput. 24, 20 (2002)), to GMRES as a new option ndeflate (Morgan's k). The default ndeflate = 0 keeps the current behavior.

Restarted GMRES(m) discards the whole Krylov subspace at every restart. When the linear map has a few eigenvalues close to zero, every cycle has to rediscover them, and the residual can stagnate for thousands of applications. GMRES-DR keeps the k harmonic Ritz vectors of smallest harmonic Ritz value across the restart, together with the residual. Their span is again a Krylov subspace, so the next cycle continues from an Arnoldi factorization of length k. This needs no extra application of the map and the same m + 1 stored vectors.

Motivation

PEPSKit's implicit CTMRG gradients solve a real linear problem with GMRES(; krylovdim = 30), which stagnates at critical states. The table gives applications of the map to a relative residual of 1e-6 / 1e-8 for the operator of the critical Ising partition function, assembled densely (n = 1280 at χ = 16, n = 5120 at χ = 32). † means no convergence within 5000 applications.

χ	GMRES(; krylovdim = 30)	GMRES(; krylovdim = 30, ndeflate = 10)
16	† / †	94 / 118
32	† / †	186 / 240

Solves that converge within the first cycle are unchanged.

Demonstration

A 1000 × 1000 matrix with four eigenvalues of order 1e-3 and the rest in [1, 10]. Restarted GMRES(30) does not converge in 100 cycles, while GMRES-DR(30, 10) converges in 80 applications.

Script

using KrylovKit, LinearAlgebra, Random, Printf

# eigenvalues: four close to zero (one complex conjugate pair), the rest in [1, 10];
# mildly non-normal and rotated to a random orthonormal basis
rng = MersenneTwister(42)
N = 1000
D = diagm(range(1.0, 10.0; length = N))
D[1, 1], D[2, 2] = 1.0e-3, 2.0e-3
D[3:4, 3:4] = [1 -2; 2 1] * 1.0e-3
D += triu(randn(rng, N, N), 1) / (2 * sqrt(N))
Q = Matrix(qr(randn(rng, N, N)).Q)
A = Q * D * Q'
b = randn(rng, N)

tol = 1.0e-8 * norm(b)
for (name, alg) in (
        ("GMRES(30)", GMRES(; krylovdim = 30, maxiter = 100, tol, verbosity = 0)),
        ("GMRES-DR(30, 10)", GMRES(; krylovdim = 30, ndeflate = 10, maxiter = 100, tol, verbosity = 0)),
    )
    x, info = linsolve(A, b, zero(b), alg)
    @printf(
        "%-17s converged = %d, applications = %4d, residual = %.1e\n",
        name, info.converged, info.numops, norm(b - A * x) / norm(b)
    )
end

GMRES(30)         converged = 0, applications = 3002, residual = 3.6e-02
GMRES-DR(30, 10)  converged = 1, applications =   80, residual = 7.4e-09

Interface

x, info = linsolve(f, b, x₀, GMRES(; krylovdim = 30, ndeflate = 10, tol), a₀, a₁)
x, info = linsolve(f, b, x₀, a₀, a₁; krylovdim = 30, ndeflate = 10, rtol)

• Constraint: 0 <= ndeflate < krylovdim, checked in the keyword constructor. The docstring suggests div(krylovdim, 3).
• Cost per cycle: maxiter still counts cycles. With ndeflate > 0, every cycle after the first applies the map krylovdim - ndeflate times.
• Coverage: real and complex problems, the a₀, a₁ form, reallinsolve and every Orthogonalizer. Since it is still a GMRES, it is also accepted as alg_rrule by the eigsolve and svdsolve pullbacks (not tested with ndeflate > 0). Log messages read "GMRES-DR" when ndeflate > 0.
• Scalar type: deflated restarts need a BlasFloat scalar type.
• Compatibility: ndeflate is appended as the last field. The five-argument positional constructors GMRES(orth, maxiter, krylovdim, tol, verbosity) and GMRES{O, S}(...) are kept and set ndeflate = 0, and the type parameters are unchanged. The only incompatibility is with GMRES objects serialized by an older version.

Implementation

The algorithm changes are in src/linsolve/gmres.jl. The new field, constructors and docstring are in src/algorithms.jl, and ndeflate is passed through the keyword interface in src/linsolve/linsolve.jl. Within a cycle the code is that of GMRES, with the update of a column of R moved into a helper, _gmres_hessenbergcolumn!. The cycle can now start from a factorization of length k with a general residual vector instead of β e₁. When ndeflate > 0, a full and nonsingular cycle ends in _gmresdr_deflate!, and otherwise the restart is the one of GMRES. With G = a₀ + a₁ H and g = a₁ h_{m+1,m}:

1. Kept vectors: the harmonic Ritz values are the eigenvalues of the Hessenberg matrix G + |g|² G⁻ᴴ eₘ eₘᵀ. Its Schur form (hschur!, reordered with partitionschur!) gives an orthonormal basis for the ndeflate values of smallest modulus. In real arithmetic, complex conjugate pairs are kept together.
2. Residual direction: the new residual direction is [-conj(g) G⁻ᴴ eₘ; 1], orthogonalized against that basis. It spans the orthogonal complement of the range of the (m+1) × m Hessenberg matrix, so it contains the GMRES residual.
3. Arnoldi form: the projected Hessenberg matrix is brought back to Arnoldi form with _restorearnoldiform!, as in the Krylov–Schur eigsolve. The basis is updated with basistransform! and the factorization with shrink!.
4. Residual coordinates: these carry over from the Givens rotations, so no application of the map is needed.

No dependency is added. The restart falls back to the one of GMRES if G is singular or no vector can be kept, and if the implicit residual reaches tol but the recomputed one does not.

Tests

New testsets in test/linsolve.jl:

• Constructors: backward-compatible positional forms and argument checks.
• Full factorization: all four scalar types and vector modes, logs, reallinsolve and the a₀, a₁ form. Without a restart the result and numops equal those with ndeflate = 0.
• With restarts: a non-normal 100 × 100 matrix with four eigenvalues of order 1e-2. GMRES(10) does not converge in 50 cycles, while GMRES-DR(10, 4) converges in about 80 applications. This testset also covers the keyword interface, complex shifts, and reallinsolve (including a real-linear map with conj).

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 98.07692% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 88.29%. Comparing base (80dc683) to head (f31a682).

Files with missing lines	Patch %	Lines
src/linsolve/gmres.jl	97.97%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #178      +/-   ##
==========================================
+ Coverage   86.63%   88.29%   +1.66%     
==========================================
  Files          36       36              
  Lines        3920     3989      +69     
==========================================
+ Hits         3396     3522     +126     
+ Misses        524      467      -57     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.