adf.h 13.7 KB
Newer Older
Sebastian Wolf's avatar
Sebastian Wolf committed
1
// Xerus - A General Purpose Tensor Library
2
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
Sebastian Wolf's avatar
Sebastian Wolf 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
// 
// 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 ADF algorithm and its variants.
 */

#pragma once

27
#include "../index.h"
Sebastian Wolf's avatar
Sebastian Wolf committed
28 29 30 31 32
#include "../ttNetwork.h"
#include "../performanceData.h"
#include "../measurments.h"

namespace xerus {
33
	
Sebastian Wolf's avatar
Sebastian Wolf committed
34 35 36
	/**
	 * @brief Wrapper class for all ADF variants.
	 * @details By creating a new object of this class and modifying the member variables, the behaviour of the solver can be modified.
37
	 * This algorithm is a modified implementation of the alternating directional fitting algrothim, first introduced by Grasedyck, Kluge and Kraemer (2015).
Sebastian Wolf's avatar
Sebastian Wolf committed
38 39
	 */
    class ADFVariant {
Sebastian Wolf's avatar
Sebastian Wolf committed
40
		protected:
Perry's avatar
Perry committed
41
			
Sebastian Wolf's avatar
Sebastian Wolf committed
42 43
		template<class MeasurmentSet> 
		class InternalSolver {
44 45 46 47 48 49 50
			/*
			 * General notation for the ADF:
			 * The vector of the measured values is denoted as b
			 * The measurment operator is denoted as A. It holds by definition A(x) = b (for noiseless measurments).
			 * The Operator constructed from the x by removing the component at corePosition is denoted as E.
			 */
			
Sebastian Wolf's avatar
Sebastian Wolf committed
51
		protected:
52
			///@brief Indices for all internal functions.
Sebastian Wolf's avatar
Sebastian Wolf committed
53
			const Index r1, r2, i1;
54
			
55
			///@brief Reference to the current solution (external ownership)
Sebastian Wolf's avatar
Sebastian Wolf committed
56
			TTTensor& x;
57 58
			
			///@brief Degree of the solution.
Sebastian Wolf's avatar
Sebastian Wolf committed
59
			const size_t degree;
60 61
			
			///@brief Maximally allowed ranks.
Sebastian Wolf's avatar
Sebastian Wolf committed
62 63
			const std::vector<size_t> maxRanks;
			
64
			///@brief Reference to the measurment set (external ownership)
Sebastian Wolf's avatar
Sebastian Wolf committed
65
			const MeasurmentSet& measurments;
66 67
			
			///@brief Number of measurments (i.e. measurments.size())
Sebastian Wolf's avatar
Sebastian Wolf committed
68
			const size_t numMeasurments;
69 70
			
			///@brief The two norm of the measured values
Sebastian Wolf's avatar
Sebastian Wolf committed
71 72
			const value_t normMeasuredValues;
			
73
			///@brief Maximal allowed number of iterations (one iteration = one sweep)
Sebastian Wolf's avatar
Sebastian Wolf committed
74
			const size_t maxIterations; 
75 76
			
			///@brief The target residual norm at which the algorithm shall stop.
Sebastian Wolf's avatar
Sebastian Wolf committed
77
			const double targetResidualNorm;
78 79
			
			///@brief Minimal relative decrease of the residual norm ( (oldRes-newRes)/oldRes ) until either the ranks are increased (if allowed) or the algorithm stops.
80
			const double minimalResidualNormDecrease; 
Sebastian Wolf's avatar
Sebastian Wolf committed
81
			
82
			///@brief The current iteration.
Sebastian Wolf's avatar
Sebastian Wolf committed
83
			size_t iteration;
84 85
			
			///@brief Current residual norm. Updated at the beginning of each iteration.
Sebastian Wolf's avatar
Sebastian Wolf committed
86
			double residualNorm;
87 88

			///@brief The residual norm of the last iteration.
Sebastian Wolf's avatar
Sebastian Wolf committed
89 90
			double lastResidualNorm;
			
91
			///@brief The current residual, saved as vector (instead of a order one tensor).
Sebastian Wolf's avatar
Sebastian Wolf committed
92
			std::vector<value_t> residual;
93 94
			
			///@brief The current projected Gradient component. That is E(A^T(Ax-b))
Sebastian Wolf's avatar
Sebastian Wolf committed
95
			Tensor projectedGradientComponent;
Sebastian Wolf's avatar
Sebastian Wolf committed
96
			
Sebastian Wolf's avatar
Sebastian Wolf committed
97 98
			///@brief Ownership holder for a (degree+2)*numMeasurments array of Tensor pointers. (Not used directly)
			std::unique_ptr<Tensor*[]> forwardStackMem;
99 100 101 102 103 104 105
			
			/** @brief Array [numMeasurments][degree]. For positions smaller than the current corePosition and for each measurment, this array contains the pre-computed
			* contraction of the first _ component tensors and the first _ components of the measurment operator. These tensors are deduplicated in the sense that for each unqiue
			* part of the position only one tensor is actually stored, which is why the is an array of pointers. The Tensors at the current corePosition are used as 
			* scatch space. For convinience the underlying array (forwardStackMem) is larger, wherefore also the positions -1 and degree are allow, all poining to a {1} tensor
			* containing 1 as only entry. Note that position degree-1 must not be used.
			**/
Sebastian Wolf's avatar
Sebastian Wolf committed
106
			Tensor* const * const forwardStack;
107
			
Sebastian Wolf's avatar
Sebastian Wolf committed
108 109
			/// @brief Ownership holder for the unqiue Tensors referenced in forwardStack.
			std::unique_ptr<Tensor[]> forwardStackSaveSlots;
110 111
			
			/// @brief Vector containing for each corePosition a vector of the smallest ids of each group of unique forwardStack entries.
Sebastian Wolf's avatar
Sebastian Wolf committed
112 113
			std::vector<std::vector<size_t>> forwardUpdates;
			
114
			
Sebastian Wolf's avatar
Sebastian Wolf committed
115 116
			///@brief Ownership holder for a (degree+2)*numMeasurments array of Tensor pointers. (Not used directly)
			std::unique_ptr<Tensor*[]> backwardStackMem;
117 118 119 120 121 122 123
			
			/** @brief Array [numMeasurments][degree]. For positions larger than the current corePosition and for each measurment, this array contains the pre-computed
			* contraction of the last _ component tensors and the last _ components of the measurment operator. These tensors are deduplicated in the sense that for each unqiue
			* part of the position only one tensor is actually stored, which is why the is an array of pointers. The Tensors at the current corePosition are used as 
			* scratch space. For convinience the underlying array (forwardStackMem) is larger, wherefore also the positions -1 and degree are allow, all poining to a {1} tensor
			* containing 1 as only entry. Note that position zero must not be used.
			**/
Sebastian Wolf's avatar
Sebastian Wolf committed
124
			Tensor* const * const backwardStack;
125

Sebastian Wolf's avatar
Sebastian Wolf committed
126 127
			/// @brief Ownership holder for the unqiue Tensors referenced in backwardStack.
			std::unique_ptr<Tensor[]> backwardStackSaveSlots;
128 129
			
			/// @brief Vector containing for each corePosition a vector of the smallest ids of each group of unique backwardStack entries.
Sebastian Wolf's avatar
Sebastian Wolf committed
130 131
			std::vector<std::vector<size_t>> backwardUpdates;
			
Sebastian Wolf's avatar
Sebastian Wolf committed
132 133
            /// @brief: Norm of each rank one measurment operator
            std::unique_ptr<double[]> measurmentNorms;
134 135
			
			///@brief: Reference to the performanceData object (external ownership)
Sebastian Wolf's avatar
Sebastian Wolf committed
136 137
			PerformanceData& perfData;
			
138
			///@brief calculates the two-norm of the measured values.
Sebastian Wolf's avatar
Sebastian Wolf committed
139 140
			static double calculate_norm_of_measured_values(const MeasurmentSet& _measurments);
			
141
			///@brief Constructes either the forward or backward stack. That is, it determines the groups of partially equale measurments. Therby stetting (forward/backward)- Updates, StackMem and SaveSlot.
Sebastian Wolf's avatar
Sebastian Wolf committed
142
			void construct_stacks(std::unique_ptr< xerus::Tensor[] >& _stackSaveSlot, std::vector< std::vector< size_t > >& _updates, const std::unique_ptr<Tensor*[]>& _stackMem, const bool _forward);
Sebastian Wolf's avatar
Sebastian Wolf committed
143
			
144
			///@brief Resizes the unqiue stack tensors to correspond to the current ranks of x.
Sebastian Wolf's avatar
Sebastian Wolf committed
145 146
			void resize_stack_tensors();
			
147
			///@brief Returns a vector of tensors containing the slices of @a _component where the second dimension is fixed.
Sebastian Wolf's avatar
Sebastian Wolf committed
148
			std::vector<Tensor> get_fixed_components(const Tensor& _component);
Sebastian Wolf's avatar
Sebastian Wolf committed
149
			
150 151 152
			///@brief For each measurment sets the forwardStack at the given _corePosition to the contraction between the forwardStack at the previous corePosition (i.e. -1)
			/// and the given component contracted with the component of the measurment operator. For _corePosition == corePosition and _currentComponent == x.components(corePosition)
			/// this really updates the stack, otherwise it uses the stack as scratch space.
Sebastian Wolf's avatar
Sebastian Wolf committed
153 154
			void update_forward_stack(const size_t _corePosition, const Tensor& _currentComponent);
			
155 156 157 158 159 160
			///@brief For each measurment sets the backwardStack at the given _corePosition to the contraction between the backwardStack at the previous corePosition (i.e. +1)
			/// and the given component contracted with the component of the measurment operator. For _corePosition == corePosition and _currentComponent == x.components(corePosition)
			/// this really updates the stack, otherwise it uses the stack as scratch space.
			void update_backward_stack(const size_t _corePosition, const Tensor& _currentComponent);
			
			///@brief (Re-)Calculates the current residual, i.e. Ax-b.
Sebastian Wolf's avatar
Sebastian Wolf committed
161 162
			void calculate_residual( const size_t _corePosition );
			
163 164 165 166 167
			///@brief Calculates one internal step of calculate_projected_gradient. In particular the dyadic product of the leftStack, the rightStack and the position vector.
			template<class PositionType>
			void perform_dyadic_product(const size_t _localLeftRank, const size_t _localRightRank, const value_t* const _leftPtr,  const value_t* const _rightPtr,  value_t* const _deltaPtr, const value_t _residual, const PositionType& _position, value_t* const _scratchSpace );
	
			
168 169
			///@brief: Calculates the component at _corePosition of the projected gradient from the residual, i.e. E(A^T(b-Ax)).
			void calculate_projected_gradient(const size_t _corePosition);
Sebastian Wolf's avatar
Sebastian Wolf committed
170 171
			
			/**
Sebastian Wolf's avatar
Sebastian Wolf committed
172
			* @brief: Calculates ||P_n (A(E(A^T(b-Ax)))))|| = ||P_n (A(E(A^T(residual)))))|| =  ||P_n (A(E(gradient)))|| for each n, 
Sebastian Wolf's avatar
Sebastian Wolf committed
173 174 175 176 177
			* where P_n sets all entries equals zero except where the index at _corePosition is equals n. In case of RankOneMeasurments,
			* the calculation is not slicewise (only n=0 is set).
			*/
			std::vector<value_t> calculate_slicewise_norm_A_projGrad( const size_t _corePosition);
			
178
			///@brief Updates the current solution x. For SinglePointMeasurments the is done for each slice speratly, for RankOneMeasurments there is only one combined update.
Sebastian Wolf's avatar
Sebastian Wolf committed
179 180
			void update_x(const std::vector<value_t>& _normAProjGrad, const size_t _corePosition);
			
181
			///@brief Basically the complete algorithm, trying to reconstruct x using its current ranks.
Sebastian Wolf's avatar
Sebastian Wolf committed
182 183 184
			void solve_with_current_ranks();
			
		public:
185
			///@brief Tries to solve the reconstruction problem with the current settings.
Sebastian Wolf's avatar
Sebastian Wolf committed
186 187 188 189 190 191 192
			double solve();
			
			InternalSolver(TTTensor& _x, 
							const std::vector<size_t>& _maxRanks, 
							const MeasurmentSet& _measurments, 
							const size_t _maxIteration, 
							const double _targetResidualNorm, 
193
							const double _minimalResidualNormDecrease, 
Sebastian Wolf's avatar
Sebastian Wolf committed
194 195 196
							PerformanceData& _perfData ) : 
				x(_x),
				degree(_x.degree()),
197
				maxRanks(TTTensor::reduce_to_maximal_ranks(_maxRanks, _x.dimensions)),
Sebastian Wolf's avatar
Sebastian Wolf committed
198 199 200 201 202 203 204
				
				measurments(_measurments),
				numMeasurments(_measurments.size()),
				normMeasuredValues(calculate_norm_of_measured_values(_measurments)),
				
				maxIterations(_maxIteration),
				targetResidualNorm(_targetResidualNorm),
205
				minimalResidualNormDecrease(_minimalResidualNormDecrease),
Sebastian Wolf's avatar
Sebastian Wolf committed
206 207 208 209 210 211 212
				
				iteration(0),
				residualNorm(std::numeric_limits<double>::max()), 
				lastResidualNorm(std::numeric_limits<double>::max()),
				
				residual(numMeasurments),
				
Sebastian Wolf's avatar
Sebastian Wolf committed
213
				forwardStackMem(new Tensor*[numMeasurments*(degree+2)]),
Sebastian Wolf's avatar
Sebastian Wolf committed
214 215 216
				forwardStack(forwardStackMem.get()+numMeasurments),
				forwardUpdates(degree),
					
Sebastian Wolf's avatar
Sebastian Wolf committed
217
				backwardStackMem(new Tensor*[numMeasurments*(degree+2)]),
Sebastian Wolf's avatar
Sebastian Wolf committed
218 219 220
				backwardStack(backwardStackMem.get()+numMeasurments),
				backwardUpdates(degree),
				
Sebastian Wolf's avatar
Sebastian Wolf committed
221 222
				measurmentNorms(new double[numMeasurments]),
				
Sebastian Wolf's avatar
Sebastian Wolf committed
223 224
				perfData(_perfData) 
				{
Sebastian Wolf's avatar
Sebastian Wolf committed
225
					_x.require_correct_format();
226 227
					XERUS_REQUIRE(numMeasurments > 0, "Need at very least one measurment.");
					XERUS_REQUIRE(measurments.degree() == degree, "Measurment degree must coincide with x degree.");
Sebastian Wolf's avatar
Sebastian Wolf committed
228 229 230
				}
				
		};
231
		
Sebastian Wolf's avatar
Sebastian Wolf committed
232
	public:
Sebastian Wolf's avatar
Sebastian Wolf committed
233 234
        size_t maxIterations; ///< Maximum number of sweeps to perform. Set to 0 for infinite.
        double targetResidualNorm; ///< Target residual. The algorithm will stop upon reaching a residual smaller than this value.
235
        double minimalResidualNormDecrease; // The minimal relative decrease of the residual per step  ( i.e. (lastResidual-residual)/lastResidual ). If the avg. of the last three steps is smaller than this value, the algorithm stops.
Sebastian Wolf's avatar
Sebastian Wolf committed
236 237
        
		/// fully defining constructor. alternatively ALSVariants can be created by copying a predefined variant and modifying it
Sebastian Wolf's avatar
Sebastian Wolf committed
238
        ADFVariant(const size_t _maxIteration, const double _targetResidual, const double _minimalResidualDecrease)
239
                : maxIterations(_maxIteration), targetResidualNorm(_targetResidual), minimalResidualNormDecrease(_minimalResidualDecrease) { }
Sebastian Wolf's avatar
Sebastian Wolf committed
240 241
        
        /**
Sebastian Wolf's avatar
Sebastian Wolf committed
242 243
		* @brief Tries to reconstruct the (low rank) tensor _x from the given measurments. 
		* @param[in,out] _x On input: an initial guess of the solution, also defining the ranks. On output: The reconstruction found by the algorithm.
244
		* @param _measurments the available measurments, can be either a SinglePointMeasurementSet or RankOneMeasurementSet.
Sebastian Wolf's avatar
Sebastian Wolf committed
245 246 247
		* @param _perfData optinal performanceData object to be used.
		* @returns the residual @f$|P_\Omega(x-b)|_2@f$ of the final @a _x.
		*/
248 249
		template<class MeasurmentSet>
		double operator()(TTTensor& _x, const MeasurmentSet& _measurments, PerformanceData& _perfData) const {
250
			InternalSolver<MeasurmentSet> solver(_x, _x.ranks(), _measurments, maxIterations, targetResidualNorm, minimalResidualNormDecrease, _perfData);
251 252
			return solver.solve();
		}
Sebastian Wolf's avatar
Sebastian Wolf committed
253
		
Sebastian Wolf's avatar
Sebastian Wolf committed
254 255 256
		/**
		* @brief Tries to reconstruct the (low rank) tensor _x from the given measurments. 
		* @param[in,out] _x On input: an initial guess of the solution, may be of smaller rank. On output: The reconstruction found by the algorithm.
257
		* @param _measurments the available measurments, can be either a SinglePointMeasurementSet or RankOneMeasurementSet.
Sebastian Wolf's avatar
Sebastian Wolf committed
258 259 260 261
		* @param _maxRanks the maximal ranks the algorithm may use to decrease the resdiual.
		* @param _perfData optinal performanceData object to be used.
		* @returns the residual @f$|P_\Omega(x-b)|_2@f$ of the final @a _x.
		*/
262 263
		template<class MeasurmentSet>
		double operator()(TTTensor& _x, const MeasurmentSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData) const {
264
			InternalSolver<MeasurmentSet> solver(_x, _maxRanks, _measurments, maxIterations, targetResidualNorm, minimalResidualNormDecrease, _perfData);
265 266 267
			return solver.solve();
		}
	};
Sebastian Wolf's avatar
Sebastian Wolf committed
268
	
Sebastian Wolf's avatar
Sebastian Wolf committed
269
	/// @brief Default variant of the ADF algorithm
Sebastian Wolf's avatar
Sebastian Wolf committed
270
    extern const ADFVariant ADF;
Sebastian Wolf's avatar
Sebastian Wolf committed
271 272
}