From 08f05dede8001d084bc0400eca329e7cf7433887 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 10 Oct 2026 22:24:02 +0200 Subject: [PATCH] CultMath BoundedLeastSquares: the stopping rule's decisions are traced and pinned; Validate is the one owner of n < 0 The release choice and the pass judgement are internal functions (SelectRelease, JudgePass) and Solve reports each decision, in order and with the values it compared, to an optional internal ISolveTrace. The tests assert those decisions with exact equality against the traces the base commit dbcb2a8e produced on the same problems, and call the two functions directly on states a solve may never reach. Solve's own n < 0 check is deleted (Validate throws the same exception), and the NonFiniteTarget member doc names the n = 0 case. Co-Authored-By: Claude Opus 5.5 --- .../src/CultMath/BoundedLeastSquares.cs | 136 ++++- .../BoundedLeastSquaresBaseTraces.cs | 477 ++++++++++++++++++ .../BoundedLeastSquaresTests.cs | 89 +--- .../BoundedLeastSquaresTrace.cs | 223 ++++++++ .../BoundedLeastSquaresTraceTests.cs | 226 +++++++++ 5 files changed, 1044 insertions(+), 107 deletions(-) create mode 100644 packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresBaseTraces.cs create mode 100644 packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTrace.cs create mode 100644 packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTraceTests.cs diff --git a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs index c005c215..b19af384 100644 --- a/packages/cultmath/src/CultMath/BoundedLeastSquares.cs +++ b/packages/cultmath/src/CultMath/BoundedLeastSquares.cs @@ -30,10 +30,41 @@ public enum BoundedLeastSquaresInput /// Some entry of a is NaN or infinite. NonFiniteMatrix, - /// Some entry of b is NaN or infinite. + /// Some entry of b is NaN or infinite. With n = 0 there is no unknown, b is never read and no entry is judged. NonFiniteTarget, } +/// A decision took, for tests that pin the stopping rule by its decisions. +internal enum SolveEvent +{ + /// A step stopped at a bound: column, step length. + Step, + + /// A bound whose violation beat the best so far, judged against its rounding bound: column, violation, noise, best before. + Candidate, + + /// A bound released: column, the stall mark after the release. + Release, + + /// KKT holds and the solve is over: the tolerance. + Kkt, + + /// A stationarity miss with no release, the cost fell and the solve goes on: cost, mark before, allowance. + Miss, + + /// A stationarity miss with no release, the cost did not fall: the solve is over: cost, mark, allowance. + Stall, + + /// The iteration cap ended the solve. + Cap, +} + +/// Observes the decisions of one solve in order; never changes one. +internal interface ISolveTrace +{ + void Record(SolveEvent kind, int column, double a, double b, double c); +} + /// /// 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. @@ -182,8 +213,20 @@ public static class BoundedLeastSquares out int iterations, double kktRelativeTolerance, int maxIterations = DefaultMaxIterations) + => SolveTraced(m, n, a, b, lo, hi, x, workspace, out iterations, kktRelativeTolerance, maxIterations, null); + + // The whole solve. A trace sink only observes: every decision below is taken from the same values with or without one. + internal static BoundedLeastSquaresStatus SolveTraced( + int m, int n, + ReadOnlySpan a, ReadOnlySpan b, + ReadOnlySpan lo, ReadOnlySpan hi, + Span x, Span workspace, + out int iterations, + double kktRelativeTolerance, + int maxIterations, + ISolveTrace? trace) { - if (n < 0 || x.Length < n || workspace.Length < WorkspaceLength(n)) + if (x.Length < n || workspace.Length < WorkspaceLength(n)) throw new ArgumentException("BoundedLeastSquares: span shorter than the stated dimensions."); iterations = 0; if (Validate(m, n, a, b, lo, hi, kktRelativeTolerance) != BoundedLeastSquaresInput.Valid) @@ -284,8 +327,9 @@ public static class BoundedLeastSquares if (block >= 0) { - if (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; } + if (iterations >= maxIterations) { trace?.Record(SolveEvent.Cap, -1, 0.0, 0.0, 0.0); status = BoundedLeastSquaresStatus.IterationLimit; break; } iterations++; + trace?.Record(SolveEvent.Step, block, t, 0.0, 0.0); for (var i = 0; i < freeCount; i++) { var j = freeIndex[i]; @@ -309,30 +353,9 @@ public static class BoundedLeastSquares 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 && violation > ReleaseNoise(ata, atb, xd, j, n)) - { - worst = violation; - release = j; - } - } - if (release < 0 && !freeViolated) { status = BoundedLeastSquaresStatus.Converged; break; } - - if (release < 0) - { - var f = Objective(xd, g, atb, n, out var fScale); - if (f >= fMark - StallRelativeTolerance * fScale) { status = BoundedLeastSquaresStatus.Converged; break; } - fMark = f; - } - else - fMark = double.PositiveInfinity; - if (iterations >= maxIterations) { status = BoundedLeastSquaresStatus.IterationLimit; break; } + var release = SelectRelease(ata, atb, xd, g, state, lo, hi, kktTol, n, trace); + if (JudgePass(release, freeViolated, kktTol, ref fMark, xd, g, atb, n, trace)) { status = BoundedLeastSquaresStatus.Converged; break; } + if (iterations >= maxIterations) { trace?.Record(SolveEvent.Cap, -1, 0.0, 0.0, 0.0); status = BoundedLeastSquaresStatus.IterationLimit; break; } iterations++; // A stationarity miss re-solves the same free set from the improved point. if (release >= 0) state[release] = 0; @@ -343,6 +366,63 @@ public static class BoundedLeastSquares return status; } + // The bound to release: the one with the most violated gradient sign, the first of equal violations, and only if + // the violation exceeds both the tolerance (the starting `worst`) and the rounding bound of its own gradient. + // -1 when none qualifies. + internal static int SelectRelease( + ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, ReadOnlySpan g, + ReadOnlySpan state, ReadOnlySpan lo, ReadOnlySpan hi, double kktTol, int n, ISolveTrace? trace) + { + 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) + { + var noise = ReleaseNoise(ata, atb, x, j, n); + trace?.Record(SolveEvent.Candidate, j, violation, noise, worst); + if (violation > noise) + { + worst = violation; + release = j; + } + } + } + return release; + } + + // What a full step's outcome means for the solve: true when it is over. Nothing to release and a stationary free + // set is KKT. Nothing to release with a stationarity miss ends when the cost has not fallen below the last such + // miss's by more than rounding, else it records the cost. A release is never judged by the cost and resets the mark. + internal static bool JudgePass( + int release, bool freeViolated, double kktTol, ref double fMark, + ReadOnlySpan x, ReadOnlySpan g, ReadOnlySpan atb, int n, ISolveTrace? trace) + { + if (release < 0 && !freeViolated) + { + trace?.Record(SolveEvent.Kkt, -1, kktTol, 0.0, 0.0); + return true; + } + if (release < 0) + { + var f = Objective(x, g, atb, n, out var fScale); + var allowance = StallRelativeTolerance * fScale; + var stalled = f >= fMark - allowance; + trace?.Record(stalled ? SolveEvent.Stall : SolveEvent.Miss, -1, f, fMark, allowance); + if (stalled) return true; + fMark = f; + } + else + { + fMark = double.PositiveInfinity; + trace?.Record(SolveEvent.Release, release, fMark, 0.0, 0.0); + } + return false; + } + private static Span Take(ref Span pool, int length) { var head = pool.Slice(0, length); @@ -365,7 +445,7 @@ public static class BoundedLeastSquares } // Rounding bound of the double gradient component j: gamma_(n+1) (Higham, inner products) times the magnitude of its terms. - private static double ReleaseNoise(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, int j, int n) + internal static double ReleaseNoise(ReadOnlySpan ata, ReadOnlySpan atb, ReadOnlySpan x, int j, int n) { var s = Math.Abs(atb[j]); for (var k = 0; k < n; k++) s += Math.Abs(ata[j * n + k] * x[k]); diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresBaseTraces.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresBaseTraces.cs new file mode 100644 index 00000000..33243d78 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresBaseTraces.cs @@ -0,0 +1,477 @@ +using System.Collections.Generic; + +namespace CultMath.Tests; + +// The decisions of the base commit dbcb2a8e (BoundedLeastSquares.cs with the same events added at the same points, nothing else changed) +// on the BlsTrace cases. Head must reproduce every line. +internal static class BlsBaseTraces +{ + internal static readonly Dictionary Expected = new() + { + ["tie-pair"] = """ + Candidate 0 1 6.661338147750939E-16 1E-06 + Release 0 Infinity 0 0 + Kkt -1 1E-06 0 0 + end Converged 1 3f800000 00000000 + """, + ["tie-triple"] = """ + Candidate 0 2 1.7763568394002505E-15 2E-06 + Release 0 Infinity 0 0 + Kkt -1 2E-06 0 0 + end Converged 1 3f800000 00000000 00000000 + """, + ["tie-pair-at-hi"] = """ + Candidate 0 3 1.9984014443252818E-15 1E-06 + Release 0 Infinity 0 0 + Step 0 0.3333333333333333 0 0 + Candidate 1 2 1.3322676295501878E-15 1E-06 + Release 1 Infinity 0 0 + Step 1 0.5 0 0 + Kkt -1 1E-06 0 0 + end Converged 4 00000000 00000000 + """, + ["step-lands-on-bound"] = """ + Kkt -1 1E-06 0 0 + end Converged 0 3f800000 + """, + ["step-blocked"] = """ + Step 0 0.2 0 0 + Kkt -1 3E-06 0 0 + end Converged 1 3f800000 bf000000 + """, + ["w5-shift0-1e-12"] = """ + Candidate 1 135111.75592265668 3.900108740394563E-10 1.5002178145105112E-07 + Candidate 7 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 7 Infinity 0 0 + Step 7 0.4366038070185627 0 0 + Candidate 1 76194.22899199842 5.600810209156691E-10 1.5002178145105112E-07 + Release 1 Infinity 0 0 + Step 1 0.7573476181786473 0 0 + Candidate 2 11849.84792945953 4.782591048738879E-10 1.5002178145105112E-07 + Candidate 4 17935.626996204504 7.122810428200547E-10 11849.84792945953 + Release 4 Infinity 0 0 + Candidate 0 633.4374486887755 6.019496088675069E-10 1.5002178145105112E-07 + Candidate 3 776.2927740423111 8.68340538093788E-10 633.4374486887755 + Candidate 11 2665.4946252912787 3.120333456889181E-11 776.2927740423111 + Release 11 Infinity 0 0 + Candidate 0 0.2169954634180158 6.022639451081885E-10 1.5002178145105112E-07 + Candidate 3 0.26539156234684924 8.683991289001988E-10 0.2169954634180158 + Release 3 Infinity 0 0 + Step 4 0.06907420089782461 0 0 + Candidate 0 0.0362638246792244 7.758593413765183E-10 1.5002178145105112E-07 + Candidate 2 0.0760455647832714 6.59955603932088E-10 0.0362638246792244 + Candidate 5 0.1938752159811088 5.480900325737566E-11 0.0760455647832714 + Release 5 Infinity 0 0 + Step 11 0.03215405619865555 0 0 + Candidate 0 0.004646654699172359 7.758594326420045E-10 1.5002178145105112E-07 + Candidate 1 0.008483701902150642 9.999787174674693E-10 0.004646654699172359 + Candidate 4 0.01007973030937137 9.812607064514548E-10 0.008483701902150642 + Release 4 Infinity 0 0 + Candidate 0 0.0032877034464036115 6.980323623931817E-10 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Candidate 2 0.002053476287983358 6.004885051312699E-10 1.5002178145105112E-07 + Release 2 Infinity 0 0 + Candidate 7 0.0014518519601551816 9.294667933498287E-10 1.5002178145105112E-07 + Release 7 Infinity 0 0 + Kkt -1 1.5002178145105112E-07 0 0 + end Converged 14 3e64b5b7 3f800000 3477cb66 3e0efa90 3f48c5cd 3f44d339 00000000 3f6c163f 00000000 00000000 00000000 00000000 + """, + ["w5-shift1-1e-12"] = """ + Candidate 0 6737.644464075565 1.9448748840962782E-11 1.5002178145105112E-07 + Candidate 2 135111.75592265668 3.900108740394563E-10 6737.644464075565 + Candidate 8 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 8 Infinity 0 0 + Step 8 0.4366038070185627 0 0 + Candidate 0 5151.85545116663 2.402625547421844E-11 1.5002178145105112E-07 + Candidate 2 76194.22899199842 5.600810209156691E-10 5151.85545116663 + Release 2 Infinity 0 0 + Step 2 0.7573476181786473 0 0 + Candidate 0 3065.16054302454 3.0049666978425744E-11 1.5002178145105112E-07 + Candidate 3 11849.84792945953 4.782591048738879E-10 3065.16054302454 + Candidate 5 17935.626996204504 7.122810428200547E-10 11849.84792945953 + Release 5 Infinity 0 0 + Candidate 0 2665.4946252912787 3.120333456889181E-11 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Candidate 1 0.21699546342279064 6.022639451081885E-10 1.5002178145105112E-07 + Candidate 4 0.2653915623377543 8.683991289001988E-10 0.21699546342279064 + Release 4 Infinity 0 0 + Step 5 0.06907420086992105 0 0 + Candidate 1 0.0362638246951974 7.758593413765183E-10 1.5002178145105112E-07 + Candidate 3 0.07604556476871949 6.59955603932088E-10 0.0362638246951974 + Candidate 6 0.19387521598218882 5.480900325737565E-11 0.07604556476871949 + Release 6 Infinity 0 0 + Step 0 0.03215405619799004 0 0 + Candidate 1 0.004646654706448317 7.758594326420045E-10 1.5002178145105112E-07 + Candidate 2 0.0084837019094266 9.999787174674693E-10 0.004646654706448317 + Candidate 5 0.010079730302095413 9.812607064514548E-10 0.0084837019094266 + Release 5 Infinity 0 0 + Candidate 1 0.0032877034464036115 6.98032362558651E-10 1.5002178145105112E-07 + Release 1 Infinity 0 0 + Candidate 3 0.0020534762807074003 6.004885049847546E-10 1.5002178145105112E-07 + Release 3 Infinity 0 0 + Candidate 8 0.0014518519383273087 9.294667932478118E-10 1.5002178145105112E-07 + Release 8 Infinity 0 0 + Kkt -1 1.5002178145105112E-07 0 0 + end Converged 14 00000000 3e64b5b7 3f800000 3477cb66 3e0efa90 3f48c5cd 3f44d339 00000000 3f6c163f 00000000 00000000 00000000 + """, + ["w5-shift5-1e-12"] = """ + Candidate 0 137847.8859493429 3.97908931879844E-10 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Step 0 0.4366038070185627 0 0 + Candidate 4 5151.85545116663 2.402625547421844E-11 1.5002178145105112E-07 + Candidate 6 76194.22899199842 5.600810209156691E-10 5151.85545116663 + Release 6 Infinity 0 0 + Step 6 0.7573476181786473 0 0 + Candidate 4 3065.16054302454 3.0049666978425744E-11 1.5002178145105112E-07 + Candidate 7 11849.84792945953 4.782591048738879E-10 3065.16054302454 + Candidate 9 17935.626996204504 7.122810428200548E-10 11849.84792945953 + Release 9 Infinity 0 0 + Candidate 4 2665.4946252912787 3.120333456889181E-11 1.5002178145105112E-07 + Release 4 Infinity 0 0 + Candidate 5 0.21699546341733367 6.022639451081885E-10 1.5002178145105112E-07 + Candidate 8 0.26539156235958217 8.683991289001988E-10 0.21699546341733367 + Release 8 Infinity 0 0 + Step 9 0.06907420085879441 0 0 + Candidate 5 0.03626382468792144 7.758593413765183E-10 1.5002178145105112E-07 + Candidate 7 0.07604556476871949 6.59955603932088E-10 0.03626382468792144 + Candidate 10 0.19387521598173407 5.480900325737565E-11 0.07604556476871949 + Release 10 Infinity 0 0 + Step 4 0.03215405619816462 0 0 + Candidate 0 0.009507730177574558 1.0214105219138004E-09 1.5002178145105112E-07 + Candidate 9 0.01007973029754794 9.81260706451455E-10 0.009507730177574558 + Release 9 Infinity 0 0 + Candidate 0 0.0016777917876424908 9.191744022622724E-10 1.5002178145105112E-07 + Candidate 5 0.003287703438218159 6.980323624876057E-10 0.0016777917876424908 + Release 5 Infinity 0 0 + Candidate 0 0.0014518512944050599 9.29466768014945E-10 1.5002178145105112E-07 + Candidate 7 0.0020534762884381053 6.004885050882204E-10 0.0014518512944050599 + Release 7 Infinity 0 0 + Candidate 0 0.0014518519615194236 9.294667933556894E-10 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Kkt -1 1.5002178145105112E-07 0 0 + end Converged 14 3f6c163f 00000000 00000000 00000000 00000000 3e64b5b7 3f800000 3477cb66 3e0efa90 3f48c5cd 3f44d339 00000000 + """, + ["w5-shift7-1e-12"] = """ + Candidate 1 67022.28908552075 1.9346519011515402E-10 1.5002178145105112E-07 + Candidate 2 137847.8859493429 3.97908931879844E-10 67022.28908552075 + Release 2 Infinity 0 0 + Step 2 0.4366038070185627 0 0 + Candidate 1 37467.36481177385 2.7877783940654753E-10 1.5002178145105112E-07 + Candidate 8 76194.22899199842 5.600810209156691E-10 37467.36481177385 + Release 8 Infinity 0 0 + Step 8 0.7573476181786473 0 0 + Candidate 1 8650.132742734299 3.6196108123398433E-10 1.5002178145105112E-07 + Candidate 9 11849.84792945953 4.782591048738879E-10 8650.132742734299 + Candidate 11 17935.626996204504 7.122810428200548E-10 11849.84792945953 + Release 11 Infinity 0 0 + Candidate 0 610.3107119330159 4.058698497767296E-11 1.5002178145105112E-07 + Candidate 6 2665.4946252912787 3.120333456889181E-11 610.3107119330159 + Release 6 Infinity 0 0 + Candidate 0 0.21130910008525916 4.2238488714754606E-11 1.5002178145105112E-07 + Candidate 7 0.21699546341733367 6.022639451081885E-10 0.21130910008525916 + Candidate 10 0.26539156235958217 8.683991289001988E-10 0.21699546341733367 + Release 10 Infinity 0 0 + Step 11 0.06907420085879441 0 0 + Candidate 0 0.19387521598173407 5.480900325737565E-11 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Step 6 0.032154056198701104 0 0 + Candidate 2 0.00950773017393658 1.0214105219138004E-09 1.5002178145105112E-07 + Candidate 11 0.01007973030937137 9.81260706451455E-10 0.00950773017393658 + Release 11 Infinity 0 0 + Candidate 2 0.0016777917990111746 9.191744022547622E-10 1.5002178145105112E-07 + Candidate 7 0.003287703439127654 6.980323624818886E-10 0.0016777917990111746 + Release 7 Infinity 0 0 + Candidate 2 0.0014518512834911235 9.2946676790583E-10 1.5002178145105112E-07 + Candidate 9 0.0020534762988972943 6.004885050176478E-10 0.0014518512834911235 + Release 9 Infinity 0 0 + Candidate 2 0.0014518519601551816 9.294667935855503E-10 1.5002178145105112E-07 + Release 2 Infinity 0 0 + Kkt -1 1.5002178145105112E-07 0 0 + end Converged 14 3f44d339 00000000 3f6c163f 00000000 00000000 00000000 00000000 3e64b5b7 3f800000 3477cb66 3e0efa90 3f48c5cd + """, + ["w5-shift11-1e-12"] = """ + Candidate 0 135111.75592265668 3.900108740394563E-10 1.5002178145105112E-07 + Candidate 6 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 6 Infinity 0 0 + Step 6 0.4366038070185627 0 0 + Candidate 0 76194.22899199842 5.600810209156691E-10 1.5002178145105112E-07 + Release 0 Infinity 0 0 + Step 0 0.7573476181786473 0 0 + Candidate 1 11849.84792945953 4.782591048738879E-10 1.5002178145105112E-07 + Candidate 3 17935.626996204504 7.122810428200547E-10 11849.84792945953 + Release 3 Infinity 0 0 + Candidate 2 776.2927740423111 8.68340538093788E-10 1.5002178145105112E-07 + Candidate 10 2665.4946252912787 3.120333456889181E-11 776.2927740423111 + Release 10 Infinity 0 0 + Candidate 2 0.26539156234684924 8.683991289001988E-10 1.5002178145105112E-07 + Release 2 Infinity 0 0 + Step 3 0.06907420089782461 0 0 + Candidate 1 0.0760455647832714 6.59955603932088E-10 1.5002178145105112E-07 + Candidate 4 0.1938752159811088 5.480900325737566E-11 0.0760455647832714 + Release 4 Infinity 0 0 + Step 10 0.03215405619865555 0 0 + Candidate 0 0.008483701902150642 9.999787174674693E-10 1.5002178145105112E-07 + Candidate 3 0.01007973030937137 9.812607064514548E-10 0.008483701902150642 + Release 3 Infinity 0 0 + Candidate 1 0.001579840580234304 5.938316649352929E-10 1.5002178145105112E-07 + Candidate 6 0.0016777917771833017 9.191744021382341E-10 0.001579840580234304 + Candidate 11 0.0032877034464036115 6.980323623931817E-10 0.0016777917771833017 + Release 11 Infinity 0 0 + Candidate 1 0.002053476281616895 6.004885051443822E-10 1.5002178145105112E-07 + Release 1 Infinity 0 0 + Candidate 6 0.0014518519456032664 9.29466793204687E-10 1.5002178145105112E-07 + Release 6 Infinity 0 0 + Kkt -1 1.5002178145105112E-07 0 0 + end Converged 14 3f800000 3477cb66 3e0efa90 3f48c5cd 3f44d339 00000000 3f6c163f 00000000 00000000 00000000 00000000 3e64b5b7 + """, + ["w5-default"] = """ + Candidate 1 135111.75592265668 3.900108740394563E-10 0.1500217814510511 + Candidate 7 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 7 Infinity 0 0 + Step 7 0.4366038070185627 0 0 + Candidate 1 76194.22899199842 5.600810209156691E-10 0.1500217814510511 + Release 1 Infinity 0 0 + Step 1 0.7573476181786473 0 0 + Candidate 2 11849.84792945953 4.782591048738879E-10 0.1500217814510511 + Candidate 4 17935.626996204504 7.122810428200547E-10 11849.84792945953 + Release 4 Infinity 0 0 + Candidate 0 633.4374486887755 6.019496088675069E-10 0.1500217814510511 + Candidate 3 776.2927740423111 8.68340538093788E-10 633.4374486887755 + Candidate 11 2665.4946252912787 3.120333456889181E-11 776.2927740423111 + Release 11 Infinity 0 0 + Candidate 0 0.2169954634180158 6.022639451081885E-10 0.1500217814510511 + Candidate 3 0.26539156234684924 8.683991289001988E-10 0.2169954634180158 + Release 3 Infinity 0 0 + Step 4 0.06907420089782461 0 0 + Candidate 5 0.1938752159811088 5.480900325737566E-11 0.1500217814510511 + Release 5 Infinity 0 0 + Step 11 0.03215405619865555 0 0 + Kkt -1 0.1500217814510511 0 0 + end Converged 10 00000000 3f800000 00000000 3f122168 3f800000 3ee2570c 00000000 3f800000 00000000 00000000 00000000 00000000 + """, + ["w5-1e-9"] = """ + Candidate 1 135111.75592265668 3.900108740394563E-10 0.00015002178145105113 + Candidate 7 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 7 Infinity 0 0 + Step 7 0.4366038070185627 0 0 + Candidate 1 76194.22899199842 5.600810209156691E-10 0.00015002178145105113 + Release 1 Infinity 0 0 + Step 1 0.7573476181786473 0 0 + Candidate 2 11849.84792945953 4.782591048738879E-10 0.00015002178145105113 + Candidate 4 17935.626996204504 7.122810428200547E-10 11849.84792945953 + Release 4 Infinity 0 0 + Candidate 0 633.4374486887755 6.019496088675069E-10 0.00015002178145105113 + Candidate 3 776.2927740423111 8.68340538093788E-10 633.4374486887755 + Candidate 11 2665.4946252912787 3.120333456889181E-11 776.2927740423111 + Release 11 Infinity 0 0 + Candidate 0 0.2169954634180158 6.022639451081885E-10 0.00015002178145105113 + Candidate 3 0.26539156234684924 8.683991289001988E-10 0.2169954634180158 + Release 3 Infinity 0 0 + Step 4 0.06907420089782461 0 0 + Candidate 0 0.0362638246792244 7.758593413765183E-10 0.00015002178145105113 + Candidate 2 0.0760455647832714 6.59955603932088E-10 0.0362638246792244 + Candidate 5 0.1938752159811088 5.480900325737566E-11 0.0760455647832714 + Release 5 Infinity 0 0 + Step 11 0.03215405619865555 0 0 + Candidate 0 0.004646654699172359 7.758594326420045E-10 0.00015002178145105113 + Candidate 1 0.008483701902150642 9.999787174674693E-10 0.004646654699172359 + Candidate 4 0.01007973030937137 9.812607064514548E-10 0.008483701902150642 + Release 4 Infinity 0 0 + Candidate 0 0.0032877034464036115 6.980323623931817E-10 0.00015002178145105113 + Release 0 Infinity 0 0 + Candidate 2 0.002053476287983358 6.004885051312699E-10 0.00015002178145105113 + Release 2 Infinity 0 0 + Candidate 7 0.0014518519601551816 9.294667933498287E-10 0.00015002178145105113 + Release 7 Infinity 0 0 + Kkt -1 0.00015002178145105113 0 0 + end Converged 14 3e64b5b7 3f800000 3477cb66 3e0efa90 3f48c5cd 3f44d339 00000000 3f6c163f 00000000 00000000 00000000 00000000 + """, + ["w5-cap-3"] = """ + Candidate 1 135111.75592265668 3.900108740394563E-10 1.5002178145105112E-07 + Candidate 7 137847.8859493429 3.97908931879844E-10 135111.75592265668 + Release 7 Infinity 0 0 + Step 7 0.4366038070185627 0 0 + Candidate 1 76194.22899199842 5.600810209156691E-10 1.5002178145105112E-07 + Release 1 Infinity 0 0 + Cap -1 0 0 0 + end IterationLimit 3 00000000 00000000 00000000 00000000 00000000 00000000 00000000 3f800000 00000000 00000000 00000000 00000000 + """, + ["confirm-pass-1"] = """ + Candidate 0 1.579785345010647 1.070372468301927E-15 3.4029892871621996E-11 + Candidate 1 267440.9365187891 1.7819678819499464E-10 1.579785345010647 + Release 1 Infinity 0 0 + Miss -1 -0.013847553852146902 Infinity 1.3847553852146903E-17 + Kkt -1 3.4029892871621996E-11 0 0 + end Converged 2 bf800000 b7d8d9d8 + """, + ["confirm-pass-2"] = """ + Candidate 0 459906.13776790636 3.064812525589559E-10 9.174025280630848E-11 + Release 0 Infinity 0 0 + Miss -1 -0.009148166097177548 Infinity 9.148166097177549E-18 + Kkt -1 9.174025280630848E-11 0 0 + end Converged 2 b9511fc3 00000000 + """, + ["noise-one-heavy"] = """ + Step 2 0.024621369447167075 0 0 + Candidate 5 0.012439710770288478 7.959493583159782E-17 0 + Release 5 Infinity 0 0 + Miss -1 6.382790758681003E-20 Infinity 3.5111012909447177E-34 + Candidate 2 1.7763568394002505E-14 5.656022558163373E-13 0 + Candidate 3 1.5178830414797062E-18 6.655962418338423E-17 0 + Miss -1 4.449537197592434E-20 6.382790758681003E-20 4.4495371975924343E-35 + Candidate 2 2.842170943040401E-14 5.656022558163373E-13 0 + Candidate 3 2.6020852139652106E-18 6.655962418338423E-17 0 + Miss -1 4.3626178182815784E-20 4.449537197592434E-20 8.122105116216791E-35 + Candidate 2 2.4868995751603507E-14 5.656022558163373E-13 0 + Candidate 3 2.168404344971009E-18 6.655962418338423E-17 0 + Stall -1 4.5923297824441354E-20 4.3626178182815784E-20 5.911346354815417E-35 + end Converged 5 bcd3d632 00000000 b56215c6 bd1ed80d 00000000 3c154621 + """, + ["noise-one-light"] = """ + Step 2 0.024621369199885876 0 0 + Candidate 5 12.439710484934402 7.959493151937227E-14 0 + Release 5 Infinity 0 0 + Candidate 2 1.164670302474663E-20 5.656022762441261E-20 0 + Candidate 3 2.0650148258027912E-14 6.655962175358023E-14 0 + Miss -1 6.070552180371575E-20 Infinity 8.382049277624677E-34 + Candidate 2 4.235164736271502E-22 5.656022762441263E-20 0 + Candidate 3 2.6645352591003757E-15 6.655962175358026E-14 0 + Miss -1 2.7840962879986504E-21 6.070552180371575E-20 1.0613730890584895E-34 + Stall -1 3.341824384345746E-20 2.7840962879986504E-21 5.303670554538265E-35 + end Converged 4 b7d8ebb8 00000000 c106c1d7 b822a7fd 00000000 3718db44 + """, + ["noise-ascending"] = """ + Step 2 0.024621365926673414 0 0 + Candidate 5 1243.970980681693 7.959492843486517E-12 0 + Release 5 Infinity 0 0 + Candidate 2 8.326672684688674E-17 5.6560222323987015E-15 0 + Candidate 3 3.1086244689504383E-15 6.655962152994665E-14 0 + Miss -1 2.4960527060967552E-20 Infinity 1.0260573563230584E-34 + Candidate 3 1.5543122344752192E-15 6.655962152994665E-14 0 + Kkt -1 0 0 0 + end Converged 3 bcd3d632 00000000 b8b0a103 b822a7fd 00000000 33c3a805 + """, + ["near-singular-3"] = """ + Kkt -1 4.011347176225968E-07 0 0 + end Converged 0 cb42b27d cacc2e2d 49ac1c2a cae3caf3 4a4a8d2d 4ac85877 49e01bee 4a1115bf ca628e52 cb4f0de4 c9a80340 ca273e77 caa50c4a 498330e8 cafde285 ca5a8f4e ca576c89 ca517d14 4aa437c5 00000000 cadb0117 ca3073e4 4aab2053 4a273cf8 + """, + ["near-singular-11"] = """ + Kkt -1 8.26586030599639E-07 0 0 + end Converged 0 c9e0f987 c999e95d c9305c40 49c39403 ca2d4887 cadb1864 49752be6 ca9b41c5 ca3d7042 4a1f2341 48f59fdc 4a865715 c950dbee 4aa1e545 49f2e233 4a19c6b5 c9865f16 4a9d9500 ca24f675 4a1f7a2c ca84489a 4a7b50aa ca34c295 00000000 + """, + ["below-noise-2"] = """ + Miss -1 -0.49773013750831796 Infinity 4.994142446397116E-16 + Stall -1 -0.49747185997938986 -0.49773013750831796 4.992128425425914E-16 + end Converged 1 4a0f075c 4a068153 ca43258d c98493ab 4a2dd375 4a1d0b0b 4969536f 48b0cc93 c99ce69b c9be826d c827f0d4 c690e0d3 c9e0c560 485f7251 49c00e74 4939fadc ca264bd4 4a21483e c93e1efa 4a0a6b52 ca401661 c9df1aa4 c97ec1d5 4a134dbe + """, + ["below-noise-9"] = """ + Miss -1 -0.49865831712645864 Infinity 5.014509747560003E-16 + Stall -1 -0.49841675352897546 -0.49865831712645864 5.014635585960667E-16 + end Converged 1 c989ed22 4a0e7096 c98bf358 498c0fda 4a59e43c c82ef778 ca09dc93 c9f6326b 49feecdd c84efa98 c9155228 cac57da7 49c5dd82 49934cfe c82a770a 496d11ec 4a1a9b15 49084d80 485f9b84 4a12918f 4a0b2594 49e90cd5 c9f587a0 c92ccca9 + """, + ["boxed-conditioned"] = """ + Step 6 1.7950541723871097E-06 0 0 + Step 3 7.627008876179801E-07 0 0 + Step 7 1.2561155441158494E-05 0 0 + Step 2 2.5030117160256515E-05 0 0 + Step 4 0.00010110066197668991 0 0 + Step 0 0.0009920440120834782 0 0 + Step 1 0.004498393465636397 0 0 + Step 5 0.01208289710393302 0 0 + Candidate 3 0.14713995573185712 3.023980745569078E-16 6.053930821702194E-10 + Candidate 4 0.49326802296517885 1.1111503645474022E-15 0.14713995573185712 + Release 4 Infinity 0 0 + Step 4 0.05208634850970442 0 0 + Candidate 3 0.14559191607117922 3.023980745569078E-16 6.053930821702194E-10 + Release 3 Infinity 0 0 + Step 3 0.0045758825891475145 0 0 + Candidate 6 0.03694908965828074 8.947385886461852E-17 6.053930821702194E-10 + Candidate 7 0.04406217168160146 1.5294423693380426E-16 0.03694908965828074 + Release 7 Infinity 0 0 + Step 7 0.09930368626410971 0 0 + Candidate 6 0.035916303301579564 8.947385886461852E-17 6.053930821702194E-10 + Release 6 Infinity 0 0 + Step 6 0.007199613540050211 0 0 + Kkt -1 6.053930821702194E-10 0 0 + end Converged 16 bd4ccccd bd4ccccd 3d4ccccd bd4ccccd bd4ccccd 3d4ccccd 3d4ccccd bd4ccccd + """, + ["random-5"] = """ + Candidate 1 2.3674161062532875 2.1026878960246686E-15 2.367416106253287E-06 + Release 1 Infinity 0 0 + Candidate 2 1.1521571531355912 1.0233211195181646E-15 2.367416106253287E-06 + Release 2 Infinity 0 0 + Kkt -1 2.367416106253287E-06 0 0 + end Converged 2 00000000 3eed0d99 3e4a4f4e + """, + ["random-5-warm-hi"] = """ + Candidate 0 2.5083502579763883 2.22786256818387E-15 2.367416106253287E-06 + Candidate 2 4.168064738401132 7.758404534582419E-15 2.5083502579763883 + Release 2 Infinity 0 0 + Candidate 0 2.3851527197990436 2.11844117341457E-15 2.367416106253287E-06 + Candidate 1 3.0072430750429278 8.055628322575466E-15 2.3851527197990436 + Release 1 Infinity 0 0 + Candidate 0 1.678136149614912 1.4904843134066255E-15 2.367416106253287E-06 + Release 0 Infinity 0 0 + Step 0 0.24975695258808922 0 0 + Kkt -1 2.367416106253287E-06 0 0 + end Converged 4 00000000 3eed0d99 3e4a4f4e + """, + ["random-17"] = """ + Step 4 0.1232927301830878 0 0 + Step 2 0.7178536883498736 0 0 + Kkt -1 3.93442891781855E-06 0 0 + end Converged 2 beb31769 bee93450 bf800000 00000000 3dcccccd 00000000 bf432cd8 3c1e7a75 + """, + ["random-17-warm-hi"] = """ + Candidate 0 1.1914277145018732 1.1651190341326312E-14 3.93442891781855E-06 + Candidate 2 8.918271049237749 2.2324704556906607E-14 1.1914277145018732 + Release 2 Infinity 0 0 + Step 2 0.9608581392999689 0 0 + Candidate 0 0.12375707360068922 1.1651190341326312E-14 3.93442891781855E-06 + Candidate 1 0.46721760431143355 1.776136934146532E-14 0.12375707360068922 + Candidate 3 5.9589017943044995 1.7762048080698952E-14 0.46721760431143355 + Candidate 5 6.330922749270252 2.0454695850835226E-14 5.9589017943044995 + Release 5 Infinity 0 0 + Step 5 0.6196380883752669 0 0 + Candidate 1 1.2265685104711581 1.6243881393846015E-14 3.93442891781855E-06 + Candidate 3 4.91009654542079 1.5666114156513806E-14 1.2265685104711581 + Release 3 Infinity 0 0 + Step 3 0.4723422790325044 0 0 + Candidate 1 2.313755717662308 1.4071244908743248E-14 3.93442891781855E-06 + Release 1 Infinity 0 0 + Candidate 2 0.54389362384502 1.823003362441137E-14 3.93442891781855E-06 + Candidate 6 2.3746920633686277 6.957925871969801E-15 0.54389362384502 + Release 6 Infinity 0 0 + Candidate 7 1.5652291166895491 9.329087077889423E-15 3.93442891781855E-06 + Release 7 Infinity 0 0 + Candidate 0 0.9877280584738057 5.425941552304341E-15 3.93442891781855E-06 + Release 0 Infinity 0 0 + Kkt -1 3.93442891781855E-06 0 0 + end Converged 10 beb31769 bee93450 bf800000 00000000 3dcccccd 00000000 bf432cd8 3c1e7a75 + """, + ["random-40"] = """ + Step 0 0.6917697000772063 0 0 + Step 1 0.08710458252709695 0 0 + Step 2 0.04883630084124485 0 0 + Candidate 4 0.9459344738504001 6.568402131960035E-15 9.16596182404494E-06 + Release 4 Infinity 0 0 + Step 4 0.29527820149229345 0 0 + Kkt -1 9.16596182404494E-06 0 0 + end Converged 5 3f800000 bf800000 bf800000 00000000 3dcccccd + """, + ["random-40-warm-hi"] = """ + Candidate 1 4.718545111140088 9.652651618643102E-15 9.16596182404494E-06 + Candidate 3 5.026292374797416 1.201620793157454E-14 4.718545111140088 + Release 3 Infinity 0 0 + Step 3 0.08526897635819114 0 0 + Candidate 1 4.75316070860798 9.60653437865909E-15 9.16596182404494E-06 + Release 1 Infinity 0 0 + Candidate 2 3.678831066315437 4.901187544235657E-15 9.16596182404494E-06 + Release 2 Infinity 0 0 + Step 1 0.525739055641334 0 0 + Step 2 0.74742310322831 0 0 + Kkt -1 9.16596182404494E-06 0 0 + end Converged 6 3f800000 bf800000 bf800000 00000000 3dcccccd + """, + }; +} diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs index 4678006f..136e6e1b 100644 --- a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs @@ -942,48 +942,6 @@ 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, 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. @@ -1016,11 +974,11 @@ public sealed class BoundedLeastSquaresTests float[] bWeak = null; if (kind == "weak") { - (a, bWeak) = ConditionedProblem(m, n, Math.Pow(10.0, 6.0 + rng.NextFloat()), seed, 1.0, true); + (a, bWeak) = BlsTrace.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; + a = BlsTrace.ConditionedProblem(m, n, Math.Pow(10.0, 5.0 + 0.5 * rng.NextFloat()), seed).a; } else { @@ -1092,7 +1050,7 @@ public sealed class BoundedLeastSquaresTests const int m = 28, n = 24; for (uint seed = 1; seed <= 30; seed++) { - var (a, b) = ConditionedProblem(m, n, cond, seed); + var (a, b) = BlsTrace.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}"); @@ -1116,7 +1074,7 @@ public sealed class BoundedLeastSquaresTests var confirmed = 0; for (uint seed = 1; seed <= 60; seed++) { - var (a, b) = ConditionedProblem(m, n, 3e7, seed); + var (a, b) = BlsTrace.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}"); @@ -1139,7 +1097,7 @@ public sealed class BoundedLeastSquaresTests var confirmed = 0; for (uint seed = 1; seed <= 30; seed++) { - var (a, b) = ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: true); + var (a, b) = BlsTrace.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}"); @@ -1159,7 +1117,7 @@ public sealed class BoundedLeastSquaresTests const int m = 28, n = 24; for (uint seed = 1; seed <= 30; seed++) { - var (a, b) = ConditionedProblem(m, n, cond, seed, 1e6); + var (a, b) = BlsTrace.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}"); @@ -1553,7 +1511,7 @@ public sealed class BoundedLeastSquaresTests [InlineData(1e-9)] public void TheStopAfterAReleaseIsOnTheGradientNotTheCost(double tolerance) { - var (m, n, a, b, lo, hi) = AllocatorW5(0); + var (m, n, a, b, lo, hi) = BlsTrace.AllocatorW5(0); var cold = new float[n]; Assert.Equal(BoundedLeastSquaresStatus.Converged, SolveAt(m, n, a, b, lo, hi, cold, out _, tolerance)); Assert.Equal(0.92221445f, cold[7], 1e-4); @@ -1575,7 +1533,7 @@ public sealed class BoundedLeastSquaresTests const int n = 12; for (var shift = 0; shift < n; shift++) { - var (m, _, a, b, lo, hi) = AllocatorW5(shift); + var (m, _, a, b, lo, hi) = BlsTrace.AllocatorW5(shift); var cold = new float[n]; Assert.True(BoundedLeastSquaresStatus.Converged == SolveAt(m, n, a, b, lo, hi, cold, out _, tolerance), $"shift {shift}"); Assert.True(Math.Abs(cold[(7 + shift) % n] - 0.92221445f) <= 1e-4, $"shift {shift}: throttle 7 is {cold[(7 + shift) % n]}"); @@ -1586,33 +1544,6 @@ public sealed class BoundedLeastSquaresTests } } - // Column j of the returned problem is column (j - shift) mod 12 of Soul's W5. - private static (int m, int n, float[] a, float[] b, float[] lo, float[] hi) AllocatorW5(int shift) - { - const int m = 15, n = 12; - var w = new float[m * n]; - w[0 * n + 2] = 100f; w[0 * n + 8] = -100f; w[0 * n + 9] = 100f; - var row1 = new[] { 13.882518f, 20.86695f, 0f, 14.8965845f, 12.38046f, 22.115599f, -13.706654f, 15.85789f, 0f, 0f, -100f, 100f }; - var row2 = new[] { -187.72859f, 239.31168f, 158.87828f, -270.31165f, 235.38467f, -14.704258f, 121.61232f, 244.81305f }; - for (var j = 0; j < n; j++) w[1 * n + j] = row1[j]; - for (var j = 0; j < row2.Length; j++) w[2 * n + j] = row2[j]; - for (var k = 0; k < 4; k++) w[(3 + k) * n + 8 + k] = 0.3f; - for (var i = 0; i < 8; i++) w[(7 + i) * n + i] = 0.1f; - var b = new float[m]; - b[1] = 67.385445f; b[2] = 558.70917f; - for (var k = 0; k < 4; k++) b[3 + k] = -3f; - var lo = new float[n]; - var hi = new float[n]; - var a = new float[m * n]; - for (var j = 0; j < n; j++) - { - var to = (j + shift) % n; - hi[to] = j < 8 ? 1f : 2f; - for (var r = 0; r < m; r++) a[r * n + to] = w[r * n + j]; - } - return (m, n, a, b, lo, hi); - } - // b = 0 with the origin inside the box: the optimum is any x with Ax = 0, the gradient scale is 0 and so is the tolerance, and // every gradient component the solve leaves behind is rounding. Releasing on such a violation cycles release / // stationarity miss / re-bind until the iteration cap (the plain-batch problem at seed 141 in miniature). @@ -1670,7 +1601,7 @@ public sealed class BoundedLeastSquaresTests var worst = 0; for (uint seed = 1; seed <= 30; seed++) { - var (a, b) = ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: true); + var (a, b) = BlsTrace.ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: true); var x = new float[n]; var status = SolveAt(m, n, a, b, Fill(-1e9f, n), Fill(1e9f, n), x, out var iterations, tol); Assert.True(status == BoundedLeastSquaresStatus.Converged, $"cond {cond} tol {tol} seed {seed}: {status} after {iterations}"); @@ -1790,7 +1721,7 @@ public sealed class BoundedLeastSquaresTests foreach (var along in new[] { true, false }) for (uint seed = 1; seed <= 30; seed++) { - var (a, b) = ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: along); + var (a, b) = BlsTrace.ConditionedProblem(m, n, cond, seed, 1.0, bAlongSmallestSingularVector: along); var byDefault = new float[n]; var tight = new float[n]; var s1 = Solve(m, n, a, b, Fill(-1e9f, n), Fill(1e9f, n), byDefault, out var i1); diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTrace.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTrace.cs new file mode 100644 index 00000000..71d05c7b --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTrace.cs @@ -0,0 +1,223 @@ +using System; +using System.Collections.Generic; +using System.Globalization; +using System.Text; +using CultMath; + +namespace CultMath.Tests; + +// The solver's decisions, recorded in order through the internal ISolveTrace seam, and the problems they are recorded on. +// BlsTrace is the only place that knows how a trace is written; BoundedLeastSquaresTraceTests compares it with the trace +// the base commit (dbcb2a8e) produced for the same problems, and with hand-built decision inputs. +internal static class BlsTrace +{ + internal sealed class Recorder : ISolveTrace + { + public readonly List Lines = new(); + + public void Record(SolveEvent kind, int column, double a, double b, double c) + => Lines.Add($"{kind} {column} {F(a)} {F(b)} {F(c)}"); + } + + // Round-trip text: equal strings are equal doubles. + internal static string F(double v) => v.ToString("R", CultureInfo.InvariantCulture); + + internal sealed record Case(string Name, int M, int N, float[] A, float[] B, float[] Lo, float[] Hi, float[] X0, double Tolerance, int MaxIterations); + + // The solve of one case: every decision, then the status, the iteration count and the bits of x. + internal static string Trace(Case c) + { + var rec = new Recorder(); + var x = (float[])c.X0.Clone(); + var status = BoundedLeastSquares.SolveTraced(c.M, c.N, c.A, c.B, c.Lo, c.Hi, x, + new float[BoundedLeastSquares.WorkspaceLength(c.N)], out var iterations, c.Tolerance, c.MaxIterations, rec); + var sb = new StringBuilder(); + foreach (var line in rec.Lines) sb.Append(line).Append('\n'); + sb.Append($"end {status} {iterations}"); + foreach (var v in x) sb.Append(' ').Append(BitConverter.SingleToInt32Bits(v).ToString("x8", CultureInfo.InvariantCulture)); + return sb.ToString(); + } + + private static float[] Fill(float value, int n) + { + var v = new float[n]; + Array.Fill(v, value); + return v; + } + + private static Case Make(string name, int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[]? x0 = null, + double tol = BoundedLeastSquares.DefaultKktRelativeTolerance, int max = BoundedLeastSquares.DefaultMaxIterations) + => new(name, m, n, a, b, lo, hi, x0 ?? (float[])lo.Clone(), tol, max); + + internal static IEnumerable Cases() + { + // Equal violations: which column is released first is the tie-break (min |x0 + x1 - 1| on [0,1]^2, then three copies). + yield return Make("tie-pair", 1, 2, new[] { 1f, 1f }, new[] { 1f }, new float[2], Fill(1f, 2)); + yield return Make("tie-triple", 2, 3, Fill(1f, 6), new[] { 1f, 1f }, new float[3], Fill(1f, 3)); + yield return Make("tie-pair-at-hi", 1, 2, new[] { 1f, 1f }, new[] { -1f }, new float[2], Fill(1f, 2), Fill(1f, 2)); + + // A step that lands exactly on a bound is a full step, not a blocked one. + yield return Make("step-lands-on-bound", 1, 1, new[] { 1f }, new[] { 1f }, new[] { 0f }, new[] { 1f }, new[] { 0.5f }); + yield return Make("step-blocked", 2, 2, new[] { 1f, 0f, 0f, 1f }, new[] { 3f, -0.5f }, new[] { 0f, -1f }, new[] { 1f, 1f }, new[] { 0.5f, 0.5f }); + + // Soul's W5, with throttle 7 released from every column position in turn (column 0 included). + foreach (var shift in new[] { 0, 1, 5, 7, 11 }) + { + var (m, n, a, b, lo, hi) = AllocatorW5(shift); + yield return Make($"w5-shift{shift}-1e-12", m, n, a, b, lo, hi, new float[n], 1e-12); + } + { + var (m, n, a, b, lo, hi) = AllocatorW5(0); + yield return Make("w5-default", m, n, a, b, lo, hi, new float[n]); + yield return Make("w5-1e-9", m, n, a, b, lo, hi, new float[n], 1e-9); + yield return Make("w5-cap-3", m, n, a, b, lo, hi, new float[n], 1e-12, 3); + } + + // A release followed by a stationarity miss: the first miss is not yet a stall. + yield return Make("confirm-pass-1", 4, 2, + new[] { -0.0067761745f, -207.21786f, -0.004115388f, -751.1994f, -8.1700025E-05f, -481.7794f, 0.0038280683f, 480.12833f }, + new[] { 1.7477255f, 2.4914713f, -2.3861952f, 2.187128f }, new[] { -1f, -0.25f }, new[] { 1f, 0.1f }, tol: 1e-12); + yield return Make("confirm-pass-2", 5, 2, + new[] { 70.07762f, 5.7015634f, 111.790436f, -12.586686f, 75.29796f, 12.872999f, -402.25464f, -1.8093401f, -524.51056f, -8.649797f }, + new[] { -1.8312539f, 2.9556413f, 0.7550676f, -2.1365998f, 2.30717f }, new[] { -1f, 0f }, new[] { 1f, 1f }, tol: 1e-12); + + // The noise-level problems: columns rescaled so the gate is the only thing between a release and a spin. + { + const int m = 2, n = 6; + var a0 = new[] + { + 1.0213068f, 0.97530603f, -0.5673449f, -0.64136946f, -0.07432885f, -0.35571396f, + 0.4326824f, 1.1969734f, -0.1343108f, -0.23431005f, -1.4009786f, 0.10656489f, + }; + var lo0 = new[] { -0.44782114f, 0f, -0.008422324f, -0.03878026f, 0f, -6.6434336f }; + var hi0 = new[] { 0.0448007f, 0f, 0.005942123f, 0.0785512f, 0f, 0.23079813f }; + var x0 = new[] { -0.28347927f, 0f, 0.0051532285f, -0.03878026f, 0f, 0.23079813f }; + foreach (var (name, scale) in new (string, float[])[] + { + ("one-heavy", new[] { 1f, 1f, 1e4f, 1f, 1f, 1f }), + ("one-light", new[] { 1e3f, 1e3f, 1e-3f, 1e3f, 1e3f, 1e3f }), + ("ascending", new[] { 1f, 1e1f, 1e2f, 1e3f, 1e4f, 1e5f }), + }) + { + var a = new float[m * n]; + var lo = new float[n]; + var hi = new float[n]; + var x = new float[n]; + for (var j = 0; j < n; j++) + { + for (var i = 0; i < m; i++) a[i * n + j] = a0[i * n + j] * scale[j]; + lo[j] = lo0[j] / scale[j]; + hi[j] = hi0[j] / scale[j]; + x[j] = Math.Clamp(x0[j] / scale[j], lo[j], hi[j]); + } + yield return Make($"noise-{name}", m, n, a, new float[m], lo, hi, x); + } + } + + // Near-singular and below-the-noise-floor free sets: stationarity misses and stalls. + foreach (var seed in new uint[] { 3, 11 }) + { + var (a, b) = ConditionedProblem(28, 24, 3e7, seed); + yield return Make($"near-singular-{seed}", 28, 24, a, b, Fill(-1e9f, 24), Fill(1e9f, 24), new float[24]); + } + foreach (var seed in new uint[] { 2, 9 }) + { + var (a, b) = ConditionedProblem(28, 24, 1e7, seed, 1.0, bAlongSmallestSingularVector: true); + yield return Make($"below-noise-{seed}", 28, 24, a, b, Fill(-1e9f, 24), Fill(1e9f, 24), new float[24]); + } + { + var (a, b) = ConditionedProblem(12, 8, 1e5, 4); + yield return Make("boxed-conditioned", 12, 8, a, b, Fill(-0.05f, 8), Fill(0.05f, 8), new float[8], 1e-9); + } + + // Boxed random problems, cold and from a far warm start. + foreach (var seed in new uint[] { 5, 17, 40 }) + { + var rng = new CultMath.Random(seed * 7919u); + var n = rng.NextInt(2, 9); + var m = rng.NextInt(n, 14); + 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.NextInt(2) == 0 ? -1f : 0f; + hi[j] = rng.NextInt(2) == 0 ? 1f : 0.1f; + } + yield return Make($"random-{seed}", m, n, a, b, lo, hi, new float[n]); + yield return Make($"random-{seed}-warm-hi", m, n, a, b, lo, hi, (float[])hi.Clone()); + } + } + + // 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. + internal 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); + } + + // Column j of the returned problem is column (j - shift) mod 12 of Soul's W5. + internal static (int m, int n, float[] a, float[] b, float[] lo, float[] hi) AllocatorW5(int shift) + { + const int m = 15, n = 12; + var w = new float[m * n]; + w[0 * n + 2] = 100f; w[0 * n + 8] = -100f; w[0 * n + 9] = 100f; + var row1 = new[] { 13.882518f, 20.86695f, 0f, 14.8965845f, 12.38046f, 22.115599f, -13.706654f, 15.85789f, 0f, 0f, -100f, 100f }; + var row2 = new[] { -187.72859f, 239.31168f, 158.87828f, -270.31165f, 235.38467f, -14.704258f, 121.61232f, 244.81305f }; + for (var j = 0; j < n; j++) w[1 * n + j] = row1[j]; + for (var j = 0; j < row2.Length; j++) w[2 * n + j] = row2[j]; + for (var k = 0; k < 4; k++) w[(3 + k) * n + 8 + k] = 0.3f; + for (var i = 0; i < 8; i++) w[(7 + i) * n + i] = 0.1f; + var b = new float[m]; + b[1] = 67.385445f; b[2] = 558.70917f; + for (var k = 0; k < 4; k++) b[3 + k] = -3f; + var lo = new float[n]; + var hi = new float[n]; + var a = new float[m * n]; + for (var j = 0; j < n; j++) + { + var to = (j + shift) % n; + hi[to] = j < 8 ? 1f : 2f; + for (var r = 0; r < m; r++) a[r * n + to] = w[r * n + j]; + } + return (m, n, a, b, lo, hi); + } +} diff --git a/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTraceTests.cs b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTraceTests.cs new file mode 100644 index 00000000..6c8cddf5 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTraceTests.cs @@ -0,0 +1,226 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using CultMath; +using Xunit; + +namespace CultMath.Tests; + +// The stopping rule is a set of decisions, so it is pinned where it is decided. The solver reports each decision, in order and +// with the values it compared, through an internal trace sink (ISolveTrace); these tests assert that record with exact +// equality, never a status or a tolerance band. Two sources of truth: +// - BlsBaseTraces: the same decisions as the base commit dbcb2a8e took them on the same problems, so Solve's defaults are +// unchanged decision by decision, not only in the end; +// - hand-built inputs to the two judgement functions, for states a solve may never reach (a release after a stationarity +// miss with a finite stall mark) where only the decision itself can be asked. +public sealed class BoundedLeastSquaresTraceTests +{ + private static string Normal(string s) => s.Replace("\r\n", "\n"); + + [Fact] + public void EveryDecisionEqualsTheBaseCommitsOnTheSameProblem() + { + foreach (var c in BlsTrace.Cases()) + { + Assert.True(BlsBaseTraces.Expected.TryGetValue(c.Name, out var expected), $"{c.Name}: no base trace"); + Assert.Equal(Normal(expected!), BlsTrace.Trace(c)); + } + Assert.Equal(BlsTrace.Cases().Count(), BlsBaseTraces.Expected.Count); + } + + // Each kind of decision, and the ones the stop rule's variants differ on, occurs in the recorded problems: a case table + // that stopped exercising a decision would pin nothing about it. + [Fact] + public void TheCasesExerciseEveryKindOfDecision() + { + var lines = BlsTrace.Cases().SelectMany(c => BlsTrace.Trace(c).Split('\n')).ToList(); + foreach (var kind in Enum.GetValues()) + Assert.True(lines.Any(l => l.StartsWith(kind + " ", StringComparison.Ordinal)), $"{kind} is never recorded"); + Assert.Contains(lines, l => l.StartsWith("Release 0 ", StringComparison.Ordinal)); + Assert.Contains(lines, l => l.StartsWith("Release 11 ", StringComparison.Ordinal)); + } + + // The sink only observes: with it, without it and through the public overload the solve is the same, bit for bit. + [Fact] + public void ARecordingSinkNeverChangesTheSolve() + { + foreach (var c in BlsTrace.Cases()) + { + var plain = (float[])c.X0.Clone(); + var traced = (float[])c.X0.Clone(); + var viaPublic = (float[])c.X0.Clone(); + var ws = new float[BoundedLeastSquares.WorkspaceLength(c.N)]; + var s1 = BoundedLeastSquares.SolveTraced(c.M, c.N, c.A, c.B, c.Lo, c.Hi, plain, ws, out var i1, c.Tolerance, c.MaxIterations, null); + var s2 = BoundedLeastSquares.SolveTraced(c.M, c.N, c.A, c.B, c.Lo, c.Hi, traced, ws, out var i2, c.Tolerance, c.MaxIterations, new BlsTrace.Recorder()); + var s3 = BoundedLeastSquares.Solve(c.M, c.N, c.A, c.B, c.Lo, c.Hi, viaPublic, ws, out var i3, c.Tolerance, c.MaxIterations); + Assert.True(s1 == s2 && s2 == s3 && i1 == i2 && i2 == i3, c.Name); + for (var j = 0; j < c.N; j++) + { + Assert.Equal(BitConverter.SingleToInt32Bits(plain[j]), BitConverter.SingleToInt32Bits(traced[j])); + Assert.Equal(BitConverter.SingleToInt32Bits(plain[j]), BitConverter.SingleToInt32Bits(viaPublic[j])); + } + } + } + + private static string Judge(int release, bool freeViolated, ref double fMark, double[] x, double[] g, double[] atb, out bool over, double kktTol = 0.5) + { + var rec = new BlsTrace.Recorder(); + over = BoundedLeastSquares.JudgePass(release, freeViolated, kktTol, ref fMark, x, g, atb, x.Length, rec); + return string.Join("\n", rec.Lines); + } + + // Cost 0.5 * sum x (g - atb) = 2 and scale 0.5 * sum |x| (|g| + |atb|) = 4 for this point. + private static readonly double[] X1 = { 2.0 }; + private static readonly double[] G1 = { 3.0 }; + private static readonly double[] Atb1 = { 1.0 }; + + // A release is never judged by the cost, whichever column it is: a column-0 release with a stall mark the cost has not + // beaten goes on, and the mark is reset to +Inf. + [Theory] + [InlineData(0)] + [InlineData(1)] + [InlineData(2)] + public void AReleaseOfAnyColumnIsNotJudgedByTheCostAndResetsTheMark(int column) + { + var x = new[] { 2.0, 2.0, 2.0 }; + var g = new[] { 3.0, 3.0, 3.0 }; + var atb = new[] { 1.0, 1.0, 1.0 }; + var mark = 1.0; // below the cost 6: a stationarity miss with this mark would be a stall + var trace = Judge(column, true, ref mark, x, g, atb, out var over); + Assert.False(over); + Assert.Equal(double.PositiveInfinity, mark); + Assert.Equal($"Release {column} Infinity 0 0", trace); + } + + // The reset is to +Inf: not the cost (that would make the next equal cost a stall at once), not nothing. + [Fact] + public void AReleaseResetsTheMarkToPositiveInfinityNotToTheCostNorLeavesItStale() + { + var mark = 5.0; + Judge(2, false, ref mark, X1, G1, Atb1, out var over); + Assert.False(over); + Assert.Equal(double.PositiveInfinity, mark); + + // After the reset a miss at the cost the stale mark would have stopped on is only recorded. + var trace = Judge(-1, true, ref mark, X1, G1, Atb1, out over); + Assert.False(over); + Assert.Equal(2.0, mark); + Assert.Equal($"Miss -1 2 Infinity {BlsTrace.F(1e-15 * 4.0)}", trace); + } + + [Fact] + public void AStationarityMissIsOverOnlyWhenTheCostHasNotFallenByMoreThanTheRoundingAllowance() + { + var allowance = 1e-15 * 4.0; + + var mark = 2.0 + 3e-15; // 2 >= mark - allowance: not enough progress + var trace = Judge(-1, true, ref mark, X1, G1, Atb1, out var over); + Assert.True(over); + Assert.Equal($"Stall -1 2 {BlsTrace.F(2.0 + 3e-15)} {BlsTrace.F(allowance)}", trace); + + mark = 2.0 + 5e-15; // 2 < mark - allowance: progress, so go on and record the cost + trace = Judge(-1, true, ref mark, X1, G1, Atb1, out over); + Assert.False(over); + Assert.Equal(2.0, mark); + Assert.Equal($"Miss -1 2 {BlsTrace.F(2.0 + 5e-15)} {BlsTrace.F(allowance)}", trace); + + mark = 2.0; // equal cost is a stall + trace = Judge(-1, true, ref mark, X1, G1, Atb1, out over); + Assert.True(over); + Assert.Equal($"Stall -1 2 2 {BlsTrace.F(allowance)}", trace); + } + + // The cost and its scale are 0.5 sum x_j (g_j - atb_j) and 0.5 sum |x_j| (|g_j| + |atb_j|) over every column. + [Fact] + public void TheCostAndItsScaleCountEveryColumnAndTheirOwnSigns() + { + var x = new[] { 2.0, -1.0 }; + var g = new[] { 3.0, 1.0 }; + var atb = new[] { 1.0, 2.0 }; + var mark = double.PositiveInfinity; + var trace = Judge(-1, true, ref mark, x, g, atb, out var over); + Assert.False(over); + Assert.Equal($"Miss -1 2.5 Infinity {BlsTrace.F(1e-15 * 5.5)}", trace); + Assert.Equal(2.5, mark); + } + + [Fact] + public void NothingToReleaseAndAStationaryFreeSetIsKktAtOnce() + { + var mark = 1.0; + var trace = Judge(-1, false, ref mark, X1, G1, Atb1, out var over, kktTol: 0.25); + Assert.True(over); + Assert.Equal("Kkt -1 0.25 0 0", trace); + Assert.Equal(1.0, mark); + } + + private static string Select(double[] ata, double[] atb, double[] x, double[] g, int[] state, float[] lo, float[] hi, double kktTol, out int release) + { + var rec = new BlsTrace.Recorder(); + release = BoundedLeastSquares.SelectRelease(ata, atb, x, g, state, lo, hi, kktTol, x.Length, rec); + return string.Join("\n", rec.Lines); + } + + // Equal violations: the first column wins, and the second is never even a candidate. + [Fact] + public void OfEquallyViolatedBoundsTheFirstIsReleased() + { + var ata = new[] { 1.0, 0.0, 0.0, 1.0 }; + var atb = new[] { 0.0, 0.0 }; + var x = new[] { 0.0, 0.0 }; + var trace = Select(ata, atb, x, new[] { -1.0, -1.0 }, new[] { -1, -1 }, new[] { 0f, 0f }, new[] { 1f, 1f }, 0.1, out var release); + Assert.Equal(0, release); + Assert.Equal("Candidate 0 1 0 0.1", trace); + + // The most violated wins wherever it is, and a later equal one does not take it over. + ata = new[] { 1.0, 0, 0, 0, 1.0, 0, 0, 0, 1.0 }; + trace = Select(ata, new[] { 0.0, 0.0, 0.0 }, new[] { 0.0, 0.0, 0.0 }, new[] { -0.5, -2.0, -2.0 }, new[] { -1, -1, -1 }, new[] { 0f, 0f, 0f }, new[] { 1f, 1f, 1f }, 0.1, out release); + Assert.Equal(1, release); + Assert.Equal("Candidate 0 0.5 0 0.1\nCandidate 1 2 0 0.5", trace); + } + + [Fact] + public void ABoundIsReleasedOnlyAboveBothTheToleranceAndItsRoundingBound() + { + var ata = new[] { 1.0 }; + var lo = new[] { 0f }; + var hi = new[] { 1f }; + var x = new[] { 0.0 }; + + // violation == tolerance: not released (and not a candidate). + var trace = Select(ata, new[] { 0.0 }, x, new[] { -0.25 }, new[] { -1 }, lo, hi, 0.25, out var release); + Assert.Equal(-1, release); + Assert.Equal("", trace); + + // violation == rounding bound: a candidate that is not released. The bound here is (1 + 1) eps |atb| = 2 eps. + var noise = 2 * 2.220446049250313e-16 * 1.0; + trace = Select(ata, new[] { 1.0 }, x, new[] { -noise }, new[] { -1 }, lo, hi, 0.0, out release); + Assert.Equal(-1, release); + Assert.Equal($"Candidate 0 {BlsTrace.F(noise)} {BlsTrace.F(noise)} 0", trace); + + // One step above both: released. + var above = Math.BitIncrement(noise); + trace = Select(ata, new[] { 1.0 }, x, new[] { -above }, new[] { -1 }, lo, hi, 0.0, out release); + Assert.Equal(0, release); + Assert.Equal($"Candidate 0 {BlsTrace.F(above)} {BlsTrace.F(noise)} 0", trace); + + // A bound at hi violates with a positive gradient; a fixed column (lo == hi) is never released however steep. + trace = Select(ata, new[] { 0.0 }, new[] { 1.0 }, new[] { 4.0 }, new[] { 1 }, lo, hi, 0.5, out release); + Assert.Equal(0, release); + trace = Select(ata, new[] { 0.0 }, x, new[] { -9.0 }, new[] { -1 }, new[] { 0.5f }, new[] { 0.5f }, 0.5, out release); + Assert.Equal(-1, release); + Assert.Equal("", trace); + } + + // gamma_(n+1) = (n + 1) eps, times the magnitude of the terms of the released column's own gradient: its |atb| and its row of + // |ata x|. + [Fact] + public void TheRoundingBoundIsGammaTimesTheMagnitudeOfTheColumnsOwnTerms() + { + var ata = new[] { 2.0, 1.0, 1.0, 3.0 }; + var atb = new[] { 0.5, -4.0 }; + var x = new[] { 0.25, -2.0 }; + Assert.Equal(3 * 2.220446049250313e-16 * 10.25, BoundedLeastSquares.ReleaseNoise(ata, atb, x, 1, 2)); + Assert.Equal(3 * 2.220446049250313e-16 * 3.0, BoundedLeastSquares.ReleaseNoise(ata, atb, x, 0, 2)); + } +}