4.4.3 Newton's method in several variables

The Newton step

Gradient descent uses only the slope, so it cannot tell a narrow valley from a round bowl. Newton's method uses the curvature too. Near the current point $\mathbf{x}_k$ , lesson 4.3.1 copies $f$ by a quadratic in the step $\mathbf{h}$ :

\begin{gathered} f(\mathbf{x}_k + \mathbf{h}) \approx f(\mathbf{x}_k) + \nabla f \cdot \mathbf{h} + \tfrac{1}{2}\,\mathbf{h}^\top H\,\mathbf{h} \end{gathered}

with $\nabla f$ and $H$ taken at $\mathbf{x}_k$ . This is lesson 4.4.1's quadratic in $\mathbf{h}$ , with $A = H$ and $\mathbf{b} = -\nabla f$ .

If $H$ is positive definite, the copy's minimum is where $H\mathbf{h} = -\nabla f$ . Newton's method jumps to that minimum, then builds a new copy there.

Rule: Newton's method

\begin{gathered} \text{solve } H(\mathbf{x}_k)\,\mathbf{h}_k = -\nabla f(\mathbf{x}_k) \\[4pt] \mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{h}_k \end{gathered}

It is often written $\mathbf{x}_{k+1} = \mathbf{x}_k - H^{-1}\nabla f$ . Compute it with np.linalg.solve (lesson 3.3.1): forming $H^{-1}$ costs more and adds rounding. The curvature sets the length of the step, so there is no $\eta$ to choose.

In one variable, $H$ is the number $f''$ and $\nabla f$ is $f'$ , so the step is $h = -f'/f''$ . That is lesson 2.4.3's Newton's method for a critical point.

A quadratic takes one step

For $f = \frac{1}{2}\mathbf{x}^\top A\mathbf{x} - \mathbf{b}^\top\mathbf{x}$ the copy is exact, so one step lands on the minimum from any start. Take lesson 4.4.1's $x^2 + xy + y^2 - 3x$ , with $H = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$ . At $(0, 1)$ the gradient is $(-2, 2)$ , so the step solves

\begin{aligned} 2h_1 + h_2 &= 2 \\ h_1 + 2h_2 &= -2 \end{aligned}

which gives $\mathbf{h} = (2, -2)$ and $\mathbf{x}_1 = (2, -1)$ , the minimum. From $(0, 0)$ the step is $(2, -1)$ , and from $(5, 5)$ it is $(-3, -6)$ : both land on $(2, -1)$ . Gradient descent needed 22 steps at its best step size (lesson 4.4.2).

Lesson 2.4.3's quartic, in two variables

Take $f(x, y) = x^4 - 3x^2 + x + y^2$ of lesson 4.4.1, with $\nabla f = (4x^3 - 6x + 1,\ 2y)$ and $H = \begin{pmatrix} 12x^2 - 6 & 0 \\ 0 & 2 \end{pmatrix}$ . This $H$ is diagonal, so the Newton step splits into two one-variable steps.

For $y$ the step is $-2y/2 = -y$ : $y$ lands on 0 at once, because $y^2$ is a quadratic. For $x$ it is lesson 2.4.3's step $-f'(x)/f''(x)$ .

From $(1, 0.5)$ the $x$ values are 1.166667, 1.132616, 1.130905 and 1.130901, exactly those of lesson 2.4.3. The end is the local minimum $(1.130901, 0)$ .

H\mathbf{h} = -\nabla f

import numpy as np def newton_step(grad, hess, x): return x + np.linalg.solve(hess(x), -grad(x)) def grad(p): return np.array([2 * p[0] + p[1] - 3, p[0] + 2 * p[1]]) def hess(p): return np.array([[2.0, 1.0], [1.0, 2.0]]) for start in [(0, 1), (0, 0), (5, 5)]: x = np.array(start, dtype=float) print(start, newton_step(grad, hess, x) - x, newton_step(grad, hess, x)) def q_grad(p): return np.array([4 * p[0]**3 - 6 * p[0] + 1, 2 * p[1]]) def q_hess(p): return np.array([[12 * p[0]**2 - 6, 0.0], [0.0, 2.0]]) x = np.array([1.0, 0.5]) for k in range(4): x = newton_step(q_grad, q_hess, x) print(k + 1, np.round(x, 6))

⌘/Ctrl + Enter runs the cell

How fast: the digits double

Lesson 4.4.2's curved valley, $f(x, y) = (1 - x)^2 + 5(y - x^2)^2$ , has its minimum at $(1, 1)$ and a condition number of about 133 there. Its gradient and Hessian are

\begin{gathered} \nabla f = \begin{pmatrix} -2(1 - x) - 20x(y - x^2) \\ 10(y - x^2) \end{pmatrix} \\[6pt] H = \begin{pmatrix} 2 - 20(y - x^2) + 40x^2 & -20x \\ -20x & 10 \end{pmatrix} \end{gathered}

At $\mathbf{x}_0 = (0.5, 0)$ , $\nabla f = (1.5, -2.5)$ and $H = \begin{pmatrix} 17 & -10 \\ -10 & 10 \end{pmatrix}$ , with eigenvalues 2.9052 and 24.0948. Solving $H\mathbf{h} = -\nabla f$ gives $\mathbf{x}_1 = (0.642857, 0.392857)$ . The copy promised $f = 0.1786$ there; the true value is 0.1296.

step $k$	$\mathbf{x}_k$	$f(\mathbf{x}_k)$	error $\lVert\mathbf{x}_k - (1, 1)\rVert$	correct digits
0	$(0.5, 0)$	0.5625	1.118
1	$(0.642857, 0.392857)$	0.1296	0.704	0.2
2	$(0.939467, 0.794621)$	0.04236	0.214	0.7
3	$(0.971669, 0.943104)$	0.000808	0.0636	1.2
4	$(0.999709, 0.998632)$	$3.175 \times 10^{-6}$	0.0014	2.9
5	$(0.999998, 0.999995)$	$5.18 \times 10^{-12}$	$5.15 \times 10^{-6}$	5.3
6	$(1, 1)$		$9.12 \times 10^{-12}$	11.0

The correct digits, $-\log_{10}$ of the error, grow slowly at first and then about double: 1.2, 2.9, 5.3, 11.0. This is the quadratic convergence of lesson 2.4.3. Near a minimum where $H$ is positive definite, each error is roughly a constant times the square of the last.

From the same start, gradient descent with $\eta = 0.02$ needs 2186 steps to bring $\lVert\nabla f\rVert$ below $10^{-8}$ ; Newton needs 6. The valley's narrowness does not slow Newton down: $H$ measures the curvature in every direction and sizes the step for each. That is the rescaling lesson 4.4.2 asked for.

Each Newton step jumps to the minimum of the quadratic copy, whose level curves are the dashed ellipses. Where the Hessian is indefinite, plain Newton heads for the nearest critical point, here a saddle; a safeguarded step goes downhill instead.

(0.5, 0)

⌘/Ctrl + Enter runs the cell

When H is not positive definite

Newton's method solves $\nabla f = \mathbf{0}$ , so it heads for the nearest critical point of any kind. Where $H$ is not positive definite, the copy has no minimum, and two things go wrong.

It can land on a saddle or a maximum. For the quartic at $(0, 0.5)$ , $H = \begin{pmatrix} -6 & 0 \\ 0 & 2 \end{pmatrix}$ is indefinite. The $x$ steps are lesson 2.4.3's from 0: 0.166667, 0.169935, 0.169938. They reach the saddle $(0.169938, 0)$ of lesson 4.4.1, which is no minimum.

Its step can go uphill. When $H$ is positive definite, $\nabla f \cdot \mathbf{h} = -\nabla f^\top H^{-1}\nabla f < 0$ , so the step points downhill (lesson 4.2.1). Otherwise it may not. On the curved valley at $(-1.5, 2.5)$ , $H = \begin{pmatrix} 87 & 30 \\ 30 & 10 \end{pmatrix}$ has determinant $-30$ and eigenvalues $-0.3083$ and 97.3083.

There the Newton step goes to $(-3.166667, 7.25)$ , where $f = 55.94$ , up from 6.5625. Plain Newton does reach $(1, 1)$ after 8 steps from this start, but it passes $f = 611.5$ on the way. That was luck.

Lesson 4.4.2's explorer has a Newton's method switch. Turn it on at the valley's start $(-1.5, 2.5)$ and set 1 step, then 8. The first step leaves the picture uphill; by step 8 the path is back at $(1, 1)$ .

Rule: A safeguarded Newton step

\begin{gathered} H \text{ positive definite: solve } H\mathbf{h} = -\nabla f \\[4pt] \text{otherwise: } \mathbf{h} = -\nabla f \\[4pt] \text{halve } \mathbf{h} \text{ until } f(\mathbf{x} + \mathbf{h}) < f(\mathbf{x}) \end{gathered}

Test $H$ with np.linalg.cholesky, which fails unless $H$ is positive definite (lesson 3.7.2). Every step now lowers $f$ , so the method heads for a minimum. Near one, $H$ is positive definite and the full Newton steps take over.

With the safeguard, the quartic from $(0, 0.5)$ takes one gradient step, to $(-1, -0.5)$ . Then the Newton steps give $x = -1.5$ , $-1.333333$ , $-1.301932$ , $-1.300841$ : lesson 2.4.3's steps from $-1$ , ending at the global minimum $(-1.30084, 0)$ . The curved valley from $(-1.5, 2.5)$ takes 12 steps, all downhill.

(-1.5, 2.5)

import numpy as np def f(p): return (1 - p[0])**2 + 5 * (p[1] - p[0]**2)**2 def grad(p): x, y = p return np.array([-2 * (1 - x) - 20 * x * (y - x**2), 10 * (y - x**2)]) def hess(p): x, y = p return np.array([[2 - 20 * (y - x**2) + 40 * x**2, -20 * x], [-20 * x, 10.0]]) x = np.array([-1.5, 2.5]) h = np.linalg.solve(hess(x), -grad(x)) print(np.round(np.linalg.eigvalsh(hess(x)), 4), round(grad(x) @ h, 4), x + h, f(x), round(f(x + h), 2)) def safe_newton(f, grad, hess, x, tol=1e-10, max_steps=100): for k in range(max_steps): g = grad(x) if np.linalg.norm(g) < tol: return x, k try: np.linalg.cholesky(hess(x)) h = np.linalg.solve(hess(x), -g) except np.linalg.LinAlgError: h = -g while not f(x + h) < f(x): h = h / 2 x = x + h raise RuntimeError("no convergence") x, steps = safe_newton(f, grad, hess, np.array([-1.5, 2.5])) print(np.round(x, 6), steps) def quartic(p): return p[0]**4 - 3 * p[0]**2 + p[0] + p[1]**2 def q_grad(p): return np.array([4 * p[0]**3 - 6 * p[0] + 1, 2 * p[1]]) def q_hess(p): return np.array([[12 * p[0]**2 - 6, 0.0], [0.0, 2.0]]) x, steps = safe_newton(quartic, q_grad, q_hess, np.array([0.0, 0.5])) print(np.round(x, 6) + 0.0, steps, round(quartic(x), 4))

⌘/Ctrl + Enter runs the cell

The cost of a step

Newton takes few steps, but each one costs more. A gradient-descent step needs the gradient, $n$ numbers.

A Newton step also needs the Hessian, which has $\frac{n(n+1)}{2}$ different entries because it is symmetric. It then solves an $n \times n$ system, about $\frac{n^3}{3}$ multiplications by elimination (lesson 3.3.1).

variables $n$	gradient entries	Hessian entries	solving, about
10	10	55	333
100	100	5,050	$3.3 \times 10^5$
1,000	1,000	500,500	$3.3 \times 10^8$
1,000,000	$10^6$	$5 \times 10^{11}$	$3.3 \times 10^{17}$

With up to a few thousand variables, Newton's handful of steps wins. With millions, as in a large machine-learning model, the Hessian cannot even be stored, and gradient methods are the only choice.

Between the two sit quasi-Newton methods. The best known, BFGS, builds an estimate of $H$ from the gradients of earlier steps, so it needs only gradients and still converges fast.

scipy.optimize.minimize

minimize(f, x0) from scipy.optimize returns a result with the minimiser .x, the value .fun, the number of iterations .nit, and .success with a .message. Without bounds or constraints it uses BFGS. Given no gradient, it estimates one by finite differences, as in lesson 2.8.1; jac=grad passes yours, which saves evaluations of $f$ .

method="Newton-CG" with hess=hess takes Newton steps. It solves $H\mathbf{h} = -\nabla f$ by an iterative method instead of elimination, which suits a large $H$ .

On the curved valley from $(-1.5, 2.5)$ , BFGS takes 17 iterations: 66 evaluations of $f$ without a gradient, and 22 with one. Newton-CG takes 23. The counts can differ slightly with the SciPy version.

minimize finds a local minimum near its start. On the quartic from $(1, 0.5)$ it returns $(1.130901, 0)$ , with $f = -1.0702$ , the minimum that is not global (lesson 4.4.1).

From $(-1, 0.5)$ it returns the global one, $-3.5139$ . Case study 4.4.5 uses minimize on a yield curve.

Try it The first three results are the curved valley: BFGS with an estimated gradient, BFGS with jac, and Newton-CG with the Hessian too. Each line prints success, the answer and the counts. The last loop runs BFGS on the quartic from two starts.

import numpy as np from scipy.optimize import minimize def f(p): return (1 - p[0])**2 + 5 * (p[1] - p[0]**2)**2 def grad(p): x, y = p return np.array([-2 * (1 - x) - 20 * x * (y - x**2), 10 * (y - x**2)]) def hess(p): x, y = p return np.array([[2 - 20 * (y - x**2) + 40 * x**2, -20 * x], [-20 * x, 10.0]]) x0 = np.array([-1.5, 2.5]) r = minimize(f, x0) print(r.success, np.round(r.x, 4), r.nit, r.nfev) r = minimize(f, x0, jac=grad) print(r.success, np.round(r.x, 4), r.nit, r.nfev) r = minimize(f, x0, jac=grad, hess=hess, method="Newton-CG") print(r.success, np.round(r.x, 4), r.nit) def quartic(p): return p[0]**4 - 3 * p[0]**2 + p[0] + p[1]**2 for start in [(1, 0.5), (-1, 0.5)]: r = minimize(quartic, np.array(start, dtype=float)) print(start, np.round(r.x, 4) + 0.0, round(r.fun, 4))

⌘/Ctrl + Enter runs the cell

Common mistakes

Do not invert $H$ . Solve $H\mathbf{h} = -\nabla f$ with np.linalg.solve.
Mind the sign. The step solves $H\mathbf{h} = -\nabla f$ ; with $+\nabla f$ it goes the wrong way.
Newton finds critical points, not minima. From $(0, 0.5)$ on the quartic it lands on a saddle; check $H$ at the answer (lesson 4.4.1).
Away from the minimum, $H$ may be indefinite and the step may go uphill. Test $H$ , fall back to a gradient step, and halve until $f$ falls.
A wrong Hessian spoils the speed. Check it against central differences of the gradient (lesson 4.1.2).
Few steps is not the same as cheap. Each step forms $H$ and solves a system; with many variables, use BFGS or gradient methods.
success is not a proof. minimize reports a local minimum near its start; try other starts.

Newton's method in several variables

0 of 12 correct

1.
One Newton step on $x^2 + xy + y^2 - 3x$ from $(0, 0)$ solves $\begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\mathbf{h} = (3, 0)$ . What is $h_1$ ?
2.
Solve $H\mathbf{h} = -\nabla f$ for $H = \begin{pmatrix} 4 & 0 \\ 0 & 1 \end{pmatrix}$ and $\nabla f = (8, -3)$ . Give $h_1$
3.
For the same $H$ and $\nabla f$ , give $h_2$
4.
Solve $H\mathbf{h} = -\nabla f$ for $H = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}$ and $\nabla f = (3, 3)$ . Give $h_1$
5.
For $x^4 - 3x^2 + x + y^2$ , take one Newton step from $(1, 0.5)$ . What is the new $y$ ?
6.
Near a minimum, two Newton errors are 0.01 and 0.0003. If each error is about a constant times the square of the last, estimate the next
7.
With $H = \begin{pmatrix} 2 & 0 \\ 0 & 4 \end{pmatrix}$ and $\nabla f = (2, 4)$ , find $\nabla f \cdot \mathbf{h}$ for the Newton step
8.
At a point the Hessian has eigenvalues 5 and $-1$ . Is the plain Newton step sure to go downhill? Type 1 for yes, 0 for no
9.
Plain Newton on $x^2 - y^2$ from $(3, 2)$ : how many steps does it take to reach the saddle $(0, 0)$ ?
10.
A function of 100 variables has a symmetric Hessian. How many different entries does it have?
11.
Elimination on an $n \times n$ system takes about $n^3/3$ multiplications. About how many for $n = 300$ ?
12.
minimize on $x^4 - 3x^2 + x + y^2$ from $(1, 0.5)$ returns $f = -1.0702$ with success. Is that the global minimum? Type 1 for yes, 0 for no

How to type answers

Powers with ^: x^7, 2^(1/3), a^-2.
Multiply with * or by writing side by side: 3x^2, 2(x+1).
ln(x) is the natural log, log(x, 2) is log base 2, e is Euler's number, sqrt(x) is a square root.
sin(x), cos(x) and tan(x) take radians; write pi for π and sin(x)^2 for sin²x.
Your first answer to each question is saved: it counts towards your score and brings the question back for review.
Exact answers like ln(20)/ln(3) are best. Decimals need about 4 significant figures.

H(\mathbf{x}_k)\mathbf{h} = -\nabla f(\mathbf{x}_k)

Newton's method in several variables

The Newton step

Rule: Newton's method

A quadratic takes one step

Lesson 2.4.3's quartic, in two variables

How fast: the digits double

When H is not positive definite

Rule: A safeguarded Newton step

The cost of a step

scipy.optimize.minimize

Common mistakes

Newton's method in several variables

Summary card

Free further reading