Shark Machine Learning Library
BFGS.cpp
Go to the documentation of this file.00001 //=========================================================================== 00073 //=========================================================================== 00074 00075 00076 #include <cmath> 00077 #include <iostream> 00078 #include <float.h> 00079 #include <SharkDefs.h> 00080 #include <ReClaM/BFGS.h> 00081 00082 00083 #define DLINMIN 00084 00085 00086 //=========================================================================== 00116 double BFGS::optimize( 00117 Model& model, 00118 ErrorFunction& errorfunction, 00119 const Array<double>& input, 00120 const Array<double>& target) 00121 { 00122 bfgs(model, errorfunction, input, target, FRET); 00123 return FRET; 00124 } 00125 00126 double BFGS::optimize( 00127 Model& model, 00128 ErrorFunction& errorfunction, 00129 const Array<double>& input, 00130 const Array<double>& target, 00131 double gtol, 00132 int iterMax) 00133 { 00134 unsigned int iter; 00135 bfgs(model, errorfunction, input, target, gtol, iter, FRET, iterMax); 00136 return FRET; 00137 } 00138 00139 00140 void BFGS::bfgs( 00141 Model& model, 00142 ErrorFunction& errorfunction, 00143 const Array<double> &in, 00144 const Array<double> &target, 00145 double gtol, 00146 unsigned& iter, 00147 double& fret, 00148 unsigned iterMax) 00149 { 00150 unsigned i, j, n = BFGS_n; 00151 DEDW.resize(n, false); 00152 // part of Init 00153 00154 if (firstIteration) 00155 { 00156 fret = errorfunction.errorDerivative(model, in, target, DEDW); 00157 BFGS_g0 = DEDW; 00158 firstIteration = false; 00159 } 00160 Array<double> p(n); 00161 for (i = 0;i < n;i++)p(i) = model.getParameter(i); 00162 00163 double d = 0., err, scale; 00164 double test, temp; 00165 00166 Array< double > p1(n); 00167 Array< double > g0(n); 00168 Array< double > g1(n); 00169 Array< double > s(n); 00170 00171 Array< double > H(n, n); 00172 00173 Array< double > gamma(n); 00174 Array< double > delta(n); 00175 00176 iter = 0; 00177 00178 double r = errorfunction.errorDerivative(model, in, target, DEDW); 00179 g0 = DEDW; 00180 00181 for (err = 0., i = 0; i < n; ++i) 00182 err += g0(i) * g0(i); 00183 if (err < gtol*gtol) 00184 { 00185 fret = r; 00186 for (i = 0;i < n;i++)model.setParameter(i, p(i)); 00187 return; 00188 } 00189 00190 for (i = 0; i < n; ++i) 00191 { 00192 for (j = 0; j < n; ++j) 00193 H(i, j) = 0.; 00194 H(i, i) = 1.; 00195 } 00196 00197 for (; ;) 00198 { 00199 for (i = 0; i < n; ++i) 00200 for (s(i) = 0., j = 0; j < n; ++j) 00201 s(i) -= H(i, j) * g0(j); 00202 00203 // 00204 // line search 00205 // 00206 00207 00208 switch (lineSearchType) 00209 { 00210 case Dlinmin: 00211 // dlinmin( p1, s, fret, in, target); 00212 dlinmin(model, errorfunction, DEDW, p1, s, fret, in, target); 00213 break; 00214 case Linmin: 00215 // linmin( p1, s, fret, in, target); 00216 linmin(model, errorfunction, DEDW, p1, s, fret, in, target); 00217 break; 00218 case Cblnsrch: 00219 // cblnsrch(model, errorfunction, CG_fret, dedw, CG_xi, CG_fret, input, target); 00220 cblnsrch(model, errorfunction, fret, DEDW, BFGS_p1, fret, in, target); 00221 // cblnsrch( BFGS_p, errorfuncion.errorDerrivative(input,target,DEDW), BFGS_g0, BFGS_s, BFGS_p1, fret, in, target); 00222 break; 00223 } 00224 00225 for (i = 0;i < n;i++)model.setParameter(i, p1(i)); 00226 errorfunction.errorDerivative(model, in, target, DEDW); 00227 g1 = DEDW; 00228 for (err = 0., i = 0; i < n; ++i) 00229 err += g1(i) * g1(i); 00230 if (++iter >= iterMax || err <= gtol*gtol) 00231 break; 00232 00233 // 00234 // test for convergence 00235 // 00236 test = 0.; 00237 for (i = 0; i < n; i++) 00238 { 00239 temp = fabs(p1(i) - p(i)) / Shark::max((p(i)), 1.); 00240 if (temp > test) 00241 test = temp; 00242 } 00243 if (test < gtol) 00244 { 00245 return; 00246 } 00247 00248 for (d = 0., i = 0; i < n; ++i) 00249 d += (gamma(i) = g1(i) - g0(i)) * (delta(i) = p1(i) - p(i)); 00250 00251 for (i = 0; i < n; ++i) 00252 for (s(i) = 0., j = 0; j < n; ++j) 00253 s(i) += H(i, j) * gamma(j); 00254 00255 if (d < 1e-8) 00256 { 00257 for (i = 0; i < n; ++i) 00258 { 00259 for (j = 0; j < n; ++j) 00260 H(i, j) = 0.; 00261 H(i, i) = 1.; 00262 } 00263 } 00264 else 00265 { 00266 for (scale = 0., i = 0; i < n; ++i) 00267 scale += gamma(i) * s(i); 00268 scale = (scale / d + 1) / d; 00269 00270 for (i = 0; i < n; ++i) 00271 { 00272 g0(i) = g1(i); 00273 p(i) = p1(i); 00274 for (j = 0; j < n; ++j) 00275 H(i, j) += scale * delta(i) * delta(j) 00276 - (s(i) * delta(j) + s(j) * delta(i)) / d; 00277 } 00278 } 00279 } 00280 00281 for (i = 0; i < n; ++i)model.setParameter(i, p1(i)); 00282 00283 } 00284 00285 00286 //=========================================================================== 00321 unsigned i; 00322 void BFGS::init(Model& model) 00323 { 00324 initBfgs(model); 00325 } 00326 00327 00328 void BFGS::initBfgs(Model& model, LineSearch ls, 00329 double ax, double bx, double lambda) 00330 { 00331 firstIteration = true; 00332 BFGS_n = model.getParameterDimension(); 00333 unsigned i, j; 00334 BFGS_p.resize(BFGS_n, false); 00335 for (i = 0;i < BFGS_n;i++)BFGS_p(i) = model.getParameter(i); 00336 lineSearchType = ls; 00337 LS_lambda = lambda; 00338 LS_ax = ax; 00339 LS_bx = bx; 00340 00341 00342 00343 00344 H.resize(BFGS_n, BFGS_n, false); 00345 BFGS_p1.resize(BFGS_n, false); 00346 BFGS_g0.resize(BFGS_n, false); 00347 BFGS_g1.resize(BFGS_n, false); 00348 BFGS_s.resize(BFGS_n, false); 00349 00350 BFGS_gamma.resize(BFGS_n, false); 00351 BFGS_delta.resize(BFGS_n, false); 00352 00353 for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p(i)); 00354 for (i = 0; i < BFGS_n; ++i) 00355 { 00356 for (j = 0; j < BFGS_n; ++j) H(i, j) = 0.; 00357 H(i, i) = 1.; 00358 } 00359 } 00360 00361 void BFGS::reset() 00362 { 00363 unsigned i, j; 00364 BFGS_g0 = DEDW; 00365 BFGS_p1 = 0.; BFGS_g1 = 0.; BFGS_s = 0.; BFGS_gamma = 0.; BFGS_delta = 0.; 00366 for (i = 0; i < BFGS_n; ++i) 00367 { 00368 for (j = 0; j < BFGS_n; ++j) 00369 H(i, j) = 0.; 00370 H(i, i) = 1.; 00371 } 00372 } 00373 00374 00375 //=========================================================================== 00404 void BFGS::bfgs(Model& model, 00405 ErrorFunction& errorfunction, 00406 const Array<double> &in, 00407 const Array<double> &target, 00408 double& fret) 00409 { 00410 unsigned i, j; 00411 DEDW.resize(BFGS_n, false); 00412 // part of Init 00413 00414 if (firstIteration) 00415 { 00416 fret = errorfunction.errorDerivative(model, in, target, DEDW); 00417 BFGS_g0 = DEDW; 00418 firstIteration = false; 00419 } 00420 for (i = 0; i < BFGS_n; ++i) 00421 { 00422 for (BFGS_s(i) = 0., j = 0; j < BFGS_n; ++j) 00423 BFGS_s(i) -= H(i, j) * BFGS_g0(j); 00424 } 00425 00426 // 00427 // line search 00428 // 00429 00430 BFGS_p1 = BFGS_p; 00431 00432 switch (lineSearchType) 00433 { 00434 case Dlinmin: 00435 // dlinmin( p1, s, fret, in, target); 00436 dlinmin(model, errorfunction, DEDW, BFGS_p1, BFGS_s, fret, in, target); 00437 break; 00438 case Linmin: 00439 // linmin( p1, s, fret, in, target); 00440 linmin(model, errorfunction, DEDW, BFGS_p1, BFGS_s, fret, in, target); 00441 break; 00442 case Cblnsrch: 00443 // cblnsrch(model, errorfunction, CG_fret, dedw, CG_xi, CG_fret, input, target); 00444 cblnsrch(model, errorfunction, fret, DEDW, BFGS_p1, fret, in, target); 00445 // cblnsrch( BFGS_p, errorfuncion.errorDerrivative(input,target,DEDW), BFGS_g0, BFGS_s, BFGS_p1, fret, in, target); 00446 break; 00447 } 00448 00449 for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p1(i)); 00450 00451 errorfunction.errorDerivative(model, in, target, DEDW); 00452 BFGS_g1 = DEDW; 00453 00454 for (BFGS_d = 0., i = 0; i < BFGS_n; ++i) 00455 BFGS_d += (BFGS_gamma(i) = BFGS_g1(i) - BFGS_g0(i)) * (BFGS_delta(i) = BFGS_p1(i) - BFGS_p(i)); 00456 00457 for (i = 0; i < BFGS_n; ++i) 00458 for (BFGS_s(i) = 0., j = 0; j < BFGS_n; ++j) 00459 BFGS_s(i) += H(i, j) * BFGS_gamma(j); 00460 00461 if (BFGS_d < 1e-8) 00462 { 00463 for (i = 0; i < BFGS_n; ++i) 00464 { 00465 for (j = 0; j < BFGS_n; ++j) 00466 H(i, j) = 0.; 00467 H(i, i) = 1.; 00468 } 00469 } 00470 else 00471 { 00472 for (BFGS_scale = 0., i = 0; i < BFGS_n; ++i) 00473 BFGS_scale += BFGS_gamma(i) * BFGS_s(i); 00474 BFGS_scale = (BFGS_scale / BFGS_d + 1) / BFGS_d; 00475 00476 for (i = 0; i < BFGS_n; ++i) 00477 { 00478 BFGS_g0(i) = BFGS_g1(i); 00479 BFGS_p(i) = BFGS_p1(i); 00480 for (j = 0; j < BFGS_n; ++j) 00481 H(i, j) += BFGS_scale * BFGS_delta(i) * BFGS_delta(j) 00482 - (BFGS_s(i) * BFGS_delta(j) + BFGS_s(j) * BFGS_delta(i)) / BFGS_d; 00483 } 00484 } 00485 00486 for (i = 0; i < BFGS_n; ++i) 00487 BFGS_p(i) = BFGS_p1(i); 00488 00489 for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p(i)); 00490 } 00491 00492 00494 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) 00495 00496 00497 00498 void BFGS::dlinmin( 00499 Model& model, 00500 ErrorFunction& errorfunction, 00501 Array<double>& dedw, 00502 Array< double >& p, 00503 const Array< double >& xi, 00504 double& fret, 00505 const Array< double >&in, 00506 const Array< double >&target) 00507 { 00508 00509 SIZE_CHECK(p.ndim() == 1 && p.samedim(xi)) 00510 const double GOLD = 1.618034; 00511 const double GLIMIT = 100.; 00512 const double TINY = 1.0e-20; 00513 const unsigned ITMAX = 100; 00514 const double ZEPS = 1.0e-10; 00515 const double TOL = 2.0e-4; 00516 bool ok1, ok2; 00517 unsigned i, iter; 00518 double fa, fb, fc, fp, ax, bx, cx; 00519 double ulim, dum; 00520 double a, b, d = 0, e, fu, fv, fw, fx, q, r, tol1, tol2, u, v, _w, x, xm; 00521 double du, dv, dw, dx, d1, d2, u1, u2, olde; 00522 00523 00524 00525 unsigned n = BFGS_n; 00526 Array< double > xt(n); 00527 Array< double > generalDerivative(n); 00528 00529 ax = LS_ax; 00530 bx = LS_bx; 00531 00532 //=================================================================== 00533 00534 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*ax); 00535 fp = fa = errorfunction.error(model, in, target); 00536 00537 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*bx); 00538 fb = errorfunction.error(model, in, target); 00539 00540 if (fb > fa) 00541 { 00542 00543 dum = ax; 00544 ax = bx; 00545 bx = dum; 00546 dum = fb; 00547 fb = fa; 00548 fa = dum; 00549 } 00550 00551 cx = bx + GOLD * (bx - ax); 00552 00553 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*cx); 00554 fc = errorfunction.error(model, in, target); 00555 00556 while (fb > fc) 00557 { 00558 r = (bx - ax) * (fb - fc); 00559 q = (bx - cx) * (fb - fa); 00560 u = bx - ((bx - cx) * q - (bx - ax) * r) / 00561 (2. * SIGN(Shark::max(fabs(q - r), TINY), q - r)); 00562 ulim = bx + GLIMIT * (cx - bx); 00563 if ((bx - u) *(u - cx) > 0.) 00564 { 00565 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00566 00567 fu = errorfunction.error(model, in, target); 00568 00569 if (fu < fc) 00570 { 00571 ax = bx; 00572 bx = u; 00573 fa = fb; 00574 fb = fu; 00575 break; 00576 } 00577 else if (fu > fb) 00578 { 00579 cx = u; 00580 fc = fu; 00581 break; 00582 } 00583 u = cx + GOLD * (cx - bx); 00584 00585 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00586 fu = errorfunction.error(model, in, target); 00587 00588 } 00589 else if ((cx - u) *(u - ulim) > 0.) 00590 { 00591 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00592 fu = errorfunction.error(model, in, target); 00593 00594 if (fu < fc) 00595 { 00596 bx = cx; 00597 cx = u; 00598 u = cx + GOLD * (cx - bx); 00599 fb = fc; 00600 fc = fu; 00601 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00602 00603 fu = errorfunction.error(model, in, target); 00604 } 00605 } 00606 else if ((u - ulim) *(ulim - cx) >= 0.) 00607 { 00608 u = ulim; 00609 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00610 fu = errorfunction.error(model, in, target); 00611 } 00612 else 00613 { 00614 u = cx + GOLD * (cx - bx); 00615 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00616 fu = errorfunction.error(model, in, target); 00617 } 00618 ax = bx; 00619 bx = cx; 00620 cx = u; 00621 fa = fb; 00622 fb = fc; 00623 fc = fu; 00624 } 00625 //======================================================================= 00626 00627 e = 0.; 00628 if (ax < cx) 00629 { 00630 a = ax; 00631 b = cx; 00632 } 00633 else 00634 { 00635 a = cx; 00636 b = ax; 00637 } 00638 x = _w = v = bx; 00639 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*x); 00640 fw = fv = fx = errorfunction.errorDerivative(model, in, target, dedw); 00641 generalDerivative = dedw; 00642 00643 for (dx = 0., i = 0; i < n; ++i) 00644 dx += xi(i) * generalDerivative(i); 00645 dw = dv = dx; 00646 00647 for (iter = 0; iter < ITMAX; iter++) 00648 { 00649 xm = 0.5 * (a + b); 00650 tol2 = 2. * (tol1 = TOL * fabs(x) + ZEPS); 00651 if (fabs(x - xm) <= (tol2 - 0.5 *(b - a))) 00652 { 00653 break; 00654 } 00655 if (fabs(e) > tol1) 00656 { 00657 d1 = 2. * (b - a); 00658 d2 = d1; 00659 if (dw != dx) 00660 d1 = (_w - x) * dx / (dx - dw); 00661 if (dv != dx) 00662 d2 = (v - x) * dx / (dx - dv); 00663 u1 = x + d1; 00664 u2 = x + d2; 00665 ok1 = (a - u1) * (u1 - b) > 0. && dx * d1 <= 0.; 00666 ok2 = (a - u2) * (u2 - b) > 0. && dx * d2 <= 0.; 00667 olde = e; 00668 e = d; 00669 if (ok1 || ok2) 00670 { 00671 if (ok1 && ok2) 00672 d = (fabs(d1) < fabs(d2) ? d1 : d2); 00673 else if (ok1) 00674 d = d1; 00675 else 00676 d = d2; 00677 if (fabs(d) <= fabs(0.5 * olde)) 00678 { 00679 u = x + d; 00680 if (u - a < tol2 || b - u < tol2) 00681 d = SIGN(tol1, xm - x); 00682 } 00683 else 00684 { 00685 d = 0.5 * (e = (dx >= 0. ? a - x : b - x)); 00686 } 00687 } 00688 else 00689 { 00690 d = 0.5 * (e = (dx >= 0. ? a - x : b - x)); 00691 } 00692 } 00693 else 00694 { 00695 d = 0.5 * (e = (dx >= 0. ? a - x : b - x)); 00696 } 00697 if (fabs(d) >= tol1) 00698 { 00699 u = x + d; 00700 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00701 fu = errorfunction.errorDerivative(model, in, target, dedw); 00702 generalDerivative = dedw; 00703 } 00704 else 00705 { 00706 u = x + SIGN(tol1, d); 00707 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00708 fu = errorfunction.errorDerivative(model, in, target, dedw); 00709 generalDerivative = dedw; 00710 if (fu > fx) 00711 break; 00712 } 00713 for (du = 0., i = 0; i < n; ++i) 00714 du += xi(i) * generalDerivative(i); 00715 if (fu <= fx) 00716 { 00717 if (u >= x) 00718 a = x; 00719 else 00720 b = x; 00721 v = _w; 00722 _w = x; 00723 x = u; 00724 fv = fw; 00725 fw = fx; 00726 fx = fu; 00727 dv = dw; 00728 dw = dx; 00729 dx = du; 00730 } 00731 else 00732 { 00733 if (u < x) 00734 a = u; 00735 else 00736 b = u; 00737 if (fu <= fw || _w == x) 00738 { 00739 v = _w; 00740 _w = u; 00741 fv = fw; 00742 fw = fu; 00743 dv = dw; 00744 dw = du; 00745 } 00746 else if (fu < fv || v == x || v == _w) 00747 { 00748 v = u; 00749 fv = fu; 00750 dv = du; 00751 } 00752 } 00753 } 00754 if (iter >= ITMAX) 00755 throw SHARKEXCEPTION("BFGS: too many iterations in routine dbrent"); 00756 if (fx < fp) 00757 { 00758 fret = fx; 00759 for (i = 0; i < n; ++i) 00760 p(i) += xi(i) * x; 00761 } 00762 else 00763 fret = fp; 00764 } 00765 00766 // linmin.cpp 00767 00768 void BFGS::linmin( 00769 Model& model, 00770 ErrorFunction& errorfunction, 00771 Array<double>& dedw, 00772 Array< double >& p, 00773 const Array< double >& xi, 00774 double& fret, 00775 const Array< double >&in, 00776 const Array< double >&target) 00777 { 00778 00779 SIZE_CHECK(p.ndim() == 1 && p.samedim(xi)) 00780 const double GOLD = 1.618034; 00781 const double GLIMIT = 100.; 00782 const double TINY = 1.0e-20; 00783 const unsigned ITMAX = 100; 00784 const double CGOLD = 0.3819660; 00785 const double ZEPS = 1.0e-10; 00786 const double TOL = 2.0e-4; 00787 00788 unsigned i, iter; 00789 double fa, fb, fc, ax, bx, cx; 00790 double ulim, dum; 00791 double a, b, d = 0., e, etemp, fu, fv, fw, fx, o, q, r, tol1, tol2, u, v, _w, x, xm; 00792 unsigned n = BFGS_n; 00793 Array< double > xt(n); 00794 00795 ax = LS_ax; 00796 bx = LS_bx; 00797 00798 //=================================================================== 00799 00800 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*ax); 00801 double initial = fa = errorfunction.error(model, in, target); 00802 00803 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*bx); 00804 fb = errorfunction.error(model, in, target); 00805 00806 if (fb > fa) 00807 { 00808 dum = ax; 00809 ax = bx; 00810 bx = dum; 00811 dum = fb; 00812 fb = fa; 00813 fa = dum; 00814 } 00815 cx = bx + GOLD * (bx - ax); 00816 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*cx); 00817 fc = errorfunction.error(model, in, target); 00818 while (fb > fc) 00819 { 00820 r = (bx - ax) * (fb - fc); 00821 q = (bx - cx) * (fb - fa); 00822 u = bx - ((bx - cx) * q - (bx - ax) * r) / 00823 (2. * SIGN(Shark::max(fabs(q - r), TINY), q - r)); 00824 ulim = bx + GLIMIT * (cx - bx); 00825 if ((bx - u) *(u - cx) > 0.) 00826 { 00827 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00828 fu = errorfunction.error(model, in, target); 00829 if (fu < fc) 00830 { 00831 ax = bx; 00832 bx = u; 00833 fa = fb; 00834 fb = fu; 00835 break; 00836 } 00837 else if (fu > fb) 00838 { 00839 cx = u; 00840 fc = fu; 00841 break; 00842 } 00843 u = cx + GOLD * (cx - bx); 00844 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00845 fu = errorfunction.error(model, in, target); 00846 } 00847 else if ((cx - u) *(u - ulim) > 0.) 00848 { 00849 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00850 fu = errorfunction.error(model, in, target); 00851 if (fu < fc) 00852 { 00853 bx = cx; 00854 cx = u; 00855 u = cx + GOLD * (cx - bx); 00856 fb = fc; 00857 fc = fu; 00858 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00859 fu = errorfunction.error(model, in, target); 00860 } 00861 } 00862 else if ((u - ulim) *(ulim - cx) >= 0.) 00863 { 00864 u = ulim; 00865 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00866 fu = errorfunction.error(model, in, target); 00867 } 00868 else 00869 { 00870 u = cx + GOLD * (cx - bx); 00871 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00872 fu = errorfunction.error(model, in, target); 00873 } 00874 ax = bx; 00875 bx = cx; 00876 cx = u; 00877 fa = fb; 00878 fb = fc; 00879 fc = fu; 00880 } 00881 //======================================================================= 00882 00883 e = 0.; 00884 if (ax < cx) 00885 { 00886 a = ax; 00887 b = cx; 00888 } 00889 else 00890 { 00891 a = cx; 00892 b = ax; 00893 } 00894 x = _w = v = bx; 00895 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*x); 00896 fw = fv = fx = errorfunction.error(model, in, target); 00897 00898 for (iter = 0; iter < ITMAX; iter++) 00899 { 00900 xm = 0.5 * (a + b); 00901 tol2 = 2. * (tol1 = TOL * fabs(x) + ZEPS); 00902 if (fabs(x - xm) <= (tol2 - 0.5 *(b - a))) 00903 { 00904 break; 00905 } 00906 if (fabs(e) > tol1) 00907 { 00908 r = (x - _w) * (fx - fv); 00909 q = (x - v) * (fx - fw); 00910 o = (x - v) * q - (x - _w) * r; 00911 q = 2. * (q - r); 00912 if (q > 0.) 00913 o = -o; 00914 q = fabs(q); 00915 etemp = e; 00916 e = d; 00917 if (fabs(o) >= fabs(0.5 * q * etemp) || 00918 o <= q *(a - x) || 00919 o >= q *(b - x)) 00920 d = CGOLD * (e = (x >= xm ? a - x : b - x)); 00921 else 00922 { 00923 d = o / q; 00924 u = x + d; 00925 if (u - a < tol2 || b - u < tol2) 00926 d = SIGN(tol1, xm - x); 00927 } 00928 } 00929 else 00930 d = CGOLD * (e = (x >= xm ? a - x : b - x)); 00931 u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d)); 00932 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u); 00933 fu = errorfunction.error(model, in, target); 00934 if (fu <= fx) 00935 { 00936 if (u >= x) 00937 a = x; 00938 else 00939 b = x; 00940 v = _w; 00941 _w = x; 00942 x = u; 00943 fv = fw; 00944 fw = fx; 00945 fx = fu; 00946 } 00947 else 00948 { 00949 if (u < x) 00950 a = u; 00951 else 00952 b = u; 00953 if (fu <= fw || _w == x) 00954 { 00955 v = _w; 00956 _w = u; 00957 fv = fw; 00958 fw = fu; 00959 } 00960 else if (fu <= fv || v == x || v == _w) 00961 { 00962 v = u; 00963 fv = fu; 00964 } 00965 } 00966 } 00967 00968 if (iter >= ITMAX) 00969 throw SHARKEXCEPTION("BFGS: too many iterations in brent"); 00970 fret = fx; 00971 //======================================================================= 00972 00973 if (fx < initial) 00974 { 00975 fret = fx; 00976 for (i = 0; i < n; ++i) 00977 p(i) += xi(i) * x; 00978 } 00979 else 00980 { 00981 fret = initial; 00982 } 00983 } 00984 //====================================================================== 00985 // Definition of the cubic line search class 00986 // Armijo and Goldstein's line search algorithm 00987 00988 // author: Doug Hart, Adapted to the COOOL library by Wenceslau Gouveia 00989 00990 // Modified to fit into new classes. H. Lydia Deng, 02/21/94, 03/15/94 00991 00992 //======================================================================== 00993 00994 00995 /* 00996 00997 * cubic line search 00998 00999 * 01000 01001 * The parameter lambda controls the accuraccy of the line search. 01002 01003 * lambda = .25 is a good choice. 01004 01005 */ 01006 01007 void BFGS::cblnsrch( 01008 Model& model, 01009 ErrorFunction& errorfunction, 01010 double fold, 01011 Array< double >& g, 01012 Array< double >& p, 01013 double& f, 01014 const Array< double >&in, 01015 const Array< double >&target) 01016 { 01017 // 01018 // the maximum of performed iterations should be choosen small enough 01019 // to ensure that neither (3^iterMax) nor (0.1^(2*iterMax)) causes a 01020 // numerical overflow or underflow. 01021 // 01022 const unsigned iterMax = 20; 01023 unsigned iterNum; 01024 unsigned i, n; 01025 bool tst = false; 01026 double slope; 01027 double alpha, alpha2; 01028 double alpha_prev, alpha_prev2; 01029 double alpha_tmp = 0; 01030 double f1, f2, fprev = 0; 01031 double a, b, c; 01032 double cm11, cm12, cm21, cm22; 01033 double disc; 01034 double lambda = LS_lambda; 01035 01036 n = BFGS_n; 01037 /* 01038 01039 * dot product of search direction and gradient 01040 01041 */ 01042 01043 for (slope = 0., i = 0; i < n; i++) 01044 slope += g(i) * p(i); 01045 01046 iterNum = 0; /* iteration counter */ 01047 alpha = 1.; /* updating step */ 01048 01049 /* 01050 * updating 01051 */ 01052 01053 for (i = 0; i < n; i++) 01054 model.setParameter(i, model.getParameter(i) + alpha * p(i)); 01055 f = errorfunction.error(model, in, target); 01056 iterNum++; 01057 /* 01058 * Implementing Goldstein's test for alpha too small 01059 */ 01060 while (f < fold + (1. - lambda)*alpha*slope && iterNum < iterMax) 01061 { 01062 alpha *= 3; 01063 for (i = 0; i < n; i++) 01064 model.setParameter(i, model.getParameter(i) + alpha * p(i)); 01065 f = errorfunction.error(model, in, target); 01066 iterNum++; 01067 } 01068 if (iterNum >= iterMax) 01069 throw SHARKEXCEPTION("BFGS: alpha overflowed!"); 01070 01071 /* 01072 * Armijo's test for alpha too large 01073 */ 01074 01075 alpha_prev = alpha; /* H.L. Deng, 6/13/95 */ 01076 while (f > fold + lambda*alpha*slope && iterNum < iterMax) 01077 { 01078 alpha2 = alpha * alpha; 01079 f1 = f - fold - slope * alpha; 01080 if (tst == false) 01081 { 01082 alpha_tmp = -slope * alpha2 / (f1 * 2.); /* tentative alpha */ 01083 tst = true; 01084 } 01085 else 01086 { 01087 alpha_prev2 = alpha_prev * alpha_prev; 01088 f2 = fprev - fold - alpha_prev * slope; 01089 c = 1. / (alpha - alpha_prev); 01090 cm11 = 1. / alpha2; 01091 cm12 = -1. / alpha_prev2; 01092 cm21 = -alpha_prev / alpha2; 01093 cm22 = alpha / alpha_prev2; 01094 a = c * (cm11 * f1 + cm12 * f2); 01095 b = c * (cm21 * f1 + cm22 * f2); 01096 disc = b * b - 3. * a * slope; 01097 01098 if ((fabs(a) > FLT_MIN) && (disc > FLT_MIN)) 01099 alpha_tmp = (-b + sqrt(disc)) / (3. * a); 01100 else 01101 alpha_tmp = slope * alpha2 / (2. * f1); 01102 if (alpha_tmp >= .5 * alpha) 01103 alpha_tmp = .5 * alpha; 01104 } 01105 alpha_prev = alpha; 01106 fprev = f; 01107 if (alpha_tmp < .1 * alpha) 01108 alpha *= .1; 01109 else 01110 alpha = alpha_tmp; 01111 for (i = 0; i < n; i++) 01112 model.setParameter(i, model.getParameter(i) + alpha * p(i)); 01113 f = errorfunction.error(model, in, target); 01114 iterNum++; 01115 } 01116 if (iterNum >= iterMax) 01117 throw SHARKEXCEPTION("BFGS: alpha underflowed!"); 01118 } 01119