Repository navigation
Fix the Chebyshev iteration recurrence - #376
Merged
haampie merged 1 commit intoSep 14, 2026
Merged
Conversation
After JuliaLinearAlgebra#367 the iteration converges, but it still does not produce the Chebyshev residual polynomial r_k = T_k((μI - A)/ρ) r_0 / T_k(μ/ρ): - the iteration counter starts at 0, so the `iteration == 1` branch ran on the second step instead of the first, resetting the direction and using α = 2/μ there; - `α = 1/(μ - β)` mixed a dimensionless β with an eigenvalue; the Templates form is α = 1/(μ - β/α_prev); - the β used on the second step must be half the generic value, since the three-term recurrence starts from T_{-1} = T_1 rather than T_0. Rewrite `iterate` in the coupled two-term form of Sleijpen's lecture notes (Exercise 5.11), tracking ζ_k = T_k(μ/ρ) / T_{k-1}(μ/ρ) with ζ_{k+1} = 2μ/ρ - 1/ζ_k, β_k = 1/ζ_k² and α_k = 2/(ρ ζ_{k+1}). On a diagonal SPD matrix with spectrum in [1, 100] the residual norm after 60 iterations is now 1.2e-4 (the exact Chebyshev value), where the previous code gave 2.2e-1 and, before JuliaLinearAlgebra#367, diverged to 4e10. The existing tests did not catch this because their matrix has condition number below 3. Add a test that checks the residual history against the exact Chebyshev polynomial evaluated per eigenvalue, with and without a diagonal left preconditioner, for real and complex element types. Also fix the `maxiter` docstring, which mentioned GMRES. Fixes JuliaLinearAlgebra#375.
haampie-bot
force-pushed
the
fix-chebyshev-recurrence
branch
from
September 14, 2026 13:03
bdaac24 to
51c6e44
Compare
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.
Follow-up to #367 / #375. That fix stops the divergence, but the iteration still does not produce the Chebyshev residual polynomial
r_k = T_k((μI - A)/ρ) r_0 / T_k(μ/ρ):iteration == 1branch ran on the second step, resetting the direction and usingα = 2/μ;α = 1/(μ - β)mixes a dimensionlessβwith an eigenvalue (Templates has1/(μ - β/α_prev));βon the second step must be half the generic value, because the three-term recurrence starts fromT_{-1} = T_1, notT_0.iterateis rewritten in the coupled two-term form (Sleijpen, lecture notes Ex. 5.11): trackζ_k = T_k(μ/ρ)/T_{k-1}(μ/ρ)withζ_{k+1} = 2μ/ρ - 1/ζ_k,β_k = 1/ζ_k²,α_k = 2/(ρ ζ_{k+1}).Residual norm on
Diagonal(range(1, 100, length=200)),b = ones:The existing tests did not catch any of this because their matrix has condition number < 3. The new test compares the residual history against the exact polynomial evaluated per eigenvalue (diagonal
A, with and without a diagonal left preconditioner, real and complex eltypes); it fails on current master and passes here to ~1e-15 relative.Also fixes the
maxiterdocstring, which mentioned GMRES.