Conversation
Implement GMRES with deflated restarting (Morgan, 2002) for the MATLAB
toolbox. The solver deflates k harmonic Ritz vectors into the starting
basis of each restart cycle, keeping the subspace dimension fixed at m
(unlike GMRES-E which uses m+d). The first k columns of the Hessenberg
matrix at each restart are formed without new matrix-vector products by
reusing the previous cycle's Arnoldi factorisation.
Files added:
matlab/src/solvers/gmres_dr.m -- solver with detailed inline docs
matlab/tests/test_gmres_dr.m -- input validation, fallback dispatch,
identity/diagonal, Embree 3x3,
Sherman1 and Sherman4 sparse tests
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Interactive script that runs all four KrySBAS solvers (GMRES-E, LGMRES, PD-GMRES, GMRES-DR) on a user-configurable set of test matrices and produces one figure per matrix showing log(relative residual) vs restart cycle with a distinct colour/line style per solver. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
- Add GMRES(m) (built-in) as a grey reference curve. Its per-inner- iteration resvec is downsampled to one point per restart cycle via unique([1:m:end, end]) to handle early convergence correctly. - Build all legend labels dynamically with solver parameters and CPU time using LaTeX strings (e.g. GMRES-DR (m=27, k=3, t=0.12 s)). - Set Interpreter latex on legend, xlabel, ylabel, and title. - Move legend to eastoutside; ylabel now shows ||r||/||r_0|| in math. - Expose alpha_p / alpha_d in CONFIGURATION and pass to pd_gmres so the legend reflects the actual PD coefficients used. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Set explicit figure size (1150x520 px) so the eastoutside legend has enough horizontal room to display full parameter + timing labels. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Labels are too wide for the eastoutside box. Placing the legend below the plot (southoutside, NumColumns=2) gives each entry half the figure width, which is enough for even the PD-GMRES label with all parameters. Figure height increased to 620px to accommodate the taller layout. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
MATLAB's built-in gmres() returns absolute residual norms in resvec. Divide by resvec(1) so GMRES(m) is plotted on the same relative scale as the KrySBAS solvers (all starting at 1.0). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
eigs() can return complex eigenvectors for real matrices when harmonic Ritz values come in conjugate pairs. For a real linear system the imaginary components are pure numerical noise. Stripping them with real() after the sort step (both in cycle 1 and in the restart loop) keeps Ek, Yk, Vk, and the solution x real throughout. Fixes test_embree_3x3_toy_example (x had O(1e-5) imaginary parts). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
GMRES-DR(m=2, k=1) uses a 2-dimensional subspace per cycle on this 3x3 system, so it needs multiple restarts to converge. GMRES-E(m=2, d=1) uses m+d=3=n in its first cycle (full space), giving an exact solution in one shot. The accumulated restarts leave GMRES-DR's solution accurate to ~1e-6 (the solver tol), which is outside assertElementsAlmostEqual's default sqrt(eps) ~ 1.49e-8 tolerance. Tightening tol to 1e-12 ensures the solver runs until the solution error is well within that threshold. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
GMRES-DR(m=2, k=1) uses a 2-dimensional subspace per restart cycle. For the 3x3 Embree matrix this requires multiple restarts, and the solution converges to the solver tolerance (1e-6), not to machine precision. Revert tol to 1e-6 and pass an explicit 'relative' 1e-4 tolerance to assertElementsAlmostEqual so the test checks correctness at the accuracy the solver actually achieves. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Collaborator
Author
|
@gusespinola @jccf19 : The PR is ready for code-review. As soon as this is merged, I will work on the Julia version of the solver. |
The generalized eigenvalue problem G*Ek=lambda*Fsq*Ek can have singular
G when the hybrid Hessenberg (free columns + Arnoldi) is ill-conditioned
in restart cycles. The fix uses a dense/sparse switcher controlled by
eig_dense_thresh (default 200, the subspace size s -- always s<=m,
never n):
s <= 200 : eig() -- LAPACK, all eigenvalues, robust against singular G
(returns Inf for degenerate directions, filtered by helper)
s > 200 : eigs() -- ARPACK, k eigenvalues, preferred for large subspaces;
falls back to eig() via try-catch if eigs fails
A local helper gmres_dr_pick_ritz selects the k smallest-|lambda| finite
eigenvectors, with a further fallback to eig(F) (non-generalized) if
fewer than k finite pairs survive after filtering.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
jccf19
approved these changes
Jun 20, 2026
ROOT CAUSE: After enough restart cycles, harmonic Ritz values converge toward true eigenvalues. When two Ritz values are nearly equal (clustered spectrum), their Ritz vectors become almost parallel. The QR of Yk then produces a near-zero diagonal in Rk_qr, making Ek_norm = Ek / Rk_qr blow up. This corrupts T = Hprev*Ek_norm, H_def, and the subsequent least-squares solve and eigs call. FIX (rank-revealing truncation in Step 1): After [Vk, Rk_qr] = qr(Yk, 0), count k_eff = the number of diagonals of Rk_qr that exceed max(diag)*sqrt(eps). Retain only those k_eff columns; the remaining k - k_eff are numerically dependent on the accepted set and are silently dropped for this cycle. The Arnoldi extension runs for m - k_eff steps so the total subspace stays at m. All k references inside the cycle body (hgs, H_def indexing, V_cycle assembly, phi, Arnoldi loop start) become k_eff. Step 8 still requests k Ritz vectors for the NEXT cycle so k_eff is recomputed fresh each restart. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
The previous fix used plain qr(Yk, 0) (no pivoting) and then truncated to the first k_eff columns. This is wrong: without column pivoting a near-zero diagonal can appear anywhere in Rk_qr, so Rk_qr(1:k_eff,1:k_eff) can still be singular even after truncation, causing Ek/Rk_qr to blow up. Replace with column-pivoting QR: [Vk, Rk_qr, Pmat] = qr(Yk, 0). Column pivoting guarantees |Rk_qr(1,1)| >= |Rk_qr(2,2)| >= ... >= |Rk_qr(k,k)|, so keeping the first k_eff columns (where k_eff = number of diagonals above diag_R(1)*sqrt(eps)) always selects the most independent subset and leaves Rk_qr(1:k_eff,1:k_eff) well-conditioned. Ek_norm is updated accordingly to account for the column permutation: Ek_norm = Ek*Pmat(:,1:k_eff)/Rk_qr. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
…mpat In Octave, [Q, R, P] = qr(A, 0) returns P as a permutation VECTOR when the 0 flag is passed, not as a permutation matrix. Using Pmat(:, 1:k_eff) on a 1×k row vector gives a 1×k_eff row vector, so Ek * (1×k_eff) fails with a dimension mismatch. In MATLAB the same call returns a matrix. Fix: immediately after the QR call, normalize Pmat to a row vector: - if isvector(Pmat): already a vector, ensure row orientation - else: extract permutation vector from the matrix form Then use column indexing Ek(:, perm(1:k_eff)) instead of the matrix multiply Ek * Pmat(:, 1:k_eff), which works identically in both runtimes. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Replace the free-column formula and eigs(F,G) approach with two more
numerically robust techniques from the reference GMRES-DR implementation:
1. INCREMENTAL QR (qrupdate_gs, new utility):
The thin QR of the growing Hessenberg H is updated column-by-column
via double-pass modified Gram-Schmidt. This gives a cheap residual-
norm estimate at every inner step and avoids recomputing QR from
scratch each cycle.
2. SCHUR-BASED HARMONIC RITZ:
The harmonic Ritz problem is solved via the explicit harmonic matrix
Ht = H_sq + h_{m+1,m}^2 * (H_sq' \ e_m) * e_m'
followed by schur() + ordschur() reordering. This replaces the
generalized eigenproblem eigs(F,G,k,'LM') which failed whenever
G = Rs'*Rs was ill-conditioned. Complex conjugate pairs in the real
Schur form are handled by setting keep = k+1 when T(k+1,k) ~= 0,
avoiding the need to force real() on complex eigenvectors.
3. THICK RESTART via orthogonal basis rotation (Wu & Simon 1999):
At each restart, Pk (the k Schur vectors) and rc (residual direction)
are combined into an orthonormal Pkp1 via thin QR with sign correction.
H and V are then rotated: H <- Pkp1'*H*Pk, V <- V*Pkp1, vr <- Pkp1'*rc.
This pure orthogonal transformation replaces the error-prone free-column
formula and its associated state (Hprev, g_res, Ek_norm, g_tilde).
Files:
matlab/src/utils/qrupdate_gs.m new incremental QR utility
matlab/src/solvers/gmres_dr.m full rewrite
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
When the Hessenberg H gains one row between consecutive calls (e.g. grows from (j)x(j-1) to (j+1)xj as the Arnoldi loop progresses), q_in has the old row count while n_rows = j+1. The previous q_out(:, 1:n_prev) = q_in assignment failed because of the row-count mismatch. Use n_q_rows = size(q_in, 1) and copy only the existing rows; the new bottom row stays zero, which is correct since the new Arnoldi vector is orthogonal to the existing Q columns by construction. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Add reference [3] in gmres_dr.m and a References section in qrupdate_gs.m crediting the GMRES-SDR repository by the NLA Group, University of Manchester (https://github.com/nla-group/GMRES-SDR), which is the source of the Schur-based thick restart and incremental QR update adopted in this implementation. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Eight tests across four categories:
Full one-shot factorisation:
- correctness (A = Q*R, Q'*Q = I, R upper triangular)
- positive diagonals of R
Incremental column-by-column update:
- incremental build matches full one-shot call
- orthogonality of Q is preserved at every step
Growing row count (the GMRES-DR use case):
- (j+1)×j submatrix passed at step j: verifies the row-growth fix
(q_in has one fewer row than the new matrix)
- least-squares residual from incremental QR matches pinv reference
Edge case:
- single-column input: Q is a unit vector, R equals column norm
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Primary path: eigs(F, G, k, 'LM') -- matches GMRES-E, fast and accurate
when G = R'*R is well-conditioned (rcond > 1e-14). The k eigenvectors
are sorted by ascending |lambda|, stripped of complex noise via real(),
and QR-orthonormalised to form Pk for the thick restart.
Fallback path: when rcond(G) < 1e-14 the generalized eigenproblem is
ill-conditioned. We compute the explicit harmonic matrix
Ht = H_sq + h_{m+1,m}^2 * (H_sq' \ e_m) * e_m'
and use the real Schur decomposition + ordschur, which is always
numerically stable. The 2x2 block check (keep = k+1) applies only to
this path since eigs already returns real arithmetic.
G is free to compute (reuses r_qr from the Arnoldi QR already done this
cycle). The rcond threshold 1e-14 was suggested by J.C. Cabral.
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
Add an optional 6th output argument stats (only populated when requested
via nargout >= 6):
[x, flag, relresvec, kdvec, time, stats] = gmres_dr(...)
stats.n_eigs : number of restart cycles that used the eigs primary
path (rcond(G) > 1e-14 -- well-conditioned case)
stats.n_schur : number of restart cycles that fell back to the Schur
decomposition (rcond(G) <= 1e-14 -- ill-conditioned G)
Counters are initialised to 0 at algorithm setup, incremented in the
respective if/else branches, and assembled into the stats struct at
every return point (including the fallback dispatches, where both
counts are 0 since no harmonic Ritz step is performed).
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
compare_solvers.m: - Add clear/clc at start for a clean workspace - Switch default matrix to sherman5 with m=50, maxit=600 - Add isfield guard for Problem.b (some matrices have no b field) - Request stats output from gmres_dr to expose eigs/Schur counters - Use m_dr = m + k + d for GMRES-DR to match GMRES-E subspace size data/: - Add ex40.mat and wang4.mat (additional SuiteSparse test matrices) Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
MATLAB Implementation of GMRES-DR as in https://epubs.siam.org/doi/10.1137/S1064827599364659 with its corresponding unit tests.
Also, added a script
compare_solvers.mthat allows for easy comparison of solvers.