htswanalysis/MACS/lib/gsl/gsl-1.11/multimin/vector_bfgs.c

   1 /* multimin/vector_bfgs.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* vector_bfgs.c -- Limited memory Broyden-Fletcher-Goldfarb-Shanno method */
  21
  22 /* Modified by Brian Gough to use single iteration structure */
  23
  24 #include <config.h>
  25 #include <gsl/gsl_multimin.h>
  26 #include <gsl/gsl_blas.h>
  27
  28 #include "directional_minimize.c"
  29
  30 typedef struct
  31 {
  32   int iter;
  33   double step;
  34   double max_step;
  35   double tol;
  36   gsl_vector *x1;
  37   gsl_vector *dx1;
  38   gsl_vector *x2;
  39   double g0norm;
  40   double pnorm;
  41   gsl_vector *p;
  42   gsl_vector *x0;
  43   gsl_vector *g0;
  44   gsl_vector *dx0;
  45   gsl_vector *dg0;
  46 }
  47 vector_bfgs_state_t;
  48
  49 static int
  50 vector_bfgs_alloc (void *vstate, size_t n)
  51 {
  52   vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;
  53
  54   state->x1 = gsl_vector_calloc (n);
  55
  56   if (state->x1 == 0)
  57     {
  58       GSL_ERROR ("failed to allocate space for x1", GSL_ENOMEM);
  59     }
  60
  61   state->dx1 = gsl_vector_calloc (n);
  62
  63   if (state->dx1 == 0)
  64     {
  65       gsl_vector_free (state->x1);
  66       GSL_ERROR ("failed to allocate space for dx1", GSL_ENOMEM);
  67     }
  68
  69   state->x2 = gsl_vector_calloc (n);
  70
  71   if (state->x2 == 0)
  72     {
  73       gsl_vector_free (state->dx1);
  74       gsl_vector_free (state->x1);
  75       GSL_ERROR ("failed to allocate space for x2", GSL_ENOMEM);
  76     }
  77
  78   state->p = gsl_vector_calloc (n);
  79
  80   if (state->p == 0)
  81     {
  82       gsl_vector_free (state->x2);
  83       gsl_vector_free (state->dx1);
  84       gsl_vector_free (state->x1);
  85       GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
  86     }
  87
  88   state->x0 = gsl_vector_calloc (n);
  89
  90   if (state->x0 == 0)
  91     {
  92       gsl_vector_free (state->p);
  93       gsl_vector_free (state->x2);
  94       gsl_vector_free (state->dx1);
  95       gsl_vector_free (state->x1);
  96       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
  97     }
  98
  99   state->g0 = gsl_vector_calloc (n);
 100
 101   if (state->g0 == 0)
 102     {
 103       gsl_vector_free (state->x0);
 104       gsl_vector_free (state->p);
 105       gsl_vector_free (state->x2);
 106       gsl_vector_free (state->dx1);
 107       gsl_vector_free (state->x1);
 108       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 109     }
 110
 111   state->dx0 = gsl_vector_calloc (n);
 112
 113   if (state->dx0 == 0)
 114     {
 115       gsl_vector_free (state->g0);
 116       gsl_vector_free (state->x0);
 117       gsl_vector_free (state->p);
 118       gsl_vector_free (state->x2);
 119       gsl_vector_free (state->dx1);
 120       gsl_vector_free (state->x1);
 121       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 122     }
 123
 124   state->dg0 = gsl_vector_calloc (n);
 125
 126   if (state->dg0 == 0)
 127     {
 128       gsl_vector_free (state->dx0);
 129       gsl_vector_free (state->g0);
 130       gsl_vector_free (state->x0);
 131       gsl_vector_free (state->p);
 132       gsl_vector_free (state->x2);
 133       gsl_vector_free (state->dx1);
 134       gsl_vector_free (state->x1);
 135       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 136     }
 137
 138   return GSL_SUCCESS;
 139 }
 140
 141 static int
 142 vector_bfgs_set (void *vstate, gsl_multimin_function_fdf * fdf,
 143                  const gsl_vector * x, double *f, gsl_vector * gradient,
 144                  double step_size, double tol)
 145 {
 146   vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;
 147
 148   state->iter = 0;
 149   state->step = step_size;
 150   state->max_step = step_size;
 151   state->tol = tol;
 152
 153   GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
 154
 155   /* Use the gradient as the initial direction */
 156
 157   gsl_vector_memcpy (state->x0, x);
 158   gsl_vector_memcpy (state->p, gradient);
 159   gsl_vector_memcpy (state->g0, gradient);
 160
 161   {
 162     double gnorm = gsl_blas_dnrm2 (gradient);
 163     state->pnorm = gnorm;
 164     state->g0norm = gnorm;
 165   }
 166
 167   return GSL_SUCCESS;
 168 }
 169
 170 static void
 171 vector_bfgs_free (void *vstate)
 172 {
 173   vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;
 174
 175   gsl_vector_free (state->dg0);
 176   gsl_vector_free (state->dx0);
 177   gsl_vector_free (state->g0);
 178   gsl_vector_free (state->x0);
 179   gsl_vector_free (state->p);
 180   gsl_vector_free (state->x2);
 181   gsl_vector_free (state->dx1);
 182   gsl_vector_free (state->x1);
 183 }
 184
 185 static int
 186 vector_bfgs_restart (void *vstate)
 187 {
 188   vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;
 189
 190   state->iter = 0;
 191   return GSL_SUCCESS;
 192 }
 193
 194 static int
 195 vector_bfgs_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
 196                      gsl_vector * x, double *f,
 197                      gsl_vector * gradient, gsl_vector * dx)
 198 {
 199   vector_bfgs_state_t *state = (vector_bfgs_state_t *) vstate;
 200
 201   gsl_vector *x1 = state->x1;
 202   gsl_vector *dx1 = state->dx1;
 203   gsl_vector *x2 = state->x2;
 204   gsl_vector *p = state->p;
 205   gsl_vector *g0 = state->g0;
 206   gsl_vector *x0 = state->x0;
 207
 208   double pnorm = state->pnorm;
 209   double g0norm = state->g0norm;
 210
 211   double fa = *f, fb, fc;
 212   double dir;
 213   double stepa = 0.0, stepb, stepc = state->step, tol = state->tol;
 214
 215   double g1norm;
 216   double pg;
 217
 218   if (pnorm == 0.0 || g0norm == 0.0)
 219     {
 220       gsl_vector_set_zero (dx);
 221       return GSL_ENOPROG;
 222     }
 223
 224   /* Determine which direction is downhill, +p or -p */
 225
 226   gsl_blas_ddot (p, gradient, &pg);
 227
 228   dir = (pg >= 0.0) ? +1.0 : -1.0;
 229
 230   /* Compute new trial point at x_c= x - step * p, where p is the
 231      current direction */
 232
 233   take_step (x, p, stepc, dir / pnorm, x1, dx);
 234
 235   /* Evaluate function and gradient at new point xc */
 236
 237   fc = GSL_MULTIMIN_FN_EVAL_F (fdf, x1);
 238
 239   if (fc < fa)
 240     {
 241       /* Success, reduced the function value */
 242       state->step = stepc * 2.0;
 243       *f = fc;
 244       gsl_vector_memcpy (x, x1);
 245       GSL_MULTIMIN_FN_EVAL_DF (fdf, x1, gradient);
 246       return GSL_SUCCESS;
 247     }
 248
 249 #ifdef DEBUG
 250   printf ("got stepc = %g fc = %g\n", stepc, fc);
 251 #endif
 252
 253   /* Do a line minimisation in the region (xa,fa) (xc,fc) to find an
 254      intermediate (xb,fb) satisifying fa > fb < fc.  Choose an initial
 255      xb based on parabolic interpolation */
 256
 257   intermediate_point (fdf, x, p, dir / pnorm, pg,
 258                       stepa, stepc, fa, fc, x1, dx1, gradient, &stepb, &fb);
 259
 260   if (stepb == 0.0)
 261     {
 262       return GSL_ENOPROG;
 263     }
 264
 265   minimize (fdf, x, p, dir / pnorm,
 266             stepa, stepb, stepc, fa, fb, fc, tol,
 267             x1, dx1, x2, dx, gradient, &(state->step), f, &g1norm);
 268
 269   gsl_vector_memcpy (x, x2);
 270
 271   /* Choose a new direction for the next step */
 272
 273   state->iter = (state->iter + 1) % x->size;
 274
 275   if (state->iter == 0)
 276     {
 277       gsl_vector_memcpy (p, gradient);
 278       state->pnorm = g1norm;
 279     }
 280   else
 281     {
 282       /* This is the BFGS update: */
 283       /* p' = g1 - A dx - B dg */
 284       /* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */
 285       /* B = dx.g/dx.dg */
 286
 287       gsl_vector *dx0 = state->dx0;
 288       gsl_vector *dg0 = state->dg0;
 289
 290       double dxg, dgg, dxdg, dgnorm, A, B;
 291
 292       /* dx0 = x - x0 */
 293       gsl_vector_memcpy (dx0, x);
 294       gsl_blas_daxpy (-1.0, x0, dx0);
 295
 296       /* dg0 = g - g0 */
 297       gsl_vector_memcpy (dg0, gradient);
 298       gsl_blas_daxpy (-1.0, g0, dg0);
 299
 300       gsl_blas_ddot (dx0, gradient, &dxg);
 301       gsl_blas_ddot (dg0, gradient, &dgg);
 302       gsl_blas_ddot (dx0, dg0, &dxdg);
 303
 304       dgnorm = gsl_blas_dnrm2 (dg0);
 305
 306       if (dxdg != 0)
 307         {
 308           B = dxg / dxdg;
 309           A = -(1.0 + dgnorm * dgnorm / dxdg) * B + dgg / dxdg;
 310         }
 311       else
 312         {
 313           B = 0;
 314           A = 0;
 315         }
 316
 317       gsl_vector_memcpy (p, gradient);
 318       gsl_blas_daxpy (-A, dx0, p);
 319       gsl_blas_daxpy (-B, dg0, p);
 320
 321       state->pnorm = gsl_blas_dnrm2 (p);
 322     }
 323
 324   gsl_vector_memcpy (g0, gradient);
 325   gsl_vector_memcpy (x0, x);
 326   state->g0norm = gsl_blas_dnrm2 (g0);
 327
 328 #ifdef DEBUG
 329   printf ("updated directions\n");
 330   printf ("p: ");
 331   gsl_vector_fprintf (stdout, p, "%g");
 332   printf ("g: ");
 333   gsl_vector_fprintf (stdout, gradient, "%g");
 334 #endif
 335
 336   return GSL_SUCCESS;
 337 }
 338
 339 static const gsl_multimin_fdfminimizer_type vector_bfgs_type = {
 340   "vector_bfgs",                /* name */
 341   sizeof (vector_bfgs_state_t),
 342   &vector_bfgs_alloc,
 343   &vector_bfgs_set,
 344   &vector_bfgs_iterate,
 345   &vector_bfgs_restart,
 346   &vector_bfgs_free
 347 };
 348
 349 const gsl_multimin_fdfminimizer_type
 350   * gsl_multimin_fdfminimizer_vector_bfgs = &vector_bfgs_type;