gradient_descent.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_OPTIMIZER_GRADIENT_DESCENT_H
#define pRC_ALGORITHMS_OPTIMIZER_GRADIENT_DESCENT_H
 
#include <prc/algorithms/optimizer/line_search/fixed.hpp>
#include <prc/core/functors/eval.hpp>
#include <prc/core/functors/scalar_product.hpp>
 
namespace pRC::Optimizer
{
    template<class LS = LineSearch::Fixed>
    class GradientDescent
    {
    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 GradientDescent(LS const &lineSearch,
            Size const maxIterations = defaultMaxIterations())
            : mLineSearch(lineSearch)
            , mMaxIterations(maxIterations)
        {
        }
 
        constexpr GradientDescent(
            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("Gradient Descent initial f(x) =", f);
 
            if(projectedGradientConverged(g, tolerance))
            {
                return x;
            }
 
            auto alpha = identity<TX>();
            for(Index iteration = 0;;)
            {
                auto const f0 = f;
                auto const g0 = g;
                auto const x0 = x;
 
                Tensor p = -g;
                auto d = -norm<2, 1>(g)();
 
                alpha = lineSearch()(x, f, g, d, function, p, alpha);
 
                callback(x);
 
                if(valueDiverged(f))
                {
                    Logging::info("Gradient Descent diverged at f(x) =", f);
                    x = x0;
                    break;
                }
 
                /*if(valueIncreased(f0, f))
                {
                    Logging::info("Gradient Descent made no further progress at
                f(x) =", f); x = x0; break;
                }*/
 
                if(++iteration; !(iteration < maxIterations()))
                {
                    Logging::info(
                        "Gradient Descent max iterations reached at f(x) =", f);
                    break;
                }
 
                if(valueConverged(f0, f, tolerance))
                {
                    Logging::info("Gradient Descent converged at f(x) =", f);
                    break;
                }
 
                if(projectedGradientConverged(g, tolerance))
                {
                    Logging::info("Gradient Descent converged at f(x) =", f);
                    break;
                }
 
                Logging::info("Gradient Descent current f(x) =", f);
            }
 
            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_GRADIENT_DESCENT_H