KBK 2004-06-14: Last week, I found myself dealing with a fairly complex numerical problem that involved minimising a function f(u) with respect to a parameter u. f(u) was evaluated by an external C program; each run took over a minute to complete. I was therefore in need of a procedure that could vary u within bounds determined in advance and minimise f(u). The expense of running the external program meant that I wanted to do this in very few steps. I therefore wound up implementing Brent's method, a method that converges very rapidly. In case someone else finds it useful, I package it up here.
The method is based on parabolic interpolation. If f(x)=fx, f(v)=fv, and f(w)=fw are the three smallest values of the function found so far, it fits a parabola to them and finds the point where the derivative of that parabola vanishes; that point becomes the next place to evaluate the function.
The method has a fair amount of bookkeeping to make sure that parabolic interpolation is actually right. It maintains an interval [a,b] where the minimum is known to lie, and rejects the parabolic interpolation if it would give rise to a point outside the interval. It also keeps track of the size of the steps taken, and rejects parabolic interpolation if two consecutive steps fail to shrink the interval enough.
If parabolic interpolation is rejected, the interval where the minimum is known to be is subdivided in the Golden Ratio, and the point where it is split is used as the next candidate point. The latter process is guaranteed to converge (because it always brackets the minimum).
The algorithm is due to R. P. Brent (Algorithms for Minimization Without Derivatives. Englewood Cliffs, New Jersey: Prentice Hall, 1975) - see Chapter 5 for the details.
#----------------------------------------------------------------------
#
# min_bound_1d --
#
# Find a local minimum of a function between two given
# abscissae. Derivative of f is not required.
#
# Usage:
# min_bound_1d f x1 x2 ?-option value?,,,
#
# Parameters:
# f - Function to minimize. Must be expressed as a Tcl
# command, to which will be appended the value at which
# to evaluate the function.
# x1 - Lower bound of the interval in which to search for a
# minimum
# x2 - Upper bound of the interval in which to search for a minimum
#
# Options:
# -relerror value
# Gives the tolerance desired for the returned
# abscissa. Default is 1.0e-7. Should never be less
# than the square root of the machine precision.
# -maxiter n
# Constrains minimize_bound_1d to evaluate the function
# no more than n times. Default is 100. If convergence
# is not achieved after the specified number of iterations,
# an error is thrown.
# -guess value
# Gives a point between x1 and x2 that is an initial guess
# for the minimum. f(guess) must be at most f(x1) or
# f(x2).
# -abserror value
# Gives the desired absolute error for the returned
# abscissa. Default is 1.0e-10.
# -trace boolean
# A true value causes a trace to the standard output
# of the function evaluations.
#
# Results:
# Returns a two-element list comprising the abscissa at which
# the function reaches a local minimum within the interval,
# and the value of the function at that point.
#
# Side effects:
# Whatever side effects arise from evaluating the given function.
#
#----------------------------------------------------------------------
proc min_bound_1d { f x1 x2 args } {
set phim1 0.6180339887498949
set twomphi 0.3819660112501051
array set params {
-relerror 1.0e-7
-abserror 1.0e-10
-maxiter 100
-trace 0
}
set params(-guess) [expr { $phim1 * $x1 + $twomphi * $x2 }]
if { ( [llength $args] % 2 ) != 0 } {
return -code error -errorcode [list min_bound_1d wrongNumArgs] \
"wrong \# args, should be\
\"[lreplace [info level 0] 1 end f x1 x2 ?-option value...?]\""
}
foreach { key value } $args {
if { ![info exists params($key)] } {
return -code error -errorcode [list min_bound_1d badoption $key] \
"unknown option \"$key\",\
should be -abserror, -initial, -maxiter, -relerror,\
or -trace"
}
set params($key) $value
}
# a and b presumably bracket the minimum of the function. Make sure
# they're in ascending order.
if { $x1 < $x2 } {
set a $x1; set b $x2
} else {
set b $x1; set a $x2
}
set x $params(-guess); # Best abscissa found so far
set w $x; # Second best abscissa found so far
set v $x; # Most recent earlier value of w
set e 0.0; # Distance moved on the step before
# last.
# Evaluate the function at the initial point
set s $f; lappend s $x; set fx [uplevel 1 $s]
if { $params(-trace) } {
puts stdout [format {x=%.17g f(x)=%.17g (initial evaluation)} $x $fx]
}
set fw $fx
set fv $fx
for { set iter 0 } { $iter < $params(-maxiter) } { incr iter } {
# Find the midpoint of the current interval
set xm [expr { 0.5 * ( $a + $b ) }]
# Compute the current tolerance for x, and twice its value
set tol [expr { $params(-relerror) * abs($x) + $params(-abserror) }]
set tol2 [expr { $tol + $tol }]
if { abs( $x - $xm ) <= $tol2 - 0.5 * ($b - $a) } {
return [list $x $fx]
}
set golden 1
if { abs($e) > $tol } {
# Use parabolic interpolation to find a minimum determined
# by the evaluations at x, v, and w. The size of the step
# to take will be $p/$q.
set r [expr { ( $x - $w ) * ( $fx - $fv ) }]
set q [expr { ( $x - $v ) * ( $fx - $fw ) }]
set p [expr { ( $x - $v ) * $q - ( $x - $w ) * $r }]
set q [expr { 2. * ( $q - $r ) }]
if { $q > 0 } {
set p [expr { - $p }]
} else {
set q [expr { - $q }]
}
set olde $e
set e $d
# Test if parabolic interpolation results in less than half
# the movement of the step two steps ago.
if { abs($p) < abs( .5 * $q * $olde )
&& $p > $q * ( $a - $x )
&& $p < $q * ( $b - $x ) } {
set d [expr { $p / $q }]
set u [expr { $x + $d }]
if { ( $u - $a ) < $tol2 || ( $b - $u ) < $tol2 } {
if { $xm-$x < 0 } {
set d [expr { - $tol }]
} else {
set d $tol
}
}
set golden 0
}
}
# If parabolic interpolation didn't come up with an acceptable
# result, use Golden Section instead.
if { $golden } {
if { $x >= $xm } {
set e [expr { $a - $x }]
} else {
set e [expr { $b - $x }]
}
set d [expr { $twomphi * $e }]
}
# At this point, d is the size of the step to take. Make sure
# that it's at least $tol.
if { abs($d) >= $tol } {
set u [expr { $x + $d }]
} elseif { $d < 0 } {
set u [expr { $x - $tol }]
} else {
set u [expr { $x + $tol }]
}
# Evaluate the function
set s $f; lappend s $u; set fu [uplevel 1 $s]
if { $params(-trace) } {
puts stdout [format {x=%.17g f(x)=%.17g (%s)} $u $fu \
[expr { $golden ?
"golden section" :
"parabolic interpolation" }]]
}
if { $fu <= $fx } {
# We've the best abscissa so far.
if { $u >= $x } {
set a $x
} else {
set b $x
}
set v $w
set fv $fw
set w $x
set fw $fx
set x $u
set fx $fu
} else {
if { $u < $x } {
set a $u
} else {
set b $u
}
if { $fu <= $fw || $w == $x } {
# We've the second-best abscissa so far
set v $w
set fv $fw
set w $u
set fw $fu
} elseif { $fu <= $fv || $v == $x || $v == $w } {
# We've the third-best so far
set v $u
set fv $fu
}
}
}
return -code error -errorcode [list min_bound_1d noconverge $iter] \
"[lindex [info level 0] 0] failed to converge after $iter steps."
}The following demonstration script gives a feel for how the program works. Note that the easy case - the function is a parabola - finds the solution on the fourth iteration (and takes two more steps to "check its work."
One awkward function, included to show how the method deals with pathological behaviour, is
f(x)=-1/( 0.01 + abs(x-5) ).
The trouble with this function is that the parabolic interpolation that the method uses to accelerate convergence fails with it; the function has an undefined derivative at its minimum, and is concave downward everywhere else. The method keeps attempting to use parabolic interpolation, and every so often backs off and simply subdivides the interval into two parts (in the ration of the Golden Section) and tries again. Even with this struggle, it evaluates the function only 25 times to come within one part in ten million of the correct answer.
if { [info exists argv0] && [string equal $argv0 [info script]] } {
set ::tcl_precision 17
puts "find the minimum of ( (x + 3) * (x - 1) ) -- easy"
proc f x {
expr { ($x + 3) * ($x - 1) }
}
puts [min_bound_1d f -10. 10. -trace 1]
puts ""
puts "find the first minumum of cos(x)"
proc f x { expr { cos($x) } }
puts [min_bound_1d f 0 6.28318 -trace 1]
puts ""
puts "find the minimum of -exp( -(x-3)**2 / 2 )"
proc f x {
set t [expr { $x - 3. }]
set v [expr { - exp ( - $t * $t / 2. ) }]
}
puts [min_bound_1d f 0. 30. -trace 1]
puts ""
puts "find the first minimum of cos(x)/x"
proc f x { expr { cos($x) / $x } }
puts [min_bound_1d f 0 6.28318 -trace 1]
puts ""
puts "find (difficult!) the minimum of ( -1 / ( 0.01 + abs(x-5) ) )"
proc f x { expr { -1.0 / ( 0.01 + abs( $x - 5. ) ) } }
puts [min_bound_1d f 0 20. -trace 1]
}