cultmath: BoundedLeastSquares primal active-set solver
Dense small bounded LS, zero allocation, caller workspace, C#-only. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens
parent
de57d8dc9c
commit
60c3d6fb0d
3 changed files with 565 additions and 0 deletions
233
packages/cultmath/src/CultMath/BoundedLeastSquares.cs
Normal file
233
packages/cultmath/src/CultMath/BoundedLeastSquares.cs
Normal file
|
|
@ -0,0 +1,233 @@
|
|||
using System.Runtime.InteropServices;
|
||||
|
||||
namespace CultMath;
|
||||
|
||||
/// <summary>Outcome of <see cref="BoundedLeastSquares.Solve"/>.</summary>
|
||||
public enum BoundedLeastSquaresStatus
|
||||
{
|
||||
/// <summary>The KKT conditions hold: every free variable has zero gradient and every bound has a gradient pointing into it.</summary>
|
||||
Converged,
|
||||
|
||||
/// <summary>The iteration cap was reached. <c>x</c> holds the last feasible iterate, which is never worse than the clamped start.</summary>
|
||||
IterationLimit,
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Dense bounded least squares: <c>min ||A x - b||^2</c> subject to <c>lo <= x <= hi</c>, for small problems
|
||||
/// (tens of rows and columns). C#-only, like the intercept helpers in <see cref="math"/>; there is no HLSL mirror.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Primal active set. Each iteration solves the free variables by Cholesky on the normal equations, takes the
|
||||
/// step to the first bound it would cross, and, once a full step lands, releases the bound with the most
|
||||
/// violated gradient sign. Free columns that are linearly dependent on earlier free columns (a zero column, a
|
||||
/// duplicate) get a zero step instead of a failed pivot, so a singular free set never throws. The normal
|
||||
/// equations square the condition number and run in float: callers that need a minimum-norm answer on
|
||||
/// rank-deficient problems add a small ridge row to <c>A</c> and <c>b</c>. The solver allocates nothing and is
|
||||
/// deterministic.
|
||||
/// </remarks>
|
||||
public static class BoundedLeastSquares
|
||||
{
|
||||
/// <summary>Default iteration cap for <see cref="Solve"/>.</summary>
|
||||
public const int DefaultMaxIterations = 100;
|
||||
|
||||
private const float PivotRelativeTolerance = 1e-6f;
|
||||
private const float KktRelativeTolerance = 1e-5f;
|
||||
|
||||
/// <summary>Floats of caller-supplied workspace needed for <paramref name="n"/> columns.</summary>
|
||||
public static int WorkspaceLength(int n) => 2 * n * n + 7 * n;
|
||||
|
||||
/// <summary>
|
||||
/// Minimises <c>||A x - b||^2</c> subject to <c>lo <= x <= hi</c>. <paramref name="x"/> is the warm start
|
||||
/// on entry (clamped to the box first) and the result on exit.
|
||||
/// </summary>
|
||||
/// <param name="m">Rows of <paramref name="a"/>.</param>
|
||||
/// <param name="n">Columns of <paramref name="a"/>, and length of <paramref name="lo"/>, <paramref name="hi"/> and <paramref name="x"/>.</param>
|
||||
/// <param name="a">Row-major <c>m * n</c> matrix.</param>
|
||||
/// <param name="b">Length <c>m</c> target.</param>
|
||||
/// <param name="lo">Lower bounds; each must be <c><= hi</c>.</param>
|
||||
/// <param name="hi">Upper bounds.</param>
|
||||
/// <param name="x">Warm start in, solution out.</param>
|
||||
/// <param name="workspace">At least <see cref="WorkspaceLength"/> floats; contents on entry are ignored.</param>
|
||||
/// <param name="iterations">Active-set iterations used.</param>
|
||||
/// <param name="maxIterations">Iteration cap.</param>
|
||||
public static BoundedLeastSquaresStatus Solve(
|
||||
int m, int n,
|
||||
ReadOnlySpan<float> a, ReadOnlySpan<float> b,
|
||||
ReadOnlySpan<float> lo, ReadOnlySpan<float> hi,
|
||||
Span<float> x, Span<float> workspace,
|
||||
out int iterations,
|
||||
int maxIterations = DefaultMaxIterations)
|
||||
{
|
||||
if (m < 0 || n < 0 || a.Length < m * n || b.Length < m || lo.Length < n || hi.Length < n
|
||||
|| x.Length < n || workspace.Length < WorkspaceLength(n))
|
||||
throw new ArgumentException("BoundedLeastSquares: span shorter than the stated dimensions.");
|
||||
|
||||
var ata = workspace.Slice(0, n * n);
|
||||
var chol = workspace.Slice(n * n, n * n);
|
||||
var off = 2 * n * n;
|
||||
var atb = workspace.Slice(off, n);
|
||||
var g = workspace.Slice(off + n, n);
|
||||
var p = workspace.Slice(off + 2 * n, n);
|
||||
var y = workspace.Slice(off + 3 * n, n);
|
||||
// One flag per column: 0 free, -1 at lo, +1 at hi. Free-list and dependency flags follow.
|
||||
var ints = MemoryMarshal.Cast<float, int>(workspace.Slice(off + 4 * n, 3 * n));
|
||||
var state = ints.Slice(0, n);
|
||||
var freeIndex = ints.Slice(n, n);
|
||||
var dependent = ints.Slice(2 * n, n);
|
||||
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
var l = lo[j];
|
||||
var h = hi[j];
|
||||
var v = x[j];
|
||||
v = !(v > l) ? l : (v < h ? v : h);
|
||||
x[j] = v;
|
||||
state[j] = v == l ? -1 : (v == h ? 1 : 0);
|
||||
}
|
||||
|
||||
var atbMax = 0f;
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
var sb = 0f;
|
||||
for (var r = 0; r < m; r++)
|
||||
sb += a[r * n + j] * b[r];
|
||||
atb[j] = sb;
|
||||
atbMax = MathF.Max(atbMax, MathF.Abs(sb));
|
||||
for (var k = 0; k <= j; k++)
|
||||
{
|
||||
var s = 0f;
|
||||
for (var r = 0; r < m; r++)
|
||||
s += a[r * n + j] * a[r * n + k];
|
||||
ata[j * n + k] = s;
|
||||
ata[k * n + j] = s;
|
||||
}
|
||||
}
|
||||
var kktTol = KktRelativeTolerance * (1f + atbMax);
|
||||
|
||||
iterations = 0;
|
||||
while (true)
|
||||
{
|
||||
Gradient(ata, atb, x, g, n);
|
||||
|
||||
var freeCount = 0;
|
||||
for (var j = 0; j < n; j++)
|
||||
if (state[j] == 0)
|
||||
freeIndex[freeCount++] = j;
|
||||
SolveFreeSet(ata, chol, g, p, y, freeIndex.Slice(0, freeCount), dependent, n);
|
||||
|
||||
var t = 1f;
|
||||
var block = -1;
|
||||
for (var i = 0; i < freeCount; i++)
|
||||
{
|
||||
var j = freeIndex[i];
|
||||
var pj = p[j];
|
||||
if (pj == 0f) continue;
|
||||
var tj = pj > 0f ? (hi[j] - x[j]) / pj : (lo[j] - x[j]) / pj;
|
||||
if (tj < t)
|
||||
{
|
||||
t = tj;
|
||||
block = j;
|
||||
}
|
||||
}
|
||||
|
||||
if (block >= 0)
|
||||
{
|
||||
if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit;
|
||||
iterations++;
|
||||
for (var i = 0; i < freeCount; i++)
|
||||
{
|
||||
var j = freeIndex[i];
|
||||
x[j] += t * p[j];
|
||||
}
|
||||
if (p[block] > 0f) { x[block] = hi[block]; state[block] = 1; }
|
||||
else { x[block] = lo[block]; state[block] = -1; }
|
||||
continue;
|
||||
}
|
||||
|
||||
for (var i = 0; i < freeCount; i++)
|
||||
{
|
||||
var j = freeIndex[i];
|
||||
var v = x[j] + p[j];
|
||||
x[j] = v < lo[j] ? lo[j] : (v > hi[j] ? hi[j] : v);
|
||||
}
|
||||
Gradient(ata, atb, x, g, n);
|
||||
|
||||
var release = -1;
|
||||
var worst = kktTol;
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
if (lo[j] == hi[j]) continue;
|
||||
var violation = state[j] < 0 ? -g[j] : (state[j] > 0 ? g[j] : 0f);
|
||||
if (violation > worst)
|
||||
{
|
||||
worst = violation;
|
||||
release = j;
|
||||
}
|
||||
}
|
||||
if (release < 0) return BoundedLeastSquaresStatus.Converged;
|
||||
if (iterations >= maxIterations) return BoundedLeastSquaresStatus.IterationLimit;
|
||||
iterations++;
|
||||
state[release] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
private static void Gradient(ReadOnlySpan<float> ata, ReadOnlySpan<float> atb, ReadOnlySpan<float> x, Span<float> g, int n)
|
||||
{
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
var s = -atb[j];
|
||||
for (var k = 0; k < n; k++)
|
||||
s += ata[j * n + k] * x[k];
|
||||
g[j] = s;
|
||||
}
|
||||
}
|
||||
|
||||
// Newton step over the free set: solves (A^T A)_FF p_F = -g_F by Cholesky. A free column whose pivot
|
||||
// collapses (dependent on earlier free columns) is dropped from the factor and gets p = 0.
|
||||
private static void SolveFreeSet(
|
||||
ReadOnlySpan<float> ata, Span<float> l, ReadOnlySpan<float> g, Span<float> p, Span<float> y,
|
||||
ReadOnlySpan<int> free, Span<int> dependent, int n)
|
||||
{
|
||||
var k = free.Length;
|
||||
for (var i = 0; i < k; i++)
|
||||
{
|
||||
var fi = free[i];
|
||||
var d = ata[fi * n + fi];
|
||||
for (var c = 0; c < i; c++)
|
||||
d -= l[i * n + c] * l[i * n + c];
|
||||
if (d <= PivotRelativeTolerance * ata[fi * n + fi] || d <= 0f)
|
||||
{
|
||||
dependent[i] = 1;
|
||||
for (var c = 0; c < i; c++) l[i * n + c] = 0f;
|
||||
l[i * n + i] = 1f;
|
||||
for (var r = i + 1; r < k; r++) l[r * n + i] = 0f;
|
||||
continue;
|
||||
}
|
||||
dependent[i] = 0;
|
||||
var lii = MathF.Sqrt(d);
|
||||
l[i * n + i] = lii;
|
||||
for (var r = i + 1; r < k; r++)
|
||||
{
|
||||
var s = ata[free[r] * n + fi];
|
||||
for (var c = 0; c < i; c++)
|
||||
s -= l[r * n + c] * l[i * n + c];
|
||||
l[r * n + i] = s / lii;
|
||||
}
|
||||
}
|
||||
|
||||
for (var i = 0; i < k; i++)
|
||||
{
|
||||
var s = dependent[i] != 0 ? 0f : -g[free[i]];
|
||||
for (var c = 0; c < i; c++)
|
||||
s -= l[i * n + c] * y[c];
|
||||
y[i] = s / l[i * n + i];
|
||||
}
|
||||
for (var i = k - 1; i >= 0; i--)
|
||||
{
|
||||
var s = y[i];
|
||||
for (var r = i + 1; r < k; r++)
|
||||
s -= l[r * n + i] * p[free[r]];
|
||||
p[free[i]] = s / l[i * n + i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,2 @@
|
|||
fileFormatVersion: 2
|
||||
guid: 4ba991e1b6c54e099bc7ea16541734c5
|
||||
|
|
@ -0,0 +1,330 @@
|
|||
using System;
|
||||
using CultMath;
|
||||
using Xunit;
|
||||
|
||||
namespace CultMath.Tests;
|
||||
|
||||
// Fixtures are deliberately non-degenerate: A is non-square and non-identity, and b is outside A's range,
|
||||
// so the unconstrained optimum has a real residual and the box actually bites.
|
||||
public sealed class BoundedLeastSquaresTests
|
||||
{
|
||||
private static double Cost(int m, int n, float[] a, float[] b, float[] x)
|
||||
{
|
||||
var c = 0.0;
|
||||
for (var r = 0; r < m; r++)
|
||||
{
|
||||
var s = -(double)b[r];
|
||||
for (var j = 0; j < n; j++) s += (double)a[r * n + j] * x[j];
|
||||
c += s * s;
|
||||
}
|
||||
return c;
|
||||
}
|
||||
|
||||
// Gradient of 0.5*||Ax-b||^2, in double.
|
||||
private static double[] Gradient(int m, int n, float[] a, float[] b, float[] x)
|
||||
{
|
||||
var g = new double[n];
|
||||
for (var r = 0; r < m; r++)
|
||||
{
|
||||
var s = -(double)b[r];
|
||||
for (var j = 0; j < n; j++) s += (double)a[r * n + j] * x[j];
|
||||
for (var j = 0; j < n; j++) g[j] += a[r * n + j] * s;
|
||||
}
|
||||
return g;
|
||||
}
|
||||
|
||||
private static double AtbScale(int m, int n, float[] a, float[] b)
|
||||
{
|
||||
var mx = 0.0;
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
var s = 0.0;
|
||||
for (var r = 0; r < m; r++) s += (double)a[r * n + j] * b[r];
|
||||
mx = Math.Max(mx, Math.Abs(s));
|
||||
}
|
||||
return 1 + mx;
|
||||
}
|
||||
|
||||
// Asserts feasibility and the KKT sign rule: free columns have ~zero gradient, columns at lo have
|
||||
// gradient >= 0, columns at hi have gradient <= 0.
|
||||
private static void AssertKkt(int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[] x)
|
||||
{
|
||||
var g = Gradient(m, n, a, b, x);
|
||||
var tol = 1e-3 * AtbScale(m, n, a, b);
|
||||
for (var j = 0; j < n; j++)
|
||||
{
|
||||
Assert.True(x[j] >= lo[j] && x[j] <= hi[j], $"column {j} outside its box: {x[j]}");
|
||||
if (lo[j] == hi[j]) continue;
|
||||
if (x[j] == lo[j]) Assert.True(g[j] >= -tol, $"column {j} at lo with gradient {g[j]}");
|
||||
else if (x[j] == hi[j]) Assert.True(g[j] <= tol, $"column {j} at hi with gradient {g[j]}");
|
||||
else Assert.True(Math.Abs(g[j]) <= tol, $"column {j} free with gradient {g[j]}");
|
||||
}
|
||||
}
|
||||
|
||||
// Projected gradient with step 1/L, L bounded by the Frobenius norm squared of A.
|
||||
private static double ReferenceCost(int m, int n, float[] a, float[] b, float[] lo, float[] hi)
|
||||
{
|
||||
var frob = 0.0;
|
||||
foreach (var v in a) frob += (double)v * v;
|
||||
var step = 1.0 / Math.Max(frob, 1e-12);
|
||||
var x = new float[n];
|
||||
for (var j = 0; j < n; j++) x[j] = Math.Clamp(0f, lo[j], hi[j]);
|
||||
var xd = new double[n];
|
||||
for (var j = 0; j < n; j++) xd[j] = x[j];
|
||||
for (var it = 0; it < 200000; it++)
|
||||
{
|
||||
for (var j = 0; j < n; j++) x[j] = (float)xd[j];
|
||||
var g = Gradient(m, n, a, b, x);
|
||||
for (var j = 0; j < n; j++) xd[j] = Math.Clamp(xd[j] - step * g[j], lo[j], hi[j]);
|
||||
}
|
||||
for (var j = 0; j < n; j++) x[j] = (float)xd[j];
|
||||
return Cost(m, n, a, b, x);
|
||||
}
|
||||
|
||||
private static BoundedLeastSquaresStatus Solve(int m, int n, float[] a, float[] b, float[] lo, float[] hi, float[] x,
|
||||
out int iterations, int maxIterations = BoundedLeastSquares.DefaultMaxIterations)
|
||||
=> BoundedLeastSquares.Solve(m, n, a, b, lo, hi, x, new float[BoundedLeastSquares.WorkspaceLength(n)], out iterations, maxIterations);
|
||||
|
||||
private static float[] Fill(float value, int n)
|
||||
{
|
||||
var v = new float[n];
|
||||
Array.Fill(v, value);
|
||||
return v;
|
||||
}
|
||||
|
||||
private static readonly float[] A54 =
|
||||
{
|
||||
1.0f, 0.3f, -0.2f,
|
||||
0.2f, 1.1f, 0.4f,
|
||||
-0.5f, 0.6f, 0.9f,
|
||||
0.7f, -0.4f, 0.5f,
|
||||
0.1f, 0.8f, -0.6f,
|
||||
};
|
||||
|
||||
private static readonly float[] B5 = { 0.9f, -0.4f, 1.3f, 0.2f, -0.7f };
|
||||
|
||||
[Fact]
|
||||
public void InteriorOptimumEqualsNormalEquations()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, A54, B5, Fill(-100f, n), Fill(100f, n), x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
|
||||
// Independent solve of (A^T A) x = A^T b by Gauss-Jordan in double.
|
||||
var mat = new double[n, n + 1];
|
||||
for (var i = 0; i < n; i++)
|
||||
{
|
||||
for (var j = 0; j < n; j++)
|
||||
for (var r = 0; r < m; r++) mat[i, j] += (double)A54[r * n + i] * A54[r * n + j];
|
||||
for (var r = 0; r < m; r++) mat[i, n] += (double)A54[r * n + i] * B5[r];
|
||||
}
|
||||
for (var c = 0; c < n; c++)
|
||||
{
|
||||
var piv = mat[c, c];
|
||||
for (var j = 0; j <= n; j++) mat[c, j] /= piv;
|
||||
for (var i = 0; i < n; i++)
|
||||
{
|
||||
if (i == c) continue;
|
||||
var f = mat[i, c];
|
||||
for (var j = 0; j <= n; j++) mat[i, j] -= f * mat[c, j];
|
||||
}
|
||||
}
|
||||
for (var i = 0; i < n; i++) Assert.Equal(mat[i, n], x[i], 1e-4);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void KktHoldsOnRandomBoxedProblems()
|
||||
{
|
||||
var anyActive = false;
|
||||
for (uint seed = 1; seed <= 150; seed++)
|
||||
{
|
||||
var rng = new CultMath.Random(seed * 7919u);
|
||||
var n = rng.NextInt(2, 17);
|
||||
var m = rng.NextInt(n, 25);
|
||||
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++)
|
||||
{
|
||||
switch (rng.NextInt(3))
|
||||
{
|
||||
case 0: lo[j] = -1f; hi[j] = 1f; break;
|
||||
case 1: lo[j] = 0f; hi[j] = 1f; break;
|
||||
default: lo[j] = -0.25f; hi[j] = 0.1f; break;
|
||||
}
|
||||
}
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, a, b, lo, hi, x, out _);
|
||||
Assert.True(status == BoundedLeastSquaresStatus.Converged, $"seed {seed} ({m}x{n}) hit the cap");
|
||||
AssertKkt(m, n, a, b, lo, hi, x);
|
||||
for (var j = 0; j < n; j++) anyActive |= x[j] == lo[j] || x[j] == hi[j];
|
||||
|
||||
var cs = Cost(m, n, a, b, x);
|
||||
var cr = ReferenceCost(m, n, a, b, lo, hi);
|
||||
Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6, $"seed {seed} ({m}x{n}): solver cost {cs}, reference {cr}");
|
||||
}
|
||||
Assert.True(anyActive);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void AllBoundsActive()
|
||||
{
|
||||
const int m = 4, n = 3;
|
||||
var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f };
|
||||
var b = new[] { 10f, -10f, 10f, 1f };
|
||||
var lo = Fill(-1f, n);
|
||||
var hi = Fill(1f, n);
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, a, b, lo, hi, x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
Assert.Equal(new[] { 1f, -1f, 1f }, x);
|
||||
AssertKkt(m, n, a, b, lo, hi, x);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void ZeroColumn()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var a = (float[])A54.Clone();
|
||||
for (var r = 0; r < m; r++) a[r * n + 1] = 0f;
|
||||
var lo = Fill(-1f, n);
|
||||
var hi = Fill(1f, n);
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, a, B5, lo, hi, x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
AssertKkt(m, n, a, B5, lo, hi, x);
|
||||
var cs = Cost(m, n, a, B5, x);
|
||||
var cr = ReferenceCost(m, n, a, B5, lo, hi);
|
||||
Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void DuplicateColumns()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var a = (float[])A54.Clone();
|
||||
for (var r = 0; r < m; r++) a[r * n + 2] = a[r * n + 0];
|
||||
var lo = Fill(-1f, n);
|
||||
var hi = new[] { 0.3f, 1f, 0.3f };
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, a, B5, lo, hi, x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
AssertKkt(m, n, a, B5, lo, hi, x);
|
||||
var cs = Cost(m, n, a, B5, x);
|
||||
var cr = ReferenceCost(m, n, a, B5, lo, hi);
|
||||
Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void WarmStartMatchesColdStart()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var lo = new[] { -0.2f, -1f, 0f };
|
||||
var hi = new[] { 0.5f, 0.3f, 1f };
|
||||
var cold = new float[n];
|
||||
Solve(m, n, A54, B5, lo, hi, cold, out _);
|
||||
AssertKkt(m, n, A54, B5, lo, hi, cold);
|
||||
|
||||
var starts = new[] { new[] { 0.5f, -1f, 1f }, new[] { 9f, 9f, -9f }, new[] { 0.1f, 0.1f, 0.5f }, (float[])cold.Clone() };
|
||||
foreach (var start in starts)
|
||||
{
|
||||
var status = Solve(m, n, A54, B5, lo, hi, start, out var iterations);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
for (var j = 0; j < n; j++) Assert.Equal(cold[j], start[j], 1e-4);
|
||||
if (ReferenceEquals(start, starts[^1])) Assert.Equal(0, iterations);
|
||||
}
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void AsymmetricBounds()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var lo = new[] { -1f, 0f, -1f };
|
||||
var hi = new[] { 1f, 1f, 1f };
|
||||
// Column 1 wants to go negative (b has a strong component against it) but its floor is 0.
|
||||
var b = new[] { 0.2f, -2.5f, 0.4f, -1.5f, -2.0f };
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, A54, b, lo, hi, x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
AssertKkt(m, n, A54, b, lo, hi, x);
|
||||
Assert.Equal(0f, x[1]);
|
||||
Assert.True(x[0] < 0f || x[2] < 0f, "a symmetric-bounds column should use its negative side");
|
||||
var cs = Cost(m, n, A54, b, x);
|
||||
var cr = ReferenceCost(m, n, A54, b, lo, hi);
|
||||
Assert.True(Math.Abs(cs - cr) <= 1e-4 * cr + 1e-6);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void IterationCapReturnsFeasibleIterateAndLimitStatus()
|
||||
{
|
||||
const int m = 4, n = 3;
|
||||
var a = new[] { 1.0f, 0.2f, 0.1f, 0.1f, 1.0f, 0.3f, -0.2f, 0.2f, 1.0f, 0.5f, 0.5f, 0.4f };
|
||||
var b = new[] { 10f, -10f, 10f, 1f };
|
||||
var lo = Fill(-1f, n);
|
||||
var hi = Fill(1f, n);
|
||||
var x = new float[n];
|
||||
var start = Cost(m, n, a, b, x);
|
||||
var status = Solve(m, n, a, b, lo, hi, x, out var iterations, maxIterations: 1);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.IterationLimit, status);
|
||||
Assert.Equal(1, iterations);
|
||||
for (var j = 0; j < n; j++) Assert.InRange(x[j], lo[j], hi[j]);
|
||||
Assert.True(Cost(m, n, a, b, x) < start);
|
||||
|
||||
var zero = new float[n];
|
||||
Assert.Equal(BoundedLeastSquaresStatus.IterationLimit, Solve(m, n, a, b, lo, hi, zero, out iterations, maxIterations: 0));
|
||||
Assert.Equal(0, iterations);
|
||||
Assert.Equal(new float[n], zero);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void StartOutsideTheBoxIsClampedFirst()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var lo = Fill(-0.1f, n);
|
||||
var hi = Fill(0.1f, n);
|
||||
var x = new[] { 50f, -50f, float.NaN };
|
||||
Solve(m, n, A54, B5, lo, hi, x, out _, maxIterations: 0);
|
||||
Assert.Equal(new[] { 0.1f, -0.1f, -0.1f }, x);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void FixedColumnStaysPut()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var lo = new[] { -1f, 0.25f, -1f };
|
||||
var hi = new[] { 1f, 0.25f, 1f };
|
||||
var x = new float[n];
|
||||
var status = Solve(m, n, A54, B5, lo, hi, x, out _);
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
Assert.Equal(0.25f, x[1]);
|
||||
AssertKkt(m, n, A54, B5, lo, hi, x);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void WarmCallAllocatesNothing()
|
||||
{
|
||||
const int m = 5, n = 3;
|
||||
var lo = new[] { -0.2f, -1f, 0f };
|
||||
var hi = new[] { 0.5f, 0.3f, 1f };
|
||||
var x = new float[n];
|
||||
var workspace = new float[BoundedLeastSquares.WorkspaceLength(n)];
|
||||
BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, x, workspace, out _);
|
||||
|
||||
var before = GC.GetAllocatedBytesForCurrentThread();
|
||||
var status = BoundedLeastSquares.Solve(m, n, A54, B5, lo, hi, x, workspace, out _);
|
||||
var after = GC.GetAllocatedBytesForCurrentThread();
|
||||
Assert.Equal(BoundedLeastSquaresStatus.Converged, status);
|
||||
Assert.Equal(0L, after - before);
|
||||
}
|
||||
|
||||
[Fact]
|
||||
public void ShortSpansThrowArgumentException()
|
||||
{
|
||||
var x = new float[3];
|
||||
Assert.Throws<ArgumentException>(() => BoundedLeastSquares.Solve(5, 3, A54, B5, Fill(-1f, 3), Fill(1f, 3), x, new float[4], out _));
|
||||
}
|
||||
}
|
||||
Loading…
Reference in a new issue