GameCult/CultLib
Watch
9
0
Fork
You've already forked CultLib
 0
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions

CultMath BoundedLeastSquares: release on the gradient against its rounding, never by the cost

Browse source 
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 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-10 08:51:18 +02:00
commit dbcb2a8ea8
2 changed files with 89 additions and 13 deletions
Download patch file Download diff file Expand all files Collapse all files

41
packages/cultmath/src/CultMath/BoundedLeastSquares.cs
View file

@ -34,11 +34,13 @@ public enum BoundedLeastSquaresStatus
/// Stopping. KKT termination uses a tolerance of <see cref="DefaultKktRelativeTolerance"/> (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 <c>kktRelativeTolerance</c>, 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
/// <see cref="BoundedLeastSquaresStatus.Converged"/>, so an optimal point whose double gradient sits above the
/// tolerance costs one confirming iteration rather than the iteration cap.
/// a smaller <c>kktRelativeTolerance</c>, 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
/// <see cref="BoundedLeastSquaresStatus.Converged"/>, 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.
/// </para>
/// <para>
/// 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
/// <param name="x">Warm start in, solution out.</param>
/// <param name="workspace">At least <see cref="WorkspaceLength"/> floats; contents on entry are ignored.</param>
/// <param name="iterations">Active-set iterations used.</param>
/// <param name="kktRelativeTolerance">Finite and greater than zero, else <see cref="BoundedLeastSquaresStatus.InvalidInput"/>. 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.</param>
/// <param name="kktRelativeTolerance">Finite and greater than zero, else <see cref="BoundedLeastSquaresStatus.InvalidInput"/>. 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.</param>
/// <param name="maxIterations">Iteration cap.</param>
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<double> ata, ReadOnlySpan<double> atb, ReadOnlySpan<double> 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<double> ata, ReadOnlySpan<double> atb, ReadOnlySpan<double> x, Span<double> g, int n)
{
for (var j = 0; j < n; j++)

61
packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs
View file

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