core/vnl/algo/vnl_lbfgsb.cxx
Go to the documentation of this file.
00001 // This is core/vnl/algo/vnl_lbfgsb.cxx
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005 //:
00006 // \file
00007 //
00008 // \author Brad King, Kitware Inc.
00009 // \date   28 Aug 07
00010 //
00011 //-----------------------------------------------------------------------------
00012 
00013 #include "vnl_lbfgsb.h"
00014 #include <vcl_cstring.h>
00015 #include <vcl_iostream.h>
00016 
00017 #include <vnl/algo/vnl_netlib.h> // setulb_()
00018 
00019 //----------------------------------------------------------------------------
00020 vnl_lbfgsb::vnl_lbfgsb(): f_(0)
00021 {
00022   init_parameters();
00023 }
00024 
00025 //----------------------------------------------------------------------------
00026 vnl_lbfgsb::vnl_lbfgsb(vnl_cost_function& f): f_(&f)
00027 {
00028   init_parameters();
00029 }
00030 
00031 //----------------------------------------------------------------------------
00032 void vnl_lbfgsb::init_parameters()
00033 {
00034   long n = this->f_->get_number_of_unknowns();
00035   this->bound_selection_.set_size(n);
00036   this->bound_selection_.fill(0);
00037   this->max_corrections_ = 5;
00038   this->convergence_factor_ = 1e+7;
00039   this->projected_gradient_tolerance_ = 1e-5;
00040 }
00041 
00042 //----------------------------------------------------------------------------
00043 bool vnl_lbfgsb::minimize(vnl_vector<double>& x)
00044 {
00045   // Basic setup.
00046   long n = this->f_->get_number_of_unknowns();
00047   long m = this->max_corrections_;
00048 
00049   // Function and gradient.
00050   double f = 0;
00051   vnl_vector<double> gradient(n);
00052 
00053   // Working space.
00054   vnl_vector<double> wa((2*m+4)*n + 12*m*m + 12*m);
00055   vnl_vector<long> iwa(3*n);
00056   char csave[60];
00057   long lsave[4];
00058   long isave[44];
00059   double dsave[29];
00060 
00061   // Task communication.
00062   char task[61]="START                                                       ";
00063 
00064   // Verbosity level inside lbfgs implementation.
00065   // (-1 no o/p, 0 start and end, 1 every iter)
00066   long const iprint = trace ? 1 : -1;
00067 
00068   // Initialize iteration.
00069   this->num_evaluations_ = 0;
00070   this->num_iterations_ = 0;
00071 
00072   // TODO: Deal with verbose_, check_derivatives_, trace, xtol,
00073   // maxfev, ftol, gtol, epsfcn members of vnl_nonlinear_minimizer.
00074 
00075   // Track the best position found.
00076   vnl_vector<double> x_best(x);
00077 
00078   bool ok = true;
00079   for (;;)
00080   {
00081     // Call the L-BFGS-B code.
00082     v3p_netlib_setulb_(
00083       &n,
00084       &m,
00085       x.data_block(),
00086       this->lower_bound_.data_block(),
00087       this->upper_bound_.data_block(),
00088       this->bound_selection_.data_block(),
00089       &f, gradient.data_block(),
00090       &this->convergence_factor_,
00091       &this->projected_gradient_tolerance_,
00092       wa.data_block(),
00093       iwa.data_block(),
00094       task,
00095       &iprint,
00096       csave, lsave, isave, dsave
00097       );
00098 
00099     // Check the current task.
00100     if (vcl_strncmp("FG", task, 2) == 0)
00101     {
00102       // Evaluate the function and gradient.
00103       this->f_->compute(x, &f, &gradient);
00104 
00105       if (this->num_evaluations_ == 0)
00106       {
00107         x_best = x;
00108         this->start_error_ = f;
00109         this->end_error_ = f;
00110       }
00111       else if (f < this->end_error_)
00112       {
00113         x_best = x;
00114         this->end_error_ = f;
00115       }
00116       this->report_eval(f);
00117     }
00118     else if (vcl_strncmp("NEW_X", task, 5) == 0)
00119     {
00120       // dsave[12] = the infinity norm of the projected gradient
00121       this->inf_norm_projected_gradient_ = dsave[12];
00122 
00123       // Iteration.a
00124       if (this->report_iter())
00125       {
00126         this->failure_code_ = FAILED_USER_REQUEST;
00127         ok = false;
00128         break;
00129       }
00130     }
00131     else if (vcl_strncmp("ERROR", task, 5) == 0)
00132     {
00133       // some error
00134       this->failure_code_ = ERROR_FAILURE;
00135       ok = false;
00136       break;
00137     }
00138     else if (vcl_strncmp("CONVERGENCE", task, 11) == 0)
00139     {
00140       // convergence has been reached
00141       if (f < this->end_error_)
00142       {
00143         x_best = x;
00144         this->end_error_ = f;
00145       }
00146 
00147       if (vcl_strncmp("CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH",
00148                       task, 47) == 0)
00149       {
00150         // function tolerance reached
00151         this->failure_code_ = CONVERGED_FTOL;
00152       }
00153       else if (vcl_strncmp("CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL",
00154                            task, 48) == 0)
00155       {
00156         // gradient tolerance reached
00157         this->failure_code_ = CONVERGED_GTOL;
00158       }
00159       else
00160       {
00161         this->failure_code_ = ERROR_FAILURE;
00162         if (trace)
00163         {
00164           vcl_cerr << "Unknown convergence type: " << task << vcl_endl;
00165         }
00166       }
00167       break;
00168     }
00169     else
00170     {
00171       // unknown task
00172       this->failure_code_ = ERROR_FAILURE;
00173       if (trace)
00174       {
00175         vcl_cerr << "Unknown failure with task: " << task << vcl_endl;
00176       }
00177       ok = false;
00178       break;
00179     }
00180 
00181     if (this->num_evaluations_ > this->get_max_function_evals())
00182     {
00183       this->failure_code_ = TOO_MANY_ITERATIONS;
00184       ok = false;
00185       break;
00186     }
00187   }
00188 
00189   // Store the best known position no matter the outcome.
00190   x = x_best;
00191 
00192   return ok;
00193 }