THE CLUSTER PROBLEM IN

MULTIVARIATE GLOBAL OPTIMIZATION

Kaisheng Du and R. Baker Kearfott

University of Southwestern Louisiana

Abstract. We consider branch and bound methods for enclosing all unconstrained global minimizers of a nonconvex nonlinear twice-continuously differentiable objective function. In particular, we consider bounds obtained with interval arithmetic, with the “midpoint test,” but no acceleration procedures. Unless the lower bound is exact, the algorithm without acceleration procedures in general gives an undesirable cluster of boxes around each minimizer. In a previous paper, we analyzed this problem for univariate objective functions. In this paper, we generalize that analysis to multidimensional objective functions. As in the univariate case, the results show that the problem is highly related to the behavior of the objective function near the global minimizers and to the order of the corresponding interval extension.

1. Introduction and Basic Concepts

Our underlying problem is:

(1) find all global minimizers to f(x) subject to x ∈ X,

where X ⊂ ℝ^{m} is a compact right parallelepiped with faces parallel to the axes. We will refer to X as a box. We denote the global minimum as f^{∗} and the set of global minimizers as X^{∗}. Interval branch and bound procedures for unconstrained nonconvex optimization, i.e. for rigorously enclosing the solution set X^{∗} of (1), are competitive with stochastic methods, such as Monte Carlo methods, and methods involving heuristics, such as simulated annealing or the tunneling method. Our analysis deals with algorithms similar to Algorithm 3, p. 111 of [10]. Also, as in [10], we will use interval arithmetic to obtain the bounds.

This paper deals with the phenomenon of clusters of small boxes around global minimizers that such algorithms produce cannot eliminate. We refer to this phenomenon as the cluster problem. Though the algorithms in which the cluster problem occurs apply to global, nonconvex, unconstrained optimization, the cluster problem is essentially a local phenomenon dealing with the portion of the domain

Typeset by AMS-TEX

Key words and phrases. branch and bound principle, inclusion function, interval extensions, midpoint test, global optimization, order of an interval extension, nonconvex optimization.

This work was partially funded by National Science Foundation grant # CCR-9203730

around minimizers in which f is convex. However, the overall algorithm’s efficiency depends on how it handles these local regions; an understanding of this local behavior aids in prediction of efficiency for the overall algorithm.

We introduced and analyzed the phenomenon for the univariate case in [6]. Here, we generalize those results to the multidimensional case. The algorithm presented here is simplified, to facilitate a clear theoretical analysis. Practical algorithms will have more features to improve efficiency. However, their running will generate additional phenomena, such as possible, but not certain, elimination of boxes through convergence to critical points, which, together, are more difficult to analyze.

For introductions to interval arithmetic, see e.g. [1], [7], or [9].

Throughout, we denote non-interval quantities (both points in multidimensional space and scalars) by lower case, interval quantities by upper case boldface, and vectors by lower case marked letters such as p⃗. We will occasionally use a lower case bold letter to denote a non-interval function value that has been bounded using interval arithmetic.

Definition 1. Let X = X_{1} × X_{2} × ⋯ × X_{m}, be an m-dimensional box, where each X_{i} an interval. Then the width of X is

w(X) = max_{1≤i≤m}{w(X_{i})}.

This differs from the componentwise definition, or diameter, used by some authors in some contexts. However, with this definition of width, the order of an inclusion function is formally identical to the corresponding order in the onedimensional case.

Definition 2. Let F(X) denote an interval extension of f evaluated over a box X. Let f(X) denote the exact range of f over X. If there is a constant K, independent of the box X such that

(2)  w(F(X)) − w(f(X)) ≤ Kw(X)^{α},

for all boxes X with w(X) sufficiently small, then we say that F is an order α inclusion function for f. When α is 1 or 2, we call the inclusion first order or second order, respectively.

The following algorithm is Algorithm 3, p. 111 of [10]. We use the following notation.

(1) mid(X): the “midpoint” of a box X, defined componentwise;

(2) ub(X): the upper bound of an interval X;

(3) lb(X): the lower bound of an interval X.

Algorithm 1.

0. Input the original box X, the inclusion function F of f, and additional parameters used in the termination criteria.

1. Set Y := X.

2. Calculate F(X) and f̃:= ub(F(c)) where c := mid(Y).

3. Set y := lb(F(Y)).

4. Initialize the list L := {(Y, y)}.

5. Choose a coordinate direction i_{0} parallel to which Y has an edge of maximum length, i.e., i_{0} ∈ {i : w(Y_{i}) = w(Y)}.

6. Bisect Y normal to direction i_{0} to obtain boxes V_{1} and V_{2} such that Y = V_{1} ∪ V_{2}.

7. Calculate F(V_{1}) and F(V_{2}).

8. Set v_{i} := lb(F(V_{i})) for i = 1, 2.

9. Enter the pairs (V_{1}, v_{1}) and (V_{2}, v_{2}) at the end of the list.

10. Choose a pair (Ỹ, ỹ) of the list which satisfies ỹ ≤ z for all pairs (Z, z) of the list.

11. Discard all pairs {Z, z} from the list that satisfy z > f̃ (midpoint test).

12. If the termination criteria hold go to 15.

13. Denote the first pair of the list by (Y, y). Then set c := midY and f̃ := min(f̃, ub(F(c))).

14. Go to 5.

15. End.

Generally, we will terminate the algorithm if all boxes in the list have widths at most a pre-specified tolerance ϵ.

Step 11 is crucial for our purposes, since we are studying the power of the midpoint test to discard boxes that do not contain global minimizers.

2. The analysis

Our multi-dimensional analysis appears in this section. It is analogous to the one-dimensional case in [6]. In one dimension, we formulated when the algorithm would reject an interval I in terms of the number of intervals generated by the algorithm between an optimizer and I. This number was compared with the excess width of F(I) over the exact range f(I) of f on I. That analysis generalizes to multivariate objective functions. However, we must count the number of boxes in a particular direction between an optimizer and the box in question.

Terminology, general assumptions and a sufficient condition for a box to be rejected.

Assumption. Suppose that x^{∗} ∈ X^{∗} is a particular minimizer that is a strict local minimizer in the sense that there is a neighborhood N of x^{∗} such that f(x) > f(x^{∗}) for x ≠ x^{∗} and x ∈ N.

Our analysis is a local analysis of the behavior of the algorithm in such a neighborhood of x^{∗}. However, the overall algorithm in a region containing many minimizers and sub-regions where f is not convex behaves according to our analysis. This is because nonconvex regions correspond to large function values (local maxima), and are rapidly rejected in Step 11. (We illustrate this with Example 2 in §3 below.)

For simplicity in the analysis, we assume a uniform subdivision of mesh ϵ each coordinate direction in a neighborhood of x^{∗}, as illustrated in Figure 1. This does not seem to diminish the predictive power of our analytical results: The subdivision may have more boxes in one direction from the minimum than in others; however, this non-uniformity does not qualitatively change the clustering behavior. The actual algorithm does produce a uniform subdivision on certain steps.

For small ϵ and near a selected local minimizer x^{∗}, we will represent a typical point x = x_{δ} in one of the boxes X^{(1)} corresponding to the list L in Algorithm 1 as

x_{δ} = x^{∗} + (nϵ + δ)p⃗,

where n is the number of boxes^{1} through which a ray from x^{∗} to the box X^{(1)} passes, p⃗ is a vector of unit length with respect to the infinity norm and δ ∈ [0, ϵ) is a parameter indicating the point’s position within the box. The set

R_{n} = \{x^{∗} + (nϵ + δ)p⃗ | p⃗ ∈ ℝ^{m}, \left∥{ \atop p}\right∥_{∞} = 1, ϵ fixed, n fixed, 0 ≤ δ ≤ ϵ\}

is the set illustrated in Figure 1. Assuming that a uniform subdivision of mesh ϵ has been achieved, the ring R_{n} consists of a union of boxes, each of which is of distance nϵ to (n + 1)ϵ in ∥∘∥_{∞} from the optimizer x^{∗}. We will call both this set of boxes and the set consisting of their union “the ring R_{n}.” In our analysis we assume that R_{n} will be left in the list upon completion of Algorithm 1 provided there exists at least one box in R_{n} that is left in the list. This is valid, since our goal is to obtain an upper bound on the number of boxes left in the list.

Figure 1. Reducing the analysis to one dimension. n = 3. The ring R_{n} is delineated in bold lines.

1or parts of boxes, if x^{∗} is not at a vertex of a sub-box

Our analysis will proceed as in the one-dimensional case: We take an arbitrary box X^{(1)} ⊂ R_{n} and consider whether the algorithm will reject it. Let x^{+} denote the current midpoint used in the midpoint test. Then, directly from Step 11 of Algorithm 1, a box X^{(1)} will be rejected provided

lb(F(X^{(1)})) > f(x^{+})

i.e.

−lb(F(X^{(1)})) < −f(x^{+})

which is equivalent to

(3)  lb(f(X^{(1)})) − lb(F(X^{(1)})) < lb(f(X^{(1)})) − f(x^{+}).

But

lb(f(X^{(1)})) − lb(F(X^{(1)}))
= lb(f(X^{(1)})) − ub(f(X^{(1)})) + ub(f(X^{(1)})) − lb(F(X^{(1)}))
≤ lb(f(X^{(1)})) − ub(f(X^{(1)})) + ub(F(X^{(1)})) − lb(F(X^{(1)}))

= w(F(X^{(1)}) − w(f(X^{(1)})) (4)

Combining (3) and (4), we see that a sufficient condition for a box to be rejected is

(5)  w(F(X^{(1)}) − w(f(X^{(1)})) ≤ lb(f(X^{(1)})) − f(x^{+}).

Now, following the pattern of representation of points in R_{n}, we let x_{δ_{0}} = x^{∗} + (nϵ+δ_{0})p⃗_{0} be a point in X^{(1)} with smallest function value: f(x_{δ_{0}}) = min_{x∈X^{(1)}} f(x). Then (5) becomes

(6)  w(F(X^{(1)})) − w(f(X^{(1)})) < f(x_{δ_{0}}) − f(x^{+}).

We will think of f̃(t) = f(x^{∗} + tp⃗) as a univariate function of t. Suppose that nϵ is sufficiently small to ensure that f̃(t) is monotonic in the interval [0, (n + 1)ϵ] for every p⃗ ∈ ℝ^{m} with \left∥{ \atop p}\right∥_{∞} = 1. (This is possible whenever x^{∗} is a strict local minimizer.) Then a p⃗ can be found such that, with x_{1} = x^{∗} + nϵp⃗, f(x_{δ_{0}}) = f(x_{1}). (See Figure 1.) From (6), we thus see that

(6′) w(F(X^{(1)})) − w(f(X^{(1)})) < f(x^{∗} + nϵp⃗) − f(x^{+})

∀p⃗ ∈ ℝ^{m} with \left∥{ \atop p}\right∥_{∞} = 1

is a sufficient condition that the box X^{(1)} be rejected in Step 11 of Algorithm 1.

The theory.

Consistent with the assumptions above, the theory reflects local behavior of the algorithm. That is, it describes the behavior of the algorithm, given a sufficiently small initial box containing a single global minimum. As stated above, in practical contexts with practical algorithms, the monotonicity test (checking that the gradient does not contain zero), combined with checking the interval objective function values (i.e. the midpoint test) effectively eliminate boxes in the nonconvex regions between optima. Thus, the theory accurately reflects behavior of such practical algorithms. We do not prove this here^{2}, but illustrate it in §3 with a numerical example.

Our first theorem gives an upper bound on the number of boxes in the cluster around a global minimizer, when we assume that the minimizer occurs at the vertex of one of the boxes produced by the algorithm. The second theorem, though proven analogously, generalizes the first theorem to the case when the global minimizers occur in the interior of the boxes in question.

Suppose that the objective function f : X ⊂ ℝ^{m} → ℝ has a Taylor expansion in some neighborhood of x^{∗}, and suppose that f attains its global minimum at an interior point x^{∗} (i.e., ∥g(x^{∗})∥ = 0, where g is the gradient of f).

Theorem 1. Assume that the global minimizer x^{∗} for which the analysis holds occurs at a vertex of one of the boxes produced by the algorithm, so that x^{+} = x^{∗}. In addition to the general assumptions above, assume that the Hessian matrix G(x) of the objective function f is positive definite in the set of sub-boxes C of the initial box X in which the cluster occurs. Denote the eigenvalues of G by λ_{1} ≤ λ_{2} ≤ ⋯ ≤ λ_{m}. Set

λ_{1,0} = min_{x∈C} λ_{1} > 0.

Also assume that the interval extension F of f is of order α. Then, upon completion of the algorithm (i.e. when each box in the list has width at most ϵ), the maximum number of boxes corresponding to x^{∗} left in the list will be

(7)  N = \{2⌊\sqrt{\frac{2K}{λ_{1,0}}} · \sqrt{ϵ^{α−2}}⌋ + 1\}^{m} ,

where ⌊r⌋ stands for the integer part of the positive real number r, and where we assume that each box in the list has width ϵ.

Proof. By assumption, the objective function f has a Taylor expansion about x^{∗} of the form

f(x) = f(x^{∗}) + g(x^{∗})^{T} (x − x^{∗}) + (x − x^{∗})^{T} \frac{G(ξ)}{2}(x − x^{∗}).

where g is the gradient of f and ξ is a point of the form ξ = x^{∗} + θx for some θ (0 ≤ θ ≤ 1). When x = x^{∗} + nϵp⃗ and g(x^{∗}) = 0, the above becomes

(8)  f(x^{∗} + nϵp⃗) − f(x^{∗}) = \frac{1}{2}(nϵ)^{2}p^{T} G(ξ)p⃗,

2We will rigorously analyze the global behavior of the algorithm in a separate work, using the results presented here and similar techniques.

where ξ = x^{∗} + tp⃗ for some t ∈ (0, nϵ). Since F is of order α,

(9)  w (F(X^{(1)})) − w (f(X^{(1)})) ≤ Kw (X^{(1)})^{α} ,

where K does not depend on ϵ. Since G = G(ξ) is positive definite, a set of eigenvectors forms an orthonormal basis for ℝ^{m}. Denote these eigenvectors by {v⃗_{i}}_{i=1}^{m}, so that p⃗ = ∑_{i=1}^{m} α_{i}v⃗_{i} for some set {α_{i}}_{i=1}^{m}. Then

Gp⃗ = ∑_{i=1}^{m} α_{i}λ_{i}v⃗_{i},

so

(10)  p^{T} Gp⃗ = (∑_{i=1}^{m} α_{i}v⃗_{i}) · (∑_{i=1}^{m} α_{i}λ_{i}v⃗_{i})

= ∑_{i=1}^{m} α_{i}^{2}λ_{i}

≥ ∥p⃗∥_{2}^{2}λ_{1,0} ≥ λ_{1,0},

since

min_{∥p⃗∥_{∞}=1} ∥p⃗∥_{2}^{2} = 1.

Combining (8) and (10) then gives

(11)  f(x^{∗} + nϵp⃗) − f(x^{∗}) ≥ \frac{1}{2}(nϵ)^{2}λ_{1,0}.

Now, (9) and (11) combined with the sufficient condition in (6′) show that

(12)  Kw(X^{(1)})^{α} < \frac{1}{2}(nϵ)^{2}λ_{1,0}

is a sufficient condition that the algorithm reject X^{(1)}. Since (12) does not contain an explicit reference to p⃗, it is also a sufficient condition that any box in R_{n} be rejected. Solving the above inequality for n with w(X^{(1)}) = ϵ, we get

(13)  n > \sqrt{\frac{2K}{λ_{1,0}}} · \sqrt{ϵ^{α−2}}

Thus the maximum index n for a box X^{(1)} to stay in the list is

N_{1} = ⌊\sqrt{\frac{2K}{λ_{1,0}}} · \sqrt{ϵ^{α−2}}⌋ + 1.

The ring R_{n} will thus be rejected provided n obeys (13). However, without loss of generality assume that Algorithm 1 has produced a uniform subdivision of boxes of width ϵ in each direction^{3}. Then set of all possible boxes inside the ring R_{n} defined by (13) is a cube of side length 2nϵ. The conclusion then follows. □

3Note that each coordinate width is the coordinate width of the original box divided by a power of two. Thus, if a coordinate width is larger than a corresponding coordinate width of an adjacent box, we could subdivide the box in that direction. We could thus obtain a uniform width and more boxes, not altering the conclusion of the theorem.

Corollary 1. Suppose λ_{1,0} > 0. If:

(1) α < 2, then there may exist a severe cluster, i.e. the number of boxes in the list associated with x^{∗} may increase without bound as ϵ becomes small;

(2) α = 2, then the cluster is not serious, but there may always be a constant number N > 1 of boxes in the list associated with x^{∗}, no matter how small ϵ is;

(3) α > 2, then there is no cluster, i.e. for sufficiently small ϵ the optima obtained by use of the midpoint test are contained in single boxes, unless the optima lie on boundaries of boxes.

In practical situations, where x^{∗} is not both the best approximate minimizer and a vertex of one of the boxes, we feel that the conclusions of Theorem 1 still hold qualitatively. The following theorem proves this, under a technical assumption.

Theorem 2. Let the setting be as in Theorem 1, but do not assume that x^{∗} occurs at a vertex, nor that f(x^{∗}) = f(x^{+}). However, do assume that the current midpoint value f(x^{+}) occurs at a vertex x^{v}. Let λ_{1,0} be as in Theorem 1, and define

λ_{m,0} = max_{x∈C} λ_{m},

where C is as in Theorem 1. Then, under the other assumptions in Theorem 1, the number of boxes corresponding to x^{∗} left in the list will not exceed

(7′)  Ñ = \{\sqrt{\frac{2Kϵ^{α−2}}{λ_{1,0}} + \frac{m}{4} \frac{λ_{m,0}}{λ_{1,0}}} + 1\}^{m} .

Proof of Theorem 2. The proof is analogous to the proof of Theorem 1, but with an extra term appearing in the analogues of (8) and (12). In particular, let x^{v} denote the vertex nearest the minimizer x^{∗}. Then x^{v} takes the place of x^{∗} in the analogues to (8) and (11). A relation between f(x^{ν}) and f(x^{∗}) is given by

(14)  f(x^{v}) = f(x^{∗}) + \frac{δ^{2}}{2} p_{v}^{T} G(ξ_{v})p⃗_{v}

for some p⃗_{v} with \left∥{ \atop p_{v}}\right∥_{∞} = 1 and δ ≤ ϵ/2. However,

p_{v}^{T} G(ξ_{v})p⃗_{v} ≤ mλ_{m,0}

since \left∥{ \atop p}\right∥_{∞} ≤ \sqrt{m}∥p⃗∥_{2}, ∀p⃗ ∈ ℝ^{m}, so

(15)  f(x^{v}) ≤ f(x^{∗}) + \frac{ϵ^{2}}{8} mλ_{m,0}.

Thus, (8), (11) and (15) give

f(x^{∗} + nϵp⃗) − f(x^{v}) = f(x^{∗} + nϵp⃗) − f(x^{∗}) + f(x^{∗}) − f(x^{v})

≥ \frac{1}{2}(nϵ)^{2}λ_{1,0} − \frac{ϵ^{2}}{8} mλ_{m,0}. (11′)

Combining (11′) with (9) and the sufficient condition (6′) now gives

(12′)  Kw(X^{(1)})^{α} < \frac{1}{2}(nϵ)^{2}λ_{1,0} − \frac{ϵ^{2}}{8} mλ_{m,0}

as a sufficient condition that the algorithm reject X^{(1)}. As in Theorem 1, we solve this equation for n to obtain

(13′)  n^{2} > \frac{2Kϵ^{α−2}}{λ_{1,0}} + \frac{m}{4} \frac{λ_{m,0}}{λ_{1,0}} .

The result now follows analogously to that of Theorem 1. □

Remarks.

Remark 1. In the corollary we say “there may ⋯ ” rather than “there must” because the theorem gives an upper bound, and not a precise value, for the number of boxes. In problems in which the lower bound given by the interval extension F is exact, there will be no cluster even though F is not of high order. Note, however, that in the experiments in the next section, the bounds correspond very closely with the actual results.

Remark 2. If the Hessian matrix is not positive definite, then λ_{1,0} = 0 and the above discussion is no longer valid. We need to have a better interval extension (at least of order 3) to obtain conclusions as in the above corollary. The midpoint test will thus be ineffective in this case.

Remark 3. The requirements for the proofs of Theorem 1 and Theorem 2 to hold are essentially that the objective function have two continuous derivatives in some neighborhood of a global minimum and that the Hessian matrix be nonsingular at the global minimum in question (in the local analysis). Thus the conclusions apply to a large range of functions. In fact, even these restrictions are somewhat technical: some extension of the analysis to less smooth functions and to functions whose derivative vanishes at the minimum in question are given in [6] for the onedimensional case.

Remark 4. When G is constant (10) becomes an equality for p⃗ = v⃗_{1}. Thus, the bound in Theorem 1 is, in this sense, sharp.

Remark 5. Of course, in practice roundoff error could cause actual algorithms to behave contrary to this theory for extremely small ϵ. However, see the results in the next section.

Remark 6. The assumption on f(x^{+}) in Theorem 2 is for convenience only, and is not crucial. For the more general case, we simply base the Taylor expansion at a current midpoint. Refer to our analysis of the univariate case in [6].

3. Numerical results

We implemented Algorithm 1 in ACRITH-XSC (also called FORTRAN-SC; see [2]) on an IBM 3090. We terminated the procedure when all of the boxes in the final list had the same width ϵ. (ϵ is largest number in the computation less than the specified tolerance ϵ′.) Below, we give numbers of intervals in the final list, together with the corresponding stopping tolerance ϵ′ (the first lines in the tables) and the given initial intervals (the first columns in the tables).

Example 1. We consider the function

f(x_{1}, x_{2}) = (x_{1} + 1)^{2} + (x_{2} − 1)^{2}

on the initial box [−10, 20] × [−3, 10]. Here x^{∗} = (−1, 1) and f^{∗} = 0.

(A) \underline{Interval extension of order 1}. We used the extension

F(X) = X_{1}(X_{1} − 2) + X_{2}(X_{2} − 2) + 2,

with convergence order 1. The number of boxes left in the final list are shown with corresponding values of ϵ′ in the first row of Table 1.

(B) \underline{Interval extension of order 2}. We used the mean value form

F(X) = f(c) + \frac{f′(c)}{2}(X − c),

where c is the midpoint of X. The number of boxes left in the final list are shown with corresponding values of ϵ′ in the first row of Table 1.

initial  ϵ′_{1} = 10^{−1} ϵ′_{2} = \frac{1}{2} · 10^{−1}  ϵ′_{3} = \frac{1}{4} · 10^{−1} ϵ′_{4} = \frac{1}{8} · 10^{−1} ϵ′_{5} = \frac{1}{16} · 10^{−1}^{α = 1  119  236  467  927  1057} α = 2  13  12  14  12  11

Table 1. Number of boxes, example 1.

Example 2. Six hump camel back function,

f(x_{1}, x_{2}) = 4x_{1}^{2} − 2.1x_{1}^{4} + \frac{1}{3}x_{1}^{6} + x_{1}x_{2} − 4x_{2}^{2} + 4x_{2}^{4}

on the initial box [−10, 20] × [−3, 10]. This function has two global minima and a total of fifteen critical points (obtained with INTBIS [5]) in this region. However, the results below show that the algorithm behaves over the entire region as predicted by the local theory of §2.

(A) \underline{Interval extension of order 1}. We used the extension

F(X) = 4X_{1}^{2} − 2.1X_{1}^{4} + \frac{1}{3}X_{1}^{6} + X_{1}X_{2} − 4X_{2}^{2} + 4X_{2}^{4}

with convergence order 1. The number of boxes left in the final list are shown with corresponding values of ϵ′ in the first row of Table 2.

(B) \underline{Interval extension of order 2}. We use the mean value form for the second order interval extension. The number of boxes left in the final list are shown with corresponding ϵ′ in the second row of Table 2.

initial  ϵ′_{1} = 10^{−1} ϵ′_{2} = \frac{1}{2} · 10^{−1}  ϵ′_{3} = \frac{1}{4} · 10^{−1} ϵ′_{4} = \frac{1}{8} · 10^{−1} ϵ′_{5} = \frac{1}{16} · 10^{−1}^{α = 1  105  149  291  571  1140} α = 2  35  23  29  29  27

Table 2. Number of boxes, example 2.

On the other hand, (7) provides the following inferences.

(1) If the extension is of order 1 (α = 1), the number of boxes in the list should increase without bound as we decrease the size of the boxes. When α = 2 the number of boxes in the list should be approximately constant as the final box size decreases.

(2) In the two dimensional case, the number of boxes in the list will be proportional to N^{2}. When α = 1, we have N_{1}^{2}/N_{2}^{2} ≈ ϵ_{1}/ϵ_{2}, where N_{1} corresponds to ϵ_{1} and N_{2} corresponds to ϵ_{2}.

The results in tables 1 and 2 closely match these inferences.

4. Conclusions and future work

The above results, both theoretical and numerical, show that interval extensions of order at least 2 should be used, if possible, if only the midpoint test is used to discard boxes. This is especially true when the initial interval contains the global minimizer. Additionally, we have the following possibilities for improvement and further work.

1. Incorporate acceleration devices such as an interval Newton method. (Note that this is standard practice in interval arithmetic-based branch and bound algorithms; see, for example [4] for one of the earlier explanations of this process.)

2. Use higher order interval extensions, when possible. When the Hessian matrix is positive definite, an interval extension of order greater than two together with the simple midpoint test should result in an efficient algorithm with no cluster, even without acceleration procedures.

3. When the Hessian matrix is no longer positive definite at the optimum (but is still nonnegative definite), a higher order extension may be necessary to avoid a cluster. Our investigations in the immediate future will involve delineating classes of

multivariate objective functions for which higher-order interval extensions are readily available. In particular, we will see when the arguments in [3] and in §2.4 of [8] generalize.

References

1. Alefeld, Götz, and Herzberger, Jürgen: Introduction to Interval Computations. Academic Press, New York (1983).

2. Bleher, J. H., Rump, S. M., Kulisch, U., Metzger, M., Ullrich, C., and Walter, W.: Fortran-SC — A study of a Fortran extension for engineering and scientific computation with access to ACRITH. Computing \underline{39}, 93–110 (1987).

3. Cornelius, H. and Lohner, R.: Computing the range of values of real functions with accuracy higher than second order. Computing \underline{33}, 331–347 (1984).

4. Hansen, E. R.: Global optimization using interval analysis– the multidimensional case. Numer. Math. \underline{34}, 247–270 (1980).

5. Kearfott, R. B. and Novoa, M.: INTBIS, A Portable Interval Newton / Bisection Package (Algorithm 681) ACM Trans. Math. Software \underline{16}, 152–157 6. Kearfott, R. B. and Du K.: The cluster problem in global optimization, the

univariate case. Computing Supplement \underline{9}, 117–127 (1992).

7. Moore, Ramon E.: Methods and Applications of Interval Analysis. SIAM, Philadelphia (1979).

8. Neumaier, A.: Interval Methods for Systems of Equations. Cambridge University Press, Cambridge, England (1990).

9. Ratschek, H. and Rokne, J.: Computer Methods for the Range of Functions. Horwood, Chichester, England (1984).

10. Ratschek, H. and Rokne, J.: New Computer Methods for Global Optimization. Wiley, New York (1988).

Department of Mathematics, U.S.L. Box 4-1010, Lafayette, LA 70504-1010