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/algorithms/optimizer/line_search/bracketing.hpp>
#include <prc/core/functors/eval.hpp>
#include <prc/core/functors/scalar_product.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;
        }
 
        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;
        }
 
    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 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::debug("BFGS initial f(x) =", f);
 
            if(projectedGradientConverged(g, tolerance))
            {
                return x;
            }
 
            Tensor H = expand(makeSeries<Index, 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;
                auto const x0 = x;
 
                alpha = lineSearch()(x, f, g, d, function, p, alpha);
 
                callback(x);
 
                if(valueDiverged(f))
                {
                    Logging::info("BFGS diverged at f(x) =", f);
                    x = x0;
                    break;
                }
 
                if(valueIncreased(f0, f))
                {
                    Logging::info("BFGS made no further progress at f(x) =", f);
                    x = x0;
                    break;
                }
 
                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);
            }
 
            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_BFGS_H