651 lines
24 KiB
Go
651 lines
24 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package optim
|
||
|
||
import (
|
||
"math"
|
||
"strings"
|
||
"testing"
|
||
|
||
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
||
)
|
||
|
||
// Line-search and convergence-exit pins for the optimisers. Each test
|
||
// names the defect it pins; the numeric fixtures are hand-derived
|
||
// optima.
|
||
|
||
// stiffQuadratic is f = k·(x−3)², whose minimum is x = 3 with value 0.
|
||
func stiffQuadratic(k float64) func(*core.Array) (float64, error) {
|
||
return func(p *core.Array) (float64, error) {
|
||
d := p.FloatAt(0) - 3
|
||
return k * d * d, nil
|
||
}
|
||
}
|
||
|
||
// TestLBFGSLineSearchStiffQuadraticConverges pins the line-search
|
||
// budget: the step that reduces a quadratic is ≈ 1/L for a curvature
|
||
// L, so with 20 halvings from a unit step the first trial is never
|
||
// acceptable once L ≳ 1e6 and L-BFGS used to return the start point as
|
||
// a converged answer (x = 0, value 9k).
|
||
func TestLBFGSLineSearchStiffQuadraticConverges(t *testing.T) {
|
||
for _, k := range []float64{1e6, 1e8, 1e12} {
|
||
point, value, err := MinimiseLBFGS(stiffQuadratic(k), nil, mustFloats(t, []float64{0}), LBFGSOptions{})
|
||
if err != nil {
|
||
t.Errorf("k=%g: MinimiseLBFGS: %v", k, err)
|
||
continue
|
||
}
|
||
if got := point.FloatAt(0); math.Abs(got-3) > 1e-6 {
|
||
t.Errorf("k=%g: x = %.12g, want 3", k, got)
|
||
}
|
||
if value > 1e-6 {
|
||
t.Errorf("k=%g: value = %.12g, want 0", k, value)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestLBFGSLineSearchMixedUnitsConverges is the same failure through an
|
||
// ordinary two-parameter fit: the y coordinate's curvature sets the
|
||
// first step and takes x down with it when the step is capped at a
|
||
// unit. The optimum of (x−3)² + K·(y−5)² is (3, 5) with value 0.
|
||
func TestLBFGSLineSearchMixedUnitsConverges(t *testing.T) {
|
||
for _, k := range []float64{1e6, 1e8} {
|
||
f := func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-3, p.FloatAt(1)-5
|
||
return dx*dx + k*dy*dy, nil
|
||
}
|
||
grad := func(p *core.Array) (*core.Array, error) {
|
||
return core.FromFloats([]float64{2 * (p.FloatAt(0) - 3), 2 * k * (p.FloatAt(1) - 5)}, 2)
|
||
}
|
||
for _, g := range []struct {
|
||
name string
|
||
fn func(*core.Array) (*core.Array, error)
|
||
}{{"finite differences", nil}, {"analytic gradient", grad}} {
|
||
point, value, err := MinimiseLBFGS(f, g.fn, mustFloats(t, []float64{0, 0}, 2), LBFGSOptions{})
|
||
if err != nil {
|
||
t.Errorf("K=%g (%s): %v", k, g.name, err)
|
||
continue
|
||
}
|
||
if math.Abs(point.FloatAt(0)-3) > 1e-6 || math.Abs(point.FloatAt(1)-5) > 1e-6 {
|
||
t.Errorf("K=%g (%s): point = (%.12g, %.12g), want (3, 5)",
|
||
k, g.name, point.FloatAt(0), point.FloatAt(1))
|
||
}
|
||
if value > 1e-6 {
|
||
t.Errorf("K=%g (%s): value = %.12g, want 0", k, g.name, value)
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestLBFGSLineSearchStallIsAnError pins the silence: an exhausted
|
||
// line search used to break out of the iteration and hand the start
|
||
// point back as converged. The objective here is flat below zero and
|
||
// jumps at it, so the finite-difference stencil just below the jump
|
||
// reports a large gradient while no step along it can reduce the
|
||
// value: the search stalls and must say so.
|
||
func TestLBFGSLineSearchStallIsAnError(t *testing.T) {
|
||
f := func(p *core.Array) (float64, error) {
|
||
if p.FloatAt(0) >= 0 {
|
||
return 1, nil
|
||
}
|
||
return 0, nil
|
||
}
|
||
start := mustFloats(t, []float64{-1e-9})
|
||
point, value, err := MinimiseLBFGS(f, nil, start, LBFGSOptions{})
|
||
if err == nil {
|
||
t.Fatalf("a stalled line search was reported as convergence: point = %v, value = %g",
|
||
floatsOf(point), value)
|
||
}
|
||
if !strings.Contains(err.Error(), "line search") {
|
||
t.Fatalf("error = %v, want a line-search refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestLBFGSBoxOptimaFixtures checks the box-constrained answers against
|
||
// hand-derived optima: (x−3)²+(y+2)² on −1 ≤ x ≤ 1, y free reaches the
|
||
// upper wall at (1, −2) with value 4; (x+5)²+(y−5)² on 0 ≤ x, y ≤ 2
|
||
// pins both coordinates at (0, 2) with value 34; (x+y−3)²+x² on
|
||
// x, y ≥ 1 bottoms out at the corner (1, 2) with value 1.
|
||
func TestLBFGSBoxOptimaFixtures(t *testing.T) {
|
||
inf := math.Inf(1)
|
||
cases := []struct {
|
||
name string
|
||
f func(*core.Array) (float64, error)
|
||
lower []float64
|
||
upper []float64
|
||
wantX []float64
|
||
wantValue float64
|
||
}{
|
||
{
|
||
name: "upper wall, y free",
|
||
f: func(p *core.Array) (float64, error) {
|
||
return (p.FloatAt(0)-3)*(p.FloatAt(0)-3) + (p.FloatAt(1)+2)*(p.FloatAt(1)+2), nil
|
||
},
|
||
lower: []float64{-1, -inf}, upper: []float64{1, inf},
|
||
wantX: []float64{1, -2}, wantValue: 4,
|
||
},
|
||
{
|
||
name: "both coordinates pinned",
|
||
f: func(p *core.Array) (float64, error) {
|
||
return (p.FloatAt(0)+5)*(p.FloatAt(0)+5) + (p.FloatAt(1)-5)*(p.FloatAt(1)-5), nil
|
||
},
|
||
lower: []float64{0, 0}, upper: []float64{inf, 2},
|
||
wantX: []float64{0, 2}, wantValue: 34,
|
||
},
|
||
{
|
||
name: "coupled bowl at the corner",
|
||
f: func(p *core.Array) (float64, error) {
|
||
x, y := p.FloatAt(0), p.FloatAt(1)
|
||
return (x+y-3)*(x+y-3) + x*x, nil
|
||
},
|
||
lower: []float64{1, 1}, upper: []float64{inf, inf},
|
||
wantX: []float64{1, 2}, wantValue: 1,
|
||
},
|
||
}
|
||
for _, c := range cases {
|
||
point, value, err := MinimiseLBFGS(c.f, nil, mustFloats(t, []float64{0, 0}, 2),
|
||
LBFGSOptions{Lower: c.lower, Upper: c.upper})
|
||
if err != nil {
|
||
t.Errorf("%s: %v", c.name, err)
|
||
continue
|
||
}
|
||
for i, want := range c.wantX {
|
||
if math.Abs(point.FloatAt(i)-want) > 1e-4 {
|
||
t.Errorf("%s: x[%d] = %.12g, want %.12g", c.name, i, point.FloatAt(i), want)
|
||
}
|
||
}
|
||
if math.Abs(value-c.wantValue) > 1e-6 {
|
||
t.Errorf("%s: value = %.12g, want %.12g", c.name, value, c.wantValue)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestLBFGSProjectionOntoBindingWall pins the box projection with a
|
||
// box that actually binds: (x−3)² on x ≤ 1 has its constrained minimum
|
||
// on the wall at x = 1 with value 4, and a start above the wall is
|
||
// projected onto it rather than refused.
|
||
func TestLBFGSProjectionOntoBindingWall(t *testing.T) {
|
||
f := func(p *core.Array) (float64, error) {
|
||
d := p.FloatAt(0) - 3
|
||
return d * d, nil
|
||
}
|
||
for _, start := range []float64{0, 5, -10} {
|
||
point, value, err := MinimiseLBFGS(f, nil, mustFloats(t, []float64{start}),
|
||
LBFGSOptions{Upper: []float64{1}})
|
||
if err != nil {
|
||
t.Fatalf("start %g: %v", start, err)
|
||
}
|
||
if math.Abs(point.FloatAt(0)-1) > 1e-6 {
|
||
t.Errorf("start %g: x = %.12g, want 1 on the wall", start, point.FloatAt(0))
|
||
}
|
||
if math.Abs(value-4) > 1e-6 {
|
||
t.Errorf("start %g: value = %.12g, want 4", start, value)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestMinimiseLevelSetStallIsNotConvergence pins the level-set trap: the
|
||
// value spread was the only convergence test, so a simplex whose
|
||
// vertices happened to lie on one level set stopped while spanning the
|
||
// space.
|
||
// For (x−1)² + (y−2)² from (0, 0) all three vertices of the stalled
|
||
// simplex sit on the circle of radius √0.5 about (1, 2), so the
|
||
// spread is zero and the answer used to be (1.5, 1.5) with value 0.5.
|
||
func TestMinimiseLevelSetStallIsNotConvergence(t *testing.T) {
|
||
f := func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-1, p.FloatAt(1)-2
|
||
return dx*dx + dy*dy, nil
|
||
}
|
||
start := mustFloats(t, []float64{0, 0}, 2)
|
||
point, value, err := Minimise(f, start, MinimiseOptions{})
|
||
if err != nil {
|
||
t.Fatalf("Minimise: %v", err)
|
||
}
|
||
if math.Abs(point.FloatAt(0)-1) > 1e-4 || math.Abs(point.FloatAt(1)-2) > 1e-4 {
|
||
t.Fatalf("point = (%v), want (1, 2)", floatsOf(point))
|
||
}
|
||
if value > 1e-8 {
|
||
t.Fatalf("value = %g, want 0", value)
|
||
}
|
||
// Neither a larger budget nor a tighter tolerance may rescue the
|
||
// stall: the loop exits at the top-of-loop test either way.
|
||
for _, opts := range []MinimiseOptions{
|
||
{MaxIterations: 20000},
|
||
{Tolerance: 1e-20},
|
||
} {
|
||
point, value, err := Minimise(f, start, opts)
|
||
if err != nil {
|
||
t.Fatalf("%+v: %v", opts, err)
|
||
}
|
||
if value > 1e-8 {
|
||
t.Fatalf("%+v: value = %g at (%v), want 0 at (1, 2)", opts, value, floatsOf(point))
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestMinimiseHitRateOnShiftedBowls sweeps starts over a grid: every
|
||
// one of them must find the minimum of the axis-aligned bowl, the
|
||
// rotated bowl and the 3-D sphere. The level-set stall used to return
|
||
// a non-minimal point for 5 of 49 axis-bowl starts and 1 of 49 rotated
|
||
// starts with the default options.
|
||
func TestMinimiseHitRateOnShiftedBowls(t *testing.T) {
|
||
bowls := []struct {
|
||
name string
|
||
f func(*core.Array) (float64, error)
|
||
dim int
|
||
}{
|
||
{"axis bowl", func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-1, p.FloatAt(1)-2
|
||
return dx*dx + dy*dy, nil
|
||
}, 2},
|
||
{"rotated bowl", func(p *core.Array) (float64, error) {
|
||
u := (p.FloatAt(0) - 1) + (p.FloatAt(1) - 2)
|
||
v := (p.FloatAt(0) - 1) - (p.FloatAt(1) - 2)
|
||
return u*u + 3*v*v, nil
|
||
}, 2},
|
||
{"3-D sphere", func(p *core.Array) (float64, error) {
|
||
s := 0.0
|
||
for i, c := range []float64{1, 2, 3} {
|
||
d := p.FloatAt(i) - c
|
||
s += d * d
|
||
}
|
||
return s, nil
|
||
}, 3},
|
||
}
|
||
for _, b := range bowls {
|
||
for x := -3.0; x <= 3; x++ {
|
||
for y := -3.0; y <= 3; y++ {
|
||
start := []float64{x, y}
|
||
if b.dim == 3 {
|
||
start = []float64{x, y, x - y}
|
||
}
|
||
point, value, err := Minimise(b.f, mustFloats(t, start, b.dim), MinimiseOptions{})
|
||
if err != nil {
|
||
t.Fatalf("%s from %v: %v", b.name, start, err)
|
||
}
|
||
if value > 1e-6 {
|
||
t.Errorf("%s from %v: value = %g at (%v), want 0", b.name, start, value, floatsOf(point))
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestMinimiseToleranceScaleIsDocumented pins the documented absolute
|
||
// tolerance: an objective whose values are ~1e-14 already counts as
|
||
// flat (the default spread test is 1e-10·max(1, |f|)), so Minimise
|
||
// reports the start point and a nil error, and rescaling the objective
|
||
// to O(1), the documented remedy, resolves the minimum (1, 2).
|
||
func TestMinimiseToleranceScaleIsDocumented(t *testing.T) {
|
||
const scale = 1e-14
|
||
tiny := func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-1, p.FloatAt(1)-2
|
||
return scale * (dx*dx + dy*dy), nil
|
||
}
|
||
start := mustFloats(t, []float64{0, 0}, 2)
|
||
point, value, err := Minimise(tiny, start, MinimiseOptions{})
|
||
if err != nil {
|
||
t.Fatalf("Minimise: %v", err)
|
||
}
|
||
// The start's own value is 5e-14; the run stops on the value spread
|
||
// long before the minimum is reached, which the documentation now
|
||
// warns about.
|
||
if math.Abs(point.FloatAt(0)-1) <= 0.1 || math.Abs(point.FloatAt(1)-2) <= 0.1 {
|
||
t.Fatalf("tiny objective: point = (%v), the documented absolute tolerance claims (1, 2) is left unfound",
|
||
floatsOf(point))
|
||
}
|
||
if value > 5e-14 {
|
||
t.Fatalf("tiny objective: value = %g, want a value no larger than the start's 5e-14", value)
|
||
}
|
||
unit := func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-1, p.FloatAt(1)-2
|
||
return dx*dx + dy*dy, nil
|
||
}
|
||
point, value, err = Minimise(unit, start, MinimiseOptions{})
|
||
if err != nil {
|
||
t.Fatalf("Minimise rescaled: %v", err)
|
||
}
|
||
if math.Abs(point.FloatAt(0)-1) > 1e-4 || math.Abs(point.FloatAt(1)-2) > 1e-4 || value > 1e-8 {
|
||
t.Fatalf("rescaled objective: point = (%v), value = %g, want (1, 2) and 0", floatsOf(point), value)
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedComplexMatrixRefused pins the dtype guard on
|
||
// the constraint matrix: a complex A used to reach FloatAt and panic
|
||
// with an indexing error instead of the family's dtype refusal.
|
||
func TestMinimiseConstrainedComplexMatrixRefused(t *testing.T) {
|
||
A, err := core.FromComplexes([]complex128{1, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
start := mustFloats(t, []float64{5, -3}, 2)
|
||
if _, _, err := MinimiseConstrained(bowlAt(1, 2), nil, start,
|
||
LinearConstraints{A: A, Lower: []float64{3}, Upper: []float64{3}}, LBFGSOptions{}); err == nil {
|
||
t.Fatal("expected a complex constraint matrix to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("error = %v, want a complex-input refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedComplexGradientRefused pins the callback guard
|
||
// in the constrained wrapper: a complex gradient used to be sliced with
|
||
// RawFloats (nil for a complex payload) and panic.
|
||
func TestMinimiseConstrainedComplexGradientRefused(t *testing.T) {
|
||
grad := func(*core.Array) (*core.Array, error) {
|
||
g, err := core.FromComplexes([]complex128{1, 1}, 2)
|
||
return g, err
|
||
}
|
||
A, err := core.FromFloats([]float64{1, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
start := mustFloats(t, []float64{5, -3}, 2)
|
||
if _, _, err := MinimiseConstrained(bowlAt(1, 2), grad, start,
|
||
LinearConstraints{A: A, Lower: []float64{3}, Upper: []float64{3}}, LBFGSOptions{}); err == nil {
|
||
t.Fatal("expected a complex gradient to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("error = %v, want a complex-input refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedShortGradientRefused pins the length contract
|
||
// in the constrained wrapper: a gradient one element short used to
|
||
// panic in the slice before MinimiseLBFGS could report it.
|
||
func TestMinimiseConstrainedShortGradientRefused(t *testing.T) {
|
||
grad := func(*core.Array) (*core.Array, error) { return mustFloats(t, []float64{1}), nil }
|
||
A, err := core.FromFloats([]float64{1, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
start := mustFloats(t, []float64{5, -3}, 2)
|
||
if _, _, err := MinimiseConstrained(bowlAt(1, 2), grad, start,
|
||
LinearConstraints{A: A, Lower: []float64{3}, Upper: []float64{3}}, LBFGSOptions{}); err == nil {
|
||
t.Fatal("expected a short gradient to be refused")
|
||
} else if !strings.Contains(err.Error(), "gradient callback") {
|
||
t.Fatalf("error = %v, want a callback-length refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestLBFGSComplexGradientRefused pins the dtype guard on the L-BFGS
|
||
// gradient callback, which used to dereference the nil int payload of a
|
||
// complex array.
|
||
func TestLBFGSComplexGradientRefused(t *testing.T) {
|
||
grad := func(*core.Array) (*core.Array, error) {
|
||
g, err := core.FromComplexes([]complex128{1, 1}, 2)
|
||
return g, err
|
||
}
|
||
if _, _, err := MinimiseLBFGS(bowlAt(1, 2), grad, mustFloats(t, []float64{5, -3}, 2), LBFGSOptions{}); err == nil {
|
||
t.Fatal("expected a complex gradient to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("error = %v, want a complex-input refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestLevenbergMarquardtComplexPayloadRefused pins both callback guards
|
||
// of the LM fitter: a complex residual and a complex analytic Jacobian
|
||
// used to panic in FloatAt.
|
||
func TestLevenbergMarquardtComplexPayloadRefused(t *testing.T) {
|
||
p0 := mustFloats(t, []float64{0, 0}, 2)
|
||
complexResidual := func(*core.Array) (*core.Array, error) {
|
||
r, err := core.FromComplexes([]complex128{1, 2, 3, 4}, 4)
|
||
return r, err
|
||
}
|
||
if _, _, err := LevenbergMarquardt(complexResidual, p0, LMOptions{}); err == nil {
|
||
t.Fatal("expected a complex residual to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("residual: error = %v, want a complex-input refusal", err)
|
||
}
|
||
residual := func(p *core.Array) (*core.Array, error) {
|
||
return core.FromFloats([]float64{p.FloatAt(0) - 1, p.FloatAt(1) - 2}, 2)
|
||
}
|
||
complexJacobian := func(*core.Array) (*core.Array, error) {
|
||
j, err := core.FromComplexes([]complex128{1, 0, 0, 1}, 2, 2)
|
||
return j, err
|
||
}
|
||
if _, _, err := LevenbergMarquardt(residual, p0, LMOptions{Jacobian: complexJacobian}); err == nil {
|
||
t.Fatal("expected a complex Jacobian to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("Jacobian: error = %v, want a complex-input refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestFindRootSystemComplexResidualRefused pins the dtype guard on the
|
||
// root-system residual, which used to dereference the nil int payload
|
||
// of a complex array.
|
||
func TestFindRootSystemComplexResidualRefused(t *testing.T) {
|
||
f := func(*core.Array) (*core.Array, error) {
|
||
r, err := core.FromComplexes([]complex128{1, 2}, 2)
|
||
return r, err
|
||
}
|
||
if _, _, err := FindRootSystem(f, mustFloats(t, []float64{1, 1}, 2), RootSystemOptions{}); err == nil {
|
||
t.Fatal("expected a complex residual to be refused")
|
||
} else if !strings.Contains(err.Error(), "complex") {
|
||
t.Fatalf("error = %v, want a complex-input refusal", err)
|
||
}
|
||
}
|
||
|
||
// TestLevenbergMarquardtChi2MatchesReturnedPoint pins the reported fit
|
||
// quality: the relative-improvement break published the new, lower χ²
|
||
// while the parameters were still the old ones, so the answer looked
|
||
// 99.9 % better than the point that came back.
|
||
func TestLevenbergMarquardtChi2MatchesReturnedPoint(t *testing.T) {
|
||
xs := []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9}
|
||
ys := make([]float64, len(xs))
|
||
for i, x := range xs {
|
||
ys[i] = 3*math.Exp(-0.5*x) + 0.5 + 0.001*math.Sin(7*x)
|
||
}
|
||
residual := func(p *core.Array) (*core.Array, error) {
|
||
a, b, c := p.FloatAt(0), p.FloatAt(1), p.FloatAt(2)
|
||
out := core.New(core.Float, len(xs))
|
||
for i := range xs {
|
||
out.RawFloats()[i] = ys[i] - (a*math.Exp(-b*xs[i]) + c)
|
||
}
|
||
return out, nil
|
||
}
|
||
point, chi2, err := LevenbergMarquardt(residual, mustFloats(t, []float64{2, 0.3, 0.1}, 3),
|
||
LMOptions{Tolerance: 0.05})
|
||
if err != nil {
|
||
t.Fatalf("LevenbergMarquardt: %v", err)
|
||
}
|
||
r, err := residual(point)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
actual := 0.0
|
||
for i := range r.Len() {
|
||
actual += r.FloatAt(i) * r.FloatAt(i)
|
||
}
|
||
if math.Abs(chi2-actual) > 1e-9*math.Max(1, actual) {
|
||
t.Fatalf("reported χ² = %.14g, χ² at the returned point = %.14g", chi2, actual)
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedExactFixtures checks the augmented Lagrangian
|
||
// against hand-derived optima: min (x−1)² + (y−2)² subject to x+y = 5
|
||
// projects to (2, 3) with value 2; min (x−10)² + (y−10)² subject to
|
||
// x+y = 1 and x ≤ −4 bottoms out at (−4, 5) with value 221; and the
|
||
// degenerate row min x+3 subject to x = 1 reaches 4.
|
||
func TestMinimiseConstrainedExactFixtures(t *testing.T) {
|
||
bowlPeak, err := core.FromFloats([]float64{1, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
boxCut, err := core.FromFloats([]float64{1, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
cases := []struct {
|
||
name string
|
||
f func(*core.Array) (float64, error)
|
||
A *core.Array
|
||
lower []float64
|
||
upper []float64
|
||
boxUpper []float64
|
||
start []float64
|
||
wantX []float64
|
||
wantValue float64
|
||
}{
|
||
{
|
||
name: "equality cuts the unconstrained minimum",
|
||
f: func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-1, p.FloatAt(1)-2
|
||
return dx*dx + dy*dy, nil
|
||
},
|
||
A: bowlPeak,
|
||
lower: []float64{5}, upper: []float64{5},
|
||
start: []float64{0, 0}, wantX: []float64{2, 3}, wantValue: 2,
|
||
},
|
||
{
|
||
name: "row cut by a box wall",
|
||
f: func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-10, p.FloatAt(1)-10
|
||
return dx*dx + dy*dy, nil
|
||
},
|
||
A: boxCut,
|
||
lower: []float64{1}, upper: []float64{1},
|
||
boxUpper: []float64{-4, math.Inf(1)},
|
||
start: []float64{0, 0}, wantX: []float64{-4, 5}, wantValue: 221,
|
||
},
|
||
}
|
||
for _, c := range cases {
|
||
point, value, err := MinimiseConstrained(c.f, nil, mustFloats(t, c.start, len(c.start)),
|
||
LinearConstraints{A: c.A, Lower: c.lower, Upper: c.upper},
|
||
LBFGSOptions{Tolerance: 1e-10, Upper: c.boxUpper})
|
||
if err != nil {
|
||
t.Errorf("%s: %v", c.name, err)
|
||
continue
|
||
}
|
||
for i, want := range c.wantX {
|
||
if math.Abs(point.FloatAt(i)-want) > 1e-4 {
|
||
t.Errorf("%s: x[%d] = %.12g, want %.12g", c.name, i, point.FloatAt(i), want)
|
||
}
|
||
}
|
||
if math.Abs(value-c.wantValue) > 1e-4*math.Max(1, c.wantValue) {
|
||
t.Errorf("%s: value = %.12g, want %.12g", c.name, value, c.wantValue)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedBadlyScaledRow pins the badly scaled row: min
|
||
// x² + y² subject to a·x + y = 1 has the closed form x = a/(a²+1),
|
||
// y = 1/(a²+1) with value 1/(a²+1). The row keeps the caller's scale,
|
||
// so the penalty's gradient at the start is mu·a and its curvature
|
||
// mu·a², which put the step that reduces the augmented Lagrangian
|
||
// below the old 20-halving line search for every a ≥ 1000: the inner
|
||
// solve took the silent stall exit and the outer loop reported "40
|
||
// rounds left the worst row violation at 1" without moving. The rows
|
||
// the long line search can reach must now reach the closed form; rows
|
||
// beyond its reach must be refused with the stall diagnostic, never
|
||
// returned as a converged answer.
|
||
func TestMinimiseConstrainedBadlyScaledRow(t *testing.T) {
|
||
f := func(p *core.Array) (float64, error) {
|
||
return p.FloatAt(0)*p.FloatAt(0) + p.FloatAt(1)*p.FloatAt(1), nil
|
||
}
|
||
rowResidual := func(a float64, point *core.Array) float64 {
|
||
return math.Abs(a*point.FloatAt(0) + point.FloatAt(1) - 1)
|
||
}
|
||
for _, a := range []float64{1e3, 1e4, 1e6, 1e8} {
|
||
A, err := core.FromFloats([]float64{a, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
point, value, err := MinimiseConstrained(f, nil, mustFloats(t, []float64{0, 0}, 2),
|
||
LinearConstraints{A: A, Lower: []float64{1}, Upper: []float64{1}}, LBFGSOptions{})
|
||
if err != nil {
|
||
t.Errorf("a=%g: %v", a, err)
|
||
continue
|
||
}
|
||
wantX, wantY, wantValue := a/(a*a+1), 1/(a*a+1), 1/(a*a+1)
|
||
if got := point.FloatAt(0); math.Abs(got-wantX) > 1e-6*wantX {
|
||
t.Errorf("a=%g: x = %.12g, want %.12g", a, got, wantX)
|
||
}
|
||
if got := point.FloatAt(1); math.Abs(got-wantY) > 1e-6*wantY {
|
||
t.Errorf("a=%g: y = %.12g, want %.12g", a, got, wantY)
|
||
}
|
||
if math.Abs(value-wantValue) > 1e-6*wantValue {
|
||
t.Errorf("a=%g: value = %.12g, want %.12g", a, value, wantValue)
|
||
}
|
||
if res := rowResidual(a, point); res > 1e-6 {
|
||
t.Errorf("a=%g: row residual = %g, want ≤ 1e-6", a, res)
|
||
}
|
||
}
|
||
// Past the line search's reach the answer is an error, not the start
|
||
// point dressed as convergence: a returned point must satisfy the
|
||
// row it claims to solve.
|
||
for _, a := range []float64{1e9, 1e12} {
|
||
A, err := core.FromFloats([]float64{a, 1}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
point, _, err := MinimiseConstrained(f, nil, mustFloats(t, []float64{0, 0}, 2),
|
||
LinearConstraints{A: A, Lower: []float64{1}, Upper: []float64{1}}, LBFGSOptions{})
|
||
if err != nil {
|
||
if !strings.Contains(err.Error(), "line search") {
|
||
t.Errorf("a=%g: error = %v, want the inner line search's stall diagnostic", a, err)
|
||
}
|
||
continue
|
||
}
|
||
if res := rowResidual(a, point); res > 1e-6 {
|
||
t.Errorf("a=%g: returned (%v) as converged with a row residual of %g, want the stall reported",
|
||
a, floatsOf(point), res)
|
||
}
|
||
}
|
||
}
|
||
|
||
// TestMinimiseConstrainedToleranceIsNotTheFeasibilityThreshold pins the
|
||
// decoupling: opts.Tolerance is the inner solver's projected-gradient
|
||
// tolerance and used to double as the absolute row-feasibility
|
||
// threshold (feasibleAt = max(Tolerance, 1e-10)), so a looser inner
|
||
// solve bought a looser row. For min (x−3)² + y² subject to 1 ≤ x ≤ 2
|
||
// (answer (2, 0), value 1) the old coupling returned (2.0033, 0) at
|
||
// Tolerance = 1e-2, a row violation of 3.3e-3, and a tolerance below
|
||
// 1e-10 hard-failed the solve with "40 rounds left the worst row
|
||
// violation at 1.3e-9". The row is now judged against the fixed 1e-10
|
||
// at every inner tolerance.
|
||
func TestMinimiseConstrainedToleranceIsNotTheFeasibilityThreshold(t *testing.T) {
|
||
f := func(p *core.Array) (float64, error) {
|
||
dx := p.FloatAt(0) - 3
|
||
return dx*dx + p.FloatAt(1)*p.FloatAt(1), nil
|
||
}
|
||
A, err := core.FromFloats([]float64{1, 0}, 1, 2)
|
||
if err != nil {
|
||
t.Fatal(err)
|
||
}
|
||
for _, tol := range []float64{0, 1e-8, 1e-10, 1e-12, 1e-4, 1e-2} {
|
||
point, value, err := MinimiseConstrained(f, nil, mustFloats(t, []float64{0, 0}, 2),
|
||
LinearConstraints{A: A, Lower: []float64{1}, Upper: []float64{2}},
|
||
LBFGSOptions{Tolerance: tol})
|
||
if err != nil {
|
||
t.Errorf("Tolerance=%g: %v", tol, err)
|
||
continue
|
||
}
|
||
if math.Abs(point.FloatAt(0)-2) > 1e-6 || math.Abs(point.FloatAt(1)) > 1e-6 {
|
||
t.Errorf("Tolerance=%g: point = (%v), want (2, 0)", tol, floatsOf(point))
|
||
}
|
||
if violation := math.Abs(point.FloatAt(0) - 2); violation > 1e-8 {
|
||
t.Errorf("Tolerance=%g: row violation = %g, want ≤ 1e-8", tol, violation)
|
||
}
|
||
if math.Abs(value-1) > 1e-6 {
|
||
t.Errorf("Tolerance=%g: value = %.12g, want 1", tol, value)
|
||
}
|
||
}
|
||
}
|
||
|
||
// bowlAt returns the separable bowl centred at (cx, cy), the fixture
|
||
// the constrained tests share.
|
||
func bowlAt(cx, cy float64) func(*core.Array) (float64, error) {
|
||
return func(p *core.Array) (float64, error) {
|
||
dx, dy := p.FloatAt(0)-cx, p.FloatAt(1)-cy
|
||
return dx*dx + dy*dy, nil
|
||
}
|
||
}
|
||
|
||
// floatsOf copies an array's values out for readable failure messages.
|
||
func floatsOf(a *core.Array) []float64 {
|
||
out := make([]float64, a.Len())
|
||
for i := range out {
|
||
out[i] = a.FloatAt(i)
|
||
}
|
||
return out
|
||
}
|