Skip to content

Solve the Newton step without forming the normal equations. Fixes the original report and mesr's assembly from #1354 - #4

Open
BoykoNeov wants to merge 2 commits into
masterfrom
fix-1354-least-squares
Open

Solve the Newton step without forming the normal equations. Fixes the original report and mesr's assembly from #1354#4
BoykoNeov wants to merge 2 commits into
masterfrom
fix-1354-least-squares

Conversation

@BoykoNeov

@BoykoNeov BoykoNeov commented Jul 27, 2026

Copy link
Copy Markdown
Owner

Fixes the original report and @mesr's assembly from solvespace#1354 — sketches that solve fine at small sizes and report "incompatible constraints" once the dimensions get large.

It is not the redundancy test either

Same dead end as solvespace#1247, and worth stating so nobody re-walks it: rankOk is true, and the rank/redundancy branch is not what rejects these sketches. The rank decision that bites is the one inside the Newton step solve.

Root cause: the normal equations square the condition number

System::SolveLeastSquares() computed the minimum-norm step as x = A'(AA')⁻¹B, forming AA' explicitly and factoring it with Eigen's rank-revealing SparseQR.

Forming AA' squares the Jacobian's condition number. Eigen's pivot cutoff is 20(m+n)·maxColNorm·ε, and applying it to AA' makes it behave, as far as A is concerned, as though the tolerance were ~√ε ≈ 2e-7 rather than ~ε. SparseQR::_solve_impl then zeroes the component of the step it considers rank-deficient — so one residual can never be driven to zero, however many iterations run, and the solver reports a set of entirely compatible constraints as incompatible.

The transition is exactly as sharp as that implies. With one length pinned at 4 m and the other swept:

lenA 5th pivot threshold rank result
2880 mm 3.657e-07 3.65574e-07 5/5 converges
2881 mm below 3.65574e-07 4/5 residual pinned at 0.0190797, ‖X‖ → 1e-13

A pivot landing 0.03 % from the cutoff is also why the threshold moves by orders of magnitude when unrelated entities are deleted, which is what @ruevs observed on the issue.

The fix

Factor A' directly: A' = QR, hence AA' = P R'R P', so z = P R⁻¹R'⁻¹P'B is two triangular solves and the rank decision is made on pivots that scale like A instead of like .

  • m > n keeps the old normal-equations path — Eigen's sparse QR wants tall matrices — and is byte-identical to master there.
  • Zero-size systems now return a zero step instead of multiplying an uninitialised vector.
  • SolveLinearSystem, TestRank, CalculateRank and NewtonSolve are untouched.

The threshold moves from 2.88 m to 76.9 km.

Redundancy detection gets better, not worse: models with genuinely duplicate constraints degrade from REDUNDANT_OKAY to REDUNDANT_DIDNT_CONVERGE above ~4 m on master, and with this change stay REDUNDANT_OKAY with correct candidate lists out to 100 m.

Regression test

test/constraint/large_dimensions/perpendicular_4m — two 4 m lines, perpendicular, with a point-on-line constraint. It checks the geometry (both magnitudes 4000.0, dot product 0, point-on-line distance 0), not just that a result code came back OKAY.

Please do not canonicalize this fixture. It is @mesr's file exactly as uploaded, and it has to stay that way: re-saving from SolveSpace writes out the solved parameter values, so the sketch loads already at its solution, the first Newton step is the zero step, and the bug cannot occur. That is the opposite of this suite's usual convention, so there's a comment saying so next to the CHECK_LOAD as well as here. It is also why @mesr's assembly needs the part re-saved to reproduce — linked files load already-solved entities — and the second commit adds a resave command to test/debugtool.cpp for exactly that.

Verification

  • Full suite passes in Debug and Release: 263 cases / 932 checks (master is 262 / 925). I checked Release specifically because this failure turns on a pivot landing 0.03 % from a threshold — precisely the kind of thing that can pass /Od and fail /O2.
  • With src/system.cpp reverted, the new test fails at IsSolvedOkay().

Scope, and what this does not fix

Per @ruevs, take this one after the constraint-deletion fix (#2, upstream solvespace#1744). They are independent — either order builds and passes — but that one is the more severe bug and the smaller diff.

Please fetch and fast-forward this branch rather than using the merge button on my fork — that keeps it a clean fast-forward for upstream.

🤖 Generated with Claude Code

BoykoNeov and others added 2 commits July 27, 2026 23:38
SolveLeastSquares() computed the minimum norm Newton step as
x = A'*(A*A')^-1*B, forming A*A' explicitly and factoring it with a
rank-revealing sparse QR. That squares the condition number of the
Jacobian. Our equations mix dimensionless quantities with lengths and
with areas, so the spread of magnitudes in A already grows with the
physical size of the sketch; squaring it pushes the smallest pivot of
A*A' below the threshold Eigen uses to call a column linearly
dependent, which is proportional to the largest column norm. Eigen's
rank-truncated solve then silently zeroes that component of the step,
so the residual of one equation can never be driven to zero, Newton's
method stalls, and a perfectly solvable sketch is reported as having
incompatible constraints.

On the file from the bug report, two 4 m lines constrained
perpendicular with a point on line, the transition is exact: with the
second line pinned at 4 m and the first at 2880 mm, the fifth pivot of
A*A' is 3.657e-7 against a threshold of 3.6557e-7 and the sketch
solves; at 2881 mm the pivot falls just below the threshold, the rank
drops from 5 to 4, and the solve fails, with the step norm collapsing
to 1e-13 while one residual stays pinned at 0.019.

Factor A' = Q*R directly instead. Then A*A' = P*R'*R*P', so the same
z = (A*A')^-1*B comes from two triangular solves against R, and the
rank decision is taken on pivots that scale like A rather than like
A*A'. In that same sweep, the longest first line that solves goes from
2.88 m to 76.9 km.

Rank determination for redundant constraints is untouched; it runs in
TestRank(), against A itself. Over-constrained sketches are in fact
reported better than before: duplicating a constraint in that file
used to be reported as redundant below about 3 m but as unsolvable
above it, and is now reported as redundant at every size, out to 100 m.

The wide case, more equations than unknowns, keeps using the normal
equations, since Eigen's sparse QR wants a matrix that is at least as
tall as it is wide. Such a system is redundant anyway, and it gives the
same result as before. A system with no equations or no unknowns now
returns a zero step, rather than multiplying an uninitialized vector by
a matrix of a mismatched size.

Also teach the debug tool to load a file and report how each group
solved, and to load a file and save it back, which is how a linked part
gets re-solved without a GUI.

Fixes the original report and mesr's assembly from solvespace#1354.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The fixture is the file from the bug report, unmodified: two 4 m lines
constrained perpendicular, with the start of the second one on the
first. It has to be the file as saved, not a canonical re-save, because
the failure depends on the stored parameter values being the ones from
before the perpendicular constraint was added; re-saving it stores the
solution, and then there is no Newton step left to take.

The test checks the geometry as well as the solve result, since the
constraints could be satisfied by more than one configuration.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@ruevs

ruevs commented Jul 28, 2026

Copy link
Copy Markdown

bug_angle_dim.zip from here is fixed by this PR.

issue.zip from here is also fixed.

If when we decide to merge this is should be AFTER #2.

@@Evil-Spirit @jwesthues @phkahler - please take a look - my linear algebra it too rusty to evaluate this.

@BoykoNeov

Copy link
Copy Markdown
Owner Author

Written by Claude Opus 5 — both this text and the code it describes; posted by @BoykoNeov.

@ruevs Thanks for testing both of those — good news, and I've corrected the PR body.

The "not claimed" bullet was mine and it was based on a headless run that could not see the symptom: bug_angle_dim.zip is under-constrained (dof = 8), so an unenforced angle just leaves the sketch sitting somewhere else — the solver returns OKAY at every length either way, with or without this change. Nothing to observe without the GUI. Your screenshots settle it.

The commit trailer currently reads "Fixes the original report and @mesr's assembly from solvespace#1354" precisely because I couldn't claim the rest. Now that you've confirmed it, say the word and I'll amend it to Fixes #1354 so that merging closes the issue.

Merge order acknowledged — #2 (upstream solvespace#1744) first. That branch now has a second commit fixing the leak its test uncovered; the details are over there, and the short version is that the leak is pre-existing and independent of the fix, measured both ways rather than assumed.

For @Evil-Spirit / @jwesthues / @phkahler, the one-paragraph version of the linear algebra:

The Newton step is under-determined (m equations, n > m unknowns) and we want the minimum-norm solution, which is x = A'(AA')⁻¹B. SolveLeastSquares() did exactly that literally: build AA', hand it to Eigen's rank-revealing sparse QR. But cond(AA') = cond(A)², and Eigen's pivot threshold 20(m+n)·maxColNorm·ε is applied to the squared matrix — so in terms of A the effective tolerance is about √ε ≈ 2e-7 instead of ε ≈ 2e-16. Past a certain sketch size a perfectly good pivot falls under it, SparseQR::_solve_impl zeroes that component of the step, and one residual can never be driven to zero: "incompatible constraints" on a sketch that is not. This change factors A' = QR directly and gets the same x from two triangular solves — AA' = P R'R P', so z = P R⁻¹ R'⁻¹ P' B — with the rank decision now made on pivots that scale like A rather than like . Same answer in exact arithmetic; the size at which it breaks moves from 2.88 m to 76.9 km. Nothing else in the rank/redundancy machinery is touched.

The sharpest evidence that this is the mechanism, rather than a plausible story: with the second dimension pinned at 4 m, lenA = 2880 gives a 5th pivot of 3.657e-07 against a threshold of 3.65574e-07 — rank 5 of 5, solves — and lenA = 2881 puts it under, rank 4 of 5, residual pinned at 0.0190797. One millimetre either side of the cliff.

@BoykoNeov

Copy link
Copy Markdown
Owner Author

Written by Claude Opus 5 — both this text and the code it describes; posted by @BoykoNeov.

@ruevs Thanks — that settles the one thing I could not settle myself.

bug_angle_dim.zip. I could not reproduce a failure with it headlessly: it reports OKAY at every length I tried, dof=8, and the symptom I saw was a silently unenforced angle rather than an error. That is why the commit trailer names the original report and @mesr's assembly instead of saying "Fixes solvespace#1354". If what you are seeing is the GUI behaviour — where the previous solution feeds the next solve's initial guess, so the failure depends on how you got there — then that is consistent with my measurement rather than contradicting it, and I am happy to take it as fixed. Confirm and I will amend the trailer to Fixes #1354. That means force-pushing b473a13f, so I would rather you say the word than have the head move under a branch you have already read.

issue.zip from #5080650753 I have not run at all. I can verify it before the trailer changes if you would prefer the claim to rest on something other than my word for it.

Merge order noted — after #2, which is now upstream solvespace#1744 and green, including the ASan leak (fixed in 81f473ff; it was a pre-existing ownership bug in Group::Clear(), not the tag clearing).

One more while you are here: of these PRs, #3 (issue solvespace#1247) is the only one nobody has looked at. It is not linear algebra — a backtracking line search whose fallback is "if no damped step reduces the residual, take the full step exactly as before" — so it is a much smaller thing to check than this one. I have left a short note on it.

@ruevs

ruevs commented Jul 29, 2026

Copy link
Copy Markdown

Don't bother with testing issue.zip - I've done it already and it works.

This will wait until @jwesthues or @Evil-Spirit reviews it.

You can amend the commit message to say "Fixes solvespace#1354" if you want; but don't rebase - #3 which is now solvespace#1748 will probably get merged first and other changes on master may happen.

@ruevs

ruevs commented Aug 4, 2026

Copy link
Copy Markdown

@phkahler

phkahler commented Aug 4, 2026

Copy link
Copy Markdown

@ruevs I followed all the links in these solver related PRs and wrote down all the issues referenced. I don't know if they'll be be fixed, but the numbers were 226, 1105, 246, 1354, and 1247.

@ruevs

ruevs commented Aug 4, 2026

Copy link
Copy Markdown

@ruevs I followed all the links in these solver related PRs and wrote down all the issues referenced. I don't know if they'll be be fixed, but the numbers were 226, 1105, 246, 1354, and 1247.

@phkahler I tested the models from all the solver issues.

Here is the status:

After/if I merge #1749 (shall I do it?) - I'll rebase this and open a PR for it in the main repository. But let's finish the NURBS fix first solvespace#1746.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants