Repository navigation
solvers: drop the relative-residual guard, it fired on rounding noise (#21) - #53
Merged
Merged
Conversation
#21 suspected this check was unreachable: partialPivLu on an exactly singular jacobian divides by a zero pivot, so dx.allFinite() catches it first, and the backward-error test would only ever fire for a NEARLY singular matrix that still yields finite garbage. The proposal was to construct that case or drop the check. Measured, on the 2x2 system this adds, sweeping the perturbation: cond ~1e8 back_err/|R| = 3.8e-12 below linear_tolerance, no fire cond 1e12 back_err/|R| = 4.7e-08 fires; dx finite at 1e9 cond 1e14 back_err/|R| = 7.1e-06 fires; dx finite at 1e11 cond 1e15 back_err/|R| = 2.3e-05 fires; dx finite at 9e11 cond >=1e16 dx = inf/nan allFinite catches it first So the window is real and about four decades wide. Between roughly 1e12 and 1e15 the linear solve returns a finite answer that does not solve the system, and this check is the only thing that notices. The check stays. The test asserts more than "reports failure", because an unguarded solve also ends up not converged -- it just exhausts max_iter instead. What the guard has to do is decline to APPLY the step: with it removed, the iterate runs to -1.0010e+10 on the first update rather than staying at the seed. Both guards are now pinned and neither covers the other. Removing the backward-error check fails the new test while SingularJacobianReportsFailure still passes; removing allFinite fails SingularJacobianReportsFailure while the new test still passes. No production code changed. 301/301. Second-to-last item of #21.
5 tasks
Review of #53 found three wrong figures in my own comments and one test that could not fail. Fixing both, plus the production comment the review showed was half the story. The comment table was wrong in two ways. "Condition ~2/eps" is ~4/eps -- sigma max ~2, sigma min ~eps/2 -- so every condition number I quoted was understated fourfold. Worse, the load-bearing row was simply false: I wrote that beyond cond 1e16 dx goes non-finite and allFinite catches it. Measured, at eps = 2e-16 the condition number is 2.6e16, dx is FINITE at 4.5e12, and the backward-error check is what fires. The boundary between the two guards is not a condition number at all -- it is fl(1 + eps) == 1.0 making the pivot exactly zero, which happens between eps = 1.5e-16 and 1e-16. The table now includes that row and the eps = 1e-13 the test actually runs, which it previously omitted. "The iterate runs to -1.0010e+10 on the first update" mis-attributed a real number: -1.0010e10 is where it ends after 50 iterations; the first update gives -1.0008e10. The assertion is now EXPECT_EQ(x, 0.0) rather than a 1e3 bound -- the solve cold-starts and the guard fires on iteration zero, so the iterate is exactly the seed, and an approximate bound would also admit a damped step. The tests could not tell WHICH guard stopped the solve: allFinite and the backward-error check produce an identical observable. A new test raises linear_tolerance past the measured 5.5e-07 and changes nothing else, so the same step must now be APPLIED. Verified by mutation -- rejecting any large but finite dx, as allFinite would if it could, fails it. eps is a parameter now rather than buried in the material, so the window can be swept instead of asserted in prose, which is how the wrong numbers survived. The switch case also moved to match the enum order; near_singular was appended, so decoupled keeps its value and #52's "mode": 3 still means what it says. One test is NOT here, on purpose. I wrote one for the '* rn' scaling -- drop it and an absolute threshold would abort a healthy solve whose residual is in stress units -- and it was vacuous. Inside the solver, J = [1 1; 1 1+1e-8] with R = (-1e9, -1.001e9) gives a backward error of EXACTLY 0, so neither form fires and it passed either way; a standalone program with the same inputs gives 1.19e-07. Backward-stable LU lands on either side of zero by rounding. The gap is recorded in the test file rather than papered over with a test that cannot fail. The production comment claimed the two guards differ because "finiteness alone would not" catch the near-singular case. True, but it omits the converse, which is the sharper half: once dx is inf, back_err is NaN, and NaN > tol is FALSE, so the backward-error check would pass an infinite step. They are independent because of IEEE comparison semantics. Also recorded: -ffast-math defeats BOTH legs at once, folding J*dx - R to zero and making allFinite() true for inf. Finally, linear_tolerance's doc called it a backward-error tolerance. It is not, in the textbook sense -- that form is normalized by |J|*|dx|, and since partialPivLu is backward-stable it sits at ~1e-17 even at cond 1e15 and would never fire. Switching to it would silently disable the guard. Renamed to what it is: a relative-residual tolerance. 302/302.
A third review of #53 found that the guard's discrimination is scale dependent, and the previous comment invited the wrong reading of it. Measured on ONE system at fixed cond ~4e13, varying only the problem scale, unmutated code: scale rel. residual fires forward error of the step 1e0 5.5e-07 yes 8.0e-04 1e2 2.4e-07 yes 8.0e-04 1e3 1.1e-16 NO 8.0e-04 1e6 1.2e-16 NO 8.0e-04 The step is equally wrong at every scale. The check fires only where the constants leave sub-ULP noise in the residual -- at scale 1e3, 1.001*1000 is exactly representable and the noise disappears. So the guard is not detecting ill-conditioning and never was; it detects whether J*dx reproduces R, which is a different question that happens to coincide at some scales. Both the header comment and the test now say so, and the test records that its constants are load-bearing for that reason rather than arbitrary. Two smaller things from the same review, folded in because the file was open. The parse-time pair check prevents undefined behaviour, not merely a worse message: without it, names[1] on a one-element entry is an out-of-bounds vector read, and the ["x"] and [[]] cases abort the binary with heap corruption instead of failing cleanly. The comment framed it purely as diagnostics quality, which invites someone to relax it. The zero_blocks message assertion on "zero_blocks" is weak on its own -- json_reader_registry wraps every reader failure as "failed to read parameter 'zero_blocks'", so it passes for an nlohmann type error too. Only the "pair" assertion discriminates, and that one is deliberately coupled to the converter's wording; both now say as much at the assertion site. Also from #52: the JSON document now derives its mode from system_mode::decoupled rather than a literal 3, and the honest sub-case asserts converged() rather than only that nothing threw. 302/302.
Member
Author
|
Review findings (design, tmech expressions instead of component loops, modern C++, comments, tests, correctness). Severity: must = wrong result, hang, build break or unreachable feature; should = design, duplication, accuracy; nit = style. backward-error guard — needs changes
|
This was referenced Oct 17, 2026
The guard compared |J*dx - R| against linear_tolerance*|R|. A partial- pivoting LU is backward stable, so that residual is rounding noise: on the near-singular fixture it is 5.5e-7 without fused multiply-add and exactly 0 with it (-march=native on x86, every AArch64 build). The check fired on one toolchain and not the other, with the same Jacobian. Without it a near-singular solve still fails, by exhausting max_iter instead of stopping on the first step; #50 already keeps a failed iterate from seeding the next solve. allFinite() stays and catches the zero pivot of an exactly singular Jacobian deterministically. Removes the two round-off dependent tests added earlier on this branch, the near_singular fixture mode, and the parameter. Nothing outside the solver header set linear_tolerance.
petlenz
changed the title
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
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.
#21 asked whether the backward-error guard in
vector_newtonwas reachable. This PR first built a case where it fires, then found that the case is an artefact, and now removes the guard.What the guard measured
The check compared
|J*dx - R|againstlinear_tolerance * |R|.partialPivLuis backward stable, so that residual is rounding noise. On the near-singular fixture[1 1; 1 1+p]withp = 1e-13, row two evaluatesdx1 + (1+p)*dx2 - R2with|dx| ~ 1e10:-march=native(FMA)Without FMA the product
(1+p)*dx2rounds at a spacing of ~1.9e-6 and the cancellation againstdx1exposes that error. With FMA the product and sum are one operation with one rounding of an O(1) result, and the residual vanishes. The step is the correct solution of the ill-conditioned system either way; the guard saw a rounding mode, not a bad Jacobian.Consequence: the two tests added earlier on this branch passed on CI and failed 2/31 under
-march=native, and would fail on any AArch64 build. No value oflinear_tolerancechanges that: zero is below every positive threshold.What changes
linear_toleranceparameter are removed. Nothing outside the solver header set it.allFinite()stays. An exactly singular Jacobian divides by a zero pivot and leavesinf/nanindx; that is deterministic andSingularJacobianReportsFailurecovers it.near_singularfixture mode are gone.zero_blocksJSON tests stay: themodeliteral is spelled through the enum, and a solve that must be accepted is also asserted to converge.Without the guard a near-singular solve still fails, by exhausting
max_iterrather than stopping on the first step. #50 already keeps a failed iterate from seeding the next solve. A conditioning estimate (rcond()or a pivot ratio) was considered and rejected as not worth it for systems of this size.682/682 on main merged with this branch; 29/29 in
test_vector_newtonunder both the CI build and-march=native.