bfgs.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_OPTIMIZER_BFGS_H
#define pRC_ALGORITHMS_OPTIMIZER_BFGS_H
 
#include <prc/config.hpp>
#include <prc/algorithms/optimizer/line_search/bracketing.hpp>
#include <prc/core/basic/functions/copy.hpp>
#include <prc/core/functors/eval.hpp>
#include <prc/core/tensor/functions/norm.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/functions/delta.hpp>
#include <prc/core/value/limits.hpp>
 
namespace pRC::Optimizer
{
    template<class LS = LineSearch::Bracketing>
    class BFGS
    {
    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;
        }
 
    public:
        constexpr BFGS(LS const &lineSearch,
            Size const maxIterations = defaultMaxIterations())
            : mLineSearch(lineSearch)
            , mMaxIterations(maxIterations)
        {
        }
 
        constexpr BFGS(Size const maxIterations = defaultMaxIterations())
            : mMaxIterations(maxIterations)
        {
        }
 
        constexpr auto &lineSearch() const
        {
            return mLineSearch;
        }
 
        constexpr auto maxIterations() const
        {
            return mMaxIterations;
        }
 
        template<class XX, class RX = RemoveReference<XX>,
            class TX = typename RX::Type, class VX = typename TX::Value,
            If<IsTensorish<RX>> = 0,
            class RXE = RemoveConstReference<ResultOf<Eval, XX>>, class FF,
            If<IsInvocable<FF, RXE const &, RXE &>> = 0,
            If<IsFloat<ResultOf<FF, RXE const &, RXE &>>> = 0, class FC,
            If<IsInvocable<FC, RXE &>> = 0, class VT = VX,
            If<All<IsFloat<VX>, IsFloat<VT>>> = 0,
            If<IsInvocable<LS, RXE &, ResultOf<FF, RXE const &, RXE &> &, RXE &,
                VX &, FF, RXE const &>> = 0>
        inline constexpr auto operator()(XX &&x0, FF &&function, FC &&callback,
            VT const &tolerance = NumericLimits<VT>::tolerance()) const
        {
            decltype(auto) x =
                copy<!(!IsReference<XX>() && !IsConst<RX>())>(eval(x0));
 
            RXE g;
            auto f = function(x, g);
 
            Logging::debug("BFGS initial f(x) =", f);
 
            if(projectedGradientConverged(g, tolerance))
            {
                return x;
            }
 
            Tensor H = expand(makeSeries<Index, typename RXE::Dimension{}>(),
                [&](auto const... seq)
                {
                    return identity<
                        Tensor<TX, RXE::size(seq)..., RXE::size(seq)...>>();
                });
 
            auto alpha = identity<TX>();
            for(Index iteration = 0;;)
            {
                Tensor p = -H * g;
 
                auto d = scalarProduct(p, g)();
                if(d > zero())
                {
                    H = identity();
                    p = -g;
                    d = -norm<2, 1>(g)();
                }
 
                auto const f0 = f;
                auto const g0 = g;
 
                alpha = lineSearch()(x, f, g, d, function, p, alpha);
 
                callback(x);
 
                if(++iteration; !(iteration < maxIterations()))
                {
                    Logging::info("BFGS max iterations reached at f(x) =", f);
                    break;
                }
 
                if(valueConverged(f0, f, tolerance))
                {
                    Logging::info("BFGS converged at f(x) =", f);
                    break;
                }
 
                if(projectedGradientConverged(g, tolerance))
                {
                    Logging::info("BFGS converged at f(x) =", f);
                    break;
                }
 
                Logging::debug("BFGS current f(x) =", f);
 
                Tensor const s = alpha * p;
                Tensor const y = g - g0;
                Tensor const rho = rcp(scalarProduct(y, s));
 
                Tensor const V =
                    identity<decltype(H)>() - rho * tensorProduct(y, s);
                H = eval(transpose(V) * H) * V + rho * tensorProduct(s, s);
            }
 
            if constexpr(IsReference<decltype(x)>())
            {
                return forward<XX>(x0);
            }
            else
            {
                return x;
            }
        }
 
    private:
        LS const mLineSearch;
        Size const mMaxIterations;
    };
}
#endif // pRC_ALGORITHMS_OPTIMIZER_BFGS_H