mmin_constr_spg.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*------------------------------------------------------------------------*
24  * Open Optimization Library - Spectral Projected Gradient Method
25  *
26  * spg/spg.c
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 2 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
41  *
42  * Ricardo Biloti
43  * since: May 3rd, 2005
44  *
45  * $Id: spg.c,v 1.4 2005/05/10 20:24:27 biloti Exp $
46  *------------------------------------------------------------------------*/
47 #ifndef O2SCL_OOL_MMIN_SPG_H
48 #define O2SCL_OOL_MMIN_SPG_H
49 
50 /** \file mmin_constr_spg.h
51  \brief File defining \ref o2scl::mmin_constr_spg
52 */
53 
54 #include <o2scl/multi_funct.h>
55 #include <o2scl/mmin_constr.h>
56 #include <gsl/gsl_math.h>
57 
58 #ifndef DOXYGEN_NO_O2NS
59 namespace o2scl {
60 #endif
61 
62  /** \brief Constrained minimization by the spectral projected
63  gradient method (OOL)
64 
65  This class applies a non-monotone line search strategy
66  to the classical projected gradient method.
67 
68  As in \ref Birgin00, this class applies a nonmonotone Armijo
69  sufficient decrease condition for accepting trial points as an
70  improvement over the classical spectral projected gradient
71  method. This method may be competitive with large problems
72  because it has low memory requirements.
73 
74  Default template arguments
75  - \c func_t - \ref multi_funct11
76  - \c dfunc_t - \ref grad_funct11
77  - \c vec_t - \ref boost::numeric::ublas::vector < double >
78 
79  \future There is some memory allocation which isn't
80  deallocated until the destructor, and should be handled
81  a bit more elegantly.
82  */
83  template<class func_t=multi_funct11, class dfunc_t=grad_funct11,
85  public mmin_constr<func_t,dfunc_t,ool_hfunct11,vec_t> {
86 
87 #ifndef DOXYGEN_INTERNAL
88 
89  protected:
90 
91  /// A convenient typedef for the unused Hessian product type
93 
94  /// Armijo parameter
95  double alpha;
96 
97  /// Temporary vector
98  vec_t xx;
99 
100  /// Temporary vector
101  vec_t d;
102 
103  /// Temporary vector
104  vec_t s;
105 
106  /// Temporary vector
107  vec_t y;
108 
109  /// Temporary vector
110  vec_t fvec;
111 
112  /// Non-monotone parameter
113  size_t m;
114 
115  /// Desc
116  int tail;
117 
118  /// Line search
119  int line_search() {
120 
121  double lambda, lambdanew;
122  double fmax, faux, fxx, dTg;
123  short armijo_failure;
124  size_t ii;
125 
126  // Saving the previous gradient and compute
127  // d = P(x - alpha * g) - x
128  for(size_t i=0;i<this->dim;i++) {
129  y[i]=this->gradient[i];
130  d[i]=this->x[i];
131  }
132  for(size_t i=0;i<this->dim;i++) {
133  d[i]-=alpha*this->gradient[i];
134  }
135  proj(d);
136  for(size_t i=0;i<this->dim;i++) {
137  d[i]-=this->x[i];
138  }
139 
140  dTg=0.0;
141  for(size_t i=0;i<this->dim;i++) dTg+=d[i]*this->gradient[i];
142 
143  lambda=1;
144  armijo_failure=1;
145  while (armijo_failure) {
146 
147  /* x trial */
148  for(size_t i=0;i<this->dim;i++) xx[i]=this->x[i];
149  for(size_t i=0;i<this->dim;i++) xx[i]+=lambda*d[i];
150  fxx=(*this->func)(this->dim,xx);
151 
152  fmax=0;
153  for(ii=0;ii<m;ii++) {
154  faux=fvec[ii]+gamma*lambda*dTg;
155  fmax=GSL_MAX(fmax,faux);
156  }
157 
158  if (fxx>fmax) {
159  armijo_failure=1;
160  lambdanew=-lambda*lambda*dTg/(2*(fxx-this->f-lambda*dTg));
161  double dtmp;
162  if (sigma2*lambda<lambdanew) dtmp=sigma2*lambda;
163  else dtmp=lambdanew;
164  if (sigma1*lambda>dtmp) lambda=sigma1*lambda;
165  else lambda=dtmp;
166  } else {
167  armijo_failure=0;
168  }
169  }
170 
171  // st->s = x_{k+1} - x_k
172  for(size_t i=0;i<this->dim;i++) s[i]=xx[i];
173  for(size_t i=0;i<this->dim;i++) s[i]-=this->x[i];
174 
175  // M->x = x_{k+1}
176  for(size_t i=0;i<this->dim;i++) this->x[i]=xx[i];
177  (*this->dfunc)(this->dim,this->x,this->gradient);
178  this->f=fxx;
179 
180  // Infinite norm of g1 = d <- [P(x - g) - x]
181  for(size_t i=0;i<this->dim;i++) d[i]=this->x[i];
182  for(size_t i=0;i<this->dim;i++) d[i]-=this->gradient[i];
183  proj(d);
184  for(size_t i=0;i<this->dim;i++) d[i]-=this->x[i];
185 
186  double maxval=fabs(d[0]);
187  for(size_t i=1;i<this->dim;i++) {
188  if (fabs(d[i])>maxval) {
189  maxval=fabs(d[i]);
190  }
191  }
192  this->size=fabs(maxval);
193 
194  // st->y = - (g(x_{k+1}) - g(x_k))
195  for(size_t i=0;i<this->dim;i++) y[i]-=this->gradient[i];
196 
197  m++;
198  if (M<m) m=M;
199  tail++;
200  tail=tail % M;
201  fvec[tail]=this->f;
202 
203  return 0;
204  }
205 
206  /// Project into feasible region
207  int proj(vec_t &xt) {
208  size_t ii, n=this->dim;
209 
210  for(ii=0;ii<n;ii++) {
211  double dtmp;
212  if (xt[ii]<this->U[ii]) dtmp=xt[ii];
213  else dtmp=this->U[ii];
214  if (this->L[ii]>dtmp) xt[ii]=this->L[ii];
215  else xt[ii]=dtmp;
216  //xt[ii]=GSL_MAX(this->L[ii],GSL_MIN(xt[ii],this->U[ii]));
217  }
218  return 0;
219  }
220 
221 #endif
222 
223  public:
224 
225  mmin_constr_spg() {
226  fmin=-1.0e99;
227  tol=1.0e-4;
228  alphamin=1.0e-30;
229  alphamax=1.0e30;
230  gamma=1.0e-4;
231  sigma1=0.1;
232  sigma2=0.9;
233  M=10;
234  }
235 
236  ~mmin_constr_spg() {
237  fvec.clear();
238  }
239 
240  /** \brief Minimum function value (default \f$ 10^{-99} \f$)
241 
242  If the function value is below this value, then the algorithm
243  assumes that the function is not bounded and exits.
244  */
245  double fmin;
246 
247  /// Tolerance on infinite norm (default \f$ 10^{-4} \f$)
248  double tol;
249 
250  /// Lower bound to spectral step size (default \f$ 10^{-30} \f$)
251  double alphamin;
252 
253  /// Upper bound to spectral step size (default \f$ 10^{30} \f$)
254  double alphamax;
255 
256  /// Sufficient decrease parameter (default \f$ 10^{-4} \f$)
257  double gamma;
258 
259  /// Lower bound to the step length reduction (default 0.1)
260  double sigma1;
261 
262  /// Upper bound to the step length reduction (default 0.9)
263  double sigma2;
264 
265  /// Monotonicity parameter (M=1 forces monotonicity) (default 10)
266  size_t M;
267 
268  /// Allocate memory
269  virtual int allocate(const size_t n) {
270  if (this->dim>0) free();
271  xx.resize(n);
272  d.resize(n);
273  s.resize(n);
274  y.resize(n);
276  }
277 
278  /// Free previously allocated memory
279  virtual int free() {
280  xx.clear();
281  d.clear();
282  s.clear();
283  y.clear();
284  return 0;
285  }
286 
287  /// Set the function, the initial guess, and the parameters
288  virtual int set(func_t &fn, dfunc_t &dfn, vec_t &init) {
289 
290  int ret=mmin_constr<func_t,dfunc_t,hfunc_t,
291  vec_t>::set(fn,dfn,init);
292 
293  // Turn initial guess into a feasible point
294  proj(this->x);
295 
296  // Evaluate function and gradient
297  this->f=(*this->func)(this->dim,this->x);
298  (*this->dfunc)(this->dim,this->x,this->gradient);
299 
300  /* Infinite norm of d <- g1 = [P(x - g) - x] */
301  o2scl::vector_copy(this->dim,this->x,d);
302  for(size_t i=0;i<this->dim;i++) {
303  d[i]-=this->gradient[i];
304  }
305  proj(d);
306  for(size_t i=0;i<this->dim;i++) d[i]-=this->x[i];
307 
308  double maxval=fabs(d[0]);
309  for(size_t i=1;i<this->dim;i++) {
310  if (fabs(d[i])>maxval) {
311  maxval=fabs(d[i]);
312  }
313  }
314  this->size=fabs(maxval);
315 
316  /* alpha_0 */
317  maxval=fabs(this->gradient[0]);
318  for(size_t i=1;i<this->dim;i++) {
319  if (fabs(this->gradient[i])>maxval) {
320  maxval=fabs(this->gradient[i]);
321  }
322  }
323  alpha=1.0/fabs(maxval);
324 
325  // FIXME: at the moment, this memory allocation is
326  // left hanging. It's taken care of the destructor,
327  // but this isn't so elegant.
328  fvec.resize(M);
329  m=1;
330  tail=0;
331  fvec[tail]=this->f;
332 
333  return 0;
334  }
335 
336  /// Restart the minimizer
337  virtual int restart() {
338 
339  proj(this->x);
340 
341  // Evaluate function and gradient
342  this->f=(*this->func)(this->dim,this->x);
343  (*this->dfunc)(this->dim,this->x,this->gradient);
344 
345  // alpha_0
346  double maxval=fabs(this->gradient[0]);
347  for(size_t i=1;i<this->dim;i++) {
348  if (fabs(this->gradient[i])>maxval) {
349  maxval=fabs(this->gradient[i]);
350  }
351  }
352  alpha=1.0/fabs(maxval);
353 
354  // Infinite norm of d <- g1 = [P(x - g) - x]
355  o2scl::vector_copy(this->dim,this->x,d);
356  for(size_t i=0;i<this->dim;i++) {
357  d[i]-=this->gradient[i];
358  }
359  proj(d);
360  for(size_t i=0;i<this->dim;i++) d[i]-=this->x[i];
361 
362  maxval=fabs(d[0]);
363  for(size_t i=1;i<this->dim;i++) {
364  if (fabs(d[i])>maxval) {
365  maxval=fabs(d[i]);
366  }
367  }
368  this->size=fabs(maxval);
369 
370  m=1;
371  tail=0;
372  fvec[tail]=this->f;
373 
374  return 0;
375  }
376 
377  /// Perform an iteration
378  virtual int iterate() {
379  double t;
380 
381  line_search();
382 
383  double bk=0.0;
384  for(size_t i=0;i<this->dim;i++) bk+=s[i]*y[i];
385 
386  if (bk>=0.0) {
387  alpha=alphamax;
388  } else {
389  double ak=0.0;
390  for(size_t i=0;i<this->dim;i++) ak+=s[i]*s[i];
391  ak=-ak/bk;
392  double dtmp;
393  if (alphamin>ak) dtmp=alphamin;
394  else dtmp=ak;
395  if (alphamax<dtmp) alpha=alphamax;
396  else alpha=dtmp;
397  //alpha=GSL_MIN(alphamax,GSL_MAX(alphamin,ak));
398  }
399 
400  return 0;
401  }
402 
403  /// See if we're finished
404  virtual int is_optimal() {
405  if (this->size>tol && this->f>fmin) {
406  return GSL_CONTINUE;
407  }
408  return GSL_SUCCESS;
409  }
410 
411  /// Return string denoting type ("mmin_constr_spg")
412  const char *type() { return "mmin_constr_spg"; }
413 
414 #ifndef DOXYGEN_INTERNAL
415 
416  private:
417 
421  (const mmin_constr_spg<func_t,dfunc_t,vec_t>&);
422 
423 #endif
424 
425  };
426 
427 #ifndef DOXYGEN_NO_O2NS
428 }
429 #endif
430 
431 #endif
432 
Class for automatically computing gradients [abstract base].
Definition: mmin.h:54
virtual int free()
Free previously allocated memory.
vec_t d
Temporary vector.
vec_t s
Temporary vector.
vec_t xx
Temporary vector.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
vec_t y
Temporary vector.
virtual int iterate()
Perform an iteration.
const char * type()
Return string denoting type ("mmin_constr_spg")
size_t m
Non-monotone parameter.
ool_hfunct11 hfunc_t
A convenient typedef for the unused Hessian product type.
double tol
Tolerance on infinite norm (default )
int line_search()
Line search.
double gamma
Sufficient decrease parameter (default )
double sigma2
Upper bound to the step length reduction (default 0.9)
Constrained multidimensional minimization (OOL) [abstract base].
Definition: mmin_constr.h:78
virtual int restart()
Restart the minimizer.
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:127
virtual int is_optimal()
See if we&#39;re finished.
size_t M
Monotonicity parameter (M=1 forces monotonicity) (default 10)
double alphamin
Lower bound to spectral step size (default )
virtual int allocate(const size_t n)
Allocate memory.
virtual int set(func_t &fn, dfunc_t &dfn, vec_t &init)
Set the function, the initial guess, and the parameters.
std::function< int(size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> grad_funct11
Array of multi-dimensional functions typedef.
Definition: mmin.h:42
double fmin
Minimum function value (default )
double alphamax
Upper bound to spectral step size (default )
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ool_hfunct11
Hessian product function for o2scl::mmin_constr.
Definition: mmin_constr.h:66
virtual int allocate(const size_t n)
Allocate memory.
Definition: mmin_constr.h:257
double sigma1
Lower bound to the step length reduction (default 0.1)
int proj(vec_t &xt)
Project into feasible region.
Constrained minimization by the spectral projected gradient method (OOL)
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct11
Multi-dimensional function typedef.
Definition: multi_funct.h:45
vec_t fvec
Temporary vector.
double alpha
Armijo parameter.

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).