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

   1 /* multimin/vector_bfgs2.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   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
  18  * 02110-1301, USA.
  19  */
  20
  21 /* vector_bfgs2.c -- Fletcher's implementation of the BFGS method,
  22    from R.Fletcher, "Practical Method's of Optimization", Second
  23    Edition, ISBN 0471915475.  Algorithms 2.6.2 and 2.6.4. */
  24
  25 /* Thanks to Alan Irwin irwin@beluga.phys.uvic.ca. for suggesting this
  26    algorithm and providing sample fortran benchmarks */
  27
  28 #include <config.h>
  29 #include <gsl/gsl_multimin.h>
  30 #include <gsl/gsl_blas.h>
  31
  32 #include "linear_minimize.c"
  33 #include "linear_wrapper.c"
  34
  35 typedef struct
  36 {
  37   int iter;
  38   double step;
  39   double g0norm;
  40   double pnorm;
  41   double delta_f;
  42   double fp0;                   /* f'(0) for f(x-alpha*p) */
  43   gsl_vector *x0;
  44   gsl_vector *g0;
  45   gsl_vector *p;
  46   /* work space */
  47   gsl_vector *dx0;
  48   gsl_vector *dg0;
  49   gsl_vector *x_alpha;
  50   gsl_vector *g_alpha;
  51   /* wrapper function */
  52   wrapper_t wrap;
  53   /* minimization parameters */
  54   double rho;
  55   double sigma;
  56   double tau1;
  57   double tau2;
  58   double tau3;
  59   int order;
  60 }
  61 vector_bfgs2_state_t;
  62
  63 static int
  64 vector_bfgs2_alloc (void *vstate, size_t n)
  65 {
  66   vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
  67
  68   state->p = gsl_vector_calloc (n);
  69
  70   if (state->p == 0)
  71     {
  72       GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
  73     }
  74
  75   state->x0 = gsl_vector_calloc (n);
  76
  77   if (state->x0 == 0)
  78     {
  79       gsl_vector_free (state->p);
  80       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
  81     }
  82
  83   state->g0 = gsl_vector_calloc (n);
  84
  85   if (state->g0 == 0)
  86     {
  87       gsl_vector_free (state->x0);
  88       gsl_vector_free (state->p);
  89       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
  90     }
  91
  92   state->dx0 = gsl_vector_calloc (n);
  93
  94   if (state->dx0 == 0)
  95     {
  96       gsl_vector_free (state->g0);
  97       gsl_vector_free (state->x0);
  98       gsl_vector_free (state->p);
  99       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 100     }
 101
 102   state->dg0 = gsl_vector_calloc (n);
 103
 104   if (state->dg0 == 0)
 105     {
 106       gsl_vector_free (state->dx0);
 107       gsl_vector_free (state->g0);
 108       gsl_vector_free (state->x0);
 109       gsl_vector_free (state->p);
 110       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 111     }
 112
 113   state->x_alpha = gsl_vector_calloc (n);
 114
 115   if (state->x_alpha == 0)
 116     {
 117       gsl_vector_free (state->dg0);
 118       gsl_vector_free (state->dx0);
 119       gsl_vector_free (state->g0);
 120       gsl_vector_free (state->x0);
 121       gsl_vector_free (state->p);
 122       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 123     }
 124
 125   state->g_alpha = gsl_vector_calloc (n);
 126
 127   if (state->g_alpha == 0)
 128     {
 129       gsl_vector_free (state->x_alpha);
 130       gsl_vector_free (state->dg0);
 131       gsl_vector_free (state->dx0);
 132       gsl_vector_free (state->g0);
 133       gsl_vector_free (state->x0);
 134       gsl_vector_free (state->p);
 135       GSL_ERROR ("failed to allocate space for g0", GSL_ENOMEM);
 136     }
 137
 138   return GSL_SUCCESS;
 139 }
 140
 141 static int
 142 vector_bfgs2_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_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
 147
 148   state->iter = 0;
 149   state->step = step_size;
 150   state->delta_f = 0;
 151
 152   GSL_MULTIMIN_FN_EVAL_F_DF (fdf, x, f, gradient);
 153
 154   /* Use the gradient as the initial direction */
 155
 156   gsl_vector_memcpy (state->x0, x);
 157   gsl_vector_memcpy (state->g0, gradient);
 158   state->g0norm = gsl_blas_dnrm2 (state->g0);
 159
 160   gsl_vector_memcpy (state->p, gradient);
 161   gsl_blas_dscal (-1 / state->g0norm, state->p);
 162   state->pnorm = gsl_blas_dnrm2 (state->p);     /* should be 1 */
 163   state->fp0 = -state->g0norm;
 164
 165   /* Prepare the wrapper */
 166
 167   prepare_wrapper (&state->wrap, fdf,
 168                    state->x0, *f, state->g0,
 169                    state->p, state->x_alpha, state->g_alpha);
 170
 171   /* Prepare 1d minimisation parameters */
 172
 173   state->rho = 0.01;
 174   state->sigma = tol;
 175   state->tau1 = 9;
 176   state->tau2 = 0.05;
 177   state->tau3 = 0.5;
 178   state->order = 3;  /* use cubic interpolation where possible */
 179
 180   return GSL_SUCCESS;
 181 }
 182
 183 static void
 184 vector_bfgs2_free (void *vstate)
 185 {
 186   vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
 187
 188   gsl_vector_free (state->x_alpha);
 189   gsl_vector_free (state->g_alpha);
 190   gsl_vector_free (state->dg0);
 191   gsl_vector_free (state->dx0);
 192   gsl_vector_free (state->g0);
 193   gsl_vector_free (state->x0);
 194   gsl_vector_free (state->p);
 195 }
 196
 197 static int
 198 vector_bfgs2_restart (void *vstate)
 199 {
 200   vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
 201
 202   state->iter = 0;
 203   return GSL_SUCCESS;
 204 }
 205
 206 static int
 207 vector_bfgs2_iterate (void *vstate, gsl_multimin_function_fdf * fdf,
 208                       gsl_vector * x, double *f,
 209                       gsl_vector * gradient, gsl_vector * dx)
 210 {
 211   vector_bfgs2_state_t *state = (vector_bfgs2_state_t *) vstate;
 212   double alpha = 0.0, alpha1;
 213   gsl_vector *x0 = state->x0;
 214   gsl_vector *g0 = state->g0;
 215   gsl_vector *p = state->p;
 216
 217   double g0norm = state->g0norm;
 218   double pnorm = state->pnorm;
 219   double delta_f = state->delta_f;
 220   double pg, dir;
 221   int status;
 222
 223   double f0 = *f;
 224
 225   if (pnorm == 0.0 || g0norm == 0.0 || state->fp0 == 0)
 226     {
 227       gsl_vector_set_zero (dx);
 228       return GSL_ENOPROG;
 229     }
 230
 231   if (delta_f < 0)
 232     {
 233       double del = GSL_MAX_DBL (-delta_f, 10 * GSL_DBL_EPSILON * fabs(f0));
 234       alpha1 = GSL_MIN_DBL (1.0, 2.0 * del / (-state->fp0));
 235     }
 236   else
 237     {
 238       alpha1 = fabs(state->step);
 239     }
 240
 241   /* line minimisation, with cubic interpolation (order = 3) */
 242
 243   status = minimize (&state->wrap.fdf_linear, state->rho, state->sigma,
 244                      state->tau1, state->tau2, state->tau3, state->order,
 245                      alpha1,  &alpha);
 246
 247   if (status != GSL_SUCCESS)
 248     {
 249       return status;
 250     }
 251
 252   update_position (&(state->wrap), alpha, x, f, gradient);
 253
 254   state->delta_f = *f - f0;
 255
 256   /* Choose a new direction for the next step */
 257
 258   {
 259     /* This is the BFGS update: */
 260     /* p' = g1 - A dx - B dg */
 261     /* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */
 262     /* B = dx.g/dx.dg */
 263
 264     gsl_vector *dx0 = state->dx0;
 265     gsl_vector *dg0 = state->dg0;
 266
 267     double dxg, dgg, dxdg, dgnorm, A, B;
 268
 269     /* dx0 = x - x0 */
 270     gsl_vector_memcpy (dx0, x);
 271     gsl_blas_daxpy (-1.0, x0, dx0);
 272
 273     gsl_vector_memcpy (dx, dx0);  /* keep a copy */
 274
 275     /* dg0 = g - g0 */
 276     gsl_vector_memcpy (dg0, gradient);
 277     gsl_blas_daxpy (-1.0, g0, dg0);
 278
 279     gsl_blas_ddot (dx0, gradient, &dxg);
 280     gsl_blas_ddot (dg0, gradient, &dgg);
 281     gsl_blas_ddot (dx0, dg0, &dxdg);
 282
 283     dgnorm = gsl_blas_dnrm2 (dg0);
 284
 285     if (dxdg != 0)
 286       {
 287         B = dxg / dxdg;
 288         A = -(1.0 + dgnorm * dgnorm / dxdg) * B + dgg / dxdg;
 289       }
 290     else
 291       {
 292         B = 0;
 293         A = 0;
 294       }
 295
 296     gsl_vector_memcpy (p, gradient);
 297     gsl_blas_daxpy (-A, dx0, p);
 298     gsl_blas_daxpy (-B, dg0, p);
 299   }
 300
 301   gsl_vector_memcpy (g0, gradient);
 302   gsl_vector_memcpy (x0, x);
 303   state->g0norm = gsl_blas_dnrm2 (g0);
 304   state->pnorm = gsl_blas_dnrm2 (p);
 305
 306   /* update direction and fp0 */
 307
 308   gsl_blas_ddot (p, gradient, &pg);
 309   dir = (pg >= 0.0) ? -1.0 : +1.0;
 310   gsl_blas_dscal (dir / state->pnorm, p);
 311   state->pnorm = gsl_blas_dnrm2 (p);
 312   gsl_blas_ddot (p, g0, &state->fp0);
 313
 314   change_direction (&state->wrap);
 315
 316   return GSL_SUCCESS;
 317 }
 318
 319 static const gsl_multimin_fdfminimizer_type vector_bfgs2_type = {
 320   "vector_bfgs2",               /* name */
 321   sizeof (vector_bfgs2_state_t),
 322   &vector_bfgs2_alloc,
 323   &vector_bfgs2_set,
 324   &vector_bfgs2_iterate,
 325   &vector_bfgs2_restart,
 326   &vector_bfgs2_free
 327 };
 328
 329 const gsl_multimin_fdfminimizer_type
 330   * gsl_multimin_fdfminimizer_vector_bfgs2 = &vector_bfgs2_type;