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