Gauss-Newton linearized least-squares step with QR/SVD rank safeguards

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.