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   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   } else {
00177     if  (ls.get_alpha() < 1)  // Reduce Delta
00178       Delta = .35 * Delta;
00179     else if  (red && (ls.get_slope_ratio() > .7))  // Increase Delta
00180       Delta = 3*Delta;
00181 
00182     //  Update ||g||
00183     ng = max(abs(g)); // norm(g,inf);
00184 
00185     if (trace) { // store trace
00186       x_values(no_iter) = x;
00187       F_values(no_iter) = F;
00188       ng_values(no_iter) = ng;
00189       Delta_values(no_iter) = Delta;
00190     }
00191 
00192     no_iter++;
00193     h = x - xp;
00194     nh = norm(h);
00195 
00196     //if  (nh == 0)
00197     //  found = 4;
00198     //else {
00199     y = g - gp;
00200     yh = dot(y,h);
00201     if  (yh > std::sqrt(eps) * nh * norm(y)) {
00202       //  Update  D
00203       v = D*y;
00204       yv = dot(y,v);
00205       a = (1 + yv/yh)/yh;
00206       w = (a/2)*h - v/yh;
00207       D += outer_product(w,h) + outer_product(h,w); //D = D + w*h' + h*w';
00208     }  // update D
00209     //  Check stopping criteria
00210     double thrx = stop_epsilon_2*(stop_epsilon_2 + norm(x));
00211     if (ng <= stop_epsilon_1)
00212       stop = true; // stop = 1, stop by small gradient
00213     else if  (nh <= thrx)
00214       stop = true; // stop = 2, stop by small x-step
00215     else if  (no_feval >= max_evaluations)
00216       stop = true; // stop = 3, number of function evaluations exeeded
00217     else
00218       Delta = std::max(Delta, 2*thrx);
00219     //} found =4
00220     }  // Nonzero h
00221       } // nofail
00222     }  // iteration
00223 
00224     //  Set return values
00225     x_end = x;
00226     finished = true;
00227 
00228     if (trace) { // trim size of trace output
00229       x_values.set_size(no_iter, true);
00230       F_values.set_size(no_iter, true);
00231       ng_values.set_size(no_iter, true);
00232       Delta_values.set_size(no_iter, true);
00233     }
00234 
00235     return true;
00236   }
00237 
00238   bool Newton_Search::search(vec &xn)
00239   {
00240     bool state = search();
00241     xn = get_solution();
00242     return state;
00243   }
00244 
00245   bool Newton_Search::search(const vec &x0, vec &xn)
00246   {
00247     set_start_point(x0);
00248     bool state = search();
00249     xn = get_solution();
00250     return state;
00251   }
00252 
00253   vec Newton_Search::get_solution()
00254   {
00255     it_assert(finished, "Newton_Search: search is not run yet");
00256     return x_end;
00257   }
00258 
00259   double Newton_Search::get_function_value()
00260   {
00261     if (finished)
00262       return F;
00263     else
00264       it_warning("Newton_Search::get_function_value, search has not been run");
00265 
00266     return 0.0;
00267   }
00268 
00269   double Newton_Search::get_stop_1()
00270   {
00271     if (finished)
00272       return ng;
00273     else
00274       it_warning("Newton_Search::get_stop_1, search has not been run");
00275 
00276     return 0.0;
00277   }
00278 
00279   double Newton_Search::get_stop_2()
00280   {
00281     if (finished)
00282       return nh;
00283     else
00284       it_warning("Newton_Search::get_stop_2, search has not been run");
00285 
00286     return 0.0;
00287   }
00288 
00289   int Newton_Search::get_no_iterations()
00290   {
00291     if (finished)
00292       return no_iter;
00293     else
00294       it_warning("Newton_Search::get_no_iterations, search has not been run");
00295 
00296     return 0;
00297   }
00298 
00299   int Newton_Search::get_no_function_evaluations()
00300   {
00301     if (finished)
00302       return no_feval;
00303     else
00304       it_warning("Newton_Search::get_no_function_evaluations, search has not been run");
00305 
00306     return 0;
00307   }
00308 
00309 
00310   void Newton_Search::get_trace(Array<vec> & xvalues, vec &Fvalues, vec &ngvalues, vec &dvalues)
00311   {
00312     if (finished) {
00313       if (trace) { // trim size of trace output
00314   xvalues = x_values;
00315   Fvalues = F_values;
00316   ngvalues = ng_values;
00317   dvalues = Delta_values;
00318       } else
00319   it_warning("Newton_Search::get_trace, trace is not enabled");
00320     } else
00321       it_warning("Newton_Search::get_trace, search has not been run");
00322   }
00323 
00324   //================================== Line_Search =============================================
00325 
00326   Line_Search::Line_Search()
00327   {
00328     method = Soft;
00329 
00330     if (method == Soft) {
00331       stop_rho = 1e-3;
00332       stop_beta = 0.99;
00333     }
00334 
00335     max_iterations = 10;
00336     max_stepsize = 10;
00337 
00338     f = NULL;
00339     df_dx = NULL;
00340     no_feval = 0;
00341     init = false;
00342     finished = false;
00343     trace = false;
00344   }
00345 
00346   void Line_Search::set_function(double(*function)(const vec&))
00347   {
00348     // Add checks to see that function is OK???
00349     f = function;
00350   }
00351 
00352   void Line_Search::set_gradient(vec(*gradient)(const vec&))
00353   {
00354     // Add checks to see that function is OK???
00355     df_dx = gradient;
00356   }
00357 
00358 
00359   void Line_Search::set_stop_values(double rho, double beta)
00360   {
00361     // test input values???
00362     stop_rho = rho;
00363     stop_beta = beta;
00364   }
00365 
00366 
00367   void Line_Search::set_start_point(const vec &x, double F, const vec &g, const vec &h)
00368   {
00369     // check values ???
00370     x_start = x;
00371     F_start = F;
00372     g_start = g;
00373     h_start = h;
00374     n = x.size();
00375 
00376     finished = false;
00377     init = true;
00378   }
00379 
00380   void Line_Search::get_solution(vec &xn, double &Fn, vec &gn)
00381   {
00382     it_assert(finished, "Line_Search: search is not run yet");
00383 
00384     xn = x_end;
00385     Fn = F_end;
00386     gn = g_end;
00387   }
00388 
00389   bool Line_Search::search()
00390   {
00391     it_assert(f != NULL, "Line_Search: Function pointer is not set");
00392     it_assert(df_dx != NULL, "Line_Search: Gradient function pointer is not set");
00393 
00394     it_assert(init, "Line_search: Starting point is not set");
00395 
00396     // Default return values and simple checks
00397     x_end = x_start; F_end = F_start; g_end = g_start;
00398 
00399     // add some checks???
00400     finished = false;
00401 
00402     vec g;
00403 
00404     // return parameters
00405     no_feval = 0;
00406     slope_ratio = 1;
00407 
00408 
00409 
00410     // Check descent condition
00411     double dF0 = dot(h_start,g_end);
00412 
00413     if (trace) { // prepare structures to store trace data
00414       alpha_values.set_size(max_iterations);
00415       F_values.set_size(max_iterations);
00416       dF_values.set_size(max_iterations);
00417       alpha_values(0) = 0;
00418       F_values(0) = F_end;
00419       dF_values(0) = dF0;
00420     }
00421 
00422 
00423     if  (dF0 >= -10*eps*norm(h_start)*norm(g_end)) { // not significantly downhill
00424       if (trace) { // store trace
00425   alpha_values.set_size(1, true);
00426   F_values.set_size(1, true);
00427   dF_values.set_size(1, true);
00428       }
00429       return false;
00430     }
00431 
00432     // Finish initialization
00433     double F0 = F_start, slope0, slopethr;
00434 
00435     if (method == Soft) {
00436       slope0 = stop_rho*dF0; slopethr = stop_beta*dF0;
00437     } else { // exact line search
00438       slope0 = 0;  slopethr = stop_rho*std::abs(dF0);
00439     }
00440 
00441     // Get an initial interval for am
00442     double a = 0, Fa = F_end, dFa = dF0;
00443     bool stop = false;
00444     double b = std::min(1.0, max_stepsize), Fb, dFb;
00445 
00446 
00447     while  (!stop) {
00448       Fb = f(x_start+b*h_start);
00449       g = df_dx(x_start+b*h_start);
00450       // check if these values are OK if not return false???
00451       no_feval++;
00452 
00453       dFb = dot(g,h_start);
00454       if (trace) { // store trace
00455   alpha_values(no_feval) = b;
00456   F_values(no_feval) = Fb;
00457   dF_values(no_feval) = dFb;
00458       }
00459 
00460       if  (Fb < F0 + slope0*b) { // new lower bound
00461   alpha = b;
00462   slope_ratio = dFb/dF0; // info(2);
00463 
00464   if (method == Soft) {
00465     a = b;  Fa = Fb;  dFa = dFb;
00466   }
00467 
00468   x_end = x_start + b*h_start;  F_end = Fb;  g_end = g;
00469 
00470   if  ( (dFb < std::min(slopethr,0.0)) && (no_feval < max_iterations) && (b < max_stepsize) ) {
00471     // Augment right hand end
00472     if (method == Exact) {
00473       a = b;  Fa = Fb;  dFa = dFb;
00474     }
00475     if  (2.5*b >= max_stepsize)
00476       b = max_stepsize;
00477     else
00478       b = 2*b;
00479   } else
00480     stop = true;
00481       } else
00482   stop = true;
00483     } // phase 1: expand interval
00484 
00485 
00486 
00487     if (stop)  // OK so far.  Check stopping criteria
00488       stop = (no_feval >= max_iterations)
00489   || (b >= max_stepsize && dFb < slopethr)
00490   || (a > 0 && dFb >= slopethr);
00491     // Commented by ediap 2006-07-17: redundant check
00492     //  || ( (method == Soft) && (a > 0 & dFb >= slopethr) );  // OK
00493 
00494 
00495     if (stop && trace) {
00496   alpha_values.set_size(no_feval, true);
00497   F_values.set_size(no_feval, true);
00498   dF_values.set_size(no_feval, true);
00499     }
00500 
00501     // Refine interval
00502     while  (!stop) {
00503 
00504       double c, Fc, dFc;
00505 
00506       //c = interpolate(xfd,n);
00507       double C = Fb-Fa - (b-a)*dFa;
00508       if (C >= 5*n*eps*b) {
00509   double A = a - 0.5*dFa*(sqr(b-a)/C);
00510   c = std::min(std::max(a+0.1*(b-a), A), b-0.1*(b-a));  // % Ensure significant resuction
00511       } else
00512   c = (a+b)/2;
00513 
00514       Fc = f(x_start+c*h_start);
00515       g = df_dx(x_start+c*h_start);
00516       dFc = dot(g,h_start);
00517       // check these values???
00518       no_feval++;
00519 
00520       if (trace) { // store trace
00521   alpha_values(no_feval) = c;
00522   F_values(no_feval) = Fc;
00523   dF_values(no_feval) = dFc;
00524       }
00525 
00526       if (method == Soft) {
00527   // soft line method
00528   if  (Fc < F0 + slope0*c) { // new lower bound
00529     alpha = c;
00530     slope_ratio = dFc/dF0;
00531 
00532     x_end = x_start + c*h_start;  F_end = Fc;  g_end = g;
00533     a = c; Fa = Fc; dFa = dFc; // xfd(:,1) = xfd(:,3);
00534     stop = (dFc > slopethr);
00535   } else { // new upper bound
00536     b = c; Fb = Fc; dFb = dFc; // xfd(:,2) = xfd(:,3);
00537   }
00538 
00539       } else { // Exact line search
00540   if  (Fc < F_end) { // better approximant
00541     alpha = c;
00542     slope_ratio = dFc/dF0;
00543     x_end = x_start + c*h_start;  F_end = Fc;  g_end = g;
00544   }
00545   if  (dFc < 0) { // new lower bound
00546     a = c; Fa = Fc; dFa = dFc; // xfd(:,1) = xfd(:,3);
00547   } else { //new upper bound
00548     b = c; Fb = Fc; dFb = dFc; // xfd(:,2) = xfd(:,3);
00549   }
00550   stop = (std::abs(dFc) <= slopethr) | ((b-a) < stop_beta*b);
00551       }
00552 
00553       stop = (stop | (no_feval >= max_iterations));
00554     } // refine
00555 
00556     finished = true;
00557 
00558     if (trace) { // store trace
00559       alpha_values.set_size(no_feval+1, true);
00560       F_values.set_size(no_feval+1, true);
00561       dF_values.set_size(no_feval+1, true);
00562     }
00563 
00564     return true;
00565   }
00566 
00567   bool Line_Search::search(vec &xn, double &Fn, vec &gn)
00568   {
00569     bool state = search();
00570     get_solution(xn, Fn, gn);
00571     return state;
00572   }
00573 
00574   bool Line_Search::search(const vec &x, double F, const vec &g, const vec &h,
00575          vec &xn, double &Fn, vec &gn)
00576   {
00577     set_start_point(x, F, g, h);
00578     bool state = search();
00579     get_solution(xn, Fn, gn);
00580     return state;
00581   }
00582 
00583 
00584   double Line_Search::get_alpha()
00585   {
00586     if (finished)
00587       return alpha;
00588     else
00589       it_warning("Line_Search::get_alpha, search has not been run");
00590 
00591     return 0.0;
00592   }
00593 
00594   double Line_Search::get_slope_ratio()
00595   {
00596     if (finished)
00597       return slope_ratio;
00598     else
00599       it_warning("Line_Search::get_slope_raio, search has not been run");
00600 
00601     return 0.0;
00602   }
00603 
00604   int Line_Search::get_no_function_evaluations()
00605   {
00606     if (finished)
00607       return no_feval;
00608     else
00609       it_warning("Line_Search::get_no_function_evaluations, search has not been run");
00610 
00611     return 0;
00612   }
00613 
00614 
00615   void Line_Search::set_max_iterations(int value)
00616   {
00617     it_assert(value > 0, "Line_Search, max iterations must be > 0");
00618     max_iterations = value;
00619   }
00620 
00621   void Line_Search::set_max_stepsize(double value)
00622   {
00623     it_assert(value > 0, "Line_Search, max stepsize must be > 0");
00624     max_stepsize = value;
00625   }
00626 
00627   void Line_Search::set_method(const Line_Search_Method &search_method)
00628   {
00629     method = search_method;
00630 
00631     if (method == Soft) {
00632       stop_rho = 1e-3;
00633       stop_beta = 0.99;
00634     } else { // exact line search
00635       method = Exact;
00636       stop_rho = 1e-3;
00637       stop_beta = 1e-3;
00638     }
00639   }
00640 
00641 
00642   void Line_Search::get_trace(vec &alphavalues, vec &Fvalues, vec &dFvalues)
00643   {
00644     if (finished) {
00645       if (trace) { // trim size of trace output
00646   alphavalues = alpha_values;
00647   Fvalues = F_values;
00648   dFvalues = dF_values;
00649       } else
00650   it_warning("Line_Search::get_trace, trace is not enabled");
00651     } else
00652       it_warning("Line_Search::get_trace, search has not been run");
00653   }
00654 
00655   // =========================== functions ==============================================
00656 
00657   vec fminunc(double(*function)(const vec&), vec(*gradient)(const vec&), const vec &x0)
00658   {
00659     Newton_Search newton;
00660     newton.set_functions(function, gradient);
00661 
00662     vec xn;
00663     newton.search(x0, xn);
00664 
00665     return xn;
00666   }
00667 
00668 
00669 
00670 } // namespace itpp
SourceForge Logo

Generated on Sun Sep 14 18:52:35 2008 for IT++ by Doxygen 1.5.6