[Toon-members] TooN/optimization conjugate_gradient.h

[Edward Rosten]

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