---

gnn_conjugate_gradient.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_conjugate_gradient.c
 *  @brief Conjugate Gradient Trainer Implementation.
 *
 *  @date   : 07-09-03 12:08, 13-09-03 18:41
 *  @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.
 */



/**
 * @defgroup gnn_conjugate_gradient_doc gnn_conjugate_gradient : Conjugate Gradient Descent Algorithm.
 * @ingroup gnn_trainer_doc
 *
 * This trainer implements the Conjugate Gradient minimization procedure.
 * The conjugate gradients method takes one step after another in special
 * conjugate directions, minimizing the error function along the current
 * line.
 *
 * The conjugate directions are built using a sequential procedure, given
 * by the formula:
 *
 * \f[ d_{new} = g_{new} - \beta_{new} d_{old} \f]
 *
 * where \f$d_{new}\f$ is the new search direction, which is conjugate to
 * the old \f$d_{old}\f$, \f$g_{new}\f$ and \f$g_{old}\f$ are the new and
 * the previous computed gradients with respect to the function's parameters.
 * The \f$\beta\f$ is a coefficient which can be computed in different ways.
 * There are three forms:
 * - \ref gnn_conjugate_gradient_hestenes_stiefel : The original form,
 *   discovered by Hestenes and Stiefel.
 * - \ref gnn_conjugate_gradient_fletcher_reeves :
 *   The Fletcher-Reeves form.
 * - \ref gnn_conjugate_gradient_polak_ribiere :
 *   The Polak-Ribière form.
 *
 * Informally, two directions are said to be conjugate, if they are "ortogonal"
 * in a sense considering the second-order Taylor approximation of the error
 * surface.
 *
 * After choosing the direction, the algorithm seeks the minimum along this
 * line. This is performed, in practice, using a line-search procedure.
 * You can install any of the available line-search procedures for the
 * training algoritm. Please refer to (\ref gnn_line_search) for further
 * details.
 */



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

#include <math.h>
#include "gnn_conjugate_gradient.h"



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

static int
gnn_conjugate_gradient_reset (gnn_trainer *trainer);

static int
gnn_conjugate_gradient_train (gnn_trainer *trainer);

static void
gnn_conjugate_gradient_destroy (gnn_trainer *trainer);




/******************************************/
/* Static Implementation                  */
/******************************************/

/**
 * @brief The trainer's "reset" implementation.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param  trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns 0 if suceeded.
 */
static int
gnn_conjugate_gradient_reset (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cgtrainer;

    assert (trainer != NULL);

    cgtrainer = (gnn_conjugate_gradient *) trainer;

    /* reset iteration counter */
    cgtrainer->iteration = 0;

    /* clear optimization information */
    gsl_vector_set_zero (cgtrainer->buf);
    gsl_vector_set_zero (cgtrainer->gnew);
    gsl_vector_set_zero (cgtrainer->gold);

    return 0;
}

/**
 * @brief The trainer's "train" implementation.
 * @ingroup gnn_conjugate_gradient_descent_doc
 *
 * @param  trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns 0 if succeeded.
 */
static int
gnn_conjugate_gradient_train (gnn_trainer *trainer)
{
    double alpha;
    double beta;
    double ax, bx, cx;
    double fa, fb, fc;
    size_t s, n;
    gnn_conjugate_gradient *cg;

    /* get view */
    cg = (gnn_conjugate_gradient *) trainer;

    /* process minibatch */
    gnn_trainer_batch_process (trainer);
    
    /* copy gradient */
    gsl_vector_memcpy (cg->gnew, gnn_trainer_batch_get_dw (trainer));

    /* check if the algorithm should be restarted */
    if (cg->iteration % cg->restart == 0)
    {
        /* initialize direction */
        gsl_vector_memcpy (cg->line->d, cg->gnew);
    }
    else
    {
        /* compute beta */
        beta = cg->beta (trainer);
        if (beta < 0.0)
        {
            /* beta is negative, restart direction. */
            gsl_vector_memcpy (cg->line->d, cg->gnew);
        }
        else
        {
            /* beta is positive, compute new direction */
            gsl_vector_scale (cg->line->d, beta);
            gsl_vector_sub (cg->line->d, cg->gnew);
        }
    }

    /* perform line search over the current minibatch's patterns */
    ax = 0.0;
    bx = cg->step;
    s = gnn_trainer_get_pattern_index (trainer);
    n = gnn_trainer_batch_get_size (trainer);
    
    gnn_line_search_bracket (cg->line, s, n, &ax, &bx, &cx, &fa, &fb, &fc);
    cg->alpha (cg->line, s, n, ax, bx, cx, &alpha, cg->tol);

    /* build new point */
    gsl_vector_memcpy (cg->buf, cg->line->d);
    gsl_vector_scale (cg->buf, alpha);

    /* sum to parameters */
    gsl_vector_add (cg->line->w, cg->buf);

    /* update parameters */
    gnn_node_param_set (trainer->node, cg->line->w);

    /* store old gradient */
    gsl_vector_memcpy (cg->gold, cg->gnew);

    /* move to next minibatch */
    gnn_trainer_batch_next (trainer);

    /* update iteration counter */
    cg->iteration++;

    return 0;
}

/**
 * @brief The trainers "destroy" implementation.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 */
static void
gnn_conjugate_gradient_destroy (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cgtrainer;

    assert (trainer != NULL);

    cgtrainer = (gnn_conjugate_gradient *) trainer;

    if (cgtrainer->buf != NULL)
        gsl_vector_free (cgtrainer->buf);

    if (cgtrainer->gnew != NULL)
        gsl_vector_free (cgtrainer->gnew);

    if (cgtrainer->gold != NULL)
        gsl_vector_free (cgtrainer->gold);

    if (cgtrainer->line != NULL)
        gnn_line_destroy (cgtrainer->line);
        
    return;
}



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

/**
 * @brief Creates a new conjugate gradient descent trainer.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function creates a new conjugate gradients trainer
 * (\ref gnn_conjugate_gradient).
 *
 * @param  node A pointer to a \ref gnn_node.
 * @param  crit A pointer to a \ref gnn_criterion.
 * @param  data A pointer to a \ref gnn_dataset.
 * @return Returns a pointer to a new \ref gnn_conjugate_gradient trainer.
 */
gnn_trainer *
gnn_conjugate_gradient_new (gnn_node *node,
                            gnn_criterion *crit,
                            gnn_dataset *data)
{
    int status;
    size_t l;
    gnn_trainer *trainer;
    gnn_conjugate_gradient *cgtrainer;

    /* allocate memory for the trainer */
    cgtrainer =
            (gnn_conjugate_gradient *) malloc (sizeof (gnn_conjugate_gradient));
    if (cgtrainer == NULL)
    {
        GSL_ERROR_VAL ("couldn't allocate memory for gnn_conjugate_gradient",
                       GSL_ENOMEM, NULL);
    }

    /* get view as gnn_trainer */
    trainer = (gnn_trainer *) cgtrainer;

    /* initialize */
    status = gnn_trainer_init (trainer,
                               "gnn_conjugate_gradient",
                               node,
                               crit,
                               data,
                               gnn_conjugate_gradient_reset,
                               gnn_conjugate_gradient_train,
                               gnn_conjugate_gradient_destroy);
    if (status)
    {
        GSL_ERROR_VAL ("couldn't initialize gnn_conjugate_gradient",
                       GSL_EFAILED, NULL);
    }

    /* set fields */
    cgtrainer->step      = GNN_CONJUGATE_GRADIENT_STEP;
    cgtrainer->tol       = GNN_CONJUGATE_GRADIENT_TOL;
    cgtrainer->iteration = 0;
    cgtrainer->restart   = GNN_CONJUGATE_GRADIENT_RESTART;

    cgtrainer->alpha = GNN_CONJUGATE_GRADIENT_ALPHA;
    cgtrainer->beta  = GNN_CONJUGATE_GRADIENT_BETA;

    l = gnn_node_param_get_size (node);
    cgtrainer->gnew = gsl_vector_alloc (l);
    cgtrainer->gold = gsl_vector_alloc (l);
    cgtrainer->buf  = gsl_vector_alloc (l);
    cgtrainer->line = gnn_line_new (trainer->grad, NULL);

    if (   cgtrainer->buf  == NULL
        || cgtrainer->gnew == NULL
        || cgtrainer->gold == NULL )
    {
        gnn_trainer_destroy (trainer);
        GSL_ERROR_VAL ("couldn't allocate memory for gnn_conjugate_gradient",
                       GSL_ENOMEM, NULL);
    }

    return trainer;
}



/**
 * @brief The Polak-Ribière form for the \f$\beta\f$ coefficient.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function returns the Polak-Ribière form of the \f$\beta\f$ coefficient
 * for the evaluation of the new conjugate direction:
 *
 * \f[ d_{new} = g_{new} - \beta_{new} d_{old} \f]
 *
 * where
 *
 * \f[ \beta_{new} = \frac{g_{new}^T (g_{new} - g_{old})}
 *                        {     g_{old}^T - g_{old}     }
 * \f]
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns the \f$beta\f$ coefficient.
 */
double
gnn_conjugate_gradient_polak_ribiere (gnn_trainer *trainer)
{
    double beta_u;
    double beta_l;
    gnn_conjugate_gradient *cgtrainer;

    /* get view */
    cgtrainer = (gnn_conjugate_gradient *) trainer;

    /* compute beta */
    gsl_vector_memcpy (cgtrainer->buf, cgtrainer->gnew);
    gsl_vector_sub (cgtrainer->buf, cgtrainer->gold);
    gsl_blas_ddot (cgtrainer->gnew, cgtrainer->buf, &beta_u);
    gsl_blas_ddot (cgtrainer->gold, cgtrainer->gold, &beta_l);

    return ( beta_u / beta_l );
}

/**
 * @brief The Hestenes-Stiefel form for the \f$\beta\f$ coefficient.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function returns the Hestenes-Stiefel form of the \f$\beta\f$
 * coefficient for the evaluation of the new conjugate direction:
 *
 * \f[ d_{new} = g_{new} - \beta_{new} d_{old} \f]
 *
 * where
 *
 * \f[ \beta_{new} = \frac{ g_{new}^T (g_{new}   - g_{old}) }
 *                        { d_{old}^T (g_{new}^T - g_{old}) }
 * \f]
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns the \f$beta\f$ coefficient.
 */
double
gnn_conjugate_gradient_hestenes_stiefel (gnn_trainer *trainer)
{
    double beta_u;
    double beta_l;
    gnn_conjugate_gradient *cgtrainer;

    /* get view */
    cgtrainer = (gnn_conjugate_gradient *) trainer;

    /* compute beta */
    gsl_vector_memcpy (cgtrainer->buf, cgtrainer->gnew);
    gsl_vector_sub (cgtrainer->buf, cgtrainer->gold);
    
    gsl_blas_ddot (cgtrainer->gnew, cgtrainer->buf, &beta_u);
    gsl_blas_ddot (cgtrainer->line->d, cgtrainer->buf, &beta_l);

    return ( beta_u / beta_l );
}

/**
 * @brief The Fletcher-Reeves form for the \f$\beta\f$ coefficient.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function returns the Fletcher-Reeves form of the \f$\beta\f$
 * coefficient for the evaluation of the new conjugate direction:
 *
 * \f[ d_{new} = g_{new} - \beta_{new} d_{old} \f]
 *
 * where
 *
 * \f[ \beta_{new} = \frac{ g_{new}^T g_{new} }
 *                        { d_{old}^T g_{old} }
 * \f]
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns the \f$beta\f$ coefficient.
 */
double
gnn_conjugate_gradient_fletcher_reeves (gnn_trainer *trainer)
{
    double beta_u;
    double beta_l;
    gnn_conjugate_gradient *cgtrainer;

    /* get view */
    cgtrainer = (gnn_conjugate_gradient *) trainer;

    /* compute beta */
    gsl_blas_ddot (cgtrainer->gnew, cgtrainer->gnew, &beta_u);
    gsl_blas_ddot (cgtrainer->gold, cgtrainer->gold, &beta_l);

    return ( beta_u / beta_l );
}

/**
 * @brief Sets the precision tolerance for the line search procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @param tol A stricly positive real value.
 * @return Returns 0 if succeeded.
 */
int
gnn_conjugate_gradient_set_tol (gnn_trainer *trainer, double tol)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    /* check value */
    if (tol <= 0.0)
    {
        GSL_ERROR ("tolerance should be stricly greater than zero",
                   GSL_EINVAL);
    }

    /* set value */
    cg = (gnn_conjugate_gradient *) trainer;
    cg->tol = tol;

    return 0;
}

/**
 * @brief Gets the tolerance for the line search procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns the tolerance's value.
 */
double
gnn_conjugate_gradient_get_tol (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    cg = (gnn_conjugate_gradient *) trainer;

    return cg->tol;
}

/**
 * @brief Sets the initial step for the interval bracketing procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @param step A stricly positive real value.
 * @return Returns 0 if succeeded.
 */
int
gnn_conjugate_gradient_set_step (gnn_trainer *trainer, double step)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    /* check value */
    if (step <= 0.0)
    {
        GSL_ERROR ("step should be stricly greater than zero",
                   GSL_EINVAL);
    }

    /* set value */
    cg = (gnn_conjugate_gradient *) trainer;
    cg->step = step;

    return 0;
}

/**
 * @brief Gets the initial step for the interval bracketing procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return The trainer's internal step.
 */
double
gnn_conjugate_gradient_get_step (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    cg = (gnn_conjugate_gradient *) trainer;

    return cg->step;
}

/**
 * @brief Sets the number of iterations before restarting.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @param restart A stricly positive integer.
 * @return Returns 0 if succeeded.
 */
int
gnn_conjugate_gradient_set_restart (gnn_trainer *trainer, size_t restart)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    /* check value */
    if (restart <= 0.0)
    {
        GSL_ERROR ("restart iteration should be stricly greater than zero",
                   GSL_EINVAL);
    }

    /* set value */
    cg = (gnn_conjugate_gradient *) trainer;
    cg->restart = restart;

    return 0;
}

/**
 * @brief Gets the number of iterations before reinitializing the direction.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function returns the number of iterations executed by the conjugate
 * gradients trainer before reinitializing the search direction.
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns the number of iterations.
 */
size_t
gnn_conjugate_gradient_get_restart (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    cg = (gnn_conjugate_gradient *) trainer;

    return cg->restart;
}

/**
 * @brief Sets the line search procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function sets a new line search procedure used by the conjugate
 * gradients trainer.
 *
 * \code
 * gnn_trainer *trainer;
 * trainer = gnn_conjugate_gradient_new (node, crit, data);
 *
 * // use the Golden-Section line search
 * gnn_conjugate_gradient_set_line_search (trainer, gnn_line_search_golden);
 * \endcode
 *
 * Please refer to (\ref gnn_line_search_doc) for the available line search
 * procedures.
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @param lsearch A pointer to a line search procedure.
 * @return Returns 0 if succeeded.
 */
int
gnn_conjugate_gradient_set_line_search (gnn_trainer *trainer,
                                        gnn_line_search_type lsearch)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);
    assert (lsearch != NULL);

    /* set value */
    cg = (gnn_conjugate_gradient *) trainer;
    cg->alpha = lsearch;

    return 0;
}

/**
 * @brief Gets the installed line search procedure.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return Returns a pointer to the installed line-search procedure.
 */
gnn_line_search_type
gnn_conjugate_gradient_get_alpha (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    cg = (gnn_conjugate_gradient *) trainer;

    return cg->alpha;
}

/**
 * @brief Sets the form of the \f$beta\f$ coefficient.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function sets a new form for evaluating the \f$beta\f$ coefficient
 * used by the conjugate gradients method to build the new search direction.
 *
 * \code
 * gnn_trainer *trainer;
 * trainer = gnn_conjugate_gradient_new (node, crit, data);
 *
 * // use the Fletcher-Reeves form
 * gnn_conjugate_gradient_set_beta (trainer,
 *                                  gnn_conjugate_gradient_fletcher_reeves);
 * \endcode
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @param beta A pointer to the a \f$\beta\f$ evaluation procedure.
 * @return Returns 0 if succeeded.
 */
int
gnn_conjugate_gradient_set_beta (gnn_trainer *trainer,
                                 gnn_conjugate_gradient_beta beta)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);
    assert (beta != NULL);

    /* set value */
    cg = (gnn_conjugate_gradient *) trainer;
    cg->beta = beta;

    return 0;
}

/**
 * @brief Gets the \f$\beta\f$ evaluation function.
 * @ingroup gnn_conjugate_gradient_doc
 *
 * This function returns a pointer to the installed \f$\beta\f$ evaluation
 * procedure.
 *
 * @param trainer A pointer to a \ref gnn_conjugate_gradient.
 * @return A pointer to a \f$\beta\f$ function.
 */
gnn_conjugate_gradient_beta
gnn_conjugate_gradient_get_beta (gnn_trainer *trainer)
{
    gnn_conjugate_gradient *cg;

    assert (trainer != NULL);

    cg = (gnn_conjugate_gradient *) trainer;

    return cg->beta;
}

---