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()
{