Notebooks · lib/minpack/lmder

MINPACK's lmder, Levenberg-Marquardt least squares with a trust region

A port of lmder.f and lmder1.f from netlib MINPACK (lib/minpack/ref), with the Fortran's arithmetic in the Fortran's order, so that the iterates match a gfortran build of the original to the last bit where the problem's own functions do. The explanation is in lib/minpack/lmder.intent; the checks are in lib/minpack/mgh_check.viv and lib/minpack/check.sh.

The function. fcn(data, x, fvec, fjac, iflag) is a def: for iflag 1 it writes the m residuals at x into fvec, for 2 the m x n Jacobian into fjac, for 0 it may print; it returns iflag, or a negative number to stop lmder. data is any array the function needs (a problem number, measurements).

Differences from the Fortran, none of which changes an operation:

indices count from 0, and ipvt holds 0-based column numbers;
fjac is an m x n array instead of a leading-dimension array;
the work arrays wa1 ... wa4 are allocated here instead of passed in;
info is returned, and nfev and njev are written into counts[0] and counts[1];
every trial step writes one row into hist (the Fortran has no such record; it is what lets a test compare iterates, and a plot draw them). The columns are the H_ constants below. A hist with fewer rows than trial steps keeps the first rows; pass zeros(1, H_X + n) to keep none worth reading.

Used with import minpack.lmder.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 32–534

from minpack.lmpar import lmpar
from minpack.qrfac import qrfac
from minpack.enorm import enorm, EPSMCH

# The columns of one row of hist: one row per trial step.
H_ITER = 0
H_NFEV = 1
H_DELTA_IN = 2
H_PAR_LM = 3
H_PNORM = 4
H_FNORM1 = 5
H_ACTRED = 6
H_PRERED = 7
H_RATIO = 8
H_DELTA = 9
H_PAR = 10
H_ACCEPT = 11
H_FNORM = 12
H_GNORM = 13
H_X = 14

# The thresholds of the trust region, from the data statement of lmder.f
# (the Fortran's names in parentheses). The halves (p5 there: `temp = 0.5`,
# `delta = pnorm / 0.5`, `par = 0.5 * par`, `/ 0.5` in prered and
# `0.5 * ratio`) stay written as 0.5, as the formulas they are part of.
#
# A trial step whose ratio of actual to predicted reduction is at most this
# shrinks the radius (p25).
LM_RATIO_SHRINK = 0.25
# A ratio at least this expands it to twice the step's length (p75).
LM_RATIO_EXPAND = 0.75
# A ratio at least this accepts the step: x moves to x + p (p0001).
LM_RATIO_ACCEPT = 0.0001
# The smallest factor t a shrink multiplies the radius by; the shrunk
# radius is t * min(delta, pnorm / LM_SHRINK_MIN), so it is also at most 10 t
# times the step's length (p1).
LM_SHRINK_MIN = 0.1
# A trial norm fnorm1 at least 1 / LM_RISE_LIMIT = 10 times fnorm counts as
# actred = -1 without squaring the quotient (which could overflow), and
# shrinks the radius by the full LM_SHRINK_MIN (p1).
LM_RISE_LIMIT = 0.1
# lmder1's settings (lmder1.f): the evaluations of f allowed per variable,
# plus one, and the initial radius as a multiple of |D x| (factor).
LMDER1_MAXFEV_PER = 100.0
LMDER1_FACTOR = 100.0

def lmder(fcn, data: Arr, x: Arr, fvec: Arr, fjac: Arr, ftol, xtol, gtol, maxfev, diag: Arr, mode, factor, nprint, ipvt: Arr, qtf: Arr, counts: Arr, hist: Arr) {
    m = len(fvec)
    n = len(x)
    step setup {
        wa1 = zeros(n)
        wa2 = zeros(n)
        wa3 = zeros(n)
        wa4 = zeros(m)
        epsmch = EPSMCH
        info = 0
        iflag = 0
        nfev = 0
        njev = 0
        nrow = 0

from minpack.lmpar import lmpar
from minpack.qrfac import qrfac
from minpack.enorm import enorm, EPSMCH

# The columns of one row of hist: one row per trial step.
H_ITER = 0
H_NFEV = 1
H_DELTA_IN = 2
H_PAR_LM = 3
H_PNORM = 4
H_FNORM1 = 5
H_ACTRED = 6
H_PRERED = 7
H_RATIO = 8
H_DELTA = 9
H_PAR = 10
H_ACCEPT = 11
H_FNORM = 12
H_GNORM = 13
H_X = 14

# The thresholds of the trust region, from the data statement of lmder.f
# (the Fortran's names in parentheses). The halves (p5 there: `temp = 0.5`,
# `delta = pnorm / 0.5`, `par = 0.5 * par`, `/ 0.5` in prered and
# `0.5 * ratio`) stay written as 0.5, as the formulas they are part of.
#
# A trial step whose ratio of actual to predicted reduction is at most this
# shrinks the radius (p25).
LM_RATIO_SHRINK = 0.25
# A ratio at least this expands it to twice the step's length (p75).
LM_RATIO_EXPAND = 0.75
# A ratio at least this accepts the step: x moves to x + p (p0001).
LM_RATIO_ACCEPT = 0.0001
# The smallest factor t a shrink multiplies the radius by; the shrunk
# radius is t * min(delta, pnorm / LM_SHRINK_MIN), so it is also at most 10 t
# times the step's length (p1).
LM_SHRINK_MIN = 0.1
# A trial norm fnorm1 at least 1 / LM_RISE_LIMIT = 10 times fnorm counts as
# actred = -1 without squaring the quotient (which could overflow), and
# shrinks the radius by the full LM_SHRINK_MIN (p1).
LM_RISE_LIMIT = 0.1
# lmder1's settings (lmder1.f): the evaluations of f allowed per variable,
# plus one, and the initial radius as a multiple of |D x| (factor).
LMDER1_MAXFEV_PER = 100.0
LMDER1_FACTOR = 100.0

def lmder(fcn, data: Arr, x: Arr, fvec: Arr, fjac: Arr, ftol, xtol, gtol, maxfev, diag: Arr, mode, factor, nprint, ipvt: Arr, qtf: Arr, counts: Arr, hist: Arr) {
    m = len(fvec)
    n = len(x)
    step setup {
        wa1 = zeros(n)
        wa2 = zeros(n)
        wa3 = zeros(n)
        wa4 = zeros(m)
        epsmch = EPSMCH
        info = 0
        iflag = 0
        nfev = 0
        njev = 0
        nrow = 0