Replace Eigen with TNL in gdcpp.h and examples (69756326) · Commits · Jakub Klinkovský / TNL gradient descent

examples/ackley.cpp

+10 −8

Original line number	Diff line number	Diff line
		@@ -8,9 +8,10 @@ struct Ackley
		Ackley()
		{ }

		double operator()(const Eigen::VectorXd &xval, Eigen::VectorXd &) const
		template <typename Vector>
		double operator()(const Vector &xval, Vector &) const
		{
		assert(xval.size() == 2);
		assert(xval.getSize() == 2);
		double x = xval(0);
		double y = xval(1);
		// Calculate ackley function, but no gradient. Let gradien be estimated
		@@ -23,6 +24,9 @@ struct Ackley

		int main()
		{
		// Define the vector type
		using Vector = TNL::Containers::Vector<double, TNL::Devices::Host, gdc::Index>;

		// Create optimizer object with Ackley functor as objective.
		//
		// You can specify a StepSize functor as template parameter.
		@@ -34,8 +38,7 @@ int main()
		// You can additionally specify a FiniteDifferences functor as template
		// parameter. There are Forward-, Backward- and CentralDifferences
		// available. (Default is CentralDifferences)
		gdc::GradientDescent<double, Ackley,
		gdc::WolfeBacktracking<double>> optimizer;
		gdc::GradientDescent<double, Ackley, gdc::WolfeBacktracking<double>> optimizer;

		// Set number of iterations as stop criterion.
		// Set it to 0 or negative for infinite iterations (default is 0).
		@@ -60,8 +63,7 @@ int main()
		optimizer.setVerbosity(4);

		// Set initial guess.
		Eigen::VectorXd initialGuess(2);
		initialGuess << -2.7, 2.2;
		Vector initialGuess = {-2.7, 2.2};

		// Start the optimization
		auto result = optimizer.minimize(initialGuess);
		@@ -73,7 +75,7 @@ int main()
		std::cout << "Final fval: " << result.fval << std::endl;

		// do something with final x-value
		std::cout << "Final xval: " << result.xval.transpose() << std::endl;
		std::cout << "Final xval: " << result.xval << std::endl;

		return 0;
		}

examples/paraboloid.cpp

+9 −7

Original line number	Diff line number	Diff line
		@@ -3,7 +3,8 @@
		// Implement an objective functor.
		struct Paraboloid
		{
		double operator()(const Eigen::VectorXd &xval, Eigen::VectorXd &gradient)
		template <typename Vector>
		double operator()(const Vector &xval, Vector &gradient) const
		{
		// compute gradient explicitly
		// if gradient calculation is omitted, then the optimizer uses
		@@ -18,6 +19,9 @@ struct Paraboloid

		int main()
		{
		// Define the vector type
		using Vector = TNL::Containers::Vector<double, TNL::Devices::Host, gdc::Index>;

		// Create optimizer object with Paraboloid functor as objective.
		//
		// You can specify a StepSize functor as template parameter.
		@@ -29,8 +33,7 @@ int main()
		// You can additionally specify a FiniteDifferences functor as template
		// parameter. There are Forward-, Backward- and CentralDifferences
		// available. (Default is CentralDifferences)
		gdc::GradientDescent<double, Paraboloid,
		gdc::ConstantStepSize<double>> optimizer;
		gdc::GradientDescent<double, Paraboloid, gdc::ConstantStepSize<double>> optimizer;

		// Set number of iterations as stop criterion.
		// Set it to 0 or negative for infinite iterations (default is 0).
		@@ -58,8 +61,7 @@ int main()
		optimizer.setVerbosity(4);

		// Set initial guess.
		Eigen::VectorXd initialGuess(2);
		initialGuess << 2, 2;
		Vector initialGuess = {2, 2};

		// Start the optimization
		auto result = optimizer.minimize(initialGuess);
		@@ -71,7 +73,7 @@ int main()
		std::cout << "Final fval: " << result.fval << std::endl;

		// do something with final x-value
		std::cout << "Final xval: " << result.xval.transpose() << std::endl;
		std::cout << "Final xval: " << result.xval << std::endl;

		return 0;
		}

include/gdcpp.h

+53 −54

Original line number	Diff line number	Diff line
		@@ -8,15 +8,16 @@
		#ifndef GDCPP_GDCPP_H_
		#define GDCPP_GDCPP_H_

		#include <Eigen/Geometry>
		#include <TNL/Containers/Vector.h>
		#include <limits>
		#include <iostream>
		#include <iomanip>
		#include <functional>
		#include <cassert>

		namespace gdc
		{
		typedef long int Index;
		using Index = long int;

		/** Functor to compute forward differences.
		* Computes the gradient of the objective f(x) as follows:
		@@ -29,8 +30,8 @@ namespace gdc
		class ForwardDifferences
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &)> Objective;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &)>;
		private:
		Scalar eps_;
		Index threads_;
		@@ -66,11 +67,11 @@ namespace gdc
		{
		assert(objective_);

		gradient.resize(xval.size());
		gradient.resize(xval.getSize());
		#pragma omp parallel for num_threads(threads_)
		for(Index i = 0; i < xval.size(); ++i)
		for(Index i = 0; i < xval.getSize(); ++i)
		{
		Vector xvalN = xval;
		Vector xvalN(xval);
		xvalN(i) += eps_;
		Scalar fvalN = objective_(xvalN);

		@@ -90,8 +91,8 @@ namespace gdc
		class BackwardDifferences
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &)> Objective;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &)>;
		private:
		Scalar eps_;
		Index threads_;
		@@ -127,11 +128,11 @@ namespace gdc
		{
		assert(objective_);

		gradient.resize(xval.size());
		gradient.resize(xval.getSize());
		#pragma omp parallel for num_threads(threads_)
		for(Index i = 0; i < xval.size(); ++i)
		for(Index i = 0; i < xval.getSize(); ++i)
		{
		Vector xvalN = xval;
		Vector xvalN(xval);
		xvalN(i) -= eps_;
		Scalar fvalN = objective_(xvalN);
		gradient(i) = (fval - fvalN) / eps_;
		@@ -150,8 +151,8 @@ namespace gdc
		struct CentralDifferences
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &)> Objective;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &)>;
		private:
		Scalar eps_;
		Index threads_;
		@@ -187,12 +188,12 @@ namespace gdc
		{
		assert(objective_);

		Vector fvals(xval.size() * 2);
		Vector fvals(xval.getSize() * 2);
		#pragma omp parallel for num_threads(threads_)
		for(Index i = 0; i < fvals.size(); ++i)
		for(Index i = 0; i < fvals.getSize(); ++i)
		{
		Index idx = i / 2;
		Vector xvalN = xval;
		Vector xvalN(xval);
		if(i % 2 == 0)
		xvalN(idx) += eps_ / 2;
		else
		@@ -201,8 +202,8 @@ namespace gdc
		fvals(i) = objective_(xvalN);
		}

		gradient.resize(xval.size());
		for(Index i = 0; i < xval.size(); ++i)
		gradient.resize(xval.getSize());
		for(Index i = 0; i < xval.getSize(); ++i)
		gradient(i) = (fvals(i * 2) - fvals(i * 2 + 1)) / eps_;
		}
		};
		@@ -211,7 +212,7 @@ namespace gdc
		template<typename Scalar>
		struct NoCallback
		{
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;

		bool operator()(const Index,
		const Vector &,
		@@ -227,9 +228,9 @@ namespace gdc
		class ConstantStepSize
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &, Vector &)> Objective;
		typedef std::function<void(const Vector &, const Scalar, Vector &)> FiniteDifferences;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &, Vector &)>;
		using FiniteDifferences = std::function<void(const Vector &, const Scalar, Vector &)>;
		private:
		Scalar stepSize_;
		public:
		@@ -277,9 +278,9 @@ namespace gdc
		class BarzilaiBorwein
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &, Vector &)> Objective;
		typedef std::function<void(const Vector &, const Scalar, Vector &)> FiniteDifferences;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &, Vector &)>;
		using FiniteDifferences = std::function<void(const Vector &, const Scalar, Vector &)>;

		enum class Method
		{
		@@ -302,8 +303,8 @@ namespace gdc
		{
		auto sk = xval - lastXval_;
		auto yk = gradient - lastGradient_;
		Scalar num = sk.dot(sk);
		Scalar denom = sk.dot(yk);
		Scalar num = TNL::dot(sk, sk);
		Scalar denom = TNL::dot(sk, yk);

		if(denom == 0)
		return 1;
		@@ -316,8 +317,8 @@ namespace gdc
		{
		auto sk = xval - lastXval_;
		auto yk = gradient - lastGradient_;
		Scalar num = sk.dot(yk);
		Scalar denom = yk.dot(yk);
		Scalar num = TNL::dot(sk, yk);
		Scalar denom = TNL::dot(yk, yk);

		if(denom == 0)
		return 1;
		@@ -355,7 +356,7 @@ namespace gdc
		const Vector &gradient)
		{
		Scalar stepSize = 0;
		if(lastXval_.size() == 0)
		if(lastXval_.getSize() == 0)
		{
		stepSize = constStep_;
		}
		@@ -369,9 +370,6 @@ namespace gdc
		case Method::Inverse:
		stepSize = inverseStep(xval, gradient);
		break;
		default:
		assert(false);
		break;
		}
		}

		@@ -395,9 +393,9 @@ namespace gdc
		class ArmijoBacktracking
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &, Vector &)> Objective;
		typedef std::function<void(const Vector &, const Scalar, Vector &)> FiniteDifferences;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &, Vector &)>;
		using FiniteDifferences = std::function<void(const Vector &, const Scalar, Vector &)>;
		protected:
		Scalar decrease_;
		Scalar cArmijo_;
		@@ -411,7 +409,7 @@ namespace gdc
		{
		gradient.resize(0);
		Scalar fval = objective_(xval, gradient);
		if(gradient.size() == 0)
		if(gradient.getSize() == 0)
		finiteDifferences_(xval, fval, gradient);
		return fval;
		}
		@@ -517,7 +515,7 @@ namespace gdc
		xvalN = xval - stepSize * gradient;
		fvalN = evaluateObjective(xvalN, gradientN);

		armijoCondition = fvalN <= fval - cArmijo_ * stepSize * gradient.dot(gradient);
		armijoCondition = fvalN <= fval - cArmijo_ * stepSize * TNL::dot(gradient, gradient);
		secondCondition = computeSecondCondition(stepSize, fval, fvalN, gradient, gradientN);

		++iterations;
		@@ -541,9 +539,9 @@ namespace gdc
		class WolfeBacktracking : public ArmijoBacktracking<Scalar>
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &, Vector &)> Objective;
		typedef std::function<void(const Vector &, const Scalar, Vector &)> FiniteDifferences;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &, Vector &)>;
		using FiniteDifferences = std::function<void(const Vector &, const Scalar, Vector &)>;
		protected:
		Scalar cWolfe_;

		@@ -553,7 +551,7 @@ namespace gdc
		const Vector &gradient,
		const Vector &gradientN)
		{
		return gradient.dot(gradientN) <= cWolfe_ * gradient.dot(gradient);
		return TNL::dot(gradient, gradientN) <= cWolfe_ * TNL::dot(gradient, gradient);
		}
		public:
		WolfeBacktracking()
		@@ -602,9 +600,9 @@ namespace gdc
		class DecreaseBacktracking
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		typedef std::function<Scalar(const Vector &, Vector &)> Objective;
		typedef std::function<void(const Vector &, const Scalar, Vector &)> FiniteDifferences;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;
		using Objective = std::function<Scalar(const Vector &, Vector &)>;
		using FiniteDifferences = std::function<void(const Vector &, const Scalar, Vector &)>;
		private:
		Scalar decrease_;
		Scalar minStep_;
		@@ -698,7 +696,7 @@ namespace gdc
		class GradientDescent
		{
		public:
		typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector;
		using Vector = TNL::Containers::Vector<Scalar, TNL::Devices::Host, Index>;

		struct Result
		{
		@@ -722,7 +720,7 @@ namespace gdc
		{
		gradient.resize(0);
		Scalar fval = objective_(xval, gradient);
		if(gradient.size() == 0)
		if(gradient.getSize() == 0)
		finiteDifferences_(xval, fval, gradient);
		return fval;
		}
		@@ -733,11 +731,11 @@ namespace gdc
		ss1 << std::fixed << std::showpoint << std::setprecision(6);
		std::stringstream ss2;
		ss2 << '[';
		for(Index i = 0; i < vec.size(); ++i)
		for(Index i = 0; i < vec.getSize(); ++i)
		{
		ss1 << vec(i);
		ss2 << std::setfill(' ') << std::setw(10) << ss1.str();
		if(i != vec.size() - 1)
		if(i != vec.getSize() - 1)
		ss2 << ' ';
		ss1.str("");
		}
		@@ -824,12 +822,13 @@ namespace gdc
		[this](const Vector &xval, const Scalar fval, Vector &gradient)
		{ this->finiteDifferences_(xval, fval, gradient); });

		Vector xval = initialGuess;
		Vector xval(initialGuess);
		Vector gradient;
		Scalar fval;
		Scalar gradientLen = minGradientLen_ + 1;
		Scalar stepSize;
		Vector step = Vector::Zero(xval.size());
		Vector step;
		step.resize(xval.getSize(), 0);
		Scalar stepLen = minStepLen_ + 1;
		bool callbackResult = true;

		@@ -841,11 +840,11 @@ namespace gdc
		{
		xval -= step;
		fval = evaluateObjective(xval, gradient);
		gradientLen = gradient.norm();
		gradientLen = TNL::l2Norm(gradient);
		// update step according to step size and momentum
		stepSize = stepSize_(xval, fval, gradient);
		step = momentum_ * step + (1 - momentum_) * stepSize * gradient;
		stepLen = step.norm();
		stepLen = TNL::l2Norm(step);
		// evaluate callback an save its result
		callbackResult = callback_(iterations, xval, fval, gradient);