Newton, quasi-Newton and Gauss–Newton methods

Before this lesson: lesson 02 (condition number, strong convexity), lesson 03 (rates of gradient descent, Armijo and Wolfe conditions), linear algebra 07 (positive definite matrices, Cholesky), linear algebra 12 (least squares), calculus 03 (Jacobian, Hessian, Taylor's theorem). Time: about 95 minutes, plus the notebook. This is a difficult lesson; take a break after section 3.

---

Why this matters #

Gradient descent needs about $\kappa\log(1/\varepsilon)$ steps (lesson 03). Second-order methods remove most of this dependence on the condition number $\kappa$. They are the standard tools when the number of parameters is moderate and the gradients are exact (full batch):

• logistic regression and other generalized linear models (scikit-learn's LogisticRegression uses L-BFGS by default; Newton's method for logistic regression is iteratively reweighted least squares);
• fitting physical models to data (scipy.optimize.curve_fit uses Levenberg–Marquardt);
• inverse problems with partial differential equations (PDEs), physics-informed neural networks (PINNs), and weather data assimilation, where the main method (incremental 4D-Var) is a Gauss–Newton method.

In deep learning they come back as approximations. The Gauss–Newton and Fisher matrices of this lesson are the starting point of K-FAC and Shampoo (lesson 12). "Why don't we train neural networks with Newton's method?" is also a standard interview question.

1. Predict first #

2. The picture #

Notation for this lesson: $x^\star$ is a minimizer and $e_k = x_k - x^\star$ is the error. We write the Hessian as $\nabla^2 f(x)$. The letters $B_k$ and $H_k$ are reserved for quasi-Newton methods: $B_k$ approximates the Hessian and $H_k$ approximates its inverse (as in Nocedal and Wright).

Newton in one dimension: fit a parabola, jump to its bottom #

At the current point $x_k$, Taylor's theorem gives a parabola with the same value, slope and curvature as $f$:

If $f''(x_k) > 0$, the parabola opens upward. Its minimum is where $m_k'(x) = f'(x_k) + f''(x_k)(x - x_k) = 0$:

This is also Newton's root-finding method applied to the equation $f'(x) = 0$: replace $f'$ by its tangent line and find where the line crosses zero. If $f$ is itself a parabola, the model is exact and one step lands on the minimum.

Interactive figure: loading…

Newton on a quadratic: the perfect preconditioner #

Take $f(x) = \tfrac12 x^\top A x$ with $A$ symmetric positive definite. Then $\nabla f(x) = Ax$ and $\nabla^2 f(x) = A$, so the Newton step from any $x$ is

One step, for any condition number. To see why, write the error in the eigenvectors $q_i$ of $A$ (eigenvalues $\lambda_i$), as in lesson 03. Gradient descent multiplies the component along $q_i$ by $1 - \eta\lambda_i$: one step size $\eta$ must serve all directions. Newton multiplies the gradient component along $q_i$ by $1/\lambda_i$: a large step in a flat direction and a small step in a steep one. Every component becomes exactly 0.

A second way to see it: change variables to $z = A^{1/2}x$. Then $f = \tfrac12\|z\|^2$, whose level sets are circles, and gradient descent with $\eta = 1$ jumps to the centre. Mapped back to $x$, this is exactly the Newton step. Newton's method is gradient descent in the coordinates where the function is round.

Interactive figure: loading…

The race on 2D surfaces #

On a non-quadratic function the model is only correct near $x_k$, so Newton needs several steps. But near a minimizer the model gets better and better, and the error shrinks very fast. The figure compares four methods. Damped Newton shortens the step with a line search when the full step does not decrease $f$ enough. BFGS (section 3) does not use the Hessian at all; it learns it from gradients.

Interactive figure: loading…

3. Definitions and main results #

The Newton step #

The quadratic model of $f$ at $x$ is

If $\nabla^2 f(x) \succ 0$, the model is strictly convex and its minimizer solves the linear system

$\Delta x_{\rm nt}$ is the Newton step, and pure Newton sets $x_{k+1} = x_k + \Delta x_{\rm nt}$. Solve the system with a Cholesky factorization $\nabla^2 f = LL^\top$ (about $d^3/3$ operations; see linear algebra 07). Never form the inverse matrix. If Cholesky fails, the Hessian is not positive definite, and that is useful information.

The number $\delta(x) = \sqrt{\nabla f(x)^\top [\nabla^2 f(x)]^{-1} \nabla f(x)}$ is the Newton decrement (Boyd and Vandenberghe write $\lambda(x)$). The model predicts the decrease $f(x) - \min m = \delta(x)^2/2$. Stopping when $\delta(x)^2/2 < \varepsilon$ is the standard rule.

Local quadratic convergence #

Why this is true (intuition). The Newton step is exact for a quadratic. A smooth function differs from its quadratic model by a term of order $\|e\|^3$ (Taylor), so after division by the Hessian the new error is of order $\|e\|^2$. Squaring a small number makes it much smaller.

Theorem 1 (local quadratic convergence). Let $f$ be twice differentiable with an $M$-Lipschitz Hessian, $\|\nabla^2 f(x) - \nabla^2 f(y)\| \le M\|x - y\|$ (spectral norm). Let $\nabla f(x^\star) = 0$ and $\nabla^2 f(x^\star) \succeq \mu I$ with $\mu > 0$. If $\|e_0\| \le \mu/(2M)$, then pure Newton satisfies, for every $k$,

So $q_k = (M/\mu)\|e_k\|$ satisfies $q_{k+1} \le q_k^2$ and $q_0 \le \tfrac12$: $q_k \le (1/2)^{2^k}$.

Proof of Theorem 1 (6 steps)

1. Since $\nabla f(x^\star) = 0$, the Newton step gives $e_{k+1} = e_k - [\nabla^2 f(x_k)]^{-1}\big(\nabla f(x_k) - \nabla f(x^\star)\big)$.
2. Along the segment from $x^\star$ to $x_k = x^\star + e_k$ (fundamental theorem of calculus): $\nabla f(x_k) - \nabla f(x^\star) = \int_0^1 \nabla^2 f(x^\star + t e_k)\, e_k\, dt$.
3. Put 2 into 1 and write $e_k = \int_0^1 e_k\,dt$: $e_{k+1} = [\nabla^2 f(x_k)]^{-1} \int_0^1 \big[\nabla^2 f(x_k) - \nabla^2 f(x^\star + t e_k)\big]\, e_k\, dt$.
4. The two points $x_k$ and $x^\star + te_k$ are $(1 - t)\|e_k\|$ apart. By the Lipschitz condition, the bracket has norm at most $M(1 - t)\|e_k\|$.
5. So $\|e_{k+1}\| \le \|[\nabla^2 f(x_k)]^{-1}\|\; M\|e_k\|^2 \int_0^1 (1 - t)\,dt = \dfrac{M}{2\,\lambda_{\min}(\nabla^2 f(x_k))}\,\|e_k\|^2$.
6. If $\|e_k\| \le \mu/(2M)$, Weyl's inequality (linear algebra 17) gives $\lambda_{\min}(\nabla^2 f(x_k)) \ge \mu - \|\nabla^2 f(x_k) - \nabla^2 f(x^\star)\| \ge \mu - M\|e_k\| \ge \mu/2$. Then step 5 gives $\|e_{k+1}\| \le (M/\mu)\|e_k\|^2$, and $(M/\mu)\|e_k\| \le \tfrac12$ gives $\|e_{k+1}\| \le \tfrac12\|e_k\|$. So $x_{k+1}$ is again in the ball, and induction finishes the proof. $\blacksquare$

In words: the number of correct digits roughly doubles at every step, once you are close enough. The theorem is local. It says nothing about a start far away (worked example (b) shows a bad one).

Damped Newton, and what goes wrong without positive curvature #

Damped Newton sets $x_{k+1} = x_k + t_k\,\Delta x_{\rm nt}$, where $t_k$ comes from backtracking from $t = 1$ with the Armijo condition (lesson 03). For strongly convex $f$ with a Lipschitz Hessian, Boyd and Vandenberghe (§9.5.3) prove two phases:

1. Damped phase. Far from $x^\star$, each step decreases $f$ by at least a fixed amount $\gamma > 0$.
2. Quadratic phase. Once $\|\nabla f\|$ is small enough, the line search always accepts $t = 1$, and Theorem 1 applies.

The total number of steps is at most $(f(x_0) - f^\star)/\gamma + \log_2\log_2(\varepsilon_0/\varepsilon)$. The second term is about 6 for any accuracy you can reach in double precision.

If $\nabla^2 f(x)$ is not positive definite (a nonconvex $f$), the model has no minimum. Its stationary point can be a maximum or a saddle, so $\Delta x_{\rm nt}$ can point uphill or towards a saddle point. In the 1D figure, $f''(0.3) < 0$ and Newton moves towards the maximum. In 2D, Newton is attracted to saddle points (lesson 11 studies saddles in neural networks). Common fixes:

• use $-\nabla f$ when the Newton step is not a descent direction;
• add a multiple of the identity, $\nabla^2 f + \tau I \succ 0$ (Nocedal and Wright §3.4);
• replace each eigenvalue $\lambda_i$ of the Hessian by $|\lambda_i|$ (the "saddle-free Newton" idea, used for the damped Newton path in the 2D figure);
• use a trust region: minimize the model only inside a ball where you trust it (Levenberg–Marquardt below is an example).

Affine invariance #

Let $x = Ty$ with an invertible matrix $T$, and $g(y) = f(Ty)$. The chain rule gives $\nabla g(y) = T^\top\nabla f(x)$ and $\nabla^2 g(y) = T^\top\nabla^2 f(x)\,T$. So the Newton step in $y$ is

The step in $y$ is the step in $x$, written in the new coordinates. An Armijo line search uses only values of $f$ and the slope $\nabla f^\top \Delta x$, which are also the same. So the iterates satisfy $y_k = T^{-1}x_k$: Newton's method is affine invariant. Its speed does not depend on how you scale or rotate the variables, and so not on $\kappa$. Gradient descent is not: its step in $y$ is $-T^\top\nabla f$, which is $-TT^\top\nabla f$ in $x$. This is why feature scaling matters for gradient descent and not for Newton.

Quasi-Newton methods: the secant condition and BFGS #

Computing and factorizing $\nabla^2 f$ costs $O(d^3)$ per step. But every step gives, almost for free, the pair

By Taylor, $y_k \approx \nabla^2 f\, s_k$: the change of the gradient shows how the Hessian acts on one direction. A quasi-Newton method keeps a matrix $B_k \approx \nabla^2 f$ (or $H_k \approx [\nabla^2 f]^{-1}$) and asks the new matrix to reproduce the observed change. This is the secant condition:

In 1D it fixes $B_{k+1} = \big(f'(x_{k+1}) - f'(x_k)\big)/(x_{k+1} - x_k)$: the secant method. In $d$ dimensions it gives $d$ equations for $d(d+1)/2$ unknowns, so we choose the smallest change that satisfies it.

Derivation of BFGS. Change $B$ by a symmetric rank-2 matrix built from the two vectors we have: $B^+ = B + a\, yy^\top + b\,(Bs)(Bs)^\top$. The secant condition $B^+s = y$ reads

It holds when $a\,(y^\top s) = 1$ and $b\,(s^\top Bs) = -1$. So

The Sherman–Morrison–Woodbury formula inverts this rank-2 change. With $\rho = 1/(y^\top s)$, the BFGS inverse update (Broyden, Fletcher, Goldfarb, Shanno, 1970) is

Check of the secant condition: $(I - \rho\, ys^\top)\,y = y - \rho\,y\,(s^\top y) = 0$, so the first term vanishes on $y$ and $H^+y = \rho\, s\,(s^\top y) = s$. The method then steps along $-H_k\nabla f(x_k)$ with a line search. Each update costs $O(d^2)$ and no linear system is solved.

Proposition 2 (BFGS keeps positive definiteness). If $H \succ 0$ and $y^\top s > 0$, then $H^+ \succ 0$.

Proof.

1. Take $z \ne 0$ and set $w = (I - \rho\, ys^\top)\,z = z - \rho\,(s^\top z)\, y$.
2. Then $z^\top H^+ z = w^\top H w + \rho\,(s^\top z)^2$. Both terms are $\ge 0$, because $H \succ 0$ and $\rho > 0$.
3. If $s^\top z \ne 0$, the second term is $> 0$. If $s^\top z = 0$, then $w = z \ne 0$ and the first term is $> 0$. In both cases $z^\top H^+ z > 0$. $\blacksquare$

Where $y^\top s > 0$ comes from. For strictly convex $f$ it always holds. In general the Wolfe curvature condition gives it (lesson 03): $\nabla f(x_{k+1})^\top p_k \ge c_2 \nabla f(x_k)^\top p_k$ with $c_2 < 1$, so $y^\top s = \alpha_k\big(\nabla f(x_{k+1}) - \nabla f(x_k)\big)^\top p_k \ge (c_2 - 1)\,\alpha_k\nabla f(x_k)^\top p_k > 0$, because $\nabla f(x_k)^\top p_k < 0$. This is why BFGS runs with a Wolfe line search.

Two facts about BFGS (Nocedal and Wright, ch. 6). Near a minimizer with $\nabla^2 f(x^\star) \succ 0$ it converges superlinearly: $\|e_{k+1}\| \le c_k\|e_k\|$ with $c_k \to 0$. On a strictly convex quadratic with exact line searches it stops after at most $d$ steps, with $H_d = A^{-1}$.

Interactive figure: loading…

L-BFGS: BFGS with a short memory #

For $d = 10^6$, the matrix $H_k$ has $10^{12}$ entries (8 TB). L-BFGS (limited-memory BFGS) stores only the $m$ newest pairs $(s_i, y_i)$, with $\rho_i = 1/(y_i^\top s_i)$. It defines $H_k$ as the result of the BFGS updates with these pairs (oldest first), applied to $H^0 = \gamma I$, where $\gamma = s^\top y/y^\top y$ for the newest pair (a scalar estimate of the inverse curvature). The two-loop recursion computes $p = -H_k g$ without forming $H_k$:

1. $q \leftarrow g$. For $i$ from newest to oldest: $\alpha_i \leftarrow \rho_i\, s_i^\top q$, then $q \leftarrow q - \alpha_i y_i$.
2. $r \leftarrow \gamma\, q$.
3. For $i$ from oldest to newest: $\beta \leftarrow \rho_i\, y_i^\top r$, then $r \leftarrow r + s_i\,(\alpha_i - \beta)$.
4. $p = -r$.

Why the two loops compute $H_k g$

One BFGS update is $H^+ = V^\top H V + \rho\, ss^\top$ with $V = I - \rho\, ys^\top$. Unroll $m$ updates: $H_k g$ is a sum of terms of the form $V_{\text{newest}}^\top \cdots V^\top_{i+1}\,\rho_i s_i s_i^\top\,V_{i+1}\cdots V_{\text{newest}}\, g$, plus $V^\top_{\text{newest}}\cdots V^\top_{\text{oldest}}\,H^0\,V_{\text{oldest}}\cdots V_{\text{newest}}\,g$. The first loop applies the $V_i$ to $g$, newest first, and stores the scalars $\alpha_i = \rho_i s_i^\top(\cdot)$ that the terms need. The second loop applies the $V_i^\top$ in the opposite order and adds the $\rho_i s_i s_i^\top$ terms through $s_i\,\alpha_i$. Nocedal and Wright (Algorithm 7.4) give the full proof.

Each loop costs $m$ dot products and $m$ vector updates: $O(md)$ time per step and $2md$ stored numbers. Typical memories are $m = 3$ to $20$ (SciPy's L-BFGS-B uses 10, torch.optim.LBFGS keeps 100). L-BFGS is the default method for smooth, full-batch problems of any size.

Gauss–Newton for nonlinear least squares #

Many problems have the form $F(x) = \tfrac12\|r(x)\|^2 = \tfrac12\sum_i r_i(x)^2$, with residuals $r: \R^d \to \R^n$ (model minus data) and Jacobian $J = \partial r/\partial x$ of size $n \times d$. The chain rule (calculus 08) gives

Linearize the residual, not the gradient: $r(x + h) \approx r + Jh$. Then

$L$ is a linear least-squares problem (linear algebra 12). Its minimizer is the Gauss–Newton step

This is the Newton equation with the residual-curvature term $S$ dropped. Three properties:

• It needs only first derivatives of $r$.
• $J^\top J \succeq 0$ always. If $J$ has full column rank, then $h_{\rm gn}^\top\nabla F = h_{\rm gn}^\top J^\top r = -h_{\rm gn}^\top J^\top J h_{\rm gn} < 0$: the step is always a descent direction. A Newton step need not be.
• Compute it as the least-squares solution of $Jh \approx -r$ (QR or lstsq), not by forming $J^\top J$, because $\kappa(J^\top J) = \kappa(J)^2$.

When is dropping $S$ safe? Near $x^\star$, $\nabla F(x_k) \approx (J^\top J + S)\,e_k$, so

• Zero residual, $r(x^\star) = 0$: then $S(x^\star) = 0$, and convergence is quadratic, like Newton.
• Small residuals or nearly linear $r$: $S$ is small and the linear rate is fast.
• Large residuals and curved $r$: the rate is the spectral radius of $(J^\top J)^{-1}S(x^\star)$. If it is larger than 1, Gauss–Newton fails even close to $x^\star$ (worked example (d)).

Levenberg–Marquardt: damping as a trust region #

The Levenberg–Marquardt (LM) step adds a damping term (here $\lambda > 0$ is the damping parameter, not an eigenvalue):

• $\lambda \to 0$: $h_{\rm lm} \to h_{\rm gn}$, fast near a small-residual solution.
• $\lambda \to \infty$: $h_{\rm lm} \approx -J^\top r/\lambda = -\nabla F/\lambda$, a short gradient step, safe far from the solution.
• For every $\lambda > 0$ the matrix is positive definite, so $h_{\rm lm}$ is a descent direction.
• Trust-region view. $h_{\rm lm}$ minimizes the model $L(h)$ among all $h$ with $\|h\| \le \|h_{\rm lm}\|$ (Nocedal and Wright, Theorem 4.1). A large $\lambda$ means a small trusted ball; $\lambda$ is the Lagrange multiplier of the constraint (lesson 08).

Choosing $\lambda$: the gain ratio. Compare the actual decrease with the decrease the model predicted:

If $\varrho > 0$, accept the step and set $\lambda \leftarrow \lambda\max\{\tfrac13,\ 1 - (2\varrho - 1)^3\}$ and $\nu \leftarrow 2$. Otherwise reject it and set $\lambda \leftarrow \lambda\nu$, $\nu \leftarrow 2\nu$. Start with $\lambda_0 = \tau\max_i (J^\top J)_{ii}$; $\tau = 10^{-3}$ is a common choice (Madsen, Nielsen and Tingleff, Algorithm 3.16, Nielsen's rule). A good model ($\varrho$ near 1) shrinks $\lambda$; a bad model grows it. Marquardt's variant uses $\operatorname{diag}(J^\top J)$ in place of $I$, which makes the step independent of the scale of each parameter.

Interactive figure: loading…

Generalized Gauss–Newton, the Fisher matrix, and Hessian-free methods #

For a loss $\ell(z)$ of a model output $z = f(\theta)$, the chain rule gives

The generalized Gauss–Newton (GGN) matrix $G = J_f^\top\,\nabla^2\ell\,J_f$ keeps the first term. It is positive semidefinite when $\ell$ is convex in $z$. For the squared error it is $J^\top J$. When $\ell$ is the negative log-likelihood of an exponential family with natural parameters $z$ (softmax cross-entropy is one), it equals the Fisher information matrix (Martens, arXiv 1412.1193). K-FAC and Shampoo approximate this matrix layer by layer (lesson 12).

You never need the full matrix to use it. The conjugate gradient method (linear algebra 14) solves $Gp = -\nabla \ell$ with products $Gv$ only, and such products cost about two gradient evaluations (calculus 08). Stopping CG early gives truncated Newton or Newton-CG; with the GGN it is the Hessian-free optimization of Martens (2010). Lesson 11 uses the same Hessian-vector products to measure sharpness.

Why not Newton for neural networks? #

• Cost. For $d = 10^8$ parameters, the Hessian has $10^{16}$ entries, and a Cholesky factorization needs about $d^3/3 \approx 3 \cdot 10^{23}$ operations per step.
• Nonconvexity. The Hessian of a network loss has negative eigenvalues; pure Newton is attracted to saddle points.
• Noise. Minibatch curvature is noisy; a quasi-Newton pair $(s, y)$ from two different minibatches measures noise more than curvature.

What deep learning uses instead: diagonal scaling (Adam, lesson 06), structured approximations of the GGN or Fisher (K-FAC, Shampoo), and orthogonalized updates (Muon), all in lesson 12.

Answers to section 1. (1) One step in both cases. The quadratic model is exact, and Newton is affine invariant, so the substitution $x_2 = 10z_2$ changes nothing. Gradient descent with the best step would need about $\ln 10^6/\ln(1001/999) \approx 6{,}908$ steps at $\kappa = 1000$. (2) Towards the maximum $x = 0$: $f''(0.3) = 3 \cdot 0.09 - 1 = -0.73 < 0$, the model parabola opens downward, and its stationary point is its maximum. The iterates go $0.3 \to -0.074 \to 0.0008 \to \dots \to 0$. (3) Problems with zero residual at the solution, $r(x^\star) = 0$, where $S(x^\star) = 0$. If $r$ is linear, $S = 0$ everywhere and one step solves the problem.

4. Worked example #

(a) Newton on a 2D quadratic. $f = \tfrac12(x_1^2 + 10x_2^2)$ and $x_0 = (1, 1)$. Then $\nabla f(x_0) = (1, 10)$ and $\nabla^2 f = \operatorname{diag}(1, 10)$. The Newton step solves $\operatorname{diag}(1, 10)\,\Delta x = -(1, 10)$, so $\Delta x = (-1, -1)$ and $x_1 = (0, 0) = x^\star$.

(b) Newton in 1D: the quadratic rate. $f(x) = x^4/4 - x$, so $f'(x) = x^3 - 1$, $f''(x) = 3x^2$ and $x^\star = 1$. The Newton map is

From $x_0 = 2$:

$k$	$x_k$	$e_k = x_k - 1$	$e_k / e_{k-1}^2$
0	2	1
1	1.416667	0.417	0.417
2	1.110534	0.111	0.637
3	1.010637	$1.06\cdot10^{-2}$	0.871
4	1.000112	$1.12\cdot10^{-4}$	0.986
5	1.0000000124	$1.24\cdot10^{-8}$	1.000

The ratio in the last column tends to $f'''(x^\star)/(2f''(x^\star)) = 6/(2 \cdot 3) = 1$, so $e_{k+1} \approx e_k^2$: the exponent doubles at each step.

A bad start. From $x_0 = 0.1$, the curvature $f''(0.1) = 0.03$ is tiny. The pure step goes to $x_1 = (2 \cdot 0.001 + 1)/(3 \cdot 0.01) = 33.4$, where $f \approx 3.1 \cdot 10^5$ (it was $-0.1$). Far out, $x^+ \approx 2x/3$, so the iterates come back slowly (22.3, 14.8, 9.9, …) and reach $1.0001$ only at $k = 12$. An Armijo line search would reject the first jump: this is the damped phase.

(c) BFGS learns $A^{-1}$ in two steps. Same quadratic, $A = \operatorname{diag}(1, 10)$, $x_0 = (1, 1)$, $H_0 = I$, exact line searches $\alpha = g^\top g/(g^\top A g)$ for the direction $-g$.

• $k = 0$: $g_0 = (1, 10)$, $p_0 = -(1, 10)$, $\alpha_0 = 101/1001 \approx 0.1009$, so $x_1 \approx (0.8991, -0.0090)$. Then $y_0^\top s_0 \approx 10.19 > 0$, and the update gives $H_1 \approx \begin{pmatrix} 1.009 & -0.00009 \\ -0.00009 & 0.1000 \end{pmatrix}$. It is already correct in the steep direction, which the first step explored.
• $k = 1$: the step $-H_1 g_1$ with $\alpha_1 \approx 0.991$ gives $x_2 = (0, 0)$ up to rounding, and $H_2 = \operatorname{diag}(1, 0.1) = A^{-1}$.

(d) Gauss–Newton, Newton and LM with a nonzero residual (Madsen et al., Example 3.3). $r(x) = (x + 1,\ a x^2 + x - 1)$ with one parameter $x$. For $a < 1$, $x^\star = 0$ is a local minimizer with $r(x^\star) = (1, -1)$ and $F(x^\star) = 1$: the residual is not zero. At $x^\star$:

• $J = (1,\ 2ax + 1)^\top = (1, 1)^\top$, so $J^\top J = 2$;
• only $r_2$ is curved, with $r_2'' = 2a$, so $S = r_2\, r_2'' = -2a$ and $F''(0) = 2 - 2a$;
• the Gauss–Newton error factor is $-S/(J^\top J) = a$.

Three cases:

• $a = 0$: $r$ is linear, and one Gauss–Newton step solves the problem.
• $a = 0.5$, $x_0 = 0.1$: $r = (1.1, -0.895)$, $J = (1, 1.1)^\top$, $J^\top r = 1.1 - 0.9845 = 0.1155$ and $J^\top J = 2.21$. Gauss–Newton: $x_1 = 0.1 - 0.1155/2.21 = 0.0477$, then $0.0233$, $0.0115$ (factor $\to 0.5$, linear). Newton uses $F''(0.1) = 2.21 + (-0.895)(1) = 1.315$: $x_1 = 0.0122$, then $2.2 \cdot 10^{-4}$, $7 \cdot 10^{-8}$ (quadratic).
• $a = -2$: factor $-2$. From $0.1$, Gauss–Newton jumps around ($-0.303$, $0.137$, $-0.468$, …) and never converges, although it starts close to $x^\star$.

One LM step for $a = 0.5$, $x_0 = 0.1$, $\lambda = 1$: $h = -0.1155/(2.21 + 1) = -0.0360$, so $x_1 = 0.0640$. The predicted decrease is $\tfrac12 h(\lambda h - J^\top r) = \tfrac12(-0.0360)(-0.0360 - 0.1155) = 0.00273$. The actual decrease is $F(0.1) - F(0.0640) = 0.00333$. So $\varrho = 1.22 > 0$: accept the step. The new damping is $\lambda \leftarrow \max\{\tfrac13,\ 1 - (2.44 - 1)^3\} = \max\{\tfrac13, -1.99\} = \tfrac13$.

5. In code #

The code checks examples (b), (c) and (d) and shows what SciPy's quasi-Newton methods learn. Fraction gives exact rational arithmetic, so the table of (b) has no rounding.

from fractions import Fraction
import numpy as np
from scipy.optimize import minimize, least_squares, rosen, rosen_der, rosen_hess

# (b) Newton on f(x) = x^4/4 - x: the map x+ = (2x^3 + 1)/(3x^2), in exact arithmetic
x, errors = Fraction(2), []
for k in range(6):
    errors.append(abs(x - 1))
    x = (2 * x**3 + 1) / (3 * x**2)
print([f"{float(e):.2e}" for e in errors])
# ['1.00e+00', '4.17e-01', '1.11e-01', '1.06e-02', '1.12e-04', '1.24e-08']
print([round(float(errors[k + 1] / errors[k] ** 2), 3) for k in range(5)])
# [0.417, 0.637, 0.871, 0.986, 1.0]

# (c) BFGS on f = ½(x1² + 10 x2²) with exact line searches: H2 = A^{-1}
A = np.diag([1.0, 10.0])
x, H = np.array([1.0, 1.0]), np.eye(2)
for k in range(2):
    g = A @ x
    p = -H @ g
    alpha = -(g @ p) / (p @ A @ p)          # exact line search on a quadratic
    s = alpha * p
    y = A @ (x + s) - g
    rho = 1.0 / (y @ s)
    V = np.eye(2) - rho * np.outer(y, s)
    H = V.T @ H @ V + rho * np.outer(s, s)  # BFGS inverse update
    x = x + s
print(np.allclose(H, np.linalg.inv(A)), np.abs(x).max() < 1e-12)   # True True

# Quasi-Newton on Rosenbrock from (-1.2, 1)
x0 = np.array([-1.2, 1.0])
for method in ["BFGS", "L-BFGS-B", "CG"]:
    res = minimize(rosen, x0, jac=rosen_der, method=method)
    print(f"{method:8s} iterations={res.nit:3d} gradient calls={res.njev:3d}")
# BFGS 32/39, L-BFGS-B 36/44, CG (nonlinear conjugate gradient) 36/77 (SciPy 1.18)

# BFGS learns the curvature: its H at the end vs the true inverse Hessian at (1, 1)
res = minimize(rosen, x0, jac=rosen_der, method="BFGS")
print(np.round(res.hess_inv, 3))                           # [[0.502 1.003] [1.003 2.01 ]]
print(np.round(np.linalg.inv(rosen_hess(np.ones(2))), 3))  # [[0.5   1.   ] [1.    2.005]]

# (d) Gauss-Newton factor and LM, Madsen et al. Example 3.3; local minimizer x* = 0
def residual(x, a):
    return np.array([x[0] + 1, a * x[0] ** 2 + x[0] - 1])

def jac(x, a):
    return np.array([[1.0], [2 * a * x[0] + 1]])

for a in [0.0, 0.5, -2.0]:
    J, r = jac(np.zeros(1), a), residual(np.zeros(1), a)
    JtJ = (J.T @ J)[0, 0]
    F2 = JtJ + r[1] * 2 * a                 # F'' = JᵀJ + Σ r_i r_i''
    sol = least_squares(residual, [0.1], jac=jac, args=(a,), method="lm")
    print(f"a={a:4.1f}  JᵀJ={JtJ:.0f}  F''={F2:.0f}  GN factor={1 - F2 / JtJ:+.1f}  LM x={sol.x[0]:.0e}")
# a= 0.0  JᵀJ=2  F''=2  GN factor=+0.0  LM x=-1e-16    r is linear: one step
# a= 0.5  JᵀJ=2  F''=1  GN factor=+0.5  LM x=4e-05     linear rate; LM stops on its F tolerance
# a=-2.0  JᵀJ=2  F''=6  GN factor=-2.0  LM x=1e-07     Gauss-Newton diverges, LM converges

The iteration counts on Rosenbrock can change a little between SciPy versions. The final BFGS matrix is close to the true inverse Hessian: BFGS has learned the curvature from gradients alone.

6. Where it shows up #

• Logistic regression, generalized linear models, curve fitting. Newton's method for logistic regression is iteratively reweighted least squares (IRLS). scikit-learn's LogisticRegression uses L-BFGS by default. scipy.optimize.minimize uses BFGS by default (L-BFGS-B when there are bounds), and curve_fit uses Levenberg–Marquardt.
• Physics-informed neural networks and other full-batch problems. PINN losses are badly conditioned because of the differential operators. Rathore et al. (ICML 2024) found that Adam followed by L-BFGS works better than either one alone. In PyTorch: torch.optim.LBFGS(params, line_search_fn="strong_wolfe"), with full batches.
• Inverse problems with PDEs. Adjoint gradients plus L-BFGS or truncated Gauss–Newton; seismic full-waveform inversion is the classic case (Virieux and Operto, 2009).
• Weather data assimilation. Strong-constraint 4D-Var finds the initial state $x_0$ that minimizes $\tfrac12\|x_0 - x_b\|^2_{B^{-1}} + \tfrac12\sum_k \|y_k - \mathcal{H}_k(\mathcal{M}_k(x_0))\|^2_{R_k^{-1}}$ ($x_b$ a prior "background" state, $\mathcal{M}_k$ the forecast model up to time $t_k$, $\mathcal{H}_k$ the observation operator, $y_k$ the observations, $B$ and $R_k$ error covariances). This is nonlinear least squares. Incremental 4D-Var linearizes the model in an outer loop and solves the quadratic problem for the increment with conjugate gradients in an inner loop: a truncated Gauss–Newton method (Gratton, Lawless and Nichols, SIAM J. Optim. 2007).
• Deep learning. Not plain Newton, but its approximations: the GGN and Fisher in Hessian-free optimization and K-FAC, Kronecker factors in Shampoo, diagonal scaling in Adam (lessons 06 and 12). Laplace approximations of the posterior also use the GGN.

7. Common mistakes #

• Inverting the Hessian. Solve the linear system with Cholesky. If Cholesky fails, the Hessian is not positive definite, and the Newton step may not be a descent direction.
• Calling Newton "globally quadratic". The quadratic rate is local. Far away you need damping or a trust region (worked example (b), from $x_0 = 0.1$).
• Pure Newton on a nonconvex function. It finds stationary points, including saddles and maxima. Check the curvature or use a modified Hessian.
• BFGS with a line search that does not enforce the curvature condition. Without $y^\top s > 0$, $H$ loses positive definiteness.
• L-BFGS on noisy minibatch gradients. Gradient differences then measure noise, not curvature. Use it with full batches (or large, fixed batches).
• Forming $J^\top J$ in a badly conditioned problem. It squares the condition number of $J$. Use QR or lstsq on $Jh \approx -r$.

8. Check yourself #

1. Gauss–Newton on the Rosenbrock residuals $r(x) = (1 - x_1,\ 10(x_2 - x_1^2))$ from $(-1.2, 1)$. Compute the first two steps. What happens to $F = \tfrac12\|r\|^2$?
   Answer

   $J = \begin{pmatrix} -1 & 0 \\ -20x_1 & 10 \end{pmatrix}$ is square and invertible, so $h = -J^{-1}r$. Step 1: $r = (2.2, -4.4)$ gives $x_1 = (1, -3.84)$. Step 2: $r = (0, -48.4)$ gives $x_2 = (1, 1)$, the solution. But $F$ rose from 12.1 to 1171.3 at step 1: LM or a line search would not take this step.
2. Why do Newton's method on $f$ and on $g(y) = f(Ty)$ give corresponding iterates, and what does this mean for feature scaling?
   Answer

   The Newton step in $y$ is $-(T^\top\nabla^2 f\,T)^{-1}T^\top\nabla f = T^{-1}\Delta x_{\rm nt}$, and the Armijo test uses only values and slopes, which are the same. So $y_k = T^{-1}x_k$. Rescaling the features does not change Newton's speed; it does change the speed of gradient descent.
3. A BFGS step on a nonconvex $f$ gives $y^\top s = -0.3$. What do implementations do, and what prevents the case?
   Answer

   They skip the update (keep $H_k$) or use a damped update. A line search that satisfies the Wolfe curvature condition guarantees $y^\top s > 0$, so the case does not happen.
4. In LM, how does $\lambda$ change after a step with $\varrho = 0.9$, and after one with $\varrho = -0.2$?
   Answer

   $\varrho = 0.9$: accept, and $\lambda \leftarrow \lambda\max\{1/3,\ 1 - 0.8^3\} = 0.488\,\lambda$ (the model was good). $\varrho = -0.2$: reject, $\lambda \leftarrow \nu\lambda$ (for example $2\lambda$) and $\nu \leftarrow 2\nu$ (the model was bad).
5. L-BFGS with $m = 10$ on $d = 10^6$ parameters: memory and cost per direction, compared with BFGS?
   Answer

   L-BFGS stores $2md = 2\cdot10^7$ numbers (160 MB in float64) and does about $4md = 4\cdot10^7$ multiplications per direction. BFGS would store $d^2 = 10^{12}$ numbers (8 TB).
6. Write the Hessian of $F = \tfrac12\|r\|^2$. Which term does Gauss–Newton drop, and when does it vanish?
   Answer

   $\nabla^2 F = J^\top J + \sum_i r_i \nabla^2 r_i$. Gauss–Newton drops $\sum_i r_i\nabla^2 r_i$. It vanishes at the solution when the residuals are zero there, and everywhere when $r$ is linear.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the Newton step at $x$?

   Show answer

2. What are the assumptions of local quadratic convergence of Newton's method?

   Show answer

3. Under quadratic convergence with constant about 1, the error is $10^{-3}$. What is it after the next step?

   Show answer

4. What does "Newton's method is affine invariant" mean?

   Show answer

5. State the secant condition.

   Show answer

6. Why does BFGS need $y^\top s > 0$, and what guarantees it?

   Show answer

7. Memory and time per step of L-BFGS with memory $m$ in dimension $d$?

   Show answer

8. What matrix does Gauss–Newton use, and what does it drop?

   Show answer

9. When does Gauss–Newton converge quadratically?

   Show answer

10. Write the Levenberg–Marquardt step and its two limits.

    Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Why not Newton for a large neural network?" — $d^2$ memory and $O(d^3)$ solves, an indefinite Hessian (saddles attract Newton), and noisy minibatch curvature. People use Hessian-vector products (Newton-CG), GGN or Fisher approximations (K-FAC, Shampoo) or Adam instead.
• "Newton vs BFGS vs L-BFGS?" — Exact Hessian, quadratic local rate, $O(d^3)$ per step; learned inverse Hessian, superlinear, $O(d^2)$; $m$ pairs, $O(md)$, the default for smooth full-batch problems.
• "Why does BFGS need a Wolfe line search?" — The curvature condition gives $y^\top s > 0$, which keeps $H$ positive definite, so every direction is a descent direction.
• "Gauss–Newton vs Newton for least squares?" — Gauss–Newton drops $\sum r_i\nabla^2 r_i$, needs only $J$, always gives a descent direction, and is quadratic for zero residuals; with large residuals it can be slow or fail, and Levenberg–Marquardt makes it safe.
• "What is 4D-Var, as an optimizer?" — Nonlinear least squares in the initial state, solved by incremental (truncated Gauss–Newton) iterations with adjoint gradients.

10. Practice #

The notebook for this lesson runs in Colab. It has four exercises, each with a check cell and a folded solution:

1. Newton's method with backtracking for $\ell_2$-regularized logistic regression on the breast cancer data set (standardized features); compare the number of iterations with gradient descent.
2. Implement the BFGS inverse update; check the secant condition and positive definiteness; run BFGS with a Wolfe line search on Rosenbrock and compare with scipy.optimize.minimize(method="BFGS").
3. Implement the L-BFGS two-loop recursion and check it against the dense BFGS matrix built from the same pairs.
4. Gauss–Newton and Levenberg–Marquardt for an exponential curve fit; compare with scipy.optimize.least_squares(method="lm").

11. Watch and read #

Visually Explained: Newton's Method in OptimizationVisually Explained · 11 min

The quadratic model, a numerical example, and "the good, the bad and the ugly" of Newton's method; compare with the 1D figure above.

Lecture 14: Newton's MethodRyan Tibshirani · CMU 10-725 Convex Optimization · 78 min

Damped Newton with backtracking, affine invariance, and the two-phase convergence analysis of section 3.

Harvard AM205 video 4.9 - Quasi-Newton methodsChris Rycroft · Harvard AM205 · 25 min

The secant condition and the BFGS update, as in section 3.

Harvard AM205 video 1.8 - Nonlinear least squaresChris Rycroft · Harvard AM205 · 27 min

Gauss–Newton and Levenberg–Marquardt for curve fitting, as in the exponential-fit figure.

• Boyd and Vandenberghe, Convex Optimization, free PDF: §9.5 (Newton step, Newton decrement, affine invariance, the two-phase analysis).
• Lieven Vandenberghe, ECE236C lectures (slides): "Newton's method" and "Quasi-Newton methods" (BFGS, positive definiteness, superlinear convergence, L-BFGS).
• Nocedal and Wright, Numerical Optimization, 2nd edition (Springer, 2006): ch. 3.4 (modified Newton), ch. 6 (BFGS), ch. 7.2 (L-BFGS, the two-loop recursion), ch. 10.3 (Gauss–Newton and Levenberg–Marquardt).
• Madsen, Nielsen and Tingleff, Methods for Non-Linear Least Squares Problems (2004), free: §3.1–3.2 (Gauss–Newton, Example 3.3, Levenberg–Marquardt, Algorithm 3.16).
• James Martens, Deep learning via Hessian-free optimization (ICML 2010), and New insights and perspectives on the natural gradient method (the GGN and the Fisher).
• SciPy, scipy.optimize.least_squares (methods lm and trf) and scipy.optimize.minimize.