lbfgs.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_OPTIMIZER_LBFGS_H
#define pRC_ALGORITHMS_OPTIMIZER_LBFGS_H
 
#include <prc/algorithms/optimizer/line_search/bracketing.hpp>
#include <prc/core/container/deque.hpp>
#include <prc/core/functors/eval.hpp>
#include <prc/core/functors/scalar_product.hpp>
 
namespace pRC::Optimizer
{
    template<class LS = LineSearch::Bracketing, Size M = 5>
    class LBFGS
    {
    private:
        static constexpr Size defaultMaxIterations()
        {
            return 1000;
        }
 
        template<class G, class T>
        static constexpr auto projectedGradientConverged(G const &g,
            T const &tolerance)
        {
            auto const infNorm = norm<2, 0>(g)();
 
            return infNorm <= tolerance * identity<T>(1e-3);
        }
 
        template<class F, class T>
        static constexpr auto valueConverged(F const &f0, F const &f,
            T const &tolerance)
        {
            auto const scale = max(abs(f0), abs(f), identity<T>());
 
            return delta(f0, f) <= tolerance * scale;
        }
 
        template<class F>
        static constexpr auto valueDiverged(F const &f)
        {
            return !isFinite(f);
        }
 
        template<class F>
        static constexpr auto valueIncreased(F const &f0, F const &f)
        {
            return f > f0;
        }
 
        template<class S, class Y, class R, class H>
        static constexpr auto resetHistory(S &s, Y &y, R &rho, H &H0)
        {
            s.clear();
            y.clear();
            rho.clear();
            H0 = identity();
 
            return;
        }
 
        template<class S, class Y, class R, class H, class G>
        static constexpr auto applyHessianMatrix(S const &s, Y const &y,
            R const &rho, H const &H0, Size const size, G const &g)
        {
            using T = ResultOf<Mul, typename R::Type,
                typename ResultOf<ScalarProduct, typename S::Type, G>::Type>;
            Deque<T, M> alpha;
 
            Tensor q = g;
            for(Index i = 0; i < size; ++i)
            {
                alpha.pushFront(rho.back(i) * scalarProduct(s.back(i), q)());
                q -= alpha.back(i) * y.back(i);
            }
 
            q *= H0;
 
            for(Index i = 0; i < size; ++i)
            {
                auto const beta = rho.front(i) * scalarProduct(y.front(i), q)();
                q += s.front(i) * (alpha.front(i) - beta);
            }
 
            return q;
        }
 
    public:
        constexpr LBFGS(LS const &lineSearch,
            Size const maxIterations = defaultMaxIterations())
            : mLineSearch(lineSearch)
            , mMaxIterations(maxIterations)
        {
        }
 
        constexpr LBFGS(Size const maxIterations = defaultMaxIterations())
            : mMaxIterations(maxIterations)
        {
        }
 
        constexpr auto &lineSearch() const
        {
            return mLineSearch;
        }
 
        constexpr auto maxIterations() const
        {
            return mMaxIterations;
        }
 
        template<class XX, class FF, class FC,
            IsTensorish RX = RemoveReference<XX>,
            IsTensor RXE = RemoveConstReference<ResultOf<Eval, XX>>,
            IsFloat VX = Value<RX>, IsFloat VT = VX>
            requires IsInvocable<FF, RXE const &, RXE &> &&
            IsFloat<ResultOf<FF, RXE const &, RXE &>> &&
            IsInvocable<FC, RXE const &> &&
            IsInvocable<LS, RXE &, ResultOf<FF, RXE const &, RXE &> &, RXE &,
                typename ResultOf<ScalarProduct, XX, XX>::Type &, FF,
                RXE const &>
        inline constexpr auto operator()(XX &&x0, FF &&function, FC &&callback,
            VT const &tolerance = NumericLimits<VT>::tolerance()) const
        {
            using TX = typename RX::Type;
 
            decltype(auto) x =
                copy<!(!IsReference<XX> && !IsConst<RX>)>(eval(x0));
 
            RXE g;
            auto f = function(x, g);
 
            Logging::info("L-BFGS initial f(x) =", f);
 
            if(projectedGradientConverged(g, tolerance))
            {
                return x;
            }
 
            Deque<RXE, M> s;
            Deque<RXE, M> y;
            Deque<TX, M> rho;
            TX H0 = identity();
 
            auto alpha = identity<TX>();
            for(Index iteration = 0;;)
            {
                Tensor p = -applyHessianMatrix(s, y, rho, H0, s.size(), g);
 
                auto d = scalarProduct(p, g)();
                if(d > zero())
                {
                    resetHistory(s, y, rho, H0);
                    p = -g;
                    d = -norm<2, 1>(g)();
                }
 
                auto const f0 = f;
                auto const g0 = g;
                auto const x0 = x;
 
                alpha = lineSearch()(x, f, g, d, function, p, alpha);
 
                callback(x);
 
                if(valueDiverged(f))
                {
                    Logging::info("L-BFGS diverged at f(x) =", f);
                    x = x0;
                    break;
                }
 
                if(valueIncreased(f0, f))
                {
                    Logging::info("L-BFGS made no further progress at f(x) =",
                        f);
                    x = x0;
                    break;
                }
 
                if(++iteration; !(iteration < maxIterations()))
                {
                    Logging::info("L-BFGS max iterations reached at f(x) =", f);
                    break;
                }
 
                if(valueConverged(f0, f, tolerance))
                {
                    Logging::info("L-BFGS converged at f(x) =", f);
                    break;
                }
 
                if(projectedGradientConverged(g, tolerance))
                {
                    Logging::info("L-BFGS converged at f(x) =", f);
                    break;
                }
 
                Logging::info("L-BFGS current f(x) =", f);
 
                s.pushBack(alpha * p);
                y.pushBack(g - g0);
                rho.pushBack(rcp(scalarProduct(y.back(), s.back()))());
                H0 = rcp(rho.back() * norm<2, 1>(y.back())());
            }
 
            callback(x);
 
            if constexpr(IsReference<decltype(x)>)
            {
                return forward<XX>(x0);
            }
            else
            {
                return x;
            }
        }
 
    private:
        LS const mLineSearch;
        Size const mMaxIterations;
    };
}
#endif // pRC_ALGORITHMS_OPTIMIZER_LBFGS_H