---

gnn_line_search.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_line_search.c
 *  @brief Line Search Procedures Implementation.
 *
 *  @date   : 15-09-03 21:55, 01-10-03 09:50
 *  @author : Pedro Ortega C. <peortega@dcc.uchile.cl>
 *  Copyright  2003  Pedro Ortega C.
 ****************************************************************************/
/*
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Library General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */



/**
 * @brief Line Search Procedures.
 * @defgroup gnn_line_search_doc Line Search Procedures for Nodes.
 * @ingroup  libgnn_node
 *
 * This module contains a set of line search procedures for a \ref gnn_node.
 * A line search procedure is a one-dimensional minimization algorithm, with
 * tries to find a minimum in a given direction.
 *
 * This module tries to factorize the common features available in trainers
 * that need line search procedures during training. It isn't implemented yet.
 *
 *
 */



/******************************************/
/* Include Files                          */
/******************************************/

#include <math.h>
#include "gnn_utilities.h"
#include "gnn_line_search.h"



/******************************************/
/* Static Declaration                     */
/******************************************/




/******************************************/
/* Public Interface                       */
/******************************************/

/**
 * @brief Bracketing procedure.
 * @ingroup gnn_line_search
 *
 * Given a \ref gnn_line buffer, and given distinct initial points \a ax and
 * \a bx, this routine searches in the downhill direction (defined by the
 * function associated to \a line as evaluated at the initial points) and
 * returns new points \a ax, \a bx, \a cx that bracket a minimum of the
 * function. Also returned are the function values at the three points,
 * \a fa, \a fb, and \a fc.
 *
 * [Adapted from Numerical Recipes in C.]
 *
 * @param line A pointer to a \ref gnn_line buffer.
 * @param s    The index of the first pattern in the minibatch to be evaluated.
 * @param n    The size of the minibatch to be evaluated.
 * @param ax   A pointer to an initial point.
 * @param bx   A pointer to the second initial point.
 * @param cx   A pointer where the middle point should be placed in.
 * @param fa   A pointer where the value \f$<E(a)>\f$ should be placed in.
 * @param fb   A pointer where the value \f$<E(b)>\f$ should be placed in.
 * @param fc   A pointer where the value \f$<E(c)>\f$ should be placed in.
 * @return Returns 0 if succeeded.
 */
int
gnn_line_search_bracket (gnn_line *line, size_t s, size_t n,
                         double *ax, double *bx, double *cx,
                         double *fa, double *fb, double *fc)
{
    double ulim;
    double u;
    double r;
    double q;
    double fu;
    double dum;

    /* evaluate function */
    gnn_line_pats (line, *ax, gnnLineE, s, n);
    *fa = GNN_LINE_E (line);
    gnn_line_pats (line, *bx, gnnLineE, s, n);
    *fb = GNN_LINE_E (line);

    /* swith roles of a and b, so that we can   */
    /* go downhill in the direction from a to b */
    if (*fb > *fa)
    {
        GNN_SHIFT3 (dum, *ax, *bx, dum);
        GNN_SHIFT3 (dum, *fb, *fa, dum);
    }

    /* First guess for c */
    *cx = *bx + GNN_GOLDEN_OLD * (*bx - *ax);
    gnn_line_pats (line, *cx, gnnLineE, s, n);
    *fc = GNN_LINE_E (line);

    /* Keep returning here until we bracket */
    while (*fb > *fc)
    {
        /* Compute u by parabolic extrapolation from  */
        /* a, b, c. GNN_GOLDEN_TINY is use to prevent */
        /* any posible division by zero.              */
        r = (*bx - *ax) * (*fb - *fc);
        q = (*bx - *cx) * (*fb - *fa);
        u = *bx
            - ((*bx - *cx) * q - (*bx - *ax) * r) /
              ( 2.0 * GNN_USE_SIGN ( GNN_MAX (fabs (q - r), GNN_GOLDEN_TINY),
                                     q - r ) );
        ulim = *bx + GNN_GOLDEN_LIMIT * (*cx - *bx);

        /* We won't go farther than this. */
        /* Test various possibilities:    */
        if ((*bx - u) * (u - *cx) > 0.0)
        {
            /* Parabolic u is between b and c: try it. */
            gnn_line_pats (line, u, gnnLineE, s, n);
            fu = GNN_LINE_E (line);

            if (fu < *fc)
            {
                /* Got a minimum between b and c. */
                *ax = *bx;
                *bx =  u;
                *fa = *fb;
                *fb =  fu;
                return 0;
            }
            else if (fu > *fb)
            {
                /* Got a minimum between a und u. */
                *cx = u;
                *fc = fu;
                return 0;
            }

            /* Parabolic fit was no use. Use default magnification. */
            u  = *cx + GNN_GOLDEN_OLD * (*cx - *bx);
            gnn_line_pats (line, u, gnnLineE, s, n);
            fu = GNN_LINE_E (line);
        }
        else if ((*cx - u) * (u - ulim) > 0.0)
        {
            /* Parabolic fit is between c and its allowed limit. */
            gnn_line_pats (line, u, gnnLineE, s, n);
            fu = GNN_LINE_E (line);
            
            if (fu < *fc)
            {
                GNN_SHIFT3 (*bx, *cx, u, *cx + GNN_GOLDEN_OLD * (*cx - *bx));
                GNN_SHIFT2 (*fb, *fc, fu);
                gnn_line_pats (line, u, gnnLineE, s, n);
                fu = GNN_LINE_E (line);
            }
        }
        else if ((u - ulim) * (ulim - *cx) >= 0.0)
        {
            /* Limit parabolic u to maximum allowed value. */
            u  = ulim;
            gnn_line_pats (line, u, gnnLineE, s, n);
            fu = GNN_LINE_E (line);
        }
        else
        {
            /* Reject parabolic u, use default magnification. */
            u  = *cx + GNN_GOLDEN_OLD * (*cx - *bx);
            gnn_line_pats (line, u, gnnLineE, s, n);
            fu = GNN_LINE_E (line);
        }
            
        /* Eliminate oldest point and continue. */
        GNN_SHIFT3 (*ax, *bx, *cx, u);
        GNN_SHIFT3 (*fa, *fb, *fc, fu);
    }

    return 0;
}

/**
 * @brief Bracketing procedure.
 * @ingroup gnn_line_search
 *
 * Given a \ref gnn_line buffer \a line, and given a bracketing triplet
 * of abscissas \a ax, \a bx, \a cx (such that \a bx is between \a ax
 * and \a cx, and \a f(bx) is less than both \a f(ax) and \a f(cx)),
 * this routine performs a golden section search for the minimum error,
 * isolating it to a fractional precision of about \a tol. The abscissa
 * of the minimum is returned as \a xmin, and the minimum function value
 * is returned as \em golden, the returned function value.
 *
 * [Adapted from Numerical Recipes in C.]
 *
 * @param line A pointer to a \ref gnn_line buffer.
 * @param s    The index of the first pattern in the minibatch to be evaluated.
 * @param n    The size of the minibatch to be evaluated.
 * @param ax   Initial point of the bracketing triplet.
 * @param bx   Middle point of the bracketing triplet.
 * @param cx   Last point of the bracketing triplet.
 * @param xmin The location of the minimum.
 * @param tol  Fractional precision tolerance.
 * @return Returns the minimum function value.
 */
double
gnn_line_search_golden (gnn_line *line, size_t s, size_t n,
                        double ax, double bx, double cx,
                        double *xmin, double tol)
{
    double f1;
    double f2;
    double x0;
    double x1;
    double x2;
    double x3;

    x0 = ax;
    
    /* At any given time we will keep track */
    /* of four points, x0, x1, x2, x3.      */
    x3 = cx;

    /* Make x0 to x1 the smaller segment,     */
    /* and fill in the new point to be tried. */
    if (fabs (cx - bx) > fabs (bx - ax))
    {
        x1 = bx;
        x2 = bx + GNN_GOLDEN_C * (cx - bx);
    }
    else
    {
        x2 = bx;
        x1 = bx - GNN_GOLDEN_C * (bx - ax);
    }

    /* The initial function evaluations. Note that we never    */
    /* need to evaluate the function at the original endpoints. */
    gnn_line_pats (line, x1, gnnLineE, s, n);
    f1 = GNN_LINE_E (line);
    gnn_line_pats (line, x2, gnnLineE, s, n);
    f2 = GNN_LINE_E (line);

    while (fabs (x3 - x0) > tol * (fabs (x1) + fabs (x2)))
    {
        if (f2 < f1)
        {   /* One possible outcome, its housekeeping, */
            /* and a new function evaluation.          */
            GNN_SHIFT3 (x0, x1, x2, GNN_GOLDEN_R * x1 + GNN_GOLDEN_C * x3);
            f1 = f2;
            gnn_line_pats (line, x2, gnnLineE, s, n);
            f2 = GNN_LINE_E (line);
        }
        else
        {
            /* The other outcome, and its new function evaluation. */
            GNN_SHIFT3 (x3, x2, x1, GNN_GOLDEN_R * x2 + GNN_GOLDEN_C * x0);
            f2 = f1;
            gnn_line_pats (line, x1, gnnLineE, s, n);
            f1 = GNN_LINE_E (line);
        }
        
        /* Back to see if we are done. */
    }

    /* We are done. Output the best of the two current values. */
    if (f1 < f2)
    {
        *xmin = x1;
        return f1;
    }
    else
    {
        *xmin = x2;
        return f2;
    }
    
    return 0;
}

/**
 * @brief Bracketing procedure.
 * @ingroup gnn_line_search
 *
 * Given a \ref gnn_line buffer \a line, and given a bracketing triplet
 * of abscissas \a ax, \a bx, \a cx (such that \a bx is between \a ax and
 * \a cx, and \c f(bx) is less than both \c f(ax) and \c f(cx)), this
 * routine isolates the minimum to a fractional precision of about \a tol
 * using Brent’s method. The abscissa of the minimum is returned as \a xmin,
 * and the minimum function value is returned as \c brent, the
 * returned function value.
 *
 * @param line A pointer to a \ref gnn_line buffer.
 * @param s    The index of the first pattern in the minibatch to be evaluated.
 * @param n    The size of the minibatch to be evaluated.
 * @param ax   Initial point of the bracketing triplet.
 * @param bx   Middle point of the bracketing triplet.
 * @param cx   Last point of the bracketing triplet.
 * @param xmin The location of the minimum.
 * @param tol  Fractional precision tolerance.
 * @return Returns the minimum function value.
 */
double
gnn_line_search_brent (gnn_line *line, size_t s, size_t n,
                       double ax, double bx, double cx,
                       double *xmin, double tol)
{
    size_t iter;
    double a, b, d;
    double e, etemp;
    
    double u,  v,  w,  x,  xm;
    double fu, fv, fw, fx;
    
    double p, q, r;
    double tol1, tol2;
    
    /* This will be the distance moved on the step before last.*/
    e = 0.0;
    
    /* a and b must be in ascending order, */
    /* but input abscissas need not be.    */
    a = (ax < cx ? ax : cx);
    b = (ax > cx ? ax : cx);
    x = w = v = bx;
    
    /* Initializations... */
    gnn_line_pats (line, x, gnnLineE, s, n);
    fw = fv = fx = GNN_LINE_E (line);
    
    /* Main program loop. */
    for (iter=1; iter<=GNN_BRENT_ITMAX; iter++)
    {
        xm   = 0.5 * (a + b);
        tol2 = 2.0 * (tol1 = tol * fabs (x) + GNN_BRENT_ZEPS);
        
        /*  Test for done here. */
        if (fabs (x - xm) <= (tol2 - 0.5 * (b - a)))
        {
            *xmin = x;
            return fx;
        }
        
        /* Construct a trial parabolic fit. */
        if (fabs(e) > tol1)
        {
            r = (x - w) * (fx - fv);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w) * r;
            q = 2.0 * (q - r);

            if (q > 0.0)
                p = -p;
                
            q     = fabs (q);
            etemp = e;
            e     = d;

            if (    fabs (p) >= fabs (0.5 * q * etemp)
                 || p <= q * (a - x)
                 || p >= q * (b - x) )
            {
                d = GNN_BRENT_CGOLD * (e = (x >= xm ? a - x : b - x));
            }
            /* The above conditions determine the acceptability  */
            /* of the parabolic fit. Here we take the golden     */
            /* section step into the larger of the two segments. */
            else
            {
                /* Take the parabolic step. */
                d = p / q;
                u = x + d;
                if (u - a < tol2 || b - u < tol2)
                    d = GNN_USE_SIGN(tol1, xm-x);
            }
        }
        else
        {
            d = GNN_BRENT_CGOLD * (e = (x >= xm ? a - x : b - x));
        }
        u = (fabs (d) >= tol1 ? x + d : x + GNN_USE_SIGN(tol1, d));

        /* This is the one function evaluation per iteration. */
        gnn_line_pats (line, u, gnnLineE, s, n);
        fu = GNN_LINE_E (line);

        /*  Now decide what to do with our function evaluation. */
        if (fu <= fx)
        {
            if (u >= x)
                a = x;
            else
                b = x;

            /*  Housekeeping follows: */
            GNN_SHIFT3 (v, w, x, u);
            GNN_SHIFT3 (fv, fw, fx, fu);
        }
        else
        {
            if (u < x)
                a = u;
            else
                b = u;

            if (fu <= fw || w == x)
            {
                v = w;
                w = u;
                fv = fw;
                fw = fu;
            }
            else if (fu <= fv || v == x || v == w)
            {
                v  = u;
                fv = fu;
            }
        }
        /* Done with housekeeping. Back for another iteration. */
    }
    GSL_ERROR ("too many iterations in brent", GSL_EFAILED);
    
    /* Never get here. */
    *xmin = x;
    return fx;
}

/**
 * @brief Simple line search procedure.
 * @ingroup gnn_line_search_doc
 * @todo Not implemented yet.
 *
 * @param line A pointer to a \ref gnn_line buffer.
 * @param s    The index of the first pattern in the minibatch to be evaluated.
 * @param n    The size of the minibatch to be evaluated.
 * @param ax   Initial point of the bracketing triplet.
 * @param bx   Middle point of the bracketing triplet.
 * @param cx   Last point of the bracketing triplet.
 * @param xmin The location of the minimum.
 * @param tol  Fractional precision tolerance.
 * @return Returns the minimum function value.
 */
double
gnn_line_search_charalambous (gnn_line *line, size_t s, size_t n,
                              double ax, double bx, double cx,
                              double *xmin, double tol)
{
    return 0.0;
}


/**
 * @brief Simple line search procedure.
 * @ingroup gnn_line_search_doc
 *
 * This function implements a simple but robust line search routine. It returns
 * the \f$\alpha\f$, where
 *
 * \f[ F(x,w +\alpha d) \f]
 *
 * is minimized.
 *
 * The strategy is the following:
 * - Start with a fixed \f$\delta\f$
 * - Compute \f$F(x, w)\f$, \f$F(x, w+\delta d)\f$, \f$F(x, w+2\delta d)\f$
 * - If \f$F(x, w) \geq F(x, w+\delta d) \leq F(x, w+2\delta d)\f$, then
 *   proceed to interpolate.
 * - If \f$F(x, w) \geq F(x, w+\delta) > F(x, w+2\delta)\f$, then
 *   \f$d \leftarrow 2d\f$ and retry.
 * - If \f$F(x, w) < F(x, w+\delta)\f$, then \f$d \leftarrow d/2\f$ and retry.
 *
 * The interpolation is performed by a quadratic polynom:
 * \f[ P(\alpha) = a_0 + a_1 \alpha + a_2 \alpha^2 \f]
 * where the coefficients are
 *
 * \f[ \left [ \begin{array}{c}
 *               a_0 \\ a_1 \\ a_2 \\
 *             \end{array}
 *     \right ] =
 *     \left [ \begin{array}{ccc}
 *               1 & 0 & 0 \\
 *               -3/(2\delta)  & 2/\delta    & - 1/(2\delta) \\
 *               1/(2\delta^2) & -1/\delta^2 & 1/(2\delta^2) \\
 *             \end{array}
 *     \right ]
 *     \left [ \begin{array}{c}
 *                F(x, w) \\ F(x, w+\delta d) \\ F(x,w+2\delta d) \\
 *             \end{array}
 *     \right ]
 * \f]
 * which gives the solution for \f$\alpha\f$:
 * \f[ \alpha = -\frac{a_1}{2 a_2} \f]
 *
 * @param line A pointer to a \ref gnn_line buffer.
 * @param s    The index of the first pattern in the minibatch to be evaluated.
 * @param n    The size of the minibatch to be evaluated.
 * @param ax   Initial point of the bracketing triplet.
 * @param bx   Middle point of the bracketing triplet.
 * @param cx   Last point of the bracketing triplet.
 * @param xmin The location of the minimum.
 * @param tol  Fractional precision tolerance.
 * @return Returns the minimum function value.
 */
double
gnn_line_search_simple (gnn_line *line, size_t s, size_t n,
                        double ax, double bx, double cx,
                        double *xmin, double tol)
{
    return 0.0;
}

---