# [Matrix-SIG] Nonlinear optimization routines anyone?

Ryszard Czerminski ryszard@moldyn.com
Tue, 16 Mar 1999 08:09:08 -0500 (EST)

```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

```