GradientDescent.hpp

#ifndef _LDAPLUSPLUS_OPTIMIZATION_GRADIENT_DESCENT_HPP_
#define _LDAPLUSPLUS_OPTIMIZATION_GRADIENT_DESCENT_HPP_

#include <functional>
#include <memory>

#include <Eigen/Core>

namespace ldaplusplus {
namespace optimization {


template <typename ProblemType, typename ParameterType>
class LineSearch
{
    public:
        typedef typename ParameterType::Scalar Scalar;

        virtual Scalar search(
            const ProblemType &problem,
            Eigen::Ref<ParameterType> x0,
            const ParameterType &grad_x0,
            const ParameterType &direction
        ) = 0;

        virtual ~LineSearch(){};
};


template <typename ProblemType, typename ParameterType>
class ConstantLineSearch : public LineSearch<ProblemType, ParameterType>
{
    public:
        typedef typename ParameterType::Scalar Scalar;

        ConstantLineSearch(Scalar alpha) : alpha_(alpha) {}

        Scalar search(
            const ProblemType &problem,
            Eigen::Ref<ParameterType> x0,
            const ParameterType &grad_x0,
            const ParameterType &direction
        ) {
            x0 -= alpha_ * direction;

            return problem.value(x0);
        }

    private:
        Scalar alpha_;
};

template <typename ProblemType, typename ParameterType>
class ArmijoLineSearch : public LineSearch<ProblemType, ParameterType>
{
    public:
        typedef typename ParameterType::Scalar Scalar;

        ArmijoLineSearch(Scalar beta=0.001, Scalar tau=0.5) : beta_(beta),
                                                              tau_(tau)
        {}

        Scalar search(
            const ProblemType &problem,
            Eigen::Ref<ParameterType> x0,
            const ParameterType &grad_x0,
            const ParameterType &direction
        ) {
            ParameterType x_copy(x0.rows(), x0.cols());
            Scalar value_x0 = problem.value(x0);
            Scalar decrease = beta_ * (grad_x0.array() * direction.array()).sum();
            Scalar value = value_x0;
            Scalar a = 1.0/tau_;

            while (value > value_x0 - a * decrease) {
                a *= tau_;
                x_copy = x0 - a * direction;
                value = problem.value(x_copy);
            }

            x0 -= a * direction;

            return value;
        }

    private:
        Scalar beta_;
        Scalar tau_;
};


template <typename ProblemType, typename ParameterType>
class GradientDescent
{
    public:
        typedef typename ParameterType::Scalar Scalar;

        GradientDescent(
            std::shared_ptr<LineSearch<ProblemType, ParameterType> > line_search,
            std::function<bool(Scalar, Scalar, size_t)> progress
        ) : line_search_(line_search),
            progress_(progress)
        {}

        void minimize(const ProblemType &problem, Eigen::Ref<ParameterType> x0) {
            // allocate memory for the gradient
            ParameterType grad(x0.rows(), x0.cols());

            // Keep the value in this variable
            Scalar value = problem.value(x0);

            // And the iterations in this one
            size_t iterations = 0;

            // Whether we stop or not is decided by someone else
            while (progress_(value, grad.template lpNorm<Eigen::Infinity>(), iterations++)) {
                problem.gradient(x0, grad);
                value = line_search_->search(problem, x0, grad, grad);
            }
        }

    private:
        std::shared_ptr<LineSearch<ProblemType, ParameterType> > line_search_;
        std::function<bool(Scalar, Scalar, size_t)> progress_;
};


}  // namespace optimization
}  // namespace ldaplusplus

#endif // _LDAPLUSPLUS_OPTIMIZATION_GRADIENT_DESCENT_HPP_