[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN optimization/golden_section.h optimization...
|
From: |
Edward Rosten |
|
Subject: |
[Toon-members] TooN optimization/golden_section.h optimization... |
|
Date: |
Tue, 07 Apr 2009 14:33:06 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 09/04/07 14:33:06
Modified files:
optimization : golden_section.h
Added files:
optimization : conjugate_gradient.h
test : cg_test.cc cg_view.gnuplot
Log message:
Added conjugate gradient optimization and test code.
Updated golden_section to make it work better with ConjugateGradient.
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/optimization/golden_section.h?cvsroot=toon&r1=1.1&r2=1.2
http://cvs.savannah.gnu.org/viewcvs/TooN/optimization/conjugate_gradient.h?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/test/cg_test.cc?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/test/cg_view.gnuplot?cvsroot=toon&rev=1.1
Patches:
Index: optimization/golden_section.h
===================================================================
RCS file: /cvsroot/toon/TooN/optimization/golden_section.h,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -b -r1.1 -r1.2
--- optimization/golden_section.h 3 Apr 2009 17:21:46 -0000 1.1
+++ optimization/golden_section.h 7 Apr 2009 14:33:06 -0000 1.2
@@ -4,27 +4,47 @@
#include <limits>
#include <cmath>
#include <cstdlib>
+#include <iomanip>
namespace TooN
{
using std::numeric_limits;
+
+ /// golden_section_search performs a golden section search line
minimization
+ /// on the functor provided. The inputs a, b, c must bracket the
minimum, and
+ /// must be in order, so that \f$ a < b < c \f$ and \f$ f(a) > f(b) <
f(c) \f$.
+ /// @param a The most negative point along the line.
+ /// @param b The central point.
+ /// @param c The most positive point along the line.
+ /// @param func The functor to minimize
+ /// @param maxiterations Maximum number of iterations
+ /// @param tolerance Tolerance at which the search should be stopped.
+ /// @return The minima position is returned as the first element of the
vector,
+ /// and the minimal value as the second element.
+ template<class Functor, class Precision> Vector<2, Precision>
golden_section_search(Precision a, Precision b, Precision c, const Functor&
func, int maxiterations, Precision tol =
sqrt(numeric_limits<Precision>::epsilon()))
+ {
+ return golden_section_search(a, b, c, func(b), func,
maxiterations, tol);
+ }
+
/// golden_section_search performs a golden section search line
minimization
/// on the functor provided. The inputs a, b, c must bracket the
minimum, and
/// must be in order, so that \f$ a < b < c \f$ and \f$ f(a) > f(b) <
f(c) \f$.
/// @param a The most negative point along the line.
/// @param b The central point.
+ /// @param fb The value of the function at the central point (\f$b\f$).
/// @param c The most positive point along the line.
/// @param func The functor to minimize
- /// @param tolerance
+ /// @param maxiterations Maximum number of iterations
+ /// @param tolerance Tolerance at which the search should be stopped.
/// @return The minima position is returned as the first element of the
vector,
/// and the minimal value as the second element.
- template<class Functor, class Precision> Vector<2, Precision>
golden_section_search(Precision a, Precision b, Precision c, const Functor&
func, Precision tol = sqrt(numeric_limits<Precision>::epsilon()))
+ template<class Functor, class Precision> Vector<2, Precision>
golden_section_search(Precision a, Precision b, Precision c, Precision fb,
const Functor& func, int maxiterations, Precision tol =
sqrt(numeric_limits<Precision>::epsilon()))
{
using std::abs;
//The golden ratio:
const Precision g = (3.0 - sqrt(5))/2;
- Precision x1, x2;
+ Precision x1, x2, fx1, fx2;
//Perform an initial iteration, to get a 4 point
//bracketing. This is rather more convenient than
@@ -33,20 +53,24 @@
{
x1 = b - g*(b-a);
x2 = b;
+
+ fx1 = func(x1);
+ fx2 = fb;
}
else
{
x2 = b + g * (c-b);
x1 = b;
+
+ fx1 = fb;
+ fx2 = func(x2);
}
//We now have an ordered list of points a x1 x2 c
- Precision fx1 = func(x1);
- Precision fx2 = func(x2);
-
//Termination condition from NR in C
- while(abs(c-a) > tol * (abs(x2)+abs(x1)))
+ int itnum = 1; //We've already done one iteration.
+ while(abs(c-a) > tol * (abs(x2)+abs(x1)) && itnum <
maxiterations)
{
if(fx1 > fx2)
{
Index: optimization/conjugate_gradient.h
===================================================================
RCS file: optimization/conjugate_gradient.h
diff -N optimization/conjugate_gradient.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ optimization/conjugate_gradient.h 7 Apr 2009 14:33:06 -0000 1.1
@@ -0,0 +1,303 @@
+#include <TooN/optimization/golden_section.h>
+#include <utility>
+#include <cmath>
+#include <cassert>
+#include <cstdlib>
+
+namespace TooN{
+
+ namespace Internal{
+
+
+ ///Turn a multidimensional function in to a 1D function by specifying a
+ ///point and direction.
+ template<int Size, typename Precision, typename Func> struct LineSearch
+ {
+ const Vector<Size, Precision>& start, direction;
+ const Func& f;
+
+ ///@param s Start point.
+ ///@param d direction
+ ///@param func Function.
+ LineSearch(const Vector<Size, Precision>& s, const Vector<Size,
Precision>& d, const Func& func)
+ :start(s),direction(d),f(func)
+ {}
+
+ ///@param x Position to evaluate function
+ ///@return \f$f(\vec{d} + x\vec{s})\f$
+ Precision operator()(Precision x) const
+ {
+ return f(start + x * direction);
+ }
+ };
+
+ ///Bracket a 1D function by searching forward from zero. The assumption
+ ///is that a minima exists in \f$f(x),\ x>0\f$, and this function
searches
+ ///for a bracket using exponentially growning or shrinking steps.
+ ///@param a_val The value of the function at zero.
+ ///@param func Function to bracket
+ ///@param initial_lambda Initial stepsize
+ ///@return <code>m[i][0]</code> contains the values of \f$x\f$ for the
bracket, in increasing order,
+ /// and <code>m[i][1]</code> contains the corresponding values
of \f$f(x)\f$.
+ template<typename Precision, typename Func> Matrix<3,2,Precision>
bracket_minimum_forward(Precision a_val, const Func& func, Precision
initial_lambda=1)
+ {
+ //Get a, b, c to bracket a minimum along a line
+ Precision a, b, c, b_val, c_val;
+
+ a=0;
+
+ //Search forward in steps of lambda
+ Precision lambda=initial_lambda;
+ b = lambda;
+ b_val = func(b);
+
+ if(b_val < a_val) //We've gone downhill, so keep searching
until we go back up
+ {
+ for(;;)
+ {
+ lambda *= 2;
+ c = lambda;
+ c_val = func(c);
+
+ if(c_val > b_val) // we have a bracket
+ break;
+ else
+ {
+ a = b;
+ a_val = b_val;
+ b=c;
+ b_val=c_val;
+
+ }
+ }
+ }
+ else //We've overshot the minimum, so back up
+ {
+ c = b;
+ c_val = b_val;
+ //Here, c_val > a_val
+
+ for(;;)
+ {
+ lambda *= .5;
+ b = lambda;
+ b_val = func(b);
+
+ if(b_val < a_val)// we have a bracket
+ break;
+ else //Contract the bracket
+ {
+ c = b;
+ c_val = b_val;
+ }
+ }
+ }
+
+ Matrix<3,2> ret;
+ ret[0] = makeVector(a, a_val);
+ ret[1] = makeVector(b, b_val);
+ ret[2] = makeVector(c, c_val);
+
+ return ret;
+ }
+
+}
+
+
+/** This class provides a nonlinear conjugate-gradient optimizer. The following
+code snippet will perform an optimization on the Rosenbrock Bananna function in
+two dimensions:
+
address@hidden
+double Rosenbrock(const Vector<2>& v)
+{
+ return sq(1 - v[0]) + 100 * sq(v[1] - sq(v[0]));
+}
+
+Vector<2> RosenbrockDerivatives(const Vector<2>& v)
+{
+ double x = v[0];
+ double y = v[1];
+
+ Vector<2> ret;
+ ret[0] = -2+2*x-400*(y-sq(x))*x;
+ ret[1] = 200*y-200*sq(x);
+
+ return ret;
+}
+
+int main()
+{
+ ConjugateGradient<2> cg(makeVector(0,0), Rosenbrock,
RosenbrockDerivatives);
+
+ while(!cg.iterate(Rosenbrock, RosenbrockDerivatives))
+ cout << "y_" << iteration << " = " << cg.y << endl;
+
+ cout << "Optimal value: " << cg.y << endl;
+}
address@hidden
+
+Linesearch is currently performed using golden-section search and conjugate
+vector updates are performed using the Polak-Ribiere equations. There many
+tunable parameters, and the internals are readily accessible, so alternative
+termination conditions etc can easily be substituted. However, ususally these
+will not be necessary.
+*/
+template<int Size, class Precision=double> struct ConjugateGradient
+{
+ const int size; ///< Dimensionality of the space.
+ Vector<Size> g; ///< Gradient vector used by the next call to
iterate()
+ Vector<Size> h; ///< Conjugate vector to be searched along in the
next call to iterate()
+ Vector<Size> old_g; ///< Gradient vector used to compute $h$ in the
last call to iterate()
+ Vector<Size> old_h; ///< Conjugate vector searched along in the last
call to iterate()
+ Vector<Size> x; ///< Current position (best known point)
+ Vector<Size> old_x; ///< Previous best known point (not set at
construction)
+ Precision y; ///< Function at \f$x\f$
+ Precision old_y; ///< Function at old_x
+
+ Precision tolerance; ///< Tolerance used to determine if the
optimization is complete. Defaults to square root of machine precision.
+ Precision epsilon; ///< Additive term in tolerance to prevent
excessive iterations if \f$x_\mathrm{optimal} = 0\f$. Known as \c ZEPS in
numerical recipies.
+ int max_iterations; ///< Maximum number of iterations. Defaults
to \c size\f$*100\f$
+
+ Precision bracket_initial_lambda;///< Initial stepsize used in
bracketing the minimum for the line search. Defaults to 1.
+ Precision linesearch_tolerance; ///< Tolerance used to determine if the
linesearch is complete. Defaults to square root of machine precision.
+ int linesearch_max_iterations; ///< Maximum number of iterations in
the linesearch. Defaults to 100.
+
+ int iterations; ///< Number of iterations performed
+
+ template<class Func, class Deriv> ConjugateGradient(const Vector<Size>&
start, const Func& func, const Deriv& deriv)
+ : size(start.size()),
+ g(size),h(size),old_g(size),old_h(size),x(start),old_x(size)
+ {
+ using std::numeric_limits;
+
+ x = start;
+
+ //Start with the conjugate direction aligned with
+ //the gradient
+ g = deriv(x);
+ h = g;
+
+ y = func(x);
+ old_y = y;
+
+ tolerance = sqrt(numeric_limits<Precision>::epsilon());
+ epsilon = 1e-20;
+ max_iterations = size * 100;
+
+ bracket_initial_lambda = 1;
+
+ linesearch_tolerance =
sqrt(numeric_limits<Precision>::epsilon());
+ linesearch_max_iterations=100;
+
+ iterations=0;
+ }
+
+
+ ///Perform a linesearch from the current point (x) along the current
+ ///conjugate vector (h). /The linesearch does not make use of
derivatives.
+ ///You probably do not want to use /this function. See iterate()
instead.
+ ///This function updates:
+ /// - x
+ /// - old_c
+ /// - y
+ /// - old_y
+ /// - iterations
+ /// Note that the conjugate direction and gradient are not updated.
+ /// @param func Functor returning the function value at a given point.
+ template<class Func> void find_next_point(const Func& func)
+ {
+ Internal::LineSearch<Size, Precision, Func> line(x, -h, func);
+
+ //Always search in the conjugate direction (h)
+
+ //First bracket a minimum.
+ Matrix<3,2,Precision> bracket =
Internal::bracket_minimum_forward(y, line, bracket_initial_lambda);
+
+ double a = bracket[0][0];
+ double b = bracket[1][0];
+ double c = bracket[2][0];
+
+ double a_val = bracket[0][1];
+ double b_val = bracket[1][1];
+ double c_val = bracket[2][1];
+
+ //We should have a bracket here
+ assert(a < b && b < c);
+ assert(a_val > b_val && b_val < c_val);
+
+ //Find the real minimum
+ Vector<2, Precision> m = golden_section_search(a, b, c, b_val,
line, linesearch_max_iterations);
+
+ assert(m[0] >= a && m[0] <= c);
+ assert(m[1] <= b_val);
+
+ //Update the current position and value
+ old_y = y;
+ old_x = x;
+
+ x -= m[0] * h;
+ y = m[1];
+
+ iterations++;
+ }
+
+ ///Check to see it iteration should stop. You probably do not want to
use
+ ///this function. See iterate() instead. This function updates nothing.
+ bool finished()
+ {
+ using std::abs;
+ return iterations > max_iterations || 2*abs(y - old_y) <=
tolerance * (abs(y) + abs(old_y) + epsilon);
+ }
+
+ ///After an iteration, update the gradient and conjugate using the
+ ///Polak-Ribiere equations. /@param deriv Functor to compute
derivatives at
+ ///the specified point. This function updates:
+ ///- g
+ ///- old_g
+ ///- h
+ ///- old_h
+ template<class Deriv> void update_vectors_PR(const Deriv& deriv)
+ {
+ //Update the position, gradient and conjugate directions
+ old_g = g;
+ old_h = h;
+
+ g = deriv(x);
+ //Precision gamma = (g * g - oldg*g)/(oldg * oldg);
+ Precision gamma = (g * g - old_g*g)/(old_g * old_g);
+ h = g + gamma * old_h;
+ }
+
+ ///Use this function to iterate over the optimization. Note that after
+ ///iterate returns false, g, h, old_g and old_h will not have been
+ ///updated.
+ ///This function updates:
+ /// - x
+ /// - old_c
+ /// - y
+ /// - old_y
+ /// - iterations
+ /// - g*
+ /// - old_g*
+ /// - h*
+ /// - old_h*
+ /// *'d variables not updated on the last iteration.
+ ///@param func Functor returning the function value at a given point.
+ ///@param deriv Functor to compute derivatives at the specified point.
+ ///@return Whether to continue.
+ template<class Func, class Deriv> bool iterate(const Func& func, const
Deriv& deriv)
+ {
+ find_next_point(func);
+
+ if(!finished())
+ {
+ update_vectors_PR(deriv);
+ return 1;
+ }
+ else
+ return 0;
+ }
+};
+
+}
Index: test/cg_test.cc
===================================================================
RCS file: test/cg_test.cc
diff -N test/cg_test.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ test/cg_test.cc 7 Apr 2009 14:33:06 -0000 1.1
@@ -0,0 +1,64 @@
+#include <TooN/optimization/golden_section.h>
+#include <TooN/optimization/conjugate_gradient.h>
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include <cstdlib>
+
+using namespace TooN;
+using namespace std;
+
+double sq(double d)
+{
+ return d*d;
+}
+
+double Rosenbrock(const Vector<2>& v)
+{
+ return sq(1 - v[0]) + 100 * sq(v[1] - sq(v[0]));
+}
+
+Vector<2> RosenbrockDerivatives(const Vector<2>& v)
+{
+ double x = v[0];
+ double y = v[1];
+
+ Vector<2> ret;
+ ret[0] = -2+2*x-400*(y-sq(x))*x;
+ ret[1] = 200*y-200*sq(x);
+
+ return ret;
+}
+
+
+double Spiral(const Vector<2>& v)
+{
+ double x = v[0];
+ double y = v[1];
+
+ return sin(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))+2.0*x*x+2.0*y*y;
+}
+
+
+Vector<2> SpiralDerivatives(const Vector<2>& v)
+{
+ double x = v[0];
+ double y = v[1];
+
+ double dx
=2.0*(10.0*cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*x*x*x+10.0*cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*x*y*y-cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*y*sqrt(x*x+y*y)+2.0*x*x*x*sqrt(x*x+y*y)+2.0*x*sqrt(x*x+y*y)*y*y)/sqrt(pow(x*x+y*y,3.0));
+ double dy =
2.0*(10.0*cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*y*x*x+10.0*cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*y*y*y+cos(20.0*sqrt(x*x+y*y)+2.0*atan(y/x))*x*sqrt(x*x+y*y)+2.0*y*sqrt(x*x+y*y)*x*x+2.0*y*y*y*sqrt(x*x+y*y))/sqrt(pow(x*x+y*y,3.0));
+
+ return makeVector(dx, dy);
+}
+
+int main()
+{
+ ConjugateGradient<2> cg(makeVector(1.5, 1.5), Spiral,
SpiralDerivatives);
+ cg.bracket_initial_lambda=1e-7; //Prevent the bracket from jumping over
spiral arms for a prettier display
+
+ cout << cg.x << " " << cg.y << endl;
+ while(cg.iterate(Spiral, SpiralDerivatives))
+ cout << cg.x << " " << cg.y << endl;
+ cout << cg.x << " " << cg.y << endl;
+}
+
Index: test/cg_view.gnuplot
===================================================================
RCS file: test/cg_view.gnuplot
diff -N test/cg_view.gnuplot
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ test/cg_view.gnuplot 7 Apr 2009 14:33:06 -0000 1.1
@@ -0,0 +1,19 @@
+####
+#### Compile cg_test.cc to cg_test
+####
+
+set multiplot
+set nocontour
+set pm3d
+unset surface
+set view 0,0
+set isosample 100
+
+splot [-2:2] [-2:2] sin(20*sqrt(x**2 + y**2) + atan(y/x)*2) + 1*(x**2 + y**2)
+unset pm3d
+set surface
+splot [-2:2] [-2:2] '< ./cg_test' using 1:2:3 with linespoints title
'Conjugate Gradient path'
+unset multiplot
+
+
+
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN optimization/golden_section.h optimization...,
Edward Rosten <=