---

gnn_gcomm.c

Go to the documentation of this file.

/***************************************************************************
 *  @file gnn_gcomm.c
 *  @brief gnn_gcomm Implementation.
 *
 *  @date   : 28-09-03 21:09
 *  @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_gcomm_doc gnn_gcomm : Generalized Committee Convergence Node.
 * @ingroup gnn_atomic_doc
 *
 * This node type implements the function given by:
 *
 * \f[ y_j = \sum_k \alpha_k x_{km+j}, \qquad 0 \leq j \leq m-1 \f]
 *
 * where \f$k\f$ and \f$j\f$ are uniquely combined in order to run over all
 * input indexes \f$i = km+j = 0,\ldots, n-1\f$, and the coefficients
 * \f$\alpha_k\f$ are:
 *
 * \f[ \alpha_k = \frac{w_k^2}{\sum_{k'} w_{k'}^2 } \f]
 *
 * Note that this assures that
 *
 * \f[ \sum_k \alpha_k = 1,\qquad 0 \leq \alpha_k \leq 1 \f]
 *
 * In other words, this function slices the input vector \f$x\f$ into several
 * subvectors of size \f$m \times 1\f$ and performs a weighted sum of them.
 *
 * If the input vector's
 * size isn't a multiple of the output vector's size,
 * that is, if \f$n\f$ can't be written as \f$km\f$, then the last vector slice
 * will consider only the remaining terms, e.g.
 * if \f$x=[x_0, \ldots, x_5]^T\f$ and \f$y\f$ is of size
 * \f$4 \times 1\f$ then
 *
 * \f[ y = \left [ \begin{array}{c}
 *                     \alpha_0 x_0 + \alpha_1 x_4 \\
 *                     \alpha_0 x_1 + \alpha_1 x_5 \\
 *                     \alpha_0 x_2 \\
 *                     \alpha_0 x_3 \\
 *         \end{array} \right ]
 * \f]
 *
 * Schematically, this idea can be depicted as
 *
 * <img src="images/gnn_gcomm1.png">
 *
 * The function's gradients are:
 * \f[ \frac{\partial E}{\partial x_i}
 *     = \frac{\partial E}{\partial x_{km+j}}
 *     = \alpha_k \frac{\partial E}{\partial y_j}
 * \f]
 * and for the \f$w_k\f$:
 * \f[ \frac{\partial E}{\partial w_k}
 *     = \sum_j \frac{\partial E}{\partial y_j}
 *              \sum_{k'} \frac{\partial \alpha_{k'}}{\partial w_k} x_{k'm+j}
 * \f]
 * and
 * \f[ \frac{\partial \alpha_{k'}}{\partial w_k}
 *     = \left \{ \begin{array}{ccl}
 *          \frac{2w_k \sum_i w_i^2 - w_k^2 2w_k}
 *               { \left ( \sum_i w_i^2 \right )^2 } & , & k = k' \\
 *         -\frac{w_{k'}^2 2w_k}
 *               { \left ( \sum_i w_i^2 \right )^2 } & , & k \neq k' \\
 *       \end{array} \right \}
 *     = 2 w_k \left ( \delta^{k'}_k - \frac{\alpha_{k'}}{\Sigma} \right )
 * \f]
 * where \f$\Sigma = \sum_i w_i^2 \f$.
 */


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

#include <gsl/gsl_blas.h>
#include "gnn_gcomm.h"



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

static int
gnn_gcomm_f (gnn_node *node,
             const gsl_vector *x,
             const gsl_vector *w,
             gsl_vector *y);

static int
gnn_gcomm_dx (gnn_node *node,
              const gsl_vector *x,
              const gsl_vector *w,
              const gsl_vector *dy,
              gsl_vector *dx);

static int
gnn_gcomm_dw (gnn_node *node,
              const gsl_vector *x,
              const gsl_vector *w,
              const gsl_vector *dy,
              gsl_vector *dw);

static void
gnn_gcomm_destroy (gnn_node *node);



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

/**
 * @brief Computes the output.
 * @ingroup gnn_gcomm_doc
 *
 * @param   node A pointer to the \ref gnn_gcomm node.
 * @param x The input vector \f$x\f$.
 * @param w The parameter vector \f$w\f$.
 * @param y The output buffer vector \f$y\f$.
 * @return Returns 0 on success.
 */
static int
gnn_gcomm_f (gnn_node *node,
             const gsl_vector *x,
             const gsl_vector *w,
             gsl_vector *y)
{
    size_t k;
    gnn_gcomm *gnode;

    /* get views */
    gnode = (gnn_gcomm *) node;

    /* copy input */
    gsl_vector_memcpy (gnode->x, x);

    /* evaluate */
    gnode->sumw = gsl_blas_dnrm2 (w);
    for (k=0; k<gnode->rep; ++k)
    {
        double wk;
        double alphak;
        
        /* compute w_k^2 */
        wk = gsl_vector_get (w, k);
        wk = wk * wk;
        
        /* compute alpha_k*/
        alphak = wk / gnode->sumw;
        gsl_vector_set (gnode->alpha, k, alphak);
        gsl_blas_daxpy (alphak, gnode->X[k], y);
    }

    return 0;
}

/**
 * @brief Computes \f$ \frac{\partial E}{\partial x} \f$.
 * @ingroup gnn_gcomm_doc
 *
 * @param    node A pointer to the \ref gnn_gcomm node.
 * @param x  The input vector \f$x\f$.
 * @param w  The parameter vector \f$w\f$.
 * @param dy The vector \f$\frac{\partial E}{\partial y}\f$.
 * @param dx The output buffer vector \f$\frac{\partial E}{\partial x}\f$.
 * @return Returns 0 on success.
 */
static int
gnn_gcomm_dx (gnn_node *node,
              const gsl_vector *x,
              const gsl_vector *w,
              const gsl_vector *dy,
              gsl_vector *dx)
{
    size_t k;
    gnn_gcomm *gnode;

    /* get views */
    gnode = (gnn_gcomm *) node;

    /* evaluate */
    for (k=0; k<gnode->rep; ++k)
    {
        double alphak;
        
        alphak = gsl_vector_get (gnode->alpha, k);
        
        gsl_vector_memcpy (gnode->X[k], dy);
        gsl_vector_scale (gnode->X[k], alphak);
    }

    /* copy buffer to dx */
    gsl_vector_memcpy (dx, gnode->x);

    return 0;
}

/**
 * @brief Computes \f$ \frac{\partial E}{\partial w} \f$.
 * @ingroup gnn_gcomm_doc
 *
 * @param    node A pointer to the \ref gnn_gcomm node.
 * @param x  The input vector \f$x\f$.
 * @param w  The parameter vector \f$w\f$.
 * @param dy The vector \f$\frac{\partial E}{\partial y}\f$.
 * @param dx The output buffer vector \f$\frac{\partial E}{\partial w}\f$.
 * @return Returns 0 on success.
 */
static int
gnn_gcomm_dw (gnn_node *node,
              const gsl_vector *x,
              const gsl_vector *w,
              const gsl_vector *dy,
              gsl_vector *dw)
{
    size_t k;
    size_t j;
    size_t kp;
    gnn_gcomm *gnode;

    /* get views */
    gnode = (gnn_gcomm *) node;

    /* copy input vector into sliced buffers */
    gsl_vector_set_zero (gnode->x);
    gsl_vector_memcpy (gnode->x, x);

    /* evaluate */
    for (k=0; k<dw->size; ++k)
    {
        double dwk;
        double wk;
        double Sj;
        
        wk  = gsl_vector_get (w, k);
        dwk = gsl_vector_get (dw, k);
        Sj  = 0.0;
        
        for (j=0; j<dy->size; ++j)
        {
            double dyj;
            double Skp;

            Skp = 0.0;
            dyj = gsl_vector_get (dy, j);
            
            for (kp=0; kp<w->size; ++kp)
            {
                size_t i;
                double alphakp;
                double dalphakp;
                double xi;
                
                i = kp * dy->size + j;
                xi = gsl_vector_get (gnode->x, i);
                
                alphakp  = gsl_vector_get (gnode->alpha, kp);
                dalphakp = (k == kp)? 1.0 - alphakp / gnode->sumw
                                    :     - alphakp / gnode->sumw;

                Skp += dalphakp * xi;
            }

            Sj += dyj * 2 * wk * Skp;
        }
        
        dwk += Sj;
        gsl_vector_set (dw, k, Sj);
    }

    return 0;
}

/**
 * @brief Computes the output.
 * @ingroup gnn_gcomm_doc
 *
 * @param node A pointer to the \ref gnn_gcomm node.
 */
static void
gnn_gcomm_destroy (gnn_node *node)
{
    gnn_gcomm *gnode;

    gnode = (gnn_gcomm *) node;

    if (gnode->alpha != NULL)
        gsl_vector_free (gnode->alpha);
    if (gnode->xbuf != NULL)
        gsl_vector_free (gnode->xbuf);
    if (gnode->X_view != NULL)
        free (gnode->X_view);
    if (gnode->X != NULL)
        free (gnode->X);
}



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

/**
 * @brief Creates an generalized committee node.
 * @ingroup gnn_gcomm_doc
 *
 * This function creates a node of the \ref gnn_gcomm type. The coefficients
 * are all initialized at the same value. This is equivalent to set the
 * \f$a_k\f$ all to \f$\frac{1}{m}\f$.
 *
 * @param input_size  The input size \f$n\f$.
 * @param output_size The output size \f$m\f$.
 * @return A pointer to a new \ref gnn_gcomm node.
 */
gnn_node *
gnn_gcomm_new (size_t input_size, size_t output_size)
{
    int status;
    size_t k;
    gsl_vector *w;
    gnn_node *node;
    gnn_gcomm *gnode;
    size_t param_size;

    /* check sizes */
    if (output_size < 1)
    {
        GSL_ERROR_VAL ("output size should be strictly positive",
                       GSL_EINVAL, NULL);
    }
    if (input_size < output_size)
    {
        GSL_ERROR_VAL ("input size should be greater than the output size",
                       GSL_EINVAL, NULL);
    }

    /* allocate the node */
    gnode = (gnn_gcomm *) malloc (sizeof (gnn_gcomm));
    if (gnode == NULL)
    {
        GSL_ERROR_VAL ("could not allocate memory for gnn_gcomm node",
                       GSL_ENOMEM, NULL);
    }

    /* get view as a gnn_node */
    node = (gnn_node *) gnode;

    /* Initialize the node */
    status = gnn_node_init (node,
                            "gnn_gcomm",
                            gnn_gcomm_f,
                            gnn_gcomm_dx,
                            gnn_gcomm_dw,
                            gnn_gcomm_destroy);
    if (status)
    {
        GSL_ERROR_VAL ("could not initialize gnn_gcomm node",
                       GSL_EFAILED, NULL);
    }

    /* compute amount of replications */
    gnode->rep = (input_size-1) / output_size + 1;
    
    /* set sizes */
    param_size = gnode->rep;
    status = gnn_node_set_sizes (node, input_size, output_size, param_size);
    if (status)
    {
        GSL_ERROR_VAL ("could not set sizes for gnn_gcomm node",
                       GSL_EFAILED, NULL);
    }

    /* create alpha buffer */
    gnode->alpha = gsl_vector_calloc (param_size);
    
    /* create full input vector buffer and its view */
    gnode->xbuf  = gsl_vector_calloc (gnode->rep * output_size);

    gnode->X_view =
        (gsl_vector_view *) malloc (sizeof (gsl_vector_view) * gnode->rep);

    gnode->X =
        (gsl_vector **) malloc (sizeof (gsl_vector *) * gnode->rep);

    if (gnode->xbuf == NULL || gnode->X_view == NULL || gnode->X == NULL)
    {
        gnn_node_destroy (node);
        GSL_ERROR_VAL ("could not allocate buffers for gnn_gcomm node",
                       GSL_ENOMEM, NULL);
    }

    /* create input vector slices : views and pointers */
    gnode->x_view = gsl_vector_subvector (gnode->xbuf, 0, input_size);
    gnode->x      = &(gnode->x_view.vector);

    for (k=0; k<gnode->rep; ++k)
    {
        gnode->X_view[k] =
            gsl_vector_subvector (gnode->xbuf, k * output_size, output_size);
        gnode->X[k] = &(gnode->X_view[k].vector);

    }

    /* set initial parameter vector */
    gsl_vector_set_all (gnode->alpha, 1.0);
    gnn_node_param_set (node, gnode->alpha);

    return node;
}

---