bracketing.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_BRACKETING_H
#define pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_BRACKETING_H
 
#include <prc/config.hpp>
#include <prc/core/basic/common.hpp>
#include <prc/core/functors/scalar_product.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/functions/cast.hpp>
#include <prc/core/value/limits.hpp>
#include <prc/core/value/type_traits.hpp>
 
namespace pRC::Optimizer::LineSearch
{
    struct Bracketing
    {
    private:
        static constexpr Bool defaultStrongWolfeConditions()
        {
            return true;
        }
 
        static constexpr Float<> defaultDecreasingStepScale()
        {
            return 0.5;
        }
 
        static constexpr Float<> defaultIncreasingStepScale()
        {
            return 2.1;
        }
 
        static constexpr Size defaultMaxIterations()
        {
            return 5;
        }
 
        static constexpr Float<> defaultC1()
        {
            return 1e-4;
        }
 
        static constexpr Float<> defaultC2()
        {
            return 0.9;
        }
 
    public:
        constexpr Bracketing(
            Bool const strongWolfeConditions = defaultStrongWolfeConditions(),
            Float<> const decreasingStepScale = defaultDecreasingStepScale(),
            Float<> const increasingStepScale = defaultIncreasingStepScale(),
            Size const maxIterations = defaultMaxIterations(),
            Float<> const c1 = defaultC1(), Float<> const c2 = defaultC2())
            : mStrongWolfeConditions(strongWolfeConditions)
            , mDecreasingStepScale(decreasingStepScale)
            , mIncreasingStepScale(increasingStepScale)
            , mMaxIterations(maxIterations)
            , mC1(c1)
            , mC2(c2)
        {
        }
 
        constexpr auto strongWolfeConditions() const
        {
            return mStrongWolfeConditions;
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) decreasingStepScale() const
        {
            return cast<T>(mDecreasingStepScale);
        }
 
        template<class T = Float<>>
        constexpr decltype(auto) increasingStepScale() const
        {
            return cast<T>(mIncreasingStepScale);
        }
 
        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 X, If<IsTensor<X>> = 0,
            class T = typename X::Type::Value, If<IsFloat<T>> = 0, class F,
            If<IsInvocable<F, X const &, X &>> = 0,
            If<IsFloat<ResultOf<F, X const &, X &>>> = 0, class FC,
            If<IsInvocable<FC, X const &>> = 0,
            If<IsConvertible<ResultOf<FC, X const &>, X>> = 0>
        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-Bracketing: Minimum alpha < 0.");
                }
 
                if(alphaMax < zero<T>())
                {
                    Logging::error("LS-Bracketing: Maximum alpha < 0.");
                }
 
                if(alphaMax < alphaMin)
                {
                    Logging::error(
                        "LS-Bracketing: 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 low = alphaMin;
            auto high = alphaMax;
 
            for(Index iteration = 0;; ++iteration,
                      alpha = min(decreasingStepScale<T>() * (low + high),
                          increasingStepScale<T>() * alpha))
            {
                if(iteration > maxIterations())
                {
                    Logging::debug("Line search: Max iterations reached.");
                    break;
                }
 
                x = constraint(x0 + alpha * p);
                f = function(x, g);
 
                if(f > f0 + c1<T>() * alpha * d0)
                {
                    high = alpha;
                    Logging::info("First Wolfe condition not satisfied, trying again");
                    continue;
                }
 
                d = scalarProduct(g, p)();
                if(d < c2<T>() * d0)
                {
                    low = alpha;
                    Logging::info("Second Wolfe condition not satisfied, trying again");
                    continue;
                }
 
                if(strongWolfeConditions())
                {
                    if(abs(d) > abs(c2<T>() * d0))
                    {
                        high = alpha;
                        Logging::info("Strong Wolfe condition not satisfied, trying again");
                        continue;
                    }
                }
 
                break;
            }
 
            return alpha;
        }
 
        template<class X, If<IsTensor<X>> = 0,
            class T = typename X::Type::Value, If<IsFloat<T>> = 0, class F,
            If<IsInvocable<F, X const &, X &>> = 0,
            If<IsFloat<ResultOf<F, X const &, X &>>> = 0>
        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:
        Bool const mStrongWolfeConditions;
        Float<> const mDecreasingStepScale;
        Float<> const mIncreasingStepScale;
        Size const mMaxIterations;
        Float<> const mC1;
        Float<> const mC2;
    };
}
#endif // pRC_ALGORITHMS_OPTIMIZER_LINE_SEARCH_BRACKETING_H