---

gnn_hessian.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_hessian.c
 *  @brief Hessian Matrix Evaluation Routines Implementation.
 *
 *  @date   : 19-09-03 15:37
 *  @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 Evaluation of Hessian Matrix.
 * @defgroup gnn_hessian_doc Hessian Matrix Evaluation Routines.
 * @ingroup libgnn_node
 *
 * The routines presented in this module deal with the computation of the
 * second derivatives of the error, given by
 *
 * \f[ \frac{\partial^2 E}{\partial w_i \partial w_j} \f]
 *
 * where both indexes run over all parameters \f$i,j=1, \ldots, l\f$. These
 * elements form the Hessian Matrix, which plays a very important role in
 * several aspects of neural computing, like, per example:
 * - Some parameter optimization algorithms use the information of the
 *   error surface's second-order properties.
 * - Pruning algorithms use the Hessian in order to determine the relevance
 *   of the parameters.
 * Although the Hessian Matrix does hold important information about a
 * modelling function, its computation can be very expensive, both in
 * time and memory requirements.
 *
 * It is important to note that some of these algorithms assume special
 * error functions. Nevertheless, sometimes this assumption does not
 * heavilly contradict the situation encountered in e.g. pruning procedures
 * using another cost function.
 */



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

#include <gsl/gsl_blas.h>
#include "gnn_evaluation.h"
#include "gnn_hessian.h"
#include <math.h>



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

/**
 * @brief Fast multiplication by the Hessian matrix.
 * @ingroup gnn_hessian_doc
 *
 * Some applications do not require the Hessian explicitally. Instead,
 * what really is needed is a fast method for computing \f$Hv\f$, where
 * \f$H\f$ is the Hessian matrix (which is assumed to be symmetric, so that
 * \f$Hv = v^TH\f$) and
 * \f$v\f$ is a vector. This numerical method computes this multiplication
 * using central differences:
 *
 * \f[ v^T H
 *     = \frac{1}{2\epsilon}
 *       \left [
 *           \frac{\partial E}{\partial w}(w + \epsilon v) -
 *           \frac{\partial E}{\partial w}(w - \epsilon v)
 *       \right ]
 * \f]
 *
 * where \f$\epsilon\f$ is a small value that gives the required precision,
 * since the error resulting from this approximation is \f$O(\epsilon^2)\f$.
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  v       A pointer to the vector \f$v\f$.
 * @param  vH      A pointer to a vector of size \f$l \times 1\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting multiplication
 *                 will be stored in this vector.
 * @param  eps     The required precision \f$\epsilon\f$.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_vector_mult (gnn_line *line,
                    const gsl_vector *v, gsl_vector *vH, double eps)
{
    size_t j;
    gnn_grad *grad;

    assert (line != NULL);
    assert (v    != NULL);
    assert (vH   != NULL);
    assert (v    != vH);
    
    /* get gnn_grad */
    grad = gnn_line_get_grad (line);

    /* check sizes */
    if (grad->l != v->size)
    {
        GSL_ERROR ("the vector's size should match the number of parameters",
                   GSL_EINVAL);
    }
    if (grad->l != vH->size)
    {
        GSL_ERROR ("the result's size should match the number of parameters",
                   GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* check precision */
    if (eps <= 0.0)
    {
        GSL_ERROR ("eps should be stricly positive", GSL_EINVAL);
    }

    /* set direction */
    gnn_line_set_direction (line, v);

    /* compute gradient step forward */
    gnn_line_all (line, eps, gnnLineDE);
    gsl_vector_memcpy (vH, GNN_GRAD_DW (grad));

    /* compute gradient step backward */
    gnn_line_all (line, -eps, gnnLineDE);
    gsl_vector_sub (vH, GNN_GRAD_DW (grad));

    /* finish approximation */
    gsl_vector_scale (vH, 1 / (2*eps));

    /* restore original parameters */
    gnn_line_all (line, 0.0, gnnLineE);

    return 0;
}

/**
 * @brief Finite differences diagonal Hessian matrix approximation.
 * @ingroup gnn_hessian_doc
 *
 * This numerical efficient method computes the diagonal values of the Hessian
 * matrix using the central differentes method, given by
 *
 * \f[ \frac{\partial^2 E}{\partial w_i^2}
 *     = \frac{1}{2\epsilon}
 *       \left [
 *           \frac{\partial E}{\partial w_i}(w_i + \epsilon) -
 *           \frac{\partial E}{\partial w_i}(w_i - \epsilon)
 *       \right ]
 * \f]
 *
 * where \f$\epsilon\f$ is a small value that gives the required precision,
 * since the error resulting from this approximation is \f$O(\epsilon^2)\f$.
 *
 * The same as above written in matrix notation is:
 *
 * \f[ H = \frac{1}{2\epsilon}
 *       \left [ g(w + \epsilon I) - g(w - \epsilon I) \right ]
 * \f]
 *
 * where \f$H\f$ is the Hessian matrix and
 * \f$g(w_0) = \frac{\partial E}{\partial w}(w_0)\f$. Note that this
 * requires a single forward and backward step for evaluation, since
 * the assumption of the diagonal Hessian is equivalent to an uncorrelated
 * parameter space.
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  H       A pointer to a vector of size \f$l \times 1\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting diagonal terms of the
 *                 Hessian will be stored in this vector.
 * @param  eps     The required precision \f$\epsilon\f$.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_diagonal (gnn_line *line, gsl_vector *H, double eps)
{
    gnn_grad   *grad;
    gsl_vector *d;
    
    assert (line != NULL);
    assert (H    != NULL);
    
    /* get gnn_grad */
    grad = gnn_line_get_grad (line);
    
    /* check sizes */
    if (grad->l != H->size)
    {
        GSL_ERROR ("the Hessian's size should be l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* check precision */
    if (eps <= 0.0)
    {
        GSL_ERROR ("eps should be stricly positive", GSL_EINVAL);
    }

    /* compute direction */
    d = (gsl_vector *) gnn_line_get_direction (line);
    gsl_vector_set_all (d, 1.0);

    /* compute gradient step forward */
    gnn_line_all (line, eps, gnnLineDE);
    gsl_vector_memcpy (H, GNN_GRAD_DW (grad));

    /* compute gradient step backward */
    gnn_line_all (line, -eps, gnnLineDE);
    gsl_vector_sub (H, GNN_GRAD_DW (grad));

    /* finish approximation */
    gsl_vector_scale (H, 1 / (2*eps));

    /* restore original parameters */
    gnn_line_all (line, 0.0, gnnLineE);

    return 0;
}

/**
 * @brief Initializes the Levenberg-Marquadt Hessian approximation.
 * @ingroup gnn_hessian_doc
 *
 * This function initializes the Hessian for the Levenberg-Marquadt
 * approximation.
 *
 * For further details, see \ref gnn_hessian_levenberg_marquadt.
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  H       A pointer to a matrix of size \f$l \times l\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting Hessian will be
 *                 stored in this matrix.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_levenberg_marquadt_init (gnn_grad *grad, gsl_matrix *H)
{
    assert (grad != NULL);
    assert (H    != NULL);

    /* check sizes */
    if (grad->l != H->size1 || grad->l != H->size2)
    {
        GSL_ERROR ("the Hessian's size should be l x l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* clear Hessian */
    gsl_matrix_set_zero (H);

    return 0;
}

/**
 * @brief Levenberg-Marquadt Hessian approximation.
 * @ingroup gnn_hessian_doc
 *
 * This function computes the k-th iteration of the Levenberg-Marquadt Hessian
 * approximation, using the k-th pattern.
 *
 * The Levenberg-Marquadt approximation (also known as the outer product
 * approximation) of the Hessian matrix is given by
 *
 * \f[ \frac{\partial^2 E}{\partial w_i \partial w_j}
 *     = \sum_{k=1}^P \frac{\partial y^k}{\partial w_i}
 *                    \frac{\partial y^k}{\partial w_j}
 * \f]
 *
 * where the sum runs over all pattern outputs \f$y^k\f$, \f$k=1,\ldots,P\f$.
 * It is very important to note that this approximation is only valid for
 * an error function of the form
 * \f[ E = \frac{1}{2} \sum_{k=1}^P (y^k - t^k)^2, \f]
 * that is, for an sum-of-squares error function (and its equivalent forms).
 *
 * The previous formula can be written in matrix form, giving the sequential
 * procedure for building the Hessian:
 *
 * \f[ H_P = \sum_{k=1}^P \frac{\partial E}{\partial w}^k
 *                        \frac{\partial E}{\partial w}^k
 * \f]
 *
 * which is the way used by \ref libgnn.
 *
 * @warning Evaluation??? Shouldn't it be handled completely by a gnn_eval
 *          or gnn_grad??
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  H       A pointer to a matrix of size \f$l \times l\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting Hessian will be
 *                 stored in this matrix.
 * @param  k       The index of the pattern to be used.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_levenberg_marquadt (gnn_grad *grad, gsl_matrix *H, size_t k)
{
    double p;
    gsl_vector *x;
    gsl_vector *t;
    gsl_matrix_view G;

    assert (grad != NULL);
    assert (H    != NULL);

    /* check sizes */
    if (grad->l != H->size1 || grad->l != H->size2)
    {
        GSL_ERROR ("the Hessian's size should be l x l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* check batch bounds */
    if (k < 0 || k >= grad->P)
    {
        GSL_ERROR ("pattern index out of bounds", GSL_EINVAL);
    }

    /* get pattern */
    gnn_dataset_get (grad->data, k, &x, &t, &p);

    /* compute gradient */
    gnn_node_evaluate_init (grad->node);
    gnn_node_evaluate_f (grad->node, x, grad->y);
    grad->e = gnn_criterion_evaluate_e (grad->crit, grad->y, t);
    gnn_criterion_evaluate_dy (grad->crit, grad->dy);
    gnn_node_evaluate_dx (grad->node, grad->dy, grad->dx);
    gnn_node_evaluate_dw (grad->node, grad->dw);

    /* scale dw gradient */
    gsl_vector_scale (grad->dw, p);

    /* get matrix view from gradient */
    G = gsl_matrix_view_vector (grad->dw, grad->l, 1);

    /* compute g*gT and sum to Hessian */
    gsl_blas_dgemm (CblasNoTrans, CblasTrans,
                    1.0, &(G.matrix), &(G.matrix), 1.0, H);

    return 0;
}

/**
 * @brief Levenberg-Marquadt Hessian Inverse approximation initialization.
 * @ingroup gnn_hessian_doc
 *
 * This function initializes the outer product inverse approximation routine.
 * For further details, please refer to \ref gnn_hessian_inverse.
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  H       A pointer to a matrix of size \f$l \times l\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting inverse Hessian will be
 *                 stored in this matrix.
 * @param  alpha   A small value for \f$\alpha\f$.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_inverse_init (gnn_grad *grad, gsl_matrix *H, double alpha)
{
    assert (grad != NULL);
    assert (H    != NULL);

    /* check sizes */
    if (grad->l != H->size1 || grad->l != H->size2)
    {
        GSL_ERROR ("the Hessian's size should be l x l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* clear Hessian */
    gsl_matrix_set_identity (H);
    gsl_matrix_scale (H, alpha);

    return 0;
}

/**
 * @brief Levenberg-Marquadt Hessian Inverse approximation.
 * @ingroup gnn_hessian_doc
 *
 * This function computes the k-th iteration of the outer-product inverse
 * Hessian approximation, using the k-th pattern.
 *
 * The outer product approximation used in \ref gnn_hessian_levenberg_marquadt
 * can be used to develop an efficient procedure for computing the inverse
 * of the Hessian matrix.
 *
 * Recalling the formula for the outer product approximation
 *
 * \f[ H_P = \sum_{k=1}^P g^k (g^k)^T \f]
 *
 * where \f$g^k = \frac{\partial E}{\partial w}^k\f$, it is easy to see that
 * this leads to a sequential form for computing the Hessian:
 *
 * \f[ H_{N+1} = H_{N} + g^{N+1} (g^{N+1})^T \f]
 *
 * Using the following matrix identity
 *
 * \f[ (A + BC)^{-1} = A^{-1} - A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1} \f]
 *
 * where \f$I\f$ is the identity matrix, and setting \f$H_{N} = A\f$,
 * \f$g^{N+1}=B\f$ and \f$(g^{N+1})^T=C\f$, then the previous formula
 * can be rewritten in order to obtain an iterative procedure for building
 * the inverse of the Hessian:
 *
 * \f[
 * H^{-1}_{N+1} = H^{-1}_N -
 *     \frac{ H^{-1}_N g^{N+1} (g^{N+1})^T H^{-1}_N }
 *          {    1 + (g^{N+1})^T H^{-1}_N g^{N+1}   }
 * \f]
 *
 * The procedure can be initialized by \f$ H_0 = \alpha I \f$, where
 * \f$\alpha\f$ is a small value (because the procedure actually seeks
 * to approximate \f$ H + \alpha I \f$).
 *
 * It is important to note that, as in the case for
 * \ref gnn_hessian_levenberg_marquadt, this approximation is only valid
 * for a sum-of-squares error function.
 *
 * @warning Evaluation??? Shouldn't it be handled completely by a gnn_eval
 *          or gnn_grad??
 *
 * @param  grad    A pointer to a \ref gnn_grad buffer.
 * @param  H       A pointer to a matrix of size \f$l \times l\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting inverse Hessian will be
 *                 stored in this matrix.
 * @param  k       The index of the pattern to be used.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_inverse (gnn_grad *grad, gsl_matrix *H, size_t k)
{
    double p;
    gsl_vector *x;
    gsl_vector *t;
    gsl_matrix_view Hg;

    assert (grad != NULL);
    assert (H    != NULL);

    /* check sizes */
    if (grad->l != H->size1 || grad->l != H->size2)
    {
        GSL_ERROR ("the Hessian's size should be l x l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* check batch bounds */
    if (k < 0 || k >= grad->P)
    {
        GSL_ERROR ("pattern index out of bounds", GSL_EINVAL);
    }

    /* get pattern */
    gnn_dataset_get (grad->data, k, &x, &t, &p);

    /* compute gradient */
    gnn_node_evaluate_init (grad->node);
    gnn_node_evaluate_f (grad->node, x, grad->y);
    grad->e = gnn_criterion_evaluate_e (grad->crit, grad->y, t);
    gnn_criterion_evaluate_dy (grad->crit, grad->dy);
    gnn_node_evaluate_dx (grad->node, grad->dy, grad->dx);
    gnn_node_evaluate_dw (grad->node, grad->dw);

    /* scale dw gradient */
    gsl_vector_scale (grad->dw, p);

    /* COMPUTE FORMULA */
    {
        double v;
        
        /* h = H^-1 * g */
        gsl_blas_dgemv (CblasNoTrans, 1.0, H, grad->dw, 0.0, grad->h);
        
        /* get matrix view from vector */
        Hg = gsl_matrix_view_vector (grad->h, grad->l, 1);
        
        /* v = 1 + g^T * h */
        gsl_blas_ddot (grad->dw, grad->dw, &v);
        v = 1.0 + v;
        
        /* compute (h*h^T / v) and sum to Hessian */
        gsl_blas_dgemm (CblasNoTrans, CblasTrans,
                        1.0, &(Hg.matrix), &(Hg.matrix), -1.0 / v, H);
    }
    
    return 0;
}

/**
 * @brief Finite differences Hessian matrix approximation.
 * @ingroup gnn_hessian_doc
 *
 * This numerical efficient method computes the values of the Hessian matrix
 * using the central differentes method, given by
 *
 * \f[ \frac{\partial^2 E}{\partial w_i \partial w_j}
 *     = \frac{1}{2\epsilon}
 *       \left [
 *           \frac{\partial E}{\partial w_i}(w_j + \epsilon) -
 *           \frac{\partial E}{\partial w_i}(w_j - \epsilon)
 *       \right ]
 * \f]
 *
 * where \f$\epsilon\f$ is a small value that gives the required precision,
 * since the error resulting from this approximation is \f$O(\epsilon^2)\f$.
 *
 * @param  line    A pointer to a \ref gnn_line buffer.
 * @param  H       A pointer to a matrix of size \f$l \times l\f$,
 *                 where \f$l\f$ is the number of parameters involved
 *                 in the \ref gnn_node. The resulting Hessian will be
 *                 stored in this matrix.
 * @param  eps     The required precision \f$\epsilon\f$.
 * @return Returns 0 if succeeded.
 */
int
gnn_hessian_finite_differences (gnn_line *line, gsl_matrix *H, double eps)
{
    size_t j;
    gnn_grad *grad;

    assert (line != NULL);
    assert (H    != NULL);
    
    /* get gnn_grad */
    grad = gnn_line_get_grad (line);

    /* check sizes */
    if (grad->l != H->size1 || grad->l != H->size2)
    {
        GSL_ERROR ("the Hessian's size should be l x l, with l the number "
                   "of parameters", GSL_EINVAL);
    }
    if (grad->P == GNN_INFTY)
    {
        GSL_ERROR ("the dataset should be finite", GSL_EINVAL);
    }

    /* check precision */
    if (eps <= 0.0)
    {
        GSL_ERROR ("eps should be stricly positive", GSL_EINVAL);
    }

    /* compute Hessian */
    for (j=0; j<grad->l; ++j)
    {
        gsl_vector *d;
        gsl_vector_view Hj;
        
        /* compute direction */
        d = (gsl_vector *) gnn_line_get_direction (line);
        gsl_vector_set_all (d, 0.0);
        gsl_vector_set (d, j, 1.0);

        /* get Hessian row view */
        Hj = gsl_matrix_row (H, j);

        /* compute gradient step forward */
        gnn_line_all (line, eps, gnnLineDE);
        gsl_vector_memcpy (&(Hj.vector), GNN_GRAD_DW (grad));

        /* compute gradient step backward */
        gnn_line_all (line, -eps, gnnLineDE);
        gsl_vector_sub (&(Hj.vector), GNN_GRAD_DW (grad));

        /* finish approximation */
        gsl_vector_scale (&(Hj.vector), 1 / (2*eps));
    }

    /* recover original parameters */
    gnn_line_all (line, 0.0, gnnGradE);

    return 0;
}

---