---

gnn_bfgs.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_bfgs.c
 *  @brief gnn_bfgs 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_bfgs_doc gnn_bfgs : Broyden-Fletcher-Goldfarb-Shanno Algorithm.
 * @ingroup gnn_trainer_doc
 *
 * This trainer implements the Broyden-Fletcher-Goldfarb-Shanno (BFGS)
 * optimization algorithm. It is a so-called \em Quasi-Newton algorithm.
 * The parameters are updated using the formula:
 *
 * \f[ w^{(\tau + 1)} = w^{(\tau)} + \alpha^{\tau} G^{(\tau)} g^{(\tau)} \f]
 *
 * where \f$G^{(\tau)}\f$ is an approximation to the inverse Hessian matrix,
 * which is built step by step using the formula:
 *
 * \f[ G^{ (\tau+1) } =
 *         G^{ (\tau + 1) }
 *         - \frac{ qq^T }{ v^T G^{(\tau)} v }
 *         + \frac{ pp^T }{ p^Tv }
 *         + ( v^T G^{(\tau)} v ) uu^T
 * \f]
 *
 * where the following vectors have been defined:
 *
 * \f[ p = w^{(\tau + 1)} - w^{(\tau)}                  \f]
 * \f[ q = G^{(\tau)} v                                 \f]
 * \f[ v = g^{(\tau + 1)} - g^{(\tau)}                  \f]
 * \f[ u = \frac{p}{p^T v} - \frac{q}{v^T G^{(\tau)} v} \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.
 *
 * Although theoretically well founded and very robust, the BFGS algorithm
 * is very expensive, both in time and space. The approximation of the
 * Hessian matrix is kept in memory, which is of size \f$O(l^2)\f$ where
 * \f$l\f$ are the number of parameters, leading to a prohibitive amount
 * of storage requirement for greater Gradient Machines.
 */



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

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



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

static int
gnn_bfgs_reset (gnn_trainer *trainer);

static int
gnn_bfgs_iteration (gnn_bfgs *bf);

static int
gnn_bfgs_train (gnn_trainer *trainer);

static void
gnn_bfgs_destroy (gnn_trainer *trainer);




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

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

    assert (trainer != NULL);

    bf = (gnn_bfgs *) 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->u);
    gsl_vector_set_zero (bf->v);
    gsl_vector_set_zero (bf->p);
    gsl_vector_set_zero (bf->q);

    /* initialize Hessian approximation */
    gsl_matrix_set_identity (bf->G);

    return 0;
}

/**
 * @brief Computes the step-wise approximation of the Hessian matrix.
 * @ingroup gnn_bfgs_doc
 *
 * @param  trainer A pointer to a \ref gnn_bfgs.
 * @return Returns 0 if suceeded.
 */
static int
gnn_bfgs_iteration (gnn_bfgs *bf)
{
    double pTv;
    double vTGv;
    double IpTv;
    double IvTGv;

    /* 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 q vector */
    gsl_blas_dgemv (CblasNoTrans, 1.0, bf->G, bf->v, 0.0, bf->q);

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

    /* compute pTv */
    gsl_blas_ddot (bf->p, bf->v, &pTv);
    
    /* compute inverses */
    IpTv = 1.0 / GNN_MAX (pTv, GNN_TINY);
    IvTGv = 1.0 / GNN_MAX (vTGv, GNN_TINY);

    /* compute u vector */
    gsl_vector_memcpy (bf->u, bf->p);
    gsl_vector_scale (bf->u, 1.0 / pTv);
    gsl_blas_daxpy (- 1.0 / vTGv, bf->q, bf->u);

    /* compute Hessian approximation */
    gsl_blas_dgemm (CblasNoTrans, CblasTrans, -IvTGv, bf->Q, bf->Q, 1, bf->G);
    gsl_blas_dgemm (CblasNoTrans, CblasTrans,   IpTv, bf->P, bf->P, 1, bf->G);
    gsl_blas_dgemm (CblasNoTrans, CblasTrans,   vTGv, bf->U, bf->U, 1, bf->G);

    return 0;
}

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

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

    size_t s, n;
    gnn_bfgs *bf;

    /* get view */
    bf = (gnn_bfgs *) 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);

    /* compute new search direction and perform one */
    /* step of the Hessian matrix approximation     */
    if (bf->iteration % bf->restart == 0)
    {
        gsl_matrix_set_identity (bf->G);
    }
    else
    {
        gnn_bfgs_iteration (bf);
    }

    /* prepare for next iteration */
    gsl_vector_memcpy (bf->wold, bf->wnew);
    gsl_vector_memcpy (bf->gold, bf->gnew);

    /* build new search direction */
    gsl_blas_dgemv (CblasNoTrans, 1.0, bf->G, bf->gnew, 0.0, 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_bfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_bfgs.
 */
static void
gnn_bfgs_destroy (gnn_trainer *trainer)
{
    gnn_bfgs *bf;

    assert (trainer != NULL);

    bf = (gnn_bfgs *) trainer;

    if (bf->G != NULL)
        gsl_matrix_free (bf->G);
    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->u != NULL)
        gsl_vector_free (bf->u);
    if (bf->v != NULL)
        gsl_vector_free (bf->v);
    if (bf->p != NULL)
        gsl_vector_free (bf->p);
    if (bf->q != NULL)
        gsl_vector_free (bf->q);
    if (bf->line != NULL)
        gnn_line_destroy (bf->line);

    return;
}



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

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

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

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

    /* initialize */
    status = gnn_trainer_init (trainer,
                               "gnn_bfgs",
                               node,
                               crit,
                               data,
                               gnn_bfgs_reset,
                               gnn_bfgs_train,
                               gnn_bfgs_destroy);
    if (status)
    {
        GSL_ERROR_VAL ("couldn't initialize gnn_bfgs",
                       GSL_EFAILED, NULL);
    }

    /* set fields */
    bf->step      = GNN_BFGS_STEP;
    bf->tol       = GNN_BFGS_TOL;
    bf->iteration = 0;
    bf->restart   = GNN_BFGS_RESTART;

    bf->alpha = GNN_BFGS_ALPHA;

    /* allocate memory for all needed buffers */
    l = gnn_node_param_get_size (node);
    bf->G    = gsl_matrix_alloc (l, l);
    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->u    = gsl_vector_alloc (l);
    bf->v    = gsl_vector_alloc (l);
    bf->p    = gsl_vector_alloc (l);
    bf->q    = gsl_vector_alloc (l);
    bf->line = gnn_line_new (trainer->grad, NULL);

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

    /* build needed matrix views of the vectors */
    bf->Uview = gsl_matrix_view_vector (bf->u, l, 1);
    bf->Vview = gsl_matrix_view_vector (bf->v, l, 1);
    bf->Pview = gsl_matrix_view_vector (bf->p, l, 1);
    bf->Qview = gsl_matrix_view_vector (bf->q, l, 1);
    bf->U     = &(bf->Uview.matrix);
    bf->V     = &(bf->Vview.matrix);
    bf->P     = &(bf->Pview.matrix);
    bf->Q     = &(bf->Qview.matrix);

    return trainer;
}



/**
 * @brief Sets the precision tolerance for the line search procedure.
 * @ingroup gnn_bfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_bfgs.
 * @param tol A stricly positive real value.
 * @return Returns 0 if succeeded.
 */
int
gnn_bfgs_set_tol (gnn_trainer *trainer, double tol)
{
    gnn_bfgs *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_bfgs *) trainer;
    bf->tol = tol;

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_bfgs *) trainer;

    return bf->tol;
}

/**
 * @brief Sets the initial step for the interval bracketing procedure.
 * @ingroup gnn_bfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_bfgs.
 * @param step A stricly positive real value.
 * @return Returns 0 if succeeded.
 */
int
gnn_bfgs_set_step (gnn_trainer *trainer, double step)
{
    gnn_bfgs *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_bfgs *) trainer;
    bf->step = step;

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_bfgs *) trainer;

    return bf->step;
}

/**
 * @brief Sets the number of iterations before restarting.
 * @ingroup gnn_bfgs_doc
 *
 * @param trainer A pointer to a \ref gnn_bfgs.
 * @param restart A stricly positive integer.
 * @return Returns 0 if succeeded.
 */
int
gnn_bfgs_set_restart (gnn_trainer *trainer, size_t restart)
{
    gnn_bfgs *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_bfgs *) trainer;
    bf->restart = restart;

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_bfgs *) trainer;

    return bf->restart;
}

/**
 * @brief Sets the line search procedure.
 * @ingroup gnn_bfgs_doc
 *
 * This function sets a new line search procedure used by the BFGS trainer.
 *
 * \code
 * gnn_trainer *trainer;
 * trainer = gnn_bfgs_new (node, crit, data);
 *
 * // use the Golden-Section line search
 * gnn_bfgs_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_bfgs.
 * @param lsearch A pointer to a line search procedure.
 * @return Returns 0 if succeeded.
 */
int
gnn_bfgs_set_line_search (gnn_trainer *trainer, gnn_line_search_type lsearch)
{
    gnn_bfgs *bf;

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

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

    return 0;
}

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

    assert (trainer != NULL);

    bf = (gnn_bfgs *) trainer;

    return bf->alpha;
}

---