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
Generated on Fri Jun 8 00:27:14 2007 for IT++ by Doxygen 1.5.2