als.cpp 20.1 KB
Newer Older
Fuchsi*'s avatar
Fuchsi* committed
1
// Xerus - A General Purpose Tensor Library
Sebastian Wolf's avatar
Sebastian Wolf committed
2
// Copyright (C) 2014-2016 Benjamin Huber and Sebastian Wolf. 
Fuchsi*'s avatar
Fuchsi* committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
// 
// 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.

20
/**
Ben Huber's avatar
Ben Huber committed
21
22
23
* @file
* @brief Implementation of the ALS variants.
*/
24

25
26
#include <xerus/misc/math.h>

27
#include <xerus/algorithms/als.h>
28
#include <xerus/basic.h>
29
#include <xerus/index.h>
30
31
#include <xerus/indexedTensorList.h>
#include <xerus/tensorNetwork.h>
32
#include <xerus/misc/internal.h>
Ben Huber's avatar
Ben Huber committed
33

34
#include <xerus/indexedTensorMoveable.h>
35
#include <xerus/indexedTensor_tensor_factorisations.h>
Baum's avatar
Baum committed
36
37
38

namespace xerus {

Ben Huber's avatar
Ben Huber committed
39
40
41
42
	// -------------------------------------------------------------------------------------------------------------------------
	//                                       local solvers
	// -------------------------------------------------------------------------------------------------------------------------
	
Sebastian Wolf's avatar
Sebastian Wolf committed
43
	void ALSVariant::lapack_solver(const TensorNetwork& _A, std::vector<Tensor>& _x, const TensorNetwork& _b, const ALSAlgorithmicData& _data) {
Ben Huber's avatar
Ben Huber committed
44
45
46
		Tensor A(_A);
		Tensor b(_b);
		Tensor x;
Sebastian Wolf's avatar
Sebastian Wolf committed
47
		
48
		xerus::solve(x, A, b);
Sebastian Wolf's avatar
Sebastian Wolf committed
49
		
50
51
		Index i,j,k,l;
		if (_data.direction == Increasing) {
Sebastian Wolf's avatar
Sebastian Wolf committed
52
53
54
			for (size_t p = 0; p+1 < _data.ALS.sites; ++p) {
				Tensor U, S;
// 				calculate_svd(U, S, x, x, 2, _data.targetRank[_data.currIndex+p], EPSILON); TODO
55
56
57
58
59
60
61
				(U(i^2,j), S(j,k), x(k,l&1)) = SVD(x(i^2,l&2), _data.targetRank[_data.currIndex+p]);
				_x[p] = std::move(U);
				x(j,l&1) = S(j,k) * x(k,l&1);
			}
			_x.back() = std::move(x);
		} else {
			// direction: decreasing index
Sebastian Wolf's avatar
Sebastian Wolf committed
62
63
64
			for (size_t p = _data.ALS.sites-1; p>0; --p) {
				Tensor S, Vt;
// 				calculate_svd(x, S, Vt, x, x.degree()-1, _data.targetRank[_data.currIndex+p-1], EPSILON); TODO
65
66
67
68
69
70
				(x(i&1,j), S(j,k), Vt(k,l&1)) = SVD(x(i&2,l^2), _data.targetRank[_data.currIndex+p-1]);
				_x[p] = std::move(Vt);
				x(i&1,k) = x(i&1,j) * S(j,k);
			}
			_x[0] = std::move(x);
		}
Ben Huber's avatar
Ben Huber committed
71
72
	}
	
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
	void ALSVariant::ASD_solver(const TensorNetwork &_A, std::vector<Tensor> &_x, const TensorNetwork &_b, const ALSAlgorithmicData &_data) {
		// performs a single gradient step, so
		// x = x + alpha * P( A^t (b - Ax) )    or for SPD: x = x + alpha * P( b - Ax )
		// where the projection P is already part of the stacks 
		// and alpha is the exact stepsize for quadratic functionals
		REQUIRE(_data.ALS.sites == 1, "ASD only defined for single site alternation at the moment");
		Tensor grad;
		Index i,j,k;
		grad(i&0) = _b(i&0) - _A(i/2,j/2) * _x[0](j&0);
		value_t alpha;
		if (_data.ALS.assumeSPD) {
			// stepsize alpha = <y,y>/<y,Ay>
			alpha = misc::sqr(frob_norm(grad)) / value_t(grad(i&0) * _A(i/2,j/2) * grad(j&0));
		} else {
			grad(i&0) = _A(j/2,i/2) * grad(j&0);
			// stepsize alpha = <y,y>/<Ay,Ay>
			alpha = frob_norm(grad) / frob_norm(_A(i/2,j/2) * grad(j&0));
		}
		_x[0] += alpha * grad;
	}
	
Ben Huber's avatar
Ben Huber committed
94
95
96
97
98
99
100
	// -------------------------------------------------------------------------------------------------------------------------
	//                                       helper functions
	// -------------------------------------------------------------------------------------------------------------------------
	
	/**
	 * @brief Finds the range of notes that need to be optimized and orthogonalizes @a _x properly
	 * @details finds full-rank nodes (these can wlog be set to identity and need not be optimized)
Sebastian Wolf's avatar
Sebastian Wolf committed
101
	 * requires canonicalizeAtTheEnd and corePosAtTheEnd to be set
Ben Huber's avatar
Ben Huber committed
102
103
104
105
106
	 * sets optimizedRange
	 * modifies x
	 */
	void ALSVariant::ALSAlgorithmicData::prepare_x_for_als() {
		const size_t d = x.degree();
107
108
109
110
111
		Index r1,r2,n1,cr1;
		
		size_t firstOptimizedIndex = 0;
		size_t dimensionProd = 1;
		while (firstOptimizedIndex + 1 < d) {
Ben Huber's avatar
Ben Huber committed
112
			const size_t localDim = x.dimensions[firstOptimizedIndex];
113
			size_t newDimensionProd = dimensionProd * localDim;
Ben Huber's avatar
Ben Huber committed
114
			if (x.rank(firstOptimizedIndex) < newDimensionProd) {
115
116
117
				break;
			}
			
Ben Huber's avatar
Ben Huber committed
118
			Tensor& curComponent = x.component(firstOptimizedIndex);
119
			curComponent.reinterpret_dimensions({curComponent.dimensions[0]*curComponent.dimensions[1], curComponent.dimensions[2]});
Ben Huber's avatar
Ben Huber committed
120
121
			curComponent(r1,n1,r2) = curComponent(r1,cr1) * x.get_component(firstOptimizedIndex+1)(cr1,n1,r2);
			x.set_component(firstOptimizedIndex+1, std::move(curComponent));
122
123
			
			//TODO sparse
Ben Huber's avatar
Ben Huber committed
124
			x.set_component(firstOptimizedIndex, Tensor(
125
126
127
128
				{dimensionProd, localDim, newDimensionProd},
				[&](const std::vector<size_t> &_idx){
					if (_idx[0]*localDim + _idx[1] == _idx[2]) {
						return 1.0;
129
130
					} 
					return 0.0;
131
132
				})
			);
133
			
Ben Huber's avatar
Ben Huber committed
134
			x.require_correct_format();
135
136
137
138
			
			firstOptimizedIndex += 1;
			dimensionProd = newDimensionProd;
		}
139
		
140
141
		size_t firstNotOptimizedIndex = d;
		dimensionProd = 1;
Ben Huber's avatar
Ben Huber committed
142
		while (firstNotOptimizedIndex > firstOptimizedIndex+ALS.sites) {
Ben Huber's avatar
Ben Huber committed
143
			const size_t localDim = x.dimensions[firstNotOptimizedIndex-1];
144
			size_t newDimensionProd = dimensionProd * localDim;
Ben Huber's avatar
Ben Huber committed
145
			if (x.rank(firstNotOptimizedIndex-2) < newDimensionProd) {
146
147
148
				break;
			}
			
Ben Huber's avatar
Ben Huber committed
149
			Tensor& curComponent = x.component(firstNotOptimizedIndex-1);
150
			curComponent.reinterpret_dimensions({curComponent.dimensions[0], curComponent.dimensions[1] * curComponent.dimensions[2]});
Ben Huber's avatar
Ben Huber committed
151
152
			curComponent(r1,n1,r2) = x.get_component(firstNotOptimizedIndex-2)(r1,n1,cr1) * curComponent(cr1,r2);
			x.set_component(firstNotOptimizedIndex-2, std::move(curComponent));
153
154
			
			//TODO sparse
Ben Huber's avatar
Ben Huber committed
155
			x.set_component(firstNotOptimizedIndex-1, Tensor(
156
157
				{newDimensionProd, localDim, dimensionProd},
				[&](const std::vector<size_t> &_idx){
158
					if (_idx[0] == _idx[1]*dimensionProd + _idx[2]) {
159
						return 1.0;
160
161
					} 
					return 0.0;
162
163
				})
			);
164

Ben Huber's avatar
Ben Huber committed
165
			x.require_correct_format();
166
167
168
169
170
			
			firstNotOptimizedIndex -= 1;
			dimensionProd = newDimensionProd;
		}
		
Sebastian Wolf's avatar
Sebastian Wolf committed
171
		if (canonicalizeAtTheEnd && corePosAtTheEnd < firstOptimizedIndex) {
Ben Huber's avatar
Ben Huber committed
172
			x.assume_core_position(firstOptimizedIndex);
173
		} else {
Sebastian Wolf's avatar
Sebastian Wolf committed
174
			if (canonicalizeAtTheEnd && corePosAtTheEnd >= firstNotOptimizedIndex) {
Ben Huber's avatar
Ben Huber committed
175
				x.assume_core_position(firstNotOptimizedIndex-1);
176
177
			}
			
Ben Huber's avatar
Ben Huber committed
178
			x.move_core(firstOptimizedIndex, true);
179
		}
180
		
Ben Huber's avatar
Ben Huber committed
181
		optimizedRange = std::pair<size_t, size_t>(firstOptimizedIndex, firstNotOptimizedIndex);
182
	}
Baum's avatar
Baum committed
183

184
185
186
187
	TensorNetwork ALSVariant::ALSAlgorithmicData::localOperatorSlice(size_t _pos) {
		//TODO optimization: create these networks without indices
		Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3;
		TensorNetwork res;
188
		INTERNAL_CHECK(A, "ie");
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
		if (ALS.assumeSPD) {
			res(r1,r2,r3, cr1,cr2,cr3) = x.get_component(_pos)(r1, n1, cr1) 
										* A->get_component(_pos)(r2, n1, n2, cr2) 
										* x.get_component(_pos)(r3, n2, cr3);
		} else {
			res(r1,r2,r3,r4, cr1,cr2,cr3, cr4) = x.get_component(_pos)(r1, n1, cr1) 
										* A->get_component(_pos)(r2, n2, n1, cr2) 
										* A->get_component(_pos)(r3, n2, n3, cr3) 
										* x.get_component(_pos)(r4, n3, cr4);
		}
		return res;
	}
	
	TensorNetwork ALSVariant::ALSAlgorithmicData::localRhsSlice(size_t _pos) {
		//TODO optimization: create these networks without indices
		Index cr1, cr2, cr3, r1, r2, r3, n1, n2;
		TensorNetwork res;
206
		if (ALS.assumeSPD || (A == nullptr)) {
207
208
209
210
211
212
213
214
215
216
			res(r1,r2, cr1,cr2) = b.get_component(_pos)(r1, n1, cr1) 
									* x.get_component(_pos)(r2, n1, cr2);
		} else {
			res(r1,r2,r3, cr1,cr2,cr3) = b.get_component(_pos)(r1, n1, cr1) 
												* A->get_component(_pos)(r2, n1, n2, cr2) 
												* x.get_component(_pos)(r3, n2, cr3);
		}
		return res;
	}
	
Ben Huber's avatar
Ben Huber committed
217
	void ALSVariant::ALSAlgorithmicData::prepare_stacks() {
Sebastian Wolf's avatar
Sebastian Wolf committed
218
		const size_t d = x.degree();
219
220
221
222
		Index r1,r2;
		
		Tensor tmpA;
		Tensor tmpB;
223
		if (ALS.assumeSPD || (A == nullptr)) {
224
225
226
227
228
229
			tmpA = Tensor::ones({1,1,1});
			tmpB = Tensor::ones({1,1});
		} else {
			tmpA = Tensor::ones({1,1,1,1});
			tmpB = Tensor::ones({1,1,1});
		}
Ben Huber's avatar
Ben Huber committed
230
231
		
		
232
233
234
235
		localOperatorCache.left.emplace_back(tmpA);
		localOperatorCache.right.emplace_back(tmpA);
		rhsCache.left.emplace_back(tmpB);
		rhsCache.right.emplace_back(tmpB);
Ben Huber's avatar
Ben Huber committed
236
237
		
		for (size_t i = d-1; i > optimizedRange.first + ALS.sites - 1; --i) {
238
			if (A != nullptr) {
239
240
				tmpA(r1&0) = localOperatorCache.right.back()(r2&0) * localOperatorSlice(i)(r1/2, r2/2);
				localOperatorCache.right.emplace_back(tmpA);
Ben Huber's avatar
Ben Huber committed
241
			}
242
243
			tmpB(r1&0) = rhsCache.right.back()(r2&0) * localRhsSlice(i)(r1/2, r2/2);
			rhsCache.right.emplace_back(tmpB);
Ben Huber's avatar
Ben Huber committed
244
245
		}
		for (size_t i = 0; i < optimizedRange.first; ++i) {
246
			if (A != nullptr) {
247
248
				tmpA(r2&0) = localOperatorCache.left.back()(r1&0) * localOperatorSlice(i)(r1/2, r2/2);
				localOperatorCache.left.emplace_back((tmpA));
Ben Huber's avatar
Ben Huber committed
249
			}
250
251
			tmpB(r2&0) = rhsCache.left.back()(r1&0) * localRhsSlice(i)(r1/2, r2/2);
			rhsCache.left.emplace_back(tmpB);
Ben Huber's avatar
Ben Huber committed
252
253
254
255
		}
	}
	
	void ALSVariant::ALSAlgorithmicData::choose_energy_functional() {
256
		if (A != nullptr) {
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
			if (ALS.assumeSPD) {
				residual_f = [&](){
					Index n1, n2;
					return frob_norm((*A)(n1/2,n2/2)*x(n2&0) - b(n1&0));
				};
				if (ALS.useResidualForEndCriterion) {
					energy_f = residual_f;
				} else {
					energy_f = [&](){
						Index r1,r2;
						Tensor res;
						// 0.5*<x,Ax> - <x,b>
						TensorNetwork xAx = localOperatorCache.left.back();
						TensorNetwork bx = rhsCache.left.back();
						for (size_t i=0; i<ALS.sites; ++i) {
							xAx(r2&0) = xAx(r1&0) * localOperatorSlice(currIndex+i)(r1/2, r2/2);
							bx(r2&0) = bx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
						}
						res() = 0.5*xAx(r1&0) * localOperatorCache.right.back()(r1&0)
								- bx(r1&0) * rhsCache.right.back()(r1&0);
Sebastian Wolf's avatar
Sebastian Wolf committed
277
						return res.frob_norm();
278
279
					};
				}
Ben Huber's avatar
Ben Huber committed
280
			} else {
281
282
283
				// not Symmetric pos def
				residual_f = [&](){
					Index r1,r2;
Ben Huber's avatar
Ben Huber committed
284
					Tensor res;
285
286
287
					// <Ax,Ax> - 2 * <Ax,b> + <b,b>
					TensorNetwork xAtAx = localOperatorCache.left.back();
					TensorNetwork bAx = rhsCache.left.back();
288
					for (size_t i=0; i<ALS.sites; ++i) {
289
290
						xAtAx(r2&0) = xAtAx(r1&0) * localOperatorSlice(currIndex+i)(r1/2, r2/2);
						bAx(r2&0) = bAx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
291
					}
292
293
294
					res() = xAtAx(r1&0) * localOperatorCache.right.back()(r1&0)
							- 2 * bAx(r1&0) * rhsCache.right.back()(r1&0);
					return res[0] + misc::sqr(normB);
Ben Huber's avatar
Ben Huber committed
295
				};
296
				energy_f = residual_f;
Ben Huber's avatar
Ben Huber committed
297
298
			}
		} else {
299
			// no operator A given
Ben Huber's avatar
Ben Huber committed
300
			residual_f = [&](){
Ben Huber's avatar
Ben Huber committed
301
				return frob_norm(x - b);
Ben Huber's avatar
Ben Huber committed
302
303
304
305
306
			};
			if (ALS.useResidualForEndCriterion) {
				energy_f = residual_f;
			} else {
				energy_f = [&](){
307
					Index r1,r2;
Ben Huber's avatar
Ben Huber committed
308
309
					Tensor res;
					// 0.5*<x,Ax> - <x,b> = 0.5*|x_i|^2 - <x,b>
310
					TensorNetwork bx = rhsCache.left.back();
311
					for (size_t i=0; i<ALS.sites; ++i) {
312
						bx(r2&0) = bx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
313
					}
314
315
					res() = 0.5*x.get_component(currIndex)(r1&0) * x.get_component(currIndex)(r1&0) 
							- bx(r1&0) * rhsCache.right.back()(r1&0);
Ben Huber's avatar
Ben Huber committed
316
317
318
319
320
321
					return res[0];
				};
			}
		}
	}
	
Ben Huber's avatar
Ben Huber committed
322
	ALSVariant::ALSAlgorithmicData::ALSAlgorithmicData(const ALSVariant &_ALS, const TTOperator *_A, TTTensor &_x, const TTTensor &_b) 
Ben Huber's avatar
Ben Huber committed
323
		: ALS(_ALS), A(_A), x(_x), b(_b)
324
325
		, targetRank(_x.ranks())
		, normB(frob_norm(_b))
Sebastian Wolf's avatar
Sebastian Wolf committed
326
		, canonicalizeAtTheEnd(_x.canonicalized)
327
328
329
330
331
332
		, corePosAtTheEnd(_x.corePosition)
		, lastEnergy2(1e102)
		, lastEnergy(1e101)
		, energy(1e100)
		, halfSweepCount(0)
		, direction(Increasing)
Ben Huber's avatar
Ben Huber committed
333
334
335
336
337
338
339
340
	{
		prepare_x_for_als();
		prepare_stacks();
		currIndex = optimizedRange.first;
		choose_energy_functional();
	}

	void ALSVariant::ALSAlgorithmicData::move_to_next_index() {
341
		Index r1,r2;
Ben Huber's avatar
Ben Huber committed
342
343
		Tensor tmpA, tmpB;
		if (direction == Increasing) {
344
			INTERNAL_CHECK(currIndex+ALS.sites < optimizedRange.second, "ie " << currIndex << " " << ALS.sites << " " << optimizedRange.first << " " << optimizedRange.second);
345
346
347
348
			// Move core to next position (assumed to be done by the solver if sites > 1)
			if (ALS.sites == 1) {
				x.move_core(currIndex+1, true);
			}
Ben Huber's avatar
Ben Huber committed
349
350
			
			// Move one site to the right
351
			if (A != nullptr) {
352
353
354
				localOperatorCache.right.pop_back();
				tmpA(r2&0) = localOperatorCache.left.back()(r1&0) * localOperatorSlice(currIndex)(r1/2, r2/2);
				localOperatorCache.left.emplace_back(std::move(tmpA));
Ben Huber's avatar
Ben Huber committed
355
356
			}
			
357
358
359
			rhsCache.right.pop_back();
			tmpB(r2&0) = rhsCache.left.back()(r1&0) * localRhsSlice(currIndex)(r1/2, r2/2);
			rhsCache.left.emplace_back(std::move(tmpB));
Ben Huber's avatar
Ben Huber committed
360
361
			currIndex++;
		} else {
362
			INTERNAL_CHECK(currIndex > optimizedRange.first, "ie");
363
364
365
366
			// Move core to next position (assumed to be done by the solver if sites > 1)
			if (ALS.sites == 1) {
				x.move_core(currIndex-1, true);
			}
Ben Huber's avatar
Ben Huber committed
367
368
			
			// move one site to the left
369
			if (A != nullptr) {
370
371
372
				localOperatorCache.left.pop_back();
				tmpA(r1&0) = localOperatorCache.right.back()(r2&0) * localOperatorSlice(currIndex)(r1/2, r2/2);
				localOperatorCache.right.emplace_back(std::move(tmpA));
Ben Huber's avatar
Ben Huber committed
373
374
			}
			
375
376
377
			rhsCache.left.pop_back();
			tmpB(r1&0) = rhsCache.right.back()(r2&0) * localRhsSlice(currIndex)(r1/2, r2/2);
			rhsCache.right.emplace_back(std::move(tmpB));
Ben Huber's avatar
Ben Huber committed
378
379
380
381
382
			currIndex--;
		}
	}
	
	
383
	TensorNetwork ALSVariant::construct_local_operator(ALSVariant::ALSAlgorithmicData& _data) const {
384
		INTERNAL_CHECK(_data.A, "IE");
385
386
387
		Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3, n4, x;
		TensorNetwork ATilde = _data.localOperatorCache.left.back();
		if (assumeSPD) {
Sebastian Wolf's avatar
Sebastian Wolf committed
388
			for (size_t p=0; p<sites; ++p) {
389
390
391
392
393
394
395
396
397
398
399
400
401
402
				ATilde(n1^(p+1), n2, r2, n3^(p+1), n4) = ATilde(n1^(p+1), r1, n3^(p+1)) * _data.A->get_component(_data.currIndex+p)(r1, n2, n4, r2);
			}
			ATilde(n1^(sites+1), n2, n3^(sites+1), n4) = ATilde(n1^(sites+1), r1, n3^(sites+1)) * _data.localOperatorCache.right.back()(n2, r1, n4);
		} else {
			for (size_t p=0;  p<sites; ++p) {
				ATilde(n1^(p+1),n2, r3,r4, n3^(p+1),n4) = ATilde(n1^(p+1), r1,r2, n3^(p+1))
						* _data.A->get_component(_data.currIndex+p)(r1, x, n2, r3)
						* _data.A->get_component(_data.currIndex+p)(r2, x, n4, r4);
			}
			ATilde(n1^(sites+1), n2, n3^(sites+1), n4) = ATilde(n1^(sites+1), r1^2, n3^(sites+1)) * _data.localOperatorCache.right.back()(n2, r1^2, n4);
		}
		return ATilde;
	}
	
Sebastian Wolf's avatar
Sebastian Wolf committed
403
	
404
	TensorNetwork ALSVariant::construct_local_RHS(ALSVariant::ALSAlgorithmicData& _data) const {
405
406
		Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3, n4, x;
		TensorNetwork BTilde;
407
		if (assumeSPD || (_data.A == nullptr)) {
408
409
410
411
412
413
414
415
			BTilde(n1,r1) = _data.rhsCache.left.back()(r1,n1);
			for (size_t p=0; p<sites; ++p) {
				BTilde(n1^(p+1), n2, cr1) = BTilde(n1^(p+1), r1) * _data.b.get_component(_data.currIndex+p)(r1, n2, cr1);
			}
			BTilde(n1^(sites+1),n2) = BTilde(n1^(sites+1), r1) * _data.rhsCache.right.back()(r1,n2);
		} else {
			BTilde(n1,r1^2) = _data.rhsCache.left.back()(r1^2,n1);
			for (size_t p=0; p<sites; ++p) {
416
417
418
				BTilde(n1^(p+1), n3, cr1, cr2) = BTilde(n1^(p+1), r1, r2) 
					* _data.b.get_component(_data.currIndex+p)(r1, n2, cr1)
					* _data.A->get_component(_data.currIndex+p)(r2, n2, n3, cr2);
419
420
421
422
423
424
425
			}
			BTilde(n1^(sites+1),n2) = BTilde(n1^(sites+1), r1^2) * _data.rhsCache.right.back()(r1^2,n2);
		}
		return BTilde;
	}


426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
	bool ALSVariant::check_for_end_of_sweep(ALSAlgorithmicData& _data, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData) const {
		if ((_data.direction == Decreasing && _data.currIndex==_data.optimizedRange.first) 
			|| (_data.direction == Increasing && _data.currIndex==_data.optimizedRange.second-sites)) 
		{
			LOG(ALS, "Sweep Done");
			_data.halfSweepCount += 1;
			
			_data.lastEnergy2 = _data.lastEnergy;
			_data.lastEnergy = _data.energy;
			_data.energy = _data.energy_f();
			
			if (_perfData) {
				size_t flags = _data.direction == Increasing ? FLAG_FINISHED_HALFSWEEP : FLAG_FINISHED_FULLSWEEP;
				if (!useResidualForEndCriterion) {
					_perfData.stop_timer();
					value_t residual = _data.residual_f();
					_perfData.continue_timer();
					_perfData.add(residual, _data.x, flags);
				} else {
					_perfData.add(_data.energy, _data.x, flags);
				}
			}
			
			// Conditions for loop termination
			if (_data.halfSweepCount == _numHalfSweeps 
				|| std::abs(_data.lastEnergy-_data.energy) < _convergenceEpsilon 
				|| std::abs(_data.lastEnergy2-_data.energy) < _convergenceEpsilon 
				|| (_data.optimizedRange.second - _data.optimizedRange.first<=sites)) 
			{
				// we are done! yay
Sebastian Wolf's avatar
Sebastian Wolf committed
456
				LOG(ALS, "ALS done, " << _data.energy << " " << _data.lastEnergy << " " << std::abs(_data.lastEnergy2-_data.energy) << " " << std::abs(_data.lastEnergy-_data.energy) << " < " << _convergenceEpsilon);
Sebastian Wolf's avatar
Sebastian Wolf committed
457
				if (_data.canonicalizeAtTheEnd && preserveCorePosition) {
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
					_data.x.move_core(_data.corePosAtTheEnd, true);
				}
				return true;
			}
				
			// change walk direction
			_data.direction = _data.direction == Increasing ? Decreasing : Increasing;
		} else {
			// we are not done with the sweep, just calculate data for the perfdata
			if (_perfData) {
				_perfData.stop_timer();
				value_t residual = _data.residual_f();
				_perfData.continue_timer();
				_perfData.add(residual, _data.x, 0);
			}
		}
		return false;
	}
	
477
	
Ben Huber's avatar
Ben Huber committed
478
479
480
481
482
483
	// -------------------------------------------------------------------------------------------------------------------------
	//                                   the actual algorithm
	// -------------------------------------------------------------------------------------------------------------------------
	
	
	double ALSVariant::solve(const TTOperator *_Ap, TTTensor &_x, const TTTensor &_b, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData) const {
484
		LOG(ALS, "ALS("<< sites <<") called");
485
		#ifndef XERUS_DISABLE_RUNTIME_CHECKS
Sebastian Wolf's avatar
Sebastian Wolf committed
486
487
			_x.require_correct_format();
			_b.require_correct_format();
Sebastian Wolf's avatar
Sebastian Wolf committed
488
489
490
			REQUIRE(_x.degree() > 0, "");
			REQUIRE(_x.dimensions == _b.dimensions, "");
			
491
			if (_Ap != nullptr) {
492
493
				_Ap->require_correct_format();
				REQUIRE(_Ap->dimensions.size() == _b.dimensions.size()*2, "");
494
				for (size_t i=0; i<_x.dimensions.size(); ++i) {
495
496
					REQUIRE(_Ap->dimensions[i] == _x.dimensions[i], "");
					REQUIRE(_Ap->dimensions[i+_Ap->degree()/2] == _x.dimensions[i], "");
497
				}
Sebastian Wolf's avatar
Sebastian Wolf committed
498
			}
Ben Huber's avatar
Ben Huber committed
499
500
		#endif
		
501
		if (_Ap != nullptr) {
502
			_perfData << "ALS for ||A*x - b||^2, x.dimensions: " << _x.dimensions << '\n'
503
					<< "A.ranks: " << _Ap->ranks() << '\n'
504
505
506
507
508
509
510
511
512
513
514
515
					<< "x.ranks: " << _x.ranks() << '\n'
					<< "b.ranks: " << _b.ranks() << '\n'
					<< "maximum number of half sweeps: " << _numHalfSweeps << '\n'
					<< "convergence epsilon: " << _convergenceEpsilon << '\n';
		} else {
			_perfData << "ALS for ||x - b||^2, x.dimensions: " << _x.dimensions << '\n'
					<< "x.ranks: " << _x.ranks() << '\n'
					<< "b.ranks: " << _b.ranks() << '\n'
					<< "maximum number of half sweeps: " << _numHalfSweeps << '\n'
					<< "convergence epsilon: " << _convergenceEpsilon << '\n';
		}
		_perfData.start();
516
		
Ben Huber's avatar
Ben Huber committed
517
		ALSAlgorithmicData data(*this, _Ap, _x, _b);
Sebastian Wolf's avatar
Sebastian Wolf committed
518
		
519
		data.energy = data.energy_f();
Ben Huber's avatar
Ben Huber committed
520
521
		
		if (_perfData) {
522
			_perfData.stop_timer();
Ben Huber's avatar
Ben Huber committed
523
			value_t residual = data.residual_f();
524
			_perfData.continue_timer();
Ben Huber's avatar
Ben Huber committed
525
526
			_perfData.add(residual, _x, FLAG_FINISHED_FULLSWEEP);
		}
527
		
Ben Huber's avatar
Ben Huber committed
528
529
		while (true) {
			LOG(ALS, "Starting to optimize index " << data.currIndex);
Sebastian Wolf's avatar
Sebastian Wolf committed
530
			
531
			// update current component tensor
532
			if (_Ap != nullptr) {
533
534
535
536
				std::vector<Tensor> tmpX;
				for (size_t p=0; p<sites; ++p) {
					tmpX.emplace_back(_x.get_component(data.currIndex+p));
				}
537
				localSolver(construct_local_operator(data), tmpX, construct_local_RHS(data), data);
538
539
				for (size_t p=0; p<sites; ++p) {
					_x.set_component(data.currIndex+p, std::move(tmpX[p]));
540
541
				}
			} else {
542
543
				//TODO?
				REQUIRE(sites==1, "approximation dmrg not implemented yet");
544
				_x.component(data.currIndex) = Tensor(construct_local_RHS(data));
Sebastian Wolf's avatar
Sebastian Wolf committed
545
546
			}
			
547
			if(check_for_end_of_sweep(data, _numHalfSweeps, _convergenceEpsilon, _perfData)) {
Sebastian Wolf's avatar
Sebastian Wolf committed
548
				return data.energy; // TODO residual?
Ben Huber's avatar
Ben Huber committed
549
			}
Sebastian Wolf's avatar
Sebastian Wolf committed
550
			
Ben Huber's avatar
Ben Huber committed
551
552
553
			data.move_to_next_index();
		}
	}
Benjamin Huber's avatar
Benjamin Huber committed
554
	
Sebastian Wolf's avatar
Sebastian Wolf committed
555
	
556
557
	const ALSVariant ALS(1, 0, ALSVariant::lapack_solver, false);
	const ALSVariant ALS_SPD(1, 0, ALSVariant::lapack_solver, true);
Sebastian Wolf's avatar
Sebastian Wolf committed
558
	
559
560
	const ALSVariant DMRG(2, 0, ALSVariant::lapack_solver, false);
	const ALSVariant DMRG_SPD(2, 0, ALSVariant::lapack_solver, true);
561
562
563
	
	const ALSVariant ASD(1, 0, ALSVariant::ASD_solver, false);
	const ALSVariant ASD_SPD(1, 0, ALSVariant::ASD_solver, true);
564
} // namespace xerus