more_thuente.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
// References:
//  Authors:    Jorge J. More, David J. Thuente
//  Title:      Line Search Algorithms with Guaranteed Sufficient Decrease
//  Year:       1994
//  URL:        https://doi.org/10.1145/192115.192132
 
#ifndef pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_MORE_THUENTE_H
#define pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_MORE_THUENTE_H
 
#include <prc/core/functors/scalar_product.hpp>
#include <prc/core/tensor/declarations.hpp>
#include <prc/core/value/float.hpp>
#include <prc/core/value/functions/delta.hpp>
 
namespace pRC::Optimizer::LineSearch
{
    struct MoreThuente
    {
    private:
        static constexpr Size defaultMaxIterations()
        {
            return 20;
        }
 
        static constexpr Float<> defaultC1()
        {
            return 1e-4;
        }
 
        static constexpr Float<> defaultC2()
        {
            return 0.9;
        }
 
        static constexpr Float<> defaultTrapLower()
        {
            return 1.1;
        }
 
        static constexpr Float<> defaultTrapUpper()
        {
            return 4.0;
        }
 
        static constexpr Float<> defaultDelta()
        {
            return 2.0 / 3.0;
        }
 
    public:
        constexpr MoreThuente(Size const maxIterations = defaultMaxIterations(),
            Float<> const c1 = defaultC1(), Float<> const c2 = defaultC2(),
            Float<> const trapLower = defaultTrapLower(),
            Float<> const trapUpper = defaultTrapUpper(),
            Float<> const delta = defaultDelta())
            : mMaxIterations(maxIterations)
            , mC1(c1)
            , mC2(c2)
            , mTrapLower(trapLower)
            , mTrapUpper(trapUpper)
            , mDelta(delta)
        {
        }
 
        constexpr auto maxIterations() const
        {
            return mMaxIterations;
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) c1() const
        {
            return cast<T>(mC1);
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) c2() const
        {
            return cast<T>(mC2);
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) trapLower() const
        {
            return cast<T>(mTrapLower);
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) trapUpper() const
        {
            return cast<T>(mTrapUpper);
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) delta() const
        {
            return cast<T>(mDelta);
        }
 
        template<IsTensor X, IsFloat T = Value<X>, class F, class FC>
            requires IsInvocable<F, X const &, X &> &&
            IsFloat<ResultOf<F, X const &, X &>> &&
            IsInvocable<FC, X const &> &&
            IsConvertible<ResultOf<FC, X const &>, X>
        constexpr auto operator()(X &x, ResultOf<F, X const &, X &> &f, X &g,
            typename ResultOf<ScalarProduct, X, X>::Type &d, F &&function,
            FC &&constraint, X const &p, T alpha = identity<T>(),
            T const alphaMin = zero<T>(),
            T const alphaMax = identity<T>(NumericLimits<T>::max())) const
        {
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(alphaMin < zero<T>())
                {
                    Logging::error("LS-MoreThuente: Minimum alpha < 0.");
                }
 
                if(alphaMax < zero<T>())
                {
                    Logging::error("LS-MoreThuente: Maximum alpha < 0.");
                }
 
                if(alphaMax < alphaMin)
                {
                    Logging::error(
                        "LS-MoreThuente: Minimum alpha > Maximum alpha.");
                }
 
                if(alpha <= alphaMin || alpha > alphaMax)
                {
                    Logging::error("Initial alpha not in range (min, max]");
                }
            }
 
            Tensor const x0 = x;
            auto const f0 = f;
            auto const d0 = d;
 
            auto const curvatureTest = c1<T>() * d0;
 
            T alphaLower = zero();
            auto fLower = f0;
            auto dLower = d0;
 
            T alphaUpper = zero();
            auto fUpper = f0;
            auto dUpper = d0;
 
            auto alphaLowerBound = zero<T>();
            auto alphaUpperBound = alpha + trapUpper<T>() * alpha;
 
            auto nextWidth = alphaMax - alphaMin;
            auto width = identity<T>(2) * nextWidth;
 
            Bool firstStage = true;
            Bool bracketed = false;
            for(Index iteration = 0;; ++iteration)
            {
                if(iteration > maxIterations())
                {
                    Logging::debug("Line search: Max iterations reached.");
                    break;
                }
 
                if(bracketed)
                {
                    if(alpha <= alphaLowerBound || alpha >= alphaUpperBound)
                    {
                        Logging::debug(
                            "Line search: Rounding errors prevent progress.");
                        break;
                    }
 
                    if(isZero(alphaUpperBound - alphaLowerBound))
                    {
                        Logging::debug("Line search: Tolerance is satisfied.");
                        break;
                    }
                }
 
                x = constraint(x0 + alpha * p);
                f = function(x, g);
                d = scalarProduct(p, g)();
 
                auto const sufficientDecreaseTest = f0 + alpha * curvatureTest;
 
                if(alpha == alphaMax && f <= sufficientDecreaseTest &&
                    d <= curvatureTest)
                {
                    Logging::debug("Line search: Max alpha reached.");
                    break;
                }
 
                if(alpha == alphaMin &&
                    (f > sufficientDecreaseTest || d >= curvatureTest))
                {
                    Logging::debug("Line search: Min alpha reached.");
                    break;
                }
 
                if(f <= sufficientDecreaseTest && abs(d) <= abs(c2<T>() * d0))
                {
                    break;
                }
 
                if(firstStage && f <= sufficientDecreaseTest && d >= zero())
                {
                    firstStage = false;
                }
 
                if(firstStage && f <= fLower && f > sufficientDecreaseTest)
                {
                    auto const fMod = f - alpha * curvatureTest;
                    auto fLowerMod = fLower - alphaLower * curvatureTest;
                    auto fUpperMod = fUpper - alphaUpper * curvatureTest;
                    auto const dMod = d - curvatureTest;
                    auto dLowerMod = dLower - curvatureTest;
                    auto dUpperMod = dUpper - curvatureTest;
 
                    alpha = computeAlpha(alphaLower, fLowerMod, dLowerMod,
                        alphaUpper, fUpperMod, dUpperMod, alpha, fMod, dMod,
                        bracketed, alphaLowerBound, alphaUpperBound);
 
                    fLower = fLowerMod + alphaLower * curvatureTest;
                    fUpper = fUpperMod + alphaUpper * curvatureTest;
                    dLower = dLowerMod + curvatureTest;
                    dUpper = dUpperMod + curvatureTest;
                }
                else
                {
                    alpha = computeAlpha(alphaLower, fLower, dLower, alphaUpper,
                        fUpper, dUpper, alpha, f, d, bracketed, alphaLowerBound,
                        alphaUpperBound);
                }
 
                if(bracketed)
                {
                    auto const alphaDelta = pRC::delta(alphaLower, alphaUpper);
 
                    if(alphaDelta >= delta<T>() * width)
                    {
                        alpha = mean(alphaLower, alphaUpper);
                    }
 
                    width = nextWidth;
                    nextWidth = alphaDelta;
                }
 
                if(bracketed)
                {
                    alphaLowerBound = min(alphaLower, alphaUpper);
                    alphaUpperBound = max(alphaLower, alphaUpper);
                }
                else
                {
                    alphaLowerBound =
                        alpha + trapLower<T>() * (alpha - alphaLower);
                    alphaUpperBound =
                        alpha + trapUpper<T>() * (alpha - alphaLower);
                }
 
                alpha = max(alpha, alphaMin);
                alpha = min(alpha, alphaMax);
 
                if(bracketed &&
                    ((alpha <= alphaLowerBound || alpha >= alphaUpperBound) ||
                        isZero(alphaUpperBound - alphaLowerBound)))
                {
                    alpha = alphaLower;
                }
            }
 
            return alpha;
        }
 
        template<IsTensor X, IsFloat T = Value<X>, class F>
            requires IsInvocable<F, X const &, X &> &&
            IsFloat<ResultOf<F, X const &, X &>>
        constexpr auto operator()(X &x, ResultOf<F, X const &, X &> &f, X &g,
            typename ResultOf<ScalarProduct, X, X>::Type &d, F &&function,
            X const &p, T alpha = identity<T>(), T const alphaMin = zero<T>(),
            T const alphaMax = identity<T>(NumericLimits<T>::max())) const
        {
            return operator()(
                x, f, g, d, forward<F>(function),
                [](auto &&x) -> decltype(auto)
                {
                    return forward<decltype(x)>(x);
                },
                p, alpha, alphaMin, alphaMax);
        }
 
    private:
        template<typename T>
        static constexpr auto secantMinimizer(T const &x, T const &dX,
            T const &y, T const &dY)
        {
            return y + dY / (dY - dX) * (x - y);
        }
 
        template<typename T>
        static constexpr auto quadraticMinimizer(T const &x, T const &fX,
            T const &dX, T const &y, T const &fY)
        {
            return x +
                dX / ((fX - fY) / (y - x) + dX) / identity<T>(2) * (y - x);
        }
 
        template<typename T>
        static constexpr auto cubicMinimizer(T const &x, T const &fX,
            T const &dX, T const &y, T const &fY, T const &dY)
        {
            auto const theta = identity<T>(3) * (fX - fY) / (y - x) + dX + dY;
 
            auto const s = max(abs(theta), abs(dX), abs(dY));
 
#ifdef __FAST_MATH__
            auto gamma = sqrt(square(theta) - dX * dY);
#else
            auto gamma = s * sqrt(square(theta / s) - (dX / s) * (dY / s));
#endif // __FAST_MATH__
 
            if(y < x)
            {
                gamma = -gamma;
            }
 
            auto const p = (gamma - dX) + theta;
            auto const q = (gamma - dX) + gamma + dY;
            auto const r = p / q;
 
            return x + r * (y - x);
        }
 
        template<typename T>
        static constexpr auto cubicMinimizer(T const &x, T const &fX,
            T const &dX, T const &y, T const &fY, T const &dY,
            T const &lowerBound, T const &upperBound)
        {
            auto const theta = identity<T>(3) * (fX - fY) / (y - x) + dX + dY;
 
            auto const s = max(abs(theta), abs(dX), abs(dY));
 
#ifdef __FAST_MATH__
            auto gamma = sqrt(max(zero<T>(), square(theta) - dX * dY));
#else
            auto gamma = s *
                sqrt(max(zero<T>(), square(theta / s) - (dX / s) * (dY / s)));
#endif // __FAST_MATH__
 
            if(y < x)
            {
                gamma = -gamma;
            }
 
            auto const p = (gamma - dX) + theta;
            auto const q = (gamma - dX) + gamma + dY;
            auto const r = p / q;
 
            if((r < zero()) && (gamma != zero()))
            {
                return x + r * (y - x);
            }
            else if(y < x)
            {
                return upperBound;
            }
            else
            {
                return lowerBound;
            }
        }
 
    private:
        template<typename T>
        constexpr auto computeAlpha(T &alphaLower, T &fLower, T &dLower,
            T &alphaUpper, T &fUpper, T &dUpper, T const &alpha, T const &f,
            T const &d, Bool &bracketed, T const &alphaLowerBound,
            T const &alphaUpperBound) const
        {
            if(f > fLower)
            {
                auto const alphaC =
                    cubicMinimizer(alphaLower, fLower, dLower, alpha, f, d);
 
                auto const alphaQ =
                    quadraticMinimizer(alphaLower, fLower, dLower, alpha, f);
 
                bracketed = true;
                alphaUpper = alpha;
                fUpper = f;
                dUpper = d;
 
                if(pRC::delta(alphaC, alphaLower) <
                    pRC::delta(alphaQ, alphaLower))
                {
                    return alphaC;
                }
                else
                {
                    return mean(alphaQ, alphaC);
                }
            }
 
            RemoveConstReference<decltype(alpha)> nextAlpha;
            if(d * dLower < zero())
            {
                auto const alphaC =
                    cubicMinimizer(alpha, f, d, alphaLower, fLower, dLower);
 
                auto const alphaS =
                    secantMinimizer(alphaLower, dLower, alpha, d);
 
                if(pRC::delta(alphaC, alpha) > pRC::delta(alphaS, alpha))
                {
                    nextAlpha = alphaC;
                }
                else
                {
                    nextAlpha = alphaS;
                }
 
                bracketed = true;
                alphaUpper = alphaLower;
                fUpper = fLower;
                dUpper = dLower;
            }
            else if(abs(d) < abs(dLower))
            {
                auto const alphaC = cubicMinimizer(alpha, f, d, alphaLower,
                    fLower, dLower, alphaLowerBound, alphaUpperBound);
 
                auto const alphaS =
                    secantMinimizer(alphaLower, dLower, alpha, d);
 
                if(bracketed)
                {
                    auto const trap = alpha + delta<T>() * (alphaUpper - alpha);
                    if(pRC::delta(alphaC, alpha) < pRC::delta(alphaS, alpha))
                    {
                        if(alpha > alphaLower)
                        {
                            nextAlpha = min(trap, alphaC);
                        }
                        else
                        {
                            nextAlpha = max(trap, alphaC);
                        }
                    }
                    else
                    {
                        if(alpha > alphaLower)
                        {
                            nextAlpha = min(trap, alphaS);
                        }
                        else
                        {
                            nextAlpha = max(trap, alphaS);
                        }
                    }
                }
                else
                {
                    if(pRC::delta(alphaC, alpha) > pRC::delta(alphaS, alpha))
                    {
                        nextAlpha = alphaC;
                    }
                    else
                    {
                        nextAlpha = alphaS;
                    }
 
                    nextAlpha = min(alphaUpperBound, nextAlpha);
                    nextAlpha = max(alphaLowerBound, nextAlpha);
                }
            }
            else
            {
                if(bracketed)
                {
                    nextAlpha =
                        cubicMinimizer(alpha, f, d, alphaUpper, fUpper, dUpper);
                }
                else if(alpha > alphaLower)
                {
                    nextAlpha = alphaUpperBound;
                }
                else
                {
                    nextAlpha = alphaLowerBound;
                }
            }
 
            alphaLower = alpha;
            fLower = f;
            dLower = d;
 
            return nextAlpha;
        }
 
    private:
        Size const mMaxIterations;
        Float<> const mC1;
        Float<> const mC2;
        Float<> const mTrapLower;
        Float<> const mTrapUpper;
        Float<> const mDelta;
    };
}
#endif // pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_MORE_THUENTE_H