Gauss-Newton and Levenberg-Marquardt methods

Let $n,m\ge 1$, let $U\subseteq {\mathbb{R}}^n$ be open, and let $r=(r_1,\dots ,r_m{)}^T:U\to {\mathbb{R}}^m$ be of class $C^2$ (each component has continuous partial derivatives through order two, as in $C^k$ maps and multi-index derivative notation in Euclidean space). Define $\varphi (x):=\frac{1}{2}{\sum }_{i=1}^m{r}_i(x{)}^2=\frac{1}{2}\mathrm{\parallel }r(x){\mathrm{\parallel }}_2^2,$ and let $J(x)\in {\mathbb{R}}^{m\times n}$ be the Jacobian with $i$th row $\mathrm{\nabla }{r}_i(x{)}^T$. Then $\varphi$ is of class $C^2$ and $\mathrm{\nabla }\varphi (x)=J(x{)}^{T}r(x),{\mathrm{\nabla }}^2\varphi (x)=J(x{)}^{T}J(x)+{\sum }_{i=1}^m{r}_i(x){\mathrm{\nabla }}^2{r}_i(x)\mathrm{.}$

The second term in the Hessian is the residual-curvature term. It vanishes when $r(x)=0$, but need not vanish at a nonzero-residual stationary point.

Let $U\subseteq {\mathbb{R}}^n$ be open, let $r:U\to {\mathbb{R}}^m$ be $C^1$ with $m,n\ge 1$, and let $J(x)=Dr(x)\in {\mathbb{R}}^{m\times n}$. Let $K\subseteq U$ be convex and let $x_{\ast },x\in K$. Suppose $L\ge 0$ and $\mathrm{\parallel }J(y)-J(z){\mathrm{\parallel }}_2\le L\mathrm{\parallel }y-z{\mathrm{\parallel }}_2$ for all $y,z\in K$, where the matrix norm is the Euclidean induced operator norm of The matrix norm induced by a published vector p-norm. Then $\mathrm{\parallel }r(x)-r(x_{\ast })-J(x)(x-x_{\ast }){\mathrm{\parallel }}_2\le \frac{L}{2}\mathrm{\parallel }x-x_{\ast }{\mathrm{\parallel }}_2^2\mathrm{.}$ The dimensions $m$ and $n$ need not be equal.

Let $J\in {\mathbb{R}}^{m\times n}$, $r\in {\mathbb{R}}^m$, $m,n\ge 1$, let $D\in {\mathbb{R}}^{n\times n}$ be nonsingular, let $\lambda >0$, and put $g:=J^{T}r$. Assume $g\mathrm{\ne }0$. There is a unique $p_{\lambda }$ satisfying $(J^{T}J+\lambda D^{T}D)p_{\lambda }=-g\mathrm{.}$ Set $z_{\lambda }:=D{p}_{\lambda }$, $\tilde{J}:=J{D}^{-1}$, and ${\mathrm{\Delta }}_{\lambda }:=\mathrm{\parallel }{z}_{\lambda }{\mathrm{\parallel }}_2>0$. Then $z_{\lambda }$ uniquely minimizes $q_{\lambda }(z):=\frac{1}{2}\mathrm{\parallel }r+\tilde{J}z{\mathrm{\parallel }}_2^2+\frac{\lambda }{2}\mathrm{\parallel }z{\mathrm{\parallel }}_2^2$ over ${\mathbb{R}}^n$, and globally minimizes the unregularized model $m(z):=\frac{1}{2}\mathrm{\parallel }r+\tilde{J}z{\mathrm{\parallel }}_2^2$ over the trust-region ball $\mathrm{\parallel }z{\mathrm{\parallel }}_2\le {\mathrm{\Delta }}_{\lambda }$ of Quadratic trust-region models and the actual-to-predicted reduction ratio. Its predicted reduction is positive and, in the original step coordinates, $m(0)-m(D{p}_{\lambda })=\frac{1}{2}\mathrm{\parallel }J{p}_{\lambda }{\mathrm{\parallel }}_2^2+\lambda \mathrm{\parallel }D{p}_{\lambda }{\mathrm{\parallel }}_2^2>0.$

Input a stored residual $r\in {\mathbb{R}}^m$, a stored Jacobian $J\in {\mathbb{R}}^{m\times n}$ with $m,n\ge 1$ (using the convention of The Jacobian matrix of partial derivatives and the gradient in the scalar-valued case), a factorization mode QR or SVD, an arithmetic mode (exact or computed), and a declared singular-value threshold $\tau \ge 0$. If $m<n$, immediately return unsupported-shape without requesting a factorization report; all remaining branches assume $m\ge n$. A computed QR report may also carry a verified full-rank certificate $s>0$ asserting $\mathrm{\parallel }Jv{\mathrm{\parallel }}_2\ge s\mathrm{\parallel }v{\mathrm{\parallel }}_2$ for every $v\in {\mathbb{R}}^n$. The target exact step, when it is returned, is the unique minimizer of $\frac{1}{2}\mathrm{\parallel }r+Jp{\mathrm{\parallel }}_2^2$ and hence solves $J^{T}Jp=-J^{T}r$; this target is asserted only on a full-column-rank branch. Do not form $J^{T}J$ in either solve path.

The threshold $\tau$ is used only by the SVD route. For QR mode the exact and computed routes are disjoint. Exact mode requires a supplied exact reduced factorization $J=QR$ with $Q^{T}Q=I_n$; only this exact factorization may reach the exact-step branch. Computed mode either consumes the report produced from this input pair $(J,-r)$ or forms it by setting $C=[J\mid -r]$ and, for $k=1,\dots ,n$, applying Scaled Householder reflector formation and application to rows $k,\dots ,m$ of column $k$ and all trailing columns of $C$, including its last right-hand-side column. Embed each returned active block in those rows, leaving earlier rows and columns unchanged, and propagate the first factorization failure. After a successful computed loop, the report consists of the computed upper-triangular block $\widehat{R}$, the first $n$ entries $\widehat{y}$ of the transformed last column, and the ordered stored reflector records (including identity markers). The exact orthogonal representatives of those records define the $Q$ used by the conditional backward-stability result; the rounded outputs themselves are not an exact factorization of $J$ and cannot be supplied to the exact branch. In exact mode, a zero diagonal in $R$ returns rank-deficient; otherwise solve $Rp=Q^T(-r)$ by backward substitution and return exact-step($p$). In computed mode, a zero diagonal returns rank-deficient-candidate. Otherwise initialize storage for ${\widehat{p}}_0,\dots ,{\widehat{p}}_{n-1}$ with all components unassigned, then for $i=n-1,n-2,\dots ,0$ compute $s_i$ by setting $s_i=0$ and, for $j=i+1,i+2,\dots ,n-1$ in increasing order, updating $s_i\leftarrow \mathrm{fl}(s_i+\mathrm{fl}({\widehat{R}}_{ij}{\widehat{p}}_j))$; then compute $t_i=\mathrm{fl}({\widehat{y}}_i-s_i)$ and ${\widehat{p}}_i=\mathrm{fl}(t_i\mathrm{/}{\widehat{R}}_{ii})$. Any first arithmetic or range failure returns linear-solve-failure with its stage and indices. A nonzero computed diagonal without an independently verified full-rank certificate $s>0$ for the input matrix $J$ returns rank-unverified with the computed candidate but no full-rank claim. With that certificate, return computed-candidate($\widehat{p}$). A nonzero computed diagonal alone is not a full-rank certificate. The computed-QR backward-stability conclusion of Computed reduced QR least squares is backward stable for full column rank is attached only when all of that theorem's successful-execution, arithmetic, and nonsingularity hypotheses are supplied; otherwise no backward-stability conclusion is returned.

For SVD mode, the input includes an exact factorization $J=U\mathrm{\Sigma }{V}^T$, or a computed report with finite stored factors ${\widehat{U}}_1\in {\mathbb{R}}^{m\times n}$, $\widehat{\mathrm{\Sigma }}=\mathrm{diag}({\widehat{\sigma }}_1,\dots ,{\widehat{\sigma }}_n)\in {\mathbb{R}}^{n\times n}$ whose entries are nonnegative, and $\widehat{V}\in {\mathbb{R}}^{n\times n}$, together with certified intervals $[{\mathrm{\ell }}_i,u_i]$ satisfying $0\le {\mathrm{\ell }}_i\le {\sigma }_i(J)\le u_i$ for each singular value $i=1,\dots ,n$. In the exact case, if any ${\sigma }_i\le \tau$, return rank-below-threshold and no Gauss-Newton step. If every ${\sigma }_i>\tau$, return $p=-V\mathrm{diag}({\sigma }_1^{-1},\dots ,{\sigma }_n^{-1})U_1^{T}r,$ where $U_1$ consists of the corresponding $n$ left singular vectors. For a computed report, if any $u_i\le \tau$, return rank-below-threshold; otherwise, if any interval satisfies ${\mathrm{\ell }}_i\le \tau <u_i$, return rank-unverified; if every lower endpoint exceeds $\tau$ but some stored ${\widehat{\sigma }}_i=0$, return factorization-failure. In the remaining branch, for $i=1,\dots ,n$ in increasing order initialize $z_i=0$ and, for $j=1,\dots ,m$ in increasing order, set $z_i\leftarrow \mathrm{fl}(z_i+\mathrm{fl}({\widehat{U}}_{ji}{r}_j))$. After these dots, for $i=1,\dots ,n$ set $w_i=\mathrm{fl}(z_i\mathrm{/}{\widehat{\sigma }}_i)$. For each $j=1,\dots ,n$ initialize $t_j=0$ and, for $i=1,\dots ,n$ in increasing order, set $t_j\leftarrow \mathrm{fl}(t_j+\mathrm{fl}({\widehat{V}}_{ji}{w}_i))$; set ${\widehat{p}}_j=-t_j$ by exact sign change. Any first arithmetic or range failure returns linear-solve-failure with its stage and indices. Otherwise return the finite candidate $\widehat{p}=-\widehat{V}\mathrm{diag}({\widehat{\sigma }}_1^{-1},\dots ,{\widehat{\sigma }}_n^{-1}){\widehat{U}}_1^{T}r$ with status computed-candidate and make no accuracy claim beyond the supplied report. The displayed formula names the result of the preceding schedule; it is not a second computation. The certified intervals determine only the threshold branch; no computed-SVD accuracy bound is asserted. A computed SVD routine is not defined here; the SVD branch consumes its report.

If the selected factorization reports failure, return factorization-failure; every status other than exact-step is explicitly not a certified exact full-rank Gauss-Newton step.

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.

Let $m,n\ge 1$, and let $r:U\to {\mathbb{R}}^m$ be $C^2$ on an open $U\subseteq {\mathbb{R}}^n$, let $\varphi =\mathrm{\parallel }r{\mathrm{\parallel }}_2^2\mathrm{/}2$, $J=Dr$, $A=J^{T}J$, and $S(x)={\sum }_i{r}_i(x){\mathrm{\nabla }}^2{r}_i(x)$. Let $x_{\ast }\in U$ satisfy $J(x_{\ast }{)}^{T}r(x_{\ast })=0$, and let the closed ball $B=\{x:\mathrm{\parallel }x-x_{\ast }{\mathrm{\parallel }}_2\le R\}$ lie in $U$, with $R>0$. Suppose on $B$:

• $J(x)$ has full column rank and smallest singular value at least $\sigma >0$;
• $\mathrm{\parallel }J(x){\mathrm{\parallel }}_2\le M$;
• $J$ is $L$-Lipschitz in the Euclidean induced matrix norm;
• $\mathrm{\parallel }S(x){\mathrm{\parallel }}_2\le \kappa$.

Define $q=(\kappa +MLR)\mathrm{/}{\sigma }^2$ and assume $q<1$. For the exact full-rank Gauss-Newton step $p(x)=-A(x{)}^{-1}J(x{)}^{T}r(x)$, supply an exact reduced QR factorization of $J(x)$ to the exact branch of Gauss-Newton linearized least-squares step with QR/SVD rank safeguards at each $x\in B$; the full-column-rank hypothesis guarantees such a factorization. Every $x\in B$ satisfies $\mathrm{\parallel }x+p(x)-x_{\ast }{\mathrm{\parallel }}_2\le q\mathrm{\parallel }x-x_{\ast }{\mathrm{\parallel }}_2\mathrm{.}$ Consequently the exact iterates from any $x_0\in B$ stay in $B$ and converge to $x_{\ast }$ at a linear rate no larger than $q$. The curvature hypothesis controls the residual term $S$; no quadratic rate is claimed when it is nonzero.

Let $m,n\ge 1$, and let $r:U\to {\mathbb{R}}^m$ be $C^1$ on an open $U\subseteq {\mathbb{R}}^n$, and let $x_{\ast }\in U$ satisfy $r(x_{\ast })=0$. Let $B=\{x:\mathrm{\parallel }x-x_{\ast }{\mathrm{\parallel }}_2\le R\}\subseteq U$ be convex. Suppose $J=Dr$ is $L$-Lipschitz on $B$ and, for some $\sigma >0$, $\mathrm{\parallel }J(x)v{\mathrm{\parallel }}_2\ge \sigma \mathrm{\parallel }v{\mathrm{\parallel }}_2(x\in B,v\in {\mathbb{R}}^n)\mathrm{.}$ Thus every $J(x)$ has full column rank. For the exact Gauss-Newton iteration $x_{k+1}=x_k+p_k$, where $p_k$ is the unique minimizer of $\mathrm{\parallel }r(x_k)+J(x_k)p{\mathrm{\parallel }}_2^2$. When the iteration uses the exact branch of Gauss-Newton linearized least-squares step with QR/SVD rank safeguards, assume each invocation is supplied with an exact reduced QR factorization of $J(x_k)$. Set $C=L\mathrm{/}(2\sigma )$. Then every $x_k\in B$ satisfies $\mathrm{\parallel }{x}_{k+1}-x_{\ast }{\mathrm{\parallel }}_2\le C\mathrm{\parallel }{x}_k-x_{\ast }{\mathrm{\parallel }}_2^2\mathrm{.}$ If $C>0$ and $\mathrm{\parallel }{x}_0-x_{\ast }{\mathrm{\parallel }}_2\le \min \{R,1\mathrm{/}(2C)\}$, all iterates stay in the ball, converge to $x_{\ast }$, and contract at least by $1\mathrm{/}2$ per iteration after the quadratic estimate. If $C=0$ and $x_0\in B$, one exact step reaches $x_{\ast }$. This result concerns exact solves; it asserts no finite-precision convergence rate.