Go to the documentation of this file.00001
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "vnl_lbfgs.h"
00014 #include <vcl_cmath.h>
00015 #include <vcl_iostream.h>
00016 #include <vcl_iomanip.h>
00017
00018 #include <vnl/algo/vnl_netlib.h>
00019
00020
00021
00022
00023 vnl_lbfgs::vnl_lbfgs():
00024 f_(0)
00025 {
00026 init_parameters();
00027 }
00028
00029
00030
00031 vnl_lbfgs::vnl_lbfgs(vnl_cost_function& f):
00032 f_(&f)
00033 {
00034 init_parameters();
00035 }
00036
00037
00038
00039 void vnl_lbfgs::init_parameters()
00040 {
00041 memory = 5;
00042 line_search_accuracy = 0.9;
00043 default_step_length = 1.0;
00044 }
00045
00046 bool vnl_lbfgs::minimize(vnl_vector<double>& x)
00047 {
00048
00049
00050
00051 long n = f_->get_number_of_unknowns();
00052 long m = memory;
00053
00054
00055
00056
00057 v3p_netlib_lbfgs_global_t lbfgs_global;
00058 v3p_netlib_lbfgs_init(&lbfgs_global);
00059
00060 long iprint[2] = {1, 0};
00061 vnl_vector<double> g(n);
00062
00063
00064 vnl_vector<double> diag(n);
00065
00066 vnl_vector<double> w(n * (2*m+1)+2*m);
00067
00068 if (verbose_)
00069 vcl_cerr << "vnl_lbfgs: n = "<< n <<", memory = "<< m <<", Workspace = "
00070 << w.size() << "[ "<< ( w.size() / 128.0 / 1024.0) <<" MB], ErrorScale = "
00071 << f_->reported_error(1) <<", xnorm = "<< x.magnitude() << vcl_endl;
00072
00073 bool we_trace = (verbose_ && !trace);
00074
00075 if (we_trace)
00076 vcl_cerr << "vnl_lbfgs: ";
00077
00078 double best_f = 0;
00079 vnl_vector<double> best_x;
00080
00081 bool ok;
00082 this->num_evaluations_ = 0;
00083 this->num_iterations_ = 0;
00084 long iflag = 0;
00085 while (true) {
00086
00087 v3p_netlib_logical diagco = false;
00088
00089
00090 double eps = gtol;
00091 double local_xtol = 1e-16;
00092 lbfgs_global.gtol = line_search_accuracy;
00093 lbfgs_global.stpinit = default_step_length;
00094
00095
00096 double f;
00097 f_->compute(x, &f, &g);
00098 if (this->num_evaluations_ == 0) {
00099 this->start_error_ = f;
00100 best_f = f;
00101 } else if (f < best_f) {
00102 best_x = x;
00103 best_f = f;
00104 }
00105
00106 #define print_(i,a,b,c,d) vcl_cerr<<vcl_setw(6)<<i<<' '<<vcl_setw(20)<<a<<' '\
00107 <<vcl_setw(20)<<b<<' '<<vcl_setw(20)<<c<<' '<<vcl_setw(20)<<d<<'\n'
00108
00109 if (check_derivatives_)
00110 {
00111 vcl_cerr << "vnl_lbfgs: f = " << f_->reported_error(f) << ", computing FD gradient\n";
00112 vnl_vector<double> fdg = f_->fdgradf(x);
00113 if (verbose_)
00114 {
00115 int l = n;
00116 int limit = 100;
00117 int limit_tail = 10;
00118 if (l > limit + limit_tail) {
00119 vcl_cerr << " [ Showing only first " <<limit<< " components ]\n";
00120 l = limit;
00121 }
00122 print_("i","x","g","fdg","dg");
00123 print_("-","-","-","---","--");
00124 for (int i = 0; i < l; ++i)
00125 print_(i, x[i], g[i], fdg[i], g[i]-fdg[i]);
00126 if (n > limit) {
00127 vcl_cerr << " ...\n";
00128 for (int i = n - limit_tail; i < n; ++i)
00129 print_(i, x[i], g[i], fdg[i], g[i]-fdg[i]);
00130 }
00131 }
00132 vcl_cerr << " ERROR = " << (fdg - g).squared_magnitude() / vcl_sqrt(double(n)) << "\n";
00133 }
00134
00135 iprint[0] = trace ? 1 : -1;
00136 iprint[1] = 0;
00137 v3p_netlib_lbfgs_(
00138 &n, &m, x.data_block(), &f, g.data_block(), &diagco, diag.data_block(),
00139 iprint, &eps, &local_xtol, w.data_block(), &iflag, &lbfgs_global);
00140
00141 this->report_eval(f);
00142
00143 if (this->report_iter()) {
00144 failure_code_ = FAILED_USER_REQUEST;
00145 ok = false;
00146 x = best_x;
00147 break;
00148 }
00149
00150 if (we_trace)
00151 vcl_cerr << iflag << ":" << f_->reported_error(f) << " ";
00152
00153 if (iflag == 0) {
00154
00155 this->end_error_ = f;
00156 ok = true;
00157 x = best_x;
00158 break;
00159 }
00160
00161 if (iflag < 0) {
00162
00163 vcl_cerr << "vnl_lbfgs: Error. Netlib routine lbfgs failed.\n";
00164 ok = false;
00165 x = best_x;
00166 break;
00167 }
00168
00169 if (this->num_evaluations_ > get_max_function_evals()) {
00170 failure_code_ = TOO_MANY_ITERATIONS;
00171 ok = false;
00172 x = best_x;
00173 break;
00174 }
00175
00176 }
00177 if (we_trace) vcl_cerr << "done\n";
00178
00179 return ok;
00180 }