IT++ Logo

newton_search.cpp

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

Generated on Sun Jul 26 08:54:57 2009 for IT++ by Doxygen 1.5.9