From 44affbb9b65125ad6005daed50c2a8f2efd9241d Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Tue, 29 Sep 2026 03:46:36 +0200 Subject: [PATCH] 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]; } } }