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)); + } +}