mmin_constr_spg.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_funct
76  - \c dfunc_t - \ref grad_funct
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_funct, class dfunc_t=grad_funct,
85  public mmin_constr<func_t,dfunc_t,ool_hfunct,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 
422 
423 #endif
424 
425  };
426 
427 #ifndef DOXYGEN_NO_O2NS
428 }
429 #endif
430 
431 #endif
432 
o2scl::mmin_constr
Constrained multidimensional minimization (OOL) [abstract base].
Definition: mmin_constr.h:78
boost::numeric::ublas::vector< double >
o2scl::gradient
Class for automatically computing gradients [abstract base].
Definition: mmin.h:54
o2scl::mmin_constr_spg::set
virtual int set(func_t &fn, dfunc_t &dfn, vec_t &init)
Set the function, the initial guess, and the parameters.
Definition: mmin_constr_spg.h:288
o2scl::mmin_constr_spg::fmin
double fmin
Minimum function value (default )
Definition: mmin_constr_spg.h:245
o2scl::mmin_constr_spg::gamma
double gamma
Sufficient decrease parameter (default )
Definition: mmin_constr_spg.h:257
o2scl::mmin_constr_spg::m
size_t m
Non-monotone parameter.
Definition: mmin_constr_spg.h:113
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::mmin_constr_spg::s
vec_t s
Temporary vector.
Definition: mmin_constr_spg.h:104
o2scl::mmin_constr_spg::sigma2
double sigma2
Upper bound to the step length reduction (default 0.9)
Definition: mmin_constr_spg.h:263
o2scl::vector_copy
void vector_copy(const vec_t &src, vec2_t &dest)
Simple vector copy.
Definition: vector.h:124
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::dim
size_t dim
Number of parameters.
Definition: mmin_constr.h:106
o2scl::mmin_constr_spg::tail
int tail
Desc.
Definition: mmin_constr_spg.h:116
o2scl::mmin_constr_spg::xx
vec_t xx
Temporary vector.
Definition: mmin_constr_spg.h:98
o2scl::mmin_constr_spg::free
virtual int free()
Free previously allocated memory.
Definition: mmin_constr_spg.h:279
o2scl::mmin_constr_spg::sigma1
double sigma1
Lower bound to the step length reduction (default 0.1)
Definition: mmin_constr_spg.h:260
o2scl::mmin_constr_spg::proj
int proj(vec_t &xt)
Project into feasible region.
Definition: mmin_constr_spg.h:207
o2scl::mmin_constr_spg::fvec
vec_t fvec
Temporary vector.
Definition: mmin_constr_spg.h:110
o2scl::ool_hfunct
std::function< int(size_t, const boost::numeric::ublas::vector< double > &, const boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> ool_hfunct
Hessian product function for o2scl::mmin_constr.
Definition: mmin_constr.h:66
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::size
double size
Desc.
Definition: mmin_constr.h:89
o2scl::multi_funct
std::function< double(size_t, const boost::numeric::ublas::vector< double > &)> multi_funct
Multi-dimensional function typedef in src/base/multi_funct.h.
Definition: multi_funct.h:45
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::x
boost::numeric::ublas::vector< double > x
The current minimum vector.
Definition: mmin_constr.h:92
o2scl::mmin_constr_spg::allocate
virtual int allocate(const size_t n)
Allocate memory.
Definition: mmin_constr_spg.h:269
o2scl::mmin_constr::allocate
virtual int allocate(const size_t n)
Allocate memory.
Definition: mmin_constr.h:257
o2scl::mmin_constr_spg::hfunc_t
ool_hfunct hfunc_t
A convenient typedef for the unused Hessian product type.
Definition: mmin_constr_spg.h:92
o2scl::mmin_constr_spg::y
vec_t y
Temporary vector.
Definition: mmin_constr_spg.h:107
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::func
multi_funct * func
User-supplied function.
Definition: mmin_constr.h:110
o2scl::mmin_constr_spg::tol
double tol
Tolerance on infinite norm (default )
Definition: mmin_constr_spg.h:248
o2scl::grad_funct
std::function< int(size_t, boost::numeric::ublas::vector< double > &, boost::numeric::ublas::vector< double > &)> grad_funct
Array of multi-dimensional functions typedef in src/min/mmin.h.
Definition: mmin.h:42
o2scl::mmin_constr_spg::alpha
double alpha
Armijo parameter.
Definition: mmin_constr_spg.h:95
o2scl::mmin_constr_spg::alphamax
double alphamax
Upper bound to spectral step size (default )
Definition: mmin_constr_spg.h:254
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::L
boost::numeric::ublas::vector< double > L
Lower bound constraints.
Definition: mmin_constr.h:117
o2scl::mmin_constr_spg::restart
virtual int restart()
Restart the minimizer.
Definition: mmin_constr_spg.h:337
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::U
boost::numeric::ublas::vector< double > U
Upper bound constraints.
Definition: mmin_constr.h:119
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::dfunc
grad_funct * dfunc
Gradient function.
Definition: mmin_constr.h:112
o2scl::mmin_constr_spg::iterate
virtual int iterate()
Perform an iteration.
Definition: mmin_constr_spg.h:378
o2scl::mmin_constr_spg::line_search
int line_search()
Line search.
Definition: mmin_constr_spg.h:119
o2scl::mmin_constr< multi_funct, grad_funct, ool_hfunct, boost::numeric::ublas::vector< double > >::f
double f
The current function value.
Definition: mmin_constr.h:87
o2scl::mmin_constr_spg::is_optimal
virtual int is_optimal()
See if we're finished.
Definition: mmin_constr_spg.h:404
o2scl::mmin_constr_spg
Constrained minimization by the spectral projected gradient method (OOL)
Definition: mmin_constr_spg.h:84
o2scl::mmin_constr_spg::d
vec_t d
Temporary vector.
Definition: mmin_constr_spg.h:101
o2scl::mmin_constr_spg::M
size_t M
Monotonicity parameter (M=1 forces monotonicity) (default 10)
Definition: mmin_constr_spg.h:266
o2scl::mmin_constr_spg::type
const char * type()
Return string denoting type ("mmin_constr_spg")
Definition: mmin_constr_spg.h:412
o2scl::mmin_constr_spg::alphamin
double alphamin
Lower bound to spectral step size (default )
Definition: mmin_constr_spg.h:251

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