Direct methods rely solely on the objective function f . These methods are also
called zero-order, black box, pattern search, or derivative-free methods. Direct methods
do not rely on derivative information to guide them toward a local minimum
or identify when they have reached a local minimum. They use other criteria tochoose the next search direction and to judge when they have converged.

7.1 Cyclic Coordinate Searchx2

7 Direct Methods

x1

(1)

(2)

(1)

(1)

(7.1)

(2)

(7.2)

x1 = arg min f ( x1 , x2 , x3 , . . . , xn )x1

x2

Cyclic coordinate search, also known as coordinate descent or taxicab search, simply
alternates between coordinate directions for its line search. The search starts froman initial x(1) and optimizes the first input:

Figure 7.1. Cyclic coordinate descent alternates between coordinate directions.

Having solved this, it optimizes the next coordinate:
(3)

(2)

(2)

x2 = arg min f ( x1 , x2 , x3 , . . . , xn )x2This process is equivalent to doing a sequence of line searches along the set of nbasis vectors, where the ith basis vector is all zero except for the ith component,
which has value 1 (algorithm 7.1). For example, the third basis function, denotede(3) , in a four-dimensional space is:e(3) = [0, 0, 1, 0]

(7.3)

x1
Figure 7.2. Above is an example of
how cyclic coordinate search can
get stuck. Moving in either of thecoordinate directions will result
only in increasing f , but moving
diagonally, which is not allowed incyclic coordinate search, can resultin lowering f .

c ha p te r 7. d i re c t me thod sbasis(i, n) = [k == i ? 1.0 : 0.0 for k in 1 : n]Figure 7.1 shows an example of a search through a two-dimensional space.
Like steepest descent, cyclic coordinate search is guaranteed either to improveor to remain the same with each iteration. No significant improvement after a fullcycle over all coordinates indicates that the method has converged. Algorithm 7.2
provides an implementation. As figure 7.2 shows, cyclic coordinate search canfail to find even a local minimum.
function cyclic_coordinate_descent(f, x, ϵ)Δ, n = Inf, length(x)while abs(Δ) > ϵx′ = copy(x)for i in 1 : n
d = basis(i, n)x = line_search(f, x, d)end
Δ = norm(x - x′)
end
return x
end

The method can be augmented with an acceleration step to help traverse
diagonal valleys. For every full cycle starting with optimizing x(1) along e(1)and ending with x(n+1) after optimizing along e(n) , an additional line search is
conducted along the direction x(n+1) − x(1) . An implementation is provided inalgorithm 7.3 and an example search trajectory is shown in figure 7.3.

7.2 Powell’s Method

Algorithm 7.1. A function for constructing the ith basis vector oflength n.

Algorithm 7.2. The cyclic coordinate descent method takes as input the objective function f and astarting point x, and it runs until
the step size over a full cycle is lessthan a given tolerance ϵ.

x2

100

x1
original
accelerated
Figure 7.3. Adding the acceleration step to cyclic coordinate descent helps traverse valleys. Six
steps are shown for both the original and accelerated versions.

Powell’s method was first introduced by M. J. D. Powell, “An Efficient Method for Finding the
Powell’s method1 can search in directions that are not orthogonal to each other. The Minimum of a Function of Sevmethod can automatically adjust for long, narrow valleys that might otherwise eral Variables Without Calculating Derivatives,” Computer Journal,
require a large number of iterations for cyclic coordinate descent or other methods vol. 7, no. 2, pp. 155–162, 1964.
An overview is presented by W. H.
that search in axis-aligned directions.
Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical
Recipes in C: The Art of Scientific
Computing. Cambridge University
Press, 1982, vol. 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
1

7. 2. p ow e l l ’s me thod 101

The algorithm maintains a list of search directions u(1) , . . . , u(n) , which are
initially the coordinate basis vectors, u(i) = e(i) for all i. Starting at x(1) , Powell’smethod conducts a line search for each search direction in succession, updating
the design point each time:x(i+1) ← line_search( f , x(i) , u(i) ) for all i in {1, . . . , n}Algorithm 7.3. The cyclic coordinate descent method with an acceleration step takes as input theobjective function f and a starting
point x, and it runs until the step
size over a full cycle is less than agiven tolerance ϵ.

x2function cyclic_coordinate_descent_with_acceleration_step(f, x, ϵ)Δ, n = Inf, length(x)while abs(Δ) > ϵx′ = copy(x)for i in 1 : n
d = basis(i, n)x = line_search(f, x, d)end
x = line_search(f, x, x - x′) # acceleration stepΔ = norm(x - x′)
end
return x
end

(7.4)

Next, all search directions are shifted down by one index, dropping the oldest
search direction, u(1) :u(i) ← u(i+1) for all i in {1, . . . , n − 1}

(7.5)

u( n ) ← x( n +1) − x(1)

(7.6)

x1
Figure 7.4. Powell’s method starts
the same as cyclic coordinate descent but iteratively learns conjugate directions.

The last search direction is replaced with the direction from x(1) to x(n+1) , whichis the overall direction of progress over the last cycle:

and another line search is conducted along the new direction to obtain a new x(1) .
This process is repeated until convergence. Algorithm 7.4 provides an implementation. Figure 7.4 shows an example search trajectory.
Powell showed that for quadratic functions, after k full iterations the last k
directions will be mutually conjugate. Recall that n line searches along mutuallyconjugate directions will optimize an n-dimensional quadratic function. Thus, nfull iterations of Powell’s method, totaling n(n + 1) line searches, will minimize aquadratic function.
© 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

102

c ha p te r 7. d i re c t me thod sfunction powell(f, x, ϵ)n = length(x)U = [basis(i,n) for i in 1 : n]
Δ = Infwhile Δ > ϵx′ = xfor i in 1 : n
d = U[i]x′ = line_search(f, x′, d)end
for i in 1 : n-1
U[i] = U[i+1]end
U[n] = d = x′ - x
x′ = line_search(f, x′, d)Δ = norm(x′ - x)
x = x′
end
return x
end

Algorithm 7.4. Powell’s method,which takes the objective functionf, a starting point x, and a toleranceϵ.

The procedure of dropping the oldest search direction in favor of the overall
direction of progress can lead the search directions to become linearly dependent.
Without search vectors that are linearly independent, the search directions canno longer cover the full design space, and the method may not be able to find theminimum. This weakness can be mitigated by periodically resetting the search
directions to the basis vectors. One recommendation is to reset every n or n + 1
iterations.

7.3 Hooke-Jeeves
The Hooke-Jeeves method (algorithm 7.5) traverses the search space based on evaluations at small steps in each coordinate direction.2 At every iteration, the HookeJeeves method evaluates f (x) and f (x ± αe(i) ) for a given step size α in every
coordinate direction from an anchoring point x. It accepts any improvement it may
find. If no improvements are found, it will decrease the step size. The processrepeats until the step size is sufficiently small. Figure 7.5 shows a few iterationsof the algorithm.

R. Hooke and T. A. Jeeves, “Direct
Search Solution of Numerical andStatistical Problems,” Journal of theACM (JACM), vol. 8, no. 2, pp. 212–
229, 1961.
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

x1

x1

x2

x2

x2

x2

7. 4 . g e n e ra l i z e d p a tte rn se a rc h 103x1x1

One step of the Hooke-Jeeves method requires 2n function evaluations foran n-dimensional problem, which can be expensive for problems with many
dimensions. The Hooke-Jeeves method is susceptible to local minima. The method
has been proven to converge on certain classes of functions.3

7.4 Generalized Pattern Search

Figure 7.5. The Hooke-Jeeves
method, proceeding left to right.
It begins with a large step size
but then reduces it once it cannot
improve by taking a step in any
coordinate direction.
3
E. D. Dolan, R. M. Lewis, and V.
Torczon, “On the Local Convergence of Pattern Search,” SIAM
Journal on Optimization, vol. 14,
no. 2, pp. 567–583, 2003.

In contrast with the Hooke-Jeeves method, which searches in the coordinate
directions, generalized pattern search can search in arbitrary directions. A pattern P
can be constructed from a set of directions D about an anchoring point x with astep size α according to:

P = {x + αd for each d in D}(7.7)

The Hooke-Jeeves method uses 2n directions for problems in n dimensions, but
generalized pattern search can use as few as n + 1.
For generalized pattern search to converge to a local minimum, certain conditions must be met. The set of directions must be a positive spanning set, whichmeans that we can construct any point in R n using a nonnegative linear combination of the directions in D . A positive spanning set ensures that at least one of thedirections is a descent direction from a location with a nonzero gradient.4
We can determine whether a given set of directions D = {d(1) , d(2) , . . . , d(m) }in R n is a positive spanning set. First, we construct the matrix D whose columns
are the directions in D (see figure 7.6). The set of directions D is a positive spanning set if D has full row rank and if Dx = −D1 with x ≥ 0 has a solution.5 This
optimization problem is identical to the initialization phase of a linear program,which is covered in chapter 11.

Convergence guarantees for generalized pattern search require thatall sampled points fall on a scaled
lattice. Each direction must thus bea product d( j) = Gz( j) for a fixednonsingular n × n matrix G andinteger vector z. V. Torczon, “On
the Convergence of Pattern Search
Algorithms,” SIAM Journal of Optimization, vol. 7, no. 1, pp. 1–25,
1997.
4

R. G. Regis, “On the Properties of
Positive Spanning Sets and Positive
Bases,” Optimization and Engineering, vol. 17, no. 1, pp. 229–262, 2016.
5© 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

104 c ha p te r 7. d i re c t me thod sfunction hooke_jeeves(f, x, α, ϵ, γ=0.5)
y, n = f(x), length(x)while α > ϵimproved = false
x_best, y_best = x, y
for i in 1 : n
for sgn in (-1,1)x′ = x + sgn*α*basis(i, n)y′ = f(x′)
if y′ < y_best
x_best, y_best, improved = x′, y′, true
end
end
end
x, y = x_best, y_best

Algorithm 7.5. The Hooke-Jeeves
method, which takes the targetfunction f, a starting point x, a starting step size α, a tolerance ϵ, anda step decay γ. The method runs
until the step size is less than ϵ andthe points sampled along the coordinate directions do not provide an
improvement. Based on the implementation from A. F. Kaupe Jr, “Algorithm 178: Direct Search,” Communications of the ACM, vol. 6, no. 6,
pp. 313–314, 1963.

if !improved
α *= γ
end

end

end
return x

Figure 7.6. A valid pattern for generalized pattern search requires apositive spanning set. These directions are stored in the set D .

only positively spans theconeonly positively spans 1d
space

positively spans R2© 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

7. 5 . n e l d e r-me a d si mp l e x me thod 105

Figure 7.7. All previous points ingeneralized pattern search lie on ascaled lattice, or mesh. The lattice is
not explicitly constructed and neednot be axis-aligned.

x (1)x (2) = x (3)x (4) = x (5)The implementation of generalized pattern search in algorithm 7.6 contains
additional enhancements over the original Hooke-Jeeves method.6 First, the implementation is opportunistic—as soon as an evaluation improves the current best
design, it is accepted as the anchoring design point for the next iteration. Second,
the implementation uses dynamic ordering to accelerate convergence—a directionthat leads to an improvement is promoted to the beginning of the list of directions.
Figure 7.7 shows a few iterations of the algorithm.

These enhancements are presented in C. Audet and J. E. Dennis
Jr., “Mesh Adaptive Direct Search
Algorithms for Constrained Optimization,” SIAM Journal on Optimization, vol. 17, no. 1, pp. 188–217,
2006.
6

7.5 Nelder-Mead Simplex Method
The Nelder-Mead simplex method7 uses a simplex to traverse the space in searchof a minimum. A simplex is a generalization of a tetrahedron to n-dimensional
space. A simplex in one dimension is a line, and in two dimensions it is a triangle(see figure 7.8). The simplex derives its name from the fact that it is the simplest
possible polytope in any given space.

7
The original simplex method
is covered in J. A. Nelder and
R. Mead, “A Simplex Method
for Function Minimization,” The
Computer Journal, vol. 7, no. 4,
pp. 308–313, 1965. We incorporate the improvement in J. C. Lagarias, J. A. Reeds, M. H. Wright,
and P. E. Wright, “Convergence
Properties of the Nelder–Mead
Simplex Method in Low Dimensions,” SIAM Journal on Optimization, vol. 9, no. 1, pp. 112–147, 1998.

© 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

106 c ha p te r 7. d i re c t me thod sfunction generalized_pattern_search(f, x, α, D, ϵ, γ=0.5)
y, n = f(x), length(x)while α > ϵimproved = falsefor (i,d) in enumerate(D)x′ = x + α*dy′ = f(x′)
if y′ < yx, y, improved = x′, y′, true
D = pushfirst!(deleteat!(D, i), d)
break
end
end
if !improved
α *= γend
end
return x
end

The Nelder-Mead method uses a series of rules that dictate how the simplex is
updated based on evaluations of the objective function at its vertices. A flowchart
outlines the procedure in figure 7.9, and algorithm 7.7 provides an implementation. Like the Hooke-Jeeves method, the simplex can move around while roughly
maintaining its size, and it can shrink down as it approaches an optimum.
The simplex consist of the points x(1) , . . . , x(n+1) . Let xh be the vertex with thehighest function value, let xs be the vertex with the second highest function value,and let xℓ be the vertex with the lowest function value. Let x̄ be the mean of allvertices except the highest point xh . Finally, for any design point xθ , let yθ = f (xθ )be its objective function value. A single iteration then evaluates four simplex
operations:

Algorithm 7.6. Generalized pattern search, which takes the target function f, a starting point x,
a starting step size α, a set of search
directions D, a tolerance ϵ, and astep decay γ. The method runs until the step size is less than ϵ andthe points sampled along the coordinate directions do not provide an
improvement.

Figure 7.8. A simplex in two dimensions is a triangle. In order forthe simplex to be valid, it musthave a nonzero area.

Reflection. xr = x̄ + α(x̄ − xh ), reflects the highest-valued point over the centroid.
This typically moves the simplex from high regions toward lower regions. Here,α > 0 and is typically set to 1.

© 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

7. 5 . n e l d e r-me a d si mp l e x me thod

107

Initial simplex
Sort the simplex entries
Compute x̄
Compute the reflection point:
xr = x̄ + α(x̄ − xh )yr < y l ?

no

yes

yr ≥ y s ?

yes

no

yr ≥ y h ?

no

yes

Replace xh with xr

Compute the expansion point:
xe = x̄ + β(xr − x̄)Compute the contraction point:
xc = x̄ + γ(xh − x̄)y e < yr ?

yc > yh ?

no

no

yes
Replace xh with xe

Shrink by replacing allx(i) with (x(i) + xℓ )/2

Replace xh with xc

Replace xh with xr

Converged?

yes

yes

Return best point

noFigure 7.9. Flowchart for the
Nelder-Mead algorithm.

© 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

108 c ha p te r 7. d i re c t me thod s

Expansion. xe = x̄ + β(xr − x̄), like reflection, but the reflected point is sent evenfurther. This is done when the reflected point has an objective function valueless than all points in the simplex. Here, β > max(1, α) and is typically set to 2.
Contraction. xc = x̄ + γ(xh − x̄), the simplex is shrunk down by moving awayfrom the worst point. It is parameterized by γ ∈ (0, 1) which is typically set to
0.5.
Shrinkage. All points are moved toward the best point, typically halving theseparation distance.
Figure 7.10 shows the four simplex operations. Figure 7.11 shows several iterationsof the algorithm.
xs

xs

xs
xe

xr
xh

x̄
xℓ
reflection

xh

x̄
xℓ
expansion

xr

xh

xs

xc

xh

x̄
xℓ
contraction

The convergence criterion for the Nelder-Mead simplex method is unlike
Powell’s method in that it considers the variation in the function values rather
than the changes to the points in the design space. It compares the standard
deviation of the sample8 y(1) , . . . , y(n+1) to a tolerance ǫ. This value is high for a
simplex over a highly curved region, and it is low for a simplex over a flat region.
A highly curved region indicates that there is still further optimization possible.

7.6 Divided Rectangles
The divided rectangles algorithm, or DIRECT for DIvided RECTangles, is a Lipschitzian optimization approach, similar in some ways to the Shubert-Piyavskii
method described in section 3.6.9 However, it eliminates the need for specifying a
Lipschitz constant, and it can be more efficiently extended to multiple dimensions.

xℓshrinkage
Figure 7.10. The Nelder-Mead simplex operations visualized in twodimensions.

The standard deviation of the sample
is also called the uncorrected sample
standard
deviation. In our case, it isq2
n +1 ( i )1− ȳ , where ȳ isn +1 ∑ i =1 y
8

the mean of y(1) , . . . , y(n+1) .

D. R. Jones, C. D. Perttunen, andB. E. Stuckman, “Lipschitzian Optimization Without the Lipschitz
Constant,” Journal of Optimization
Theory and Application, vol. 79, no. 1,
pp. 157–181, 1993.
9© 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

7. 6 . d i vi d e d re c ta n g l e s 109

Figure 7.11. The Nelder-Meadmethod, proceeding left to rightand top to bottom.

© 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

110 c ha p te r 7. d i re c t me thod sfunction nelder_mead(f, S, ϵ; α=1.0, β=2.0, γ=0.5)
Δ, y_arr = Inf, f.(S)while Δ > ϵp = sortperm(y_arr) # sort lowest to highest
S, y_arr = S[p], y_arr[p]xl, yl = S[1], y_arr[1] # lowest
xh, yh = S[end], y_arr[end] # highest
xs, ys = S[end-1], y_arr[end-1] # second-highest
xm = mean(S[1:end-1]) # centroid
xr = xm + α*(xm - xh) # reflection point
yr = f(xr)Algorithm 7.7. The Nelder-Meadsimplex method, which takes theobjective function f, a starting simplex S consisting of a list of vectors, and a tolerance ϵ. The NelderMead parameters can be specified as well and default to recommended values.

if yr < ylxe = xm + β*(xr-xm) # expansion pointye = f(xe)S[end],y_arr[end] = ye < yr ? (xe, ye) : (xr, yr)
elseif yr ≥ ys
if yr < yh
xh, yh, S[end], y_arr[end] = xr, yr, xr, yr
end
xc = xm + γ*(xh - xm) # contraction point
yc = f(xc)
if yc > yh
for i in 2 : length(y_arr)
S[i] = (S[i] + xl)/2
y_arr[i] = f(S[i])endelse
S[end], y_arr[end] = xc, yc
end
elseS[end], y_arr[end] = xr, yr
end

end

Δ = std(y_arr, corrected=false)
end
return S[argmin(y_arr)]© 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

7. 6 . d i vi d e d re c ta n g l e s 111

The notion of Lipschitz continuity can be extended to multiple dimensions. If
f is Lipschitz continuous over a domain X with Lipschitz constant ℓ > 0, thenfor a given design x(1) and y = f (x(1) ), the circular cone
(7.8)

f (x(1) ) − ℓkx − x(1) k2

forms a lower bound of f . Given m function evaluations with design points{x(1) , . . . , x(m) }, we can construct a superposition of these lower bounds by taking
their maximum:
maximize f (x(i) ) − ℓkx − x(i) k2(7.9)
i

divided rectangles lower bound

Lipschitz lower bound
2

x2

1

0

Figure 7.12. The Lipschitz lower
bound is an intersection of cones,
which creates complicated surfaces
in multidimensional space. The divided rectangle lower bound isolates each lower-bound cone toits own hyper-rectangular region,
making it trivial to compute theminimum value in each region
given a Lipschitz constant.

−1
−2
−2

−1

0
x1

1

2

−2

−1

0
x1

1

2

The Shubert-Piyavskii method samples at a point where the bound derived
from a known Lipschitz constant is lowest. Unfortunately, a Lipschitz lower boundhas intricate geometry whose complexity increases with the dimensionality ofthe design space. The left contour plot in figure 7.12 shows such a lower boundusing five function evaluations. The right contour plot shows the approximation
made by DIRECT, which divides the region into hyper-rectangles—one centeredabout each design point. Making this assumption allows for the rapid calculationof the minimum of the lower bound.

© 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

112 c ha p te r 7. d i re c t me thod s

The DIRECT method does not assume a Lipschitz constant is known. Figures 7.13 and 7.14 show lower bounds constructed using Lipschitz continuity
and the DIRECT approximation, respectively, for several different Lipschitz constants. Notice that the location of the minimum changes as the Lipschitz constantchanges, with small values of ℓ leading to designs near the lowest function evaluations and larger values of ℓ leading to designs the farthest from previous function
evaluations.

y

4

3

−4

−2

0
x

2

4

y

4

3

−4

−2

0
x

2

4

Figure 7.13. The Lipschitz lowerf (x)bound for different Lipschitz conℓ = 0.1
stants ℓ. Not only does the estiℓ = 0.2
mated minimum change locally asℓ = 0.3
the Lipschitz constant is varied, theℓ = 0.4
region in which the minimum lies
ℓ = 0.5
can vary as well.
ℓ = 0.6ℓ = 0.7ℓ = 0.8ℓ = 0.9ℓ = 1.0
lower-bound minima
evaluation points

Figure 7.14. The DIRECT lowerf (x)bound for different Lipschitz conℓ = 0.1
stants ℓ. The lower bound is notℓ = 0.2
continuous. The minimum doesℓ = 0.3
not change locally but can changeℓ = 0.4
regionally as the Lipschitz constantℓ = 0.5changes.
ℓ = 0.6ℓ = 0.7ℓ = 0.8ℓ = 0.9ℓ = 1.0
lower-bound minima
evaluation points© 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

7. 6 . d i vi d e d re c ta n g l e s 113

7.6.1 Univariate DIRECT
In one dimension, DIRECT recursively divides intervals into thirds and thensamples the objective function at the center of the intervals as shown in figure 7.15.
This scheme is in contrast with the Shubert-Piyavskii method, where sampling
occurs at a point where the bound derived from a known Lipschitz constant islowest.
For an interval [ a, b] with center c = ( a + b)/2, the lower bound based on f (c)is:f ( x ) ≥ f (c) − ℓ| x − c|(7.10)

before division:
after division:

Figure 7.15. Center-point sampling,using the DIRECT scheme, divides
intervals into thirds.

where ℓ is a Lipschitz constant that is unknown to us. The lowest value that thebound obtains on that interval is f (c) − ℓ(b − a)/2, and it occurs at the edges of
the interval.
Even though we do not know ℓ, we can deduce that the lower bound of someintervals are lower than others. If we have two intervals of the same length, forexample, and the first one has a lower evaluation at the center point than thesecond one, then the lower bound of the first interval is lower than that of thesecond. Although this does not entail that the first interval contains a minimizer,
it is an indication that we may want to focus our search in that interval.
During our search, we will have many intervals [ a1 , b1 ], . . . , [ an , bn ] of different
widths. We can plot our intervals according to their center value and interval
width, as we do in figure 7.16. The lower bound for each interval is the verticalintercept of a line of slope ℓ through its center point. The center of the intervalwith the lowest lower bound will be the first point intersected as we shift a linewith slope ℓ upwards from below.
DIRECT splits all intervals for which a Lipschitz constant exists such thatthey have the lowest lower bound, as shown in figure 7.17. We refer to theseselected intervals as potentially optimal. Any interval may technically containthe optimum, though the selected points heuristically have the best chance ofcontaining the optimum.
One iteration of the one-dimensional DIRECT algorithm consists of identifying
the set of potentially optimal intervals and then dividing each interval into thirds.
Example 7.1 demonstrates the univariate DIRECT method.

© 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

c ha p te r 7. d i re c t me thod sinterval center value f (c)114(bi − ai )/2

slope = ℓFigure 7.16. Interval selection for
a particular Lipschitz constant ℓ.
Black dots represent DIRECT intervals and the corresponding function evaluations at their centers.
A black line of slope ℓ is drawn
through the dot belonging to the selected interval. All other dots mustlie on or above this line.

f ( ci )f (ci ) − ℓ(bi − ai )/2

interval center value f (c)interval half-width (b − a)/2

Figure 7.17. The potentially optimal intervals for the DIRECT
method form a piecewise boundary that encloses all intervals along
the lower-right. Each dot corresponds to an interval.

suboptimalpotentially optimal

interval half-width (b − a)/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

7. 6 . d i vi d e d re c ta n g l e s 115

Consider the function f ( x ) = sin( x ) + sin(2x ) + sin(4x ) + sin(8x ) on theinterval [−2, 2] with a global minimizer near −0.272. Optimization is difficult
because of multiple local minima.
The figure below shows the progression of the univariate DIRECT method,with intervals chosen for splitting rendered in blue. The left side shows theintervals overlaid on the objective function. The right side shows the intervals
as scatter points in the interval half-width versus center value space.

Example 7.1. The DIRECT method
applied to a univariate function.

y

2
0

−2

y

2
0

−2

y

2
0

−2

y

2
0

−2

y

2
0

−2

y

2
0

−2
−2

−1

0
x

1

20

0.5

1

1.5

2

interval half-width© 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

116

c ha p te r 7. d i re c t me thod s

7.6.2 Multivariate DIRECT
In multiple dimensions, we divide rectangles (or hyper-rectangles in more than
two dimensions) instead of intervals. Similar to the univariate case, we divide therectangles into thirds along the axis directions. Before commencing the divisionof the rectangles, DIRECT normalizes the search space to be the unit hypercube.
As illustrated in figure 7.18, the choice of ordering of the directions whensplitting the unit hypercube matters. DIRECT prioritizes the assignment of larger
rectangles for the points with lower function evaluations. Larger rectangles areprioritized for additional splitting.
horizontal split

vertical split

vertical split

horizontal split

sampling

When splitting a region without equal side lengths, only the longest dimensions
are split (figure 7.19). Splitting then proceeds on these dimensions in the samemanner as with a hypercube.
The set of potentially optimal intervals is obtained as it is with one dimension.
The lower bound for each hyper-rectangle can be computed based on the longest
center-to-vertex distance and center value. We can construct a diagram similar tofigure 7.16 to identify the potentially optimal rectangles.10

7.6.3 Implementation
DIRECT (algorithm 7.8) is best understood when it is broken down into subroutines. We present these subroutines below.
Normalization to the unit hypercube is done by algorithm 7.9.

Figure 7.18. Interval splitting in
multiple dimensions (using DIRECT) requires choosing an ordering for the split dimensions.

Figure 7.19. DIRECT will split only
the longest dimensions of hyperrectangles.
As an additional requirement for
selection, DIRECT also requiresthat the lower bound for the interval improve the current best value
by a nontrivial amount.
10© 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

7. 6 . d i vi d e d re c ta n g l e s 117

function direct(f, a, b, ϵ, k_max)
g = reparameterize_to_unit_hypercube(f, a, b)intervals = Intervals()n = length(a)c = fill(0.5, n)interval = Interval(c, g(c), fill(0, n))add_interval!(intervals, interval)c_best, y_best = copy(interval.c), interval.y

Algorithm 7.8. DIRECT, whichtakes the multidimensional objective function f, vector of lower
bounds a, vector of upper bounds
b, tolerance parameter ϵ, and number of iterations k_max. It returns
the best coordinate.

for k in 1 : k_max
S = get_opt_intervals(intervals, ϵ, y_best)
to_add = Interval[]
for interval in S
append!(to_add, divide(g, interval))dequeue!(intervals[vertex_dist(interval)])end
for interval in to_add
add_interval!(intervals, interval)
if interval.y < y_best
c_best, y_best = copy(interval.c), interval.y
end
end
end
end

return rev_unit_hypercube_parameterization(c_best, a, b)

rev_unit_hypercube_parameterization(x, a, b) = x.*(b-a) + afunction reparameterize_to_unit_hypercube(f, a, b)
Δ = b-areturn x->f(x.*Δ + a)end

Algorithm 7.9. A function that creates a function defined over theunit hypercube that is a reparameterized version of the function fdefined over the hypercube withlower and upper bounds a and b.

© 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

118

c ha p te r 7. d i re c t me thod s

Computing the set of potentially optimal rectangles can be done efficiently(algorithm 7.10). We use the fact that the interval width can only take on powersof one-third, and thus many of the points will share the same x coordinate. For
any given x coordinate, only the one with the lowest y value can be a potentially
optimal interval. We store the rectangular intervals according to their depth andthen in a priority queue according to their centerpoint value.
using DataStructures
struct Intervalcydepthsend
min_depth(interval) = minimum(interval.depths)
vertex_dist(interval) = norm(0.5*3.0.^(-interval.depths), 2)
const Intervals =OrderedDict{Float64,PriorityQueue{Interval, Float64}}function add_interval!(intervals, interval)
d = vertex_dist(interval)
if !haskey(intervals, d)intervals[d] = PriorityQueue{Interval, Float64}()
end
return enqueue!(intervals[d], interval, interval.y)end

Algorithm 7.10. The data structureused in DIRECT. Here, Interval
has three fields: the interval centerc, the center point value y = f(c),and the number of divisions ineach dimension depths. The function add_interval! inserts a new
Interval into the data structure.

We can use this data structure to obtain all potentially optimal intervals (algorithm 7.11). The algorithm proceeds from lowest interval width to highest interval
width. For each point, we first determine whether it is above or below the line
joining the previous two points. If it is below, we skip it. The same determinationis then made for the next point.
Finally, we need a method for dividing the intervals. This is implemented by
the divide method (algorithm 7.12).
The intervals obtained from running DIRECT in two dimensions are visualized
in figure 7.20. Two iterations of DIRECT in two dimensions are worked out inexample 7.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

7. 6 . d i vi d e d re c ta n g l e s 119

function get_opt_intervals(intervals, ϵ, y_best)stack = Interval[]for (x, pq) in intervals
if !isempty(pq)interval = DataStructures.peek(pq)[1]y = interval.y

Algorithm 7.11. A routine for obtaining the potentially optimal intervals, where intervals is of typeIntervals, ϵ is a tolerance parameter, and y_best is the best function
evaluation.

while length(stack) > 1
interval1 = stack[end]interval2 = stack[end-1]
x1, y1 = vertex_dist(interval1), interval1.y
x2, y2 = vertex_dist(interval2), interval2.yℓ = (y2 - y) / (x2 - x)
if y1 <= ℓ*(x1-x) + y + ϵ
break
end
# remove previous intervalpop!(stack)end
if !isempty(stack) && interval.y > stack[end].y + ϵ# skip new interval
continue
end
end

end

push!(stack, interval) # add new interval

end
return stack© 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

120c ha p te r 7. d i re c t me thod sfunction divide(f, interval)c, d, n = interval.c, min_depth(interval), length(interval.c)dirs = findall(interval.depths .== d)cs = [(c + 3.0^(-d-1)*basis(i,n),c - 3.0^(-d-1)*basis(i,n)) for i in dirs]
vs = [(f(C[1]), f(C[2])) for C in cs]
minvals = [min(V[1], V[2]) for V in vs]

end

7.7

Algorithm 7.12. The divide routine for dividing an interval, where
f is the objective function andinterval is the interval to be divided. It returns a list of the resulting smaller intervals.

intervals = Interval[]
depths = copy(interval.depths)
for j in sortperm(minvals)
depths[dirs[j]] += 1
C, V = cs[j], vs[j]push!(intervals, Interval(C[1], V[1], copy(depths)))
push!(intervals, Interval(C[2], V[2], copy(depths)))
end
push!(intervals, Interval(c, interval.y, copy(depths)))
return intervals

Summary• Direct methods rely solely on the objective function and do not use derivative
information.
• Cyclic coordinate search optimizes one coordinate direction at a time.
• Powell’s method adapts the set of search directions based on the direction ofprogress.
• Hooke-Jeeves searches in each coordinate direction from the current point
using a step size that is adapted over time.
• Generalized pattern search is similar to Hooke-Jeeves, but it uses fewer search
directions that positively span the design space.
• The Nelder-Mead simplex method uses a simplex to search the design space,
adaptively expanding and contracting the size of the simplex in response toevaluations of the objective function.
• The divided rectangles algorithm extends the Shubert-Piyavskii approach tomultiple dimensions and does not require specifying a valid Lipschitz constant.
© 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

7. 7. su mma ry

1

300

250

x2

200

121

Figure 7.20. The DIRECT method
after 16 iterations on the Braninfunction, appendix B.3. Each cell isbordered by white lines. The cells
are much denser around the minima of the Branin function, as the
DIRECT method procedurally increases its resolution in those regions.

150

100

50

0

0

1
x1

© 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

122

c ha p te r 7. d i re c t me thod s

Consider using DIRECT to optimize the flower function (appendix B.4)
over x1 ∈ [−1, 3], x2 ∈ [−2, 1]. The function is first normalized to the unit
hypercube such that we optimize x1′ , x2′ ∈ [0, 1]:

Example 7.2. The first two iterations of DIRECT worked out in detail.

f ( x1′ , x2′ ) = flower(4x1′ − 1, 3x2′ − 2)
The objective function is sampled at [0.5, 0.5] to obtain 0.158. We have asingle interval with center [0.5, 0.5] and side lengths [1, 1]. The interval is
divided twice, first into thirds in x1′ and then the center interval is divided
into thirds in x2′ .
1
x2′
0

0

x1′

1

We now have five intervals:
centerside lengths

vertex distancecenter value

1
2
3
4
5

[1/6, 3/6]
[5/6, 3/6]
[3/6, 3/6]
[3/6, 1/6]
[3/6, 5/6]

[1/3, 1]
[1/3, 1]
[1/3, 1/3]
[1/3, 1/3]
[1/3, 1/3]

0.527
0.527
0.236
0.236
0.236

0.500
1.231
0.158
2.029
1.861

interval center valueinterval

4
5

2

2
1

1

suboptimalpotentially optimal

3

0
0

1/6

1/2

interval half-width

We next split on the two intervals centered at [1/6, 3/6] and [3/6, 3/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

7. 8 . e x e rc i se s 123

7.8 Exercises
Exercise 7.1. Previous chapters covered methods that leverage the derivative
to descend toward a minimum. Direct methods are able to use only zero-orderinformation—evaluations of f . How many evaluations are needed to approximate
the derivative and the Hessian of an n-dimensional objective function using finite
difference methods? Why do you think it is important to have zero-order methods?
Exercise 7.2. Design an objective function and a starting point x0 such that HookeJeeves will fail to attain a reduction in the objective function. You need to choose
x0 such that it is not a local minimizer.
Exercise 7.3. Is the design point obtained using the Hooke-Jeeves method guaranteed to be within ǫ of a local minimum?
Exercise 7.4. Give an example of a concrete engineering problem where you may
not be able to compute analytical derivatives.
Exercise 7.5. State a difference between the divided rectangles algorithm in one
dimension and the Shubert-Piyavskii method.
Exercise 7.6. Suppose our search algorithm has us transition from x(k) = [1, 2, 3, 4]
to x(k+1) = [2, 2, 2, 4]. Could our search algorithm be (a) cyclic coordinate search,(b) Powell’s method, (c) both a and b, or (d) neither a nor b? Why?

© 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