Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
bfgs.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-05-02
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
30 
31 #ifndef _BFGS_HPP_
32 #define _BFGS_HPP_
33 
34 namespace Feel
35 {
36 namespace ublas = boost::numeric::ublas;
44 enum BFGSType { BFGS = 0, DFP };
45 
59 // Object representing the inverse of the Hessian
60 template <typename VECTOR>
62 {
63  typedef VECTOR vector_type;
64  typedef typename vector_type::value_type T;
65  typedef typename vector_type::value_type value_type;
66  //typedef typename number_traits<T>::magnitude_type R;
67  typedef value_type magnitude_type;
68 
69 
70  BFGSInvHessian( BFGSType v = BFGS )
71  {
72  version = v;
73  }
74 
75  template<typename VEC1, typename VEC2>
76  void hmult( const VEC1 &X, VEC2 &Y )
77  {
78  Y.assign( X );
79 
80  for ( size_type k = 0 ; k < delta.size(); ++k )
81  {
82  T xdelta = ublas::inner_prod( X, delta[k] );
83  T xzeta = ublas::inner_prod( X, zeta[k] );
84 
85  switch ( version )
86  {
87  case BFGS :
88  Y.plus_assign( rho[k]*xdelta*zeta[k] );
89  Y.plus_assign( rho[k]*( xzeta-rho[k]*tau[k]*xdelta )*delta[k] );
90  break;
91 
92  case DFP :
93  Y.plus_assign( rho[k]*xdelta*delta[k] );
94  Y.minus_assign( xzeta/tau[k]*zeta[k] );
95  break;
96  }
97  }
98  }
99 
100  void restart( void )
101  {
102  delta.clear();
103  gamma.clear();
104  zeta.clear();
105  tau.clear();
106  rho.clear();
107  }
108 
109  template<typename VECT1, typename VECT2>
110  void update( const VECT1 &deltak, const VECT2 &gammak )
111  {
112  size_type N = deltak.size();
113  size_type k = delta.size();
114 
115  vector_type Y( N );
116 
117  hmult( gammak, Y );
118 
119  delta.resize( k+1 );
120  gamma.resize( k+1 );
121  zeta.resize( k+1 );
122  tau.resize( k+1 );
123  rho.resize( k+1 );
124 
125  delta[k].resize( N );
126  gamma[k].resize( N );
127  zeta[k].resize( N );
128 
129  delta[k].assign( deltak );
130  gamma[k].assign( gammak );
131 
132  rho[k] = magnitude_type( 1 ) / ublas::inner_prod( deltak, gammak );
133 
134  if ( version == BFGS )
135  zeta[k].plus_assign( delta[k] - Y );
136 
137  else
138  zeta[k].assign( Y );
139 
140  tau[k] = ublas::inner_prod( gammak, zeta[k] );
141  }
142 
143  //
144  // Data
145  //
146  std::vector<vector_type> delta, gamma, zeta;
147  std::vector<T> tau, rho;
148  int version;
149 };
150 
151 
152 template <typename FUNCTION, typename DERIVATIVE, typename VECTOR, typename IterationBFGS>
153 void bfgs( FUNCTION f,
154  DERIVATIVE grad,
155  VECTOR &x,
156  int restart,
157  IterationBFGS& iter,
158  BFGSType version = BFGS,
159  float lambda_init = 0.001,
160  float /*print_norm*/ = 1.0 )
161 {
162 
163  //typedef typename linalg_traits<VECTOR>::value_type T;
164  //typedef typename number_traits<T>::magnitude_type R;
165 
166  typedef VECTOR vector_type;
167  typedef typename vector_type::value_type T;
168  typedef typename vector_type::value_type value_type;
169  typedef typename type_traits<T>::real_type real_type;
170  typedef value_type magnitude_type;
171 
172 
173  BFGSInvHessian<vector_type> invhessian( version );
174 
175  VECTOR r( x.size() ), d( x.size() ), y( x.size() ), r2( x.size() );
176  grad( x, r );
177  real_type lambda = lambda_init, valx = f( x ), valy;
178  int nb_restart( 0 );
179 
180  //if (iter.get_noisy() >= 1) cout << "value " << valx / print_norm << " ";
181  while ( ! iter.isFinished( r ) )
182  {
183  invhessian.hmult( r, d );
184  d *= -1;
185 
186  // Wolfe Line search
187  real_type derivative = ublas::inner_prod( r, d );
188  real_type lambda_min( 0 );
189  real_type lambda_max( 0 );
190  real_type m1 = 0.27;
191  real_type m2 = 0.57;
192  bool unbounded = true, blocked = false, grad_computed = false;
193 
194  for ( ;; )
195  {
196 
197  //add(x, scaled(d, lambda), y);
198  y = lambda*d+x;
199  valy = f( y );
200 
201 #if 0
202 
203  if ( iter.get_noisy() >= 2 )
204  {
205  cout << "Wolfe line search, lambda = " << lambda
206  << " value = " << valy /print_norm << endl;
207  }
208 
209 #endif
210 
211  if ( valy <= valx + m1 * lambda * derivative )
212  {
213  grad( y, r2 );
214  grad_computed = true;
215 
216  T derivative2 = ublas::inner_prod( r2, d );
217 
218  if ( derivative2 >= m2*derivative )
219  break;
220 
221  lambda_min = lambda;
222  }
223 
224  else
225  {
226  lambda_max = lambda;
227  unbounded = false;
228  }
229 
230  if ( unbounded )
231  lambda *= real_type( 10 );
232 
233  else
234  lambda = ( lambda_max + lambda_min ) / real_type( 2 );
235 
236  if ( valy <= real_type( 2 )*valx &&
237  ( lambda < real_type( lambda_init*1E-8 ) ||
238  ( !unbounded && lambda_max-lambda_min < real_type( lambda_init*1E-8 ) ) ) )
239  {
240  blocked = true;
241  lambda = lambda_init;
242  break;
243  }
244  }
245 
246  // Rank two update
247  ++iter;
248 
249  if ( !grad_computed )
250  grad( y, r2 );
251 
252  //gmm::add(scaled(r2, -1), r);
253  r.minus_assign( r2 );
254 
255  if ( iter.numberOfIterations() % restart == 0 || blocked )
256  {
257  //if (iter.get_noisy() >= 1) cout << "Restart\n";
258  invhessian.restart();
259 
260  if ( ++nb_restart > 10 )
261  {
262  //if (iter.get_noisy() >= 1) cout << "BFGS is blocked, exiting\n";
263  return;
264  }
265  }
266 
267  else
268  {
269  //invhessian.update(gmm::scaled(d,lambda), gmm::scaled(r,-1));
270  invhessian.update( lambda*d, -r );
271  nb_restart = 0;
272  }
273 
274  r.assign( r2 );
275  x.assign( y );
276  valx = valy;
277  //if (iter.get_noisy() >= 1) cout << "BFGS value " << valx/print_norm << "\t";
278  }
279 
280 }
281 
282 
283 template <typename FUNCTION, typename DERIVATIVE, typename VECTOR, typename IterationBFGS>
284 inline void
285 dfp( FUNCTION f,
286  DERIVATIVE grad,
287  VECTOR &x,
288  int restart,
289  IterationBFGS& iter,
290  BFGSType version = DFP )
291 {
292  bfgs( f, grad, x, restart, iter, version );
293 }
294 
295 
296 }
297 #endif
Definition: solverlinear.hpp:33
BFGSType
Definition: bfgs.hpp:44
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Definition: bfgs.hpp:61

Generated on Sun Dec 22 2013 13:10:55 for Feel++ by doxygen 1.8.5