Repository navigation
fix(lsmr): underflowed norms, Givens update scale, LSMR's own residual estimate - #391
Merged
Merged
Conversation
This was referenced Sep 23, 2026
schroedk
added this pull request to stack #393
September 23, 2026 07:24
schroedk
force-pushed
the
fix/lsmr-numerics
branch
3 times, most recently
from
September 23, 2026 13:16
61f3f1f to
8ad42ac
Compare
schroedk
removed this pull request from stack #393
September 23, 2026 13:17
schroedk
added this pull request to stack #396
September 23, 2026 13:17
Every α and β in the stream is a norm read off an unscaled squared sum: the square vanishes at magnitudes the norm itself represents, and the exit read that zero as a verdict. `α₁ = √(p̃ᵀM⁻¹p̃)` on a design with `‖b‖ = 1e100` reported `x = 0` converged at zero iterations. `alpha_from_vp` falls back to a norm-scaled product, and the fused kernel `axpy_with_norm` returns `‖y‖` itself through `norm_from_sq`, which only pays the scaled pass when the sum left the normal range. The three `underflowed_alpha` cases now solve instead of certifying `x = 0`.
`SolutionState::update` floored its Givens diagonals at an absolute `f64::EPSILON`, on the claim that they are O(1). They scale with `‖A‖`: at `‖A‖ = 1e-9` the floor dropped the whole update and left `x = 0` describing a residual the recurrence reported as solved. Their product scales with `‖A‖²` and leaves the double range near `‖A‖ ≈ 1e±154`, where `ζ / (ρρ̄)` read as zero or NaN. Only an exactly vanished diagonal now ends a chain, and each factor divides by one diagonal in turn.
`residual_estimate` returned `|φ̄_k|`, the residual of LSQR's iterate on the same Krylov space. LSQR minimizes `‖r‖` there, so it understates the residual of the `x_k` LSMR returns (0.453 against a recomputed 0.479 on the Vandermonde fixture after two steps), and a residual-tolerance stop could fire before the iterate met the tolerance. The third rotation chain of Fong & Saunders §3.4, the one SciPy's `lsmr` runs, recovers `‖r_k‖` from scalars already in hand. It takes their sign for `ᾱ`, under which LSQR's `−c α` would flip `β̂` on every other step.
schroedk
force-pushed
the
fix/lsmr-numerics
branch
from
September 23, 2026 13:21
8ad42ac to
46db20c
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
1 participant
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.
Stack 3/6 (fixes), based on main (#394 merged). Independent of the exit rework in #388.
‖b‖ = 1e100reportedx = 0converged at 0 iterations. The fused axpy kernel returns‖y‖itself, so no step holds a squared norm.(x, h, h̄)update divides by one Givens diagonal at a time without an absoluteεfloor:‖A‖ ≈ 1e-9or beyond1e±154lost the update.‖r_k‖, from Fong & Saunders' third rotation sequenceResidualChain(Q̃, §3.4), instead of LSQR's|φ̄|, which understates it, so a residual-tolerance stop could fire before the iterate met the tolerance.Each fix's test fails with the fix reverted. The third fix moves residual-tolerance stops and has not been benchmarked on its own.