[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN/optimization conjugate_gradient.h
From: |
Edward Rosten |
Subject: |
[Toon-members] TooN/optimization conjugate_gradient.h |
Date: |
Wed, 24 Jun 2009 14:11:36 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 09/06/24 14:11:36
Modified files:
optimization : conjugate_gradient.h
Log message:
Prevent conjugate_gradient from going in to an infinite loop if it
actually
reaches the optimum exactly.
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/optimization/conjugate_gradient.h?cvsroot=toon&r1=1.7&r2=1.8
Patches:
Index: conjugate_gradient.h
===================================================================
RCS file: /cvsroot/toon/TooN/optimization/conjugate_gradient.h,v
retrieving revision 1.7
retrieving revision 1.8
diff -u -b -r1.7 -r1.8
--- conjugate_gradient.h 23 Jun 2009 11:49:07 -0000 1.7
+++ conjugate_gradient.h 24 Jun 2009 14:11:36 -0000 1.8
@@ -5,7 +5,6 @@
#include <cstdlib>
namespace TooN{
-
namespace Internal{
@@ -44,10 +43,12 @@
///@param a_val The value of the function at zero.
///@param func Function to bracket
///@param initial_lambda Initial stepsize
+ ///@param zeps Minimum bracket size.
///@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$.
+ /// and <code>m[i][1]</code> contains the corresponding values
of \f$f(x)\f$. If the bracket
+ /// drops below the minimum bracket size, all zeros are returned.
///@ingroup gOptimize
- template<typename Precision, typename Func> Matrix<3,2,Precision>
bracket_minimum_forward(Precision a_val, const Func& func, Precision
initial_lambda=1)
+ template<typename Precision, typename Func> Matrix<3,2,Precision>
bracket_minimum_forward(Precision a_val, const Func& func, Precision
initial_lambda, Precision zeps)
{
//Get a, b, c to bracket a minimum along a line
Precision a, b, c, b_val, c_val;
@@ -93,6 +94,8 @@
if(b_val < a_val)// we have a bracket
break;
+ else if(lambda < zeps)
+ return Zeros;
else //Contract the bracket
{
c = b;
@@ -177,6 +180,8 @@
Precision linesearch_epsilon; ///< Additive term in tolerance to
prevent excessive iterations if \f$x_\mathrm{optimal} = 0\f$. Known as \c ZEPS
in numerical recipies. Defaults to 1e-20
int linesearch_max_iterations; ///< Maximum number of iterations in
the linesearch. Defaults to 100.
+ Precision bracket_epsilon; ///<Minimum size for initial minima
bracketing. Below this, it is assumed that the system has converged. Defaults
to 1e-20.
+
int iterations; ///< Number of iterations performed
///Initialize the ConjugateGradient class with sensible values.
@@ -229,6 +234,8 @@
linesearch_epsilon = 1e-20;
linesearch_max_iterations=100;
+ bracket_epsilon=1e-20;
+
iterations=0;
}
@@ -243,15 +250,16 @@
/// - old_y
/// - iterations
/// Note that the conjugate direction and gradient are not updated.
+ /// If bracket_minimum_forward detects a local maximum, then
essentially a zero
+ /// sized step is taken.
/// @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);
+ Matrix<3,2,Precision> bracket =
Internal::bracket_minimum_forward(y, line, bracket_initial_lambda,
bracket_epsilon);
double a = bracket[0][0];
double b = bracket[1][0];
@@ -261,6 +269,14 @@
double b_val = bracket[1][1];
double c_val = bracket[2][1];
+ old_y = y;
+ old_x = x;
+ iterations++;
+
+ //Local maximum achieved!
+ if(a==0 && b== 0 && c == 0)
+ return;
+
//We should have a bracket here
assert(a < b && b < c);
assert(a_val > b_val && b_val < c_val);
@@ -272,13 +288,8 @@
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