GameCult/CultLib
Watch
9
0
Fork
You've already forked CultLib
 0
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions

cultmath: near-collinear test uses generators that reproduce the gap on the old solver

Browse source 
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-09-29 03:44:57 +02:00
commit 8b1a1ffde0
1 changed files with 25 additions and 19 deletions
Download patch file Download diff file Expand all files Collapse all files

44
packages/cultmath/tests/CultMath.Tests/BoundedLeastSquaresTests.cs
View file

@ -657,32 +657,38 @@ public sealed class BoundedLeastSquaresTests
return best;
}
// m < n with columns that are nearly parallel (condition number of A around 1e2), b outside the reachable set.
[Fact]
public void NearCollinearColumnsReachTheExactOptimum()
{
const int m = 2, n = 7;
var worst = 0.0;
for (uint seed = 1; seed <= 60; seed++)
for (var family = 0; family < 4; family++)
{
var rng = new CultMath.Random(seed * 15485863u);
var a = new float[m * n];
for (var j = 0; j < n; j++)
for (uint seed = 1; seed <= 120; seed++)
{
a[j] = rng.NextFloat(0.5f, 1f);
a[n + j] = 0.01f * rng.NextFloat(-1f, 1f);
var rng = new CultMath.Random(seed * 15485863u + (uint)family * 977u);
var a = new float[m * n];
for (var j = 0; j < n; j++)
{
switch (family)
{
case 0: a[j] = rng.NextFloat(0.5f, 1f); a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); break;
case 1: a[j] = 1f + 0.02f * rng.NextFloat(-1f, 1f); a[n + j] = 0.01f * rng.NextFloat(-1f, 1f); break;
case 2: a[j] = rng.NextFloat(-1f, 1f); a[n + j] = 0.01f * a[j] + 0.01f * rng.NextFloat(-1f, 1f); break;
default: a[j] = rng.NextFloat(0.5f, 1f); a[n + j] = a[j] * 0.9f + 0.01f * rng.NextFloat(-1f, 1f); break;
}
}
var b = new[] { rng.NextFloat(-8f, 8f), rng.NextFloat(-8f, 8f) };
var lo = new float[n];
var hi = new float[n];
for (var j = 0; j < n; j++) { lo[j] = rng.NextFloat() < 0.5f ? -1f : 0f; hi[j] = 1f; }
var x = new float[n];
var status = Solve(m, n, a, b, lo, hi, x, out _);
Assert.True(status == BoundedLeastSquaresStatus.Converged, $"family {family} seed {seed}: {status}");
var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi);
var allowed = 1e-5 * (b[0] * b[0] + b[1] * b[1]);
Assert.True(gap <= allowed, $"family {family} seed {seed}: cost {gap} above the exact optimum (allowed {allowed})");
}
var b = new[] { 6f * (rng.NextFloat() < 0.5f ? -1f : 1f), rng.NextFloat(0.5f, 2f) };
var lo = new float[n];
var hi = new float[n];
for (var j = 0; j < n; j++) { lo[j] = -1f; hi[j] = seed % 2 == 0 ? 1f : 0.5f; }
var x = new float[n];
var status = Solve(m, n, a, b, lo, hi, x, out _);
Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed}: {status}");
AssertKkt(m, n, a, b, lo, hi, x);
var gap = Cost(m, n, a, b, x) - ExactBoxedCost(m, n, a, b, lo, hi);
var bound = 1e-5 * (b[0] * b[0] + b[1] * b[1]);
worst = Math.Max(worst, gap / bound);
Assert.True(gap <= bound, $"seed {seed}: cost {gap} above the exact optimum (allowed {bound})");
}
}
}