USAAIO Lesson 21, from Week 7, fully worked. It builds the second-order Taylor model term by term, covers the Hessian and how its eigenvalues classify minima, maxima, and saddles, and derives Newton's update by setting the model's gradient to zero. It shows why Newton is exact on a quadratic, taking one step from [5,-3] to [0.2,0.4], and why it converges quadratically, with errors falling 5.9e-1, 3.4e-2, 1.3e-3, 1.7e-6, 3.3e-12. It then races IRLS Newton for logistic regression against gradient descent - 7 iterations against 1147, reaching the same w, matched to sklearn - and covers the O(d^3) per-step cost, the saddle trap that an indefinite Hessian creates, and L-BFGS and the natural gradient at scale. Every snippet runs standalone, and every number came from real execution. The lesson runs to 61 slides.
Subject: Machine Learning · 110 slides · code lesson
Open the interactive version of this deck · Homework for this lesson
Title
USAAIO · Lesson 21 · Week 7 (Optimization)
Gradient descent knows only which way is downhill. Today we use curvature: the Hessian, Newton's method derived from the second-order Taylor model, why it is exact on a quadratic and quadratically convergent, and — with real numbers — why 7 Newton steps beat 1147 gradient-descent steps on the same logistic fit.
Objectives
f and read off its gradient and Hessian H = ∇²fH's eigenvalues: all + → min, all − → max, mixed → saddleθ ← θ − H⁻¹∇f by minimizing the quadratic model, and prove it is exact in one step on a quadraticsklearn, and say why Newton costs O(d³) — with L-BFGS / natural gradient as the scalable fixesWarm-up
Discussion prompt
Before we open Lesson 21: Second-Order Methods: without looking back, what was the main idea of Hypothesis Testing, and what could you do by the end of it that you could not do before?
Hint: One sentence for the idea, one for the skill. If the second one is blank, that is the part to revisit.
Answer:
null vs alternative hypotheses, Type I/II errors, what a p-value really means, the t/chi-squared/F test family, statistical power, and the multiple-testing problem with Bonferroni. Build a one-sample t-test from scratch and verify against scipy.
Section
Part 1 of 9
Intuition
Standing on a hillside, the gradient tells you which way is downhill and how steep it is — but not how far the bottom is. Gradient descent papers over that gap with a hand-tuned learning rate.
Curvature — how fast the slope itself changes — is the missing information. On a gently curved valley you can stride far; on a sharply curved one you must step carefully.
Second-order methods read the curvature and compute the step size and direction at once, with no learning rate to tune. The object that stores curvature is the Hessian.
Counterexample
Discussion prompt
Standing on a hillside, the gradient tells you which way is downhill and how steep it is — but not how far the bottom is. Gradient descent papers over that gap with a hand-tuned learning rate.
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Answer:
Second-order methods read the curvature and compute the step size and direction at once, with no learning rate to tune. The object that stores curvature is the Hessian.
Concept
In one variable, near a point x₀ a smooth f looks like a parabola — its second-order Taylor polynomial. Slope f' sets the tilt; f'' sets how sharply it curves.
\[ f(x) \approx f(x_0) + f'(x_0)\,(x - x_0) + \tfrac{1}{2} f''(x_0)\,(x - x_0)^2 \]
If f''(x₀) > 0 the parabola opens up and has a bottom; that bottom is Newton's next guess. Everything today is the multi-variable version of this one line.
Analogy
Discussion prompt
Explain The 1-D Taylor picture by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
In one variable, near a point x₀ a smooth f looks like a parabola — its second-order Taylor polynomial. Slope f' sets the tilt; f'' sets how sharply it curves.
Concept
In many variables the tilt is the gradient ∇f (a vector) and the curvature is the Hessian H (a matrix). The Taylor model around θ₀, writing the step s = θ − θ₀, is:
\[ m(\theta) = f(\theta_0) + \nabla f(\theta_0)^\top s + \tfrac{1}{2}\, s^\top H\, s \]
The sᵀHs term is a quadratic form — exactly the object whose sign is governed by definiteness (Lesson 13). That is why the Hessian's eigenvalues decide the shape of the bowl.
Explain it
Discussion prompt
Explain The multivariable second-order model to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
In many variables the tilt is the gradient ∇f (a vector) and the curvature is the Hessian H (a matrix). The Taylor model around θ₀, writing the step s = θ − θ₀, is:
Section
Part 2 of 9 - the matrix of curvature
Concept
The Hessian collects every second partial derivative — entry (i, j) is how the i-th component of the gradient changes as you move in the j-th direction.
\[ H_{ij} = \frac{\partial^2 f}{\partial \theta_i\, \partial \theta_j}, \qquad H \in \mathbb{R}^{d \times d} \]
Hessian — The d×d matrix of all second partial derivatives of a scalar function f: R^d -> R. It is the Jacobian of the gradient field - how the gradient itself changes from point to point.
Concept
For any twice-continuously-differentiable f, mixed partials commute (Clairaut's / Schwarz's theorem): differentiating in order i then j gives the same result as j then i.
\[ \frac{\partial^2 f}{\partial \theta_i\, \partial \theta_j} = \frac{\partial^2 f}{\partial \theta_j\, \partial \theta_i} \;\Longrightarrow\; H = H^\top \]
So H is symmetric, which means everything from Lesson 13 applies: real eigenvalues, orthogonal eigenvectors, and definiteness read straight off the eigenvalue signs.
Estimation
Predict first
Our running example is the quadratic f(x) = ½ xᵀA x − bᵀx with A = [[3,1],[1,2]], b = [1,1]. Its gradient is ∇f = Ax − b and its Hessian is the constant matrix A. Let's confirm A's eigenvalues numerically:
Commit before you compute: what does Compute a Hessian: the running quadratic come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: Both eigenvalues positive → A is positive definite
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. eigvalsh returns 1.381966 and 3.618034, both > 0.
Worked example
Our running example is the quadratic f(x) = ½ xᵀA x − bᵀx with A = [[3,1],[1,2]], b = [1,1]. Its gradient is ∇f = Ax − b and its Hessian is the constant matrix A. Let's confirm A's eigenvalues numerically:
import numpy as np
A = np.array([[3., 1.], [1., 2.]]) # the Hessian of f = 1/2 xᵀAx − bᵀx
print(np.linalg.eigvalsh(A).round(6)) # [1.381966 3.618034]
print('symmetric?', np.allclose(A, A.T))Both eigenvalues positive → A is positive definite
Why: eigvalsh returns 1.381966 and 3.618034, both > 0. A symmetric matrix with all-positive eigenvalues is positive definite, so f is a convex, upward-opening bowl with a unique minimum.
| quantity | value (verified) |
|---|---|
| eig(A) smaller | 1.381966 |
| eig(A) larger | 3.618034 |
| both > 0 ? | yes → positive definite |
| symmetric ? | True |
Concept
At a critical point (∇f = 0) the linear term vanishes, so the model is f(θ₀) + ½ sᵀHs. The sign of sᵀHs in every direction s decides the shape:
| Hessian at ∇f = 0 | sᵀHs for all s ≠ 0 | critical point is a |
|---|---|---|
| positive definite (all λ > 0) | > 0 (curves up everywhere) | local minimum |
| negative definite (all λ < 0) | < 0 (curves down everywhere) | local maximum |
| indefinite (mixed-sign λ) | changes sign | saddle point |
| semidefinite (some λ = 0) | = 0 in some direction | inconclusive (flat direction) |
A zero gradient alone only says flat. The Hessian's eigenvalues are what certify which kind of flat.
Intuition
Picture a horse's saddle or a mountain pass: from front to back it dips down into a valley, but from side to side it rises over a ridge. There is a direction that goes down and a direction that goes up.
That is exactly a mixed-sign Hessian: a positive eigenvalue is an up-curving direction, a negative one is a down-curving direction. Sitting at the center you are at a critical point that is neither a min nor a max.
Figure (svg): A saddle surface: a curve dipping down through the center from left to right, and a curve rising up through the center from front to back, crossing at the saddle point.
Missing information
Discussion prompt
f = x² + y² gives H = diag(2, 2); f = x² − y² gives H = diag(2, −2). Take the eigenvalues with eigvalsh and classify. Standalone:
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
Every direction curves up, so the bowl x² + y² has a genuine bottom at the origin.
Worked example
f = x² + y² gives H = diag(2, 2); f = x² − y² gives H = diag(2, −2). Take the eigenvalues with eigvalsh and classify. Standalone:
import numpy as np
H_min = np.array([[2., 0.], [0., 2.]]) # f = x² + y²
H_saddle = np.array([[2., 0.], [0., -2.]]) # f = x² − y²
print(np.linalg.eigvalsh(H_min)) # [2. 2.] → all + → minimum
print(np.linalg.eigvalsh(H_saddle)) # [-2. 2.] → mixed → saddleH_min eigenvalues [2, 2]: all positive → local minimum
Why: Every direction curves up, so the bowl x² + y² has a genuine bottom at the origin.
H_saddle eigenvalues [−2, 2]: mixed sign → saddle point
Why: x² + curves up (eigenvalue +2), −y² curves down (eigenvalue −2). The origin is a saddle, not a minimum.
| function | H | eigenvalues | critical point |
|---|---|---|---|
| x² + y² | diag(2, 2) | [2, 2] | minimum |
| x² − y² | diag(2, −2) | [−2, 2] | saddle |
| −x² − y² | diag(−2, −2) | [−2, −2] | maximum |
Comparison
Comparison matrix
From Read three critical points off their Hessians: refill the critical point column from what you know. The rest of the table is as it appeared.
| function | H | eigenvalues | critical point |
|---|---|---|---|
| x² + y² | diag(2, 2) | [2, 2] | minimum |
| x² − y² | diag(2, −2) | [−2, 2] | saddle |
| −x² − y² | diag(−2, −2) | [−2, −2] | maximum |
Concept
For a random symmetric H in d dimensions, a critical point is a minimum only if all d eigenvalues happen to be positive — like flipping d coins all heads. As d grows that gets exponentially unlikely.
So in a neural network with millions of parameters, the overwhelming majority of critical points are saddles, not local minima. This is why 'zero gradient' is a weak stopping signal — and why the Hessian's signs matter so much on the exam.
Section
Part 3 of 9 - derive the update
Intuition
Newton's method: replace f by its second-order model, then jump straight to that model's minimum. Repeat from the new point.
Because the model is a parabola (quadratic), its minimum is a closed-form solve — no line search, no learning rate. The whole method is 'fit a parabola, hop to its bottom, refit.'
Ranking
Put in order
Put the moves of Minimize the quadratic model into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. Use ∇(gᵀs) = g and ∇(½ sᵀHs) = Hs for symmetric H (the matrix-calculus rules from Lesson 6).
Worked example
The local model in the step s is m(s) = f₀ + gᵀs + ½ sᵀHs, where g = ∇f(θ₀). Minimize it in s by setting its gradient to zero — step by step, nothing skipped.
Differentiate the model with respect to s
Why: Use ∇(gᵀs) = g and ∇(½ sᵀHs) = Hs for symmetric H (the matrix-calculus rules from Lesson 6). The constant f₀ drops out.
\[ \nabla_s\, m(s) = g + H s \]
Set the gradient to zero at the model's minimum
Why: The model is a convex parabola when H is positive definite, so its unique minimizer is where its gradient vanishes.
\[ g + H s = 0 \;\Longrightarrow\; H s = -g \]
Solve for the Newton step s
Why: Multiply by H⁻¹ (formally; in code we SOLVE, never invert). The step is the negative gradient rescaled by the inverse curvature.
\[ s = -H^{-1} g \]
Notation
Annotate
From Minimize the quadratic model — read this one piece at a time. What is each part doing?
On: \( \nabla_s\, m(s) = g + H s \)
Concept
The next iterate is θ₀ + s, giving the update every exam wants you to know cold:
\[ \boxed{\; \theta \leftarrow \theta - H^{-1}\, \nabla f \;} \]
Compare gradient descent, θ ← θ − η∇f: Newton replaces the scalar learning rate η with the matrix H⁻¹. It rescales the step per-direction by that direction's curvature — big steps where the bowl is flat, small steps where it is sharp.
Sorting
Sort into buckets
These are the pieces of Lesson 21: Second-Order Methods, out of order. Put each one back under the part of the lesson it belongs to.
Intuition
Think of H⁻¹ as a per-direction learning rate. Gradient descent uses one scalar η for every direction, so a valley that is steep one way and shallow another forces η to be tiny — bottlenecked by the steepest direction while the shallow one crawls.
Newton stretches the space so the bowl becomes round: in the stretched coordinates every direction has the same curvature, and the gradient points straight at the minimum. That is why it needs no tuning and takes so few steps.
The price is that you must know the curvature — you have to build and factor H. The rest of the lesson is about when that price is worth paying.
Anomaly
Predict first
A student writes this, and it looks reasonable:
Curvature scales the gradient, so multiply the gradient by the Hessian: θ ← θ − H ∇f.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: This multiplies BY the curvature, so it takes BIGGER steps in sharply-curved directions — exactly backwards, and it blows up.
Solve Hs = −g; the step carries the inverse Hessian, s = −H⁻¹g.
Why: This multiplies BY the curvature, so it takes BIGGER steps in sharply-curved directions — exactly backwards, and it blows up. The model minimization gave Hs = −g, i.e. s = −H⁻¹g, with the INVERSE.
Trap
Curvature scales the gradient, so multiply the gradient by the Hessian: θ ← θ − H ∇f.
θ ← θ − H ∇f ✗
Why: This multiplies BY the curvature, so it takes BIGGER steps in sharply-curved directions — exactly backwards, and it blows up. The model minimization gave Hs = −g, i.e. s = −H⁻¹g, with the INVERSE.
Solve Hs = −g; the step carries the inverse Hessian, s = −H⁻¹g.
θ ← θ − H⁻¹ ∇f ✓
Why: H⁻¹ takes SMALLER steps where curvature is large and larger steps where it is small — the correct rescaling. In code: step = np.linalg.solve(H, g); never form H⁻¹ explicitly.
Break the constraint
Discussion prompt
The rule this trap just fixed:
H⁻¹ takes SMALLER steps where curvature is large and larger steps where it is small — the correct rescaling. In code: step = np.linalg.solve(H, g); never form H⁻¹ explicitly.
Now break it on purpose. Build a case that violates it and follow the consequences until something visibly fails. Where does the failure first show up — and would you have noticed it if you had not been looking?
Hint: The dangerous rules are the ones whose violation still produces an answer. If yours fails loudly, try to find one that fails quietly.
Answer:
This multiplies BY the curvature, so it takes BIGGER steps in sharply-curved directions — exactly backwards, and it blows up. The model minimization gave Hs = −g, i.e. s = −H⁻¹g, with the INVERSE.
Section
Part 4 of 9 - one step, any start
Concept
For a quadratic f, the second-order Taylor model is not an approximation — it is f, exactly. So jumping to the model's minimum lands you at the true minimum, from any starting point, in a single step.
\[ f(x) = \tfrac{1}{2} x^\top A x - b^\top x \;\Longrightarrow\; \nabla f = Ax - b,\quad H = A \]
The minimum solves ∇f = 0, i.e. Ax = b, i.e. x★ = A⁻¹b. Newton's step x − A⁻¹(Ax − b) = A⁻¹b gives exactly that — watch it happen.
Pattern
Predict first
The table runs: start x | [5, −3] · gradient g = Ax − b | [11, −2] · after 1 Newton step | [0.2, 0.4]
In One Newton step from a wild start, given the rows so far: what is the next one — the row where quantity is true min A⁻¹b?
Correct: true min A⁻¹b | [0.2, 0.4]
| quantity | value (verified) |
|---|---|
| start x | [5, −3] |
| gradient g = Ax − b | [11, −2] |
| after 1 Newton step | [0.2, 0.4] |
| true min A⁻¹b | [0.2, 0.4] |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. g = Ax − b = [3·5 + 1·(−3), 1·5 + 2·(−3)] − [1, 1] = [12, −1] − [1, 1] = [11, −2].
Worked example
Take A = [[3,1],[1,2]], b = [1,1], and a deliberately silly start x = [5, −3]. One Newton step must land on A⁻¹b. Standalone and runnable:
import numpy as np
A = np.array([[3., 1.], [1., 2.]])
b = np.array([1., 1.])
x = np.array([5., -3.]) # arbitrary start
g = A @ x - b # gradient = [11., -2.]
x_new = x - np.linalg.solve(A, g) # Newton step (solve, don't invert)
print(x_new.round(6)) # [0.2 0.4]
print(np.linalg.solve(A, b).round(6)) # [0.2 0.4] = true min A⁻¹bGradient at the start is [11, −2]
Why: g = Ax − b = [3·5 + 1·(−3), 1·5 + 2·(−3)] − [1, 1] = [12, −1] − [1, 1] = [11, −2]. A big gradient — we are far away.
One step lands at [0.2, 0.4] = A⁻¹b, the exact minimum
Why: x_new = [5,−3] − A⁻¹[11,−2] = [0.2, 0.4], which equals A⁻¹b exactly. Distance from the start was irrelevant: the quadratic model is exact, so one step finishes the job.
| quantity | value (verified) |
|---|---|
| start x | [5, −3] |
| gradient g = Ax − b | [11, −2] |
| after 1 Newton step | [0.2, 0.4] |
| true min A⁻¹b | [0.2, 0.4] |
Discrimination
Sort into buckets
Sort these by value (verified), from memory, without looking back at One Newton step from a wild start. Telling them apart on the spot is the skill; the table is only where the answer happens to be written down.
Fill the middle
Fill in the blanks
From Confirm the step actually solves Hs = −g — one line has had its right-hand side removed. Put it back.
import numpy as np
A = np.array([[3., 1.], [1., 2.]])
b = np.array([1., 1.])
x = np.array([5., -3.])
g = A @ x - b # gradient at the start
s = -np.linalg.solve(A, g) # Newton step s = -H⁻¹g
print('H s + g =', (A @ s + g).round(12)) # [0. 0.]
print('grad at x+s =', (A @ (x + s) - b).round(12)) # [0. 0.]
Why: x is what everything below it consumes, so the wrong expression here fails later and somewhere else. The residual A·s + g is zero to machine precision, so s is the exact solution of Hs = −g.
Worked example
The whole derivation rested on Hs = −g. Let's confirm the computed step satisfies it exactly — the residual Hs + g should be numerically zero. Same A, b, start. Standalone:
import numpy as np
A = np.array([[3., 1.], [1., 2.]])
b = np.array([1., 1.])
x = np.array([5., -3.])
g = A @ x - b # gradient at the start
s = -np.linalg.solve(A, g) # Newton step s = -H⁻¹g
print('H s + g =', (A @ s + g).round(12)) # [0. 0.]
print('grad at x+s =', (A @ (x + s) - b).round(12)) # [0. 0.]H s + g = [0, 0] — the step exactly solves the model
Why: The residual A·s + g is zero to machine precision, so s is the exact solution of Hs = −g. This is the correctness check you can always run on a Newton step.
Gradient at x + s is [0, 0] — we are at the true minimum
Why: Because f is quadratic, landing at the model's minimum lands at f's minimum: ∇f(x+s) = A(x+s) − b = 0. One step, done.
| check | value (verified) |
|---|---|
| Newton step s | [−4.8, 3.4] |
| residual H s + g | [0, 0] |
| gradient at x + s | [0, 0] |
| x + s | [0.2, 0.4] (the minimum) |
Trade off
Comparison matrix
From Confirm the step actually solves Hs = −g: every row here is a choice with a cost. Fill the value (verified) column, then say which row you would actually pick and what you give up for it.
| check | value (verified) |
|---|---|
| Newton step s | [−4.8, 3.4] |
| residual H s + g | [0, 0] |
| gradient at x + s | [0, 0] |
| x + s | [0.2, 0.4] (the minimum) |
Section
Part 5 of 9 - digits double
Intuition
Suppose your error is 0.01 — two correct decimal places. Linear convergence with factor 0.1 takes it to 0.001, then 0.0001: one new digit per step. Steady, but slow.
Quadratic convergence squares it: 0.01 → 0.0001 → 0.00000001. Two digits, then four, then eight — the count of correct digits doubles every step. A handful of iterations reaches machine precision.
That is the whole appeal of Newton near a minimum: once you are close, you finish almost instantly. The catch is 'near' — far away, the quadratic model may not be trustworthy (Part 7).
Concept
Near a minimum, Newton's error squares each step: if eₖ = ‖θₖ − θ★‖, then eₖ₊₁ ≈ C · eₖ². Squaring a small number like 10⁻² gives 10⁻⁴ — so the number of correct digits roughly doubles per iteration.
\[ e_{k+1} \approx C\, e_k^2 \qquad\text{(quadratic)} \quad\text{vs.}\quad e_{k+1} \approx \rho\, e_k \;\;(\rho < 1) \qquad\text{(linear, GD)} \]
Gradient descent is linear: error shrinks by a constant factor ρ each step, adding a fixed number of digits. Newton accelerates — near the optimum it is untouchably fast.
Faded example
Fill in the blanks
See the digits double (1-D Newton), with the scaffolding fading: two lines are gone now — fill both.
import numpy as np
xstar = (3/4)(1/3) # true minimizer ≈ 0.90856
x = 1.5
for k in range(6):
err = abs(x - xstar)
print(k, round(x, 12), format(err, '.3e'))
x = x - (4x3 - 3) / (12x2) # Newton step on f'(x)
Why: Reproducing these unaided, rather than reading them, is what tells you the method has transferred. Each error is roughly the SQUARE of the previous (up to the constant C): 3.5e−2 squared ≈ 1.2e−3 ≈ the next error; 1.3e−3 squared ≈ 1.6e−6 ≈ the next.
Worked example
Minimize f(x) = x⁴ − 3x + 1. Setting f'(x) = 4x³ − 3 = 0 gives x★ = (3/4)^(1/3) ≈ 0.90856. Newton on f' is x ← x − f'/f'' with f'' = 12x². Watch the error column. Standalone:
import numpy as np
xstar = (3/4)**(1/3) # true minimizer ≈ 0.90856
x = 1.5
for k in range(6):
err = abs(x - xstar)
print(k, round(x, 12), format(err, '.3e'))
x = x - (4*x**3 - 3) / (12*x**2) # Newton step on f'(x)Errors: 5.9e−1 → 2.0e−1 → 3.5e−2 → 1.3e−3 → 1.7e−6 → 3.3e−12
Why: Each error is roughly the SQUARE of the previous (up to the constant C): 3.5e−2 squared ≈ 1.2e−3 ≈ the next error; 1.3e−3 squared ≈ 1.6e−6 ≈ the next. Correct digits double every step.
| iter | x | error |x − x★| | ≈ prev² |
|---|---|---|---|
| 0 | 1.500000000000 | 5.914e−01 | — |
| 1 | 1.111111111111 | 2.026e−01 | 0.35 |
| 2 | 0.943240740741 | 3.468e−02 | 0.041 |
| 3 | 0.909819776353 | 1.259e−03 | 1.2e−3 |
| 4 | 0.908562039132 | 1.743e−06 | 1.6e−6 |
| 5 | 0.908560296419 | 3.343e−12 | 3.0e−12 |
Pattern
Step through it
Step through See the digits double (1-D Newton) one row at a time. What is driving the change, and what would the row after the last one be?
Section
Part 6 of 9 - IRLS vs GD
Concept
For logistic regression with predictions p = σ(Xw) and labels y, the average cross-entropy loss has a famously clean gradient and Hessian. Let W = diag(pᵢ(1 − pᵢ)):
\[ \nabla L = \tfrac{1}{n} X^\top (p - y), \qquad H = \tfrac{1}{n} X^\top W X \]
Because each pᵢ(1 − pᵢ) ≥ 0, the Hessian XᵀWX is positive semidefinite — the loss is convex, so there are no saddles or local maxima to trip over. Newton here is called IRLS (iteratively reweighted least squares).
Ranking
Put in order
Put the moves of Where that gradient comes from (chain rule) into the order they have to happen.
Why: These are the moves of the worked example in the order it makes them, and each one is set up by the one before it. Differentiate −[y log p + (1−y) log(1−p)] in p: the two terms give −y/p and +(1−y)/(1−p).
Worked example
The clean form is worth deriving once. The per-sample loss is the negative log-likelihood ℓ = −[y log p + (1−y) log(1−p)] with p = σ(z) and z = xᵀw. Chain it, one link at a time.
Derivative of the loss w.r.t. p
Why: Differentiate −[y log p + (1−y) log(1−p)] in p: the two terms give −y/p and +(1−y)/(1−p).
\[ \frac{\partial \ell}{\partial p} = -\frac{y}{p} + \frac{1-y}{1-p} \]
Use the sigmoid identity σ'(z) = p(1 − p)
Why: The sigmoid's derivative is p(1−p). Multiply by ∂ℓ/∂p and the fractions cancel beautifully.
\[ \frac{\partial \ell}{\partial z} = \left(-\frac{y}{p} + \frac{1-y}{1-p}\right) p(1-p) = p - y \]
Chain to w via z = xᵀw, then average over n samples
Why: ∂z/∂w = x, so ∂ℓ/∂w = (p − y)x. Stacking all samples and averaging gives the matrix form — exactly what the code computes.
\[ \nabla L = \tfrac{1}{n}\sum_i (p_i - y_i)\, x_i = \tfrac{1}{n} X^\top (p - y) \]
| link | result |
|---|---|
| ∂ℓ/∂p | −y/p + (1−y)/(1−p) |
| σ'(z) | p(1 − p) |
| ∂ℓ/∂z | p − y (fractions cancel) |
| ∇L (all n) | Xᵀ(p − y)/n |
Notation
Annotate
From Where that gradient comes from (chain rule) — read this one piece at a time. What is each part doing?
On: \( \nabla L = \tfrac{1}{n}\sum_i (p_i - y_i)\, x_i = \tfrac{1}{n} X^\top (p - y) \)
Concept
We generate 200 points with two features (plus a bias column, so X is 200×3) from a known w_true = [0.5, −1.5, 2.0], seed 0 so every run matches. Standalone:
import numpy as np
np.random.seed(0)
n = 200
X = np.random.randn(n, 2)
Xb = np.c_[np.ones(n), X] # bias + 2 features → 200×3
w_true = np.array([0.5, -1.5, 2.0])
p = 1 / (1 + np.exp(-(Xb @ w_true)))
y = (np.random.rand(n) < p).astype(float) # Bernoulli labels
print(Xb.shape, round(y.mean(), 2)) # (200, 3) 0.57
print(Xb[:3].round(4))Xb is (200, 3); the labels are 57% ones
Why: np.c_ prepends the ones column, giving 3 columns (bias + 2 features). With seed 0 the label mean is exactly 0.57 — reproducible for every student.
| object | value (verified, seed 0) |
|---|---|
| Xb.shape | (200, 3) |
| y.mean() | 0.57 |
| Xb[0] | [1.0, 1.7641, 0.4002] |
| Xb[1] | [1.0, 0.9787, 2.2409] |
Worked example
Run Newton on the logistic loss until the step is tiny. We add 1e−8·I to the Hessian for numerical safety (it never changes the answer here). Standalone and runnable:
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
w = np.zeros(3)
for it in range(100):
p = sig(Xb @ w)
g = Xb.T @ (p - y) / n # gradient
H = (Xb.T * (p*(1-p))) @ Xb / n + 1e-8*np.eye(3) # Hessian XᵀWX
step = np.linalg.solve(H, g) # SOLVE, not inv
w = w - step
if np.linalg.norm(step) < 1e-8:
break
print('Newton iters:', it+1) # 7
print('w =', w.round(6)) # [0.48834 -1.34714 2.672486]The loss drops to its floor in 7 iterations
Why: The per-iteration trace: loss 0.6931 → 0.4245 → 0.3756 → 0.3660 → 0.36540 → 0.365401 → 0.365401. By iteration 5 the step norm is already 7.6e−5; by iteration 7 it is below 1e−8 and the loop stops.
| Newton iter | loss | ‖step‖ | ‖grad‖ |
|---|---|---|---|
| 0 | 0.69314718 | 1.356e+00 | 3.438e−01 |
| 1 | 0.42447689 | 8.768e−01 | 9.140e−02 |
| 2 | 0.37560792 | 5.943e−01 | 2.842e−02 |
| 3 | 0.36598198 | 1.920e−01 | 5.765e−03 |
| 4 | 0.36540413 | 1.461e−02 | 3.818e−04 |
| 5 | 0.36540133 | 7.580e−05 | 1.968e−06 |
| 6 | 0.36540133 | 2.058e−09 | 5.371e−11 |
Pattern
Step through it
Step through Newton (IRLS) converges in 7 steps one row at a time. What is driving the change, and what would the row after the last one be?
Fill the middle
Fill in the blanks
From Gradient descent needs 1147 steps — one line has had its right-hand side removed. Put it back.
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
wg = np.zeros(3)
for it in range(1000000):
p = sig(Xb @ wg)
g = Xb.T @ (p - y) / n
wg = wg - 0.5 * g # fixed learning rate
if np.linalg.norm(g) < 1e-8:
break
print('GD iters:', it+1) # 1147
print('w =', wg.round(6)) # same as Newton
Why: wg is what everything below it consumes, so the wrong expression here fails later and somewhere else. The loss is essentially flat (0.3654013) by iteration ~400, yet the gradient norm only crosses 1e−8 at iteration 1147.
Worked example
Same loss, same data, same stopping tolerance ‖g‖ < 1e−8. Plain gradient descent with a fixed step 0.5. Standalone:
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
wg = np.zeros(3)
for it in range(1000000):
p = sig(Xb @ wg)
g = Xb.T @ (p - y) / n
wg = wg - 0.5 * g # fixed learning rate
if np.linalg.norm(g) < 1e-8:
break
print('GD iters:', it+1) # 1147
print('w =', wg.round(6)) # same as NewtonGD reaches the same w — but takes 1147 iterations
Why: The loss is essentially flat (0.3654013) by iteration ~400, yet the gradient norm only crosses 1e−8 at iteration 1147. GD's linear convergence adds a fixed few digits per step, so squeezing the last digits is slow.
| GD iter | loss | ‖grad‖ |
|---|---|---|
| 0 | 0.69314718 | 3.438e−01 |
| 10 | 0.45468248 | 1.245e−01 |
| 50 | 0.37549626 | 2.816e−02 |
| 100 | 0.36720375 | 1.054e−02 |
| 400 | 0.36540186 | 1.651e−04 |
| 1147 | 0.36540133 | < 1e−08 (stop) |
Scale up
Step through it
Step through Gradient descent needs 1147 steps and watch the numbers move. Now imagine the input ten times bigger: which column is the one that stops this being practical?
Estimation
Predict first
Both methods must land on the same weights, and those weights must match sklearn's logistic regression (with regularization effectively off, C = 1e12). Standalone:
Commit before you compute: what does Same answer, checked against sklearn come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: Our Newton weights equal sklearn's to 6 places
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. Both give w = [0.48834, −1.34714, 2.672486].
Worked example
Both methods must land on the same weights, and those weights must match sklearn's logistic regression (with regularization effectively off, C = 1e12). Standalone:
import numpy as np
from sklearn.linear_model import LogisticRegression
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
w = np.zeros(3)
for it in range(100):
p = sig(Xb@w); g = Xb.T@(p-y)/n
H = (Xb.T*(p*(1-p)))@Xb/n + 1e-8*np.eye(3)
w = w - np.linalg.solve(H, g)
m = LogisticRegression(C=1e12, fit_intercept=False, tol=1e-10,
max_iter=1000).fit(Xb, y)
print('our Newton:', w.round(6)) # [0.48834 -1.34714 2.672486]
print('sklearn :', m.coef_.ravel().round(6))
print('match?', np.allclose(w, m.coef_.ravel(), atol=1e-4)) # TrueOur Newton weights equal sklearn's to 6 places
Why: Both give w = [0.48834, −1.34714, 2.672486]. np.allclose returns True. Note w ≠ w_true = [0.5, −1.5, 2.0]: with only 200 noisy points the maximum-likelihood fit differs slightly from the generating weights — that is sampling noise, not a bug.
| source | w₀ | w₁ | w₂ |
|---|---|---|---|
| our Newton (7 iters) | 0.48834 | −1.34714 | 2.672486 |
| gradient descent (1147 iters) | 0.48834 | −1.34714 | 2.672486 |
| sklearn LogisticRegression | 0.48834 | −1.34714 | 2.672486 |
| w_true (generating) | 0.5 | −1.5 | 2.0 |
Discrimination
Sort into buckets
Sort these by w₀, from memory, without looking back at Same answer, checked against sklearn. Telling them apart on the spot is the skill; the table is only where the answer happens to be written down.
Worked example
Newton had no knob to tune. Gradient descent does — and it matters enormously. Sweep the learning rate and count iterations to the same ‖g‖ < 1e−8. Standalone:
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
def gd_iters(lr):
w = np.zeros(3)
for it in range(100000):
g = Xb.T @ (sig(Xb@w) - y) / n
w = w - lr*g
if np.linalg.norm(g) < 1e-8: return it+1
for lr in [0.1, 0.5, 1.0, 2.0]:
print(lr, gd_iters(lr))Same problem, 20× spread in iterations just from the rate
Why: lr=0.1 needs 5766 iterations; lr=2.0 needs only 281. A poorly chosen rate is 20× slower — and a too-large rate on a worse-conditioned problem would diverge outright. Newton sidesteps the whole issue by reading the curvature.
| learning rate | GD iterations | vs Newton (7) |
|---|---|---|
| 0.1 | 5766 | 824× |
| 0.5 | 1147 | 164× |
| 1.0 | 570 | 81× |
| 2.0 | 281 | 40× |
Pattern
Step through it
Step through GD's iteration count depends wildly on the learning rate one row at a time. What is driving the change, and what would the row after the last one be?
Section
Part 7 of 9 - why not always Newton
Concept
Each Newton step must solve the d×d system Hs = −g. A dense factorization (LU / Cholesky) is O(d³) arithmetic, and merely storing the Hessian is O(d²) memory.
\[ \text{solve } H s = -g:\quad O(d^3)\ \text{time}, \quad O(d^2)\ \text{memory} \]
Doubling d multiplies the per-step work by about 2³ = 8. For a model with d in the millions, d² alone is a matrix too big to store, let alone factor — Newton is simply infeasible at that scale.
Explain it
Discussion prompt
Explain Newton costs O(d³) per step to a student a year behind you. No notation, no jargon they have not met — and it still has to be true.
Hint: If your explanation needs a symbol they have never seen, you are describing the notation rather than the idea.
Answer:
Each Newton step must solve the d×d system Hs = −g. A dense factorization (LU / Cholesky) is O(d³) arithmetic, and merely storing the Hessian is O(d²) memory.
Concept
The honest comparison is iterations × cost-per-iteration. GD does many cheap steps; Newton does few expensive ones. Which wins depends entirely on d.
| method | cost / iter | iters (our logistic, d=3) | shines when |
|---|---|---|---|
| gradient descent | O(n·d) | 1147 | d huge (millions) |
| Newton | O(n·d² + d³) | 7 | d small–medium (≤ 1000s) |
| L-BFGS | O(m·d) | tens | d large, want curvature |
At d = 3 the d³ term is nothing and Newton's 164× fewer iterations dominate. At d = 10⁶ the d³ solve is a wall, and cheap first-order (or L-BFGS) wins despite needing far more steps. There is no universally best method — only the best for your d.
Comparison
Comparison matrix
From The real tradeoff: per-step vs total: refill the cost / iter column from what you know. The rest of the table is as it appeared.
| method | cost / iter | iters (our logistic, d=3) | shines when |
|---|---|---|---|
| gradient descent | O(n·d) | 1147 | d huge (millions) |
| Newton | O(n·d² + d³) | 7 | d small–medium (≤ 1000s) |
| L-BFGS | O(m·d) | tens | d large, want curvature |
Worked example
Time np.linalg.solve on random SPD systems as d doubles. The per-step time scales like d³ — each doubling multiplies it by roughly 8. Standalone (absolute microseconds vary by machine; the ratio is the point):
import numpy as np, time
np.random.seed(0)
for d in [50, 100, 200, 400]:
M = np.random.randn(d, d); M = M @ M.T + d*np.eye(d) # SPD
v = np.random.randn(d)
t0 = time.perf_counter()
for _ in range(20):
np.linalg.solve(M, v)
print(f'd={d:>4} {(time.perf_counter()-t0)/20*1e6:8.1f} us')Each doubling of d multiplies solve time by ~8
Why: A representative run: d=50 → ~40 µs, d=100 → ~96 µs, d=200 → ~525 µs, d=400 → ~4680 µs. Ratios 96/40 ≈ 2.4, 525/96 ≈ 5.5, 4680/525 ≈ 8.9 — trending toward the 8× that O(d³) predicts as d grows. Your exact microseconds will differ; the growth shape won't.
| d | solve time (µs, representative) | ratio to prev |
|---|---|---|
| 50 | ≈ 40 | — |
| 100 | ≈ 96 | ≈ 2.4× |
| 200 | ≈ 525 | ≈ 5.5× |
| 400 | ≈ 4680 | ≈ 8.9× |
Pattern
Step through it
Step through Watch the solve time grow ~8× per doubling one row at a time. What is driving the change, and what would the row after the last one be?
Intuition
Newton jumps to the model's critical point. If the current Hessian is indefinite (you are near a saddle), that critical point is the saddle itself — and Newton walks toward it instead of downhill.
Worse, in a negative-curvature direction H⁻¹ flips the sign of the step, sending you uphill. Raw Newton is only trustworthy where H is positive definite — near a minimum. Elsewhere it needs safeguards.
Faded example
Fill in the blanks
Raw Newton jumps straight to a saddle, with the scaffolding fading: two lines are gone now — fill both.
import numpy as np
H = np.array([[2., 0.], [0., -2.]]) # indefinite: f = x² − y²
pt = np.array([1., 1.])
g = np.array([2pt[0], -2pt[1]]) # gradient = [2, -2]
newton = pt - np.linalg.solve(H, g)
print('Newton lands at:', newton) # [0. 0.] ← the SADDLE
print('eig(H):', np.linalg.eigvalsh(H)) # [-2. 2.] indefinite
Why: Reproducing these unaided, rather than reading them, is what tells you the method has transferred. solve(H, g) = [1, 1], so newton = (1,1) − (1,1) = (0,0).
Worked example
On f(x, y) = x² − y² at the point (1, 1): gradient [2, −2], Hessian diag(2, −2) (indefinite). See where one Newton step lands. Standalone:
import numpy as np
H = np.array([[2., 0.], [0., -2.]]) # indefinite: f = x² − y²
pt = np.array([1., 1.])
g = np.array([2*pt[0], -2*pt[1]]) # gradient = [2, -2]
newton = pt - np.linalg.solve(H, g)
print('Newton lands at:', newton) # [0. 0.] ← the SADDLE
print('eig(H):', np.linalg.eigvalsh(H)) # [-2. 2.] indefiniteOne Newton step goes from (1,1) to (0,0) — the saddle
Why: solve(H, g) = [1, 1], so newton = (1,1) − (1,1) = (0,0). The origin is the SADDLE of x² − y², not a minimum (which doesn't exist — the surface runs to −∞ along y). Newton was pulled to the model's critical point regardless of type.
| quantity | value (verified) |
|---|---|
| point | [1, 1] |
| gradient | [2, −2] |
| eig(H) | [−2, 2] (indefinite) |
| Newton lands at | [0, 0] = the saddle |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
One Newton step goes from (1,1) to (0,0) — the saddle
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
On f(x, y) = x² − y² at the point (1, 1): gradient [2, −2], Hessian diag(2, −2) (indefinite). See where one Newton step lands. Standalone:
Anomaly
Predict first
A student writes this, and it looks reasonable:
The formula is θ ← θ − H⁻¹∇f, so just solve Hs = −g and step, whatever H looks like right now.
It is wrong. Say what breaks — and say it before you turn the page.
Correct: In a negative-curvature direction, H⁻¹ FLIPS the sign of the gradient, so the 'step' points UPHILL — Newton climbs toward a saddle or maximum instead of descending.
Check H first: only step raw when it is positive definite; otherwise modify it (add τI until all eigenvalues are > 0) or fall back to a gradient step.
Why: In a negative-curvature direction, H⁻¹ FLIPS the sign of the gradient, so the 'step' points UPHILL — Newton climbs toward a saddle or maximum instead of descending. On x²−y² it walked from (1,1) straight to the saddle (0,0).
Trap
The formula is θ ← θ − H⁻¹∇f, so just solve Hs = −g and step, whatever H looks like right now.
Solve and step with an indefinite H
Why: In a negative-curvature direction, H⁻¹ FLIPS the sign of the gradient, so the 'step' points UPHILL — Newton climbs toward a saddle or maximum instead of descending. On x²−y² it walked from (1,1) straight to the saddle (0,0).
Check H first: only step raw when it is positive definite; otherwise modify it (add τI until all eigenvalues are > 0) or fall back to a gradient step.
Ensure H ≻ 0 (or add τI), then solve and damp with t ≤ 1
Why: A positive-definite H guarantees the step is a descent direction; the τI shift lifts any zero/negative eigenvalue off the floor. A line-search factor t ≤ 1 keeps you inside the region where the quadratic model is trustworthy. This is what production Newton solvers actually do.
Two truths and a lie
Sort into buckets
Some of these hold up and some are the exact mistakes this lesson is built to prevent. Sort them.
x₀ a smooth f looks like a parabola — its second-order Taylor polynomial. Slope f' sets the tilt; f'' sets how sharply it curves.; The Hessian collects every second partial derivative — entry (i, j) is how the i-th component of the gradient changes as you move in the j-th direction.; At a critical point (∇f = 0) the linear term vanishes, so the model is f(θ₀) + ½ sᵀHs. The sign of sᵀHs in every direction s decides the shape:θ ← θ − H ∇f.; The formula is θ ← θ − H⁻¹∇f, so just solve Hs = −g and step, whatever H looks like right now.Concept
Two standard fixes keep Newton honest away from a minimum. Damping (line search): take θ ← θ − t·H⁻¹∇f with a step length t ≤ 1 chosen so the loss actually decreases. Modified Hessian: add τI to push all eigenvalues positive before solving — exactly what our +1e−8·I did in miniature.
Trust regions go further: only trust the quadratic model within a radius, and shrink that radius when the model over-promises. These turn raw Newton into a globally reliable method (Levenberg–Marquardt is the famous instance).
Section
Part 8 of 9 - L-BFGS & natural gradient
Concept
L-BFGS never forms or stores the d×d Hessian. It keeps only the last m (say 10) pairs of gradient- and step-differences and uses them to apply an approximate H⁻¹ to a vector — second-order-like curvature at near first-order cost and O(m·d) memory.
It is the default for medium-scale smooth problems (and scipy.optimize.minimize(method='L-BFGS-B')). You get much of Newton's speed without paying O(d³) or O(d²).
Concept
Natural gradient replaces the Hessian with the Fisher information matrix F — the expected Hessian of the log-likelihood — giving the update θ ← θ − η F⁻¹∇f. F is always positive semidefinite, so it never produces the saddle-seeking behavior raw Newton can.
For logistic regression the Fisher matrix equals the Hessian XᵀWX — natural gradient and Newton coincide there. Scalable approximations of F (like K-FAC) power second-order training of large neural nets.
Intuition
A quick decision guide. Small d and a Hessian you can write down (logistic regression, GLMs, small nets): use full Newton / IRLS — it converges in a handful of steps.
Medium d, smooth loss, but d² too big to store: use L-BFGS — near-Newton speed with O(m·d) memory. Huge d (deep nets): use first-order (Adam, SGD) or approximate second-order (K-FAC, natural gradient) — you can never afford the exact Hessian.
The through-line: second-order information is always helpful and rarely free. You buy as much of it as your d can afford.
Analogy
Discussion prompt
Explain Which method, when by analogy to something with no Machine Learning in it at all — a queue, a recipe, a map, a bank balance, whatever fits. Then say where your analogy breaks.
Hint: An analogy that never breaks is not an analogy, it is the same idea wearing a hat. Find the seam — that is the part that is actually new.
Answer:
A quick decision guide. Small d and a Hessian you can write down (logistic regression, GLMs, small nets): use full Newton / IRLS — it converges in a handful of steps.
Constraint
Discussion prompt
Run The second-order toolkit with this step confiscated:
Step: solve Hs = −g (never inv) → Newton update θ ← θ − H⁻¹∇f
Is it still possible? If it is, say what takes its place and what it costs you. If it is not, say exactly what that step was providing that nothing else does.
Hint: A step you can drop for free was never load-bearing. If you cannot drop it, name the thing that goes wrong the moment it is gone.
Answer:
f by f₀ + gᵀs + ½ sᵀHs; the Hessian H = ∇²f is symmetric∇f = 0, eigenvalues of H — all + → min, all − → max, mixed → saddleHs = −g (never inv) → Newton update θ ← θ − H⁻¹∇fO(d³)/step and unsafe far from a min → damp / trust-region for safety; L-BFGS or natural gradient at scalePattern
f by f₀ + gᵀs + ½ sᵀHs; the Hessian H = ∇²f is symmetric∇f = 0, eigenvalues of H — all + → min, all − → max, mixed → saddleHs = −g (never inv) → Newton update θ ← θ − H⁻¹∇fO(d³)/step and unsafe far from a min → damp / trust-region for safety; L-BFGS or natural gradient at scaleEdge cases
Discussion prompt
The second-order toolkit works on the cases you have just seen. Push it to the edge: what is the most degenerate input it still handles — empty, zero, one item, everything equal — and what is the first case where it stops being true? Name the case, not just "it breaks".
Hint: Try the smallest legal input, then the largest, then the one where two things collide. Methods are specified at their edges; the middle takes care of itself.
Answer:
f by f₀ + gᵀs + ½ sᵀHs; the Hessian H = ∇²f is symmetric∇f = 0, eigenvalues of H — all + → min, all − → max, mixed → saddleHs = −g (never inv) → Newton update θ ← θ − H⁻¹∇fO(d³)/step and unsafe far from a min → damp / trust-region for safety; L-BFGS or natural gradient at scaleElimination
Eliminate the wrong options
Newton's method updates the parameters by:
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: Minimizing the second-order model gives Hs = −g, so s = −H⁻¹g and θ ← θ − H⁻¹∇f. The inverse Hessian rescales the gradient by the local curvature — small steps where curvature is large, big where it is small.
Check
Read it off the model minimization Hs = −g.
Check your understanding
Newton's method updates the parameters by:
Answer: A
Why: Minimizing the second-order model gives Hs = −g, so s = −H⁻¹g and θ ← θ − H⁻¹∇f. The inverse Hessian rescales the gradient by the local curvature — small steps where curvature is large, big where it is small.
Prediction
Predict first
At a point where ∇f = 0, the Hessian has eigenvalues {+4, −1}. The point is a:
Answer it in your own words, now, with nothing to choose from. The options are on the next slide — and picking the right one off a list is an easier skill than producing it.
Correct: saddle point (indefinite Hessian)
Why: Mixed-sign eigenvalues make the Hessian indefinite: the surface curves up along the +4 eigenvector and down along the −1 eigenvector. That up-and-down is exactly a saddle point.
Check
The gradient is zero; look at the eigenvalue signs.
Check your understanding
At a point where ∇f = 0, the Hessian has eigenvalues {+4, −1}. The point is a:
Answer: A
Why: Mixed-sign eigenvalues make the Hessian indefinite: the surface curves up along the +4 eigenvector and down along the −1 eigenvector. That up-and-down is exactly a saddle point.
Prediction
Predict first
Compared with gradient descent on a well-behaved (locally convex) problem, Newton's method typically…
Answer it in your own words, now, with nothing to choose from. The options are on the next slide — and picking the right one off a list is an easier skill than producing it.
Correct: needs far fewer iterations but costs O(d³) per step
Why: Newton converges quadratically — here 7 iterations vs GD's 1147 — but each step solves a d×d system at O(d³) cost. Fast in iterations, expensive per iteration, so it shines on small-to-medium d and is replaced by L-BFGS / natural gradient at scale.
Check
Trade off iteration count against per-step price.
Check your understanding
Compared with gradient descent on a well-behaved (locally convex) problem, Newton's method typically…
Answer: A
Why: Newton converges quadratically — here 7 iterations vs GD's 1147 — but each step solves a d×d system at O(d³) cost. Fast in iterations, expensive per iteration, so it shines on small-to-medium d and is replaced by L-BFGS / natural gradient at scale.
Elimination
Eliminate the wrong options
Why compute the Newton step as np.linalg.solve(H, g) rather than np.linalg.inv(H) @ g?
3 of these 4 are wrong. Strike them one at a time, and say what rules each one out before you strike the next. The survivor is the answer.
Survives elimination: A
Why: solve factors H once (LU/Cholesky) and back-substitutes — fewer operations and better backward stability than building the full inverse and multiplying. Same math, safer and faster arithmetic.
Check
The coding section rewards numerically sound habits.
Check your understanding
Why compute the Newton step as np.linalg.solve(H, g) rather than np.linalg.inv(H) @ g?
Answer: A
Why: solve factors H once (LU/Cholesky) and back-substitutes — fewer operations and better backward stability than building the full inverse and multiplying. Same math, safer and faster arithmetic.
Section
Part 9 of 9 - the project
Concept
Implement the logistic gradient and Hessian, run Newton, race it against gradient descent, and confirm both match sklearn. You have derived every piece — now assemble it yourself on the seed-0 dataset.
| # | requirement | tool |
|---|---|---|
| 1 | gradient Xᵀ(p−y)/n and Hessian XᵀWX/n | sigmoid, W = p(1−p) |
| 2 | Newton loop via solve(H, g) until ‖step‖ < 1e−8 | np.linalg.solve |
| 3 | race gradient descent to the same tolerance | ‖grad‖-based stop |
Build rules: type every line yourself, solve the system (never inv), add a tiny 1e−8·I to H for safety, and when a shape error appears, read the shapes — don't delete the line.
Counterexample
Discussion prompt
Build rules: type every line yourself, solve the system (never inv), add a tiny 1e−8·I to H for safety, and when a shape error appears, read the shapes — don't delete the line.
That is stated as though it always holds. Do one of two things: produce a case where it fails, or say precisely what rules such a case out. "It just does" is not on the menu.
Hint: Hunt at the extremes first — zero, one, negative, empty, equal. If every extreme survives, the reason they survive is the proof.
Fill the middle
Fill in the blanks
From Milestone 1 - gradient & Hessian — one line has had its right-hand side removed. Put it back.
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
def grad_hess(w):
p = sig(Xb @ w)
g = Xb.T @ (p - y) / n
H = (Xb.T * (p*(1-p))) @ Xb / n
return g, H
g, H = grad_hess(np.zeros(3))
print(g.shape, H.shape) # (3,) (3, 3)
print(np.linalg.eigvalsh(H).round(4)) # all positive → PD
Why: p is what everything below it consumes, so the wrong expression here fails later and somewhere else. eigvalsh(H) = [0.2284, 0.2392, 0.271] — every value > 0, so H is positive definite and the loss is locally a bowl.
Worked example
Your turn: write a function returning the logistic gradient and Hessian at w. Predict the Hessian's shape and its definiteness for 3 features before you print.
Hint: p = sig(Xb@w); g = Xb.T@(p−y)/n; H = (Xb.T*(p*(1−p)))@Xb/n. H is 3×3 and positive semidefinite.
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
def grad_hess(w):
p = sig(Xb @ w)
g = Xb.T @ (p - y) / n
H = (Xb.T * (p*(1-p))) @ Xb / n
return g, H
g, H = grad_hess(np.zeros(3))
print(g.shape, H.shape) # (3,) (3, 3)
print(np.linalg.eigvalsh(H).round(4)) # all positive → PDAt w = 0 the Hessian eigenvalues are all positive
Why: eigvalsh(H) = [0.2284, 0.2392, 0.271] — every value > 0, so H is positive definite and the loss is locally a bowl. Shapes: g is (3,), H is (3, 3).
| object | value (verified) |
|---|---|
| g.shape | (3,) |
| H.shape | (3, 3) |
| eig(H) at w=0 | [0.2284, 0.2392, 0.271] |
| definiteness | positive definite (all > 0) |
Missing information
Discussion prompt
Your turn: loop the Newton step until the step is tiny. Predict roughly how many iterations before you run it (hint: quadratic convergence is fast).
What do you need to know — or decide — before the first line can be written? List everything the problem has to hand you.
Hint: Anything you would have to invent to get started is a thing the problem must supply.
Answer:
The step norm falls 1.36 → 0.88 → 0.59 → 0.19 → 0.015 → 7.6e−5 → 2.1e−9 — roughly squaring each step, the fingerprint of quadratic convergence. Below 1e−8 at iteration 7, the loop breaks.
Worked example
Your turn: loop the Newton step until the step is tiny. Predict roughly how many iterations before you run it (hint: quadratic convergence is fast).
Hint: step = np.linalg.solve(H + 1e−8*np.eye(3), g); w = w − step; stop when norm(step) < 1e−8.
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
w = np.zeros(3)
for it in range(100):
p = sig(Xb @ w)
g = Xb.T @ (p - y) / n
H = (Xb.T * (p*(1-p))) @ Xb / n + 1e-8*np.eye(3)
step = np.linalg.solve(H, g)
w = w - step
if np.linalg.norm(step) < 1e-8:
break
print('iters:', it+1) # 7
print('w =', w.round(6))Converges in 7 iterations to w = [0.48834, −1.34714, 2.672486]
Why: The step norm falls 1.36 → 0.88 → 0.59 → 0.19 → 0.015 → 7.6e−5 → 2.1e−9 — roughly squaring each step, the fingerprint of quadratic convergence. Below 1e−8 at iteration 7, the loop breaks.
| metric | value (verified) |
|---|---|
| Newton iterations | 7 |
| stop rule | ‖step‖ < 1e−8 |
| final w | [0.48834, −1.34714, 2.672486] |
Reverse engineer
Discussion prompt
Work backwards. The example finished here:
Converges in 7 iterations to w = [0.48834, −1.34714, 2.672486]
What was it asked to do, and what must it have been given? Reconstruct the problem from its answer.
Hint: Every quantity in the result had to enter somewhere. Account for each one.
Answer:
Your turn: loop the Newton step until the step is tiny. Predict roughly how many iterations before you run it (hint: quadratic convergence is fast).
Estimation
Predict first
Your turn: run GD on the same loss to the same tolerance and count iterations. Predict the ratio to Newton before you print (order of magnitude).
Commit before you compute: what does Milestone 3 - race gradient descent come out to? A rough magnitude and the right form is enough — the point is to have something concrete to be wrong about.
Correct: GD takes 1147 iterations — about 164× Newton's 7
Why: A prediction you can defend turns the computation into a check rather than a leap of faith — and an answer that contradicts it is caught on the spot. Same solution [0.48834, −1.34714, 2.672486], but GD's linear convergence needs 1147 steps to reach ‖grad‖ < 1e−8 versus Newton's 7.
Worked example
Your turn: run GD on the same loss to the same tolerance and count iterations. Predict the ratio to Newton before you print (order of magnitude).
Hint: GD loop wg = wg − 0.5*g until norm(g) < 1e−8; it takes over a thousand steps. Confirm it lands on the same w.
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
wg = np.zeros(3)
for it in range(1000000):
p = sig(Xb @ wg)
g = Xb.T @ (p - y) / n
wg = wg - 0.5 * g
if np.linalg.norm(g) < 1e-8:
break
print('GD iters:', it+1) # 1147
print('w =', wg.round(6)) # same as NewtonGD takes 1147 iterations — about 164× Newton's 7
Why: Same solution [0.48834, −1.34714, 2.672486], but GD's linear convergence needs 1147 steps to reach ‖grad‖ < 1e−8 versus Newton's 7. That ~160× gap is the whole point of using curvature.
| method | iterations | final w |
|---|---|---|
| Newton (IRLS) | 7 | [0.48834, −1.34714, 2.672486] |
| gradient descent | 1147 | [0.48834, −1.34714, 2.672486] |
| ratio | ≈ 164× | identical solution |
Trade off
Comparison matrix
From Milestone 3 - race gradient descent: every row here is a choice with a cost. Fill the final w column, then say which row you would actually pick and what you give up for it.
| method | iterations | final w |
|---|---|---|
| Newton (IRLS) | 7 | [0.48834, −1.34714, 2.672486] |
| gradient descent | 1147 | [0.48834, −1.34714, 2.672486] |
| ratio | ≈ 164× | identical solution |
Pattern
Predict first
The table runs: w | [0.48834, −1.34714, 2.672486] · Newton iterations | 7 · GD iterations (same tol) | 1147
In The full program, given the rows so far: what is the next one — the row where output is matches sklearn??
Correct: matches sklearn? | yes (atol 1e−4)
| output | value (verified) |
|---|---|
| w | [0.48834, −1.34714, 2.672486] |
| Newton iterations | 7 |
| GD iterations (same tol) | 1147 |
| matches sklearn? | yes (atol 1e−4) |
Why: The relationship between the columns, not the individual numbers, is what generates the next row. Gradient, Hessian, solve, subtract, check the step norm — the entire second-order method.
Concept
import numpy as np
np.random.seed(0)
n = 200
Xb = np.c_[np.ones(n), np.random.randn(n, 2)]
w_true = np.array([0.5, -1.5, 2.0])
y = (np.random.rand(n) < 1/(1+np.exp(-(Xb@w_true)))).astype(float)
sig = lambda z: 1/(1+np.exp(-z))
def newton_logistic(X, y, iters=100):
w = np.zeros(X.shape[1])
for it in range(iters):
p = sig(X @ w); m = len(y)
g = X.T @ (p - y) / m
H = (X.T * (p*(1-p))) @ X / m + 1e-8*np.eye(X.shape[1])
step = np.linalg.solve(H, g) # solve, never invert
w = w - step
if np.linalg.norm(step) < 1e-8:
return w, it+1
return w, iters
w, k = newton_logistic(Xb, y)
print('w =', w.round(6), ' iters =', k) # [...] 7The whole method is a dozen lines
Why: Gradient, Hessian, solve, subtract, check the step norm — the entire second-order method. Returns w = [0.48834, −1.34714, 2.672486] in 7 iterations.
| output | value (verified) |
|---|---|
| w | [0.48834, −1.34714, 2.672486] |
| Newton iterations | 7 |
| GD iterations (same tol) | 1147 |
| matches sklearn? | yes (atol 1e−4) |
Comparison
Comparison matrix
From The full program: refill the value (verified) column from what you know. The rest of the table is as it appeared.
| output | value (verified) |
|---|---|
| w | [0.48834, −1.34714, 2.672486] |
| Newton iterations | 7 |
| GD iterations (same tol) | 1147 |
| matches sklearn? | yes (atol 1e−4) |
Concept
Slides closed, out loud: explain (1) why Newton solves a quadratic in one step, (2) what the Hessian's eigenvalue signs say about a critical point, and (3) why we don't run full Newton on a million-parameter net.
Stretch: show a gradient-descent step is a first-order Taylor model and Newton's is second-order; then prove the logistic Hessian XᵀWX equals the Fisher information matrix — so Newton and natural gradient coincide here. Second-order ideas return in natural-gradient and SAM (Week 35).
Connect it up
Draw it
One page, no notation unless you need it: draw how these connect — Why curvature: the second-order model · The Hessian and definiteness · Newton's method from the model · Newton is exact on a quadratic · Quadratic convergence · Newton in practice: logistic regression. Put an arrow wherever one of them is what makes another possible, and label the arrow with why.
Recap
H = ∇²f and classify a critical point by its eigenvalue signs — all + → min, all − → max, mixed → saddleθ ← θ − H⁻¹∇f by minimizing the second-order model (Hs = −g)[5,−3] → [0.2, 0.4]) and quadratically convergent near a minimum7 iters vs GD's 1147, same w, matched to sklearn| idea | the one thing to remember |
|---|---|
| Hessian | symmetric; PD → min, indefinite → saddle |
| Newton step | solve Hs = −g → θ −= H⁻¹∇f (never inv) |
| speed | 1 step on a quadratic; quadratic convergence (digits double) |
| cost & safety | O(d³)/step; damp or trust-region far from a min |
| at scale | L-BFGS / natural gradient (Fisher = XᵀWX here) |
Want this taught 1-on-1? Alexander tutors Machine Learning — $55/session, free consultation.