127 lines
4.3 KiB
Go
127 lines
4.3 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
||
// SPDX-License-Identifier: MIT
|
||
|
||
package optim
|
||
|
||
import (
|
||
"math"
|
||
|
||
"sourcedock.dev/petrbalvin/tensor/internal/base"
|
||
)
|
||
|
||
// Broyden's quasi-Newton maintenance for FindRootSystem, enabled by
|
||
// RootSystemOptions.UseBroyden. The central-difference Jacobian is
|
||
// built once at the start and the iteration walks on its maintained
|
||
// inverse: after every accepted step the rank-one update corrects the
|
||
// inverse so that it satisfies the secant equation H·y = s for the
|
||
// step just taken.
|
||
//
|
||
// The update maintained here is the bad Broyden form in its inverse
|
||
// shape, H + (s − H·y)·yᵀ/(yᵀy) with H the maintained inverse: a
|
||
// rank-one correction along the residual-difference direction y that
|
||
// satisfies the secant equation exactly for the step just taken, and
|
||
// turns every later iteration into a single matrix-vector product,
|
||
// δ = −H·r. The good form, whose correction runs along sᵀH instead,
|
||
// satisfies the same equation with a different matrix; maintaining it
|
||
// needs the product sᵀH·y on top of the step, which spends back the
|
||
// per-iteration saving the option exists for.
|
||
|
||
// broydenInvert fills h with the inverse of the central-difference
|
||
// Jacobian jac by solving jac·h = I column-wise through the library's
|
||
// LU solver: one factorisation, n right-hand columns, the same cost
|
||
// class as the Newton path's single solve. jacWork is consumed in
|
||
// place, exactly as the Newton path consumes it, so jac stays
|
||
// pristine for the steepest-descent fallback. A singular Jacobian is
|
||
// returned as the solver's own error for the caller to fall back on.
|
||
func broydenInvert(jac, jacWork, h [][]float64) error {
|
||
n := len(jac)
|
||
rhs := make([][]float64, n)
|
||
for i := range n {
|
||
rhs[i] = make([]float64, n)
|
||
rhs[i][i] = 1
|
||
}
|
||
for i := range n {
|
||
copy(jacWork[i], jac[i])
|
||
}
|
||
sol, err := base.SolveSystem("FindRootSystem", jacWork, rhs)
|
||
if err != nil {
|
||
return err
|
||
}
|
||
// sol[k] is jac⁻¹·eₖ, the kth column of the inverse.
|
||
for i := range n {
|
||
for k := range n {
|
||
h[i][k] = sol[k][i]
|
||
}
|
||
}
|
||
return nil
|
||
}
|
||
|
||
// broydenMaintain applies the rank-one update
|
||
//
|
||
// H ← H + (s − H·y)·yᵀ/(yᵀy)
|
||
//
|
||
// to the maintained inverse, where s is the accepted step and y the
|
||
// residual change it produced, and reports whether the inverse is fit
|
||
// to carry forward, together with the running count of consecutive
|
||
// steps that failed to lower the residual infinity norm. A false
|
||
// report leaves h untouched and orders FindRootSystem to rebuild the
|
||
// Jacobian numerically before it steps again, the restart the option
|
||
// documents. Two observations order the restart, both a degradation of
|
||
// the rank-one model:
|
||
//
|
||
// - the update denominator yᵀy is zero, non-finite, or at rounding
|
||
// level against the residual's own scale (‖y‖∞ ≤ ε·max(1, ‖r‖∞)):
|
||
// the division would amplify cancellation noise into H, and a y
|
||
// that small carries no curvature information at all;
|
||
// - two consecutive accepted steps each failed to lower the
|
||
// residual infinity norm. The damping guarantees the residual sum
|
||
// of squares falls on every accepted step, so a flat infinity
|
||
// norm twice in a row means the maintained inverse has stopped
|
||
// predicting the landscape and a fresh Jacobian is cheaper than
|
||
// more crawling.
|
||
func broydenMaintain(h [][]float64, s, y, r []float64, res, resPrev float64, stalled int) (bool, int) {
|
||
den := 0.0
|
||
ynorm := 0.0
|
||
for i := range y {
|
||
den += y[i] * y[i]
|
||
if v := math.Abs(y[i]); v > ynorm {
|
||
ynorm = v
|
||
}
|
||
}
|
||
if den == 0 || math.IsNaN(den) || math.IsInf(den, 0) || ynorm <= base.EpsF*math.Max(1, normInfOfStep(r)) {
|
||
return false, 0
|
||
}
|
||
if res < resPrev {
|
||
stalled = 0
|
||
} else {
|
||
stalled++
|
||
if stalled >= 2 {
|
||
return false, 0
|
||
}
|
||
}
|
||
for i := range h {
|
||
hy := 0.0
|
||
for k := range y {
|
||
hy += h[i][k] * y[k]
|
||
}
|
||
w := (s[i] - hy) / den
|
||
for k := range y {
|
||
h[i][k] += w * y[k]
|
||
}
|
||
}
|
||
return true, stalled
|
||
}
|
||
|
||
// broydenStep writes the quasi-Newton step δ = −H·r into step: the
|
||
// whole per-iteration linear algebra the maintained inverse leaves,
|
||
// a single matrix-vector product.
|
||
func broydenStep(h [][]float64, r, step []float64) {
|
||
for i := range step {
|
||
s := 0.0
|
||
for k := range r {
|
||
s += h[i][k] * r[k]
|
||
}
|
||
step[i] = -s
|
||
}
|
||
}
|