29 #ifndef SOLVERUNCONSTRAINED_H
30 #define SOLVERUNCONSTRAINED_H
32 #include <boost/numeric/ublas/banded.hpp>
33 #include <boost/numeric/ublas/symmetric.hpp>
34 #include <boost/numeric/ublas/io.hpp>
35 #include <boost/lambda/lambda.hpp>
36 #include <boost/lambda/if.hpp>
40 #include <feel/feelcore/feel.hpp>
53 template<
class>
class Problem = problem
90 typedef typename problem_type::f_type
f_type;
100 typedef typename problem_type::symmetric_matrix_type symmetric_matrix_type;
101 typedef ublas::banded_matrix<double> banded_matrix_type;
103 typedef directionalScalingMatrix<value_type> dsm_type;
118 allow_Trust_radius_calculations(
true ),
121 rho_increase_big( 1.5 ),
123 rho_increase_small( 1.1 )
133 bool allow_Trust_radius_calculations;
180 int n_CGiter,
int n_restarts,
int n_indef,
181 int n_crosses_def,
int n_crosses_indef,
182 int n_truss_exit_def,
int n_truss_exit_indef,
196 std::vector<value_type> norm_Tgrad_fx_hstr, norm_err;
197 std::vector<int> n_CGiter_hstr, n_restarts_hstr, n_indef_hstr,
198 n_crosses_def_hstr, n_crosses_indef_hstr,
199 n_truss_exit_def_hstr, n_truss_exit_indef_hstr;
200 std::vector<value_type> Delta_hstr, ared_til_hstr, phi_til_hstr, rho_hstr;
221 M_prob.define_problem( x_definitions );
229 return M_solver_stats;
250 banded_matrix_type & _Hg,
251 symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil,
258 banded_matrix_type & _Hg,
259 symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil );
268 int &_n_restarts,
int &_n_indef,
269 int &_n_crosses_def,
int &_n_crosses_indef,
270 int &_n_truss_exit_def,
int &_n_truss_exit_indef,
317 template<
typename Data,
template<
class>
class Problem>
320 M_prob( x_definitions ),
321 M_solver_stats ( this, true ),
326 template<
typename Data,
template<
class>
class Problem>
330 M_solver_stats ( this, true ),
331 M_theta( M_prob.lowerBounds(), M_prob.upperBounds() )
336 template<
typename Data,
template<
class>
class Problem>
341 template<
typename Data,
template<
class>
class Problem>
345 M_solver_stats.clear();
353 value_type norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
356 M_prob.copy_x0_to_x( x );
362 if ( M_solver_stats.collectStats() )
363 M_solver_stats.push( norm_Tgrad_fx, 0, 0, 0, 0, 0, 0, 0, Delta, 0, 0, 0 );
374 ( iter < M_options.max_TR_iter && norm_Tgrad_fx > M_options.TR_tol ) )
377 int n_CGiter, n_restarts, n_indef,
378 n_crosses_def, n_crosses_indef, n_truss_exit_def, n_truss_exit_indef;
379 value_type _s_til_x_G_til_x_s_til, phi_til, rho, Delta_used = Delta;
381 DVLOG(2) <<
"\n===================== iter = " << iter <<
" ===========================";
383 DVLOG(2) <<
"\n -> norm_Tgrad_fx = " << norm_Tgrad_fx;
391 this->CGstep( x, Delta, stot, norm_s_til,
394 n_crosses_def, n_crosses_indef,
395 n_truss_exit_def, n_truss_exit_indef,
396 _s_til_x_G_til_x_s_til, phi_til );
402 M_prob.evaluate( x_new, __fx_new, diff_order<0>() );
404 M_prob.evaluate( x, __fx, diff_order<0>() );
407 value_type ared_til = __fx_new.value( 0 ) - __fx.value( 0 ) + 0.5 * _s_til_x_G_til_x_s_til;
409 rho = ared_til / phi_til;
413 if ( M_options.allow_Trust_radius_calculations )
416 Delta = M_options.rho_decrease * Delta;
422 if ( rho > M_options.rho_big )
423 Delta = std::max( M_options.rho_increase_big*norm_s_til, Delta );
425 else if ( rho > M_options.rho_small )
426 Delta = std::max( M_options.rho_increase_small*norm_s_til, Delta );
429 norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
435 norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
440 if ( M_solver_stats.collectStats() )
442 M_solver_stats.push( norm_Tgrad_fx,
443 n_CGiter, n_restarts, n_indef,
444 n_crosses_def, n_crosses_indef,
445 n_truss_exit_def, n_truss_exit_indef,
446 Delta_used, ared_til, phi_til, rho );
449 if ( norm_Tgrad_fx > 1e-5 )
450 M_prob.setAccuracy( std::min( 1e-1, norm_Tgrad_fx ) );
452 DVLOG(2) <<
"norm_Tgrad_fx = " << norm_Tgrad_fx <<
"\n";
456 catch ( std::exception
const& __ex )
459 M_prob.evaluate ( x, __fx_new, diff_order<2>() );
461 if ( norm_inf( __fx_new.gradient( 0 ) ) > 1e-10 )
466 lambda_LS( x, __l, __u );
468 if ( M_solver_stats.collectStats() )
469 M_solver_stats.push( x, __l, __u );
474 template<
typename Data,
template<
class>
class Problem>
void
477 M_solver_stats.show();
480 this->lambda_LS( x, lambda_l, lambda_u );
481 M_prob.print_stationary_x( x, lambda_l, lambda_u );
484 template<
typename Data,
template<
class>
class Problem>
486 SolverUnconstrained<Data,Problem>::read_options()
490 std::ifstream fin(
"t_options.inp" );
502 if ( !strcmp( str,
"Delta_init" ) ) M_options.Delta_init = num;
504 if ( !strcmp( str,
"zeta_min" ) ) M_options.zeta_min = num;
506 if ( !strcmp( str,
"CGtol" ) ) M_options.CGtol = num;
508 if ( !strcmp( str,
"deps" ) ) M_options.deps = num;
510 if ( !strcmp( str,
"max_TR_iter" ) ) M_options.max_TR_iter = num;
512 if ( !strcmp( str,
"TR_tol" ) ) M_options.TR_tol = num;
514 if ( !strcmp( str,
"allow_Trust_radius_calculations" ) ) M_options.allow_Trust_radius_calculations = num;
516 if ( !strcmp( str,
"rho_decrease" ) ) M_options.rho_decrease = num;
518 if ( !strcmp( str,
"rho_big" ) ) M_options.rho_big = num;
520 if ( !strcmp( str,
"rho_increase_big" ) ) M_options.rho_increase_big = num;
522 if ( !strcmp( str,
"rho_small" ) ) M_options.rho_small = num;
524 if ( !strcmp( str,
"rho_increase_small" ) ) M_options.rho_increase_small = num;
532 template<
typename Data,
template<
class>
class Problem>
534 SolverUnconstrained<Data,Problem>::makeCauchyStep( vector_type & _x, value_type _Delta,
536 vector_type & _Tgrad_fx,
537 banded_matrix_type & _Hg,
538 symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil,
539 vector_type & _neg_grad_fx )
541 M_prob.evaluate ( _x, __fx, diff_order<2>() );
544 _neg_grad_fx = - __fx.gradient( 0 );
546 makeStep( _x, _neg_grad_fx, _Delta, __fx, _Tgrad_fx,
547 _Hg, _Thess_fxT, _Htil );
550 template<
typename Data,
template<
class>
class Problem>
552 SolverUnconstrained<Data,Problem>::makeStep( vector_type & _x, vector_type & _s, value_type _Delta,
554 vector_type & _Tgrad_fx,
555 banded_matrix_type & _Hg,
556 symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil )
558 M_theta.update( _Delta, _x, _s );
560 _Tgrad_fx = prod( M_theta(), __fx.gradient( 0 ) );
562 banded_matrix_type __diag_grad_fx( _E_n, _E_n, 0, 0 );
564 matrix<value_type> __m( outer_prod( _Tgrad_fx,
565 scalar_vector<value_type>( _Tgrad_fx.size(), 1 ) ) );
566 __diag_grad_fx = banded_adaptor<matrix<value_type> >( __m, 0, 0 );
568 _Hg = prod( M_theta.jacobian(), __diag_grad_fx );
570 _Thess_fxT = prod( M_theta(), prod( __fx.hessian( 0 ), M_theta() ) );
571 _Htil = _Thess_fxT + _Hg;
574 template<
typename Data,
template<
class>
class Problem>
575 typename SolverUnconstrained<Data,Problem>::value_type
576 SolverUnconstrained<Data,Problem>::norm_Theta_x_grad_fx( vector_type & _x, value_type _Delta )
579 M_prob.evaluate ( _x, __fx, diff_order<1>() );
581 vector_type __Tgrad_fx( _E_n );
582 vector_type __neg_grad_fx( _E_n );
583 __neg_grad_fx = - __fx.gradient( 0 );
585 M_theta.update( _Delta, _x, __neg_grad_fx, dsm_type::NO_JACOBIAN );
587 __Tgrad_fx = prod ( M_theta(), __fx.gradient( 0 ) );
589 return norm_2( __Tgrad_fx );
592 template<
typename Data,
template<
class>
class Problem>
593 typename SolverUnconstrained<Data,Problem>::value_type
594 SolverUnconstrained<Data,Problem>::tau( vector_type
const& _s,
595 vector_type
const& _d,
596 value_type
const& _Delta )
598 value_type a = inner_prod( _d, _d );
599 GST_SMART_ASSERT( a != 0 ).error(
"a is 0 and should not be. It will cause a divide by zero." );
600 value_type b = 2.0 * inner_prod( _s, _d );
601 value_type c = inner_prod( _s, _s ) - _Delta*_Delta;
602 return ( -b + sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a );
605 template<
typename Data,
template<
class>
class Problem>
606 typename SolverUnconstrained<Data,Problem>::value_type
607 SolverUnconstrained<Data,Problem>::xi( vector_type
const& __s,
608 vector_type
const& __d,
609 value_type
const& _xi1 )
614 namespace lambda = boost::lambda;
616 value_type __xi = _xi1;
617 std::for_each( __frac.begin(), __frac.end(),
618 lambda::if_then( lambda::_1 > 0 && lambda::_1 < lambda::var( __xi ),
619 lambda::var( __xi ) = lambda::_1 ) );
623 template<
typename Data,
template<
class>
class Problem>
625 SolverUnconstrained<Data,Problem>::CGstep( vector_type & _x, value_type _Delta, vector_type & _sCG,
626 value_type &norm_s_til,
628 int &_n_restarts,
int &_n_indef,
629 int &_n_crosses_def,
int &_n_crosses_indef,
630 int &_n_truss_exit_def,
631 int &_n_truss_exit_indef,
632 value_type &_s_til_x_G_til_x_s_til,
633 value_type &phi_til )
635 bool _restart =
false;
641 _n_crosses_indef = 0;
642 _n_truss_exit_def = 0;
643 _n_truss_exit_indef = 0;
647 vector_type _neg_grad_fx ( _x.size() );
648 vector_type _Tgrad_fx ( _x.size() );
649 vector_type _Htil_d ( _x.size() );
650 vector_type _Htil_s_til ( _x.size() );
652 vector_type _s_til ( _x.size() );
653 vector_type _s_til_old ( _x.size() );
654 vector_type _s_til_eps ( _x.size() );
655 vector_type _rCG ( _x.size() );
656 vector_type _rCG_old ( _x.size() );
657 vector_type _dCG ( _x.size() );
659 value_type _alpha = 0;
666 banded_matrix_type _Hg ( _x.size(), _x.size(), 0, 0 );
668 symmetric_matrix_type _Thess_fxT ( _E_nA, _E_nA );
669 symmetric_matrix_type _Htil ( _E_nA, _E_nA );
671 DVLOG(2) <<
"\n\n[value_type SolverUnconstrained<Data,Problem>::CGstep]...\n";
674 this->makeCauchyStep( _x, _Delta, __fx, _Tgrad_fx, _Hg, _Thess_fxT, _Htil, _neg_grad_fx );
676 DVLOG(2) <<
"Trust region active (C) : " << M_theta.isTrustRegionActive() <<
"\n";
678 value_type _norm_gtil = norm_2( _Tgrad_fx );
679 _s_til = zero_vector<value_type>( _s_til.size() );
685 size_t _inner_iter = 0;
698 _s_til_eps = _s_til + M_options.deps * _dCG;
700 this->makeStep( _x, _s_til_eps, _Delta,
702 _Tgrad_fx, _Hg, _Thess_fxT, _Htil );
704 _norm_gtil = norm_2( _Tgrad_fx );
706 _Htil_s_til = prod( _Htil, _s_til );
709 _rCG = -_Tgrad_fx - _Htil_s_til;
719 _gamma = inner_prod( _dCG, prod( _Htil,_dCG ) );
724 _tau = this->tau( _s_til, _dCG, _Delta );
726 if ( ( !M_theta.isTrustRegionActive() ) || ( _CGiter == 1 ) )
731 else if ( _inner_iter == 0 )
733 _xi = this->xi( _s_til_eps, _dCG, _tau );
738 _xi = this->xi( _s_til, _dCG, _tau );
742 _s_til += _xi * _dCG;
750 else if ( M_theta.isTrustRegionActive() )
754 _n_truss_exit_indef++;
762 if ( !_restart && !_done )
764 _alpha = inner_prod( _rCG, _rCG ) / _gamma;
766 if ( ( !M_theta.isTrustRegionActive() ) || ( _CGiter == 1 ) )
769 else if ( _inner_iter == 0 )
770 _xi = this->xi( _s_til_eps, _dCG, _alpha );
773 _xi = this->xi( _s_til, _dCG, _alpha );
776 _s_til += _xi * _dCG;
778 if ( norm_2( _s_til ) > _Delta )
780 _tau = this->tau( _s_til_old, _dCG, _Delta );
782 _sCG = _s_til_old + _tau * _dCG;
792 if ( !_restart && !_done )
797 _rCG -= _alpha * prod( _Htil, _dCG );
800 else if ( M_theta.isTrustRegionActive() )
807 if ( norm_2( _rCG ) / _norm_gtil < M_options.CGtol )
811 DVLOG(2) <<
"\n\nNormal CG exit 1: ||rCG||/||g_til|| = " << norm_2( _rCG ) / _norm_gtil <<
"\n";
816 if ( _CGiter >= _x.size() )
820 DVLOG(2) <<
"\n\nNormal CG exit 2: _CGiter = " << _CGiter <<
"\n";
828 if ( !_restart && !_done )
830 _beta = inner_prod( _rCG, _rCG ) / inner_prod( _rCG_old, _rCG_old );
832 DVLOG(2) <<
"\nbeta = " << _beta <<
"\n";
834 _dCG = _rCG + _beta * _dCG;
841 _s_til_x_G_til_x_s_til = inner_prod( _sCG, prod( _Hg, _sCG ) );
843 phi_til = inner_prod( _Tgrad_fx, _sCG ) + 0.5 * inner_prod( _sCG, prod( _Htil,_sCG ) );
845 norm_s_til = norm_2( _sCG );
848 _sCG = prod( M_theta(), _sCG );
852 template<
typename Data,
template<
class>
class Problem>
854 SolverUnconstrained<Data,Problem>::lambda_LS( vector_type & _x, vector_type & _lambda_l, vector_type & _lambda_u )
857 matrix_type __A ( _E_nL, _E_nL );
858 vector_type __b ( _E_nL );
861 __A = zero_matrix<double>( _E_nL, _E_nL );
863 for (
int __i = 0; __i < _E_n; ++__i )
865 __A( __i, __i ) = 1.0 + ( M_prob.x_l( __i )-_x ( __i ) )*( M_prob.x_l( __i )-_x ( __i ) );
867 __A( __i, _E_n+__i ) = -1.0;
868 __A( _E_n+__i,__i ) = -1.0;
870 __A( ( _E_n+__i ), _E_n+__i ) = 1.0+( _x ( __i )-M_prob.x_u( __i ) )*( _x ( __i )-M_prob.x_u( __i ) );
875 M_prob.evaluate( _x, __f_x, diff_order<1>() );
877 for (
int i = 0; i < _E_n; i++ )
879 value_type
val = __f_x.gradient( 0, i );
881 __b ( _E_n+i ) = -
val;
892 ublas::vector<int> __itype( __N );
894 dgesv_( &__N,&__nrhs,__A.data().begin(),&__N,__itype.data().begin(),__b.data().begin(),&__N,&__info );
898 std::ostringstream __err;
900 __err <<
"[" << __PRETTY_FUNCTION__ <<
"] dppsv failed: " << __info <<
"\n"
901 <<
"A = " << __A <<
"\n"
902 <<
"B = " << __b <<
"\n";
904 throw std::out_of_range( __err.str() );
907 for (
int i = 0; i < _E_n; i++ )
909 _lambda_l ( i ) = __b ( i );
910 _lambda_u ( i ) = __b ( _E_n+i );
919 template<
typename Data,
template<
class>
class Problem>
923 norm_Tgrad_fx_hstr.clear();
925 n_CGiter_hstr.clear();
926 n_restarts_hstr.clear();
927 n_indef_hstr.clear();
928 n_crosses_def_hstr.clear();
929 n_crosses_indef_hstr.clear();
930 n_truss_exit_def_hstr.clear();
931 n_truss_exit_indef_hstr.clear();
933 ared_til_hstr.clear();
934 phi_til_hstr.clear();
937 template<
typename Data,
template<
class>
class Problem>
946 template<
typename Data,
template<
class>
class Problem>
949 int n_CGiter,
int n_restarts,
int n_indef,
950 int n_crosses_def,
int n_crosses_indef,
951 int n_truss_exit_def,
int n_truss_exit_indef,
958 norm_Tgrad_fx_hstr.push_back( norm_Tgrad_fx );
959 n_CGiter_hstr.push_back( n_CGiter );
960 n_restarts_hstr.push_back( n_restarts );
961 n_indef_hstr.push_back( n_indef );
962 n_crosses_def_hstr.push_back( n_crosses_def );
963 n_crosses_indef_hstr.push_back( n_crosses_indef );
964 n_truss_exit_def_hstr.push_back( n_truss_exit_def );
965 n_truss_exit_indef_hstr.push_back( n_truss_exit_indef );
966 Delta_hstr.push_back( Delta );
967 ared_til_hstr.push_back( ared_til );
968 phi_til_hstr.push_back( phi_til );
969 rho_hstr.push_back( rho );
972 template<
typename Data,
template<
class>
class Problem>
974 SolverUnconstrained<Data,Problem>::Stats::show()
const
976 M_s->problem().print_stationary_x ( M_x, M_l, M_u );
978 iter = norm_Tgrad_fx_hstr.size();
979 std::cerr <<
"\n\nScaled Trust-Region Method Statistics:\n";
980 std::cerr <<
"\nNumber of outter iterations: " << iter-1 <<
"\n";
983 std::cerr <<
" k CGiters restarts indefs crosses trustEx ||g_til|| Delta^k ared_til^k phi_til^k rho^k \n";
984 std::cerr <<
"--------------------------------------------------------------------------------------------------------\n";
986 for (
int k = 0; k < iter; k++ )
988 fprintf( stderr,
"%3d ", k );
991 std::cerr <<
" N/A N/A N/A N/A N/A";
995 std::cerr <<
" " << n_CGiter_hstr[k]
996 <<
" " << n_restarts_hstr[k]
997 <<
" " << n_indef_hstr[k]
998 <<
" " << n_crosses_def_hstr[k] + n_crosses_indef_hstr[k]
999 <<
" " << n_truss_exit_def_hstr[k] + n_truss_exit_indef_hstr[k] <<
" ";
1002 fprintf ( stderr,
" %13.3e %9f", norm_Tgrad_fx_hstr[ k ], Delta_hstr[k] );
1005 std::cerr <<
" N/A N/A N/A";
1009 fprintf ( stderr,
" %11.3e", ared_til_hstr[ k ] );
1010 fprintf ( stderr,
" %11.3e", phi_til_hstr[ k ] );
1011 fprintf ( stderr,
" %11f", rho_hstr[ k ] );
1014 std::cerr << std::endl;