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 "gsl__config.h"
29 #include "gsl_multimin.h"
30 #include "gsl_blas.h"
31
32 #include "gsl_multimin__linear_minimize.c"
33 #include "gsl_multimin__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
vector_bfgs2_alloc(void * vstate,size_t n)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
vector_bfgs2_set(void * vstate,gsl_multimin_function_fdf * fdf,const gsl_vector * x,double * f,gsl_vector * gradient,double step_size,double tol)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
vector_bfgs2_free(void * vstate)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
vector_bfgs2_restart(void * vstate)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
vector_bfgs2_iterate(void * vstate,gsl_multimin_function_fdf * fdf,gsl_vector * x,double * f,gsl_vector * gradient,gsl_vector * dx)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;
331