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
Generated on Sun Sep 14 18:57:05 2008 for IT++ by Doxygen 1.5.6