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 iterates in double and stops on progress

Browse source 
The working iterate now lives in double inside the workspace (WorkspaceLength
4n^2 + 17n); A, b and the API x stay float and x is converted in at the start
and out at the end. The gradient, the free-stationarity check and the KKT test
all run in the double domain, which removes the float noise floor that made the
free-stationarity re-solve spin to the cap on optimal points (cond >= ~20).

A fixed tolerance is no longer the only stopping rule: after a stationarity
miss or a release, a pass that lowers the double cost by no more than rounding
reports Converged.

The tolerance scale is the larger of the gradient at the origin and at the box
point nearest the origin, so it never depends on the warm start; b = 0 with a
far start is no longer judged converged by its own steepness.

Two clauses were equivalent or unobservable and are deleted: the d <= 0 pivot
test (implied by d <= tol * diag) and the dependent[] exclusion from the
free-stationarity check (the progress rule now terminates it).

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-09-29 05:36:50 +02:00
commit 43035bc827
1 changed files with 82 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -23,11 +23,29 @@ public enum BoundedLeastSquaresStatus
/// Primal active set. Each iteration solves the free variables by Cholesky on the normal equations, takes the
/// step to the first bound it would cross, and, once a full step lands, releases the bound with the most
/// violated gradient sign. Free columns that are linearly dependent on earlier free columns (a zero column, a
/// duplicate) get a zero step instead of a failed pivot, so a singular free set never throws. The normal
/// equations square the condition number, so they are formed and factored in double (inputs and <c>x</c> stay
/// float): callers that need a minimum-norm answer on rank-deficient problems add a small ridge row to
/// <c>A</c> and <c>b</c>; each free-set solve gets one step of iterative refinement. The KKT tolerance is relative to the problem's gradient scale (max |A^T b|), so
/// scaling A, b or the bounds does not change the answer. The solver allocates nothing and is deterministic.
/// duplicate) get a zero step instead of a failed pivot, so a singular free set never throws.
/// <para>
/// Precision. A, b and <c>x</c> are float at the API; everything else, including the working iterate, the
/// gradient and the KKT test, is double, carved from the caller's workspace. <c>x</c> is converted in at the start
/// and back out at the end, so a returned float <c>x</c> is the nearest float to the double solution and always
/// inside the box. Each free-set solve gets one step of iterative refinement.
/// </para>
/// <para>
/// Stopping. KKT termination uses a tolerance of 1e-6 times the problem's gradient scale (the larger of the
/// gradient magnitudes at the origin and at the box point nearest the origin), so scaling A, b or the bounds, and
/// the choice of warm start, do not change it. 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.
/// </para>
/// <para>
/// Near-singular problems. Columns whose pivot falls below 1e-12 of their own diagonal (cond(A) beyond roughly 1e6)
/// are treated as dependent and not moved, so the answer there is one optimum among many, not the minimum-norm one.
/// Callers that need the minimum-norm answer, or that expect cond(A) beyond that cutoff, add a small ridge row to
/// <c>A</c> and <c>b</c>; that stays the caller's policy. Below the cutoff the solver owns optimality: near-collinear
/// problems reach the exact optimum to about 1e-7 of |b|^2 at any absolute scale of A.
/// </para>
/// The solver allocates nothing and is deterministic.
/// <para>
/// Input contract: <c>a</c> and <c>b</c> must be finite; <c>lo</c> and <c>hi</c> must not be NaN and need
/// <c>lo &lt;= hi</c>. Infinite bounds are allowed and mean unbounded on that side (<c>lo = -Inf</c>,
@ -45,9 +63,10 @@ public static class BoundedLeastSquares
private const double PivotRelativeTolerance = 1e-12;
private const double KktRelativeTolerance = 1e-6;
private const double StallRelativeTolerance = 1e-15;
/// <summary>Floats of caller-supplied workspace needed for <paramref name="n"/> columns.</summary>
public static int WorkspaceLength(int n) => 2 * (2 * n * n + 6 * n) + 3 * n;
public static int WorkspaceLength(int n) => 2 * (2 * n * n + 7 * n) + 3 * n;
/// <summary>
/// Minimises <c>||A x - b||^2</c> subject to <c>lo &lt;= x &lt;= hi</c>. <paramref name="x"/> is the warm start
@ -76,7 +95,7 @@ public static class BoundedLeastSquares
throw new ArgumentException("BoundedLeastSquares: span shorter than the stated dimensions.");
// Double-precision vectors and matrices first, then the integer flags, all carved from the caller's floats.
var doubleCount = 2 * n * n + 6 * n;
var doubleCount = 2 * n * n + 7 * n;
var dbl = MemoryMarshal.Cast<float, double>(workspace.Slice(0, 2 * doubleCount));
var ata = dbl.Slice(0, n * n);
var chol = dbl.Slice(n * n, n * n);
@ -87,6 +106,7 @@ public static class BoundedLeastSquares
var y = dbl.Slice(off + 3 * n, n);
var r1 = dbl.Slice(off + 4 * n, n);
var q = dbl.Slice(off + 5 * n, n);
var xd = dbl.Slice(off + 6 * n, n);
// One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow.
var ints = MemoryMarshal.Cast<float, int>(workspace.Slice(2 * doubleCount, 3 * n));
var state = ints.Slice(0, n);
@ -104,14 +124,12 @@ public static class BoundedLeastSquares
if (!float.IsFinite(b[i]))
return BoundedLeastSquaresStatus.InvalidInput;
var atbMax = 0.0;
for (var j = 0; j < n; j++)
{
var sb = 0.0;
for (var r = 0; r < m; r++)
sb += (double)a[r * n + j] * b[r];
atb[j] = sb;
atbMax = Math.Max(atbMax, Math.Abs(sb));
for (var k = 0; k <= j; k++)
{
var sa = 0.0;
@ -126,6 +144,7 @@ public static class BoundedLeastSquares
// From here x is written. A warm start that is NaN or an infinite value on an unbounded side is
// moved onto the box; an infinite result (unbounded side) restarts from 0 clamped to the box.
// The working iterate lives in double; x is only the API's float view of it.
for (var j = 0; j < n; j++)
{
var l = lo[j];
@ -135,23 +154,34 @@ public static class BoundedLeastSquares
if (float.IsInfinity(v))
v = 0f < l ? l : (0f > h ? h : 0f);
x[j] = v;
xd[j] = v;
state[j] = v == l ? -1 : (v == h ? 1 : 0);
}
// Scale-free tolerance: relative to the gradient scale of the problem, with no absolute floor.
// When A^T b is exactly zero the start's gradient supplies the scale (zero only if x is already optimal).
var gScale = atbMax;
if (gScale == 0.0)
// The tolerance's scale is the gradient magnitude the problem itself carries: the larger of the
// gradient at the origin (-A^T b) and at the box point nearest the origin. It depends on A, b and the
// bounds, never on the warm start, so a far start cannot be judged converged by its own steepness, and
// b = 0 with a box that excludes the origin still gets a real scale. It is zero only when that point
// is already stationary, and then the progress rule below does the stopping.
var gScale = 0.0;
for (var j = 0; j < n; j++)
{
Gradient(ata, atb, x, g, n);
for (var j = 0; j < n; j++)
gScale = Math.Max(gScale, Math.Abs(g[j]));
q[j] = Math.Clamp(0.0, lo[j], hi[j]);
gScale = Math.Max(gScale, Math.Abs(atb[j]));
}
Gradient(ata, atb, q, g, n);
for (var j = 0; j < n; j++)
gScale = Math.Max(gScale, Math.Abs(g[j]));
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.
var fMark = double.PositiveInfinity;
BoundedLeastSquaresStatus status;
while (true)
{
Gradient(ata, atb, x, g, n);
Gradient(ata, atb, xd, g, n);
var freeCount = 0;
for (var j = 0; j < n; j++)
@ -166,7 +196,7 @@ public static class BoundedLeastSquares
var j = freeIndex[i];
var pj = p[j];
if (pj == 0.0) continue;
var tj = pj > 0f ? (hi[j] - x[j]) / pj : (lo[j] - x[j]) / pj;
var tj = pj > 0.0 ? (hi[j] - xd[j]) / pj : (lo[j] - xd[j]) / pj;
if (tj < t)
{
t = tj;
@ -176,30 +206,29 @@ public static class BoundedLeastSquares
if (block >= 0)
{
if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit;
if (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; }
iterations++;
for (var i = 0; i < freeCount; i++)
{
var j = freeIndex[i];
x[j] = (float)(x[j] + t * p[j]);
xd[j] = Math.Clamp(xd[j] + t * p[j], lo[j], hi[j]);
}
if (p[block] > 0.0) { x[block] = hi[block]; state[block] = 1; }
else { x[block] = lo[block]; state[block] = -1; }
if (p[block] > 0.0) { xd[block] = hi[block]; state[block] = 1; }
else { xd[block] = lo[block]; state[block] = -1; }
continue;
}
for (var i = 0; i < freeCount; i++)
{
var j = freeIndex[i];
var v = (float)(x[j] + p[j]);
x[j] = v < lo[j] ? lo[j] : (v > hi[j] ? hi[j] : v);
xd[j] = Math.Clamp(xd[j] + p[j], lo[j], hi[j]);
}
Gradient(ata, atb, x, g, n);
Gradient(ata, atb, xd, g, n);
// Stationarity of the free set is part of KKT: float Cholesky can leave a residual gradient.
// Stationarity of the free set is part of KKT: a Newton step on a near-singular free set can leave a residual gradient.
var freeViolated = false;
for (var i = 0; i < freeCount; i++)
if (dependent[i] == 0 && Math.Abs(g[freeIndex[i]]) > kktTol)
if (Math.Abs(g[freeIndex[i]]) > kktTol)
freeViolated = true;
var release = -1;
@ -214,15 +243,37 @@ public static class BoundedLeastSquares
release = j;
}
}
if (release < 0 && !freeViolated) return BoundedLeastSquaresStatus.Converged;
if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit;
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 (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; }
iterations++;
// A stationarity miss re-solves the same free set from the improved point.
if (release >= 0) state[release] = 0;
}
for (var j = 0; j < n; j++)
x[j] = (float)xd[j];
return status;
}
private static void Gradient(ReadOnlySpan<double> ata, ReadOnlySpan<double> atb, ReadOnlySpan<float> x, Span<double> g, int n)
// 0.5 x^T A^T A x - x^T A^T b (the cost up to a constant), from the gradient at x; scale is the magnitude of
// the terms summed, which is what the rounding error of the sum is proportional to.
private static double Objective(ReadOnlySpan<double> x, ReadOnlySpan<double> g, ReadOnlySpan<double> atb, int n, out double scale)
{
var f = 0.0;
scale = 0.0;
for (var j = 0; j < n; j++)
{
f += 0.5 * x[j] * (g[j] - atb[j]);
scale += 0.5 * Math.Abs(x[j]) * (Math.Abs(g[j]) + Math.Abs(atb[j]));
}
return f;
}
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++)
{
@ -247,7 +298,7 @@ public static class BoundedLeastSquares
var d = ata[fi * n + fi];
for (var c = 0; c < i; c++)
d -= l[i * n + c] * l[i * n + c];
if (d <= PivotRelativeTolerance * ata[fi * n + fi] || d <= 0.0)
if (d <= PivotRelativeTolerance * ata[fi * n + fi])
{
dependent[i] = 1;
for (var c = 0; c < i; c++) l[i * n + c] = 0.0;