IT++ Logo Newcom Logo

newton_search.cpp

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

Generated on Fri Jun 8 00:37:33 2007 for IT++ by Doxygen 1.5.2