Newton and Quasi-Newton Methods for Nonlinear Systems: Broyden's Method

Intuition Beginner

A single equation in a single unknown, like $x^2=2$, asks where one curve crosses zero. A system asks something harder: find the one set of numbers that makes several equations all hold at once. Picture two unknowns and two equations. Each equation carves out a curve in the plane, and you want the point where both curves meet. With three unknowns you want where three surfaces all intersect. The number of equations matches the number of unknowns, so a single crossing point is what you expect to find.

Newton's idea from the one-variable case still works. Near your current guess, replace each curved equation by its flat tangent stand-in — a straight line or a flat plane. Flat objects you can solve exactly: finding where several flat planes meet is just solving a linear system, the kind with a fixed recipe. The meeting point of the flat stand-ins becomes your next guess, and near the answer the stand-ins hug the true curves so closely that each step roughly doubles your correct digits, exactly as before.

The catch is the same one the secant method fixed in one dimension. Building those flat stand-ins needs all the slopes — a whole grid of partial derivatives called the Jacobian — and recomputing that grid every step is expensive, sometimes impossible. So we carry the secant trick upward: instead of recomputing the slope grid, watch how the equations' values change across each step and nudge a running guess of the grid to match. That running-guess version is Broyden's method, the many-equation cousin of the secant method.

Visual Beginner

Think of two unknowns. Each equation is a curve in the plane, and the solution is the dot where the two curves cross. Standing at a guess, you draw the tangent line to each curve. Two straight lines cross at one point you can compute exactly, and that crossing is your next guess. Repeat, and the guesses zip toward the true crossing of the curved equations.

Broyden's method draws the same kind of crossing but never measures the exact tangents. It keeps a cheap stand-in for the slope grid and, after each step, adjusts it using only how the two equation-values changed — the many-dimensional echo of "slope equals rise over run".

method	what it uses each step	cost per step	speed near the solution
Newton for systems	exact slope grid (Jacobian) + linear solve	expensive	digits double
Broyden	a learned slope grid, nudged from value changes	cheap (no new derivatives)	digits almost double
Newton-Krylov	Jacobian only through directional probes	medium	digits double, big problems

The takeaway: solving several equations at once is "draw the flat stand-ins, find where they meet, repeat", and Broyden makes it cheap by learning the slope grid from the value changes you already measured.

Worked example Beginner

Take two equations in two unknowns:

$$F_1(x,y)=x^2+y^2-4=0,F_2(x,y)=xy-1=0.$$

The first says the point sits on a circle of radius $2$; the second says the product of the coordinates is $1$. We hunt for a crossing near $x=2$, $y=0.5$.

The slope grid (Jacobian) collects the four partial slopes:

$$J=\left[\begin{array}{cc}2x& 2y\\ y& x\end{array}\right].$$

At the guess $(2,0.5)$ this is $J=\left[\begin{array}{cc}4& 1\\ 0.5& 2\end{array}\right]$. The equation values there are $F_1=4+0.25-4=0.25$ and $F_2=1-1=0$, so $F=(0.25,0)$.

Newton's step solves the flat system $Jp=-F$, that is

$$\left[\begin{array}{cc}4& 1\\ 0.5& 2\end{array}\right]\left[\begin{array}{c}{p}_1\\ p_2\end{array}\right]=\left[\begin{array}{c}-0.25\\ 0\end{array}\right].$$

The determinant is $4\cdot 2-1\cdot 0.5=7.5$. Solving, $p_1=(-0.25\cdot 2-1\cdot 0)/7.5=-0.5/7.5=-0.0667$ and $p_2=(4\cdot 0-0.5\cdot (-0.25))/7.5=0.125/7.5=0.0167$. The new guess is $(2-0.0667,0.5+0.0167)=(1.9333,0.5167)$.

Check it improved: $F_1=1.9333^2+0.5167^2-4=3.7376+0.2670-4=0.0046$, far smaller than $0.25$, and $F_2=1.9333\cdot 0.5167-1=-0.0011$. Both values dropped toward zero.

What this tells us: solving the flat stand-in system $Jp=-F$ for the step $p$ is the whole Newton move in many dimensions. Broyden does the same move but reuses a stored slope grid instead of the exact $J$, so it skips recomputing the four partial slopes every step.

Check your understanding Beginner

Formal definition Intermediate+

Let $F:\Omega \subseteq {\mathbb{R}}^n\to {\mathbb{R}}^n$ be continuously differentiable on an open set $\Omega$ containing a root $x^{\ast }$ with $F(x^{\ast })=0$. The Jacobian is the matrix field $J(x)=F^{\prime }(x)\in {\mathbb{R}}^{n\times n}$ with entries $J(x{)}_{ij}=\partial F_i/\partial x_j$. The root $x^{\ast }$ is regular if $J(x^{\ast })$ is invertible, the multivariable analogue of the simple-root condition $f^{\prime }(x^{\ast })\ne 0$ of the scalar theory 43.02.03 and the nonsingularity hypothesis of the inverse function theorem 02.05.04.

Newton's method for the system generates iterates by linearising $F$ at $x_k$ and solving for the root of the linear model $F(x_k)+J(x_k)(x-x_k)$:

$$J(x_k)p_k=-F(x_k),x_{k+1}=x_k+p_k.$$

Each step is a linear solve with the current Jacobian — the structural cost over the scalar case, where the division by $f^{\prime }(x_n)$ is the $n=1$ instance. After each step define the secant pair

$$s_k=x_{k+1}-x_k,y_k=F(x_{k+1})-F(x_k),$$

and require of any Jacobian model $B_{k+1}\approx J(x_{k+1})$ the secant equation

$$B_{k+1}{s}_k=y_k.$$

By the integral mean-value form of Taylor's theorem 02.05.05, $y_k=({\int }_0^{1}J(x_k+t{s}_k)dt)s_k$, so the secant equation forces $B_{k+1}$ to reproduce the average Jacobian along the last step. Unlike the scalar case $B_{k+1}=y_k/s_k$, the matrix secant equation is $n$ scalar conditions on an $n^2$-entry matrix and so is badly underdetermined for $n>1$.

Broyden's method (the "good" Broyden update) resolves the underdetermination by the least-change principle: among all matrices satisfying the secant equation, choose the one closest to $B_k$ in Frobenius norm. The minimiser is the rank-one update

$$B_{k+1}=B_k+\frac{(y_k-B_k{s}_k)s_k^{\top }}{s_k^{\top }{s}_k},$$

which changes $B_k$ only along the direction $s_k$ and leaves it unaltered on the orthogonal complement $s_k^{\perp }$. The model $B_k$ is nonsymmetric in general — a Jacobian need not be symmetric — which is the central contrast with the optimization quasi-Newton updates 44.03.04, whose model is a symmetric Hessian and where BFGS additionally enforces positive-definiteness. Here no symmetry and no definiteness are imposed; the update is a single rank-one correction rather than the rank-two BFGS form. The Broyden step then solves $B_k{p}_k=-F(x_k)$.

Applying the Sherman-Morrison rank-one inverse identity to the update gives the inverse Jacobian model $H_k=B_k^{-1}$ directly, so the step costs $O(n^2)$ with no refactorisation:

$$H_{k+1}=H_k+\frac{(s_k-H_k{y}_k)s_k^{\top }{H}_k}{s_k^{\top }{H}_k{y}_k},p_k=-H_{k}F(x_k).$$

Because pure Newton and Broyden converge only locally, both are wrapped in a globalisation: a line search on the merit function $\mathfrak{m}(x)=\frac{1}{2}\Vert F(x){\Vert }_2^2$, whose gradient is $\nabla \mathfrak{m}(x)=J(x{)}^{\top }F(x)$, so the Newton direction $p_k=-J(x_k{)}^{-1}F(x_k)$ satisfies $\nabla \mathfrak{m}(x_k{)}^{\top }{p}_k=-F(x_k{)}^{\top }F(x_k)=-\Vert F(x_k){\Vert }^2<0$ and is a descent direction for $\mathfrak{m}$. The symbols $F$, $J$, $B_k$ (Jacobian model), $H_k$ (inverse model), the secant pair $s_k,y_k$, and the merit function $\mathfrak{m}$ are recorded in _meta/NOTATION.md.

Counterexamples to common slips Intermediate+

• "Broyden is just BFGS for systems." No. BFGS 44.03.04 is a symmetric rank-two update enforcing positive-definiteness of a Hessian model; Broyden is a nonsymmetric rank-one update of a Jacobian model with no definiteness. They share only the least-change secant philosophy. Using BFGS for a general $F(x)=0$ would impose a symmetry the true Jacobian does not have.

• "The matrix secant equation determines the update." For $n=1$ it does ($B_{k+1}=y_k/s_k$), but for $n>1$ the secant equation is $n$ conditions on $n^2$ entries and leaves an $n^2-n$ dimensional affine family. Broyden selects from that family by the Frobenius least-change projection; without such a selection the model can drift away from the Jacobian.

• "Broyden converges quadratically like Newton." It converges q-superlinearly, not quadratically, because $B_k$ only approximates $J(x^{\ast })$. The rate is governed by the Dennis-Moré condition along the steps, not by a Lipschitz-Jacobian square law.

• "The merit function $\frac{1}{2}\Vert F{\Vert }^2$ has the same minimisers as $F(x)=0$." Every root of $F$ is a global minimiser of $\mathfrak{m}$, but $\mathfrak{m}$ can have spurious stationary points where $J(x{)}^{\top }F(x)=0$ with $F(x)\ne 0$ (a singular Jacobian away from any root). A line search can stall there, which is why globalisation needs a singularity-aware safeguard, not merit-function descent alone.

Key theorem with proof Intermediate+

The signature structural result is that the Broyden update is exactly the least-change matrix satisfying the secant equation. This is the property that makes the nonsymmetric secant model stable: a single step's information overwrites the model only in the direction it actually measured, $s_k$, and preserves everything the model already encoded about the orthogonal directions. It is the nonsymmetric counterpart of the variational characterisation that singles out BFGS and DFP 44.03.04.

Theorem (least-change characterisation of the Broyden update). Fix $B_k\in {\mathbb{R}}^{n\times n}$ and a secant pair $s_k\ne 0$, $y_k$. Among all matrices $B\in {\mathbb{R}}^{n\times n}$ satisfying the secant equation $B{s}_k=y_k$, the unique minimiser of the Frobenius distance $\Vert B-B_k{\Vert }_F$ is

$$B_{k+1}=B_k+\frac{(y_k-B_k{s}_k)s_k^{\top }}{s_k^{\top }{s}_k}.$$

Proof. Let $\mathcal{A}=\{B:B{s}_k=y_k\}$, an affine subspace of ${\mathbb{R}}^{n\times n}$: it is nonempty (the displayed $B_{k+1}$ lies in it, verified below) and its difference space is $\mathcal{V}=\{D:D{s}_k=0\}$. Minimising $\Vert B-B_k{\Vert }_F$ over $B\in \mathcal{A}$ is the orthogonal projection of $B_k$ onto $\mathcal{A}$ in the Frobenius inner product $⟨X,Y⟩=\mathrm{tr}(X^{\top }Y)$, so the minimiser is the unique $B_{k+1}\in \mathcal{A}$ with $B_{k+1}-B_k\perp \mathcal{V}$.

Write the candidate correction $E=B_{k+1}-B_k=\frac{(y_k-B_k{s}_k)s_k^{\top }}{s_k^{\top }{s}_k}$. First, $B_{k+1}\in \mathcal{A}$:

$$B_{k+1}{s}_k=B_k{s}_k+\frac{(y_k-B_k{s}_k)(s_k^{\top }{s}_k)}{s_k^{\top }{s}_k}=B_k{s}_k+(y_k-B_k{s}_k)=y_k.$$

Second, $E\perp \mathcal{V}$. The correction has the form $E=u{s}_k^{\top }$ with $u=(y_k-B_k{s}_k)/(s_k^{\top }{s}_k)$. For any $D\in \mathcal{V}$ (so $D{s}_k=0$),

$$⟨E,D⟩=\mathrm{tr}(E^{\top }D)=\mathrm{tr}(s_k{u}^{\top }D)=u^{\top }D{s}_k=u^{\top }(D{s}_k)=0,$$

using the cyclic property of the trace and $D{s}_k=0$. Therefore $E$ is Frobenius-orthogonal to the entire difference space $\mathcal{V}$, and $B_{k+1}=B_k+E$ is the orthogonal projection of $B_k$ onto $\mathcal{A}$. Orthogonal projection onto a nonempty closed affine subspace of a finite-dimensional inner-product space is unique, so $B_{k+1}$ is the unique minimiser. $\square$

Bridge. This least-change property is the foundational reason Broyden's nonsymmetric secant model does not deteriorate uncontrollably: each update is a Frobenius projection that overwrites only the $s_k$-direction, so the error $\Vert B_k-J(x^{\ast })\Vert$ can grow by at most an amount controlled by the step sizes — the bounded-deterioration estimate that appears again in the Advanced results as the engine of the local superlinear theorem. The argument builds toward that convergence theorem exactly as the hereditary-positivity proof of 44.03.04 built toward the global BFGS theorem; here positivity is irrelevant and the projection geometry takes its place. The construction generalises the scalar secant model $B_{k+1}=y_k/s_k$ of 43.02.03, where for $n=1$ the orthogonal complement $s_k^{\perp }$ is empty and the least-change update is forced to equal the exact divided difference. The central insight is that the matrix secant equation underdetermines the model, and the Frobenius least-change projection is the canonical way to close that gap — this is exactly the nonsymmetric analogue of the weighted-Frobenius variational origin of BFGS and DFP. Putting these together, the secant equation owns the new measurement and the least-change projection owns everything else, and the bridge is the rank-one correction $u{s}_k^{\top }$ that ties the two together so the model changes as little as the new data permit.

Exercises Intermediate+

Advanced results Master

Broyden's method solves $F(x)=0$ by carrying the scalar secant model 43.02.03 to a nonsymmetric Jacobian approximation, and the theory below organises around four facts: the bounded-deterioration estimate that powers local convergence, the resulting local q-superlinear theorem with its Dennis-Moré characterisation, the inexact-Newton/Newton-Krylov framework that solves the linearised step approximately, and the global structure supplied by merit-function line search and homotopy continuation.

Theorem 1 (least-change update and bounded deterioration). Among all matrices satisfying the secant equation $B{s}_k=y_k$, the good-Broyden update is the unique Frobenius-nearest matrix to $B_k$, namely $B_{k+1}=B_k+(y_k-B_k{s}_k)s_k^{\top }/(s_k^{\top }{s}_k)$. Writing $E_k=B_k-J(x^{\ast })$, the update obeys the bounded-deterioration estimate

$$\Vert E_{k+1}{\Vert }_F\le \Vert E_k{\Vert }_F+\frac{\Vert ({\overline{J}}_k-J(x^{\ast }))s_k\Vert }{\Vert s_k\Vert }\le \Vert E_k{\Vert }_F+\gamma \max (\Vert e_{k+1}\Vert ,\Vert e_k\Vert ),$$

where ${\overline{J}}_k={\int }_0^{1}J(x_k+t{s}_k)dt$ and $\gamma$ depends on the Jacobian's Lipschitz constant. The model error can grow, but only by a quantity that is summable along a convergent sequence, which is exactly what keeps $\{B_k\}$ uniformly invertible near a regular root ^{[Dennis, J. E. & Schnabel, R. B. — Numerical Methods for Unconstrained Optimization and Nonlinear Equations]}.

Theorem 2 (local q-superlinear convergence; Dennis-Moré). Let $F$ be $C^1$ with $J$ Lipschitz near a regular root $x^{\ast }$ ($J(x^{\ast })$ invertible). There exist $\epsilon ,\delta >0$ such that if $\Vert x_0-x^{\ast }\Vert \le \epsilon$ and $\Vert B_0-J(x^{\ast })\Vert \le \delta$, the Broyden iterates are well defined, remain in the ball, and converge to $x^{\ast }$. The convergence is q-superlinear, and it is q-superlinear if and only if

$$\underset{k\to \infty }{\lim }\frac{\Vert (B_k-J(x^{\ast }))s_k\Vert }{\Vert s_k\Vert }=0.$$

The condition requires only that $B_k$ act like the true Jacobian along the step directions, not that $B_k\to J(x^{\ast })$; Broyden's models generically do not converge to $J(x^{\ast })$, yet attain the superlinear rate. This is the nonsymmetric counterpart of the Dennis-Moré theorem for optimization quasi-Newton 44.03.04, with the symmetric Hessian replaced by the Jacobian and positive-definiteness dropped entirely ^{[Dennis, J. E. & Schnabel, R. B. — Numerical Methods for Unconstrained Optimization and Nonlinear Equations]}.

Theorem 3 (inexact-Newton and Newton-Krylov). Solving the Newton equation $J(x_k)p_k=-F(x_k)$ exactly is unnecessary for fast convergence. An inexact-Newton method accepts any step satisfying the residual bound

$$\Vert J(x_k)p_k+F(x_k)\Vert \le {\eta }_k\Vert F(x_k)\Vert ,{\eta }_k\in [0,1),$$

with forcing terms ${\eta }_k$. Local convergence is linear if ${\eta }_k\le {\eta }_{\max }<1$, q-superlinear if ${\eta }_k\to 0$, and quadratic if ${\eta }_k=O(\Vert F(x_k)\Vert )$. The practical realisation is Newton-Krylov: solve the linear Newton equation by a Krylov method — typically GMRES 43.07.03 for the nonsymmetric Jacobian — using only Jacobian-vector products, which can be obtained matrix-free by the finite-difference probe $J(x)v\approx (F(x+ϵv)-F(x))/ϵ$. The Krylov solver's stopping tolerance is the forcing term ${\eta }_k$, so the inexact-Newton theory directly tunes the inner iteration count against the outer convergence order. Newton-Krylov is the method of choice when $n$ is large and forming or factorising $J$ is infeasible ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 4 (merit-function globalisation and its limits). A line search on the merit function $\mathfrak{m}(x)=\frac{1}{2}\Vert F(x){\Vert }_2^2$ globalises Newton and Broyden: since $\nabla \mathfrak{m}=J^{\top }F$, the Newton direction satisfies $\nabla {\mathfrak{m}}^{\top }p=-\Vert F{\Vert }^2<0$, so a backtracking line search enforcing the Armijo condition on $\mathfrak{m}$ produces a sequence with $\nabla \mathfrak{m}(x_k)\to 0$. Every limit point therefore satisfies $J(x{)}^{\top }F(x)=0$. At a regular point this forces $F(x)=0$, a root; but at a point where $J(x)$ is singular the line search can converge to a spurious stationary point of $\mathfrak{m}$ with $F(x)\ne 0$. This is the structural cost of the least-squares reformulation, and it is why robust solvers combine the merit-function search with a trust-region 44.03.02 safeguard or a restart, rather than relying on descent alone ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Theorem 5 (continuation/homotopy as globalisation). When no starting point is close enough for local convergence, embed $F$ in a homotopy $H:{\mathbb{R}}^n\times [0,1]\to {\mathbb{R}}^n$ with $H(x,0)=G(x)$ for an easily solved $G$ (often $G(x)=x-x_0$) and $H(x,1)=F(x)$ — for instance the convex homotopy $H(x,\lambda )=\lambda F(x)+(1-\lambda )(x-x_0)$. If the solution set $\{(x,\lambda ):H(x,\lambda )=0\}$ is a smooth curve with ${\partial }_{x}H$ invertible along it (a transversality condition holding for almost every $x_0$ by Sard's theorem), one tracks the curve from $\lambda =0$ to $\lambda =1$ by a predictor-corrector continuation, taking a Newton or Broyden corrector at each $\lambda$-increment. Continuation converts a hard global problem into a sequence of locally easy ones, paying with the cost of the path and the need to handle turning points where ${\partial }_{x}H$ becomes singular ^{[Nocedal, J. & Wright, S. J. — Numerical Optimization (2nd ed.)]}.

Synthesis. Solving $F(x)=0$ in ${\mathbb{R}}^n$ is one construction repeated — linearise, solve, step — and its entire theory descends from the regularity hypothesis $J(x^{\ast })$ invertible, which is the foundational reason the local solution branch exists and the multivariable analogue of the simple-root condition $f^{\prime }(x^{\ast })\ne 0$ of the scalar theory 43.02.03. Newton's method takes the exact Jacobian and attains the quadratic rate; Broyden replaces it by the least-change nonsymmetric secant model and attains the superlinear rate characterised by the Dennis-Moré condition along the steps — this is exactly the relaxation that the optimization quasi-Newton units 44.03.04 perform for a symmetric Hessian, and the central insight unifying the two is that a secant model need match the true derivative only in the directions actually stepped. The bounded-deterioration estimate is dual to the hereditary positive-definiteness of BFGS: where BFGS preserves a positivity invariant, Broyden preserves a near-Jacobian invariant, and each is what keeps the model usable across the run.

Putting these together, the inexact-Newton framework reveals that even the linear solve need only be approximate, so Newton-Krylov with GMRES 43.07.03 generalises the whole scheme to large matrix-free problems, and the merit function $\frac{1}{2}\Vert F{\Vert }^2$ with continuation/homotopy supplies the globalisation that the locally-convergent core lacks. The bridge forward is explicit: the nonsymmetric secant model here is the sibling of the symmetric Hessian models that drive every large-scale smooth optimizer, the inexact-Newton tolerance is the same forcing-term idea that tunes interior-point and SQP inner solves, and the homotopy view is the entry point to predictor-corrector path-following throughout numerical continuation.

Full proof set Master

Proposition 1 (least-change characterisation of the Broyden update). Given $B_k$ and a secant pair $s_k\ne 0$, $y_k$, the unique minimiser of $\Vert B-B_k{\Vert }_F$ subject to $B{s}_k=y_k$ is $B_{k+1}=B_k+(y_k-B_k{s}_k)s_k^{\top }/(s_k^{\top }{s}_k)$.

Proof. As in the Key theorem. The constraint set $\mathcal{A}=\{B:B{s}_k=y_k\}$ is a nonempty affine subspace with difference space $\mathcal{V}=\{D:D{s}_k=0\}$. The Frobenius-nearest point of $\mathcal{A}$ to $B_k$ is the orthogonal projection, characterised by $B_{k+1}\in \mathcal{A}$ and $B_{k+1}-B_k\perp \mathcal{V}$. The displayed correction $E=u{s}_k^{\top }$ with $u=(y_k-B_k{s}_k)/(s_k^{\top }{s}_k)$ gives $B_{k+1}{s}_k=y_k$, and for any $D\in \mathcal{V}$, $⟨E,D⟩=\mathrm{tr}(s_k{u}^{\top }D)=u^{\top }D{s}_k=0$. Uniqueness of orthogonal projection onto a closed affine subspace finishes. $\square$

Proposition 2 (Sherman-Morrison inverse update). If $B_k$ is invertible with $H_k=B_k^{-1}$ and $s_k^{\top }{H}_k{y}_k\ne 0$, the Broyden-updated $B_{k+1}=B_k+(y_k-B_k{s}_k)s_k^{\top }/(s_k^{\top }{s}_k)$ is invertible with inverse $H_{k+1}=H_k+(s_k-H_k{y}_k)s_k^{\top }{H}_k/(s_k^{\top }{H}_k{y}_k)$.

Proof. Write the update as a rank-one perturbation $B_{k+1}=B_k+u{v}^{\top }$ with $u=(y_k-B_k{s}_k)/(s_k^{\top }{s}_k)$ and $v=s_k$. The Sherman-Morrison formula states that $B_k+u{v}^{\top }$ is invertible iff $1+v^{\top }{B}_k^{-1}u\ne 0$, with inverse $H_k-\frac{H_{k}u{v}^{\top }{H}_k}{1+v^{\top }{H}_{k}u}$. Compute $H_{k}u=H_k(y_k-B_k{s}_k)/(s_k^{\top }{s}_k)=(H_k{y}_k-s_k)/(s_k^{\top }{s}_k)$ and $v^{\top }{H}_{k}u=s_k^{\top }(H_k{y}_k-s_k)/(s_k^{\top }{s}_k)=(s_k^{\top }{H}_k{y}_k-s_k^{\top }{s}_k)/(s_k^{\top }{s}_k)$, so $1+v^{\top }{H}_{k}u=s_k^{\top }{H}_k{y}_k/(s_k^{\top }{s}_k)$, nonzero by hypothesis. Substituting,

$$H_{k+1}=H_k-\frac{(H_k{y}_k-s_k)/(s_k^{\top }{s}_k)s_k^{\top }{H}_k}{s_k^{\top }{H}_k{y}_k/(s_k^{\top }{s}_k)}=H_k-\frac{(H_k{y}_k-s_k)s_k^{\top }{H}_k}{s_k^{\top }{H}_k{y}_k}=H_k+\frac{(s_k-H_k{y}_k)s_k^{\top }{H}_k}{s_k^{\top }{H}_k{y}_k}.$$

This is the stated inverse update, so each Broyden step costs $O(n^2)$ (a matrix-vector product and a rank-one update) with no factorisation. $\square$

Proposition 3 (merit-function descent). For $\mathfrak{m}(x)=\frac{1}{2}\Vert F(x){\Vert }_2^2$ one has $\nabla \mathfrak{m}(x)=J(x{)}^{\top }F(x)$, and if $J(x)$ is invertible and $F(x)\ne 0$ the Newton direction $p=-J(x{)}^{-1}F(x)$ satisfies $\nabla \mathfrak{m}(x{)}^{\top }p=-\Vert F(x){\Vert }^2<0$.

Proof. Differentiate $\mathfrak{m}=\frac{1}{2}{F}^{\top }F$: by the product and chain rules $\nabla \mathfrak{m}(x)=(F^{\prime }(x){)}^{\top }F(x)=J(x{)}^{\top }F(x)$. Then $\nabla \mathfrak{m}(x{)}^{\top }p=F^{\top }J(-J^{-1}F)=-F^{\top }F=-\Vert F{\Vert }^2$, which is strictly negative when $F(x)\ne 0$. A point where $\nabla \mathfrak{m}(x)=0$ has $J(x{)}^{\top }F(x)=0$; if $J(x)$ is invertible this forces $F(x)=0$, but a singular $J(x)$ admits $F(x)\ne 0$ with $J(x{)}^{\top }F(x)=0$, the spurious stationary point. $\square$

Proposition 4 (bounded deterioration). With $E_k=B_k-J(x^{\ast })$ and $J$ Lipschitz of constant $L$ near a regular root, the Broyden update satisfies $\Vert E_{k+1}{\Vert }_F\le \Vert E_k{\Vert }_F+\frac{L}{2}\max (\Vert e_{k+1}\Vert ,\Vert e_k\Vert )$ for some adjusted constant absorbed into $L$.

Proof. Subtract $J(x^{\ast })$ from the update. Since $J(x^{\ast })s_k=J(x^{\ast })s_k$ and the secant pair satisfies $y_k={\overline{J}}_k{s}_k$ for ${\overline{J}}_k={\int }_0^{1}J(x_k+t{s}_k)dt$,

$$E_{k+1}=B_{k+1}-J(x^{\ast })=E_k+\frac{(y_k-B_k{s}_k)s_k^{\top }}{s_k^{\top }{s}_k}=E_k(I-\frac{s_k{s}_k^{\top }}{s_k^{\top }{s}_k})+\frac{({\overline{J}}_k-J(x^{\ast }))s_k{s}_k^{\top }}{s_k^{\top }{s}_k},$$

using $y_k-B_k{s}_k=({\overline{J}}_k-J(x^{\ast }))s_k-E_k{s}_k$ and regrouping. The matrix $P=s_k{s}_k^{\top }/(s_k^{\top }{s}_k)$ is the orthogonal projector onto $\mathrm{span}(s_k)$, so $\Vert I-P{\Vert }_2=1$ (a projector) and $\Vert E_k(I-P){\Vert }_F\le \Vert E_k{\Vert }_F$. The second term has Frobenius norm $\le \Vert ({\overline{J}}_k-J(x^{\ast }))s_k\Vert /\Vert s_k\Vert$ (a rank-one matrix $w{s}_k^{\top }/\Vert s_k{\Vert }^2$ has Frobenius norm $\Vert w\Vert /\Vert s_k\Vert$ with $w=({\overline{J}}_k-J(x^{\ast }))s_k$). The Lipschitz bound gives $\Vert {\overline{J}}_k-J(x^{\ast })\Vert \le L\max (\Vert e_{k+1}\Vert ,\Vert e_k\Vert )$ on the segment, so $\Vert ({\overline{J}}_k-J(x^{\ast }))s_k\Vert /\Vert s_k\Vert \le L\max (\Vert e_{k+1}\Vert ,\Vert e_k\Vert )$. Combining the two bounds gives the stated deterioration estimate. $\square$

Proposition 5 (local q-superlinear convergence under the Dennis-Moré condition). Under $x_k\to x^{\ast }$, $J(x^{\ast })$ invertible, $J$ Lipschitz near $x^{\ast }$, and unit Broyden steps $B_k{s}_k=-F(x_k)$ with $\{B_k\}$ uniformly invertible, convergence is q-superlinear iff ${\lim }_k\Vert (B_k-J(x^{\ast }))s_k\Vert /\Vert s_k\Vert =0$.

Proof. Write $J^{\ast }=J(x^{\ast })$, $e_k=x_k-x^{\ast }$. The unit step gives $B_k{s}_k=-F(x_k)$, and $y_k={\overline{J}}_k{s}_k$ with ${\overline{J}}_k\to J^{\ast }$. From $F(x_{k+1})=F(x_k)+y_k=-B_k{s}_k+{\overline{J}}_k{s}_k=({\overline{J}}_k-B_k)s_k$ and the expansion $F(x_{k+1})=J^{\ast }{e}_{k+1}+o(\Vert e_{k+1}\Vert )$ 02.05.05,

$$J^{\ast }{e}_{k+1}=({\overline{J}}_k-B_k)s_k+o(\Vert e_{k+1}\Vert )=(J^{\ast }-B_k)s_k+({\overline{J}}_k-J^{\ast })s_k+o(\Vert e_{k+1}\Vert ).$$

The Lipschitz bound makes $\Vert ({\overline{J}}_k-J^{\ast })s_k\Vert \le L\max (\Vert e_{k+1}\Vert ,\Vert e_k\Vert )\Vert s_k\Vert =o(\Vert s_k\Vert )$, so

$$\Vert J^{\ast }{e}_{k+1}\Vert \le \Vert (B_k-J^{\ast })s_k\Vert +o(\Vert s_k\Vert )+o(\Vert e_{k+1}\Vert ).$$

Near the root $\Vert s_k\Vert$ and $\Vert e_k\Vert$ are comparable (the step and error are of the same order once convergence sets in), so dividing by $\Vert e_k\Vert$ and using invertibility of $J^{\ast }$ gives $\Vert e_{k+1}\Vert /\Vert e_k\Vert \to 0$ exactly when $\Vert (B_k-J^{\ast })s_k\Vert /\Vert s_k\Vert \to 0$. The converse runs the identity backward: $\Vert e_{k+1}\Vert =o(\Vert e_k\Vert )$ forces $\Vert J^{\ast }{e}_{k+1}\Vert =o(\Vert s_k\Vert )$, whence $\Vert (B_k-J^{\ast })s_k\Vert /\Vert s_k\Vert \to 0$. The bounded-deterioration estimate of Proposition 4, summed, keeps $\{B_k\}$ uniformly invertible and supplies the hypotheses; combined with the Dennis-Moré limit it yields the q-superlinear rate. $\square$

Connections Master

• Newton and the secant method for equations 43.02.03 is the parent: this unit is its ${\mathbb{R}}^n$ continuation, where the scalar division by $f^{\prime }(x_n)$ becomes the Jacobian linear solve $J(x_k)p_k=-F(x_k)$, the simple-root condition $f^{\prime }(x^{\ast })\ne 0$ becomes the regularity $J(x^{\ast })$ invertible, and the scalar secant model $B_{k+1}=y_k/s_k$ becomes the least-change Broyden update. The scalar quadratic-convergence proof re-runs with the integral Taylor remainder and a Jacobian-Lipschitz constant, exactly as the systems theorem of that unit anticipated; Broyden is the multivariable secant method, superlinear where the scalar secant attained the golden-ratio order.

• Optimization quasi-Newton methods (BFGS, DFP, SR1) 44.03.04 are the symmetric siblings: both families resolve an underdetermined secant equation by a least-change projection, but the optimization model is a symmetric Hessian with positive-definiteness enforced (BFGS) or indefiniteness tolerated (SR1), while Broyden's model is a nonsymmetric Jacobian with neither constraint. The Dennis-Moré superlinear characterisation and the bounded-deterioration analysis are shared machinery applied to the two settings; the contrast is exactly symmetry-plus-definiteness (minimisation) versus neither (general rootfinding).

• GMRES and Krylov subspace methods 43.07.03 supply the inner linear solver for Newton-Krylov: the nonsymmetric Newton equation $J(x_k)p_k=-F(x_k)$ is solved approximately by GMRES using only Jacobian-vector products, obtainable matrix-free by a finite-difference probe, with the GMRES stopping tolerance playing the role of the inexact-Newton forcing term ${\eta }_k$. The inexact-Newton theory ties that tolerance to the attained order (linear, superlinear, or quadratic), so the Krylov solver and the outer Newton iteration are a matched pair, the large-scale realisation of the linearise-and-solve idea.

• Trust-region methods 44.03.02 and the inverse and implicit function theorems 02.05.04 are the safeguard and the foundation: the merit function $\frac{1}{2}\Vert F{\Vert }^2$ can stall at a singular-Jacobian stationary point, and a trust-region globalisation handles that failure where a pure line search cannot, while the inverse-function nonsingularity $\det J(x^{\ast })\ne 0$ is the hypothesis that makes the local solution branch well defined. Multivariable Taylor expansion 02.05.05 provides the average-Jacobian identity $y_k=({\int }_0^{1}J(x_k+t{s}_k)dt)s_k$ that interprets the secant equation and powers both the bounded-deterioration and the superlinear-convergence proofs.

Historical & philosophical context Master

The rank-one secant update for nonlinear systems is due to Charles G. Broyden, whose 1965 paper A class of methods for solving nonlinear simultaneous equations introduced both the "good" update used here and a "bad" variant updating the inverse directly, and observed that the Frobenius least-change property singles out the good update ^{[Broyden 1965]}. Broyden's name later attached also to the symmetric optimization update through the BFGS quartet, but the original Broyden method is the nonsymmetric Jacobian scheme of this unit. The local superlinear-convergence theory was established by John Dennis and Jorge Moré in the early 1970s, whose bounded-deterioration framework and superlinear characterisation cover both the Broyden and the symmetric secant families ^{[Dennis-More 1974]}, and was given its definitive textbook treatment by Dennis and Schnabel.

The inexact-Newton framework was formalised by Ron Dembo, Stanley Eisenstat, and Trond Steihaug in 1982, who proved the forcing-term characterisation of the convergence order ^{[Dembo-Eisenstat-Steihaug 1982]}, and the matrix-free Newton-Krylov realisation with GMRES was developed through the 1980s and 1990s and surveyed by C. T. Kelley ^{[Kelley 1995]}. Continuation and homotopy methods for systems trace to the degree-theoretic and path-following work of Eugene Allgower and Kurt Georg and the probability-one homotopy of Layne Watson; the transversality argument that the tracked path is generically smooth rests on Sard's theorem ^{[Allgower-Georg 1990]}.

Bibliography Master

@article{broyden1965,
  author  = {Broyden, Charles G.},
  title   = {A class of methods for solving nonlinear simultaneous equations},
  journal = {Mathematics of Computation},
  volume  = {19},
  number  = {92},
  year    = {1965},
  pages   = {577--593}
}

@book{dennisschnabel1996b,
  author    = {Dennis, John E. and Schnabel, Robert B.},
  title     = {Numerical Methods for Unconstrained Optimization and Nonlinear Equations},
  series    = {Classics in Applied Mathematics},
  volume    = {16},
  publisher = {SIAM},
  year      = {1996},
  note      = {Reprint of the 1983 Prentice-Hall edition}
}

@article{dennismore1974,
  author  = {Dennis, John E. and Mor\'{e}, Jorge J.},
  title   = {A characterization of superlinear convergence and its application to quasi-Newton methods},
  journal = {Mathematics of Computation},
  volume  = {28},
  number  = {126},
  year    = {1974},
  pages   = {549--560}
}

@article{demboeisenstatsteihaug1982,
  author  = {Dembo, Ron S. and Eisenstat, Stanley C. and Steihaug, Trond},
  title   = {Inexact Newton methods},
  journal = {SIAM Journal on Numerical Analysis},
  volume  = {19},
  number  = {2},
  year    = {1982},
  pages   = {400--408}
}

@book{kelley1995,
  author    = {Kelley, C. T.},
  title     = {Iterative Methods for Linear and Nonlinear Equations},
  series    = {Frontiers in Applied Mathematics},
  volume    = {16},
  publisher = {SIAM},
  year      = {1995}
}

@book{nocedalwright2006b,
  author    = {Nocedal, Jorge and Wright, Stephen J.},
  title     = {Numerical Optimization},
  edition   = {2},
  series    = {Springer Series in Operations Research and Financial Engineering},
  publisher = {Springer},
  year      = {2006}
}

@book{allgowergeorg1990,
  author    = {Allgower, Eugene L. and Georg, Kurt},
  title     = {Numerical Continuation Methods: An Introduction},
  series    = {Springer Series in Computational Mathematics},
  volume    = {13},
  publisher = {Springer},
  year      = {1990}
}