Global Optimization by Multilevel Coordinate Search.
Inspired by a method by Jones et al. (1993), we present a global optimization algorithm based on multilevel coordinate search. It is guaranteed to converge if the function is continuous in the neighborhood of a global minimizer. By starting a local search from certain good points, an improved convergence result is obtained. We discuss implementation details and give some numerical results.
Introduction
Problems involving global optimization (traditionally usually minimization) of a multivariate
function are widespread in the mathematical modeling of real world systems
for a broad range of applications (see, e.g., Pint' er [18]). Many problems
can be described only by nonlinear relationships, which introduces the possibility
of multiple local minima. The task of global optimization is to find a solution for
with the objective function obtains its smallest value, the global minimum. When
the objective function has a huge number of local minima, local optimization techniques
are likely to get stuck before the global minimum is reached, and some kind
of global search is needed to find the global minimum with some reliability. Our
Global Optimization Homepage on the World Wide Web (WWW) at the address
http://solon.cma.univie.ac.at/ ~ neum/glopt.html contains a large number of commented
links to information and software packages relevant to global optimization.
Algorithms for solving global minimization problems can be classified into heuristic
methods that find the global minimum only with high probability, and methods that
guarantee to find a global optimum with a required accuracy. An important class
belonging to the former type are the stochastic methods (e.g., Boender & Romeijn
[3]), which involve function evaluations at a suitably chosen random sample of points
and subsequent manipulation of the sample to find good local (and hopefully global)
minima. A number of techniques like simulated annealing (e.g., Ingber [9, 10]) and
genetic algorithms (e.g., Michalewicz [13]) use analogies to physics and biology to
approach the global optimum.
The most important class of methods of the second type are branch and bound meth-
ods. They derive their origin from combinatorial optimization (e.g., Nemhauser &
Wolsey [14]), where also global optima are wanted but the variables are discrete and
take a few values only. Branch and bound methods guarantee to find a global minimizer
with a desired accuracy after a predictable (though often exponential) number
of steps. The basic idea is that the configuration space is split recursively by branching
into smaller and smaller parts. This is not done uniformly but instead some parts
are preferred and others are eliminated. The details depend on bounding procedures.
Lower bounds on the objective allow to eliminate large portions of the configuration
space early in the computation so that only a (usually small) part of the branching
tree has to be generated and processed. The lower bounds may be obtained by using
dc-methods (e.g., Horst & Tuy [8]), techniques of interval analysis (e.g., Hansen
or majorization resp. minorization methods based on the knowledge of Lipschitz
constants (e.g., Pint' er [18]). Unlike heuristic methods, however, these methods are
only applicable if something about the analytical properties of the objective is known,
since one needs to be able to compute powerful and reliable underestimating functions.
The algorithm we are going to describe in this paper is an intermediate between purely
heuristic methods and methods that allow an assessment of the quality of the minimum
obtained; it is in spirit similar to the direct method for global optimization by Jones
et al. [11]. As the latter method, our method is guaranteed to converge in the long run
if the objective is continuous in the neighborhood of a global minimizer; no additional
smoothness properties are required. In contrast to many stochastic methods that
operate only at the global level and are therefore quite slow, our algorithm contains
local enhancements that lead to quick convergence once the global part of the algorithm
has found a point in the basin of convergence of the global optimum. Moreover, for all
control variables in our algorithm meaningful default values can be chosen that work
simultaneously for most problems. Thus fine tuning is usually not required, and since
it is a deterministic method, there is no need for multiple runs.
In this paper, we consider the bound constrained optimization problem
min f(x)
(1)
with finite or infinite bounds, where we use interval notation for rectangular boxes,
ng
for . In the case where all bounds are infinite we obtain an unconstrained
optimization problem.
In direct, a finite box is normalized to [0; 1] n and a tree of boxes is constructed. Each
box is characterized by its midpoint, and the side lengths of the boxes are always of
the form 3 \Gammak , k 2 N 0 . A disadvantage of that algorithm is that infinite box bounds
cannot be handled. Moreover, since the boundary can never be reached, it converges
more slowly than necessary in cases where the minimizer lies at the boundary. For
example, in the case of functions that are monotonous in each variable the optimizer
is at a vertex, but direct converges very slowly, especially when the bounds are wide.
Inspired by direct, we devised a global optimization algorithm based on multilevel
coordinate search (MCS). In our algorithm we remedy the above shortcomings and
allow for a more irregular splitting procedure.
In Section 2, we first define a class of algorithms for solving the global optimization
problem (1) satisfying six properties A1-A6 and prove their convergence in the long
run. In addition to the splitting procedure, an enhancement is obtained by starting
local searches from certain good points, which results in in an improved convergence
result. In Section 3, an outline of our implementation of the MCS algorithm is given,
and details are explained in Sections 4 to 6. Finally, numerical results are presented in
Section 7.

Acknowledgments

The authors gratefully acknowledge partial support of this research by the Austrian
Fond zur F-orderung der wissenschaftlichen Forschung (FWF) under grant P11516-
MAT.
2. Convergence results
In this section we consider a class of algorithms for solving the minimization problem
(1) with six properties A1-A6 defined as follows.
Each box in the tree is characterized by a base point x and its function value f(x).
Moreover, we assign to each box a level s 2 g. When a box is split, its
level is set to zero; the levels of non-split boxes are positive. The level of the initial
box is 1. The base points of boxes with level s max and their function values are
put into the base list.
The levels of MCS correspond to the side lengths of direct, i.e., the boxes with small
level are the 'large' boxes that have not been split very often yet. Only boxes of the
same level compete for being chosen for branching; cf. A3.
A2: A box is always split along a single coordinate i (to be chosen appropriately), the
splitting index, into two or more (at most p - 2) parts, and the descendants of a box
with level latter only when it is - s
The base points of the descendants of a box are chosen in a natural way such that they
differ from the base point of the parent box in (at most) one coordinate and thus this
procedure of generating new function values is a variant of the standard coordinate
search method. In contrast to our algorithm, direct splits a box along all coordinates
in each step.
A3: The branching process proceeds by a series of sweeps through the levels defining
a record. A sweep is defined by the following three steps.
Step 1. Scan the list of non-split boxes and define a record list containing for each
level pointing to a box with the lowest function value among
all boxes at level s. If there is no box at level s, set b Initialize s to the lowest
level with b s 6= 0.
Step 2. The box with label b s is a candidate for splitting. If it is split, mark it as split
in the list and insert its children. Update the record list if any of the children yields
a strict improvement of f on its level. If it is not considered to be worthwhile to split
the box, its level is increased by one and possibly b s+1 has to be updated.
Step 3. Increase s by 1. If new sweep. Else if b to Step 3, else
go to Step 2.
Clearly, each sweeps ends after at most s visits to Step 3.
Like in direct, this strategy yields a combination of global and local search. The key
to balancing global and local search is the multilevel approach. The fact that we start
with the boxes at the lowest levels (i.e., the 'large' boxes) in each sweep constitutes
the global part, and the selection of the box with lowest function value at each level
forms the local part of the algorithm.
A4: For each M ? 0 there exists an r(M) 2 (0; 1) independent of the choice of s max
such that for each box [x that has been generated by splitting the box [x; y] with
along the ith coordinate we have
This assumption ensures that the descendants of a finite box will eventually become
arbitrarily narrow after sufficiently many splits along each coordinate. The worst case
r(M) of the shrinking factor may depend on the box bounds since we do not want to
overemphasize large numbers in the case of wide box bounds.
A5: The rule for splitting boxes with infinite bounds is as follows. If the ith box
bounds are (\Gamma1; 1) and we split along the ith coordinate, the interval is always split
at one or more values contained in a fixed interval [z 1
(independent of the choice of
means of which it turns into two semibounded intervals contained in (\Gamma1; z 2
possibly additional finite ones contained in [z 1
When a box with
the ith bounds [x along the ith coordinate, we obtain a box with ith
bounds one or more boxes with finite ith bounds. Generally let [x
be the unbounded part generated by splitting the ith bounds [x
the ith coordinate. There exist functions
monotonously increasing on (0; 1), such that F 2
(which implies F m
all x 0 2 R . An analogous assumption holds for successive splits of (\Gamma1; x 0 ].
The lower bound F 1 guarantees that the unbounded descendants of the interval [x
shrink sufficiently fast. However, they should not shrink too fast since then the bounded
descendants could become arbitrarily large, and therefore we impose the upper bound
Assume that a box with level s, 1 - s ! s max , is a candidate
for splitting since b s points to this box, and let n j be the number of times coordinate
j has been split in the history of this box. If s ? H(n; minn j ), the box is in any case
eligible for splitting, and the splitting index is a coordinate i with
This assumption makes sure that along each path in the tree, each coordinate is split
inifinitely often.
Proposition 1. In each sweep, one box leaves the lowest non-empty level and no box is
added at that level. Each level will eventually become empty; in particular the splitting
procedure will come to an end when all non-split boxes have level s max .
Proof. Directly from A1-A3.
Proposition 2. Let 2 - s 0 - s max . Then the levels up to s are empty after at
most (p s defined by A2. In particular, the algorithm
finishes after at most (p smax
Proof. In the course of the algorithm, the root box [u; v] with level one can at most have
descendants of level s, 2 - s - s max \Gamma1, since the worst case is splitting a box into p
parts and advancing the levels of the descendants in steps of one. By Proposition 1, the
levels 1 to s are cleared after at most
As expected, this upper bound is exponential. However, in practice, the bound is very
pessimistic since we do not only split a box at the lowest non-empty level (or advance
its level by one) in each sweep but do that at each non-empty level. Moreover, the
levels of the descendants of the boxes will sometimes advance by two, or a box will
be advanced to the next level without being split. Also, good approximations to the
global minimizer will usually have been reached long before all levels are cleared.
Proposition 3. For a box with level
0-l-m
each coordinate has been split at least m times.
Proof. This is proved by induction over m. The statement is trivial for
be a box with level s ?
i be the number of times coordinate i has been split in the history of B 0 . Then
by induction hypothesis. Assume that minn be the parent
box of B its level. Since the levels advance in steps of one
or two and s ? 2(m 1)n, it is possible to go back n times in the history of box B 0 ,
and we have s that
was obtained by splitting B i according to A6. Let n 0
i be the number of times
coordinate i has been split in the history of box B n . Since B 0 was generated from B n
by splitting n times according to A6, we have minn 0
contradiction
to the induction hypothesis. Consequently our assumption was wrong and we must
have
Proposition 4. For each exists an m i (ffi) such that the
ith side length of the box containing a point - x 2 [u; v] with finite coordinates is less
than ffi if it has been split at least m i (ffi) times along the ith coordinate.
Proof. Consider the ith side length of the box containing - x.
Case 1. u i is finite, implies that there exists an l 2 N such that
F l
after at most l splits in the ith coordinate the box containing - x has
finite ith bounds contained in [u
and the problem is reduced to Case 2. The
case treated similarly. The case reduced to
a semibounded interval contained in [z 1
or a finite interval contained
in [z 1
after one split.
Case 2. u i and v i finite, M := max(ju i j; jv i j). Then there exists a k 2 N 0 such that
defined in A4, i.e., by A4 the ith side length of the
containing - x is less than ffi if the box has been split at least k times along the ith
coordinate.
Propositions 2-4 yield the following convergence theorem.
Theorem 5. Suppose that the global minimization problem (1) has a finite solution -
R is continuous in a neighborhood of -
x. Let the algorithm satisfy the
assumptions A1-A6, and let " ? 0. Then there exists an s 0 such that for s
after at most (p s a base point x is found with
i.e., the algorithm converges if the number of levels tends to 1.
In order to enhance the performance of the algorithm, we subject a base point of a
box of level s max to local search before putting it into the base list, which is another
local feature of our algorithm. For theoretical purposes we make here the obviously
idealized assumption that the local search algorithm reaches a local minimizer after
finitely many steps if it is started in its basin of attraction.
Theorem 6. Let the assumptions of Theorem 5 be satisfied, assume that a local search
is started from each candidate for the base list, at least when its function value is smaller
than the function value of any local minimum already in the base list, and assume that
there is an " ? 0 such that f(y) ? nonglobal local minimizer y and
for any y 2 [u; v] with sufficiently large norm. Then there exist numbers L and S such
that, for any s max - L, a global minimizer is found after at most S sweeps.
Proof. By assumption there exists a
with
Propositions 3 and 4 there exists an L such that the
side lengths of the box containing -
x are smaller than ffi when its level is at least L, and by
Proposition 2 the box containing -
x reaches level L after at most S := (p
sweeps. Let s max - L. Then after at most S sweeps the box B containing -
x is at level
L, and for its base point x we have f(x) - f(-x) Two cases are possible. If
we have found a better point before, we have already made a local search from such a
point, which finds a global minimizer by assumption. In the case that x is the current
best point, at the end of that sweep a point y with f(y) - f(x) is subjected to local
search, by means of which we find a global minimizer.
3. The MCS algorithm
There are many ways to design algorithms satisfying A1-A6 and hence guaranteeing
convergence to a global minimizer. However, trivial implementations are very slow.
A number of heuristic enhancements are needed to obtain a high quality method. In
the following we describe the most useful version that we could find. We first give an
overview over the techniques used and discuss the details of each technique in Sections
4 to 6.
Unlike in direct, the base point of a box in our algorithm is usually not the midpoint
but a point at the boundary, often but not always a vertex, and a base point can belong
to more than one box. Moreover, we also assign to each box an opposite point. The
construction is such that the base point x and the opposite point y determine the box,
and we call such a box B[x; y].
The algorithm starts with a so-called initialization procedure producing an initial tree
of boxes. Whenever a box is split along some coordinate i for the first time (either
in the initialization procedure or later), this is done at three or more given values x l
(where function values are computed) and some adaptively chosen intermediate points,
and at least four subboxes are obtained.
After the initialization procedure we start with the sweeps as described in A3. When a
box with level 1 - s ! s max has been selected for splitting by the procedure described
in A3, we use the following splitting strategy: the box is
ffl either split according to A6 (splitting by rank),
ffl or it is split along a coordinate where the maximal gain in function value is
expected according to a local separable quadratic model obtained by fitting 2n+1
function values (splitting by expected gain),
ffl or, if the expected gain is not large enough, it is not split at all but its level is
increased by one.
If the box has already been split along the coordinate chosen as splitting index, only
a single new function evaluation is needed, and we obtain two or three subboxes. The
new function value is taken at a point that is predetermined by the box bounds in the
case of splitting by rank, and in the case of splitting by expected gain, at a point where
the maximal gain in function value is expected.
Before putting a base point of a box of level s max into the base list, we first check
whether it is in the the basin of attraction of a local minimum already in the base list,
and if this is not the case, we subject it to local search. The local search algorithm
used in our program essentially consists of building a local quadratic model by triple
searches, defining a promising search direction with the aid of the quadratic model
and a modified Cholesky factorization and making a line search along this direction.
This procedure is repeated until a stopping criterion is fulfilled. Moreover, devices to
handle the case that one or more coordinates of the current point are at the boundary
are included.
4. Initialization
Before building the tree of boxes, one first evaluates f at an initial point x 0 and sets
evaluated at two or more points in [u; v] that agree
with x   in all coordinates k 6= i. Thus one has
x l
f l
and x  . The point with smallest function value is then renamed x
before repeating the procedure with the next coordinate. The numbers x l
and the indices l i are stored in an initialization list. The choice of the initialization list
is left to the user, who may incorporate in it knowledge about the likely distribution of
good points. A good known starting point can be defined as x 0 . Some possible choices
are discussed in Section 7.
From the initialization list and the corresponding list of function values, an initial tree
of boxes is constructed as follows. The root box is B[x;
point and as y one of the corners of [u; v] farthest away from x. Note that x need not
be a vertex and that some or all coordinates of y can be infinite. For
current box is split along the ith coordinate into 2L
exactly one of the x l
i as endpoints, depending on whether two, one or none of the x l
are
on the boundary, which means that in addition to x l
we have to split at
z l
. The additional splitting points are chosen as z l
is the golden section ratio and
chosen such that the part with the smaller function value gets the larger fraction of
the interval. The resulting subboxes get as base point the point x 0 obtained from the
current by changing x i to the x l
i that is a boundary point of the corresponding
ith coordinate interval, so that f(x 0
i , and as opposite point the point obtained
from y by changing y i to the other end of that interval.
The information available so far allows us to define priorities on the coordinates. For
each i, we compute the union of the ranges of the quadratic interpolant through any
three consecutive
take the difference of the upper and lower bound obtained
as a crude measure of the variability of f with the ith component. Components with
higher variability get a higher priority, and this ranking is saved in a vector - such that
the component with index i has the - i th highest estimated variability. Moreover, if the
x   obtained after splitting the ith coordinate belongs to two boxes, the one containing
the minimizer of the quadratic models is taken as current box for coordinate i + 1.
The root box gets level 1. If a box of level s is split, the boxes with the smaller
fraction of the golden section split get level s
Thus the current box for splitting in the next coordinate is in any case one with level
after finishing the initialization procedure, the first level is empty and the
non-split boxes have levels which implies that it is meaningful to take
Two examples for the set of boxes, their base points and their levels after
the initialization procedure in the two-dimensional case are shown in Figure 1 (where
6
in both cases x
It is easy to see the connection between the golden section split and the assignment of
levels. When a box B with level s is split along coordinate i according to the golden
section split and its larger part B 0 is again subjected to the golden section split along
the ith coordinate, the larger descendant of B 0 has the same ith coordinate length as
the smaller descendant of B, and these boxes both have level s 2. Moreover, the box
with the better function value gets the larger fraction of the interval and the smaller
r
r
Figure
r
r

Figure
level because then it is more likely to be split again more quickly, which was also the
strategy adopted in direct.
Any choice of x l
including the endpoints u i , v i in the list guarantees that, in the simple
case that f is monotonous in each variable, the final x   of the initialization phase is
already the global minimizer.
5. Splitting
We do not store the box bounds for each box but the information in the tree (label of
the parent box, splitting index, splitting value, a label identifying which of the many
children it is etc. This keeps the amount of storage proportional to the number of
function evaluations and allows us to recover information to build a separable quadratic
model by going back in the history of the box.
Suppose that we want to split the ith coordinate interval u t fx
for is the ith component of the base point of the box being
considered. In order to fulfil assuption A4, we may not split too close to x i . If y i is
large, we also do not want the new component x 0
i to be too large and therefore force it
to be from some smaller interval u t f- g. We choose this interval according to
where the function subint is given by the following pseudo-Matlab function:
function
else
Having selected a box with level s ! s max for splitting (by the procedure described
in A3), we recover its base point x, the opposite point y and the number n i of times
coordinate i has been split in the history of the box.
Splitting by rank. Let s ? H(n; minn i ), where H is the function defined by A6.
Then we select the splitting index i among the indices i with smallest n i as the one
with lowest - i (and hence highest variability rank).
the splitting is done according to the initialization list at x l
and at the golden section split points, as discussed in Section 4, and the new base
points and opposite points are defined as before. The boxes with the smaller fraction
of the golden section split (and thus larger function values) get level min(s
and all other ones get level s + 1.
the ith component ranges between x i and y i , and the splitting value is
chosen as z The box is split at z i and
at the golden section split point, and we obtain three parts with only one additional
function evaluation at the point x 0 obtained by changing the ith coordinate of x to
z i . The smaller fraction of the golden section split gets level min(s and the
two other parts get level s + 1. Moreover, the base point of the first child is taken to
be x, the base point of the second and third child is the point x 0 defined above, and
the opposite points are obtained by changing y i to the other end of the ith coordinate
interval of the corresponding box.
Splitting by expected gain. Let s - H(n; minn i ). In order to build a local separable
quadratic model, we need two additional points and corresponding function values for
each coordinate. Whenever we have split in the ith coordinate in the history of the
box, we obtain values that can be used for quadratic interpolation in this coordinate.
For each coordinate we take the first two points and function values found by pursuing
the history of the box back to [u; v] since these points are expected to be closest to the
base point x. For coordinates that have not yet been split, we obtain this information
from the initialization list. Let
be the local separable model for f(-) generated by interpolation at x and the 2n
additional points collected as above. For each coordinate i, we define the expected
gain - e i in function value when we evaluate at a new point obtained by changing this
coordinate in the base point. Again two cases have to be distinguished.
Case 1. In the history of the current box, coordinate i was never split, i.e.,
Then we split according to the initialization list at points where we already know the
obtainable function differences, and therefore compute the expected gain as
Case 2. If n i ? 0, the ith component ranges between x i and y i , and with the quadratic
partial correction function
at our disposal, we can calculate the maximal gain expected when changing the value of
only. For the reasons discussed above, we choose the splitting value from u t f-
where
Then we compute
with minimum achieved at - . If the expected best function value satisfies
1-i-n
best ;
where f best is the current best function value (including the function values obtained
by local optimization), we expect the box to contain a better point and split, using
as splitting index the component with minimal -
e i . The condition (2) prevents wasting
function evaluations by splitting boxes with bad base point function values; these boxes
will eventually be split by rank anyway.
In Case 1 we again split according to the initialization list, and the definition of the
new base points and opposite points and the assignment of levels are as before. In Case
2 we use z i as splitting value and the box is split at z i (if z i 6= y i ) and at the golden
section split point, and we obtain two or three parts. The larger fraction of the golden
section split gets level s + 1, the smaller fraction level min(s
the third part is larger than the smaller fraction of the golden section split, it gets level
Moreover, the base point of the first child
is taken to be x, the base point of the second and third (if z i 6= y i ) child is obtained
by changing the ith coordinate of x to z i , and the opposite points are again obtained
by changing y i to the other end of the ith coordinate interval of the box.
If (2) is violated, we do not expect any improvement and therefore do not split but
increase the level by 1.
6. Local search
The theory of local optimization provides powerful tools for the task of optimizing a
smooth function when knowledge of the gradient or even the Hessian is assumed. When
no derivative information is available, the traditional methods are based on the employment
of conjugate directions and successive line searches; cf. the direction set method
of Powell [19] and a modification due to Brent [4]. These algorithms, however, do
not allow the specification of bound constraints. Elster & Neumaier [6] developed
an algorithm for optimization of low-dimensional bound constrained functions, based
on the use of quadratic models and a restriction of the evaluation points to successively
refined grids. Hoewever, the work in that algorithm grows with the dimension n as
O(n 6 ) and hence is unsuitable for larger dimensions.
The local optimization algorithm we are going to describe in the sequel also makes use of
quadratic models and successive line searches, and devices to handle bound constraints
are incorporated. We first explain the procedure of building a local quadratic model
by triple searches. Since we do not want to start local searches from points belonging
to the domain of attraction of the same local minimum and to put several copies of
essentially the same point into the base list, we have devised a criterion to avoid that;
see Step 2 in the subsection 'base list' below. We conclude with a discussion of the
validity of assumptions A4 and A5 for our algorithm.
Triple search. We want to use
function values to construct a quadratic model
best
best best
Assume that we have three vectors x l (with componentwise inequalities).
The function values are to be taken at points x with x
as follows. If
x best denotes current best point in the triple search, denote by x (i;1) and x (i;2) the points
obtained from x best by changing its ith coordinate to the other two values in fx l
and by x ik , k ! i, the points obtained by changing the ith and kth coordinate to the
ones with the smaller q(x (i) ) resp. q(x (k) ) with the current quadratic model q. Thus
we obtain the following procedure that we are going to describe in more detail in the
sequel:
best )
compute f(x (i;1) ) and f(x (i;2) ); compute g i and G ii
store x newbest but do not update x best
compute q(x (k;1) ) and q(x (k;2) ) from the current model
compute f(x ik )
update x newbest but do not update x best
compute G ik
if x newbest 6= x best , update x best , f and g 1:i
When, for an computed approximations for g l and G lk ,
by interpolating at
points that differ only in the first
best is the current best point in the triple search, we obtain
approximations for g i and G ii by determining these numbers such that the quadratic
polynomial
best best
best
interpolates at (x (i;j)
2. At this stage we do not yet update x best , but
if min(f(x (i;1) ); f(x (i;2) best ), we store this point as x bestnew .
Assume that, in addition, we have already calculated approximations for G il , 1 - l -
1. Then we can compute
best
best
best
with the current quadratic model q, and we have q(x (i;j) 2. Let
x ik be defined as above. Then we choose G ki such that the quadratic model
interpolates at x ik , i.e., such that the equation
best
best
best
best
best
best
best
is satisfied. Again we do not update x best but update x bestnew if x ik yields a strict
improvement in function value.
Only after finishing the loop over k, we reexpand the model around the best point by
G kl (x bestnew
best
l
It is easy to see that the above method gives the unique quadratic interpolant to
f at
distinct points; in particular, we recover the exact objective function as
In a diagonal triple search we only carry out the diagonal part of the above algorithm
and take the off-diagonal elements of the Hessian from the previous iteration; thus only
2n additional function values are needed.
Coordinate search. To find x l use a coordinate search based on
a line search routine. A Matlab version of the actual line search used can be obtained
electronically from http://solon.cma.univie.ac.at/ ~ neum/ms/ls/ This univariate
line search program ls0 contains a parameter smaxls limiting the number of points
used for the line search. It is possible to feed other points in addition to the starting
point and their function values into the program, and these points are included in
smaxls.
A line search with along each coordinate. The first line search
is started with the candidate for the base list. After the line search along the first
coordinate, we can take the best point and its two nearest neighbors on both sides
(or, if such points do not exist, the two nearest neighbors on one side) as fx l
g.
The subsequent line searches are started with the current best point obtained from
the previous line search. After a line search in a coordinate i ? 1, we take the best
point, the starting point of the line search (if it is different from the best point) and, if
possible, the nearest neighbor of the best point on the other side as fx l
g. The
ith coordinate of the old x best has to be among fx l
since otherwise we would
lose all points through which the surface has been fitted previously.
Local search. Now we have all ingredients at our disposal to describe the steps of the
local search algorithm used in our implementation of MCS.
Step 1. Starting with the candidate for the base list, we make a full triple search,
where fx l are found by a coordinate search as described above. This procedure
yields a new point x, its function value f , an approximation g of the gradient and an
approximation of the Hessian G.
Step 2. Set is the modified Cholesky factorization
without scaling of [15, Section 2] and the starting value " used. Set
of the current x is u k or v k and p k points outward,
set we make a line search with ls0 along x+ffp. Let ff min and ff max be the
bounds for ff obtained from the box bounds [u; v]. If p T Gp ? 0 and \Gammag T p=(p T Gp) 2
\Gammag T p=(p T Gp); otherwise we take -
as the one
where we expect the smaller function value according to the quadratic model. The
values
ff are used as input for ls0. If G is positive definite and "
is chosen, the modified Cholesky factorization is the Cholesky factorization, the search
direction is g, and we have -
positive definite
ff is close to 1. Set old is the
function value of the current point at the beginning of Step 2 and f pred is the function
value predicted by the above quadratic model with
ff, i.e., fac is a measure of
the quality of the quadratic model, and it is 1 if the line search did not improve the
function value.
Step 3. Stop if some limit on the number of visits to Step 3 (per local search) or the
limit on function calls has been exceeded. Also stop if none of the components of the
current x are at the boundary and the last triple search was a full one and the stopping
criterion is fulfilled.
Step 4. If some components of the current x are at the boundary and the stopping
criterion is fulfilled, we make line searches with at most smaxls points along these
coordinates. If the function value was not improved by these coordinate searches, stop.
Step 5. If fac ? 0:1 or the stopping criterion was fulfilled in Step 3, we make a
full triple search, otherwise we make a diagonal triple search. We make these triple
searches only in the coordinates i such that the component x i of the current x is not
at the boundary. The set fx l
is taken to consist of x i ,
two neighbors at distance ffi and 2ffi if x i lies at the boundary)
and
", where " denotes the machine accuracy. We obtain a new point x and an
approximation of its reduced gradient g and its reduced Hessian G.
Step 6. We set only in the components of the
current x that are not at the boundary (p make a line search
along x
ff as input as before. The quantity fac is defined
as in Step 2, where f old is the function value at the current point at the beginning of
Step 6. and go to Step 3.
The stopping criterion is fulfilled if the function value was not improved in Steps 5
and 6 (resp. Steps 1 and 2) or if jgj T max(jxj; jx old is an input
parameter of our program and f 0 is a 'typical' function value, e.g., the smallest function
value found in the initialization procedure.
Base list. The local searches are only carried out at the end of each sweep. Then
all candidates for the base list that have been collected in this sweep are ordered by
ascending function value, and they are gone through as follows. Let x be a candidate
for the base list.
Step 1. We check whether we have already made a local search from this point. (Often
a point belongs to two boxes.) If this is the case, we take the next candidate for the
base list and go to Step 1.
Step 2. Let w
be the points already in the base list, and by renumbering
assume that they are sorted by their distance to x, starting with the nearest point. For
we do the following.
Step 2a. Compute the function value at x increase
i by one and go to Step 2a.
Step 2b. Compute the function value at x
increase i by one and go to Step 2a. Else if
points seem to lie in the same valley. However,
we do not discard x for local search but set x to the value x 0 or x 00 with the smaller
function value, increase i by one and go to Step 2a. Else we do not subject x to local
search because the four function values are monotonous, take the next candidate for
the base list and go to Step 1.
Step 3. Local search.
Step 4. The new point x obtained from local search is subjected to a procedure similar
to Step 2 in order to find out whether we have really found a new point, and only in
that case it is put into the base list.
If we disregard the splits according to the initialization list, A4 is fulfilled for our
algorithm with 0:001. The factor
0.9 applies to the 'regular' case where for [x 00 y. For
the 'safeguarded' case of subint for large y the worst case is at most the largest of
the three numbers 1(
and
1). Moreover, A5 is satisfied for the choice
7. Numerical results
Test functions. Jones et al. [11] gave an extensive comparison of their direct
method with various methods on seven standard test functions from Dixon & Szeg-
[5] and two test functions from Yao [24]. Since our MCS algorithm is based on important
insights from [11], we first consider the same test set to evaluate the efficiency of
MCS. For each test function, the dimensions and box bounds, used by Jones et al. but
inadvertently omitted in [11], are given in Table 1. We thank Don Jones for providing
us with the code of the test functions.
Label Test function Dimension Default box bounds
GP Goldstein-Price 2 [\Gamma2; 2] 2
C6 Six-hump camel 2 [\Gamma3; 3] \Theta [\Gamma2; 2]

Table

1. Dixon & Szeg- dimensions and box bounds
All but the last four lines of Table 2 are taken from one of the tables of results of
[11]. The fourth last line contains results for the differential evolution algorithm DE
by Storn & Price [22] from http://http.icsi.berkeley.edu/ ~ storn/code.html,
using the Matlab program devec2.m with the default values for the control parame-
ters. Since this algorithm operates only at the global level, it takes a rather long time
to find a minimum with high accuracy and therefore we used obtaining a relative error
stopping criterion. The number of function evaluations needed for convergence
was averaged over 25 runs for each test function. In the case of Hartman6, one
run did not converge after 12 000 function evaluations and we averaged only over the
remaining 24 runs. The last three lines give results for MCS, first over the same bound
constraints as in [11], then averages over strongly perturbed box bounds, and finally
for the unconstrained version. Details will be given below.
When assessing the results, we have to bear in mind that the first 11 algorithms already
appeared in the 1978 anthology edited by Dixon & Szeg- are therefore are
somewhat old. The number of function calls needed for convergence is not the only
method of assessing the quality of an algorithm, but it is an important one in the case
of most real life applications, where function evaluations are expensive.
Termination. In the presentation of test results, methods are usually compared on
the basis of their performance on problems with known solutions. The algorithm is
terminated when a function value within some tolerance of the global minimum has
been found, and we also adopt this strategy. However, in practical problems, one does
not know the solution in advance and needs a criterion that tells the program when to
stop searching for a better local minimizer. This criterion should be stringent enough
that it does not waste too many function values after the global minimum has been
found, but it should also be loose enough to ensure that in typical cases, the algorithm
does not terminate before the global minimizer has been found.
Stochastic approaches to the design of suitable stopping criteria are surveyed in Section
6 of Boender & Romeijn [3]. One of the methods proposed there consists in
stopping when the number n of local searches done is larger than a function N(w)
of the number w of different local minima found so far. The function N(w) depends
on the assumptions, and several specific implicit definitions of N(w) are given in [3].
Method S5 S7 S10 H3 H6 GP BR C6 SHU
Bremmerman [5] (a) (a) (a) (a) (a) (a) 250
Mod. Bremmerman [5] (a) (a) (a) (a) 515 300 160
Zilinskas [5] (a) (a) (a) 8641 5129
Gomulka-Branin [5] 5500 5020 4860
T-orn [5] 3679 3606 3874 2584 3447 2499 1558
Gomulka-T-orn [5] 6654 6084 6144
Gomulka-V.M. [5] 7085 6684 7352 6766 11125 1495 1318
Price [5] 3800 4900 4400 2400 7600 2500 1800
Mockus [5] 1174 1279 1209 513 1232 362 189
Belisle et al. [1] (b) 339 302 4728 1846
Boender et al. [2] 567 624 755 235 462 398 235
Kostrowicki-Piela [12] (g) (g) (g) 200 200 120 120
Perttunen [16] (c) 516 371 250 264 82 97 54 197
Perttunen-
Stuckman [17] (c) 109 109 109 140 175 113 109 96 (a)
Jones et al. [11] (c) 155 145 145 199 571 191 195 285 2967
Storn-Price [22] (d)(h)6400 6194 6251 476 7220 1018 1190 416 1371
All but the last four lines are taken from [11]; missing entries were not available
from the literature.
(a) Method converged to a local minimum.
(b) Average evaluations when converges. For H6, converged only 70 % of time.
(c) Convergence defined as obtaining a relative error ! 0:01 %.
(d) Convergence defined as obtaining a relative error
Perturbed box bounds.
(f) Unconstrained optimization problem.
(g) Global minimum not found with less than 12 000 function calls.
Average over 25 cases. For H6, average over 24 cases
one case did not converge within 12 000 function values.
(i) An asterisk indicates that the first local optimization gave the local optimum.

Table

2. Number of function evaluations for various methods
compared to MCS
This result is theoretically justified for the random multiple start method only but may
serve as a guideline also for other methods that use local searches.
However, with MCS we try to do very few local optimizations only, and this reasoning
appears inadequate. So far, we have not yet found a useful general purpose stopping
criterion for MCS. For the purposes of the numerical tests reported here, the stopping
criterion for MCS was taken as obtaining a relative error ! 0:01 % in the optimal
objective function value (which happens to be nonzero always), i.e.,
\Gamma4 , which was also the criterion used by Jones et al. [11]. For the algorithms quoted in
[11], results based on the definition of convergence used by their authors are reported.
MCS control parameter settings. We applied a Matlab version of MCS with
is the dimension of the problem, H(n; m)
to the test functions and used a simple initialization list consisting of
midpoint and boundary points, i.e.,
2:
The limit on visits to Step 3 per local search was set to 50, and the parameter fl in the
stopping criterion for local optimization was taken as (cf. the subsection 'local
search' of Section 6). Note that all examples have been run with identical parameter
settings, so that no tuning to the individual test problems was involved.
Modified bounds. We also investigated the stability of our results for MCS with
respect to random perturbations of the box bounds. Instead of the default box bounds
[u; v] given in Table 1, we employed the box bounds [u
where j is a random variable that is uniformly distributed in the interval (\Gamma0:5; 0:5),
but a value of j was only accepted in a given problem if at least one of the global
minimizers was in [u The results given in the second last line of Table 2 were
taken as an average over 25 runs with different perturbed box bounds for each test
function. For Hartman6, we obtained one outlier for which the algorithm had not
found the global minimum after 12 000 function calls, and we report a result averaged
over the 24 remaining runs.
Moreover, we applied MCS to the unconstrained optimization problem for the Dixon
test set and added the results to Table 2. In this case, we cannot use an
initialization list consisting of midpoint and boundary points any more. For infinite
or very wide bounds, we suggest to use the following safeguarded version based on the
function subint and again take l 2:
if
else
Discussion. The results show that MCS seems to be strongly competitive with existing
algorithms in the case of problems with reasonable finite bound constrained. MCS with
unperturbed box bounds wins in 7 of 9 test cases against every competing algorithm
and is only beaten once by Perttunen resp. Perttunen-Stuckman for the remaining
two test functions. MCS with perturbed box bounds still wins against all competing
algorithms for 5 test functions. Only the results for Shekel's functions seem to depend
heavily on the choice of the box bounds, but they are comparable with the results of
some other algorithms, and we do not know whether direct is stable with respect to
perturbation of the box bounds.
For unconstrained problems of dimension n - 4, the performance of MCS is less
satisfactory. The reason is that in the exploration of an unbounded domain, it is easy
to miss the region where the global minimum lies if one has already found a low-lying
nonglobal minimizer. For example, for Shekel's functions the algorithm gets caught in
the second best local minimizer (1; 1; 1; 1), which is a point in the initialization list.
Further test problems. The Dixon & Szeg- test set has been criticized for containing
mainly easier test problems. A more challenging test set was used in the first
contest on evolutionary optimization (ICEO) at the ICEC'96 conference; cf. Storn
Price [23]. This test bed contains five problems, each in a 5-dimensional and a
10-dimensional version, and on these test functions, MCS showed some limitations.
The names and default box bounds of the ICEO test functions are given in Table 3,
and the results are shown in Table 4. The first two lines in Table 4 are results, taken
from [23], of two different versions of DE.
Problem Name Box bounds
1 Sphere model [\Gamma5; 5] n
3 Shekel's foxholes [0; 10] n
4 Michalewicz's function [0; -] n
5 Langerman's function [0; 10] n

Table

3. ICEO test functions and their box bounds
We applied MCS with three different choices of the initialization list to the ICEO test
functions. MCS1 is the standard version with midpoints and boundary points. For
MCS2, we took x 1
i.e., the points are uniformly spaced but do not include the boundary points.
For MCS3, we generated an initialization list with the aid of line searches. Starting with
the absolutely smallest point in [u; v], we made line searches with ls0 with
and along each coordinate, where the best point was taken as starting point
for the next line search. The parameter nloc in ls0 determines how local or global the
line search is since the algorithm tries to find up to nloc minima within smaxls function
values. For the line searches in the local search method described in Section 6,
was taken (entirely local line search). For each coordinate i, all local minimizers found
by the line searches were put into the initialization list, and if their number was less
than three, they were supplemented with the values obtained from ls0 closest to u i
and v i .
Finally, we applied MCS to the test functions used by Storn & Price [22]. Their
names and box bounds are shown in Table 5, and their definition can be found in [22].
Problem 7 is a shifted version of ICEO2 for Problems 8 and 9 have general
constraints and were therefore not used here.
The first four lines of Table 6 are taken from [22]. ANM denotes the annealed Nelder
Mead strategy of [20] and ASA the Adaptive Simulated Annealing method by Ingber
[9, 10], and DE1 and DE3 are two different versions of DE.
For MCS, we used an initialization list consisting of midpoint and boundary points
when there was no global minimizer among the points in the initialization list; in the
latter case, a different initialization list with Problem
4 contains a random variable, and the result presented for MCS was averaged over
Problem 1 Problem 2 Problem 3 Problem 4 Problem 5
5D 10D 5D 10D 5D 10D 5D 10D 5D 10D
MCS3 26 51 26956 - 32077 - 1903 -
values.
A dash indicates that a global minimizer was not found after 100 000 function
calls.

Table

4. Number of function values for the ICEO functions
Problem Name Dimension n Box bounds
4 Quartic with random noise
8 Zimmerman's problem 2 constraints
9 Polynomial fit, 9, 17 constraints

Table

5. Test functions of Storn & Price [22], their dimensions
and box bounds
ANM 95 106 90258 -
ASA 397 11275 354 4812 1379 3581 -
discontinuous test function; 2 noisy test function.
A dash indicates that a global minimizer was not
found after 100 000 function calls.

Table

6. Number of function values for the Storn & Price functions
runs, where convergence was defined as reaching a point with function value - 15.
Problem 1 has a quadratic objective function, hence is easy for MCS, and Problem 3
is easy for MCS since the objective function is monotonous.
8. Conclusions
The multilevel coordinate search algorithm MCS, presented in this paper, has excellent
theoretical convergence properties if the function is continuous in the neighborhood of
a global minimizer. In the current implementation, our test results show that MCS
is strongly competitive with existing algorithms in the case of problems with reason-able
finite bound constrained. In our comparison, MCS outperforms the competing
algorithms almost always on the classical test problems set of Dixon & Szeg- o.
For unconstrained problems of dimension n - 4, the performance of MCS is less
since in the exploration of an unbounded domain, it is easy to miss the
region where the global minimum lies if one has already found a low-lying nonglobal
minimizer. The same problem applies for some hard test problems with a huge number
of local minima.



--R


Rinnoy Kan

Algorithms for Minimization without Derivatives
The global optimization problem: an introduction
A grid algorithm for bound constrained optimization of noisy func- tions

Deterministic Approaches
Very fast simulated re-annealing

Lipschitzian optimization without the Lipschitz constant
Diffusion equation method of global minimization: performance on standard test functions
Genetic Algorithms
Integer and Combinatorial Optimization
On satisfying second-order optimality conditions using modified Cholesky factoriza- tions

The rank transformation applied to a multiunivariate method of global optimization

An efficient method for finding the minimum of a function of several variables without calculating derivatives
Numerical Recipes in C
A multistart global minimization algorithm with dynamic search trajectories
Differential evolution - a simple and efficient adaptive scheme for global optimization over continuous spaces
Minimizing the real functions of the ICEC'96 contest by differential evolution
Dynamic tunneling algorithm for global optimization
--TR

--CTR
Rong Yan , Alexander G. Hauptmann, The combination limit in multimedia retrieval, Proceedings of the eleventh ACM international conference on Multimedia, November 02-08, 2003, Berkeley, CA, USA
H.-M. Gutmann, A Radial Basis Function Method for Global Optimization, Journal of Global Optimization, v.19 n.3, p.201-227, March 2001
evket lker Birbil , Shu-Cherng Fang , Ruey-Lin Sheu, On the Convergence of a Population-Based Global Optimization Algorithm, Journal of Global Optimization, v.30 n.2-3, p.301-318, November  2004
J. M. Gablonsky , C. T. Kelley, A Locally-Biased form of the DIRECT Algorithm, Journal of Global Optimization, v.21 n.1, p.27-37, September 2001
Y. Wu , L. Ozdamar , A. Kumar, TRIOPT: a triangulation-based partitioning algorithm for global optimization, Journal of Computational and Applied Mathematics, v.177 n.1, p.35-53, 1 May 2005
U. M. Garcia-Palomares , F. J. Gonzalez-Castao , J. C. Burguillo-Rial, A Combined Global & Local Search (CGLS) Approach to Global Optimization, Journal of Global Optimization, v.34 n.3, p.409-426, March     2006
Jaewook Lee, A novel three-phase trajectory informed search methodology for global optimization, Journal of Global Optimization, v.38 n.1, p.61-77, May       2007
A. Ismael Vaz , Lus N. Vicente, A particle swarm pattern search method for bound constrained global optimization, Journal of Global Optimization, v.39 n.2, p.197-219, October   2007
Panayiotis G. Georgiou , Chris Kyriakakis, Maximum likelihood parameter estimation under impulsive conditions, a sub-Gaussian signal approach, Signal Processing, v.86 n.10, p.3061-3075, October 2006
Michael Bartholomew-Biggs , Bruce Christianson , Ming Zuo, Optimizing Preventive Maintenance Models, Computational Optimization and Applications, v.35 n.2, p.261-279, October   2006
local optimization;coordinate search;global optimization;bound constraints