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;