IT++ Logo
newton_search.cpp
Go to the documentation of this file.
00001 
00029 #include <itpp/optim/newton_search.h>
00030 #include <itpp/base/specmat.h>
00031 #include <itpp/stat/misc_stat.h>
00032 
00033 
00034 namespace itpp
00035 {
00036 
00037 
00038 Newton_Search::Newton_Search()
00039 {
00040   method = BFGS;
00041 
00042   initial_stepsize = 1.0;
00043   stop_epsilon_1 = 1e-4;
00044   stop_epsilon_2 = 1e-8;
00045   max_evaluations = 100;
00046 
00047   f = NULL;
00048   df_dx = NULL;
00049 
00050   no_feval = 0;
00051   init = false;
00052   finished = false;
00053   trace = false;
00054 }
00055 
00056 void Newton_Search::set_function(double(*function)(const vec&))
00057 {
00058   // Add checks to see that function is OK???
00059   f = function;
00060 }
00061 
00062 void Newton_Search::set_gradient(vec(*gradient)(const vec&))
00063 {
00064   // Add checks to see that function is OK???
00065   df_dx = gradient;
00066 }
00067 
00068 void Newton_Search::set_start_point(const vec &x, const mat &D)
00069 {
00070   // check that parameters are valid???
00071   x_start = x;
00072   n = x.size();
00073   D_start = D;
00074 
00075   finished = false;
00076   init = true;
00077 }
00078 
00079 void Newton_Search::set_start_point(const vec &x)
00080 {
00081   // check that parameters are valid???
00082   x_start = x;
00083   n = x.size();
00084   D_start = eye(n);
00085 
00086   finished = false;
00087   init = true;
00088 }
00089 
00090 bool Newton_Search::search()
00091 {
00092   // Check parameters and function call ???
00093   // check that x_start is a valid point, not a NaN and that norm(x0) is not inf
00094 
00095   it_assert(f != NULL, "Newton_Search: Function pointer is not set");
00096   it_assert(df_dx != NULL, "Newton_Search: Gradient function pointer is not set");
00097 
00098   it_assert(init, "Newton_Search: Starting point is not set");
00099 
00100 
00101   F = f(x_start); // function initial value
00102   vec g = df_dx(x_start); // gradient initial value
00103   vec x = x_start;
00104   no_feval++;
00105 
00106   finished = false;
00107 
00108   // Initial inverse Hessian, D
00109   mat D = D_start;
00110 
00111 
00112   bool fst = true; // what is this???
00113 
00114   bool stop = false;
00115 
00116   // Finish initialization
00117   no_iter = 0;
00118   ng = max(abs(g)); // norm(g,inf)
00119 
00120   double Delta = initial_stepsize;
00121   nh = 0; // what is this???
00122   vec h;
00123 
00124   if (trace) { // prepare structures to store trace data
00125     x_values.set_size(max_evaluations);
00126     F_values.set_size(max_evaluations);
00127     ng_values.set_size(max_evaluations);
00128     Delta_values.set_size(max_evaluations);
00129   }
00130 
00131   Line_Search ls;
00132   ls.set_functions(f, df_dx);
00133 
00134   if (ng <= stop_epsilon_1)
00135     stop = true;
00136   else {
00137     h = zeros(n);
00138     nh = 0;
00139     ls.set_stop_values(0.05, 0.99);
00140     ls.set_max_iterations(5);
00141     ls.set_max_stepsize(2);
00142   }
00143 
00144   bool more = true; //???
00145 
00146   while (!stop && more) {
00147     vec h, w, y, v;
00148     double yh, yv, a;
00149 
00150     // Previous values
00151     vec xp = x, gp = g;
00152     // double Fp = F;           ### 2006-02-03 by ediap: Unused variable!
00153     double nx = norm(x);
00154 
00155     h = D * (-g);
00156     nh = norm(h);
00157     bool red = false;
00158 
00159     if (nh <= stop_epsilon_2*(stop_epsilon_2 + nx))  // stop criterion
00160       stop = true;
00161     else {
00162       if (fst || nh > Delta) {  // Scale to ||h|| = Delta
00163         h = (Delta / nh) * h;
00164         nh = Delta;
00165         fst = false;
00166         red = true;
00167       }
00168       //  Line search
00169       ls.set_start_point(x, F, g, h);
00170       more = ls.search(x, F, g);
00171       no_feval = no_feval + ls.get_no_function_evaluations();
00172 
00173       if (more == false) {  // something wrong in linesearch?
00174         x_end = x;
00175         return false;
00176       }
00177       else {
00178         if (ls.get_alpha() < 1)   // Reduce Delta
00179           Delta = .35 * Delta;
00180         else if (red && (ls.get_slope_ratio() > .7))   // Increase Delta
00181           Delta = 3 * Delta;
00182 
00183         //  Update ||g||
00184         ng = max(abs(g)); // norm(g,inf);
00185 
00186         if (trace) { // store trace
00187           x_values(no_iter) = x;
00188           F_values(no_iter) = F;
00189           ng_values(no_iter) = ng;
00190           Delta_values(no_iter) = Delta;
00191         }
00192 
00193         no_iter++;
00194         h = x - xp;
00195         nh = norm(h);
00196 
00197         //if  (nh == 0)
00198         //  found = 4;
00199         //else {
00200         y = g - gp;
00201         yh = dot(y, h);
00202         if (yh > std::sqrt(eps) * nh * norm(y)) {
00203           //  Update  D
00204           v = D * y;
00205           yv = dot(y, v);
00206           a = (1 + yv / yh) / yh;
00207           w = (a / 2) * h - v / yh;
00208           D += outer_product(w, h) + outer_product(h, w); //D = D + w*h' + h*w';
00209         }  // update D
00210         //  Check stopping criteria
00211         double thrx = stop_epsilon_2 * (stop_epsilon_2 + norm(x));
00212         if (ng <= stop_epsilon_1)
00213           stop = true; // stop = 1, stop by small gradient
00214         else if (nh <= thrx)
00215           stop = true; // stop = 2, stop by small x-step
00216         else if (no_feval >= max_evaluations)
00217           stop = true; // stop = 3, number of function evaluations exeeded
00218         else
00219           Delta = std::max(Delta, 2 * thrx);
00220         //} found =4
00221       }  // Nonzero h
00222     } // nofail
00223   }  // iteration
00224 
00225   //  Set return values
00226   x_end = x;
00227   finished = true;
00228 
00229   if (trace) { // trim size of trace output
00230     x_values.set_size(no_iter, true);
00231     F_values.set_size(no_iter, true);
00232     ng_values.set_size(no_iter, true);
00233     Delta_values.set_size(no_iter, true);
00234   }
00235 
00236   return true;
00237 }
00238 
00239 bool Newton_Search::search(vec &xn)
00240 {
00241   bool state = search();
00242   xn = get_solution();
00243   return state;
00244 }
00245 
00246 bool Newton_Search::search(const vec &x0, vec &xn)
00247 {
00248   set_start_point(x0);
00249   bool state = search();
00250   xn = get_solution();
00251   return state;
00252 }
00253 
00254 vec Newton_Search::get_solution()
00255 {
00256   it_assert(finished, "Newton_Search: search is not run yet");
00257   return x_end;
00258 }
00259 
00260 double Newton_Search::get_function_value()
00261 {
00262   if (finished)
00263     return F;
00264   else
00265     it_warning("Newton_Search::get_function_value, search has not been run");
00266 
00267   return 0.0;
00268 }
00269 
00270 double Newton_Search::get_stop_1()
00271 {
00272   if (finished)
00273     return ng;
00274   else
00275     it_warning("Newton_Search::get_stop_1, search has not been run");
00276 
00277   return 0.0;
00278 }
00279 
00280 double Newton_Search::get_stop_2()
00281 {
00282   if (finished)
00283     return nh;
00284   else
00285     it_warning("Newton_Search::get_stop_2, search has not been run");
00286 
00287   return 0.0;
00288 }
00289 
00290 int Newton_Search::get_no_iterations()
00291 {
00292   if (finished)
00293     return no_iter;
00294   else
00295     it_warning("Newton_Search::get_no_iterations, search has not been run");
00296 
00297   return 0;
00298 }
00299 
00300 int Newton_Search::get_no_function_evaluations()
00301 {
00302   if (finished)
00303     return no_feval;
00304   else
00305     it_warning("Newton_Search::get_no_function_evaluations, search has not been run");
00306 
00307   return 0;
00308 }
00309 
00310 
00311 void Newton_Search::get_trace(Array<vec> & xvalues, vec &Fvalues, vec &ngvalues, vec &dvalues)
00312 {
00313   if (finished) {
00314     if (trace) { // trim size of trace output
00315       xvalues = x_values;
00316       Fvalues = F_values;
00317       ngvalues = ng_values;
00318       dvalues = Delta_values;
00319     }
00320     else
00321       it_warning("Newton_Search::get_trace, trace is not enabled");
00322   }
00323   else
00324     it_warning("Newton_Search::get_trace, search has not been run");
00325 }
00326 
00327 //================================== Line_Search =============================================
00328 
00329 Line_Search::Line_Search()
00330 {
00331   method = Soft;
00332 
00333   if (method == Soft) {
00334     stop_rho = 1e-3;
00335     stop_beta = 0.99;
00336   }
00337 
00338   max_iterations = 10;
00339   max_stepsize = 10;
00340 
00341   f = NULL;
00342   df_dx = NULL;
00343   no_feval = 0;
00344   init = false;
00345   finished = false;
00346   trace = false;
00347 }
00348 
00349 void Line_Search::set_function(double(*function)(const vec&))
00350 {
00351   // Add checks to see that function is OK???
00352   f = function;
00353 }
00354 
00355 void Line_Search::set_gradient(vec(*gradient)(const vec&))
00356 {
00357   // Add checks to see that function is OK???
00358   df_dx = gradient;
00359 }
00360 
00361 
00362 void Line_Search::set_stop_values(double rho, double beta)
00363 {
00364   // test input values???
00365   stop_rho = rho;
00366   stop_beta = beta;
00367 }
00368 
00369 
00370 void Line_Search::set_start_point(const vec &x, double F, const vec &g, const vec &h)
00371 {
00372   // check values ???
00373   x_start = x;
00374   F_start = F;
00375   g_start = g;
00376   h_start = h;
00377   n = x.size();
00378 
00379   finished = false;
00380   init = true;
00381 }
00382 
00383 void Line_Search::get_solution(vec &xn, double &Fn, vec &gn)
00384 {
00385   it_assert(finished, "Line_Search: search is not run yet");
00386 
00387   xn = x_end;
00388   Fn = F_end;
00389   gn = g_end;
00390 }
00391 
00392 bool Line_Search::search()
00393 {
00394   it_assert(f != NULL, "Line_Search: Function pointer is not set");
00395   it_assert(df_dx != NULL, "Line_Search: Gradient function pointer is not set");
00396 
00397   it_assert(init, "Line_search: Starting point is not set");
00398 
00399   // Default return values and simple checks
00400   x_end = x_start;
00401   F_end = F_start;
00402   g_end = g_start;
00403 
00404   // add some checks???
00405   finished = false;
00406 
00407   vec g;
00408 
00409   // return parameters
00410   no_feval = 0;
00411   slope_ratio = 1;
00412 
00413 
00414 
00415   // Check descent condition
00416   double dF0 = dot(h_start, g_end);
00417 
00418   if (trace) { // prepare structures to store trace data
00419     alpha_values.set_size(max_iterations);
00420     F_values.set_size(max_iterations);
00421     dF_values.set_size(max_iterations);
00422     alpha_values(0) = 0;
00423     F_values(0) = F_end;
00424     dF_values(0) = dF0;
00425   }
00426 
00427 
00428   if (dF0 >= -10*eps*norm(h_start)*norm(g_end)) {  // not significantly downhill
00429     if (trace) { // store trace
00430       alpha_values.set_size(1, true);
00431       F_values.set_size(1, true);
00432       dF_values.set_size(1, true);
00433     }
00434     return false;
00435   }
00436 
00437   // Finish initialization
00438   double F0 = F_start, slope0, slopethr;
00439 
00440   if (method == Soft) {
00441     slope0 = stop_rho * dF0;
00442     slopethr = stop_beta * dF0;
00443   }
00444   else { // exact line search
00445     slope0 = 0;
00446     slopethr = stop_rho * std::abs(dF0);
00447   }
00448 
00449   // Get an initial interval for am
00450   double a = 0, Fa = F_end, dFa = dF0;
00451   bool stop = false;
00452   double b = std::min(1.0, max_stepsize), Fb = 0, dFb = 0;
00453 
00454 
00455   while (!stop) {
00456     Fb = f(x_start + b * h_start);
00457     g = df_dx(x_start + b * h_start);
00458     // check if these values are OK if not return false???
00459     no_feval++;
00460 
00461     dFb = dot(g, h_start);
00462     if (trace) { // store trace
00463       alpha_values(no_feval) = b;
00464       F_values(no_feval) = Fb;
00465       dF_values(no_feval) = dFb;
00466     }
00467 
00468     if (Fb < F0 + slope0*b) {  // new lower bound
00469       alpha = b;
00470       slope_ratio = dFb / dF0; // info(2);
00471 
00472       if (method == Soft) {
00473         a = b;
00474         Fa = Fb;
00475         dFa = dFb;
00476       }
00477 
00478       x_end = x_start + b * h_start;
00479       F_end = Fb;
00480       g_end = g;
00481 
00482       if ((dFb < std::min(slopethr, 0.0)) && (no_feval < max_iterations) && (b < max_stepsize)) {
00483         // Augment right hand end
00484         if (method == Exact) {
00485           a = b;
00486           Fa = Fb;
00487           dFa = dFb;
00488         }
00489         if (2.5*b >= max_stepsize)
00490           b = max_stepsize;
00491         else
00492           b = 2 * b;
00493       }
00494       else
00495         stop = true;
00496     }
00497     else
00498       stop = true;
00499   } // phase 1: expand interval
00500 
00501 
00502 
00503   if (stop)  // OK so far.  Check stopping criteria
00504     stop = (no_feval >= max_iterations)
00505            || (b >= max_stepsize && dFb < slopethr)
00506            || (a > 0 && dFb >= slopethr);
00507   // Commented by ediap 2006-07-17: redundant check
00508   //  || ( (method == Soft) && (a > 0 & dFb >= slopethr) );  // OK
00509 
00510 
00511   if (stop && trace) {
00512     alpha_values.set_size(no_feval, true);
00513     F_values.set_size(no_feval, true);
00514     dF_values.set_size(no_feval, true);
00515   }
00516 
00517   // Refine interval
00518   while (!stop) {
00519 
00520     double c, Fc, dFc;
00521 
00522     //c = interpolate(xfd,n);
00523     double C = Fb - Fa - (b - a) * dFa;
00524     if (C >= 5*n*eps*b) {
00525       double A = a - 0.5 * dFa * (sqr(b - a) / C);
00526       c = std::min(std::max(a + 0.1 * (b - a), A), b - 0.1 * (b - a));  // % Ensure significant resuction
00527     }
00528     else
00529       c = (a + b) / 2;
00530 
00531     Fc = f(x_start + c * h_start);
00532     g = df_dx(x_start + c * h_start);
00533     dFc = dot(g, h_start);
00534     // check these values???
00535     no_feval++;
00536 
00537     if (trace) { // store trace
00538       alpha_values(no_feval) = c;
00539       F_values(no_feval) = Fc;
00540       dF_values(no_feval) = dFc;
00541     }
00542 
00543     if (method == Soft) {
00544       // soft line method
00545       if (Fc < F0 + slope0*c) {  // new lower bound
00546         alpha = c;
00547         slope_ratio = dFc / dF0;
00548 
00549         x_end = x_start + c * h_start;
00550         F_end = Fc;
00551         g_end = g;
00552         a = c;
00553         Fa = Fc;
00554         dFa = dFc; // xfd(:,1) = xfd(:,3);
00555         stop = (dFc > slopethr);
00556       }
00557       else { // new upper bound
00558         b = c;
00559         Fb = Fc;
00560         dFb = dFc; // xfd(:,2) = xfd(:,3);
00561       }
00562 
00563     }
00564     else { // Exact line search
00565       if (Fc < F_end) {  // better approximant
00566         alpha = c;
00567         slope_ratio = dFc / dF0;
00568         x_end = x_start + c * h_start;
00569         F_end = Fc;
00570         g_end = g;
00571       }
00572       if (dFc < 0) {  // new lower bound
00573         a = c;
00574         Fa = Fc;
00575         dFa = dFc; // xfd(:,1) = xfd(:,3);
00576       }
00577       else { //new upper bound
00578         b = c;
00579         Fb = Fc;
00580         dFb = dFc; // xfd(:,2) = xfd(:,3);
00581       }
00582       stop = (std::abs(dFc) <= slopethr) | ((b - a) < stop_beta * b);
00583     }
00584 
00585     stop = (stop | (no_feval >= max_iterations));
00586   } // refine
00587 
00588   finished = true;
00589 
00590   if (trace) { // store trace
00591     alpha_values.set_size(no_feval + 1, true);
00592     F_values.set_size(no_feval + 1, true);
00593     dF_values.set_size(no_feval + 1, true);
00594   }
00595 
00596   return true;
00597 }
00598 
00599 bool Line_Search::search(vec &xn, double &Fn, vec &gn)
00600 {
00601   bool state = search();
00602   get_solution(xn, Fn, gn);
00603   return state;
00604 }
00605 
00606 bool Line_Search::search(const vec &x, double F, const vec &g, const vec &h,
00607                          vec &xn, double &Fn, vec &gn)
00608 {
00609   set_start_point(x, F, g, h);
00610   bool state = search();
00611   get_solution(xn, Fn, gn);
00612   return state;
00613 }
00614 
00615 
00616 double Line_Search::get_alpha()
00617 {
00618   if (finished)
00619     return alpha;
00620   else
00621     it_warning("Line_Search::get_alpha, search has not been run");
00622 
00623   return 0.0;
00624 }
00625 
00626 double Line_Search::get_slope_ratio()
00627 {
00628   if (finished)
00629     return slope_ratio;
00630   else
00631     it_warning("Line_Search::get_slope_raio, search has not been run");
00632 
00633   return 0.0;
00634 }
00635 
00636 int Line_Search::get_no_function_evaluations()
00637 {
00638   if (finished)
00639     return no_feval;
00640   else
00641     it_warning("Line_Search::get_no_function_evaluations, search has not been run");
00642 
00643   return 0;
00644 }
00645 
00646 
00647 void Line_Search::set_max_iterations(int value)
00648 {
00649   it_assert(value > 0, "Line_Search, max iterations must be > 0");
00650   max_iterations = value;
00651 }
00652 
00653 void Line_Search::set_max_stepsize(double value)
00654 {
00655   it_assert(value > 0, "Line_Search, max stepsize must be > 0");
00656   max_stepsize = value;
00657 }
00658 
00659 void Line_Search::set_method(const Line_Search_Method &search_method)
00660 {
00661   method = search_method;
00662 
00663   if (method == Soft) {
00664     stop_rho = 1e-3;
00665     stop_beta = 0.99;
00666   }
00667   else { // exact line search
00668     method = Exact;
00669     stop_rho = 1e-3;
00670     stop_beta = 1e-3;
00671   }
00672 }
00673 
00674 
00675 void Line_Search::get_trace(vec &alphavalues, vec &Fvalues, vec &dFvalues)
00676 {
00677   if (finished) {
00678     if (trace) { // trim size of trace output
00679       alphavalues = alpha_values;
00680       Fvalues = F_values;
00681       dFvalues = dF_values;
00682     }
00683     else
00684       it_warning("Line_Search::get_trace, trace is not enabled");
00685   }
00686   else
00687     it_warning("Line_Search::get_trace, search has not been run");
00688 }
00689 
00690 // =========================== functions ==============================================
00691 
00692 vec fminunc(double(*function)(const vec&), vec(*gradient)(const vec&), const vec &x0)
00693 {
00694   Newton_Search newton;
00695   newton.set_functions(function, gradient);
00696 
00697   vec xn;
00698   newton.search(x0, xn);
00699 
00700   return xn;
00701 }
00702 
00703 
00704 
00705 } // namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:32 for IT++ by Doxygen 1.7.4