On Stopping Criteria in Verified Nonlinear Systems or Optimization Algorithms

R. B. Kearfott and G. W. Walster

University of Louisiana at Lafayette and Sun Microsystems

Traditionally, iterative methods for nonlinear systems use heuristic domain and range stopping criteria to determine when accuracy tolerances have been met. However, such heuristics can cause stopping at points far from actual solutions, and can be unreliable due to the effects of roundoff error or inaccuracies in data.

In verified computations, rigorous determination of when a set of bounds has met a tolerance can be done analogously to the traditional approximate setting. Nonetheless, the range tolerance possibly cannot be met. If the criteria are used to determine when to stop subdivision of n-dimensional bounds into subregions, then failure of a range tolerance results in excessive, unnecessary subdivision, and could make the algorithm impractical.

On the other hand, interval techniques can detect when inaccuracies or roundoff will not permit residual bounds to be narrowed. These techniques can be incorporated into range thickness stopping criteria that complement the range stopping criteria. In this note, the issue is first introduced and illustrated with a simple example. The thickness stopping criterion is then formally introduced and analyzed. Third, inclusion of the criterion within a general verified global optimization algorithm is studied. An industrial example is presented. Finally, consequences and implications are discussed.

Categories and Subject Descriptors: G.1.5 [Numerical Analysis]: Roots of Nonlinear Equations—error analysis, iterative methods, systems of equations; G.1.6 [Numerical Analysis]: Optimization—unconstrained optimization, constrained optimization, global optimization

General Terms: Algorithms, Performance, Reliability, Verification

Additional Key Words and Phrases: stopping criteria, iterative methods, verified computations, global optimization, GlobSol

This work was partially funded by National Science Foundation grant DMS-9701540 and by a Sun Microsystems Cooperative Research grant. Name: R. B. Kearfott Address: Department of Mathematics, University of Louisiana at Lafayette, Box 4-1010, Lafayette, LA 70504-1010, USA email: rbk@usl.edu Name: G. W. Walster Address: MTV12-40, Sun Microsystems, Inc., 2550 Garcia Avenue, Mountain View, CA 94043 email: Bill.Walster@Eng.Sun.COM

2 · R. B. Kearfott

1. INTRODUCTION: BACKGROUND AND SIMPLE EXAMPLES

In a general implementation of Newton’s method, the steepest descent method, or a hybrid method to compute solutions of f(x) = 0, where f : ℝ^{n} → ℝ^{n}, a domain tolerance of a form similar to

δ_{k} = \frac{∥x^{(k+1)} − x^{(k)}∥}{max \{∥x^{(k+1)}∥, 1\}} < δ  (point domain tolerance condition) (1)

is used, for some domain stopping tolerance δ. The well-known idea is that, if x^{(∗)} solves f(x^{(∗)}) = 0 exactly, then

∥x^{(k+1)} − x^{(∗)}∥ ∼ C∥x^{(k)} − x^{(∗)}∥^{2}  (2)

for some constant C and for x^{(k)} sufficiently close to x^{(∗)}. The quadratic convergence assumption represented in (2) combined with the triangle inequality leads to

∥x^{(k+1)} − x^{(∗)}∥ ˙≤ C (\frac{∥x^{(k+1)} − x^{(k)}∥}{1 − C∥x^{(k)} − x^{(∗)}∥})^{2} ≤ 2C∥x^{(k+1)} − x^{(k)}∥^{2}  (3)

for x^{(k)} sufficiently close to x^{(∗)}. The approximation (3) justifies the domain heuristic (1).

Similarly, f may be flat (with unknown flatness) near a solution x^{(∗)}, or else merely a point x^{(k)} with a small residual ∥f(x^{(k)})∥ is desired. In such instances, the following range tolerance is appropriate.

ρ_{k} = ∥f(x^{(k+1)})∥ ≤ ϵ  (point range tolerance condtion)  (4)

for some range tolerance ϵ.

For illustration, take the simple example

Example 1. Consider

f(x) = x^{2} − 2,

with x ∈ [1, 2] and x̌ = 1.5.

Newton’s method for Example 1 converges rapidly for any starting point x^{(0)} > \sqrt{2}, and in particular for x^{(0)} = x̌ = 1.5. Since f′ (x^{(∗)}) = 2\sqrt{2} is non-singular, the domain and range tolerances are roughly equivalent, and the iterations in Table 1 are obtained^{1}.

It is seen in Table 1 that the domain tolerance (1) and the range tolerance (4) reflect actual error behavior well until the errors become small enough for roundoff error to dominate the computation. For this example, a serious problem could occur only if δ were set to be on the order of the square of the machine epsilon, or ϵ were set to be on the order of the machine epsilon^{2}.

However, uncertainties in data or computational errors can lead to loss of accuracy orders of magnitude greater than the machine epsilon. Under these conditions,

1Here, IEEE arithmetic was used. In lines with errors of 0, the machine representable approximations being subtracted were identical. 2The machine epsilon for this IEEE arithmetic machine is approximately 2.220446049250313 × 10^{−16}.

Stopping Tests for Verified Optimization · 3

k  x^{(k+1)}  δ_{k}^{2}  ρ_{k}  ∥x^{(k+1)} − x^{(∗)}∥
0 \underline{1}.50000000000000
1 \underline{1.41}666666666667 3.46 × 10^{−3} 6.94 × 10^{−3} 2.45 × 10^{−3}
2 \underline{1.41421}568627451 3.00 × 10^{−6} 6.01 × 10^{−6} 2.12 × 10^{−6}
3 \underline{1.41421356237}469 2.26 × 10^{−12} 4.51 × 10^{−12} 1.59 × 10^{−12}
4 \underline{1.4142135623731}0 1.27 × 10^{−24} 4.44 × 10^{−16} 0
5 \underline{1.414213562373}09 2.47 × 10^{−32} 4.44 × 10^{−16} 2.22 × 10^{−16}
6 \underline{1.4142135623731}0 2.47 × 10^{−32} 4.44 × 10^{−16} 0

Table 1. Stopping tolerance behavior for the point Newton method with function x^{2} − 2.

the maximum attainable accuracy in a Newton iteration can be difficult to predict with floating-point iterations. In such instances, the domain and range tolerances will not terminate the algorithm, and behavior as in iterations 4 through 6 in Table 1 will occur, although the iteration may be more erratic for more complicated functions and with more than one variable, and the actual magnitudes of the values δ_{k}^{2} and ρ_{k} may not be predictable.

Such uncertainties can be automatically handled with multidimensional interval Newton methods, as described in the texts on interval computations, such as [7], [1], [8], [3], [5]. In one dimension, interval Newton methods are of the form

x^{(k+1)} = \{x̌^{(k)} − \frac{f (x̌^{(k)})}{f′(x^{(k)})}\} ∩ x^{(k)}.  (5)

(Here, interval quantities are denoted by boldface, so that x^{(k)} is an interval, x̌^{(k)} is a representative point, often taken to be the midpoint of x^{(k)}, and f′(x^{(k)}) is an interval extension of the derivative of f over x^{(k)}. ) A domain tolerance condition corresponding to (1) is

\overline{δ}_{k} = max_{1≤i≤n} \begin{cases}w & (x_{i}^{(k)}) \\ max & \{|x_{i}^{(k)}| , 1\} \end{cases}\} < δ (interval domain tolerance condition), (6)

where w(x_{i}^{(k)}) represents the width of the i-th coordinate interval. A range tolerance condition corresponding to (4) is

0 ∈ f_{i}(x^{(k)}) and w(f_{i}(x^{(k)})) < ϵ, 1 ≤ i ≤ n,

or, equivalently,

\overline{ρ}_{k} = max_{1≤i≤n} |f_{i}(x^{(k)})| < ϵ  (interval range tolerance condition).  (7)

Application of the interval Newton method (5) to Example 1, with starting bounds x^{(0)} = [1, 2] gives Table 2. In Table 2, the library INTLIB [6] and Fortran 90 module INTERVAL ARITHMETIC [4] were used for simulated directed roundings; the intervals in the table were rounded outward for display. The actual iterates became exactly stationary at iteration number 4.

Analogously to the floating point case, \overline{δ}_{k} and \overline{ρ}_{k} are reasonable quantities to check for stopping tolerances, except that the interval versions represent mathematically rigorous domain and range error bounds, rather than heuristic values.

\begin{matrix}k & x^{(k+1)} & \overline{δ}_{k} & \overline{ρ}_{k} \\ 0 & [1.00000000000000, 2.00000000000000] & 5.00 × 10^{−1} & 2.00 × 10^{0} \\ 1 & [\underline{1}.37499999999999, \underline{1}.43750000000001] & 4.35 × 10^{−2} & 1.09 × 10^{−1} \\ 2 & [\underline{1.414}06249999999, \underline{1.414}41761363637] & 2.51 × 10^{−4} & 5.77 × 10^{−4} \\ 3 & [\underline{1.4142135}5929452, \underline{1.4142135}6594718] & 4.70 × 10^{−9} & 1.01 × 10^{−8} \\ 4 & [\underline{1.414213562373}09, \underline{1.414213562373}10] & 4.71 × 10^{−16} & 1.33 × 10^{−15} \\ 5 & [\underline{1.414213562373}09, \underline{1.414213562373}10] & 4.71 × 10^{−16} & 1.33 × 10^{−15} \\ 6 & [\underline{1.414213562373}09, \underline{1.414213562373}10] & 4.71 × 10^{−16} & 1.33 × 10^{−15} \end{matrix}

Table 2. Stopping tolerance behavior for the interval Newton method with function x^{2} − 2.

However, as in the floating point case, the tolerance conditions (6) and (7) possibly cannot be met if δ or ϵ are set too small. With inclusion-monotonic interval arithmetic, the interval Newton iteration (5) must eventually become stationary. However, multidimensional versions of (5) often become stationary far before reaching an accurate enclosure for a solution, and many iterations may be spent with only small progress towards the limiting endpoints. As an example, take

Example 2. Take

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

with x_{1}^{(0)} = x_{2}^{(0)} = [−1, 1].

In Example 2, the interval extension of the Jacobi matrix over x^{(0)} is

( \begin{matrix}[0, 3] & −1 \\ 1 & [0, 3] \end{matrix} ) ,

and it does not contain any singular matrices. Furthermore, the system resulting from applying the inverse midpoint preconditioner also does not contain any singular matrices, and there is a unique solution within x^{(0)}. However, both the interval Gauss–Seidel method and interval Gaussian elimination are stationary at x^{(0)}. In such instances, subdividing the initial bounds and applying interval techniques to sub-regions leads to success. In this case, using a starting box of x_{1}^{(0)} = x_{2}^{(0)} = [−0.2, 0.1] (and also perhaps larger starting boxes) results in convergence to narrow bounds around the true solution x^{(∗)} = (0, 0), while other subregions, such as say x_{1} = [0.1, 1], x_{2} = [0.1, 1], can be rejected with other interval techniques. (For an explanation of these terms and algorithms, see [5].)

On the other hand, the interval Newton method may become stationary at some x^{(k)} due to roundout error, as in Table 2, or due to uncertainty in the data^{3}. In such cases, it is not appropriate to subdivide x^{(k)} further. To illustrate the situation where uncertainties in the data lead to questions about stopping tolerances, consider the following illustrative example.

Example 3. Suppose the task is to

minimize ϕ(x) = \frac{1}{2} \{(x_{1} − [1, 2])^{2} + (x_{2} − [3, 4])^{2}\} .

3See the next section for an example.

Here, the intervals [1, 2] and [3, 4] represent uncertainty in the problem specification. Subject to this uncertainty, find maximum possible and minimum possible values of x_{1} and x_{2}.

If an interval Newton method, starting with initial coordinate bounds x^{(0)} = ([−10, 10], [−10, 10], is applied to the gradient system

∇ϕ = (x_{1} − [1, 2], x_{2} − [3, 4])^{T} = (0, 0)^{T} ,

then the first iteration of an interval Gauss–Seidel method gives x^{(1)} = ([1, 2], [3, 4])^{T} (to within roundout error) and all subsequent iterates are stationary at x^{(k)} = x^{(1)}, k ≥ 1. Thus, \overline{δ}_{k} = 1/2 and \overline{ρ}_{k} = 2 for k ≥ 1. In some algorithms, lack of tolerance satisfaction causes the box ([1, 2], [3, 4])^{T} to be bisected, say, into the boxes ([1, 1.5], [3, 4])^{T} and ([1.5, 2], [3, 4])^{T}. However, since the box ([1, 2], [3, 4]) represents the sharpest possible bounds on the answer in this case, such bisection is not appropriate.

Verified global optimization algorithms also sometimes use the width of the interval extension of the objective function to determine when a box should be accepted as small enough. With this criterion,

w(ϕ(x^{(1)})) = w(([1, 2] − [1, 2])^{2} + ([3, 4] − [3, 4]) ∗ 2)
= w([0, 1] + [0, 1]) = 2,

not small.

Thus, stopping due to domain or range tolerance conditions, or stopping when the interval Newton method becomes stationary, is not sufficient. A stopping criterion that determines when further subdivision will not enable convergence will be useful.

2. THE THICKNESS STOPPING CRITERION, AND CONSEQUENCES OF ITS USE

In rigorous global optimization algorithms with objective function ϕ(x), the progress is made either by comparing bounds for the objective function over a region (interval vector) x^{(0)} or by successful iteration of an interval Newton method applied to ∇ϕ(x) = 0, with x^{(0)} as initial iterate^{4}. Subdivision of x^{(0)} is done to either enable sharper bounds on the objective ϕ or to enable the interval Newton method to converge to narrower bounds, and such subdivision leads to unnecessary inefficiency if it does not accomplish either of these things. The thickness stopping criterion explained here detects when further bisection will not lead to sharper bounds on the objective function, and can also be applied similarly to constraints or gradients.

Applied to the objective function ϕ, the criterion is simply

w(ϕ(x̌^{(k)})) ≥ η_{ϕ}w(ϕ(x^{(k)})) (objective thickness condition),  (8)

where η_{ϕ} ∈ [0, 1] is termed the “objective function thickness factor.” Smaller values of η_{ϕ} give a looser tolerance, but also allow more overestimation in the the range of ϕ.

4In nonlinear systems algorithms, just a nonlinear system is considered, and in constrained algorithms, the both the constraints by themselves and the Lagrange multiplier or Fritz–John system are considered. This will be discussed in the next section.

In Example 3 with x^{(k)} = ([1, 2], [3, 4])^{T} and x̌^{(k)} = (1.5, 3.5)^{T},

w(ϕ(x̌^{(k)})) = w([−0.5, 0.5]^{2} + [−0.5, 0.5]^{2})
= 0.5

and

w(ϕ(x^{(k)})) = 2.

Thus, the optimally small box x^{(k)} = ([1, 2], [3, 4])^{T} will be accepted without further subdivision if η_{ϕ} is chosen so η_{ϕ} < 0.25.

The following observations aid in determining appropriate values of η_{ϕ}.

Theorem 1. (Upper bound on the range overestimation) Suppose that the objective thickness condition (8) is satisfied for some η_{ϕ} ∈ (0, 1). Let ϕ^{♯}(x^{(k)}) represent the actual range of ϕ over x^{(k)}. Then the range overestimation w(ϕ^{♯}(x^{(k)})) −

w(ϕ(x^{(k)})) satisfies

w(ϕ^{♯}(x^{(k)})) − w(ϕ(x^{(k)})) ≤ \frac{1 − η_{ϕ}}{η_{ϕ}} w(ϕ(x̌^{(k)})).

Furthermore, this bound is achieved for some objective functions ϕ.

Proof.

w(ϕ^{♯}(x^{(k)})) − w(ϕ(x^{(k)})) ≤ w(ϕ(x^{(k)})) − w(ϕ(x̌^{(k)}))
≤ \frac{w(ϕ(x̌^{(k)}))}{η_{ϕ}} − w(ϕ(x̌^{(k)}))
= \frac{1 − η_{ϕ}}{η_{ϕ}} w(ϕ(x̌^{(k)})).

Furthermore, for constant ϕ(x) ≡ [a, b], w(ϕ(x̌^{(k)})) = w(ϕ(x^{(k)})) = w(ϕ^{♯}(x^{(k)})) = b − a, and the inequality becomes an equation.

Theorem 1 gives guidance on how near η_{ϕ} should be to 1 to ensure a given accuracy in the range of the objective function and, hence, a given accuracy in the global optimum. Conversely, η_{ϕ} should be chosen closer to 0 for the algorithm to be more efficient. It has been seen that, in Example 3, choosing η_{ϕ} < 0.25 will totally prevent unnecessary subdivision. Using a larger value of η_{ϕ} for this example will result in some extra subdivision, but it will not be excessive for η_{ϕ} around, say, 0.5. In general, an appropriate value of η_{ϕ} to avoid unnecessary subdivision depends both on the amount of interval uncertainty in the objective and specific characteristics of the interval evaluation. (The width overestimation can, in principle, be arbitrarily large.) When the objective ϕ is an uncertain constant, any η_{ϕ} ∈ [0, 1] will do. Otherwise, the value η_{ϕ} = 0.5 seems to work well in practice.

3. STOPPING CRITERIA IN GLOBAL OPTIMIZATION ALGORITHMS

Philosophy and details of incorporation into actual global optimization algorithms are considered here.

3.1 “Thick” Parameters: What is a Solution?

Considering Example 3, one may ask why the uncertain constants, [1, 2] and [3, 4] should not be considered as variables, that is, why not consider the four-dimensional bound-constrained problem

minimize ϕ(x) = \frac{1}{2} \{(x_{1} − x_{3})^{2} + (x_{2} − x_{4})^{2}\} , x_{3} ∈ [1, 2], x_{4} ∈ [3, 4] (9)

instead? In fact, in this particular example, the GlobSol package [2], starting with x_{1} ∈ [−10, 10], x_{2} ∈ [−10, 10], considers 21 boxes to give an answer of four boxes whose union is the box x^{(∗)} = ([1, 2], [3, 4], [1, 2], [3, 4])^{T}, exactly the solution that is obtained in the two variable problem with uncertain constants. But this due only to the special form of Example 3. Consider

Example 4.

minimize ϕ(x) = (x − a)^{2} + 10a^{2} with a ∈ [−1, 1], for x ∈ [−10, 10].

With η_{ϕ} = 0.2, GlobSol processes exactly two subintervals to obtain a solution list consisting of two subintervals whose union is the exact solution set x^{(∗)} = [−1, 1]. However, if a is considered as a variable, that is, if GlobSol solves the problem

minimize ϕ(x) = (x_{1} − x_{2})^{2} + 10x_{2}^{2}  for x ∈ ([−10, 10], [−1, 1])^{T} ,  (10)

GlobSol processes two boxes total, and returns a small box centered on the unique solution (0, 0), a very different solution from that to the original formulation in Example 4. The solution set for Example 4 contains the set of all solutions to the univariate problem for a any fixed point in [1, 2], while the solution to the variant in (10) contains only the solution for that a that gives the minimum value of ϕ. One can think of the interval value of ϕ, [0, 14] in the case of GlobSol’s solution to Example 4 with η_{ϕ} = 0.2, as bounds on the minimum value of the minimum and the maximum value of the minimum, with some overestimation if the constant appears in more than one place. (In Example 4, the actual range of ϕ over the minimum values is [0, 10]. See §3.4 below.)

Efficiency considerations may also play a role in how the solution set is viewed. Consider

Example 5.

minimize ϕ(x) = a(x_{1} − a)^{2} + b(x_{2} − b)^{2}
with a ∈ [1, 2], b ∈ [3, 4], for x_{1} ∈ [−10, 10] and x_{2} ∈ [−10, 10].

With η_{ϕ} = 0.2, GlobSol processed four boxes total, and, with η_{ϕ} = 0.5, GlobSol processed 10 boxes total^{5}. In both cases, the union of the output boxes consisted of the set

{(x_{1}, x_{2}) | x_{1} ∈ [0.8408, 2.1592], x_{2} ∈ [2.8333, 4.1667], }

5GlobSol was not able to complete without considering more than 20,000 subregions when η_{ϕ} = 0.9.

to within roundout. (Note that this is an overestimate of the actual solution set x^{(∗)} = ([1, 2], [3, 4])^{T}.) If Example 5 is reformulated as

minimize ϕ(x) = x_{3}(x_{1} − x_{3})^{2} + x_{4}(x_{2} − x_{4})^{2}
with x_{3} ∈ [1, 2], x_{4} ∈ [3, 4], for x_{1} ∈ [−10, 10] and x_{2} ∈ [−10, 10],

(11)

then GlobSol was unable to complete a search in an hour of processor time, corresponding to 89954 boxes processed. Although Problem (11) is similar to Problem (9) over the region in question, both problems have a two-dimensional set of solutions in ℝ^{4}. The small but significant overestimation in Problem (11), due to x_{3} and x_{4} appearing more than once, is enough to cause GlobSol to subdivide the region into too many small sub-boxes surrounding the solution set.

Recapitulating:

—The solution set for a problem with thick constants is different from the solution

set of the corresponding problem with the constants considered as variables. —The solution set to a problem with non-zero width (“thick”) constants contains

information about the solution to a kind of minimax problem. The corresponding interval values of the objective function approximate the minimum minimum and maximum minimum, and are in general better approximations when the widths of the constants are smaller.

—Combined with the stopping criteria described here, formulating the problem with thick constants can be more efficient than formulating the problem with variables only.

3.2 Structuring the Stopping Criteria

Global optimization problems are generally of the form

minimize ϕ(x)
subject to \{ \begin{matrix}c(x) & = & 0 and \\ g(x) & ≤ & 0, \end{matrix} \}

(12)

Here ϕ : x ⊂ ℝ^{n} → ℝ, c : x → ℝ^{m_{1}}, and g : x → ℝ^{m_{2}}, where x is an interval vector x = ([\underline{x}_{1}, \overline{x}_{1}], …, [\underline{x}_{n}, \overline{x}_{n}])^{T}. In a step of a global optimization code such as GlobSol, a current box x^{(k)} is processed, roughly with the following steps.

Algorithm 1. (General Steps in verified global search)

(1) The stopping tolerances are checked.

(2) Constraint propagation techniques are applied to narrow the coordinate bounds x_{i}^{(k)}, 1 ≤ i ≤ n. (See [5, Ch. 7] and the references therein.

(3) It is checked whether the lower bound \underline{ϕ}(x^{(k)}) > \overline{f}, where \overline{f} is a previously computed rigorous upper bound on the global optimum. If it is, then reject x^{(k)}, and proceed to step 9.

(4) If consideration of x^{(k)} leads to a smaller value of \overline{f}, then step 3 is applied with the boxes x^{(∗)} from the list of previously stored boxes.

(5) It is checked whether 0 ∈ c_{i}(x^{(k)}), 1 ≤ i ≤ m_{1} and whether \overline{g}_{i} ≤ 0, 1 ≤ i ≤ m_{2}.

(6) It is checked whether 0 ∈ L, where L is a generalized Lagrange multiplier system or Fritz–John system corresponding to Problem (12).

(7) An interval Newton method is applied to the generalized Lagrange system.

(8) The stopping tolerances are checked again, and the box is either discarded, stored on a solution list, or subdivided, so the components can be stored for further processing.

(9) A new x^{(k)} is selectedfrom the list of boxes yet to be processed. If there is no such x^{(k)}, then the algorithm exits.

In steps 3 and 5, interval values of the objective, equality, and inequality constraints are used, while interval extensions of the gradients of these are needed in step 6. The sharpness of these interval bounds can possibly be improved through subdivision. Thus, application of the stopping criteria should take account of the objective ϕ, the equality constraints c, the inequality constraints g, and the gradients of these.

Different problems may benefit from different sharpnesses with which the objective and two kinds of constraints are enclosed. Also, it may be unknown what the exact effects of roundout error or thick tolerances are on the minimum achievable widths in the ranges of the objective and two kinds of constraints. Thus, a rational way of combining range tolerances, thickness conditions, gradient range tolerances, and domain tolerances is needed.

The range stopping test is structured into three parallel parts, corresponding to the objective, the equality constraints, and the inequality constraints. For simplicity, the internal structure will be explained for the objective; the internal structure of the other two is similar.

Algorithm 2. (Range and thickness criteria for the objective function)

Input: The objective range width tolerance ϵ_{ϕ} and the objective thickness factor η_{ϕ}, the present box x^{(k)}, its midpoint x̌^{(k)}, and the machine epsilon ϵ_{m}. Output: Either “x^{(k)} meets the criteria” or “x^{(k)} needs to be subdivided further.”

If w(ϕ(x^{(k)})) < max {ϵ_{ϕ}, 100ϵ_{m}} Or w(ϕ(x̌^{(k)})) > η_{ϕ}w(ϕ(x^{(k)}))

Then signal “x^{(k)} meets the criteria,”

Else signal “x^{(k)} needs to be subdivided further.” End If

End Algorithm 2

If ϕ(x^{(k)}) is thick, it cannot be made narrow, since w(ϕ(x^{(k)})) ≥ w(ϕ(x̌^{(k)})). Conversely, if ϕ(x^{(k)}) is narrow to within tolerance ϵ_{ϕ}, further subdivision will not improve the chances that \underline{ϕ}(x^{(k)}) > \overline{f} (step 3 described above in the global optimization algorithm process) can be verified. Thus, disjunction as illustrated, rather than conjunction or a user-controlled choice, is appropriate in the conditional in Algorithm 2.

In the complete stopping test, the analogue of Algorithm 2 is applied to the constraints c and the inequality constraints g, in addition to the objective ϕ. Since the importance of accuracy in the objective and constraints, as well as the efficiency with which such accuracy can be achieved, varies from problem to problem, user control of the importance of each of conditions on the objective and two types of constraints is useful. This results in the following scheme.

Algorithm 3. (Combined range width and thickness conditions)

Input:

(1) flags F_{ϕ}, F_{c}, and F_{g} indicating whether it is necessary for the range or width condition on the objective, equality, or inequality constraints, respectively, to be satisfied for convergence to be signalled,

(2) the width tolerances ϵ_{ϕ}, ϵ_{c}, and ϵ_{g} for the objective, equality, and inequality constraints, respectively,

(3) the thickness factors η_{ϕ}, η_{c}, and η_{g} for the objective, equality, and inequality constraints, respectively,

(4) the present box x^{(k)}, its midpoint x̌^{(k)}, and the machine epsilon ϵ_{m}.

Output: Either “x^{(k)} meets the criteria” or “x^{(k)} needs to be subdivided further.”

(1) Define

(a) P = “true” if and only if Algorithm 2 returns “x^{(k)} meets the criteria”.

(b) C = “true” if and only if Algorithm 2, applied with ϵ_{c} replacing ϵ_{ϕ} and with each c_{i}, 1 ≤ i ≤ m_{1}, replacing ϕ, returns “x^{(k)} meets the criteria”.

(c) G = “true” if and only if Algorithm 2, applied with ϵ_{g} replacing ϵ_{ϕ} and with each g_{i}, 1 ≤ i ≤ m_{2}, replacing ϕ, returns “x^{(k)} meets the criteria”.

(2) If

{{(F_{ϕ} ∨ F_{c} ∨ F_{g})} ∧ (P ∨ ¬F_{ϕ}) ∧ (C ∨ ¬F_{c}) ∧ (G ∨ ¬F_{g})}
∨ {{¬(F_{ϕ} ∨ F_{c} ∨ F_{g})} ∧ (P ∨ C ∨ G)}

Then signal “x^{(k)} meets the criteria,”

Else signal “x^{(k)} needs to be subdivided further.” End If End Algorithm 3

In the conditional in Algorithm 3, the user-set flags F_{ϕ}, F_{c}, and F_{g} define whether or not it is necessary for the objective, the equality constraints, or the inequality constraints, respectively, to meet the width or thickness criterion, unless all three flags are set to “false.” If all three flags are set to “false,” then the box is considered subdivided enough if at least one of the width/thickness criteria is met. This scheme, although not the only one, provides substantial flexibility with a small number of parameters. The second line of logic in step 2, handling the case when all three flags are “false,” prevents vacuously signalling “x^{(k)} meets the criteria,” and allows stopping the subdivision if any of the criteria is met. In this case, that is, when (¬F_{ϕ} ∧ ¬F_{c} ∧ ¬F_{g}), individual tolerances can be turned completely off by adjusting the width tolerances and the thickness factors. For example, in the case (¬F_{ϕ} ∧ ¬F_{c} ∧ ¬F_{g}), P can be turned off by setting ϵ_{ϕ} = 0 and η_{ϕ} = 1.

Algorithm 3 includes all conditions on the ranges of the objective and constraints, but does not include either the domain conditions or the condition on the gradient.

Also, it is useful to include the domain condition (6) for cases for which the domain tolerance δ is easier to specify than the range width tolerances and range thickness factors. Finally, a condition to ensure that the gradient of the Lagrange function is narrow is useful. This condition, which can also be viewed as a scaled domain condition, has been termed maximum smear in [5, p. 157 ff and p. 175], and justified in [9] and elsewhere, is

σ_{i} < δσ,  1 ≤ i ≤ n, where  (13)
σ_{i} = |\frac{∂ϕ}{∂x_{i}} (x^{(k)})| + ∑_{j=1}^{m_{1}} |\frac{∂c_{j}}{∂x_{i}} (x^{(k)})| + ∑_{j=1}^{m_{2}} |\frac{∂g_{j}}{∂x_{i}} (x^{(k)})| and where
σ = max_{1≤i≤n} σ_{i}.

Here, |[a, b]| = max{|a|, |b|}. In the sums for σ_{i}, it would be advantageous for the partial derivatives of the c_{j} and g_{j} to be multiplied by corresponding generalized Lagrange multipliers. However, such multipliers are in general not known when the stopping test is applied.

In the global optimization system GlobSol, the range/thickness test, scaled domain tolerance (6) and maximum smear (13) are combined alternatively. The reasoning is that it is hard to predict, in general, which conditions are easily satisfiable, and alternative satisfaction causes subdivision to stop under the weakest possible conditions, thus avoiding impractical amounts of subdivision. (If the answers given by such an algorithm are not narrow enough, it is relatively easy to refine them.) This alternative formulation results in the following overall stopping criterion implemented in GlobSol.

Algorithm 4. (Overall stopping criterion)

Input:

(1) the domain tolerance δ,

(2) flags F_{ϕ}, F_{c}, and F_{g} indicating whether it is necessary for the range or width condition on the objective, equality, or inequality constraints, respectively, to be satisfied for convergence to be signalled,

(3) the width tolerances ϵ_{ϕ}, ϵ_{c}, and ϵ_{g} for the objective, equality, and inequality constraints, respectively,

(4) the thickness factors η_{ϕ}, η_{c}, and η_{g} for the objective, equality, and inequality constraints, respectively,

(5) the present box x^{(k)}, its midpoint x̌^{(k)}, and the machine epsilon ϵ_{m}.

Output: Either “x^{(k)} meets the criteria” or “x^{(k)} needs to be subdivided further.”

If The scaled diameter is small according to (6),

Then signal “x^{(k)} meets the criteria,” and Return Else If The range conditions (Algorithm 3) are satisfied,

Then signal “x^{(k)} meets the criteria,” and Return Else If The smear diameter according to (13) is satisfied,

Then signal “x^{(k)} meets the criteria,” and Return Else

signal “x^{(k)} needs to be subdivided further.” End If

End Algorithm 4

3.3 Which Coordinate Should be Subdivided?

A simple and common way of subdividing a box x^{(k)} is to bisect a coordinate to form \underline{x}^{(k)} and \overline{x}^{(k)}. For example, if x^{(k)} = ([−1, 1], [−2, 2])^{T}, then the second coordinate may be bisected to form \underline{x}^{(k)} = ([−1, 1], [−2, 0])^{T} and \overline{x}^{(k)} = ([−1, 1], [0, 2])^{T}, or the first coordinate may be bisected to form \underline{x}^{(k)} = ([−1, 0], [−2, 2])^{T} and \overline{x}^{(k)} = ([0, 1], [−2, 2])^{T}. Ideally, the coordinate to be bisected will be chosen to maximize the chance that the overall stopping criterion as defined by Algorithm 4 will be satisfied by the resulting boxes. A good heuristic for this is to order the coordinates according to the smears σ_{i} of equation (13), after removing those terms for which ϕ, c_{j}, or g_{j} is thick or narrow according to Algorithm 2. This is what is presently done in GlobSol.

3.4 Is it Advantageous to Subdivide Thick Constants?

In Example 5, due to interval dependence in the constants a and b, the estimate for ϕ was an overestimate of the actual range. Subdividing the constants, such as considering two problems, one with a ∈ [1, 1.5] and one with a ∈ [1.5, 2], will reduce the overestimation^{6}. Furthermore, especially in more complicated problems, subdivision can reduce significantly the total amount of work involved. See §3.5 below.

There are differences between subdivision of thick constants and subdivision of variables. Typically, subdivision of the constants will not lead to rapid rejection of one half of the original region, as happens when subdividing a variable interval. Thus, subdividing a multidimensional array of thick constants potentially produces a very large number of subproblems to be considered. Furthermore, supplying different subintervals is logically equivalent to supplying a different objective function, and the midpoint test, that is, step 4 of the general optimization scheme (Algorithm 1) will reject some boxes x^{(∗)} that were stored when a previous sub-interval was being processed. Similarly, the value of \overline{f} from steps 3 and 4 is different for different subintervals. Thus, a more sophisticated algorithm would need to be supplied if thick constants were subdivided automatically.

For these reasons, GlobSol only allows for subdivision of a maximum of 10 thick constants, and an example driver is supplied that only subdivides one thick constant. Nonetheless, it is possible to treat constants as variables, as explained in §3.3, when deciding how to subdivide. Further implementation effort and research are required to determine how valuable this is.

To illustrate the advantages of subdivision of a constant, a one-dimensional version of Example 5 was tried:

minimize ϕ(x) = a(x_{1} − a)^{2},  with a ∈ [1, 2], for x_{1} ∈ [−10, 10].  (14)

6The amount by which the overestimation is reduced can be estimated by Theorem 1.

With no subdivision, GlobSol gave two answer boxes whose union was enclosed in x^{(∗)} = [0.9481, 1.4996], with objective value enclosed in ϕ(x^{(∗)}) ∈ [−2.23 × 10^{−308}, 2.213]. Subdividing a into five pieces gave 10 solution boxes whose union was enclosed in x^{(∗)} = [0.9992, 2.0004], with ϕ(x^{(∗)}) ∈ [−2.23×10^{−308}, 4.840×10^{−2}].

3.5 An Industrial Example

In [10], GlobSol and the thick parameter scheme are used to solve a maintenancescheduling problem. The first case treated there is a nine-variable mixed integer programming problem:

Example 6. (A maintenance scheduling problem)

minimize \{s_{0} + ∑_{i=1}^{8} s_{i}/x_{i+1}\} /x_{1} + \{∑_{i=1}^{8} (\frac{c_{i}}{l_{i}^{β_{i}}} ) (x_{1}x_{i+1})^{β_{i}−1}\} ,

where l = (8, 7, 9, 14, 6, 15, 3, 5)^{T}, s = (100, 105, 225, 345, 165, 500, 345, 105, 345)^{T}, β = (1.7, 1.7, 2.0, 2.0, 1.7, 2.0, 1.25, 1.75)^{T}, where the c_{i} are interval constants, where the variables x_{2} through x_{9} are integers, and where the solution is to be found within the bounds x^{(0)} = ([0.5, 20.5], [0.5, 2.5], [0.5, 2.5], [0.5, 4.5], [0.5, 4.5], [0.5, 2.5], [1.5, 5.5], [0.5, 2.5], [0.5, 3.5])^{T}.

One experiment had c = (92, [163.8, 200.2], 28, 30, 172, 30, 90, 50). (That is, only c_{2} was an interval, for this particular experiment.)

The integer conditions on variables x_{2} through x_{9} were enforced by introducing equality constraints of the form

sin(πx_{i}) = 0, 2 ≤ i ≤ 9.

(This technique was first suggested in [3].) As a consequence of the special form of these constraints, it was logical to set ϵ_{c} = 0.5 (a value that assures accuracy to within the nearest integer). Also, η_{ϕ} and η_{c} were both set to 0.5.

It was not possible for GlobSol to complete this problem within several days of processor time without the thickness stopping conditions described above. It also took an inordinate amount of time when c_{2} was not subdivided. When c_{2} was subdivided into varying numbers of intervals, the results were as in the following table.

5 c_{2} subdivisions 10 c_{2} subdivisions 20 c_{2} subdivisions
# answer boxes  411  194  82
# boxes tried  18,994  28,996  50,073
CPU sec.  42,310  19,010  27,090
objective range  [317, 327]  [317, 326]  [317, 325]
x_{1} range  [10.9, 14.4]  [11.1, 14.3]  [11.3, 14.3]
x_{2} range  1  1  1
x_{3} range  1  1  1
x_{4} range  [2, 4]  [2, 3]  [2, 3]
x_{5} range  [2, 4]  [2, 4]  [2, 3]
x_{6} range  1  1  1
x_{7} range  [3, 5]  [3, 5]  [3, 5]
x_{8} range  1  1  1
x_{9} range  [1, 2]  [1, 2]  [1, 2]

Also of interest in this problem is the range on the objective ϕ associated with particular integer values of x_{2} through x_{9}. This information is also available from the GlobSol output boxes.

4. SUMMARY

A criterion for determining when to quit subdividing interval bounds into subregions has been introduced. This criterion is based on comparing the width of interval evaluations at points with interval evaluations over the entire bounds. The criterion enables verified global optimization algorithms to handle problems with significant uncertainty in the original data. Inclusion of the criterion into a verified global optimization code has been studied, and experiments on a significant problem have been carried out.

5. ACKNOWLEDGEMENT

I wish to thank G. W. Walster for suggesting that I study the issue of thick constants, and encouraging that I do so.

REFERENCES

[1] G. Alefeld and J. Herzberger. Introduction to Interval Computations. Academic Press, New York, 1983.

[2] G. F. Corliss. Globsol entry page, 1998. http://www.mscs.mu.edu/ globsol/.

[3] E. R. Hansen. Global Optimization Using Interval Analysis. Marcel Dekker, Inc., New York, 1992.

[4] R. B. Kearfott. Algorithm 763: INTERVAL ARITHMETIC: A Fortran 90 module for an inerval data type. ACM Trans. Math. Software, 22(4):385–392, December 1996.

[5] R. B. Kearfott. Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht, Netherlands, 1996.

[6] R. B. Kearfott, M. Dawande, K.-S. Du, and C.-Y. Hu. Algorithm 737: INTLIB, a portable FORTRAN 77 interval standard function library. ACM Trans. Math. Software, 20(4):447–459, December 1994.

[7] R. E. Moore. Methods and Applications of Interval Analysis. SIAM, Philadelphia, 1979.

[8] A. Neumaier. Interval Methods for Systems of Equations. Cambridge University Press, Cambridge, England, 1990.

[9] D. Ratz and T. Csendes. On the selection of subdivision directions in interval branch-and-bound methods for global optimization. Journal of Global Optimization, 7:183–207, 1995.

[10] C. M. Rocco, A. J. Miller, and R. B. Kearfott. Uncertainty analysis of indirect-grouping maintenance strategy using an interval arithmetic approach, 1999. Preprint, C/O POBA INTERNATIONAL No. 100, P. O. Box 02-5255, Universidad Central de Venezuela, Miami, Fl 33102.