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 "gsl__config.h" 25 #include "gsl_multimin.h" 26 #include "gsl_blas.h" 27 28 #include "gsl_multimin__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; 351