27 #ifndef CDPL_MATH_BFGSMINIMIZER_HPP
28 #define CDPL_MATH_BFGSMINIMIZER_HPP
52 template <typename VA, typename VT = typename MinimizerVariableArrayTraits<VA>::ValueType,
typename FVT = VT>
75 rho(0.01), tau1(9), tau2(0.05), tau3(0.5), order(3), sigma(0.1), func(func), gradFunc(grad_func), status(
SUCCESS) {}
110 for (std::size_t i = 0; max_iter == 0 || i < max_iter; i++) {
111 status =
iterate(fValue, x, g);
116 if (g_norm >=
ValueType(0) && g0Norm <= g_norm)
119 if (delta_f >=
ValueType(0) && deltaF <= delta_f)
136 startF = gradFunc(x, g);
145 multiply(p, -1 / g0Norm);
176 alpha1 = std::min(
ValueType(1), 2 * del / -fp0);
188 updatePosition(alpha, x, f, g);
220 A = -(1 + dgnorm * dgnorm / dxdg) *
B + dgg / dxdg;
238 multiply(p, dir / pNorm);
256 void initFuncEvalCache()
291 return MinimizerVariableArrayTraits<VA>::template norm2<ValueType>(v);
296 return MinimizerVariableArrayTraits<VA>::template dot<ValueType>(v1, v2);
306 return dot(gAlpha, p);
311 if (alpha == xCacheKey)
318 axpy(alpha, p, xAlpha);
325 if (alpha == fCacheKey)
330 fAlpha = func(xAlpha);
338 if (alpha == dfCacheKey)
343 if (alpha != gCacheKey) {
344 fAlpha = gradFunc(xAlpha, gAlpha);
359 if (alpha == fCacheKey && alpha == dfCacheKey) {
365 if (alpha == fCacheKey || alpha == dfCacheKey) {
373 fAlpha = gradFunc(xAlpha, gAlpha);
390 getFDF(alpha, f_alpha, df_alpha);
397 void changeDirection()
440 ValueType r = TypeTraits<ValueType>::abs(TypeTraits<ValueType>::sqrt(disc) / (a * 2));
446 ValueType temp = -(b + sgnb * TypeTraits<ValueType>::sqrt(disc)) / 2;
480 ValueType fl = f0 + zl * (fp0 + zl * (f1 - f0 - fp0));
481 ValueType fh = f0 + zh * (fp0 + zh * (f1 - f0 - fp0));
494 if (z > zl && z < zh) {
495 ValueType f = f0 + z * (fp0 + z * (f1 - f0 - fp0));
521 return c0 + z * (c1 + z * (c2 + z * c3));
541 ValueType eta = 3 * (f1 - f0) - 2 * fp0 - fp1;
542 ValueType xi = fp0 + fp1 - 2 * (f1 - f0);
543 ValueType c0 = f0, c1 = fp0, c2 = eta, c3 = xi;
546 ValueType fmin = cubic(c0, c1, c2, c3, zl);
548 checkExtremum(c0, c1, c2, c3, zh, zmin, fmin);
551 std::size_t n = solveQuadratic(3 * c3, 2 * c2, c1, z0, z1);
554 if (z0 > zl && z0 < zh)
555 checkExtremum(c0, c1, c2, c3, z0, zmin, fmin);
557 if (z1 > zl && z1 < zh)
558 checkExtremum(c0, c1, c2, c3, z1, zmin, fmin);
561 if (z0 > zl && z0 < zh)
562 checkExtremum(c0, c1, c2, c3, z0, zmin, fmin);
585 if (order > 2 && std::isfinite(fpb))
586 z = interpolateCubic(fa, fpa * (b - a), fb, fpb * (b - a), zmin, zmax);
588 z = interpolateQuadratic(fa, fpa * (b - a), fb, zmin, zmax);
597 ValueType falpha, fpalpha, delta, alpha_next;
614 const std::size_t num_brack_iter = 100, section_iters = 100;
617 while (i++ < num_brack_iter) {
618 falpha = getF(alpha);
622 if (falpha > f0 + alpha * rho * fp0 || falpha >= falpha_prev) {
628 fpb = std::numeric_limits<ValueType>::quiet_NaN();
632 fpalpha = getDF(alpha);
636 if (TypeTraits<ValueType>::abs(fpalpha) <= -sigma * fp0) {
651 delta = alpha - alpha_prev;
656 alpha_next = interpolate(alpha_prev, falpha_prev, fpalpha_prev,
657 alpha, falpha, fpalpha, lower, upper);
660 falpha_prev = falpha;
661 fpalpha_prev = fpalpha;
667 while (i++ < section_iters) {
673 alpha = interpolate(a, fa, fpa, b, fb, fpb, lower, upper);
674 falpha = getF(alpha);
676 if ((a - alpha) * fpa <= std::numeric_limits<ValueType>::epsilon()) {
681 if (falpha > f0 + rho * alpha * fp0 || falpha >= fa) {
685 fpb = std::numeric_limits<ValueType>::quiet_NaN();
688 fpalpha = getDF(alpha);
690 if (TypeTraits<ValueType>::abs(fpalpha) <= -sigma * fp0) {
718 const std::size_t order;
749 #endif // CDPL_MATH_BFGSMINIMIZER_HPP