diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs
new file mode 100644
index 00000000..4cfbad74
--- /dev/null
+++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs
@@ -0,0 +1,347 @@
+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,
+
+ /// The input broke the contract; nothing was solved and x is exactly as passed.
+ InvalidInput,
+}
+
+///
+/// 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.
+///
+/// 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.
+///
+///
+/// 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.
+///
+///
+/// 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; 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
+/// A and b; that stays the caller's policy.
+///
+/// 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
+{
+ /// Default iteration cap for .
+ public const int DefaultMaxIterations = 100;
+
+ 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;
+
+ ///
+ /// 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; -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.
+ /// 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.");
+
+ // 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 = 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 = 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;
+
+ // 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];
+ ata[j * n + k] = sa;
+ ata[k * n + j] = sa;
+ }
+ }
+
+ // 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 = float.IsNaN(x[j]) ? l : Math.Clamp(x[j], l, h);
+ if (float.IsInfinity(v))
+ 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 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++)
+ {
+ 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));
+ }
+ 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;
+ var status = BoundedLeastSquaresStatus.IterationLimit;
+ while (true)
+ {
+ Gradient(ata, atb, xd, 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, rhs, freeIndex.Slice(0, freeCount), dependent, n);
+
+ var t = 1.0;
+ var block = -1;
+ for (var i = 0; i < freeCount; i++)
+ {
+ var j = freeIndex[i];
+ var pj = p[j];
+ if (pj == 0.0) continue;
+ var tj = pj > 0.0 ? (hi[j] - xd[j]) / pj : (lo[j] - xd[j]) / pj;
+ if (tj < t)
+ {
+ t = tj;
+ block = j;
+ }
+ }
+
+ if (block >= 0)
+ {
+ if (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; }
+ iterations++;
+ for (var i = 0; i < freeCount; i++)
+ {
+ var j = freeIndex[i];
+ xd[j] = Math.Clamp(xd[j] + t * p[j], lo[j], hi[j]);
+ }
+ 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];
+ xd[j] = Math.Clamp(xd[j] + p[j], lo[j], hi[j]);
+ }
+ Gradient(ata, atb, xd, g, n);
+
+ // 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 (Math.Abs(g[freeIndex[i]]) > kktTol)
+ freeViolated = true;
+
+ var release = -1;
+ var worst = kktTol;
+ for (var j = 0; j < n; j++)
+ {
+ if (lo[j] == hi[j]) continue;
+ // A bound at lo wants a positive gradient and one at hi a negative one; a free column has state 0.
+ var violation = state[j] * g[j];
+ if (violation > worst)
+ {
+ worst = violation;
+ release = j;
+ }
+ }
+ 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 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)
+ {
+ 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++)
+ {
+ 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. 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, 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])
+ {
+ dependent[i] = 1;
+ 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 = Math.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++)
+ rhs[i] = dependent[i] != 0 ? 0.0 : -g[free[i]];
+ Substitute(l, rhs, y, p, free, n);
+ }
+
+ // 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];
+ }
+ 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] * outByColumn[free[r]];
+ outByColumn[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..a9a950be
--- /dev/null
+++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs
@@ -0,0 +1,1141 @@
+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]}");
+ }
+ }
+
+ // 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 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++) 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(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++) 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,
+ 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 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,
+ 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);
+
+ var expected = NormalEquationSolution(m, n, A54, B5);
+ for (var i = 0; i < n; i++) Assert.Equal(expected[i], x[i], 1e-4);
+ }
+
+ [Fact]
+ public void KktHoldsOnRandomBoxedProblems()
+ {
+ var anyActive = false;
+ for (uint seed = 1; seed <= 150; seed++)
+ {
+ 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");
+ 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 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);
+ }
+
+ [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 FixedColumnStaysPutEvenWhenItsGradientWantsOut()
+ {
+ const int m = 5, n = 3;
+ 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);
+ 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);
+ }
+
+ [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 LargeInteriorOptimumEqualsNormalEquations()
+ {
+ 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);
+ }
+ }
+ }
+
+ // 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()
+ {
+ 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 _));
+ }
+
+ // 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, 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;
+ 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);
+ }
+ }
+ 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]) < 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++)
+ {
+ 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;
+ }
+
+ // m < n with columns that are nearly parallel (condition number of A around 1e2), b outside the reachable set.
+ // 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++)
+ {
+ for (uint seed = 1; seed <= 120; seed++)
+ {
+ 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;
+ }
+ }
+ 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];
+ 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})");
+ }
+ }
+ }
+
+ // 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, bool bAlongSmallestSingularVector = false)
+ {
+ 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] = bAlongSmallestSingularVector ? (float)u[n - 1][i] : rng.NextFloat(-3f, 3f);
+ 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.
+ [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);
+ }
+
+ // 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]
+ [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)}");
+ }
+
+ // 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.
+ [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),
+ (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)
+ {
+ // 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.True(x[0] == expected, $"box [{l}, {h}] start {start}: {x[0]}, expected {expected}");
+ }
+ }
+
+ // 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()
+ {
+ 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");
+ }
+ }
+}