`vp.optimize.least_squares`

Minimize a sum of squares of nonlinear residuals by MINPACK's Levenberg-Marquardt, with the residuals in Viv.

Part of the distribution vipyr; it runs from a checkout of the repository. Planned on PyPI as pip install vipyr.

vp.optimize.least_squares is scipy.optimize.least_squares(method="lm") with the solver and the residuals compiled together: lib/minpack's port of MINPACK's lmder calls your function directly, with no Python in the loop, and takes the steps the Fortran takes. The defaults are SciPy's (ftol = xtol = gtol = 1e-8, x_scale = 1, max_nfev = 100 n), so a call written for SciPy gives the same iterates.

Give the residuals as

a def of your own file, "model.viv:residuals", or k.residuals after k = vp.load("model.viv");
a Python function in py2viv's numeric subset (tools/py2viv/SUBSET.md), translated to Viv once and cached.

The function is called as fun(x, *args) and returns the m residuals as an array. Each of args is a number or a 1-D array (data, measurements); they are passed to the compiled function on every call, so one compilation serves every data set of the same form. The first call for a function writes a small driver .viv and compiles it (a few seconds); later calls, in this process or the next, load it from vipyr's cache.

The Jacobian is exact by default: the driver differentiates the residuals by reverse-mode grad, one row per residual, so no step size has to be chosen and the iteration does not depend on difference errors. jac="fd" takes MINPACK's forward differences (fdjac2), which reproduces what SciPy does when no Jacobian is given, evaluation for evaluation.

Signature

vp.optimize.least_squares(fun, x0, *, jac=None, args=(), ftol=1e-08, xtol=1e-08, gtol=1e-08, x_scale=1.0, max_nfev=None, diff_step=None, fast=False)

Parameters

fun : str | function, required: The residuals r(x, *args), a Viv def returning an array of m: 'file.viv:def', a def of a module from vp.load, or a Python function in py2viv's subset.
x0 : ndarray[float64], required: The starting point, length n <= m.
jac : str | function, default None: None: exact, by reverse-mode grad of the residuals; 'fd' (or '2-point'): MINPACK's forward differences; or a def (a name of fun's file, or a function) jac(x, *args) returning the m by n Jacobian.
args : tuple, default (): Extra arguments of fun and jac, each a number or a 1-D array.
ftol : float, default 1e-08: Stop when the relative reduction of the sum of squares is at most ftol.
xtol : float, default 1e-08: Stop when the relative change of x is at most xtol.
gtol : float, default 1e-08: Stop when the largest cosine between the residuals and a column of the Jacobian is at most gtol.
x_scale : float | ndarray[float64] | str, default 1.0: The scale of each variable (MINPACK's diag is 1 / x_scale); 'jac' scales by the Jacobian's column norms.
max_nfev : int, default None: The residual evaluations allowed; None: 100 n, or 100 n (n + 1) with jac='fd', as SciPy.
diff_step : float, default None: With jac='fd', the relative step; None: sqrt(eps).
fast : bool, default False: Compile with the fast flags (fused multiply-add): faster, last bits may differ.

Returns

A frozen LeastSquaresResult. Iterating it gives the fields below in order; it also has params, tests, log and check().

x : ndarray[float64]: The solution.
cost : float: Half the sum of squares of the residuals at x.
fun : ndarray[float64]: The residuals at x.
jac : ndarray[float64]: The m by n Jacobian at x.
grad : ndarray[float64]: The gradient of cost at x, jac^T fun.
optimality : float: The largest absolute entry of grad.
active_mask : ndarray[int64]: Zeros: method 'lm' has no bounds.
nfev : int: Residual evaluations (with jac='fd', including the n of each difference Jacobian, as SciPy counts them).
njev : int: Jacobian evaluations.
status : int: SciPy's: -1 improper input, 0 max_nfev reached, 1 gtol, 2 ftol, 3 xtol, 4 ftol and xtol.
message : str: The status in words.
success : bool: Status > 0.
info : int: MINPACK's info, 0 to 8, which status is made from.

Raises

ValueError: A tolerance below machine epsilon, fewer residuals than variables, or an array of the wrong shape; vipyr._optimize.OptimizeProblemError (a ValueError) when the function or jac cannot be made into a Viv problem, or py2viv refuses a Python function.
vipyr.VivCompileError: A SyntaxError: the function's .viv file (or py2viv's translation) does not compile.
TypeError: An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError: An ImportError: no C compiler, or the program's C does not build.

Notes

The method. Each iteration solves the trust-region problem

p min \frac{1}{2} ∥ J p + r ∥^{2} subject to ∥ D p ∥ \leq Δ

for the step $p$ , with $J$ the Jacobian of the residuals $r$ at the current $x$ and $D$ the diagonal scaling (x_scale, or the column norms of $J$ with x_scale="jac"). lmpar finds the Levenberg-Marquardt parameter $λ$ for which $p = - (J^{T} J + λ D^{2})^{- 1} J^{T} r$ has $∥ D p ∥ \approx Δ$ , from a QR factorization of $J$ with column pivoting (qrfac) and a Givens update for each $λ$ (qrsolv). The step is accepted when the actual reduction of $∥ r ∥$ is at least $1 0^{- 4}$ of the reduction the linear model predicts; the radius shrinks when the ratio is below 1/4 and doubles when it is above 3/4. This is Moré's 1978 implementation, unchanged.

Stopping. ftol: the actual and predicted relative reductions of the sum of squares are both at most ftol. xtol: the relative change of $D x$ is at most xtol. gtol: the largest cosine between $r$ and a column of $J$ is at most gtol. status says which (SciPy's numbering) and info gives MINPACK's own code.

Exactness. lib/minpack/lmder.viv keeps the Fortran's operations in the Fortran's order; on the 24 Moré-Garbow-Hillstrom cases of MINPACK's test driver it reproduces a gfortran build of netlib's lmder iterate by iterate (lib/minpack/mgh_check.viv). Through vp.optimize the residuals are Viv (or py2viv's translation of your Python), so the counts agree with SciPy's whenever the residuals round the same way; lib/minpack/bench_optimize.json lists every case.

Limits. No bounds: MINPACK's lmder has none, and neither has SciPy's method="lm" (SciPy's bounded methods, 'trf' and 'dogbox', are not ported). Dense Jacobians only. With the exact Jacobian, grad runs once per residual, so a Jacobian costs about m evaluations of the whole residual vector; for thousands of residuals jac="fd" (n evaluations) is cheaper.

References

J. J. Moré, The Levenberg-Marquardt algorithm: implementation and theory, in Numerical Analysis (Dundee 1977), Lecture Notes in Mathematics 630, 105-116, Springer (1978). doi:10.1007/BFb0067700
J. J. Moré, B. S. Garbow and K. E. Hillstrom, User Guide for MINPACK-1, Argonne National Laboratory report ANL-80-74 (1980); the Fortran is netlib's, in lib/minpack/ref. https://www.netlib.org/minpack/
J. J. Moré, B. S. Garbow and K. E. Hillstrom, Testing unconstrained optimization software, ACM Trans. Math. Softw. 7, 17-41 (1981). doi:10.1145/355934.355936
P. Virtanen et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261-272 (2020): scipy.optimize.least_squares, whose interface this follows. doi:10.1038/s41592-019-0686-2

Examples

import numpy as np
import vipyr as vp

def rosenbrock(x):
    return np.array([10.0 * (x[1] - x[0] ** 2), 1.0 - x[0]])

# translated to Viv by py2viv once, compiled with lmder, Jacobian by grad
r = vp.optimize.least_squares(rosenbrock, [-1.2, 1.0])
print(r.x, r.nfev, r.njev, r.status)
assert r.check()

[1. 1.] 21 16 1

Residuals with data passed through args: one compiled function serves any data set of the same form. Forward differences give SciPy's evaluation count.

import numpy as np
import vipyr as vp

def decay(b, t, y):
    return b[0] * np.exp(-b[1] * t) + b[2] - y

t = np.linspace(0.0, 4.0, 40)
y = 3.0 * np.exp(-1.3 * t) + 0.5 + 0.01 * np.sin(7.0 * t)
for jac in (None, "fd"):
    r = vp.optimize.least_squares(decay, [1.0, 1.0, 0.0], args=(t, y), jac=jac)
    print(jac, np.round(r.x, 6), r.nfev, r.njev, r.status)

None [3.004324 1.304273 0.501609] 7 6 4
fd [3.004324 1.304273 0.501609] 25 6 4

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference: lib/minpack/lmder.viv matches a gfortran build of netlib's lmder on 24 Moré-Garbow-Hillstrom cases (lib/minpack/mgh_check.viv); through vp.optimize, NIST's 27 StRD nonlinear regression problems from both certified starts and 18 MGH cases against SciPy 1.12's least_squares(method='lm') (lib/minpack/bench_optimize.py, lib/minpack/bench_optimize.json; python/tests/test_optimize.py).

Credits

A port of MINPACK (Argonne National Laboratory), under MINPACK's BSD-style copyright notice (1999, University of Chicago), which the port and the Fortran in lib/minpack/ref both carry. This product includes software developed by the University of Chicago, as Operator of Argonne National Laboratory.