From 60c3d6fb0d46509e32275ef2b0bcb0d828bfea93 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 01:33:08 +0200 Subject: [PATCH 01/19] cultmath: BoundedLeastSquares primal active-set solver Dense small bounded LS, zero allocation, caller workspace, C#-only. Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 233 +++++++++++++ .../src/CultMath/BoundedLeastSquares.cs.meta | 2 + .../BoundedLeastSquaresTests.cs | 330 ++++++++++++++++++ 3 files changed, 565 insertions(+) create mode 100644 packages/cultmath/src/CultMath/BoundedLeastSquares.cs create mode 100644 packages/cultmath/src/CultMath/BoundedLeastSquares.cs.meta create mode 100644 packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs new file mode 100644 index 00000000..e18a42e2 --- /dev/null +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -0,0 +1,233 @@ +using System.Runtime.InteropServices; + +namespace CultMath; + +/// Outcome of . +public enum BoundedLeastSquaresStatus +{ + /// The KKT conditions hold: every free variable has zero gradient and every bound has a gradient pointing into it. + Converged, + + /// The iteration cap was reached. x holds the last feasible iterate, which is never worse than the clamped start. + IterationLimit, +} + +/// +/// Dense bounded least squares: min ||A x - b||^2 subject to lo <= x <= hi, for small problems +/// (tens of rows and columns). C#-only, like the intercept helpers in ; there is no HLSL mirror. +/// +/// +/// 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 and run in float: callers that need a minimum-norm answer on +/// rank-deficient problems add a small ridge row to A and b. The solver allocates nothing and is +/// deterministic. +/// +public static class BoundedLeastSquares +{ + /// Default iteration cap for . + public const int DefaultMaxIterations = 100; + + private const float PivotRelativeTolerance = 1e-6f; + private const float KktRelativeTolerance = 1e-5f; + + /// Floats of caller-supplied workspace needed for columns. + public static int WorkspaceLength(int n) => 2 * n * n + 7 * n; + + /// + /// Minimises ||A x - b||^2 subject to lo <= x <= hi. is the warm start + /// on entry (clamped to the box first) and the result on exit. + /// + /// Rows of . + /// Columns of , and length of , and . + /// Row-major m * n matrix. + /// Length m target. + /// Lower bounds; each must be <= hi. + /// Upper bounds. + /// Warm start in, solution out. + /// At least floats; contents on entry are ignored. + /// Active-set iterations used. + /// Iteration cap. + public static BoundedLeastSquaresStatus Solve( + int m, int n, + ReadOnlySpan a, ReadOnlySpan b, + ReadOnlySpan lo, ReadOnlySpan hi, + Span x, Span workspace, + out int iterations, + int maxIterations = DefaultMaxIterations) + { + if (m < 0 || n < 0 || a.Length < m * n || b.Length < m || lo.Length < n || hi.Length < n + || x.Length < n || workspace.Length < WorkspaceLength(n)) + throw new ArgumentException("BoundedLeastSquares: span shorter than the stated dimensions."); + + var ata = workspace.Slice(0, n * n); + var chol = workspace.Slice(n * n, n * n); + var off = 2 * n * n; + var atb = workspace.Slice(off, n); + var g = workspace.Slice(off + n, n); + var p = workspace.Slice(off + 2 * n, n); + var y = workspace.Slice(off + 3 * n, n); + // One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow. + var ints = MemoryMarshal.Cast(workspace.Slice(off + 4 * n, 3 * n)); + var state = ints.Slice(0, n); + var freeIndex = ints.Slice(n, n); + var dependent = ints.Slice(2 * n, n); + + for (var j = 0; j < n; j++) + { + var l = lo[j]; + var h = hi[j]; + var v = x[j]; + v = !(v > l) ? l : (v < h ? v : h); + x[j] = v; + state[j] = v == l ? -1 : (v == h ? 1 : 0); + } + + var atbMax = 0f; + for (var j = 0; j < n; j++) + { + var sb = 0f; + for (var r = 0; r < m; r++) + sb += a[r * n + j] * b[r]; + atb[j] = sb; + atbMax = MathF.Max(atbMax, MathF.Abs(sb)); + for (var k = 0; k <= j; k++) + { + var s = 0f; + for (var r = 0; r < m; r++) + s += a[r * n + j] * a[r * n + k]; + ata[j * n + k] = s; + ata[k * n + j] = s; + } + } + var kktTol = KktRelativeTolerance * (1f + atbMax); + + iterations = 0; + while (true) + { + Gradient(ata, atb, x, g, n); + + var freeCount = 0; + for (var j = 0; j < n; j++) + if (state[j] == 0) + freeIndex[freeCount++] = j; + SolveFreeSet(ata, chol, g, p, y, freeIndex.Slice(0, freeCount), dependent, n); + + var t = 1f; + var block = -1; + for (var i = 0; i < freeCount; i++) + { + var j = freeIndex[i]; + var pj = p[j]; + if (pj == 0f) continue; + var tj = pj > 0f ? (hi[j] - x[j]) / pj : (lo[j] - x[j]) / pj; + if (tj < t) + { + t = tj; + block = j; + } + } + + if (block >= 0) + { + if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit; + iterations++; + for (var i = 0; i < freeCount; i++) + { + var j = freeIndex[i]; + x[j] += t * p[j]; + } + if (p[block] > 0f) { x[block] = hi[block]; state[block] = 1; } + else { x[block] = lo[block]; state[block] = -1; } + continue; + } + + for (var i = 0; i < freeCount; i++) + { + var j = freeIndex[i]; + var v = x[j] + p[j]; + x[j] = v < lo[j] ? lo[j] : (v > hi[j] ? hi[j] : v); + } + Gradient(ata, atb, x, g, n); + + var release = -1; + var worst = kktTol; + for (var j = 0; j < n; j++) + { + if (lo[j] == hi[j]) continue; + var violation = state[j] < 0 ? -g[j] : (state[j] > 0 ? g[j] : 0f); + if (violation > worst) + { + worst = violation; + release = j; + } + } + if (release < 0) return BoundedLeastSquaresStatus.Converged; + if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit; + iterations++; + state[release] = 0; + } + } + + private static void Gradient(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, Span g, int n) + { + for (var j = 0; j < n; j++) + { + var s = -atb[j]; + for (var k = 0; k < n; k++) + s += ata[j * n + k] * x[k]; + g[j] = s; + } + } + + // Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky. A free column whose pivot + // collapses (dependent on earlier free columns) is dropped from the factor and gets p = 0. + private static void SolveFreeSet( + ReadOnlySpan ata, Span l, ReadOnlySpan g, Span p, Span y, + ReadOnlySpan free, Span dependent, int n) + { + var k = free.Length; + for (var i = 0; i < k; i++) + { + var fi = free[i]; + 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 <= 0f) + { + dependent[i] = 1; + for (var c = 0; c < i; c++) l[i * n + c] = 0f; + l[i * n + i] = 1f; + for (var r = i + 1; r < k; r++) l[r * n + i] = 0f; + continue; + } + dependent[i] = 0; + var lii = MathF.Sqrt(d); + l[i * n + i] = lii; + for (var r = i + 1; r < k; r++) + { + var s = ata[free[r] * n + fi]; + for (var c = 0; c < i; c++) + s -= l[r * n + c] * l[i * n + c]; + l[r * n + i] = s / lii; + } + } + + for (var i = 0; i < k; i++) + { + var s = dependent[i] != 0 ? 0f : -g[free[i]]; + for (var c = 0; c < i; c++) + s -= l[i * n + c] * y[c]; + y[i] = s / l[i * n + i]; + } + for (var i = k - 1; i >= 0; i--) + { + var s = y[i]; + for (var r = i + 1; r < k; r++) + s -= l[r * n + i] * p[free[r]]; + p[free[i]] = s / l[i * n + i]; + } + } +} diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs.meta b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs.meta new file mode 100644 index 00000000..fb5bd819 --- /dev/null +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 4ba991e1b6c54e099bc7ea16541734c5 diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs new file mode 100644 index 00000000..653b9158 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -0,0 +1,330 @@ +using System; +using CultMath; +using Xunit; + +namespace CultMath.Tests; + +// Fixtures are deliberately non-degenerate: A is non-square and non-identity, and b is outside A's range, +// so the unconstrained optimum has a real residual and the box actually bites. +public sealed class BoundedLeastSquaresTests +{ + private static double Cost(int m, int n, float[] a, float[] b, float[] x) + { + var c = 0.0; + for (var r = 0; r < m; r++) + { + var s = -(double)b[r]; + for (var j = 0; j < n; j++) s += (double)a[r * n + j] * x[j]; + c += s * s; + } + return c; + } + + // Gradient of 0.5*||Ax-b||^2, in double. + private static double[] Gradient(int m, int n, float[] a, float[] b, float[] x) + { + var g = new double[n]; + for (var r = 0; r < m; r++) + { + var s = -(double)b[r]; + for (var j = 0; j < n; j++) s += (double)a[r * n + j] * x[j]; + for (var j = 0; j < n; j++) g[j] += a[r * n + j] * s; + } + return g; + } + + private static double AtbScale(int m, int n, float[] a, float[] b) + { + var mx = 0.0; + for (var j = 0; j < n; j++) + { + var s = 0.0; + for (var r = 0; r < m; r++) s += (double)a[r * n + j] * b[r]; + mx = Math.Max(mx, Math.Abs(s)); + } + return 1 + mx; + } + + // Asserts feasibility and the KKT sign rule: free columns have ~zero gradient, columns at lo have + // gradient >= 0, columns at hi have gradient <= 0. + private static void AssertKkt(int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[] x) + { + var g = Gradient(m, n, a, b, x); + var tol = 1e-3 * AtbScale(m, n, a, b); + for (var j = 0; j < n; j++) + { + Assert.True(x[j] >= lo[j] && x[j] <= hi[j], $"column {j} outside its box: {x[j]}"); + if (lo[j] == hi[j]) continue; + if (x[j] == lo[j]) Assert.True(g[j] >= -tol, $"column {j} at lo with gradient {g[j]}"); + else if (x[j] == hi[j]) Assert.True(g[j] <= tol, $"column {j} at hi with gradient {g[j]}"); + else Assert.True(Math.Abs(g[j]) <= tol, $"column {j} free with gradient {g[j]}"); + } + } + + // Projected gradient with step 1/L, L bounded by the Frobenius norm squared of A. + private static double ReferenceCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi) + { + var frob = 0.0; + foreach (var v in a) frob += (double)v * v; + var step = 1.0 / Math.Max(frob, 1e-12); + var x = new float[n]; + for (var j = 0; j < n; j++) x[j] = Math.Clamp(0f, lo[j], hi[j]); + var xd = new double[n]; + for (var j = 0; j < n; j++) xd[j] = x[j]; + for (var it = 0; it < 200000; it++) + { + for (var j = 0; j < n; j++) x[j] = (float)xd[j]; + var g = Gradient(m, n, a, b, x); + for (var j = 0; j < n; j++) xd[j] = Math.Clamp(xd[j] - step * g[j], lo[j], hi[j]); + } + for (var j = 0; j < n; j++) x[j] = (float)xd[j]; + return Cost(m, n, a, b, x); + } + + private static BoundedLeastSquaresStatus Solve(int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[] x, + out int iterations, int maxIterations = BoundedLeastSquares.DefaultMaxIterations) + => BoundedLeastSquares.Solve(m, n, a, b, lo, hi, x, new float[BoundedLeastSquares.WorkspaceLength(n)], out iterations, maxIterations); + + private static float[] Fill(float value, int n) + { + var v = new float[n]; + Array.Fill(v, value); + return v; + } + + private static readonly float[] A54 = + { + 1.0f, 0.3f, -0.2f, + 0.2f, 1.1f, 0.4f, + -0.5f, 0.6f, 0.9f, + 0.7f, -0.4f, 0.5f, + 0.1f, 0.8f, -0.6f, + }; + + private static readonly float[] B5 = { 0.9f, -0.4f, 1.3f, 0.2f, -0.7f }; + + [Fact] + public void InteriorOptimumEqualsNormalEquations() + { + const int m = 5, n = 3; + var x = new float[n]; + var status = Solve(m, n, A54, B5, Fill(-100f, n), Fill(100f, n), x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + + // Independent solve of (A^T A) x = A^T b by Gauss-Jordan in double. + var mat = new double[n, n + 1]; + for (var i = 0; i < n; i++) + { + for (var j = 0; j < n; j++) + for (var r = 0; r < m; r++) mat[i, j] += (double)A54[r * n + i] * A54[r * n + j]; + for (var r = 0; r < m; r++) mat[i, n] += (double)A54[r * n + i] * B5[r]; + } + for (var c = 0; c < n; c++) + { + var piv = mat[c, c]; + for (var j = 0; j <= n; j++) mat[c, j] /= piv; + for (var i = 0; i < n; i++) + { + if (i == c) continue; + var f = mat[i, c]; + for (var j = 0; j <= n; j++) mat[i, j] -= f * mat[c, j]; + } + } + for (var i = 0; i < n; i++) Assert.Equal(mat[i, n], x[i], 1e-4); + } + + [Fact] + public void KktHoldsOnRandomBoxedProblems() + { + var anyActive = false; + for (uint seed = 1; seed <= 150; seed++) + { + var rng = new CultMath.Random(seed * 7919u); + var n = rng.NextInt(2, 17); + var m = rng.NextInt(n, 25); + var a = new float[m * n]; + var b = new float[m]; + for (var i = 0; i < a.Length; i++) a[i] = rng.NextFloat(-1f, 1f); + for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) + { + switch (rng.NextInt(3)) + { + case 0: lo[j] = -1f; hi[j] = 1f; break; + case 1: lo[j] = 0f; hi[j] = 1f; break; + default: lo[j] = -0.25f; hi[j] = 0.1f; break; + } + } + var x = new float[n]; + var status = Solve(m, n, a, b, lo, hi, x, out _); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed} ({m}x{n}) hit the cap"); + AssertKkt(m, n, a, b, lo, hi, x); + for (var j = 0; j < n; j++) anyActive |= x[j] == lo[j] || x[j] == hi[j]; + + var cs = Cost(m, n, a, b, x); + var cr = ReferenceCost(m, n, a, b, lo, hi); + Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6, $"seed {seed} ({m}x{n}): solver cost {cs}, reference {cr}"); + } + Assert.True(anyActive); + } + + [Fact] + public void AllBoundsActive() + { + const int m = 4, n = 3; + var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f }; + var b = new[] { 10f, -10f, 10f, 1f }; + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + var x = new float[n]; + var status = Solve(m, n, a, b, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(new[] { 1f, -1f, 1f }, x); + AssertKkt(m, n, a, b, lo, hi, x); + } + + [Fact] + public void ZeroColumn() + { + const int m = 5, n = 3; + var a = (float[])A54.Clone(); + for (var r = 0; r < m; r++) a[r * n + 1] = 0f; + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + var x = new float[n]; + var status = Solve(m, n, a, B5, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + AssertKkt(m, n, a, B5, lo, hi, x); + var cs = Cost(m, n, a, B5, x); + var cr = ReferenceCost(m, n, a, B5, lo, hi); + Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6); + } + + [Fact] + public void DuplicateColumns() + { + const int m = 5, n = 3; + var a = (float[])A54.Clone(); + for (var r = 0; r < m; r++) a[r * n + 2] = a[r * n + 0]; + var lo = Fill(-1f, n); + var hi = new[] { 0.3f, 1f, 0.3f }; + var x = new float[n]; + var status = Solve(m, n, a, B5, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + AssertKkt(m, n, a, B5, lo, hi, x); + var cs = Cost(m, n, a, B5, x); + var cr = ReferenceCost(m, n, a, B5, lo, hi); + Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6); + } + + [Fact] + public void WarmStartMatchesColdStart() + { + const int m = 5, n = 3; + var lo = new[] { -0.2f, -1f, 0f }; + var hi = new[] { 0.5f, 0.3f, 1f }; + var cold = new float[n]; + Solve(m, n, A54, B5, lo, hi, cold, out _); + AssertKkt(m, n, A54, B5, lo, hi, cold); + + var starts = new[] { new[] { 0.5f, -1f, 1f }, new[] { 9f, 9f, -9f }, new[] { 0.1f, 0.1f, 0.5f }, (float[])cold.Clone() }; + foreach (var start in starts) + { + var status = Solve(m, n, A54, B5, lo, hi, start, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + for (var j = 0; j < n; j++) Assert.Equal(cold[j], start[j], 1e-4); + if (ReferenceEquals(start, starts[^1])) Assert.Equal(0, iterations); + } + } + + [Fact] + public void AsymmetricBounds() + { + const int m = 5, n = 3; + var lo = new[] { -1f, 0f, -1f }; + var hi = new[] { 1f, 1f, 1f }; + // Column 1 wants to go negative (b has a strong component against it) but its floor is 0. + var b = new[] { 0.2f, -2.5f, 0.4f, -1.5f, -2.0f }; + var x = new float[n]; + var status = Solve(m, n, A54, b, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + AssertKkt(m, n, A54, b, lo, hi, x); + Assert.Equal(0f, x[1]); + Assert.True(x[0] < 0f || x[2] < 0f, "a symmetric-bounds column should use its negative side"); + var cs = Cost(m, n, A54, b, x); + var cr = ReferenceCost(m, n, A54, b, lo, hi); + Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6); + } + + [Fact] + public void IterationCapReturnsFeasibleIterateAndLimitStatus() + { + const int m = 4, n = 3; + var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f }; + var b = new[] { 10f, -10f, 10f, 1f }; + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + var x = new float[n]; + var start = Cost(m, n, a, b, x); + var status = Solve(m, n, a, b, lo, hi, x, out var iterations, maxIterations: 1); + Assert.Equal(BoundedLeastSquaresStatus.IterationLimit, status); + Assert.Equal(1, iterations); + for (var j = 0; j < n; j++) Assert.InRange(x[j], lo[j], hi[j]); + Assert.True(Cost(m, n, a, b, x) < start); + + var zero = new float[n]; + Assert.Equal(BoundedLeastSquaresStatus.IterationLimit, Solve(m, n, a, b, lo, hi, zero, out iterations, maxIterations: 0)); + Assert.Equal(0, iterations); + Assert.Equal(new float[n], zero); + } + + [Fact] + public void StartOutsideTheBoxIsClampedFirst() + { + const int m = 5, n = 3; + var lo = Fill(-0.1f, n); + var hi = Fill(0.1f, n); + var x = new[] { 50f, -50f, float.NaN }; + Solve(m, n, A54, B5, lo, hi, x, out _, maxIterations: 0); + Assert.Equal(new[] { 0.1f, -0.1f, -0.1f }, x); + } + + [Fact] + public void FixedColumnStaysPut() + { + const int m = 5, n = 3; + var lo = new[] { -1f, 0.25f, -1f }; + var hi = new[] { 1f, 0.25f, 1f }; + var x = new float[n]; + var status = Solve(m, n, A54, B5, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0.25f, x[1]); + AssertKkt(m, n, A54, B5, lo, hi, x); + } + + [Fact] + public void WarmCallAllocatesNothing() + { + const int m = 5, n = 3; + var lo = new[] { -0.2f, -1f, 0f }; + var hi = new[] { 0.5f, 0.3f, 1f }; + var x = new float[n]; + var workspace = new float[BoundedLeastSquares.WorkspaceLength(n)]; + BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, x, workspace, out _); + + var before = GC.GetAllocatedBytesForCurrentThread(); + var status = BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, x, workspace, out _); + var after = GC.GetAllocatedBytesForCurrentThread(); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0L, after - before); + } + + [Fact] + public void ShortSpansThrowArgumentException() + { + var x = new float[3]; + Assert.Throws(() => BoundedLeastSquares.Solve(5, 3, A54, B5, Fill(-1f, 3), Fill(1f, 3), x, new float[4], out _)); + } +} From eec1e40c258156827cae846be66f5c6a82d30444 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 01:38:56 +0200 Subject: [PATCH 02/19] cultmath: pin that a fixed column is never released Co-Authored-By: Claude Sonnet 5 --- .../tests/CultMath.Tests/BoundedLeastSquaresTests.cs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 653b9158..6cf65fd4 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -292,15 +292,17 @@ public sealed class BoundedLeastSquaresTests } [Fact] - public void FixedColumnStaysPut() + public void FixedColumnStaysPutEvenWhenItsGradientWantsOut() { const int m = 5, n = 3; - var lo = new[] { -1f, 0.25f, -1f }; - var hi = new[] { 1f, 0.25f, 1f }; + var lo = new[] { -1f, -1f, -1f }; + var hi = new[] { 1f, -1f, 1f }; var x = new float[n]; - var status = Solve(m, n, A54, B5, lo, hi, x, out _); + var status = Solve(m, n, A54, B5, lo, hi, x, out var iterations); Assert.Equal(BoundedLeastSquaresStatus.Converged, status); - Assert.Equal(0.25f, x[1]); + Assert.Equal(-1f, x[1]); + Assert.True(Gradient(m, n, A54, B5, x)[1] < -1e-3, "fixture must pull the fixed column inward"); + Assert.Equal(0, iterations); AssertKkt(m, n, A54, B5, lo, hi, x); } From 5382b81036c0a6fd74b958d8462ea08959c2f4b3 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 02:14:39 +0200 Subject: [PATCH 03/19] cultmath: widen fixed-column fixture so no other column blocks Co-Authored-By: Claude Sonnet 5 --- .../cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 6cf65fd4..0d99b341 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -295,8 +295,8 @@ public sealed class BoundedLeastSquaresTests public void FixedColumnStaysPutEvenWhenItsGradientWantsOut() { const int m = 5, n = 3; - var lo = new[] { -1f, -1f, -1f }; - var hi = new[] { 1f, -1f, 1f }; + var lo = new[] { -100f, -1f, -100f }; + var hi = new[] { 100f, -1f, 100f }; var x = new float[n]; var status = Solve(m, n, A54, B5, lo, hi, x, out var iterations); Assert.Equal(BoundedLeastSquaresStatus.Converged, status); From c65777b40013b0cf1c015373b500dd07d2f65a61 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 02:33:51 +0200 Subject: [PATCH 04/19] cultmath: stop the projected-gradient reference once it stops moving Co-Authored-By: Claude Sonnet 5 --- .../tests/CultMath.Tests/BoundedLeastSquaresTests.cs | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 0d99b341..e14cb50c 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -75,7 +75,14 @@ public sealed class BoundedLeastSquaresTests { for (var j = 0; j < n; j++) x[j] = (float)xd[j]; var g = Gradient(m, n, a, b, x); - for (var j = 0; j < n; j++) xd[j] = Math.Clamp(xd[j] - step * g[j], lo[j], hi[j]); + var moved = 0.0; + for (var j = 0; j < n; j++) + { + var next = Math.Clamp(xd[j] - step * g[j], lo[j], hi[j]); + moved = Math.Max(moved, Math.Abs(next - xd[j])); + xd[j] = next; + } + if (moved < 1e-10) break; } for (var j = 0; j < n; j++) x[j] = (float)xd[j]; return Cost(m, n, a, b, x); From 32fb7c15034e73080fb881e82e2fdc31887a9c6e Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 02:52:51 +0200 Subject: [PATCH 05/19] cultmath: FISTA reference keeps the random KKT test fast Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 31 +++++++++++-------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index e14cb50c..2b3e870d 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -61,31 +61,36 @@ public sealed class BoundedLeastSquaresTests } } - // Projected gradient with step 1/L, L bounded by the Frobenius norm squared of A. + // Accelerated projected gradient (FISTA) with step 1/L, L bounded by the Frobenius norm squared of A. + // Independent of the solver's active-set logic; run to a fixed point. private static double ReferenceCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi) { var frob = 0.0; foreach (var v in a) frob += (double)v * v; var step = 1.0 / Math.Max(frob, 1e-12); - var x = new float[n]; - for (var j = 0; j < n; j++) x[j] = Math.Clamp(0f, lo[j], hi[j]); - var xd = new double[n]; - for (var j = 0; j < n; j++) xd[j] = x[j]; - for (var it = 0; it < 200000; it++) + var x = new double[n]; + var y = new double[n]; + var xf = new float[n]; + for (var j = 0; j < n; j++) x[j] = y[j] = Math.Clamp(0.0, lo[j], hi[j]); + var t = 1.0; + for (var it = 0; it < 100000; it++) { - for (var j = 0; j < n; j++) x[j] = (float)xd[j]; - var g = Gradient(m, n, a, b, x); + for (var j = 0; j < n; j++) xf[j] = (float)y[j]; + var g = Gradient(m, n, a, b, xf); var moved = 0.0; + var tNext = (1.0 + Math.Sqrt(1.0 + 4.0 * t * t)) / 2.0; for (var j = 0; j < n; j++) { - var next = Math.Clamp(xd[j] - step * g[j], lo[j], hi[j]); - moved = Math.Max(moved, Math.Abs(next - xd[j])); - xd[j] = next; + var next = Math.Clamp(y[j] - step * g[j], lo[j], hi[j]); + moved = Math.Max(moved, Math.Abs(next - x[j])); + y[j] = next + (t - 1.0) / tNext * (next - x[j]); + x[j] = next; } + t = tNext; if (moved < 1e-10) break; } - for (var j = 0; j < n; j++) x[j] = (float)xd[j]; - return Cost(m, n, a, b, x); + for (var j = 0; j < n; j++) xf[j] = (float)x[j]; + return Cost(m, n, a, b, xf); } private static BoundedLeastSquaresStatus Solve(int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[] x, From 986711c9d597db6cc505380f7ad5b98482dd4a5a Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:19:28 +0200 Subject: [PATCH 06/19] cultmath: tests for Stryker survivors (dirty workspace, exact caps, spans, near-duplicates) Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 215 ++++++++++++++---- 1 file changed, 173 insertions(+), 42 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 2b3e870d..a7a640b5 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -104,6 +104,54 @@ public sealed class BoundedLeastSquaresTests return v; } + private static double[] NormalEquationSolution(int m, int n, float[] a, float[] b) + { + var mat = new double[n, n + 1]; + for (var i = 0; i < n; i++) + { + for (var j = 0; j < n; j++) + for (var r = 0; r < m; r++) mat[i, j] += (double)a[r * n + i] * a[r * n + j]; + for (var r = 0; r < m; r++) mat[i, n] += (double)a[r * n + i] * b[r]; + } + for (var c = 0; c < n; c++) + { + var piv = mat[c, c]; + for (var j = 0; j <= n; j++) mat[c, j] /= piv; + for (var i = 0; i < n; i++) + { + if (i == c) continue; + var f = mat[i, c]; + for (var j = 0; j <= n; j++) mat[i, j] -= f * mat[c, j]; + } + } + var x = new double[n]; + for (var i = 0; i < n; i++) x[i] = mat[i, n]; + return x; + } + + private static (int m, int n, float[] a, float[] b, float[] lo, float[] hi) RandomProblem(uint seed) + { + var rng = new CultMath.Random(seed * 7919u); + var n = rng.NextInt(2, 17); + var m = rng.NextInt(n, 25); + var a = new float[m * n]; + var b = new float[m]; + for (var i = 0; i < a.Length; i++) a[i] = rng.NextFloat(-1f, 1f); + for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) + { + switch (rng.NextInt(3)) + { + case 0: lo[j] = -1f; hi[j] = 1f; break; + case 1: lo[j] = 0f; hi[j] = 1f; break; + default: lo[j] = -0.25f; hi[j] = 0.1f; break; + } + } + return (m, n, a, b, lo, hi); + } + private static readonly float[] A54 = { 1.0f, 0.3f, -0.2f, @@ -123,26 +171,8 @@ public sealed class BoundedLeastSquaresTests var status = Solve(m, n, A54, B5, Fill(-100f, n), Fill(100f, n), x, out _); Assert.Equal(BoundedLeastSquaresStatus.Converged, status); - // Independent solve of (A^T A) x = A^T b by Gauss-Jordan in double. - var mat = new double[n, n + 1]; - for (var i = 0; i < n; i++) - { - for (var j = 0; j < n; j++) - for (var r = 0; r < m; r++) mat[i, j] += (double)A54[r * n + i] * A54[r * n + j]; - for (var r = 0; r < m; r++) mat[i, n] += (double)A54[r * n + i] * B5[r]; - } - for (var c = 0; c < n; c++) - { - var piv = mat[c, c]; - for (var j = 0; j <= n; j++) mat[c, j] /= piv; - for (var i = 0; i < n; i++) - { - if (i == c) continue; - var f = mat[i, c]; - for (var j = 0; j <= n; j++) mat[i, j] -= f * mat[c, j]; - } - } - for (var i = 0; i < n; i++) Assert.Equal(mat[i, n], x[i], 1e-4); + var expected = NormalEquationSolution(m, n, A54, B5); + for (var i = 0; i < n; i++) Assert.Equal(expected[i], x[i], 1e-4); } [Fact] @@ -151,24 +181,7 @@ public sealed class BoundedLeastSquaresTests var anyActive = false; for (uint seed = 1; seed <= 150; seed++) { - var rng = new CultMath.Random(seed * 7919u); - var n = rng.NextInt(2, 17); - var m = rng.NextInt(n, 25); - var a = new float[m * n]; - var b = new float[m]; - for (var i = 0; i < a.Length; i++) a[i] = rng.NextFloat(-1f, 1f); - for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); - var lo = new float[n]; - var hi = new float[n]; - for (var j = 0; j < n; j++) - { - switch (rng.NextInt(3)) - { - case 0: lo[j] = -1f; hi[j] = 1f; break; - case 1: lo[j] = 0f; hi[j] = 1f; break; - default: lo[j] = -0.25f; hi[j] = 0.1f; break; - } - } + var (m, n, a, b, lo, hi) = RandomProblem(seed); var x = new float[n]; var status = Solve(m, n, a, b, lo, hi, x, out _); Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed} ({m}x{n}) hit the cap"); @@ -191,9 +204,10 @@ public sealed class BoundedLeastSquaresTests var lo = Fill(-1f, n); var hi = Fill(1f, n); var x = new float[n]; - var status = Solve(m, n, a, b, lo, hi, x, out _); + var status = Solve(m, n, a, b, lo, hi, x, out var iterations); Assert.Equal(BoundedLeastSquaresStatus.Converged, status); Assert.Equal(new[] { 1f, -1f, 1f }, x); + Assert.Equal(3, iterations); AssertKkt(m, n, a, b, lo, hi, x); } @@ -336,9 +350,126 @@ public sealed class BoundedLeastSquaresTests } [Fact] - public void ShortSpansThrowArgumentException() + public void LargeInteriorOptimumEqualsNormalEquations() { - var x = new float[3]; - Assert.Throws(() => BoundedLeastSquares.Solve(5, 3, A54, B5, Fill(-1f, 3), Fill(1f, 3), x, new float[4], out _)); + const int m = 24, n = 16; + var rng = new CultMath.Random(4242u); + var a = new float[m * n]; + var b = new float[m]; + for (var i = 0; i < a.Length; i++) a[i] = rng.NextFloat(-1f, 1f); + for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + var x = new float[n]; + var status = Solve(m, n, a, b, Fill(-1000f, n), Fill(1000f, n), x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0, iterations); + var expected = NormalEquationSolution(m, n, a, b); + for (var i = 0; i < n; i++) Assert.Equal(expected[i], x[i], 1e-3); + } + + [Fact] + public void ScaledAndNearDuplicateColumnsAreOptimal() + { + const int m = 5, n = 3; + foreach (var noise in new[] { 0f, 1e-5f }) + { + var a = (float[])A54.Clone(); + for (var r = 0; r < m; r++) a[r * n + 2] = 0.37f * a[r * n + 0] + noise * (r - 2); + var lo = Fill(-1f, n); + var hi = new[] { 0.6f, 1f, 0.6f }; + var x = new float[n]; + var status = Solve(m, n, a, B5, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + AssertKkt(m, n, a, B5, lo, hi, x); + var cs = Cost(m, n, a, B5, x); + var cr = ReferenceCost(m, n, a, B5, lo, hi); + Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6, $"noise {noise}: solver cost {cs}, reference {cr}"); + } + } + + [Fact] + public void WorkspaceContentsOnEntryAreIgnored() + { + const int m = 5, n = 3; + var zeroCol = (float[])A54.Clone(); + for (var r = 0; r < m; r++) zeroCol[r * n + 1] = 0f; + var dupCol = (float[])A54.Clone(); + for (var r = 0; r < m; r++) dupCol[r * n + 2] = dupCol[r * n + 0]; + foreach (var a in new[] { A54, zeroCol, dupCol }) + { + var lo = new[] { -1f, -1f, -1f }; + var hi = new[] { 0.3f, 1f, 0.3f }; + var clean = new float[n]; + Solve(m, n, a, B5, lo, hi, clean, out var cleanIterations); + foreach (var junk in new[] { float.NaN, 7f, -1e30f }) + { + var workspace = new float[BoundedLeastSquares.WorkspaceLength(n)]; + Array.Fill(workspace, junk); + var x = new float[n]; + BoundedLeastSquares.Solve(m, n, a, B5, lo, hi, x, workspace, out var iterations); + Assert.Equal(clean, x); + Assert.Equal(cleanIterations, iterations); + } + } + } + + [Fact] + public void WarmStartAtAVertexTakesNoIterations() + { + const int m = 4, n = 3; + var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f }; + var b = new[] { 10f, -10f, 10f, 1f }; + var x = new[] { 1f, -1f, 1f }; + var status = Solve(m, n, a, b, Fill(-1f, n), Fill(1f, n), x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0, iterations); + Assert.Equal(new[] { 1f, -1f, 1f }, x); + } + + [Fact] + public void IterationCountFollowsTheCapExactly() + { + var sawRelease = false; + for (uint seed = 1; seed <= 40; seed++) + { + var (m, n, a, b, lo, hi) = RandomProblem(seed); + var full = new float[n]; + Solve(m, n, a, b, lo, hi, full, out var needed); + var previousCost = double.MaxValue; + for (var cap = 0; cap <= needed; cap++) + { + var x = new float[n]; + var status = Solve(m, n, a, b, lo, hi, x, out var iterations, cap); + Assert.Equal(cap, iterations); + Assert.Equal(cap < needed ? BoundedLeastSquaresStatus.IterationLimit : BoundedLeastSquaresStatus.Converged, status); + for (var j = 0; j < n; j++) Assert.InRange(x[j], lo[j], hi[j]); + var cost = Cost(m, n, a, b, x); + Assert.True(cost <= previousCost * (1 + 1e-6) + 1e-9, $"seed {seed} cap {cap}: cost rose from {previousCost} to {cost}"); + previousCost = cost; + } + sawRelease |= needed > n; + } + Assert.True(sawRelease, "no random problem needed more iterations than columns, so no bound was ever released"); + } + + [Fact] + public void EachShortSpanIsRejected() + { + const int m = 5, n = 3; + var ws = new float[BoundedLeastSquares.WorkspaceLength(n)]; + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + void Expect(Action call) => Assert.Contains("BoundedLeastSquares", Assert.Throws(call).Message); + Expect(() => BoundedLeastSquares.Solve(m, n, new float[m * n - 1], B5, lo, hi, new float[n], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, n, A54, new float[m - 1], lo, hi, new float[n], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, n, A54, B5, new float[n - 1], hi, new float[n], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, n, A54, B5, lo, new float[n - 1], new float[n], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, new float[n - 1], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, new float[n], new float[ws.Length - 1], out _)); + Expect(() => BoundedLeastSquares.Solve(-1, n, A54, B5, lo, hi, new float[n], ws, out _)); + Expect(() => BoundedLeastSquares.Solve(m, -1, A54, B5, lo, hi, new float[n], ws, out _)); + // Exactly sized spans, and empty problems, are accepted. + BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, new float[n], ws, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, BoundedLeastSquares.Solve(0, 0, default, default, default, default, default, default, out _)); + Assert.Equal(BoundedLeastSquaresStatus.Converged, BoundedLeastSquares.Solve(0, n, default, default, lo, hi, new float[n], ws, out _)); } } From 6f8cb642ff999b37eb6eedfa1ca0f69bbfe4e49c Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:43:32 +0200 Subject: [PATCH 07/19] cultmath: InvalidInput status and Soul tests (red on the old solver) Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 3 + .../BoundedLeastSquaresTests.cs | 213 ++++++++++++++++++ 2 files changed, 216 insertions(+) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index e18a42e2..112c1418 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -10,6 +10,9 @@ public enum BoundedLeastSquaresStatus /// The iteration cap was reached. x holds the last feasible iterate, which is never worse than the clamped start. IterationLimit, + + /// The input broke the contract; nothing was solved and x is exactly as passed. + InvalidInput, } /// diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index a7a640b5..e55dd64e 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -472,4 +472,217 @@ public sealed class BoundedLeastSquaresTests Assert.Equal(BoundedLeastSquaresStatus.Converged, BoundedLeastSquares.Solve(0, 0, default, default, default, default, default, default, out _)); Assert.Equal(BoundedLeastSquaresStatus.Converged, BoundedLeastSquares.Solve(0, n, default, default, lo, hi, new float[n], ws, out _)); } + + // Well-conditioned enough (m = n + 6) that x is comparable across scales. + private static (int m, int n, float[] a, float[] b, float[] lo, float[] hi) ScaleProblem(uint seed) + { + var rng = new CultMath.Random(seed * 104729u); + var n = rng.NextInt(3, 11); + var m = n + 6; + var a = new float[m * n]; + var b = new float[m]; + for (var i = 0; i < a.Length; i++) a[i] = rng.NextFloat(-1f, 1f); + for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) + { + lo[j] = rng.NextFloat(-0.6f, -0.05f); + hi[j] = rng.NextFloat(0.05f, 0.6f); + } + return (m, n, a, b, lo, hi); + } + + private static float[] Scaled(float[] v, float k) + { + var r = new float[v.Length]; + for (var i = 0; i < v.Length; i++) r[i] = v[i] * k; + return r; + } + + [Fact] + public void SolutionIsScaleFree() + { + var sawActive = false; + for (uint seed = 1; seed <= 40; seed++) + { + var (m, n, a, b, lo, hi) = ScaleProblem(seed); + var x = (float[])lo.Clone(); + Assert.Equal(BoundedLeastSquaresStatus.Converged, Solve(m, n, a, b, lo, hi, x, out _)); + AssertKkt(m, n, a, b, lo, hi, x); + for (var j = 0; j < n; j++) sawActive |= x[j] == lo[j] || x[j] == hi[j]; + + foreach (var k in new[] { 1e-6f, 1e3f, 1e6f }) + { + // b and the bounds scaled: x scales with them. + var xb = Scaled(lo, k); + var status = Solve(m, n, a, Scaled(b, k), Scaled(lo, k), Scaled(hi, k), xb, out _); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed} k {k}: {status}"); + for (var j = 0; j < n; j++) + Assert.True(Math.Abs(xb[j] - x[j] * k) <= 1e-3 * Math.Abs(k) * (1 + Math.Abs(x[j])), $"seed {seed} b*{k} col {j}: {xb[j]} vs {x[j] * k}"); + + // A scaled: x scales inversely, so the bounds stay put and only b is scaled with A to keep the same box. + var xa = (float[])lo.Clone(); + var lo2 = Scaled(lo, 1f / k); + var hi2 = Scaled(hi, 1f / k); + for (var j = 0; j < n; j++) xa[j] = lo2[j]; + status = Solve(m, n, Scaled(a, k), b, lo2, hi2, xa, out _); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed} A*{k}: {status}"); + for (var j = 0; j < n; j++) + Assert.True(Math.Abs(xa[j] - x[j] / k) <= 1e-3 * Math.Abs(1 / k) * (1 + Math.Abs(x[j])), $"seed {seed} A*{k} col {j}: {xa[j]} vs {x[j] / k}"); + } + } + Assert.True(sawActive); + } + + [Fact] + public void NonFiniteOrInvertedInputIsRejectedAndLeavesXUntouched() + { + const int m = 5, n = 3; + var nan = float.NaN; + var inf = float.PositiveInfinity; + var cases = new (string name, Action corrupt)[] + { + ("A NaN", (a, b, lo, hi) => a[4] = nan), + ("A +Inf", (a, b, lo, hi) => a[7] = inf), + ("A -Inf", (a, b, lo, hi) => a[1] = -inf), + ("b NaN", (a, b, lo, hi) => b[2] = nan), + ("b +Inf", (a, b, lo, hi) => b[0] = inf), + ("lo NaN", (a, b, lo, hi) => lo[1] = nan), + ("hi NaN", (a, b, lo, hi) => hi[2] = nan), + ("lo > hi", (a, b, lo, hi) => { lo[1] = 0.5f; hi[1] = 0.25f; }), + ("lo = +Inf", (a, b, lo, hi) => { lo[0] = inf; hi[0] = inf; }), + ("hi = -Inf", (a, b, lo, hi) => { lo[0] = -inf; hi[0] = -inf; }), + }; + foreach (var (name, corrupt) in cases) + { + var a = (float[])A54.Clone(); + var b = (float[])B5.Clone(); + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + corrupt(a, b, lo, hi); + foreach (var start in new[] { new[] { 0.5f, -3f, 9f }, new[] { nan, 0f, 0f } }) + { + var x = (float[])start.Clone(); + var status = Solve(m, n, a, b, lo, hi, x, out var iterations); + Assert.True(status == BoundedLeastSquaresStatus.InvalidInput, $"{name}: {status}"); + Assert.Equal(0, iterations); + for (var j = 0; j < n; j++) + Assert.Equal(BitConverter.SingleToInt32Bits(start[j]), BitConverter.SingleToInt32Bits(x[j])); + } + } + } + + [Fact] + public void InfiniteBoundsMeanUnbounded() + { + const int m = 5, n = 3; + var inf = float.PositiveInfinity; + var x = new float[n]; + Assert.Equal(BoundedLeastSquaresStatus.Converged, Solve(m, n, A54, B5, Fill(-inf, n), Fill(inf, n), x, out _)); + var expected = NormalEquationSolution(m, n, A54, B5); + for (var j = 0; j < n; j++) Assert.Equal(expected[j], x[j], 1e-4); + + // One-sided, with an infinite warm start: the start is clamped to a finite point first. + var y = new[] { -inf, inf, float.NaN }; + var lo = new[] { -inf, 0.2f, -inf }; + var hi = new[] { inf, inf, 0.1f }; + Assert.Equal(BoundedLeastSquaresStatus.Converged, Solve(m, n, A54, B5, lo, hi, y, out _)); + AssertKkt(m, n, A54, B5, lo, hi, y); + Assert.All(y, v => Assert.True(float.IsFinite(v))); + } + + // Exact reference: every assignment of each column to {lo, hi, free}, free set solved in double. + private static double ExactBoxedCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi) + { + var best = double.MaxValue; + var total = 1; + for (var j = 0; j < n; j++) total *= 3; + var x = new double[n]; + for (var code = 0; code < total; code++) + { + var c = code; + var free = new System.Collections.Generic.List(); + for (var j = 0; j < n; j++) + { + var d = c % 3; c /= 3; + if (d == 0) x[j] = lo[j]; else if (d == 1) x[j] = hi[j]; else free.Add(j); + } + var k = free.Count; + if (k > 0) + { + var mat = new double[k, k + 1]; + for (var i = 0; i < k; i++) + { + for (var l = 0; l < k; l++) + for (var r = 0; r < m; r++) mat[i, l] += (double)a[r * n + free[i]] * a[r * n + free[l]]; + for (var r = 0; r < m; r++) + { + var fixedPart = 0.0; + for (var j = 0; j < n; j++) if (!free.Contains(j)) fixedPart += (double)a[r * n + j] * x[j]; + mat[i, k] += (double)a[r * n + free[i]] * (b[r] - fixedPart); + } + } + var ok = true; + for (var col = 0; col < k && ok; col++) + { + var piv = col; + for (var r = col + 1; r < k; r++) if (Math.Abs(mat[r, col]) > Math.Abs(mat[piv, col])) piv = r; + if (Math.Abs(mat[piv, col]) < 1e-9) { ok = false; break; } + for (var l = 0; l <= k; l++) (mat[col, l], mat[piv, l]) = (mat[piv, l], mat[col, l]); + for (var r = 0; r < k; r++) + { + if (r == col) continue; + var f = mat[r, col] / mat[col, col]; + for (var l = col; l <= k; l++) mat[r, l] -= f * mat[col, l]; + } + } + if (!ok) continue; + for (var i = 0; i < k; i++) + { + x[free[i]] = mat[i, k] / mat[i, i]; + if (x[free[i]] < lo[free[i]] - 1e-12 || x[free[i]] > hi[free[i]] + 1e-12) ok = false; + } + if (!ok) continue; + } + var cost = 0.0; + for (var r = 0; r < m; r++) + { + var s = -(double)b[r]; + for (var j = 0; j < n; j++) s += (double)a[r * n + j] * x[j]; + cost += s * s; + } + best = Math.Min(best, cost); + } + return best; + } + + [Fact] + public void NearCollinearColumnsReachTheExactOptimum() + { + const int m = 2, n = 7; + var worst = 0.0; + for (uint seed = 1; seed <= 60; seed++) + { + var rng = new CultMath.Random(seed * 15485863u); + var a = new float[m * n]; + for (var j = 0; j < n; j++) + { + a[j] = rng.NextFloat(0.5f, 1f); + a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); + } + var b = new[] { 6f * (rng.NextFloat() < 0.5f ? -1f : 1f), rng.NextFloat(0.5f, 2f) }; + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) { lo[j] = -1f; hi[j] = seed % 2 == 0 ? 1f : 0.5f; } + var x = new float[n]; + var status = Solve(m, n, a, b, lo, hi, x, out _); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed}: {status}"); + AssertKkt(m, n, a, b, lo, hi, x); + var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi); + var bound = 1e-5 * (b[0] * b[0] + b[1] * b[1]); + worst = Math.Max(worst, gap / bound); + Assert.True(gap <= bound, $"seed {seed}: cost {gap} above the exact optimum (allowed {bound})"); + } + } } From 8b1a1ffde091e552db7d34bb0d80fc1fc6abf085 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:44:57 +0200 Subject: [PATCH 08/19] cultmath: near-collinear test uses generators that reproduce the gap on the old solver Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 44 +++++++++++-------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index e55dd64e..0addfddf 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -657,32 +657,38 @@ public sealed class BoundedLeastSquaresTests return best; } + // m < n with columns that are nearly parallel (condition number of A around 1e2), b outside the reachable set. [Fact] public void NearCollinearColumnsReachTheExactOptimum() { const int m = 2, n = 7; - var worst = 0.0; - for (uint seed = 1; seed <= 60; seed++) + for (var family = 0; family < 4; family++) { - var rng = new CultMath.Random(seed * 15485863u); - var a = new float[m * n]; - for (var j = 0; j < n; j++) + for (uint seed = 1; seed <= 120; seed++) { - a[j] = rng.NextFloat(0.5f, 1f); - a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); + var rng = new CultMath.Random(seed * 15485863u + (uint)family * 977u); + var a = new float[m * n]; + for (var j = 0; j < n; j++) + { + switch (family) + { + case 0: a[j] = rng.NextFloat(0.5f, 1f); a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); break; + case 1: a[j] = 1f + 0.02f * rng.NextFloat(-1f, 1f); a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); break; + case 2: a[j] = rng.NextFloat(-1f, 1f); a[n + j] = 0.01f * a[j] + 0.01f * rng.NextFloat(-1f, 1f); break; + default: a[j] = rng.NextFloat(0.5f, 1f); a[n + j] = a[j] * 0.9f + 0.01f * rng.NextFloat(-1f, 1f); break; + } + } + var b = new[] { rng.NextFloat(-8f, 8f), rng.NextFloat(-8f, 8f) }; + var lo = new float[n]; + var hi = new float[n]; + for (var j = 0; j < n; j++) { lo[j] = rng.NextFloat() < 0.5f ? -1f : 0f; hi[j] = 1f; } + var x = new float[n]; + var status = Solve(m, n, a, b, lo, hi, x, out _); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"family {family} seed {seed}: {status}"); + var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi); + var allowed = 1e-5 * (b[0] * b[0] + b[1] * b[1]); + Assert.True(gap <= allowed, $"family {family} seed {seed}: cost {gap} above the exact optimum (allowed {allowed})"); } - var b = new[] { 6f * (rng.NextFloat() < 0.5f ? -1f : 1f), rng.NextFloat(0.5f, 2f) }; - var lo = new float[n]; - var hi = new float[n]; - for (var j = 0; j < n; j++) { lo[j] = -1f; hi[j] = seed % 2 == 0 ? 1f : 0.5f; } - var x = new float[n]; - var status = Solve(m, n, a, b, lo, hi, x, out _); - Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed}: {status}"); - AssertKkt(m, n, a, b, lo, hi, x); - var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi); - var bound = 1e-5 * (b[0] * b[0] + b[1] * b[1]); - worst = Math.Max(worst, gap / bound); - Assert.True(gap <= bound, $"seed {seed}: cost {gap} above the exact optimum (allowed {bound})"); } } } From 44affbb9b65125ad6005daed50c2a8f2efd9241d Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:46:36 +0200 Subject: [PATCH 09/19] cultmath: scale-free KKT, InvalidInput contract, free-set refinement Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 126 ++++++++++++++---- 1 file changed, 97 insertions(+), 29 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index 112c1418..5730e8a1 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -25,8 +25,18 @@ public enum BoundedLeastSquaresStatus /// 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 and run in float: callers that need a minimum-norm answer on -/// rank-deficient problems add a small ridge row to A and b. The solver allocates nothing and is -/// deterministic. +/// rank-deficient problems add a small ridge row to A and b; 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. +/// +/// Input contract: a and b must be finite; lo and hi must not be NaN and need +/// lo <= hi. Infinite bounds are allowed and mean unbounded on that side (lo = -Inf, +/// hi = +Inf); lo = +Inf or hi = -Inf is invalid. Any violation returns +/// without throwing, with iterations 0 and +/// x exactly as passed. A NaN or infinite warm start is legal: it is clamped into the box (an infinite +/// result on an unbounded side becomes 0 clamped to the box). Spans shorter than the stated dimensions are a +/// caller bug and throw . +/// /// public static class BoundedLeastSquares { @@ -37,7 +47,7 @@ public static class BoundedLeastSquares private const float KktRelativeTolerance = 1e-5f; /// Floats of caller-supplied workspace needed for columns. - public static int WorkspaceLength(int n) => 2 * n * n + 7 * n; + public static int WorkspaceLength(int n) => 2 * n * n + 9 * n; /// /// Minimises ||A x - b||^2 subject to lo <= x <= hi. is the warm start @@ -47,8 +57,8 @@ public static class BoundedLeastSquares /// Columns of , and length of , and . /// Row-major m * n matrix. /// Length m target. - /// Lower bounds; each must be <= hi. - /// Upper bounds. + /// Lower bounds; each must be <= hi; -Inf means unbounded below. + /// Upper bounds; +Inf means unbounded above. /// Warm start in, solution out. /// At least floats; contents on entry are ignored. /// Active-set iterations used. @@ -72,21 +82,24 @@ public static class BoundedLeastSquares var g = workspace.Slice(off + n, n); var p = workspace.Slice(off + 2 * n, n); var y = workspace.Slice(off + 3 * n, n); + var r1 = workspace.Slice(off + 4 * n, n); + var q = workspace.Slice(off + 5 * n, n); // One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow. - var ints = MemoryMarshal.Cast(workspace.Slice(off + 4 * n, 3 * n)); + var ints = MemoryMarshal.Cast(workspace.Slice(off + 6 * n, 3 * n)); var state = ints.Slice(0, n); var freeIndex = ints.Slice(n, n); var dependent = ints.Slice(2 * n, n); + iterations = 0; for (var j = 0; j < n; j++) - { - var l = lo[j]; - var h = hi[j]; - var v = x[j]; - v = !(v > l) ? l : (v < h ? v : h); - x[j] = v; - state[j] = v == l ? -1 : (v == h ? 1 : 0); - } + if (!(lo[j] <= hi[j]) || float.IsPositiveInfinity(lo[j]) || float.IsNegativeInfinity(hi[j])) + return BoundedLeastSquaresStatus.InvalidInput; + for (var i = 0; i < m * n; i++) + if (!float.IsFinite(a[i])) + return BoundedLeastSquaresStatus.InvalidInput; + for (var i = 0; i < m; i++) + if (!float.IsFinite(b[i])) + return BoundedLeastSquaresStatus.InvalidInput; var atbMax = 0f; for (var j = 0; j < n; j++) @@ -98,16 +111,41 @@ public static class BoundedLeastSquares atbMax = MathF.Max(atbMax, MathF.Abs(sb)); for (var k = 0; k <= j; k++) { - var s = 0f; + var sa = 0f; for (var r = 0; r < m; r++) - s += a[r * n + j] * a[r * n + k]; - ata[j * n + k] = s; - ata[k * n + j] = s; + sa += a[r * n + j] * a[r * n + k]; + if (!float.IsFinite(sa) || !float.IsFinite(sb)) + return BoundedLeastSquaresStatus.InvalidInput; + ata[j * n + k] = sa; + ata[k * n + j] = sa; } } - var kktTol = KktRelativeTolerance * (1f + atbMax); - iterations = 0; + // 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. + for (var j = 0; j < n; j++) + { + var l = lo[j]; + var h = hi[j]; + var v = x[j]; + v = !(v > l) ? l : (v < h ? v : h); + if (float.IsInfinity(v)) + v = 0f < l ? l : (0f > h ? h : 0f); + x[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 == 0f) + { + Gradient(ata, atb, x, g, n); + for (var j = 0; j < n; j++) + gScale = MathF.Max(gScale, MathF.Abs(g[j])); + } + var kktTol = KktRelativeTolerance * gScale; + while (true) { Gradient(ata, atb, x, g, n); @@ -116,7 +154,7 @@ public static class BoundedLeastSquares for (var j = 0; j < n; j++) if (state[j] == 0) freeIndex[freeCount++] = j; - SolveFreeSet(ata, chol, g, p, y, freeIndex.Slice(0, freeCount), dependent, n); + SolveFreeSet(ata, chol, g, p, y, r1, q, freeIndex.Slice(0, freeCount), dependent, n); var t = 1f; var block = -1; @@ -155,6 +193,12 @@ public static class BoundedLeastSquares } Gradient(ata, atb, x, g, n); + // Stationarity of the free set is part of KKT: float Cholesky can leave a residual gradient. + var freeViolated = false; + for (var i = 0; i < freeCount; i++) + if (dependent[i] == 0 && MathF.Abs(g[freeIndex[i]]) > kktTol) + freeViolated = true; + var release = -1; var worst = kktTol; for (var j = 0; j < n; j++) @@ -167,10 +211,11 @@ public static class BoundedLeastSquares release = j; } } - if (release < 0) return BoundedLeastSquaresStatus.Converged; + if (release < 0 && !freeViolated) return BoundedLeastSquaresStatus.Converged; if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit; iterations++; - state[release] = 0; + // A stationarity miss re-solves the same free set from the improved point. + if (release >= 0) state[release] = 0; } } @@ -185,11 +230,12 @@ public static class BoundedLeastSquares } } - // Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky. A free column whose pivot - // collapses (dependent on earlier free columns) is dropped from the factor and gets p = 0. + // Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky, then one step of iterative + // refinement with the same factor. A free column whose pivot collapses (dependent on earlier free columns) + // is dropped from the factor and gets p = 0. private static void SolveFreeSet( ReadOnlySpan ata, Span l, ReadOnlySpan g, Span p, Span y, - ReadOnlySpan free, Span dependent, int n) + Span rhs, Span dp, ReadOnlySpan free, Span dependent, int n) { var k = free.Length; for (var i = 0; i < k; i++) @@ -218,9 +264,31 @@ public static class BoundedLeastSquares } } + for (var i = 0; i < k; i++) + rhs[i] = dependent[i] != 0 ? 0f : -g[free[i]]; + Substitute(l, rhs, y, p, free, n); + for (var i = 0; i < k; i++) { - var s = dependent[i] != 0 ? 0f : -g[free[i]]; + if (dependent[i] != 0) { rhs[i] = 0f; continue; } + var s = -g[free[i]]; + for (var c = 0; c < k; c++) + if (dependent[c] == 0) + s -= ata[free[i] * n + free[c]] * p[free[c]]; + rhs[i] = s; + } + Substitute(l, rhs, y, dp, free, n); + for (var i = 0; i < k; i++) + p[free[i]] += dp[free[i]]; + } + + // Solves L L^T out = rhs; rhs is indexed by free-list position, out by column. + private static void Substitute(ReadOnlySpan l, ReadOnlySpan rhs, Span y, Span outByColumn, ReadOnlySpan free, int n) + { + var k = free.Length; + for (var i = 0; i < k; i++) + { + var s = rhs[i]; for (var c = 0; c < i; c++) s -= l[i * n + c] * y[c]; y[i] = s / l[i * n + i]; @@ -229,8 +297,8 @@ public static class BoundedLeastSquares { var s = y[i]; for (var r = i + 1; r < k; r++) - s -= l[r * n + i] * p[free[r]]; - p[free[i]] = s / l[i * n + i]; + s -= l[r * n + i] * outByColumn[free[r]]; + outByColumn[free[i]] = s / l[i * n + i]; } } } From 50f26de4a75a9e0579ef47ce3da8b2f8eb4a88e6 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:49:29 +0200 Subject: [PATCH 10/19] cultmath: BoundedLeastSquares factors the normal equations in double Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 87 ++++++++++--------- 1 file changed, 45 insertions(+), 42 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index 5730e8a1..25900bfe 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -24,9 +24,9 @@ public enum BoundedLeastSquaresStatus /// 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 and run in float: callers that need a minimum-norm answer on -/// rank-deficient problems add a small ridge row to A and b; 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 +/// equations square the condition number, so they are formed and factored in double (inputs and x stay +/// float): callers that need a minimum-norm answer on rank-deficient problems add a small ridge row to +/// A and b; 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. /// /// Input contract: a and b must be finite; lo and hi must not be NaN and need @@ -43,11 +43,11 @@ public static class BoundedLeastSquares /// Default iteration cap for . public const int DefaultMaxIterations = 100; - private const float PivotRelativeTolerance = 1e-6f; - private const float KktRelativeTolerance = 1e-5f; + private const double PivotRelativeTolerance = 1e-12; + private const double KktRelativeTolerance = 1e-6; /// Floats of caller-supplied workspace needed for columns. - public static int WorkspaceLength(int n) => 2 * n * n + 9 * n; + public static int WorkspaceLength(int n) => 2 * (2 * n * n + 6 * n) + 3 * n; /// /// Minimises ||A x - b||^2 subject to lo <= x <= hi. is the warm start @@ -75,17 +75,20 @@ public static class BoundedLeastSquares || x.Length < n || workspace.Length < WorkspaceLength(n)) throw new ArgumentException("BoundedLeastSquares: span shorter than the stated dimensions."); - var ata = workspace.Slice(0, n * n); - var chol = workspace.Slice(n * n, n * n); + // 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 dbl = MemoryMarshal.Cast(workspace.Slice(0, 2 * doubleCount)); + var ata = dbl.Slice(0, n * n); + var chol = dbl.Slice(n * n, n * n); var off = 2 * n * n; - var atb = workspace.Slice(off, n); - var g = workspace.Slice(off + n, n); - var p = workspace.Slice(off + 2 * n, n); - var y = workspace.Slice(off + 3 * n, n); - var r1 = workspace.Slice(off + 4 * n, n); - var q = workspace.Slice(off + 5 * n, n); + var atb = dbl.Slice(off, n); + var g = dbl.Slice(off + n, n); + var p = dbl.Slice(off + 2 * n, n); + var y = dbl.Slice(off + 3 * n, n); + var r1 = dbl.Slice(off + 4 * n, n); + var q = dbl.Slice(off + 5 * n, n); // One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow. - var ints = MemoryMarshal.Cast(workspace.Slice(off + 6 * n, 3 * n)); + var ints = MemoryMarshal.Cast(workspace.Slice(2 * doubleCount, 3 * n)); var state = ints.Slice(0, n); var freeIndex = ints.Slice(n, n); var dependent = ints.Slice(2 * n, n); @@ -101,20 +104,20 @@ public static class BoundedLeastSquares if (!float.IsFinite(b[i])) return BoundedLeastSquaresStatus.InvalidInput; - var atbMax = 0f; + var atbMax = 0.0; for (var j = 0; j < n; j++) { - var sb = 0f; + var sb = 0.0; for (var r = 0; r < m; r++) - sb += a[r * n + j] * b[r]; + sb += (double)a[r * n + j] * b[r]; atb[j] = sb; - atbMax = MathF.Max(atbMax, MathF.Abs(sb)); + atbMax = Math.Max(atbMax, Math.Abs(sb)); for (var k = 0; k <= j; k++) { - var sa = 0f; + var sa = 0.0; for (var r = 0; r < m; r++) - sa += a[r * n + j] * a[r * n + k]; - if (!float.IsFinite(sa) || !float.IsFinite(sb)) + sa += (double)a[r * n + j] * a[r * n + k]; + if (!double.IsFinite(sa) || !double.IsFinite(sb)) return BoundedLeastSquaresStatus.InvalidInput; ata[j * n + k] = sa; ata[k * n + j] = sa; @@ -138,11 +141,11 @@ public static class BoundedLeastSquares // 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 == 0f) + if (gScale == 0.0) { Gradient(ata, atb, x, g, n); for (var j = 0; j < n; j++) - gScale = MathF.Max(gScale, MathF.Abs(g[j])); + gScale = Math.Max(gScale, Math.Abs(g[j])); } var kktTol = KktRelativeTolerance * gScale; @@ -156,13 +159,13 @@ public static class BoundedLeastSquares freeIndex[freeCount++] = j; SolveFreeSet(ata, chol, g, p, y, r1, q, freeIndex.Slice(0, freeCount), dependent, n); - var t = 1f; + var t = 1.0; var block = -1; for (var i = 0; i < freeCount; i++) { var j = freeIndex[i]; var pj = p[j]; - if (pj == 0f) continue; + if (pj == 0.0) continue; var tj = pj > 0f ? (hi[j] - x[j]) / pj : (lo[j] - x[j]) / pj; if (tj < t) { @@ -178,9 +181,9 @@ public static class BoundedLeastSquares for (var i = 0; i < freeCount; i++) { var j = freeIndex[i]; - x[j] += t * p[j]; + x[j] = (float)(x[j] + t * p[j]); } - if (p[block] > 0f) { x[block] = hi[block]; state[block] = 1; } + if (p[block] > 0.0) { x[block] = hi[block]; state[block] = 1; } else { x[block] = lo[block]; state[block] = -1; } continue; } @@ -188,7 +191,7 @@ public static class BoundedLeastSquares for (var i = 0; i < freeCount; i++) { var j = freeIndex[i]; - var v = x[j] + p[j]; + var v = (float)(x[j] + p[j]); x[j] = v < lo[j] ? lo[j] : (v > hi[j] ? hi[j] : v); } Gradient(ata, atb, x, g, n); @@ -196,7 +199,7 @@ public static class BoundedLeastSquares // Stationarity of the free set is part of KKT: float Cholesky can leave a residual gradient. var freeViolated = false; for (var i = 0; i < freeCount; i++) - if (dependent[i] == 0 && MathF.Abs(g[freeIndex[i]]) > kktTol) + if (dependent[i] == 0 && Math.Abs(g[freeIndex[i]]) > kktTol) freeViolated = true; var release = -1; @@ -204,7 +207,7 @@ public static class BoundedLeastSquares for (var j = 0; j < n; j++) { if (lo[j] == hi[j]) continue; - var violation = state[j] < 0 ? -g[j] : (state[j] > 0 ? g[j] : 0f); + var violation = state[j] < 0 ? -g[j] : (state[j] > 0 ? g[j] : 0.0); if (violation > worst) { worst = violation; @@ -219,7 +222,7 @@ public static class BoundedLeastSquares } } - private static void Gradient(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, Span g, int n) + private static void Gradient(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, Span g, int n) { for (var j = 0; j < n; j++) { @@ -234,8 +237,8 @@ public static class BoundedLeastSquares // refinement with the same factor. A free column whose pivot collapses (dependent on earlier free columns) // is dropped from the factor and gets p = 0. private static void SolveFreeSet( - ReadOnlySpan ata, Span l, ReadOnlySpan g, Span p, Span y, - Span rhs, Span dp, ReadOnlySpan free, Span dependent, int n) + ReadOnlySpan ata, Span l, ReadOnlySpan g, Span p, Span y, + Span rhs, Span dp, ReadOnlySpan free, Span dependent, int n) { var k = free.Length; for (var i = 0; i < k; i++) @@ -244,16 +247,16 @@ 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 <= 0f) + if (d <= PivotRelativeTolerance * ata[fi * n + fi] || d <= 0.0) { dependent[i] = 1; - for (var c = 0; c < i; c++) l[i * n + c] = 0f; - l[i * n + i] = 1f; - for (var r = i + 1; r < k; r++) l[r * n + i] = 0f; + for (var c = 0; c < i; c++) l[i * n + c] = 0.0; + l[i * n + i] = 1.0; + for (var r = i + 1; r < k; r++) l[r * n + i] = 0.0; continue; } dependent[i] = 0; - var lii = MathF.Sqrt(d); + var lii = Math.Sqrt(d); l[i * n + i] = lii; for (var r = i + 1; r < k; r++) { @@ -265,12 +268,12 @@ public static class BoundedLeastSquares } for (var i = 0; i < k; i++) - rhs[i] = dependent[i] != 0 ? 0f : -g[free[i]]; + rhs[i] = dependent[i] != 0 ? 0.0 : -g[free[i]]; Substitute(l, rhs, y, p, free, n); for (var i = 0; i < k; i++) { - if (dependent[i] != 0) { rhs[i] = 0f; continue; } + if (dependent[i] != 0) { rhs[i] = 0.0; continue; } var s = -g[free[i]]; for (var c = 0; c < k; c++) if (dependent[c] == 0) @@ -283,7 +286,7 @@ public static class BoundedLeastSquares } // Solves L L^T out = rhs; rhs is indexed by free-list position, out by column. - private static void Substitute(ReadOnlySpan l, ReadOnlySpan rhs, Span y, Span outByColumn, ReadOnlySpan free, int n) + private static void Substitute(ReadOnlySpan l, ReadOnlySpan rhs, Span y, Span outByColumn, ReadOnlySpan free, int n) { var k = free.Length; for (var i = 0; i < k; i++) From 43035bc827ac42d9928d15ddea57ef088f466fb4 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 05:36:50 +0200 Subject: [PATCH 11/19] cultmath: BoundedLeastSquares iterates in double and stops on progress 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 --- .../src/CultMath/BoundedLeastSquares.cs | 113 +++++++++++++----- 1 file changed, 82 insertions(+), 31 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index 25900bfe..dfda5516 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -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 x stay -/// float): callers that need a minimum-norm answer on rank-deficient problems add a small ridge row to -/// A and b; 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. +/// +/// Precision. A, b and x 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. x is converted in at the start +/// and back out at the end, so a returned float x is the nearest float to the double solution and always +/// inside the box. Each free-set solve gets one step of iterative refinement. +/// +/// +/// 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 +/// , so an optimal point whose double gradient sits above the +/// tolerance costs one confirming iteration rather than the iteration cap. +/// +/// +/// 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 +/// A and b; 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. +/// +/// The solver allocates nothing and is deterministic. /// /// Input contract: a and b must be finite; lo and hi must not be NaN and need /// lo <= hi. Infinite bounds are allowed and mean unbounded on that side (lo = -Inf, @@ -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; /// Floats of caller-supplied workspace needed for columns. - 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; /// /// Minimises ||A x - b||^2 subject to lo <= x <= hi. 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(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(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 ata, ReadOnlySpan atb, ReadOnlySpan x, Span 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 x, ReadOnlySpan g, ReadOnlySpan 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 ata, ReadOnlySpan atb, ReadOnlySpan x, Span 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; From 6982adfaafde22e370b9759bfd03032b5b9b31a4 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 05:36:50 +0200 Subject: [PATCH 12/19] cultmath: BoundedLeastSquares tests pin iteration counts and the b = 0 start - n=24 fully free problems at cond 1e2/1e3/1e4, and A x 1e6 with unscaled bounds, converge in at most two iterations (a regression to cap-spinning fails on the count, not on a clock). - b = 0 from a far start reaches the optimum instead of the start. - Warm-start clamp table, free duplicate column gets no step, parallel columns at absolute scales 1e3 and 1e-3, near-singular free sets confirm once and stop, near-collinear optimum at A x 1e6, and zero allocation on the invalid, converged and cap paths (first and second call). Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 253 +++++++++++++++++- 1 file changed, 251 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 0addfddf..03177d80 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -658,8 +658,12 @@ public sealed class BoundedLeastSquaresTests } // m < n with columns that are nearly parallel (condition number of A around 1e2), b outside the reachable set. - [Fact] - public void NearCollinearColumnsReachTheExactOptimum() + // The absolute scale of A must not matter: the KKT tolerance is relative, and the working iterate is double, so a + // 1e6 A with untouched bounds reaches the same optimum to the same allowance. + [Theory] + [InlineData(1f)] + [InlineData(1e6f)] + public void NearCollinearColumnsReachTheExactOptimum(float scale) { const int m = 2, n = 7; for (var family = 0; family < 4; family++) @@ -678,6 +682,7 @@ public sealed class BoundedLeastSquaresTests default: a[j] = rng.NextFloat(0.5f, 1f); a[n + j] = a[j] * 0.9f + 0.01f * rng.NextFloat(-1f, 1f); break; } } + for (var i = 0; i < a.Length; i++) a[i] *= scale; var b = new[] { rng.NextFloat(-8f, 8f), rng.NextFloat(-8f, 8f) }; var lo = new float[n]; var hi = new float[n]; @@ -691,4 +696,248 @@ public sealed class BoundedLeastSquaresTests } } } + + // A = U diag(s) V^T with s geometric from 1 down to 1/cond, so cond(A) is the requested value (float rounding of A + // perturbs it slightly). m x n with m >= n; b is generic, so it has a real residual when m > n. + private static (float[] a, float[] b) ConditionedProblem(int m, int n, double cond, uint seed, double aScale = 1.0) + { + var rng = new CultMath.Random(seed * 2654435761u + 17u); + double[][] Orthonormal(int dim, int count) + { + var q = new double[count][]; + for (var k = 0; k < count; k++) + { + var v = new double[dim]; + for (var i = 0; i < dim; i++) v[i] = rng.NextFloat(-1f, 1f); + for (var pass = 0; pass < 2; pass++) + for (var l = 0; l < k; l++) + { + var dot = 0.0; + for (var i = 0; i < dim; i++) dot += v[i] * q[l][i]; + for (var i = 0; i < dim; i++) v[i] -= dot * q[l][i]; + } + var norm = 0.0; + for (var i = 0; i < dim; i++) norm += v[i] * v[i]; + norm = Math.Sqrt(norm); + for (var i = 0; i < dim; i++) v[i] /= norm; + q[k] = v; + } + return q; + } + var u = Orthonormal(m, n); + var w = Orthonormal(n, n); + var a = new float[m * n]; + for (var k = 0; k < n; k++) + { + var s = aScale * Math.Pow(cond, -(double)k / Math.Max(1, n - 1)); + for (var r = 0; r < m; r++) + for (var c = 0; c < n; c++) + a[r * n + c] += (float)(u[k][r] * s * w[k][c]); + } + var b = new float[m]; + for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + return (a, b); + } + + // Unsaturated allocator-shaped ticks: 24 columns, every column free at the optimum. The solve is one Newton step; + // at most one more may confirm it. Before the working iterate moved to double, the float gradient's noise floor + // exceeded the tolerance from cond ~1e3 and the solver spun to the 100-iteration cap on optimal points. + [Theory] + [InlineData(1e2)] + [InlineData(1e3)] + [InlineData(1e4)] + public void UnsaturatedIllConditionedProblemsConvergeInAtMostTwoIterations(double cond) + { + const int m = 28, n = 24; + for (uint seed = 1; seed <= 30; seed++) + { + var (a, b) = ConditionedProblem(m, n, cond, seed); + var x = new float[n]; + var status = Solve(m, n, a, b, Fill(-1e9f, n), Fill(1e9f, n), x, out var iterations); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"cond {cond} seed {seed}: {status} after {iterations}"); + Assert.True(iterations <= 2, $"cond {cond} seed {seed}: {iterations} iterations"); + var bb = 0.0; + foreach (var v in b) bb += (double)v * v; + var expected = NormalEquationSolution(m, n, a, b); + var xs = new float[n]; + for (var j = 0; j < n; j++) xs[j] = (float)expected[j]; + Assert.True(Cost(m, n, a, b, x) - Cost(m, n, a, b, xs) <= 1e-6 * bb, $"cond {cond} seed {seed}: cost above the exact optimum"); + } + } + + // Near the pivot cutoff (cond(A) ~ 3e7, cond(A^T A) ~ 1e15) a Newton step leaves a free gradient above the + // tolerance on some problems. Those take one confirming iteration and stop because it lowers nothing; the rest + // are stationary at once. Never the cap, and never spinning. + [Fact] + public void NearSingularFreeSetsConfirmOnceAndStop() + { + const int m = 28, n = 24; + var confirmed = 0; + for (uint seed = 1; seed <= 60; seed++) + { + var (a, b) = ConditionedProblem(m, n, 3e7, seed); + var x = new float[n]; + var status = Solve(m, n, a, b, Fill(-1e9f, n), Fill(1e9f, n), x, out var iterations); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed}: {status}"); + Assert.True(iterations <= 3, $"seed {seed}: {iterations} iterations"); + if (iterations > 0) confirmed++; + } + Assert.InRange(confirmed, 1, 20); + } + + // Scale-freeness was only claimed for scaling A, b and the bounds together; A alone at 1e6 with the bounds left + // alone must also converge, and as fast. + [Theory] + [InlineData(1e2)] + [InlineData(1e4)] + public void LargeAWithUnscaledBoundsConvergesInAtMostTwoIterations(double cond) + { + const int m = 28, n = 24; + for (uint seed = 1; seed <= 30; seed++) + { + var (a, b) = ConditionedProblem(m, n, cond, seed, 1e6); + var x = new float[n]; + var status = Solve(m, n, a, b, Fill(-1e3f, n), Fill(1e3f, n), x, out var iterations); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"cond {cond} seed {seed}: {status} after {iterations}"); + Assert.True(iterations <= 2, $"cond {cond} seed {seed}: {iterations} iterations"); + } + } + + // b = 0 must not borrow its tolerance from the start point: a far start's own gradient is enormous, so that + // scale judged it converged after one iteration at cost 5001. The optimum here is ~1e-9. + [Fact] + public void ZeroTargetFromAFarStartIsNotJudgedConvergedByItsOwnSteepness() + { + var a = new[] { 1f, 1f, 1f, 1.0001f }; + var b = new float[2]; + var lo = new[] { 0.5f, -1e6f }; + var hi = new[] { 1e6f, 1e6f }; + var x = new[] { 1e6f, 1e6f }; + var status = Solve(2, 2, a, b, lo, hi, x, out _); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + AssertKkt(2, 2, a, b, lo, hi, x); + Assert.True(Cost(2, 2, a, b, x) < 1e-6, $"cost {Cost(2, 2, a, b, x)}"); + Assert.True(x[0] >= 0.5f); + + // The origin inside the box: the answer is the origin, from any start. + var lo2 = new[] { -1e6f, -1e6f }; + var y = new[] { 1e6f, -1e6f }; + Assert.Equal(BoundedLeastSquaresStatus.Converged, Solve(2, 2, a, b, lo2, hi, y, out _)); + Assert.True(Cost(2, 2, a, b, y) < 1e-9, $"cost {Cost(2, 2, a, b, y)}"); + } + + [Fact] + public void WarmStartClampIsExact() + { + var inf = float.PositiveInfinity; + var cases = new (float lo, float hi, float start, float expected)[] + { + (-1f, 1f, 0.25f, 0.25f), + (-1f, 1f, 5f, 1f), + (-1f, 1f, -5f, -1f), + (-1f, 1f, float.NaN, -1f), + (-inf, inf, float.NaN, 0f), + (-inf, inf, inf, 0f), + (-inf, inf, -inf, 0f), + (0.2f, inf, inf, 0.2f), + (0.2f, inf, -inf, 0.2f), + (-inf, -2f, -inf, -2f), + (-inf, -2f, 3f, -2f), + (-inf, 3f, float.NaN, 0f), + (-2f, inf, float.NaN, -2f), + }; + foreach (var (l, h, start, expected) in cases) + { + var x = new[] { start, start, start }; + // No rows: every column is dependent and takes a zero step, so x is exactly the clamped start. + Solve(0, 3, default, default, Fill(l, 3), Fill(h, 3), x, out var iterations); + Assert.Equal(0, iterations); + Assert.Equal(new[] { expected, expected, expected }, x); + } + } + + // A duplicate column that is free and interior is stepped by neither the factor nor the refinement: the columns + // share one gradient, and moving both along it would overshoot. + [Fact] + public void FreeDuplicateColumnGetsNoStep() + { + const int m = 5, n = 3; + var a = (float[])A54.Clone(); + for (var r = 0; r < m; r++) a[r * n + 2] = a[r * n + 0]; + var x = new float[n]; + var status = Solve(m, n, a, B5, Fill(-100f, n), Fill(100f, n), x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0, iterations); + Assert.Equal(0f, x[2]); + var reduced = new float[m * 2]; + for (var r = 0; r < m; r++) { reduced[r * 2] = a[r * n]; reduced[r * 2 + 1] = a[r * n + 1]; } + var xr = NormalEquationSolution(m, 2, reduced, B5); + Assert.Equal(xr[0], x[0], 1e-4); + Assert.Equal(xr[1], x[1], 1e-4); + } + + // Columns 1000x apart in scale and exactly parallel: the pivot tolerance is relative to the column's own + // diagonal, so it catches the duplicate at any absolute scale. + [Theory] + [InlineData(1e3f)] + [InlineData(1e-3f)] + public void ParallelColumnsAreDependentAtAnyAbsoluteScale(float scale) + { + const int m = 5, n = 3; + var a = (float[])A54.Clone(); + for (var r = 0; r < m; r++) a[r * n + 2] = 0.37f * a[r * n + 0]; + for (var i = 0; i < a.Length; i++) a[i] *= scale; + var lo = Fill(-1000f / scale, n); + var hi = Fill(1000f / scale, n); + var x = new float[n]; + var status = Solve(m, n, a, B5, lo, hi, x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.True(iterations <= 2, $"{iterations} iterations"); + var reduced = new float[m * 2]; + for (var r = 0; r < m; r++) { reduced[r * 2] = a[r * n]; reduced[r * 2 + 1] = a[r * n + 1]; } + var xr = NormalEquationSolution(m, 2, reduced, B5); + var xs = new float[n]; + xs[0] = (float)xr[0]; + xs[1] = (float)xr[1]; + var bb = 0.0; + foreach (var v in B5) bb += (double)v * v; + Assert.True(Cost(m, n, a, B5, x) - Cost(m, n, a, B5, xs) <= 1e-6 * bb, "cost above the reduced optimum"); + } + + // No path allocates: not the first InvalidInput, not a cap hit, not a solve that releases bounds. Each path is + // measured on its first call and again after, so a one-time runtime cost shows up as a difference between the two + // rather than as a solver allocation. + [Fact] + public void EveryPathAllocatesNothing() + { + const int m = 4, n = 3; + var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f }; + var b = new[] { 10f, -10f, 10f, 1f }; + var lo = Fill(-1f, n); + var hi = Fill(1f, n); + var workspace = new float[BoundedLeastSquares.WorkspaceLength(n)]; + var x = new float[n]; + var bad = (float[])a.Clone(); + bad[4] = float.NaN; + + var paths = new (string name, float[] a, int cap, BoundedLeastSquaresStatus expected)[] + { + ("invalid", bad, 100, BoundedLeastSquaresStatus.InvalidInput), + ("converged", a, 100, BoundedLeastSquaresStatus.Converged), + ("limit", a, 1, BoundedLeastSquaresStatus.IterationLimit), + }; + foreach (var (name, matrix, cap, expected) in paths) + { + var bytes = new long[2]; + for (var call = 0; call < 2; call++) + { + Array.Clear(x); + var before = GC.GetAllocatedBytesForCurrentThread(); + var status = BoundedLeastSquares.Solve(m, n, matrix, b, lo, hi, x, workspace, out _, cap); + bytes[call] = GC.GetAllocatedBytesForCurrentThread() - before; + Assert.Equal(expected, status); + } + Assert.True(bytes[0] == 0 && bytes[1] == 0, $"{name}: first call {bytes[0]} B, second {bytes[1]} B"); + } + } } From 4072fce0e764be0b3e9fd466f4e17eb103a07167 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 09:29:06 +0200 Subject: [PATCH 13/19] cultmath: BoundedLeastSquares drops the refinement step it cannot observe Stryker left the refinement step (removed, negated, or subtracted) unobserved, and it is unobservable for a reason: refinement with a same-precision residual does not improve forward error, and Cholesky is backward stable, so the residual the KKT test reads is already at rounding level. Measured on n=24 problems at cond(A) 1e3..1e6, the solver with and without the step gives the same iteration counts and the same x error against a Gram-Schmidt QR reference (both limited by cond(A^T A) * eps). The refinement pass, its dp buffer and the scratch vector are deleted; WorkspaceLength returns to 4n^2 + 15n. The doc states where the normal equations stop being optimal (cost gap below 1e-6 |b|^2 through cond(A) 3e5, about 2e-5 at 1e6, independent of the scale of A) and that a ridge row beyond that is the caller's policy. status is now initialised: the unassigned local made Stryker discard every mutant inside Solve as a compile error, so 255 mutants were never measured. Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 52 +++++++------------ .../BoundedLeastSquaresTests.cs | 4 +- 2 files changed, 22 insertions(+), 34 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index dfda5516..e3a4454b 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -28,7 +28,7 @@ public enum BoundedLeastSquaresStatus /// Precision. A, b and x 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. x is converted in at the start /// and back out at the end, so a returned float x is the nearest float to the double solution and always -/// inside the box. Each free-set solve gets one step of iterative refinement. +/// inside the box. /// /// /// Stopping. KKT termination uses a tolerance of 1e-6 times the problem's gradient scale (the larger of the @@ -39,11 +39,12 @@ public enum BoundedLeastSquaresStatus /// tolerance costs one confirming iteration rather than the iteration cap. /// /// -/// 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 -/// A and b; 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. +/// Near-singular problems. The normal equations square the condition number, so the cost gap grows with cond(A) +/// (measured: below 1e-6 of |b|^2 through cond(A) 3e5, about 2e-5 at 1e6), independent of the absolute scale of A. +/// The solver owns optimality up to that range. 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 beyond that range, or that need the minimum-norm answer, add a small ridge row to +/// A and b; that stays the caller's policy. /// /// The solver allocates nothing and is deterministic. /// @@ -66,7 +67,7 @@ public static class BoundedLeastSquares private const double StallRelativeTolerance = 1e-15; /// Floats of caller-supplied workspace needed for columns. - public static int WorkspaceLength(int n) => 2 * (2 * n * n + 7 * n) + 3 * n; + public static int WorkspaceLength(int n) => 2 * (2 * n * n + 6 * n) + 3 * n; /// /// Minimises ||A x - b||^2 subject to lo <= x <= hi. is the warm start @@ -95,7 +96,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 + 7 * n; + var doubleCount = 2 * n * n + 6 * n; var dbl = MemoryMarshal.Cast(workspace.Slice(0, 2 * doubleCount)); var ata = dbl.Slice(0, n * n); var chol = dbl.Slice(n * n, n * n); @@ -104,9 +105,8 @@ public static class BoundedLeastSquares var g = dbl.Slice(off + n, n); var p = dbl.Slice(off + 2 * n, n); 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); + var rhs = dbl.Slice(off + 4 * n, n); + var xd = dbl.Slice(off + 5 * n, n); // One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow. var ints = MemoryMarshal.Cast(workspace.Slice(2 * doubleCount, 3 * n)); var state = ints.Slice(0, n); @@ -166,10 +166,10 @@ public static class BoundedLeastSquares var gScale = 0.0; for (var j = 0; j < n; j++) { - q[j] = Math.Clamp(0.0, lo[j], hi[j]); + y[j] = Math.Clamp(0.0, lo[j], hi[j]); gScale = Math.Max(gScale, Math.Abs(atb[j])); } - Gradient(ata, atb, q, g, n); + Gradient(ata, atb, y, g, n); for (var j = 0; j < n; j++) gScale = Math.Max(gScale, Math.Abs(g[j])); var kktTol = KktRelativeTolerance * gScale; @@ -178,7 +178,7 @@ public static class BoundedLeastSquares // 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; + var status = BoundedLeastSquaresStatus.IterationLimit; while (true) { Gradient(ata, atb, xd, g, n); @@ -187,7 +187,7 @@ public static class BoundedLeastSquares for (var j = 0; j < n; j++) if (state[j] == 0) freeIndex[freeCount++] = j; - SolveFreeSet(ata, chol, g, p, y, r1, q, freeIndex.Slice(0, freeCount), dependent, n); + SolveFreeSet(ata, chol, g, p, y, rhs, freeIndex.Slice(0, freeCount), dependent, n); var t = 1.0; var block = -1; @@ -284,12 +284,13 @@ public static class BoundedLeastSquares } } - // Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky, then one step of iterative - // refinement with the same factor. A free column whose pivot collapses (dependent on earlier free columns) - // is dropped from the factor and gets p = 0. + // Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky. A free column whose pivot collapses + // (dependent on earlier free columns) is dropped from the factor and gets p = 0. No refinement step: Cholesky is + // backward stable, so the residual, which is the gradient the KKT test reads, is already at rounding level, and + // refining in the same precision does not improve the forward error. private static void SolveFreeSet( ReadOnlySpan ata, Span l, ReadOnlySpan g, Span p, Span y, - Span rhs, Span dp, ReadOnlySpan free, Span dependent, int n) + Span rhs, ReadOnlySpan free, Span dependent, int n) { var k = free.Length; for (var i = 0; i < k; i++) @@ -321,19 +322,6 @@ public static class BoundedLeastSquares for (var i = 0; i < k; i++) rhs[i] = dependent[i] != 0 ? 0.0 : -g[free[i]]; Substitute(l, rhs, y, p, free, n); - - for (var i = 0; i < k; i++) - { - if (dependent[i] != 0) { rhs[i] = 0.0; continue; } - var s = -g[free[i]]; - for (var c = 0; c < k; c++) - if (dependent[c] == 0) - s -= ata[free[i] * n + free[c]] * p[free[c]]; - rhs[i] = s; - } - Substitute(l, rhs, y, dp, free, n); - for (var i = 0; i < k; i++) - p[free[i]] += dp[free[i]]; } // Solves L L^T out = rhs; rhs is indexed by free-list position, out by column. diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 03177d80..f310bd9d 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -856,8 +856,8 @@ public sealed class BoundedLeastSquaresTests } } - // A duplicate column that is free and interior is stepped by neither the factor nor the refinement: the columns - // share one gradient, and moving both along it would overshoot. + // A duplicate column that is free and interior gets no step: the columns share one gradient, and moving both + // along it would overshoot. [Fact] public void FreeDuplicateColumnGetsNoStep() { From d808941b87c5dda1ce48e5845ac2b4207c2dcd89 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 09:44:04 +0200 Subject: [PATCH 14/19] cultmath: BoundedLeastSquares sheds redundant checks and scratch arithmetic Measured with Stryker once Solve was actually mutated (21 survivors), most of them were code that decided nothing: - The finite-input loops over A and b duplicated the non-finite test on the accumulated A^T b, which is exactly equivalent (floats cannot overflow a double product, and NaN/infinity never cancel to a finite value). One check on each A^T b entry remains. - The tolerance scale is the gradient at the box point nearest the origin, computed inline; the max with |A^T b| duplicated it whenever the origin is in the box, and the scratch vector for it is gone. - Bound violation is state * gradient, which is zero for a free column; the three-way ternary duplicated the free-stationarity check. - The warm-start clamp is Math.Clamp plus a NaN test. - Workspace slices come off a cursor (Take) instead of hand-summed offsets. Co-Authored-By: Claude Sonnet 5 --- .../src/CultMath/BoundedLeastSquares.cs | 81 ++++++++++--------- 1 file changed, 41 insertions(+), 40 deletions(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index e3a4454b..f03fea3b 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -31,9 +31,9 @@ public enum BoundedLeastSquaresStatus /// inside the box. /// /// -/// 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 +/// 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. 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. @@ -98,80 +98,73 @@ public static class BoundedLeastSquares // 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 dbl = MemoryMarshal.Cast(workspace.Slice(0, 2 * doubleCount)); - var ata = dbl.Slice(0, n * n); - var chol = dbl.Slice(n * n, n * n); - var off = 2 * n * n; - var atb = dbl.Slice(off, n); - var g = dbl.Slice(off + n, n); - var p = dbl.Slice(off + 2 * n, n); - var y = dbl.Slice(off + 3 * n, n); - var rhs = dbl.Slice(off + 4 * n, n); - var xd = dbl.Slice(off + 5 * n, n); + var ata = Take(ref dbl, n * n); + var chol = Take(ref dbl, n * n); + var atb = Take(ref dbl, n); + var g = Take(ref dbl, n); + var p = Take(ref dbl, n); + var y = Take(ref dbl, n); + var rhs = Take(ref dbl, n); + var xd = Take(ref dbl, n); // One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow. var ints = MemoryMarshal.Cast(workspace.Slice(2 * doubleCount, 3 * n)); - var state = ints.Slice(0, n); - var freeIndex = ints.Slice(n, n); - var dependent = ints.Slice(2 * n, n); + var state = Take(ref ints, n); + var freeIndex = Take(ref ints, n); + var dependent = Take(ref ints, n); iterations = 0; for (var j = 0; j < n; j++) if (!(lo[j] <= hi[j]) || float.IsPositiveInfinity(lo[j]) || float.IsNegativeInfinity(hi[j])) return BoundedLeastSquaresStatus.InvalidInput; - for (var i = 0; i < m * n; i++) - if (!float.IsFinite(a[i])) - return BoundedLeastSquaresStatus.InvalidInput; - for (var i = 0; i < m; i++) - if (!float.IsFinite(b[i])) - return BoundedLeastSquaresStatus.InvalidInput; + // A float is never large enough to overflow a double product or a sum of tens of them, so a non-finite + // A^T b entry is exactly a non-finite entry of A or b (NaN and infinity never cancel to a finite value). 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]; + if (!double.IsFinite(sb)) + return BoundedLeastSquaresStatus.InvalidInput; atb[j] = sb; for (var k = 0; k <= j; k++) { var sa = 0.0; for (var r = 0; r < m; r++) sa += (double)a[r * n + j] * a[r * n + k]; - if (!double.IsFinite(sa) || !double.IsFinite(sb)) - return BoundedLeastSquaresStatus.InvalidInput; ata[j * n + k] = sa; ata[k * n + j] = sa; } } - // 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. + // From here x is written. A warm start that is NaN is moved to the lower bound, and an infinite result + // (an 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]; var h = hi[j]; - var v = x[j]; - v = !(v > l) ? l : (v < h ? v : h); + var v = float.IsNaN(x[j]) ? l : Math.Clamp(x[j], l, h); if (float.IsInfinity(v)) - v = 0f < l ? l : (0f > h ? h : 0f); + v = Math.Clamp(0f, l, h); x[j] = v; xd[j] = v; state[j] = v == l ? -1 : (v == h ? 1 : 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. + // The tolerance's scale is the gradient magnitude at the box point nearest the origin (the origin itself + // when it is inside the box, where the gradient is -A^T b). 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++) { - y[j] = Math.Clamp(0.0, lo[j], hi[j]); - gScale = Math.Max(gScale, Math.Abs(atb[j])); + var s = -atb[j]; + for (var k = 0; k < n; k++) + s += ata[j * n + k] * Math.Clamp(0.0, lo[k], hi[k]); + gScale = Math.Max(gScale, Math.Abs(s)); } - Gradient(ata, atb, y, 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 @@ -236,7 +229,8 @@ public static class BoundedLeastSquares for (var j = 0; j < n; j++) { if (lo[j] == hi[j]) continue; - var violation = state[j] < 0 ? -g[j] : (state[j] > 0 ? g[j] : 0.0); + // 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) { worst = violation; @@ -259,6 +253,13 @@ public static class BoundedLeastSquares return status; } + private static Span Take(ref Span pool, int length) + { + var head = pool.Slice(0, length); + pool = pool.Slice(length); + return head; + } + // 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 x, ReadOnlySpan g, ReadOnlySpan atb, int n, out double scale) From 0a806249307ab932ba667e0ec4ab387893b364f0 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 09:44:04 +0200 Subject: [PATCH 15/19] cultmath: BoundedLeastSquares tests observe the warm-start clamp and the workspace The full-step clamp repairs an out-of-box start on a free column, so the warm-start clamp was invisible to the old test (Stryker: five survivors). The new table runs one column with b beyond the box and a cap of 0, where the solver returns the clamped start untouched, across finite, infinite and NaN starts and finite and infinite bounds. Workspace junk (including NaN doubles) is now checked on random boxed problems that release bounds and use the tolerance, with bit-identical x and equal iteration counts. Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 46 ++++++++++++++++--- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index f310bd9d..bed7bdae 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -412,6 +412,29 @@ public sealed class BoundedLeastSquaresTests } } + // Junk in the workspace, including NaN doubles, must not reach any decision: random boxed problems (which release + // bounds and use the tolerance) give bit-identical x and the same iteration count. + [Fact] + public void WorkspaceJunkNeverReachesTheSolveOnRandomProblems() + { + for (uint seed = 1; seed <= 25; seed++) + { + var (m, n, a, b, lo, hi) = RandomProblem(seed); + var clean = new float[n]; + Solve(m, n, a, b, lo, hi, clean, out var cleanIterations); + foreach (var junk in new[] { float.NaN, 7f, -1e30f }) + { + var workspace = new float[BoundedLeastSquares.WorkspaceLength(n)]; + Array.Fill(workspace, junk); + var x = new float[n]; + BoundedLeastSquares.Solve(m, n, a, b, lo, hi, x, workspace, out var iterations); + Assert.Equal(cleanIterations, iterations); + for (var j = 0; j < n; j++) + Assert.Equal(BitConverter.SingleToInt32Bits(clean[j]), BitConverter.SingleToInt32Bits(x[j])); + } + } + } + [Fact] public void WarmStartAtAVertexTakesNoIterations() { @@ -826,6 +849,9 @@ public sealed class BoundedLeastSquaresTests Assert.True(Cost(2, 2, a, b, y) < 1e-9, $"cost {Cost(2, 2, a, b, y)}"); } + // One column, A = [1], with b beyond the box: the first Newton step would cross a bound, so with a cap of 0 the + // solver returns the clamped start untouched. That makes the clamp itself the thing under test, for every + // combination of a finite, infinite or NaN start with finite and infinite bounds. [Fact] public void WarmStartClampIsExact() { @@ -836,23 +862,29 @@ public sealed class BoundedLeastSquaresTests (-1f, 1f, 5f, 1f), (-1f, 1f, -5f, -1f), (-1f, 1f, float.NaN, -1f), - (-inf, inf, float.NaN, 0f), - (-inf, inf, inf, 0f), - (-inf, inf, -inf, 0f), + (0.5f, 2f, inf, 2f), + (0.5f, 2f, -inf, 0.5f), (0.2f, inf, inf, 0.2f), (0.2f, inf, -inf, 0.2f), + (0.2f, inf, float.NaN, 0.2f), (-inf, -2f, -inf, -2f), (-inf, -2f, 3f, -2f), (-inf, 3f, float.NaN, 0f), + (-inf, 3f, -inf, 0f), (-2f, inf, float.NaN, -2f), + (-inf, -2f, float.NaN, -2f), + (-inf, 3f, 1f, 1f), }; foreach (var (l, h, start, expected) in cases) { - var x = new[] { start, start, start }; - // No rows: every column is dependent and takes a zero step, so x is exactly the clamped start. - Solve(0, 3, default, default, Fill(l, 3), Fill(h, 3), x, out var iterations); + // b beyond the box on the side the start is not on the bound of, so the step crosses a bound (or the + // start already sits on one). + var b = new[] { float.IsInfinity(h) ? l - 10f : h + 10f }; + if (float.IsInfinity(l) && float.IsInfinity(h)) b[0] = 0f; + var x = new[] { start }; + Solve(1, 1, new[] { 1f }, b, new[] { l }, new[] { h }, x, out var iterations, maxIterations: 0); Assert.Equal(0, iterations); - Assert.Equal(new[] { expected, expected, expected }, x); + Assert.True(x[0] == expected, $"box [{l}, {h}] start {start}: {x[0]}, expected {expected}"); } } From 854f01e4c7e6cb7be155fe370625acb9187871d0 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 10:37:27 +0200 Subject: [PATCH 16/19] cultmath: BoundedLeastSquares tests pin the tolerance scale and the zero gradient Two tolerance behaviours had no observer: the scale taken at the box point nearest the origin (a 6e-5 release on a 1000-scale problem, in either column position), and the strict comparisons when the scale and the gradient are both exactly zero (b = 0 at the origin takes no iteration). Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index bed7bdae..dd9ba3f9 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -849,6 +849,45 @@ public sealed class BoundedLeastSquaresTests Assert.True(Cost(2, 2, a, b, y) < 1e-9, $"cost {Cost(2, 2, a, b, y)}"); } + // The tolerance scale is the gradient at the box point nearest the origin. Here that point is 1000 and the + // wanted release is a gradient of about 6e-5 against a 1000-scale problem: only a scale taken from the box point + // (about 6e-5, tolerance 6e-11) releases it; a scale taken from A^T b (about 1000, tolerance 1e-3) would not. + // Run with the excluded column in each position so the scale must come from the right entry. + [Fact] + public void ToleranceScaleIsTheGradientAtTheBoxPointNearestTheOrigin() + { + var a = new[] { 1f, 0f, 0f, 1f }; + foreach (var excluded in new[] { 0, 1 }) + { + var b = new float[2]; + var lo = new[] { -1e6f, -1e6f }; + var hi = new[] { 1e6f, 1e6f }; + b[excluded] = 1000.0001f; + lo[excluded] = 1000f; + hi[excluded] = 2000f; + var x = new float[2]; + x[excluded] = 1000f; + var status = Solve(2, 2, a, b, lo, hi, x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.True(iterations >= 1, $"column {excluded}: no release"); + Assert.Equal(b[excluded], x[excluded]); + Assert.Equal(0f, x[1 - excluded]); + } + } + + // b = 0 with the origin in the box and the start at the origin: the gradient is exactly zero and so is the + // tolerance. Nothing is violated, so the solve is over without an iteration. + [Fact] + public void ZeroTargetAtTheOriginTakesNoIterations() + { + const int m = 5, n = 3; + var x = new float[n]; + var status = Solve(m, n, A54, new float[m], Fill(-1f, n), Fill(1f, n), x, out var iterations); + Assert.Equal(BoundedLeastSquaresStatus.Converged, status); + Assert.Equal(0, iterations); + Assert.Equal(new float[n], x); + } + // One column, A = [1], with b beyond the box: the first Newton step would cross a bound, so with a cap of 0 the // solver returns the clamped start untouched. That makes the clamp itself the thing under test, for every // combination of a finite, infinite or NaN start with finite and infinite bounds. From 6328f153453a8f7a7062cd653ea6c30fed75c1df Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 10:52:39 +0200 Subject: [PATCH 17/19] cultmath: BoundedLeastSquares test drives the tolerance below the noise floor With b along the smallest left singular vector, |A^T b| is tiny, so the KKT tolerance (1e-6 of it) sits below the gradient's rounding noise at cond(A) >= 1e6 and can never be met. The solver must stop on the progress rule, in at most six passes, never at the cap. This is the path the progress rule exists for; the earlier fixtures reached it only by accident. Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 27 +++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index dd9ba3f9..3888b004 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -722,7 +722,7 @@ public sealed class BoundedLeastSquaresTests // A = U diag(s) V^T with s geometric from 1 down to 1/cond, so cond(A) is the requested value (float rounding of A // perturbs it slightly). m x n with m >= n; b is generic, so it has a real residual when m > n. - private static (float[] a, float[] b) ConditionedProblem(int m, int n, double cond, uint seed, double aScale = 1.0) + private static (float[] a, float[] b) ConditionedProblem(int m, int n, double cond, uint seed, double aScale = 1.0, bool bAlongSmallestSingularVector = false) { var rng = new CultMath.Random(seed * 2654435761u + 17u); double[][] Orthonormal(int dim, int count) @@ -758,7 +758,7 @@ public sealed class BoundedLeastSquaresTests a[r * n + c] += (float)(u[k][r] * s * w[k][c]); } var b = new float[m]; - for (var i = 0; i < m; i++) b[i] = rng.NextFloat(-3f, 3f); + for (var i = 0; i < m; i++) b[i] = bAlongSmallestSingularVector ? (float)u[n - 1][i] : rng.NextFloat(-3f, 3f); return (a, b); } @@ -808,6 +808,29 @@ public sealed class BoundedLeastSquaresTests Assert.InRange(confirmed, 1, 20); } + // b along the smallest left singular vector makes the optimum lie along the weakest direction: |A^T b| is tiny, so + // the tolerance (1e-6 of it) sits below the rounding noise of the gradient (~1e-16 |A^T A| |x|) for cond(A) >= 1e6. + // The tolerance is then unreachable, and the solver must stop on the progress rule, after the handful of Newton + // passes that still lower the cost, never at the 100-iteration cap. + [Theory] + [InlineData(1e6)] + [InlineData(1e7)] + public void ToleranceBelowTheNoiseFloorStopsOnProgressNotOnTheCap(double cond) + { + const int m = 28, n = 24; + var confirmed = 0; + for (uint seed = 1; seed <= 30; seed++) + { + var (a, b) = ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: true); + var x = new float[n]; + var status = Solve(m, n, a, b, Fill(-1e9f, n), Fill(1e9f, n), x, out var iterations); + Assert.True(status == BoundedLeastSquaresStatus.Converged, $"cond {cond} seed {seed}: {status} after {iterations}"); + Assert.True(iterations <= 6, $"cond {cond} seed {seed}: {iterations} iterations"); + if (iterations >= 1) confirmed++; + } + Assert.True(confirmed >= 15, $"cond {cond}: only {confirmed} of 30 reached the tolerance floor"); + } + // Scale-freeness was only claimed for scaling A, b and the bounds together; A alone at 1e6 with the bounds left // alone must also converge, and as fast. [Theory] From c8227c6347e648c9a04fb671813178716e4fea8b Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 13:14:21 +0200 Subject: [PATCH 18/19] cultmath: BoundedLeastSquares tests drive stalled passes to the stall rule Batches of small boxed problems whose passes end with the KKT test failing and the cost unmoved (dependent free columns, releases the factor drops, weak-direction noise-floor problems). They must report Converged, and land within the documented gap of an exact active-set enumeration. Pins the stall comparison and the cost scale in Objective, which no test reached before. Co-Authored-By: Claude Sonnet 5 --- .../BoundedLeastSquaresTests.cs | 108 +++++++++++++++++- 1 file changed, 106 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 3888b004..a9a950be 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -616,8 +616,15 @@ public sealed class BoundedLeastSquaresTests } // Exact reference: every assignment of each column to {lo, hi, free}, free set solved in double. - private static double ExactBoxedCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi) + private static double ExactBoxedCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi, double pivotTol = 1e-9, double ridge = 0.0) { + var maxDiag = 0.0; + for (var j = 0; j < n; j++) + { + var d = 0.0; + for (var r = 0; r < m; r++) d += (double)a[r * n + j] * a[r * n + j]; + maxDiag = Math.Max(maxDiag, d); + } var best = double.MaxValue; var total = 1; for (var j = 0; j < n; j++) total *= 3; @@ -646,12 +653,13 @@ public sealed class BoundedLeastSquaresTests mat[i, k] += (double)a[r * n + free[i]] * (b[r] - fixedPart); } } + for (var i = 0; i < k; i++) mat[i, i] += ridge * maxDiag; var ok = true; for (var col = 0; col < k && ok; col++) { var piv = col; for (var r = col + 1; r < k; r++) if (Math.Abs(mat[r, col]) > Math.Abs(mat[piv, col])) piv = r; - if (Math.Abs(mat[piv, col]) < 1e-9) { ok = false; break; } + if (Math.Abs(mat[piv, col]) < pivotTol) { ok = false; break; } for (var l = 0; l <= k; l++) (mat[col, l], mat[piv, l]) = (mat[piv, l], mat[col, l]); for (var r = 0; r < k; r++) { @@ -762,6 +770,102 @@ public sealed class BoundedLeastSquaresTests return (a, b); } + // A pass can end with the KKT test still failing and the cost exactly where it was: a dependent free column keeps a + // gradient above the tolerance, or a bound release lands on a column the factor drops. The solver must call that + // converged; if the stall rule mis-reads the cost scale or the progress, it spins to the iteration cap instead. + // Each kind is a fixed-seed batch of small boxed problems with warm starts on and off the bounds and pinned columns: + // rank-deficient (duplicate and zero columns), nearly collinear (1e-2 to 1e-5), integer-tied, wide, cond(A) 1e5 to + // 3e5, and "weak", where b lies along the weakest singular direction of a cond(A) 1e6 to 1e7 matrix so the gradient + // tolerance sits below the rounding noise with interior columns nonzero. Every kind must converge; all but "weak" + // (past the documented condition range) must also land within the documented cost gap, relative to |b|^2, of an + // exact enumeration of the active sets (a 1e-15 ridge keeps singular free sets solvable). + [Theory] + [InlineData("plain", 300)] + [InlineData("rankdef", 300)] + [InlineData("collinear", 300)] + [InlineData("integer-tie", 300)] + [InlineData("wide", 300)] + [InlineData("cond", 300)] + [InlineData("weak", 6000)] + public void StalledPassesConvergeInsteadOfSpinningToTheCap(string kind, int seeds) + { + var failures = new System.Collections.Generic.List(); + for (uint seed = 1; seed <= seeds; seed++) + { + var rng = new CultMath.Random(seed * 2246822519u + (uint)kind.Length * 3266489917u); + double Gauss() => Math.Sqrt(-2.0 * Math.Log(1.0 - rng.NextFloat() * 0.999999)) * Math.Cos(2.0 * Math.PI * rng.NextFloat()); + var n = rng.NextInt(2, 8); + var m = rng.NextInt(2, 12); + if (kind == "wide") m = Math.Max(1, n - rng.NextInt(1, 4)); + if (kind == "cond" || kind == "weak") m = n + rng.NextInt(0, 5); + float[] a; + float[] bWeak = null; + if (kind == "weak") + { + (a, bWeak) = ConditionedProblem(m, n, Math.Pow(10.0, 6.0 + rng.NextFloat()), seed, 1.0, true); + } + else if (kind == "cond") + { + a = ConditionedProblem(m, n, Math.Pow(10.0, 5.0 + 0.5 * rng.NextFloat()), seed).a; + } + else + { + a = new float[m * n]; + for (var i = 0; i < m; i++) + { + var row = new double[n]; + for (var j = 0; j < n; j++) row[j] = kind == "integer-tie" ? rng.NextInt(-2, 3) : Gauss(); + for (var j = 1; j < n; j++) + { + if (kind == "rankdef" && j % 3 == 1) row[j] = row[j - 1]; + if (kind == "rankdef" && j % 3 == 2) row[j] = 0.0; + if (kind == "collinear" && j % 2 == 1) row[j] = row[j - 1] + Math.Pow(10.0, -rng.NextInt(2, 6)) * Gauss(); + } + for (var j = 0; j < n; j++) a[i * n + j] = (float)row[j]; + } + } + var lo = new float[n]; + var hi = new float[n]; + var x = new float[n]; + var xt = new double[n]; + for (var j = 0; j < n; j++) + { + var w = Math.Pow(10.0, rng.NextInt(-2, 2)); + lo[j] = (float)(-w * rng.NextFloat()); + hi[j] = (float)(w * rng.NextFloat()); + if (rng.NextInt(6) == 0) hi[j] = lo[j]; + else if (rng.NextInt(8) == 0) lo[j] = hi[j] = 0f; + else if (rng.NextInt(8) == 0) { lo[j] = (float)(0.5 + rng.NextFloat()); hi[j] = lo[j] + 0.5f; } + if (kind == "weak") { var c = (float)Math.Pow(10.0, rng.NextInt(3, 9)); lo[j] = rng.NextInt(2) == 0 ? -c : 0f; hi[j] = c; } + if (kind == "integer-tie") { lo[j] = -1f; hi[j] = rng.NextInt(4) == 0 ? -1f : 1f; } + x[j] = rng.NextInt(3) == 0 ? lo[j] : rng.NextInt(2) == 0 ? hi[j] : lo[j] + (hi[j] - lo[j]) * rng.NextFloat(); + xt[j] = (lo[j] + hi[j]) / 2.0 + (hi[j] - lo[j]) * 3.0 * (rng.NextFloat() - 0.5); + } + var noise = Math.Pow(10.0, rng.NextInt(-2, 4)); + var mode = rng.NextInt(4); + var b = new float[m]; + for (var i = 0; i < m; i++) + { + var s = 0.0; + for (var j = 0; j < n; j++) s += (double)a[i * n + j] * xt[j]; + b[i] = bWeak != null ? bWeak[i] : mode == 0 ? 0f : (float)(s + noise * Gauss() * (mode == 3 ? 1.0 : 0.01)); + } + var status = Solve(m, n, a, b, lo, hi, x, out var iterations); + if (status != BoundedLeastSquaresStatus.Converged) + { + failures.Add($"seed {seed}: {status} after {iterations}"); + continue; + } + if (kind == "weak") continue; + var bb = 0.0; + foreach (var v in b) bb += (double)v * v; + if (bb == 0.0) bb = 1.0; + var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi, 1e-300, 1e-15); + if (gap > 1e-5 * bb) failures.Add($"seed {seed}: cost gap {gap / bb:E2} of |b|^2"); + } + Assert.True(failures.Count == 0, $"{kind}: {failures.Count} of {seeds} failed; first: {string.Join("; ", failures.GetRange(0, Math.Min(3, failures.Count)))}"); + } + // Unsaturated allocator-shaped ticks: 24 columns, every column free at the optimum. The solve is one Newton step; // at most one more may confirm it. Before the working iterate moved to double, the float gradient's noise floor // exceeded the tolerance from cond ~1e3 and the solver spun to the 100-iteration cap on optimal points. From 09507901d00a8fe4a35bc6cabcb1fb58b5b08045 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 13:14:21 +0200 Subject: [PATCH 19/19] cultmath: BoundedLeastSquares doc states the cost gap is relative to |b|^2 Co-Authored-By: Claude Sonnet 5 --- packages/cultmath/src/CultMath/BoundedLeastSquares.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index f03fea3b..4cfbad74 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -40,7 +40,7 @@ public enum BoundedLeastSquaresStatus /// /// /// Near-singular problems. The normal equations square the condition number, so the cost gap grows with cond(A) -/// (measured: below 1e-6 of |b|^2 through cond(A) 3e5, about 2e-5 at 1e6), independent of the absolute scale of A. +/// (measured: below 1e-6 of |b|^2 through cond(A) 3e5, about 2e-5 at 1e6; the gap is relative to |b|^2), independent of the absolute scale of A. /// The solver owns optimality up to that range. 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 beyond that range, or that need the minimum-norm answer, add a small ridge row to