[Matrix-SIG] Nonlinear optimization routines anyone?

[Ryszard Czerminski]

On 15 Mar 1999, Janne Sinkkonen wrote:

> David Ascher <da@ski.org> writes:
> 
> > I suspect that if I knew enough about optimization I could code up
> > the subset that I need, but I'm not sure I'd trust my own code to do
> > this...
> 
> Coding CG or some such with a simple line search is easy enough, but
> getting the line search robust may be more difficult. And robust stuff
> is hard to code also because it may "work" (although suboptimally) even
> with grave errors in the code. :)
> 
> -- 
> Janne

I have coded this up few month ago (based on NR).
It is simple line search (without derivatives though).
It seems to work OK.

Ryszard

def SIGN(a,b):
    if b < 0: return -abs(a)
    else    : return  abs(a)

def SWAP(a,b): return b,a

def mnbrak(ax,bx,cx,fa,fb,fc,func):
#
# brakets minimum of 1D function
# i.e. returns a,b,c values such that fb < min(fa,fb)
#
   GOLD = 1.618034; GLIMIT = 100.; TINY = 1.0e-20; dum = 0.

   fa = func(ax); fb = func(bx)

   if fb > fa:
      ax,bx = SWAP(ax,bx)
      fa,fb = SWAP(fa,fb)

   cx = bx + GOLD*(bx-ax)  # first guess for c
   fc = func(cx)

   #print 'mnbrak> a,b,c,fa,fb,fc = ',ax,bx,cx,fa,fb,fc

   #iter = 0
   while fb > fc:
      #iter = iter + 1
      #print 'mnbrak> iter,a,b,c,fa,fb,fc = ',iter,ax,bx,cx,fa,fb,fc
      r = (bx-ax)*(fb-fc)
      q = (bx-cx)*(fb-fa)
      u = bx - ((bx-cx)*q - (bx-ax)*r) / (2*SIGN(max(abs(q-r),TINY),q-r))
      ulim = bx + GLIMIT*(cx-bx)

      if (bx-u)*(u-cx) > 0:
         fu = func(u)
         if fu < fc:
            ax = bx; bx = u
            fa = fb; fb = fu
            return ax,bx,cx,fa,fb,fc
         elif fu > fb:
            cx = u; fc = fu
            return ax,bx,cx,fa,fb,fc

         u = cx + GOLD*(cx-bx)
         fu = func(u)

      elif (cx-u)*(u-ulim) > 0:
         fu = func(u)
         if fu < fc:
            bx = cx; cx = u; u = cx + GOLD*(cx-bx)
            fb = fc; fc = fu; fu = func(u)
            
      else:
         u = cx + GOLD*(cx-bx)
         fu = func(u)

      ax = bx; bx = cx; cx = u
      fa = fb; fb = fc; fc = fu

   return ax,bx,cx,fa,fb,fc

def brent(ax,bx,cx,f,tol,maxiter):
#
# minimizes 1D function
#
   ZEPS = 1.0e-10
   CGOLD = 0.381966

   a = min(ax,cx)
   b = max(ax,cx)
   x = w = v = bx
   fx = fw = fv = f(x)

   e = 0.

   iter = 0

   while iter < maxiter:
      iter = iter + 1
      xm = (a+b)/2
      tol1 = tol*abs(x) + ZEPS
      tol2 = 2*tol1

      if abs(x-xm) < (tol2-(b-a)/2): # test for convergence
         return x, fx

      if abs(e) > tol:
         r = (x-w)*(fx-fv)
         q = (x-v)*(fx-fw)
         p = (x-v)*q - (x-w)*r
         q = 2*(q-r)

         if q > 0: p = -p

         q = abs(q)
         etemp = e
         e = d

         if abs(p) >= abs(q*etemp/2) or p <= q*(a-x) or p >= q*(b-x):
            if x >= xm: e = a-x
            else      : e = b-x
            d = CGOLD*e
         else:
            d = p/q
            u = x + d
            if u-a < tol2 or b-u < tol2: d = SIGN(tol1,xm-x)

      else:
         if x >= xm: e = a-x
         else      : e = b-x
         d = CGOLD*e

      if abs(d) >= tol1: u = x + d
      else             : u = x + SIGN(tol1,d)
      fu = f(u)

      if fu < fx:

         if u >= x : a = x
         else      : b = x

         v  = w;  w =  x;  x =  u
         fv = w; fw = fx; fx = fu

      else:

         if u < x : a = u
         else     : b = u

         if fu <= fw or w == x:
            v = w; w = u; fv = fw; fw = fu
         elif fu <= fv or v == x or v == w:
            v = u; fv = fu

   print 'brent> error: iter, maxiter = ',iter, maxiter
   return x, fx