6 Second-Order Methods
The previous chapter focused on optimization methods that involve first-orderapproximations of the objective function using the gradient. This chapter focuseson leveraging second-order approximations that use the second derivative in univariate optimization or the Hessian in multivariate optimization to direct the
search. This additional information can help improve the local model used for
informing the selection of directions and step lengths in descent algorithms.

6.1 Newton’s Method
Knowing the function value and gradient for a design point can help determine
the direction to travel, but this first-order information does not directly help
determine how far to step to reach a local minimum. Second-order information,on the other hand, allows us to make a quadratic approximation of the objectivefunction and approximate the right step size to reach a local minimum as shownin figure 6.1. As we have seen with quadratic fit search in chapter 3, we cananalytically obtain the location where a quadratic approximation has a zero
gradient. We can then use that location as the next iteration to approach a local
minimum.
In univariate optimization, the quadratic approximation about a point x (k)comes from the second-order Taylor expansion:q ( x ) = f ( x (k) ) + ( x − x (k) ) f ′ ( x (k) ) +( x − x (k) )2 ′′ (k)f (x )
2(6.1)

Figure 6.1. A comparison of firstorder and second-order approximations. Bowl-shaped quadratic
approximations have unique locations where the derivative is zero.

f

f

88 c ha p te r 6 . se c on d -ord e r me thod sxx

Setting the derivative to zero and solving for the root yields the update equation
for Newton’s method:∂q( x ) = f ′ ( x (k) ) + ( x − x (k) ) f ′′ ( x (k) ) = 0∂xf ′ ( x (k) )x ( k +1) = x ( k ) −f ′′ ( x (k) )

(6.2)
(6.3)

This update is shown in figure 6.2.

f

f′

Figure 6.2. Newton’s method canbe interpreted as a root-finding
method applied to f ′ that iteratively improves a univariate design
point by taking the tangent line at( x, f ′ ( x )), finding the intersection
with the x-axis, and using that x
value as the next design point.

0x ( k +1)xx (k)x ( k +1)xx (k)The update rule in Newton’s method involves dividing by the second derivative.
The update is undefined if the second derivative is zero, which occurs whenthe quadratic approximation is a line. Instability also occurs when the second
derivative is very close to zero, in which case the next iterate will lie very far fromthe current design point, far from where the local quadratic approximation isvalid. Poor local approximations can lead to poor performance with Newton’s
method. Figure 6.3 shows three kinds of failure cases.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

6 . 1. n e w ton ’s me thod

Oscillation

89

Negative f ′′fOvershoot

x

x (k)

x ( k +1)x ( k +1)xx

Newton’s method does tend to converge quickly when in a bowl-like region
that is sufficiently close to a local minimum. It has quadratic convergence, meaning
the difference between the minimum and the iterate is approximately squared
with every iteration. This rate of convergence holds for Newton’s method startingfrom x (1) within a distance δ of a root x ∗ if1• f ′′ ( x ) 6= 0 for all points in I,• f ′′′ ( x ) is continuous on I, and
•

x (k)

′′′ (1)f ′′′ ( x ∗ )1 f (x )2 | f ′′ ( x (1) ) | < c | f ′′ ( x ∗ ) | for some c < ∞Figure 6.3. Examples of failurecases with Newton’s method.

The final condition enforces sufficient closeness, ensuring that thefunction is sufficiently approximated by the Taylor expansion.
J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, 3rd ed.
Springer, 2002.

1

for an interval I = [ x ∗ − δ, x ∗ + δ]. The final condition guards against overshoot.
Newton’s method can be extended to multivariate optimization (algorithm 6.1).
The multivariate second-order Taylor expansion at x(k) is:1
f ( x ) ≈ q (x ) = f ( x( k ) ) + ( g( k ) ) ⊤ (x − x( k ) ) + (x − x( k ) ) ⊤ H( k ) (x − x( k ) )
2

(6.4)

where g(k) and H(k) are the gradient and Hessian at x(k) , respectively.
We evaluate the gradient and set it to zero:∇ q (x ) = g( k ) + H( k ) (x − x( k ) ) = 0(6.5)

We then solve for the next iterate, thereby obtaining Newton’s method in multivariate form:x( k +1) = x( k ) − (H( k ) ) −1 g( k )(6.6)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

90

c hap te r 6 . se c on d -ord e r me thod s

If f is quadratic and its Hessian is positive definite, then the update converges
to the global minimum in one step. For general functions, Newton’s method
is often terminated once x ceases to change by more than a given tolerance.2
Example 6.1 shows how Newton’s method can be used to minimize a function.
With x(1) = [9, 8], we will use Newton’s method to minimize Booth’s function:f (x) = ( x1 + 2x2 − 7)2 + (2x1 + x2 − 5)2

Termination conditions for descent methods are given in chapter 4.

2

Example 6.1. Newton’s method
used to minimize Booth’s function;
see appendix B.2.

The gradient of Booth’s function is:∇ f (x) = [10x1 + 8x2 − 34, 8x1 + 10x2 − 38]
The Hessian of Booth’s function is:"10
H(x) =
8

8
10#

The first iteration of Newton’s method yields:
x

(2)

=x

(1)



− H

(1)

" # "
10
9
−
=
8
8

" # "
10
9
−
g =
8
8
# " #
# −1 "
1
120
8
=
3
114
10

 −1

(1)

8
10

# −1 "

10 · 9 + 8 · 8 − 34
8 · 9 + 10 · 8 − 38#

The gradient at x(2) is zero, so we have converged after a single iteration. The
Hessian is positive definite everywhere, so x(2) is the global minimum.
Newton’s method can also be used to supply a descent direction to line searchor can be modified to use a step factor.3 Smaller steps toward the minimum orline searches along the descent direction can increase the method’s robustness.
The descent direction is:4d(k) = −(H(k) )−1 g(k)

(6.7)

3

See chapter 5.

The descent direction given by
Newton’s method is similar to thenatural gradient or covariant gradient. S. Amari, “Natural Gradient Works Efficiently in Learning,”
Neural Computation, vol. 10, no. 2,
pp. 251–276, 1998.

4© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

6 . 2. se c a n t me thod 91

Algorithm 6.1. Newton’s method,which takes the gradient of thefunction ∇f, the Hessian of the objective function H, an initial point x,
a step size tolerance ϵ, and a maximum number of iterations k_max.

function newtons_method(∇f, H, x, ϵ, k_max)k, Δ = 1, fill(Inf, length(x))
while norm(Δ) > ϵ && k ≤ k_max
Δ = H(x) \ ∇f(x)x -= Δk += 1end
return x
end

6.2 Secant Method
Newton’s method for univariate function minimization requires the first andsecond derivatives f ′ and f ′′ . In many cases, f ′ is known but the second derivative
is not. The secant method (algorithm 6.2) applies Newton’s method using estimatesof the second derivative and thus only requires f ′ . This property makes the secant
method more convenient to use in practice.
The secant method uses the last two iterates to approximate the second derivative:f ′ ( x ( k ) ) − f ′ ( x ( k −1) )f ′′ ( x (k) ) ≈(6.8)
x ( k ) − x ( k −1)
This estimate is substituted into Newton’s method:x ( k +1) ← x ( k ) −x ( k ) − x ( k −1)
f ′ ( x (k) )f ′ ( x ( k ) ) − f ′ ( x ( k −1) )(6.9)

The secant method requires an additional initial design point. It suffers fromthe same problems as Newton’s method and may take more iterations to converge
due to approximating the second derivative.

6.3 Quasi-Newton Methods
Just as the secant method approximates f ′′ in the univariate case, quasi-Newton
methods approximate the inverse Hessian. Quasi-Newton method updates havethe form:x( k +1) ← x( k ) − α ( k ) Q( k ) g( k )(6.10)

where α(k) is a scalar step factor and Q(k) approximates the inverse of the Hessianat x(k) .

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

92

c ha p te r 6 . se c on d -ord e r me thod s

Algorithm 6.2. The secant method
for univariate function minimization. The inputs are the first derivative f′ of the target function, twoinitial points x0 and x1, and thedesired tolerance ϵ. The final xcoordinate is returned.

function secant_method(f′, x0, x1, ϵ)
g0 = f′(x0)Δ = Infwhile abs(Δ) > ϵ
g1 = f′(x1)Δ = (x1 - x0)/(g1 - g0)*g1
x0, x1, g0 = x1, x1 - Δ, g1
end
return x1
end

These methods typically set Q(1) to the identity matrix, and they then apply updates to reflect information learned with each iteration. To simplify the equations
for the various quasi-Newton methods, we define the following:γ( k +1) ≡ g( k +1) − g( k )δ( k +1)≡x( k +1)−x(k)(6.11)(6.12)The Davidon-Fletcher-Powell (DFP) method (algorithm 6.3) uses:5
Q ← Q−Qγγ⊤ Q δδ⊤
+ ⊤
γ⊤ Qγ
δ γ

(6.13)

where all terms on the right hand side are evaluated at the same iteration.
The update for Q in the DFP method has three properties:
1. Q remains symmetric and positive definite.
2. If f (x) = 21 x⊤ Ax + b⊤ x + c, then Q = A−1 . Thus the DFP has the same
convergence properties as the conjugate gradient method.
3. For high-dimensional problems, storing and updating Q can be significant
compared to other methods like the conjugate gradient method.
An alternative to DFP, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method(algorithm 6.4), uses:6!
!
γ⊤ Qγ δδ⊤δγ⊤ Q + Qγδ⊤+ 1+ ⊤(6.14)Q ← Q−
δ⊤ γ
δ γ
δ⊤ γ

5
The original concept was presented in a technical report, W. C.
Davidon, “Variable Metric Method
for Minimization,” Argonne National Laboratory, Tech. Rep. ANL5990, 1959. It was later published:
W. C. Davidon, “Variable Metric
Method for Minimization,” SIAM
Journal on Optimization, vol. 1, no. 1,
pp. 1–17, 1991. The method was
modified by R. Fletcher and M. J. D.
Powell, “A Rapidly Convergent Descent Method for Minimization,”
The Computer Journal, vol. 6, no. 2,
pp. 163–168, 1963.

R. Fletcher, Practical Methods of Optimization, 2nd ed. Wiley, 1987.

6© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

6 . 3. q u a si -n e w ton me thod s 93

mutable struct DFP <: DescentMethodQ
endfunction init!(M::DFP, f, ∇f, x)
m = length(x)M.Q = Matrix(1.0I, m, m)
return M
endfunction step!(M::DFP, f, ∇f, x)
Q, g = M.Q, ∇f(x)x′ = line_search(f, x, -Q*g)g′ = ∇f(x′)δ = x′ - x
γ = g′ - g
Q[:] = Q - Q*γ*γ'*Q/(γ'*Q*γ) + δ*δ'/(δ'*γ)
return x′
end

Algorithm 6.3. The DavidonFletcher-Powell descent method.

mutable struct BFGS <: DescentMethodQ
endfunction init!(M::BFGS, f, ∇f, x)
m = length(x)M.Q = Matrix(1.0I, m, m)
return M
endfunction step!(M::BFGS, f, ∇f, x)
Q, g = M.Q, ∇f(x)x′ = line_search(f, x, -Q*g)g′ = ∇f(x′)δ = x′ - x
γ = g′ - g
Q[:] = Q - (δ*γ'*Q + Q*γ*δ')/(δ'*γ) +(1 + (γ'*Q*γ)/(δ'*γ))[1]*(δ*δ')/(δ'*γ)
return x′
end

Algorithm 6.4. The BroydenFletcher-Goldfarb-Shanno descentmethod.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

94

c ha p te r 6 . se c on d -ord e r me thod s

BFGS does better than DFP with approximate line search but still uses ann × n dense matrix. For very large problems where space is a concern, the Limitedmemory BFGS method (algorithm 6.5), or L-BFGS, can be used to approximate
BFGS.7 L-BFGS stores the last m values for δ and γ rather than the full inverse
Hessian, where i = 1 indexes the oldest value and i = m indexes the most recent.
The process for computing the descent direction d at x begins by computing
q(m) = ∇ f (x). The remaining vectors q(i) for i from m − 1 down to 1 are computed
using

⊤
δ( i +1) q( i +1)γ( i +1)(6.15)q( i ) = q( i +1) −⊤γ( i +1) δ( i +1)

7
J. Nocedal, “Updating QuasiNewton Matrices with Limited
Storage,” Mathematics of Computation, vol. 35, no. 151, pp. 773–782,
1980.

These vectors are used to compute another m + 1 vectors, starting withz(0) =γ( m ) ⊙ δ( m ) ⊙ q( m )⊤γ( m ) γ( m )(6.16)and proceeding with z(i) for i from 1 to m according to



⊤

⊤
( i −1)( i −1)( i −1)( i −1)δqzγ
z( i ) = z( i −1) + δ( i −1)  ⊤ ( i −1) − ⊤ ( i −1) γ( i −1) δγ( i −1) δ(6.17)The descent direction is d = −z(m) .
For minimization, the inverse Hessian Q must remain positive definite. The
initial Hessian is often set to the diagonal of⊤γ(1) δ(1)Q(1) =
⊤γ(1) γ(1)(6.18)

Computing the diagonal for the above expression and substituting the result intoz(1) = Q(1) q(1) results in the equation for z(1) .
The quasi-Newton methods discussed in this section are compared in figure 6.4.
They often perform quite similarly.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

6 . 4 . su mma ry 95

Figure 6.4. Several quasi-Newton
methods compared on the Rosenbrock function; see appendix B.6.
All methods have nearly identical
updates, with L-BFGS noticeably
deviating only when its history, m,
is 1.

x2

DFP
BFGS
L-BFGS (m = 3)
L-BFGS (m = 2)
L-BFGS (m = 1)

x1

6.4 Summary• Incorporating second-order information in descent methods often speeds convergence.
• Newton’s method is a root-finding method that leverages second-order information to quickly descend to a local minimum.
• The secant method and quasi-Newton methods approximate Newton’s method
when the second-order information is not directly available.

6.5

Exercises

Exercise 6.1. What advantage does second-order information provide about
convergence that first-order information lacks?
Exercise 6.2. When finding roots in one dimension, when would we use Newton’s
method instead of the bisection method?
Exercise 6.3. Apply Newton’s method to f ( x ) = x2 from a starting point of yourchoice. How many steps do we need to converge?

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

96 c ha p te r 6 . se c on d -ord e r me thod smutable struct LimitedMemoryBFGS <: DescentMethod
m
δs
γs
qs
end
function init!(M::LimitedMemoryBFGS, f, ∇f, x)
M.δs = []M.γs = []M.qs = []return M
endfunction step!(M::LimitedMemoryBFGS, f, ∇f, x)
δs, γs, qs, g = M.δs, M.γs, M.qs, ∇f(x)m = length(δs)
if m > 0q = g
for i in m : -1 : 1
qs[i] = copy(q)q -= (δs[i]⋅q)/(γs[i]⋅δs[i])*γs[i]endz = (γs[m] .* δs[m] .* q) / (γs[m]⋅γs[m])for i in 1 : m
z += δs[i]*(δs[i]⋅qs[i] - γs[i]⋅z)/(γs[i]⋅δs[i])endx′ = line_search(f, x, -z)else
x′ = line_search(f, x, -g)
end
g′ = ∇f(x′)
push!(δs, x′ - x); push!(γs, g′ - g)push!(qs, zeros(length(x)))
while length(δs) > M.m
popfirst!(δs); popfirst!(γs); popfirst!(qs)end
return x′
end

Algorithm 6.5.
The Limitedmemory BFGS descent method,which avoids storing the approximate inverse Hessian. The
parameter m determines the history
size. The LimitedMemoryBFGS typealso stores the step differences
δs, the gradient changes γs, andstorage vectors qs.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

6 . 5 . e x e rc i se s 97

Exercise 6.4. Apply Newton’s method to f ( x ) = 21 x⊤ Hx starting from x(1) =[1, 1]. What have you observed? Use H as follows:#"
10(6.19)
H=
0 1000
Next, apply gradient descent to the same optimization problem by steppingwith the unnormalized gradient. Do two steps of the algorithm. What have you
observed? Finally, apply the conjugate gradient method. How many steps do you
need to converge?
Exercise 6.5. Compare Newton’s method and the secant method on f ( x ) = x2 + x4 ,with x (1) = −3 and x (0) = −4. Run each method for 10 iterations. Make twoplots:
1. Plot f vs. the iteration for each method.
2. Plot f ′ vs. x. Overlay the progression of each method, drawing lines from( x (i) , f ′ ( x (i) )) to ( x (i+1) , 0) to ( x (i+1) , f ′ ( x (i+1) )) for each transition.
What can we conclude about this comparison?
Exercise 6.6. Give an example of a sequence of points x (1) , x (2) , . . . and a function
f such that f ( x (1) ) > f ( x (2) ) > · · · and yet the sequence does not converge to a
local minimum. Assume f is bounded from below.
Exercise 6.7. What is the advantage of a Quasi-Newton method over Newton’s
method?
Exercise 6.8. Give an example where the BFGS update does not exist. What would
you do in this case?
Exercise 6.9. Suppose we have a function f (x) = ( x1 + 1)2 + ( x2 + 3)2 + 4. If westart at the origin, what is the resulting point after one step of Newton’s method?
Exercise 6.10. In this problem we will derive the optimization problem fromwhich the Davidon-Fletcher-Powell update is obtained. Start with a quadratic
approximation at x(k) :
⊤



⊤ 
 1
x − x( k ) H( k ) x − x( k )f ( k ) (x ) = y ( k ) + g( k )x − x( k ) +2© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

98 c ha p te r 6 . se c on d -ord e r me thod s

where y(k) , g(k) , and H(k) are the objective function value, the true gradient, anda positive definite Hessian approximation at x(k) .
The next iterate is chosen using line search to obtain: −1
x( k +1) ← x( k ) − α ( k ) H( k )g( k )

We can construct a new quadratic approximation f (k+1) at x(k+1) . The approximation should enforce that the local function evaluation is correct:f ( k +1) ( x( k +1) ) = y ( k +1)
and that the local gradient is correct:∇ f ( k +1) (x( k +1) ) = g( k +1)
and that the previous gradient is correct:∇ f ( k +1) (x( k ) ) = g( k )Show that updating the Hessian approximation to obtain H(k+1) requires:8
H( k +1) δ( k +1) = γ( k +1)
Then, show that in order for H(k+1) to be positive definite, we require:9δ( k +1)⊤γ( k +1) > 0

Finally, assuming that the curvature condition is enforced, explain why one
then solves the following optimization problem to obtain H(k+1) :10
minimize

H − H( k )subject to

H = H⊤

H

Hδ(k+1) = γ(k+1)

8This condition is called the secant
equation. The vectors δ and γ are
defined in equation (6.11).

This condition is called the curvature condition. It can be enforced
using the Wolfe conditions during
line search.

9

The Davidon-Fletcher-Powell update is obtained by solving suchan optimization problem to obtainan analytical solution and thenfinding the corresponding update
equation for the inverse Hessian
approximation.

10

where H − H(k) is a matrix norm that defines a distance between H and H(k) .

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com