Finite-budget Levenberg-Marquardt with a trust-ratio safeguard

This algorithm performs one finite-budget LM trial update. Inputs are $x\in U$, a $C^1$ residual function $r:U\to {\mathbb{R}}^m$ on open $U$ with $m,n\ge 1$ and $J=Dr$, the exact values $r(x)$ and $J(x)$, a fixed nonsingular $D\in {\mathbb{R}}^{n\times n}$, an arithmetic mode (exact or computed), a floor ${\lambda }_{\min }>0$, current damping $\lambda \ge {\lambda }_{\min }$, factors $\gamma >1$, $0<{\eta }_1<{\eta }_2<1$, a trust radius $\mathrm{\Delta }>0$, and an attempt cap $N\ge 1$. In computed mode the input values are finite stored numbers treated as the exact data for this update; only arithmetic performed below is rounded. A computed run may include a certificate that the exact scalar product $J(x{)}^{T}r(x)$ is zero. Put $\varphi (x)=\mathrm{\parallel }r(x){\mathrm{\parallel }}_2^2\mathrm{/}2$ and $m_0(p)=\mathrm{\parallel }r(x)+J(x)p{\mathrm{\parallel }}_2^2\mathrm{/}2$. In exact mode the step is the unique minimizer of $\frac{1}{2}\mathrm{\parallel }r(x)+J(x)p{\mathrm{\parallel }}_2^2+\frac{\lambda }{2}\mathrm{\parallel }Dp{\mathrm{\parallel }}_2^2,$ which is equivalent to solving $(J(x{)}^{T}J(x)+\lambda D^{T}D)p=-J(x{)}^{T}r(x)\mathrm{.}$ In computed mode the returned vector is only a candidate for this same subproblem; no exact minimizer claim is made.

In computed mode use Operationwise real and complex floating point arithmetic with rounded square roots and fused multiply-add. Evaluate every dot product by initializing $s=0$ and, for indices in increasing order, computing $u=\mathrm{fl}(a_i{b}_i)$ followed by $s=\mathrm{fl}(s+u)$. Matrix-vector products and squared norms use this schedule row by row and coordinate by coordinate. Each other scalar arithmetic operation uses $\mathrm{fl}$, with copies and sign changes exact; the first failed operation returns the status assigned to its stage. Function evaluation returns either a finite stored residual vector or failure.

For attempts $j=1,\dots ,N$, do the following with the current damping. After a damping change, construct a new augmented matrix and obtain a report for that actual matrix; do not reuse an earlier damping value's factorization or report.

1. In exact arithmetic, if $J(x{)}^{T}r(x)=0$, return stationary before calling the linear solver. In a computed run, evaluate $J(x{)}^{T}r(x)$ by the stated dot schedule and return stationary before the solve only with the supplied certificate that this exact gradient is zero; if evaluating the product fails, return gradient-evaluation-failure; if the computed product is a rounded zero without that certificate, return gradient-unverified and do not form a ratio. Otherwise perform the augmented construction and linear solve afresh inside this attempt at its current damping: Form $A_{\lambda }=[J(x);\sqrt{\lambda }D]$ and $b_{\lambda }=[r(x);0]$, then call Gauss-Newton linearized least-squares step with QR/SVD rank safeguards in QR mode on $(A_{\lambda },b_{\lambda })$. Exact mode uses the exact positive square root and requires a supplied exact reduced factorization of $A_{\lambda }$ for each call; absence of that factorization returns linear-solve-failure. This matrix has full column rank because $D$ is nonsingular and $\lambda >0$. In computed mode, copy the stored $J(x)$ and $r(x)$ into the top blocks, compute $t=\mathrm{fl}(\sqrt{\lambda })$, then fill the bottom block entries $\mathrm{fl}(t{D}_{ij})$ in increasing row and column order and fill the appended residual entries with exact zeros. Any first arithmetic or range failure in this construction returns linear-solve-failure with its stage and indices. The call receives a computed QR report for this actual stored $A_{\lambda }$. A computed-candidate is usable only when that report includes the low-level procedure's independently verified full-rank certificate $s_{\lambda }>0$ for the actual stored matrix, namely $\mathrm{\parallel }{A}_{\lambda }v{\mathrm{\parallel }}_2\ge s_{\lambda }\mathrm{\parallel }v{\mathrm{\parallel }}_2$ for every $v\in {\mathbb{R}}^n$. Pass this bound as the computed report's certificate field; it is not inferred from a nonzero computed diagonal or from positive definiteness of the ideal exact system alone. A rank-unverified report and every other status besides exact-step in exact mode or computed-candidate in computed mode are unusable and return linear-solve-failure with the report. A computed-candidate remains a candidate and carries no convergence guarantee. If the solver status is not exact-step in exact mode or computed-candidate in computed mode, return linear-solve-failure with its report; in particular, rank-unverified is a failure, not a usable trial.
2. If exact mode returns $p=0$ after the nonzero-gradient check, return linear-solve-failure with the inconsistent exact solve report. If computed mode returns $p=0$, return zero-step-unverified and do not compute a ratio. For every remaining candidate $p$ (the zero candidates have already returned a status), compute $\mathrm{pred}=m_0(0)-m_0(p)$. If this computation fails, return model-evaluation-failure; if $\mathrm{pred}\le 0$, return nonpositive-predicted-decrease. In either case do not compute a ratio.
3. If computing $x+p$ or $Dp$ fails in the selected arithmetic, return trial-evaluation-failure. If $x+p\mathrm{\notin }U$ or $\mathrm{\parallel }Dp{\mathrm{\parallel }}_2>\mathrm{\Delta }$, reject the trial and compute the damping update ${\lambda }^{+}=\gamma \lambda$ in the selected arithmetic. If that update fails, is nonfinite, or is outside the representable range, return damping-range-failure with the current $x,\lambda$ and the trial report. Otherwise set $\lambda \leftarrow {\lambda }^{+}$ and continue.
4. Otherwise evaluate $r(x+p)$ and compute $\mathrm{ared}=\varphi (x)-\varphi (x+p)$ and $\rho =\mathrm{ared}\mathrm{/}\mathrm{pred}$. A failed trial evaluation or nonfinite computed result returns trial-evaluation-failure. If $\rho <{\eta }_1$, reject and compute ${\lambda }^{+}=\gamma \lambda$ in the selected arithmetic. If that damping update fails, is nonfinite, or is outside the representable range, return damping-range-failure with the current $x,\lambda$ and the trial report. Otherwise set $\lambda \leftarrow {\lambda }^{+}$. If this was attempt $N$, return rejection-budget-exhausted; otherwise continue.
5. If $\rho \ge {\eta }_1$, first determine the next damping. If also $\rho \ge {\eta }_2$ and $\mathrm{\parallel }Dp{\mathrm{\parallel }}_2=\mathrm{\Delta }$, compute $q=\lambda \mathrm{/}\gamma$ in exact mode and $q=\mathrm{fl}(\lambda \mathrm{/}\gamma )$ in computed mode. If this division fails, has a nonfinite or nonpositive result, or is outside the representable range, return damping-range-failure with the current $x,\lambda$ and the trial report, without accepting that trial. Otherwise set ${\lambda }_{\mathrm{n}\mathrm{e}\mathrm{x}\mathrm{t}}=\max \{{\lambda }_{\min },q\}$. In every other accepted case set ${\lambda }_{\mathrm{n}\mathrm{e}\mathrm{x}\mathrm{t}}=\lambda$. Only after the next damping is available, accept and return $x+p$, $p$, $\rho$, and ${\lambda }_{\mathrm{n}\mathrm{e}\mathrm{x}\mathrm{t}}$.

If all $N$ attempts are rejected, return rejection-budget-exhausted. The equality in the boundary test is exact for exact arithmetic and is the equality test on stored values in a computed run. No claim is made that every finite budget succeeds.