steepestDescent.h 6.66 KB
Newer Older
Benjamin Huber's avatar
Benjamin Huber committed
1
// Xerus - A General Purpose Tensor Library
2
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
Benjamin Huber's avatar
Benjamin Huber committed
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

/**
 * @file
 * @brief Header file for the steepest descent algorithms.
 */

#pragma once

#include "../ttNetwork.h"
Benjamin Huber's avatar
Benjamin Huber committed
28
#include "../performanceData.h"
29
#include "retractions.h"
Benjamin Huber's avatar
Benjamin Huber committed
30 31 32

namespace xerus {

33
	
34
	void line_search(TTTensor &_x, value_t &_alpha, const TTTangentVector &_direction, value_t _derivative, value_t &_residual,
35 36 37 38 39
					 std::function<void(TTTensor &, const TTTangentVector &)> _retraction,
					 std::function<value_t()> _calculate_residual,
					 value_t _changeInAlpha = 0.5
	);
	
Benjamin Huber's avatar
Benjamin Huber committed
40 41 42 43 44 45 46
	/**
	* @brief Wrapper class for all steepest descent variants
	* @note only implemented for TTTensors at the moment.
	* @details By creating a new object of this class and modifying the member variables, the behaviour of the solver can be modified.
	*/
	class SteepestDescentVariant {
	protected:
Benjamin Huber's avatar
Benjamin Huber committed
47
		double solve(const TTOperator *_Ap, TTTensor &_x, const TTTensor &_b, size_t _numSteps, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const;
Benjamin Huber's avatar
Benjamin Huber committed
48 49 50 51
	
	public:
		size_t numSteps; ///< maximum number of steps to perform. set to 0 for infinite
		value_t convergenceEpsilon; ///< default value for the change in the residual at which the algorithm assumes it is converged
52
		bool assumeSymmetricPositiveDefiniteOperator; ///< calculates the gradient as b-Ax instead of A^T(b-Ax)
53
		TTOperator *preconditioner;
Benjamin Huber's avatar
Benjamin Huber committed
54
		
55
		TTRetractionII retraction; ///< the retraction to project from point + tangent vector to a new point on the manifold
Benjamin Huber's avatar
Benjamin Huber committed
56 57 58 59
		
		// TODO preconditioner
		
		/// fully defining constructor. alternatively SteepestDescentVariant can be created by copying a predefined variant and modifying it
60 61
		SteepestDescentVariant(size_t _numSteps, value_t _convergenceEpsilon, bool _symPosOp, TTRetractionII _retraction)
				: numSteps(_numSteps), convergenceEpsilon(_convergenceEpsilon),
62
				  assumeSymmetricPositiveDefiniteOperator(_symPosOp), preconditioner(nullptr), retraction(_retraction)
Benjamin Huber's avatar
Benjamin Huber committed
63 64 65
		{ }
		
		/// definition using only the retraction. In the following an operator() including either convergenceEpsilon or numSteps must be called or the algorithm will never terminate
66 67
		SteepestDescentVariant(TTRetractionII _retraction)
				: numSteps(0), convergenceEpsilon(0.0), assumeSymmetricPositiveDefiniteOperator(false), preconditioner(nullptr), retraction(_retraction)
Benjamin Huber's avatar
Benjamin Huber committed
68 69 70 71 72 73 74 75 76 77 78
		{ }
		
		/**
		* call to solve @f$ A\cdot x = b@f$ for @f$ x @f$ (in a least-squares sense)
		* @param _A operator to solve for
		* @param[in,out] _x in: initial guess, out: solution as found by the algorithm
		* @param _b right-hand side of the equation to be solved
		* @param _convergenceEpsilon minimum change in residual / energy under which the algorithm terminates
		* @param _perfData vector of performance data (residuals after every microiteration)
		* @returns the residual @f$|Ax-b|@f$ of the final @a _x
		*/
Benjamin Huber's avatar
Benjamin Huber committed
79
		double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const {
Benjamin Huber's avatar
Benjamin Huber committed
80 81 82 83 84 85 86 87 88 89 90 91
			return solve(&_A, _x, _b, numSteps, _convergenceEpsilon, _perfData);
		}
		
		/**
		* call to solve @f$ A\cdot x = b@f$ for @f$ x @f$ (in a least-squares sense)
		* @param _A operator to solve for
		* @param[in,out] _x in: initial guess, out: solution as found by the algorithm
		* @param _b right-hand side of the equation to be solved
		* @param _numHalfSweeps maximum number of half-sweeps to perform
		* @param _perfData vector of performance data (residuals after every microiteration)
		* @returns the residual @f$|Ax-b|@f$ of the final @a _x
		*/
Benjamin Huber's avatar
Benjamin Huber committed
92
		double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, size_t _numSteps, PerformanceData &_perfData = NoPerfData) const {
Benjamin Huber's avatar
Benjamin Huber committed
93 94 95 96 97 98 99 100 101 102 103
			return solve(&_A, _x, _b, _numSteps, convergenceEpsilon, _perfData);
		}
		
		/**
		* call to solve @f$ A\cdot x = b@f$ for @f$ x @f$ (in a least-squares sense)
		* @param _A operator to solve for
		* @param[in,out] _x in: initial guess, out: solution as found by the algorithm
		* @param _b right-hand side of the equation to be solved
		* @param _perfData vector of performance data (residuals after every microiteration)
		* @returns the residual @f$|Ax-b|@f$ of the final @a _x
		*/
Benjamin Huber's avatar
Benjamin Huber committed
104
		double operator()(const TTOperator &_A, TTTensor &_x, const TTTensor &_b, PerformanceData &_perfData = NoPerfData) const {
Benjamin Huber's avatar
Benjamin Huber committed
105 106 107 108 109 110 111 112 113 114 115
			return solve(&_A, _x, _b, numSteps, convergenceEpsilon, _perfData);
		}
		
		/**
		* call to minimze @f$ \|x - b\|^2 @f$ for @f$ x @f$
		* @param[in,out] _x in: initial guess, out: solution as found by the algorithm
		* @param _b right-hand side of the equation to be solved
		* @param _convergenceEpsilon minimum change in residual / energy under which the algorithm terminates
		* @param _perfData vector of performance data (residuals after every microiteration)
		* @returns the residual @f$|x-b|@f$ of the final @a _x
		*/
Benjamin Huber's avatar
Benjamin Huber committed
116
		double operator()(TTTensor &_x, const TTTensor &_b, value_t _convergenceEpsilon, PerformanceData &_perfData = NoPerfData) const {
Benjamin Huber's avatar
Benjamin Huber committed
117 118 119 120 121 122 123 124 125 126 127
			return solve(nullptr, _x, _b, numSteps, _convergenceEpsilon, _perfData);
		}
		
		/**
		* call to minimze @f$ \|x - b\|^2 @f$ for @f$ x @f$
		* @param[in,out] _x in: initial guess, out: solution as found by the algorithm
		* @param _b right-hand side of the equation to be solved
		* @param _numHalfSweeps maximum number of half-sweeps to perform
		* @param _perfData vector of performance data (residuals after every microiteration)
		* @returns the residual @f$|x-b|@f$ of the final @a _x
		*/
Benjamin Huber's avatar
Benjamin Huber committed
128
		double operator()(TTTensor &_x, const TTTensor &_b, size_t _numSteps, PerformanceData &_perfData = NoPerfData) const {
Benjamin Huber's avatar
Benjamin Huber committed
129 130 131
			return solve(nullptr, _x, _b, _numSteps, convergenceEpsilon, _perfData);
		}
		
Benjamin Huber's avatar
Benjamin Huber committed
132 133
		double operator()(TTTensor &_x, const TTTensor &_b, PerformanceData &_perfData = NoPerfData) const {
			return solve(nullptr, _x, _b, numSteps, convergenceEpsilon, _perfData);
Benjamin Huber's avatar
Benjamin Huber committed
134 135 136 137
		}
	};
	
	/// default variant of the steepest descent algorithm using the lapack solver
Sebastian Wolf's avatar
Sebastian Wolf committed
138
	extern const SteepestDescentVariant SteepestDescent;
Benjamin Huber's avatar
Benjamin Huber committed
139 140
}