<- previous    index    next ->

Lecture 17, Optimization, finding minima

A difficult numerical problem is the finding of a global minima
of an unknown function. This type of problem is called an
optimization problem because the global minima is typically
the optimum solution.

In general we are given a numerically computable function of
some number of parameters  v = f(x_1, x_2, ... , x_n)
and must find the values of x_1, x_2, ... , x_n that gives
the smallest value of v. Or, by taking the absolute value,
find the values of x_1, x_2, ... , x_n that give the value
of v closest to zero.

Generally the problem is bounded and there are given maximum and
minimum values for each parameter. There are typically many
places where local minima exists. Thus, the general solution
must include a global search then a local search to find the
local minima. There is no general guaranteed optimal solution.

First consider a case of only one variable on a non differentiable
function, y = f(x) where x has bounds xmin and xmax. There may be
many local minima, valleys that are not the deepest.

  *     *       *      *         *
   *   * *     * *    *  *      * *   *
    * *   *   *   *  *  * *    *   * * 
     *     * *     **      *  *     *
            *                *
                            *

Consider evaluating f(x) at some initial point x0 and x0+dx.
If f(x0) < f(x0+dx) you might move to x0-dx.
if f(x0) > f(x0+dx) you might move to x0+dx+dx.
The above may be very bad choices!

Here are the cases you should consider:
Compute  yl=f(x-dx)  y=f(x)  yh=f(x+dx)  for some dx

      yh      y            yh       y      yl         yl
    y           yh    yl         yl             yh       y
 yl        yl            y            yh      y            yh

 case 1    case 2     case 3     case 4    case 5     case6

For your next three points, always keep best x:
case 1 x=x-dx                possibly dx=2*dx
case 2 x=x-dx  dx=dx/2 
case 3 dx=dx/2
case 4 x=x+dx  dx=dx/2
case 5 dx=dx/2
case 6 x=x+dx                possibly dx=2*dx
Then loop.

There could be a local minima, thus when dx gets small enough,
remember the best x and use another global search value to
look for a better optimum. Some heuristics may be needed to
increase dx. This is one of many possible algorithms.

Another algorithm that is useful for large areas in two
dimensions for z=f(x,y) is:
Use a small dx and dy to evaluate a preferred direction.
Use an expanding search, doubling dx and dy until no more
progress is made. Then use a contracting search, halving dx and dy
to find the local minima on that direction.
Repeat until the dx and dy are small enough.
The numbers indicate a possible order of evaluation of the points
(in one dimension).

   1
     2
        
         3



                                          5
                              6
                         7
                    4
                       8
                      9



The pseudo derivatives are used to find the preferred direction:
(After finding the best case from above, make positive dx and dy best.)
  z=f(x,y)
  zx=f(x+dx,y)                    zx < z
  zy=f(x,y+dy)                    zy < z
  r=sqrt(((z-zx)^2+(z-zy)^2))
  dx=dx*(z-zx)/r
  dy=dy*(z-zy)/r

This method has worked well on the spiral trough.


The really tough problems have many discontinuities.
I demonstrated a function that was everywhere discontinuous.
The function was f(x)=x^2-1 with f(x)=1 if the bottom bit
of x^2-1 is a one.


A sample program that works for some functions of three
floating point parameters is shown below. Then, a more general
program with a variable number of parameters is presented
with a companion crude global search program.
 
Three parameter optimization:

optm3.h
optm3.c
test_optm3.c
test_optm3_c.out

N parameter optimization:

optmn.h
optmn.c
test_optmn.c
test_optmn_c.out


In MatLab use "fminsearch" see the help file.
Each search is from one staring point.
You need nested loops to try many starting points.
I got error warnings that I ignored, OK to leave them in your output.

An interesting test case is a spiral trough:




test_spiral.f90
test_spiral_f90.out
spiral.f90
spiral.c


    <- previous    index    next ->

Other links

Go to top