10 YEARS OF MATHEMATICA

Interval Methods for
Global Optimization

Roman E. Maeder

MathConsult Dr. R. Mäder
Samstagernstrasse 58a
8832 Wollerau
Switzerland
mailto:maeder@mathconsult.ch
http://www.mathconsult.ch/

The use of interval arithmetic for constrained and unconstrained global optimization. Interval methods allow the localization of global minima even in the presence of discontinuities and singularities. Contraints can be nonlinear.

A presentation given at the 1998 Worldwide Mathematica Conference, Chicago, June 1998.

© 1998 MathConsult Dr. R. Mäder. All rights reserved.

Interval Arithmetic

To appear in The Mathematica Journal, Vol. 7, issue 2.

Guaranteed Inclusion

For an interval I, J = f(I) is guaranteed to include the images f(x) for all x ∈ I.

Interval arithmetic is closely related to Mathematica's ordinary arbitrary-precision arithmetic but its emphasis is on obtaining global results, not on high-precision computation.

Sin[Interval[{0.1,  1.4}]]

Sin[Interval[{1.2,  1.8}]]

Tan[Interval[{1, 2}]]

Function[x, Sin[x] Cos[x]][Interval[{0, Pi}]]

Function[x, Sin[2  x]/2][Interval[{0, Pi}]]

Interval Plots

Because of the inclusions principle, we know that the graph of f lies in the rectangle given by I and J = f(I). For an interval plot, we refine I until either I or J are small enough (adaptive refinement).

Example sin((1/x)).

f[x_] := Sin[1/x]

Necessary auxiliary functions:

rect[{xi_, yi_}] := 
      Rectangle[{Min[xi], Min[yi]}, {Max[xi], Max[yi]}]

\(split[i_Interval] := 
\ \ 
      With[{m = \((Min[i] + Max[i])\)/2}, 
        \ {Interval[{Min[i], m}], Interval[{m, Max[i]}]}\ ]\)

Init:

i = Interval[{0.02, 0.2}]

fp = Plot[f[x], {x, Min[i], Max[i]}, DisplayFunction -> Identity];

The list of domain intervals

dom = {i};

Loop

Apply f to all domain intervals

vals = Transpose[{dom, f/@dom}]

Construct the rectangles

rs = rect/@vals;

Display them together with an ordinary plot of f

Show[Graphics[{RGBColor[0.5, 0.7, 1], rs}], fp];

Split the domain intervals in half

dom = Flatten[split/@dom]

And repeat...

Further Examples of Interval Plots

Needs["Intervals`"]

IntervalPlot[Sin[1/x], {x, -0.1, 0.1}, Tolerance -> 0.002];

IntervalPlot[1/x, {x, -3.2, 3.21}];

Plot[x  +  Abs[x - 1]^-0.01, {x, 0, 2}];

IntervalPlot[x  +  Abs[x - 1]^-0.01, {x, 0, 2}];

Priority Queues

A priority queue is a container for linearly ordered data that supports these operations

A priority queue is often implemented using a heap.
A heap is a binary tree where the item in each node is larger than the items in the left and right subtrees. This condition guarantees that the largest of all items is in the root of the tree.
Note that the condition does not require the tree to be completely sorted; therefore, we can hope to find more efficient implementations of the heap operations than are necessary for complete sorting.

A suitable implementation keeps the binary tree in an array such that the root is at position 1 and the left and right subtrees of the node at position i are at positions 2  i and 2  i + 1, respectively. The number of items currently in the heap is kept in the variable n.

To insert an item, first put it at position i = n + 1, then compare it with the item at position j = ⌊i/2⌋ and exchange the two if necessary. If an exchange was necessary, set i = j and continue until no exchange is performed or j < 1; the latter condition happens if the new elements turns out to be larger than all other items in the heap.

The largest item is always at position 1. To remove it, put the last item at position n in its place and then compare it with the larger of the two descendants, exchanging it if necessary until it has found its place.

All you ever wanted to know about heaps can, of course, be found in Knuth's The Art of Computer Programming, Vol 3.

The Code

An Example

Needs["PriorityQueue`"]

We create a new queue that sorts its elements into the default decreasing order.

q  =  MakeQueue[]

Now, we insert one item, which is, of course, the largest one.

EnQueue[q,  17]

Using Scan[], we can insert the elements of a list one by one into the queue.

Scan[EnQueue[q,  #]&,  {1,  21,  6,  8}];  q

We retrieve and remove the largest element.

DeQueue[q]

The next largst element has risen to the top of the heap.

q

In a While loop we can retrieve the elements one by one, in order (so, a heap could be used for sorting).

While[! EmptyQueue[q], Print[DeQueue[q]]];

For minimization, we will put data of the form {domain, res}, where domain is a list of intervals and res is a single interval, into the queue, sorted by the smallest value in res. Therefore, our predicate will be Min[#1[[2]]] < Min[#2[[2]]]&.

Global Minimization

Introduction

Local minima

Examples: FindMinimum.

Global techniques

Examples:
Global Optimization by Loehle Enterprises (MathLink interface to external optimizer)
LGO by János Pintér

Interval methods

Algorithm

The domain rectangle is subdivided into (2^d) intervals (d is the dimension of the problem) until either

The Program

Examples

Needs["Intervals`"]

A simple isolated minimum at the origin.

IntervalMinimum[x^2  +  y^2, {x, -1, 1}, {y, -1, 1}, 
      Tolerance -> 10^-6]

Show[RectangleGraphColor[%[[2]]], AspectRatio -> Automatic];

The minimum is on the border of the domain.

IntervalMinimum[x^2  +  y^2, {x, 0.1, 1}, {y, 0.1, 1}, 
      Tolerance -> 10^-6]

Squashing the function leads to a non-square collection of intervals.

IntervalMinimum[0.1  x^2  +  y^2, {x, -1, 1}, {y, -1, 1}, 
      Tolerance -> 0.1]

Show[RectangleGraphColor[%[[2]]], AspectRatio -> Automatic];

Here, the minimum is -∞, lying at the border. In one-dimensional problems, there is no combinatorial explosion, allowing us to use a rather small value for the tolerance.

IntervalMinimum[Log[x], {x, 0.0, 1}, Tolerance -> 10^-10]

This function has an isolated minimum at the origin and a one-dimensional set of minima on a circle.

IntervalMinimum[
      -Cos[Sqrt[x^2 + y^2]], {x, -7, 7}, {y, -7, 7}, 
      Tolerance -> 10^-2];

Show[RectangleGraph[ %[[2]]], AspectRatio -> Automatic];

The example from the section on interval plotting.

\(min = IntervalMinimum[Sin[1/x], {x, 0, 2}, Tolerance -> 10^\(-4\)]; 
    
{min[\([1]\)], Length[min[\([2]\)]]}\)

Show[RectangleGraphColor[ min[[2]]], AspectRatio -> 0.2];

An example proposed by János Pintér.

min = 
      IntervalMinimum[
        x^2 + y^2 + Log[1 + 20  Abs[Sin[x^2 + y^2]]] + 
          10  Abs[Sin[5  x^2 + 15  y^2]], {x, -4, 10}, {y, -4, 10}, 
        Tolerance -> 10^-3]

Plot3D[
      x^2 + y^2 + Log[1 + 20  Abs[Sin[x^2 + y^2]]] + 
        10  Abs[Sin[5  x^2 + 15  y^2]], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 80];

Discussion

János D. Pintér, Global Optimization in Action. Kluwer Academic Publishers. ISBN 0-7923-3757-3.

Reiner Horst and Panos Pardalos, Handbook of Global Optimization. Kluwer Academic Publishers, 1995.

Constrained Minima

Method

Note that the domain borders are a simple form of rectangular constraints.

An arbitrary constraint is a predicate that tells whether a point of the domain is feasible. When applied to a rectangle, the predicate can return three possible results

The main problem is finding a feasible point that allows discarding rectangles known not to contain the minimum. As slong as we don't know a feasible point, we cannot discard any rectangles.

Examples

A simple constraint.

IntervalMinimumConstrained[x,  {x, 0, 1}, Abs[x] > 1/2]

Search constrained to an annular region:

min = 
      IntervalMinimumConstrained[Sin[x y]^2, 
         {x, -3, 3}, {y, -3, 3}, 0.2 < x^2 + y^2 < 7, 
        Tolerance -> 10^-2];

{First[min],  Length[min[[2]]]}

Show[RectangleGraphColor[min[[2]]], AspectRatio -> Automatic];

Equality constraints do not work well with interval methods. No rectangle will ever be completely inside the feasible region.

min = 
      IntervalMinimumConstrained[2  x + y,  {x, -3, 3}, {y, -3, 3}, 
        x^2 + y^2 == 1, Tolerance -> 0.02];

{First[min],  Length[min[[2]]]}

The presence of ∞ alerts to the fact that no feasible point was found.

Show[RectangleGraphColor[min[[2]]], AspectRatio -> Automatic];

If we know a feasible point, we might as well tell the program. Not that it helps much in this case...

min = 
      IntervalMinimumConstrained[2  x + y,  {x, -3, 3}, {y, -3, 3}, 
        x^2 + y^2 == 1, Tolerance -> 0.02,  FeasiblePoint -> {0, 1}];

Conclusions

Please look for this notebook and the packages
Intervals.m
PriorityQueue.m
on my web site, http://www.mathconsult.ch/math/