---

gnn_lmbfgs.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_lmbfgs.c
 *  @brief gnn_lmbfgs Implementation.
 *
 *  @date   : 05-10-03 01:05
 *  @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_lmbfgs_doc gnn_lmbfgs : Limited Memory BFGS Algorithm.
 * @ingroup gnn_trainer_doc
 *
 * This trainer implements the Limited Memory Broyden-Fletcher-Goldfarb-Shanno
 * (LMBFGS) optimization algorithm. It is a so-called \em Quasi-Newton
 * method. The parameters are updated using the formula:
 *
 * \f[ w^{(\tau + 1)} = w^{(\tau)} + \alpha^{\tau} d^{(\tau)} \f]
 *
 * where \f$d^{(\tau)}\f$ is a search direction given by:
 *
 * \f[ d^{(\tau+1)} = -g^{(\tau+1)} + Ap + Bp \f]
 *
 * and where the following vectors have been defined:
 *
 * \f[ p = w^{(\tau + 1)} - w^{(\tau)}                  \f]
 * \f[ v = g^{(\tau + 1)} - g^{(\tau)}                  \f]
 *
 * and \f$g^{(\tau)}\f$ is the \f$\tau\f$-th gradient
 * \f$\frac{\partial E}{\partial w}\f$. The step size \f$\alpha^{(\tau)}\f$
 * is found using a line-minimization, which can be one of those available
 * in \ref gnn_line_search. The coefficients \f$A\f$ and \f$B\f$ are defined
 * as
 *
 * \f[ A = -\left ( 1 + \frac{ v^T v }{ p^T v } \right )
 *             \frac{ p^T g^{(\tau + 1)} }{ p^T v }
 *         + \frac{ v^T g^{(\tau + 1)} }{ p^T v }
 * \f]
 * \f[ B = \frac{ p^T g^{(\tau + 1)} }{ p^T v } \f]
 *
 * The difference with the BFGS method lies in the fact that the expensive
 * approximation and storage of the inverse Hessian matrix has been replaced
 * by the unit matrix. The resulting algorithm takes only \f$O(l)\f$ space,
 * where \f$l\f$ is the number of parameters of the Gradient Machine. The
 * algorithm is comparable to the conjugate gradients algorithm
 * (refer to \ref gnn_conjugate_gradient).
 */



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

#include <math.h>
#include <gsl/gsl_blas.h>
#include "gnn_utilities.h"
#include "gnn_lmbfgs.h"



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

static int
gnn_lmbfgs_reset (gnn_trainer *trainer);

static int
gnn_lmbfgs_iteration (gnn_lmbfgs *bf, double *A, double *B);

static int
gnn_lmbfgs_train (gnn_trainer *trainer);

static void
gnn_lmbfgs_destroy (gnn_trainer *trainer);




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

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

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

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

    /* clear optimization information */
    gsl_vector_set_zero (bf->wnew);
    gsl_vector_set_zero (bf->wold);
    gsl_vector_set_zero (bf->gnew);
    gsl_vector_set_zero (bf->gold);
    gsl_vector_set_zero (bf->v);
    gsl_vector_set_zero (bf->p);

    return 0;
}

/**
 * @brief Computes the A and B coefficients.
 * @ingroup gnn_lmbfgs_doc
 *
 * @param  trainer A pointer to a \ref gnn_lmbfgs.
 * @param  A       A pointer where to store the \f$A\f$ coefficient.
 * @param  B       A pointer where to store the \f$B\f$ coefficient.
 * @return Returns 0 if suceeded.
 */
static int
gnn_lmbfgs_iteration (gnn_lmbfgs *bf, double *A, double *B)
{
    double pTv;
    double pTg;
    double vTg;
    double vTv;
    
    /* compute p vector */
    gsl_vector_memcpy (bf->p, bf->wnew);
    gsl_vector_sub (bf->p, bf->wold);

    /* compute v vector */
    gsl_vector_memcpy (bf->v, bf->gnew);
    gsl_vector_sub (bf->v, bf->gold);

    /* compute pTv and its inverse */
    gsl_blas_ddot (bf->p, bf->v, &pTv);
    pTv = GNN_MAX (pTv, GNN_TINY);

    /* compute vTv */
    gsl_blas_ddot (bf->v, bf->v, &vTv);

    /* compute pTg */
    gsl_blas_ddot (bf->p, bf->gnew, &vTv);

    /* compute vTg */
    gsl_blas_ddot (bf->v, bf->gnew, &vTv);

    /* compute coefficients */
    *A = -(1.0 + vTv / pTv) * pTg / pTv + vTg / pTv;
    *B =  pTg / pTv;

    return 0;
}

/**
 * @brief The trainer's "train" implementation.
 * @ingroup gnn_lmbfgs_descent_doc
 *
 * @param  trainer A pointer to a \ref gnn_lmbfgs.
 * @return Returns 0 if succeeded.
 */
static int
gnn_lmbfgs_train (gnn_trainer *trainer)
{
    double A;
    double B;
    double alpha;

    double ax, bx, cx;
    double fa, fb, fc;

    size_t s, n;
    gnn_lmbfgs *bf;

    /* get view */
    bf = (gnn_lmbfgs *) trainer;

    /* process minibatch */
    gnn_trainer_batch_process (trainer);

    /* copy gradient */
    gsl_vector_memcpy (bf->gnew, gnn_trainer_batch_get_dw (trainer));

    /* copy current parameter vector */
    gnn_node_param_get (trainer->node, bf->wnew);

    /* check if the algorithm should be restarted */
    if (bf->iteration % bf->restart == 0)
    {
        /* only take the negative gradient search direction at start */
        A = 0.0;
        B = 0.0;
    }
    else
    {
        /* compute new search direction */
        gnn_lmbfgs_iteration (bf, &A, &B);
    }
    
    /* prepare for next iteration */
    gsl_vector_memcpy (bf->wold, bf->wnew);
    gsl_vector_memcpy (bf->gold, bf->gnew);

    /* build new search direction */
    gsl_vector_memcpy (bf->line->d, bf->gnew);
    gsl_vector_scale (bf->line->d, -1.0);
    gsl_blas_daxpy (A, bf->p, bf->line->d);
    gsl_blas_daxpy (B, bf->v, bf->line->d);

    /* perform line search over the current minibatch's patterns */
    ax = 0.0;
    bx = bf->step;
    s = gnn_trainer_get_pattern_index (trainer);
    n = gnn_trainer_batch_get_size (trainer);

    gnn_line_search_bracket (bf->line, s, n, &ax, &bx, &cx, &fa, &fb, &fc);
    bf->alpha (bf->line, s, n, ax, bx, cx, &alpha, bf->tol);

    /* build new origin */
    gsl_vector_scale (bf->line->d, alpha);
    gsl_vector_add (bf->line->w, bf->line->d);

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

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

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

    return 0;
}

/**
 * @brief The trainers "destroy" implementation.
 * @ingroup gnn_lmbfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_lmbfgs.
 */
static void
gnn_lmbfgs_destroy (gnn_trainer *trainer)
{
    gnn_lmbfgs *bf;

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

    if (bf->wnew != NULL)
        gsl_vector_free (bf->wnew);
    if (bf->wold != NULL)
        gsl_vector_free (bf->wold);
    if (bf->gnew != NULL)
        gsl_vector_free (bf->gnew);
    if (bf->gold != NULL)
        gsl_vector_free (bf->gold);
    if (bf->v != NULL)
        gsl_vector_free (bf->v);
    if (bf->p != NULL)
        gsl_vector_free (bf->p);
    if (bf->line != NULL)
        gnn_line_destroy (bf->line);

    return;
}



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

/**
 * @brief Creates a new LMBFGS trainer.
 * @ingroup gnn_lmbfgs_doc
 *
 * This function creates a new LMBFGS trainer (\ref gnn_lmbfgs).
 *
 * @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_lmbfgs trainer.
 */
gnn_trainer *
gnn_lmbfgs_new (gnn_node *node, gnn_criterion *crit, gnn_dataset *data)
{
    int status;
    size_t l;
    gnn_trainer *trainer;
    gnn_lmbfgs *bf;

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

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

    /* initialize */
    status = gnn_trainer_init (trainer,
                               "gnn_lmbfgs",
                               node,
                               crit,
                               data,
                               gnn_lmbfgs_reset,
                               gnn_lmbfgs_train,
                               gnn_lmbfgs_destroy);
    if (status)
    {
        GSL_ERROR_VAL ("couldn't initialize gnn_lmbfgs",
                       GSL_EFAILED, NULL);
    }

    /* set fields */
    bf->step      = GNN_LMBFGS_STEP;
    bf->tol       = GNN_LMBFGS_TOL;
    bf->iteration = 0;
    bf->restart   = GNN_LMBFGS_RESTART;

    bf->alpha = GNN_LMBFGS_ALPHA;

    /* allocate memory for all needed buffers */
    l = gnn_node_param_get_size (node);
    bf->wnew = gsl_vector_alloc (l);
    bf->wold = gsl_vector_alloc (l);
    bf->gnew = gsl_vector_alloc (l);
    bf->gold = gsl_vector_alloc (l);
    bf->v    = gsl_vector_alloc (l);
    bf->p    = gsl_vector_alloc (l);
    bf->line = gnn_line_new (trainer->grad, NULL);

    if (   bf->wnew == NULL
        || bf->wold == NULL
        || bf->gnew == NULL
        || bf->gold == NULL
        || bf->v    == NULL
        || bf->p    == NULL
        || bf->line == NULL )
    {
        gnn_trainer_destroy (trainer);
        GSL_ERROR_VAL ("couldn't allocate memory for gnn_lmbfgs",
                       GSL_ENOMEM, NULL);
    }

    return trainer;
}



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

    assert (trainer != NULL);

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

    /* set value */
    bf = (gnn_lmbfgs *) trainer;
    bf->tol = tol;

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

    return bf->tol;
}

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

    assert (trainer != NULL);

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

    /* set value */
    bf = (gnn_lmbfgs *) trainer;
    bf->step = step;

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

    return bf->step;
}

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

    assert (trainer != NULL);

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

    /* set value */
    bf = (gnn_lmbfgs *) trainer;
    bf->restart = restart;

    return 0;
}

/**
 * @brief Gets the number of iterations before reinitializing the direction.
 * @ingroup gnn_lmbfgs_doc
 *
 * This function returns the number of iterations executed by the LMBFGS
 * trainer before reinitializing the search direction.
 *
 * @param trainer A pointer to a \ref gnn_lmbfgs.
 * @return Returns the number of iterations.
 */
size_t
gnn_lmbfgs_get_restart (gnn_trainer *trainer)
{
    gnn_lmbfgs *bf;

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

    return bf->restart;
}

/**
 * @brief Sets the line search procedure.
 * @ingroup gnn_lmbfgs_doc
 *
 * This function sets a new line search procedure used by the LMBFGS trainer.
 *
 * \code
 * gnn_trainer *trainer;
 * trainer = gnn_lmbfgs_new (node, crit, data);
 *
 * // use the Golden-Section line search
 * gnn_lmbfgs_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_lmbfgs.
 * @param lsearch A pointer to a line search procedure.
 * @return Returns 0 if succeeded.
 */
int
gnn_lmbfgs_set_line_search (gnn_trainer *trainer, gnn_line_search_type lsearch)
{
    gnn_lmbfgs *bf;

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

    /* set value */
    bf = (gnn_lmbfgs *) trainer;
    bf->alpha = lsearch;

    return 0;
}

/**
 * @brief Gets the installed line search procedure.
 * @ingroup gnn_lmbfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_lmbfgs.
 * @return Returns a pointer to the installed line-search procedure.
 */
gnn_line_search_type
gnn_lmbfgs_get_alpha (gnn_trainer *trainer)
{
    gnn_lmbfgs *bf;

    assert (trainer != NULL);

    bf = (gnn_lmbfgs *) trainer;

    return bf->alpha;
}

---