From dbcb2a8ea8084b3ade14482b57969e922fc0bc9e Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 10 Oct 2026 08:51:18 +0200 Subject: [PATCH] CultMath BoundedLeastSquares: release on the gradient against its rounding, never by the cost A bound is released only when its violation exceeds both kktTol and the rounding bound of its double gradient (Higham gamma_(n+1)); the cost-progress rule now judges only passes that release nothing, and a release resets fMark. Fixes a cold start stopping Converged with a violated bound (Soul W5, finding aetheria-release:finding:cut-thrust-allocator-core.s1.solver-converged-early). Tests: the W5 problem (cold == re-solve at 1e-12 and 1e-9) and the seed-141 noise case. Builds; verified on Yggdrasil in the following report. Co-Authored-By: Claude Sonnet 5.5 --- .../src/CultMath/BoundedLeastSquares.cs | 41 +++++++++---- .../BoundedLeastSquaresTests.cs | 61 +++++++++++++++++++ 2 files changed, 89 insertions(+), 13 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index 0f08ec5c..2d31620e 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -34,11 +34,13 @@ public enum BoundedLeastSquaresStatus /// Stopping. KKT termination uses a tolerance of (1e-6) times the problem's /// gradient scale (the largest gradient component at the box point nearest the origin), so scaling A, b or the /// bounds, and the choice of warm start, do not change it. A caller whose rows span many decades of weight can pass -/// a smaller kktRelativeTolerance, which changes only this tolerance: the progress rule below and the pivot -/// threshold are fixed. The tolerance is only a resolution: once a stationarity miss or a -/// bound release has been followed by a pass that lowers the cost by no more than rounding, the solver reports -/// , so an optimal point whose double gradient sits above the -/// tolerance costs one confirming iteration rather than the iteration cap. +/// a smaller kktRelativeTolerance, which changes only this tolerance: the rounding bound and progress rule +/// below and the pivot threshold are fixed. A bound is released only when its violation exceeds both the tolerance +/// and the rounding bound of its double gradient, so a tolerance below the noise cannot cycle on a release. The +/// tolerance is otherwise only a resolution: once a stationarity miss that releases nothing has been followed by a +/// pass that lowers the cost by no more than rounding, the solver reports +/// , so a free set whose residual is rounding costs one confirming +/// iteration rather than the iteration cap. A pass that releases a bound is never judged by the cost. /// /// /// Near-singular problems. The normal equations square the condition number, so the cost gap grows with cond(A) @@ -111,7 +113,7 @@ public static class BoundedLeastSquares /// Warm start in, solution out. /// At least floats; contents on entry are ignored. /// Active-set iterations used. - /// Finite and greater than zero, else . Tighten it when rows span many decades of weight; a value below the gradient's rounding noise is harmless because the progress rule ends the solve. + /// Finite and greater than zero, else . Tighten it when rows span many decades of weight; releases below the gradient's rounding bound are not taken, and free-set residuals at rounding level end on the progress rule. /// Iteration cap. public static BoundedLeastSquaresStatus Solve( int m, int n, @@ -200,9 +202,9 @@ public static class BoundedLeastSquares } var kktTol = kktRelativeTolerance * gScale; - // Cost at the last point from which the solver chose to continue past a stationarity miss or a - // release. A pass that ends without lowering it has nothing left to gain: the tolerance is a - // resolution, not the stopping rule, and a noise floor above it must not spin to the cap. + // Cost at the last point from which the solver continued past a stationarity miss with no release; a + // release resets it. A stationarity miss that ends without lowering it has nothing left to gain: the + // tolerance is a resolution, not the stopping rule, and a noise floor above it must not spin to the cap. var fMark = double.PositiveInfinity; var status = BoundedLeastSquaresStatus.IterationLimit; while (true) @@ -264,7 +266,7 @@ public static class BoundedLeastSquares if (lo[j] == hi[j]) continue; // A bound at lo wants a positive gradient and one at hi a negative one; a free column has state 0. var violation = state[j] * g[j]; - if (violation > worst) + if (violation > worst && violation > ReleaseNoise(ata, atb, xd, j, n)) { worst = violation; release = j; @@ -272,9 +274,14 @@ public static class BoundedLeastSquares } if (release < 0 && !freeViolated) { status = BoundedLeastSquaresStatus.Converged; break; } - var f = Objective(xd, g, atb, n, out var fScale); - if (f >= fMark - StallRelativeTolerance * fScale) { status = BoundedLeastSquaresStatus.Converged; break; } - fMark = f; + if (release < 0) + { + var f = Objective(xd, g, atb, n, out var fScale); + if (f >= fMark - StallRelativeTolerance * fScale) { status = BoundedLeastSquaresStatus.Converged; break; } + fMark = f; + } + else + fMark = double.PositiveInfinity; if (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; } iterations++; // A stationarity miss re-solves the same free set from the improved point. @@ -307,6 +314,14 @@ public static class BoundedLeastSquares return f; } + // Rounding bound of the double gradient component j: gamma_(n+1) (Higham, inner products) times the magnitude of its terms. + private static double ReleaseNoise(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, int j, int n) + { + var s = Math.Abs(atb[j]); + for (var k = 0; k < n; k++) s += Math.Abs(ata[j * n + k] * x[k]); + return (n + 1) * 2.220446049250313e-16 * s; + } + private static void Gradient(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, Span g, int n) { for (var j = 0; j < n; j++) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 3523d240..63bc913e 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -1331,6 +1331,67 @@ public sealed class BoundedLeastSquaresTests Assert.True(spread >= 0.5, $"default tolerance: starts differ by only {spread}"); } + // Soul's W5 (finding aetheria-release:finding:cut-thrust-allocator-core.s1.solver-converged-early): an allocator + // matrix (3 force rows, 4 slack pulls, 8 ridge rows; 8 throttles in [0,1], 4 slacks in [0,2]). From x = 0 the + // last release lowered the cost by less than the stall band while throttle 7 still sat at its upper bound + // with a gradient ten thousand times the tolerance, and the solver reported Converged there. + [Theory] + [InlineData(1e-12)] + [InlineData(1e-9)] + public void TheStopAfterAReleaseIsOnTheGradientNotTheCost(double tolerance) + { + const int m = 15, n = 12; + var a = new float[m * n]; + a[0 * n + 2] = 100f; a[0 * n + 8] = -100f; a[0 * n + 9] = 100f; + var row1 = new[] { 13.882518f, 20.86695f, 0f, 14.8965845f, 12.38046f, 22.115599f, -13.706654f, 15.85789f, 0f, 0f, -100f, 100f }; + var row2 = new[] { -187.72859f, 239.31168f, 158.87828f, -270.31165f, 235.38467f, -14.704258f, 121.61232f, 244.81305f }; + for (var j = 0; j < n; j++) a[1 * n + j] = row1[j]; + for (var j = 0; j < row2.Length; j++) a[2 * n + j] = row2[j]; + for (var k = 0; k < 4; k++) a[(3 + k) * n + 8 + k] = 0.3f; + for (var i = 0; i < 8; i++) a[(7 + i) * n + i] = 0.1f; + var b = new float[m]; + b[1] = 67.385445f; b[2] = 558.70917f; + for (var k = 0; k < 4; k++) b[3 + k] = -3f; + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) hi[j] = j < 8 ? 1f : 2f; + + var cold = new float[n]; + Assert.Equal(BoundedLeastSquaresStatus.Converged, SolveAt(m, n, a, b, lo, hi, cold, out _, tolerance)); + Assert.Equal(0.92221445f, cold[7], 1e-4); + + var again = (float[])cold.Clone(); + Assert.Equal(BoundedLeastSquaresStatus.Converged, SolveAt(m, n, a, b, lo, hi, again, out var iterations, tolerance)); + Assert.Equal(0, iterations); + for (var j = 0; j < n; j++) Assert.Equal(cold[j], again[j], 1e-6); + } + + // b = 0 with the origin inside the box: the optimum is any x with Ax = 0, the gradient scale is 0 and so is the tolerance, and + // every gradient component the solve leaves behind is rounding. Releasing on such a violation cycles release / + // stationarity miss / re-bind until the iteration cap (the plain-batch problem at seed 141 in miniature). + [Fact] + public void ANoiseLevelViolationIsNeverReleased() + { + const int m = 2, n = 6; + var a = new[] + { + 1.0213068f, 0.97530603f, -0.5673449f, -0.64136946f, -0.07432885f, -0.35571396f, + 0.4326824f, 1.1969734f, -0.1343108f, -0.23431005f, -1.4009786f, 0.10656489f, + }; + var b = new float[m]; + var lo = new[] { -0.44782114f, 0f, -0.008422324f, -0.03878026f, 0f, -6.6434336f }; + var hi = new[] { 0.0448007f, 0f, 0.005942123f, 0.0785512f, 0f, 0.23079813f }; + var x = new[] { -0.28347927f, 0f, 0.0051532285f, -0.03878026f, 0f, 0.23079813f }; + Assert.Equal(BoundedLeastSquaresStatus.Converged, Solve(m, n, a, b, lo, hi, x, out var iterations)); + Assert.True(iterations <= 4, $"{iterations} iterations"); + for (var i = 0; i < m; i++) + { + var r = 0.0; + for (var j = 0; j < n; j++) r += a[i * n + j] * x[j]; + Assert.Equal(0.0, r, 1e-5); + } + } + [Fact] public void AnInvalidToleranceIsRejectedAndLeavesXUntouched() {