36 namespace ublas = boost::numeric::ublas;
60 template <
typename VECTOR>
63 typedef VECTOR vector_type;
64 typedef typename vector_type::value_type T;
65 typedef typename vector_type::value_type value_type;
67 typedef value_type magnitude_type;
75 template<
typename VEC1,
typename VEC2>
76 void hmult(
const VEC1 &X, VEC2 &Y )
80 for (
size_type k = 0 ; k < delta.size(); ++k )
82 T xdelta = ublas::inner_prod( X, delta[k] );
83 T xzeta = ublas::inner_prod( X, zeta[k] );
88 Y.plus_assign( rho[k]*xdelta*zeta[k] );
89 Y.plus_assign( rho[k]*( xzeta-rho[k]*tau[k]*xdelta )*delta[k] );
93 Y.plus_assign( rho[k]*xdelta*delta[k] );
94 Y.minus_assign( xzeta/tau[k]*zeta[k] );
109 template<
typename VECT1,
typename VECT2>
110 void update(
const VECT1 &deltak,
const VECT2 &gammak )
125 delta[k].resize( N );
126 gamma[k].resize( N );
129 delta[k].assign( deltak );
130 gamma[k].assign( gammak );
132 rho[k] = magnitude_type( 1 ) / ublas::inner_prod( deltak, gammak );
134 if ( version == BFGS )
135 zeta[k].plus_assign( delta[k] - Y );
140 tau[k] = ublas::inner_prod( gammak, zeta[k] );
146 std::vector<vector_type> delta, gamma, zeta;
147 std::vector<T> tau, rho;
152 template <
typename FUNCTION,
typename DERIVATIVE,
typename VECTOR,
typename IterationBFGS>
153 void bfgs( FUNCTION f,
159 float lambda_init = 0.001,
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;
175 VECTOR r( x.size() ), d( x.size() ), y( x.size() ), r2( x.size() );
177 real_type lambda = lambda_init, valx = f( x ), valy;
181 while ( ! iter.isFinished( r ) )
183 invhessian.hmult( r, d );
187 real_type derivative = ublas::inner_prod( r, d );
188 real_type lambda_min( 0 );
189 real_type lambda_max( 0 );
192 bool unbounded =
true, blocked =
false, grad_computed =
false;
203 if ( iter.get_noisy() >= 2 )
205 cout <<
"Wolfe line search, lambda = " << lambda
206 <<
" value = " << valy /print_norm << endl;
211 if ( valy <= valx + m1 * lambda * derivative )
214 grad_computed =
true;
216 T derivative2 = ublas::inner_prod( r2, d );
218 if ( derivative2 >= m2*derivative )
231 lambda *= real_type( 10 );
234 lambda = ( lambda_max + lambda_min ) / real_type( 2 );
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 ) ) ) )
241 lambda = lambda_init;
249 if ( !grad_computed )
253 r.minus_assign( r2 );
255 if ( iter.numberOfIterations() % restart == 0 || blocked )
258 invhessian.restart();
260 if ( ++nb_restart > 10 )
270 invhessian.update( lambda*d, -r );
283 template <
typename FUNCTION,
typename DERIVATIVE,
typename VECTOR,
typename IterationBFGS>
292 bfgs( f, grad, x, restart, iter, version );