The Mathematica Journal
Feature Articles
New Products
New Publications
News Bulletins
Write Us
About the Journal
Staff and Contributors
Back Issues
Download this Issue

Global Minima

Minimization is an important numerical technique. Mathematica's built-in FindMinimum[expr,{x,x0},{y,y0},...] implements a typical method that uses the derivatives of the function expr to follow the steepest descent of the surface defined by expr to arrive at a local minimum. This local minimum need not be the global minimum, however.

Global techniques need to search the whole valid range in an attempt to find the global minimum or minima. Methods that evaluate expr for a number of sample values do not work well with rapidly varying functions and might miss singular points just as sampled function plots may miss such phenomena.

Interval methods work in much the same way as the interval plotting method described in the previous sections. Adaptive subdivision of the domain is used to zoom into the places where the global minima could lie.

An efficient algorithm must try to exclude as many regions of the domain as early as possible. If the domain has already been subdivided into a collection of (hyper) rectangles, we apply the function to each of them, giving us a collection of range intervals. Define minmax as the minimum of the maxima of these intervals. The global minimum must lie in a rectangle whose corresponding range intersects minmax . Because there is at least one range interval whose maximum is equal to minmax we know that there exists a value of the object function that is less than or equal to minmax. Any range interval lying completely above minmax cannot contain a smaller value.

There are two obvious ways to choose the next domain rectangle to be subdivided. We could use the one with the smallest maximal value of the object function (which is equal to minmax ) in the hope of finding one subinterval that will decrease this value. However, depending on the function, we might spend a lot of effort zooming into an area that will eventually contain only a local minimum. An alternative is to choose the range with the smallest lower bound of the range interval, hoping that it will contain the global minimum. This method is the one we chose; it leads to the criterion mentioned at the end of the previous section for sorting the priority queue that holds the rectangles.

The stopping criterion is the same as with interval plots: either the domain is small enough, or the range is narrow, because the function is rather flat in this domain. At the end (when the priority queue becomes empty), we have to go once more through the list of final rectangles to remove any that are now no longer minimal, because the value of minmax could have changed since the rectangle was put on this list. We return a list [Graphics:../Images/index_gr_37.gif], giving the interval containing the minimal value of the object function and the list of domains where it could lie. Listing 4 shows the code, which is also part of Intervals.m.

IntervalMinimum[expr_, ranges:{_Symbol,_,_}.., opts___?OptionQ] :=
    With[{vars = First/@{ranges},n = Length[{ranges}],
          tol = Tolerance /. {opts} /. Options[IntervalMinimum],
          dom = N[Interval /@ Rest /@ {ranges}]},
    With[{f = Function@@Hold[vars, expr]},
    Module[{div = MakeQueue[Min[#1[[2]]] < Min[#2[[2]]]&], fs,
            minmax = Infinity, fin = h[], new},
      EnQueue[div, {dom, f@@dom}];
        new = DeQueue[div];   (* get smallest one *)
        If[ Min[new[[2]]] >= minmax,
            Continue[] ]; (* no longer in the game *)
        (* split it *)
        new = split[new[[1]]];
        fs = interval[Apply[f, new, {1}]];
        new = Transpose[{new, fs}];
        (* update minmax *)
        minmax = Min[Max /@ fs, minmax];
        (* stopping criteria *)
         If[ Min[#[[2]]] < minmax, (* requeue *)
          If[ length[#[[2]]] < tol || length[#[[1,1]]] < tol,
            fin = h[#, fin], (* no further division *)
            EnQueue[div, #]  (* divide further *)
         ]]&, new ];
      fin = List@@Flatten[fin, Infinity, h];
      (* throw away non-minima *)
      fin = Select[fin, Min[#[[2]]]<minmax&];
      {new, fs} = Transpose[fin];
      {Interval[{Min[fs], minmax}], new}

Listing 4: Global Minimization.

This simple example with a minimum of 0 at 0 shows how the domain is subdivided and a number of intervals are returned.


Here is a simple two-dimensional example. The number of rectangles is too large to show them in full.


However, we can show them in a picture. The auxiliary function RectangleGraph produces a graphic of a list of rectangles. It was also used for interval plotting.



In the next example, we get seemingly trivial domains. Note that the minimum lies at the border of the domain. This is a case that cannot be handled well with methods based on steepest descent.


The input form shows that the numbers differ.

Interval[{0.09999999999999953, 0.10000000020954718}]

Here is a case where the minimum is equal to [Graphics:../Images/index_gr_48.gif]:


The following function has one-dimensional sets of minima; consequently, we get a large number of intervals back.


The graph shows them clearly.



The function [Graphics:../Images/index_gr_55.gif] has infinitely many disjoint minima. At some point near zero they will no longer be separated by the intervals.


The union of all domain intervals shows the isolating intervals found. The first one contains infinitely many minima.


Culioli [2] mentions a family of functions with minima that are difficult to find using steepest descent methods. The particular behavior of the derivatives also makes them hard test cases for interval methods, which rely on Taylor series and thus also on derivatives to find enclosing rectangles.

Here is a four-dimensional test case.


We have not made much progress; the range where the minimum could be is still rather large. The enclosing rectangle of the domains leads to the same conclusion.


Decreasing the tolerance further will quickly lead to an explosive growth of the number of intervals generated and Mathematica might well run out of memory even on a larger computer.

Converted by Mathematica      May 8, 2000

[The Mathematica Programmer Index] [Prev Page][Next Page]