fullTensor_solve.cxx 7.28 KB
Newer Older
Baum's avatar
Baum 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. 
Baum's avatar
Baum 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
21
22

#include<xerus.h>

23
#include "../../include/xerus/test/test.h"
24
#include "../../include/xerus/misc/internal.h"
25
using namespace xerus;
Baum's avatar
Baum committed
26

27
static misc::UnitTest tensor_solve("Tensor", "solve_Ax_equals_b", [](){
Baum's avatar
Baum committed
28
29
    Index i,j,k;
          
Sebastian Wolf's avatar
Sebastian Wolf committed
30
    Tensor A1({4,2,2});
Baum's avatar
Baum committed
31
32
33
34
35
    A1[{0,0,0}] = 1;
    A1[{1,0,1}] = 1;
    A1[{2,1,0}] = 1;
    A1[{3,1,1}] = 1;
    
Sebastian Wolf's avatar
Sebastian Wolf committed
36
    Tensor b1({4});
37
38
39
40
    b1[0] = 73;
    b1[1] = -73;
    b1[2] = 128;
    b1[3] = 93;
Baum's avatar
Baum committed
41
    
Sebastian Wolf's avatar
Sebastian Wolf committed
42
    Tensor x1({2,2});
Baum's avatar
Baum committed
43
44
45
46
    
    
    x1(k,i) = (b1(j)/A1(j,k,i));
    
47
48
49
50
    TEST((b1[0] - x1[{0,0}]) < 1e-14);
    TEST((b1[1] - x1[{0,1}]) < 1e-14);
    TEST((b1[2] - x1[{1,0}]) < 1e-14);
    TEST((b1[3] - x1[{1,1}]) < 1e-14);
Baum's avatar
Baum committed
51
    
Sebastian Wolf's avatar
Sebastian Wolf committed
52
    Tensor A2({4,2,2});
Baum's avatar
Baum committed
53
54
55
56
57
    A2[{0,0,0}] = 1;
    A2[{1,0,1}] = 1;
    A2[{2,1,0}] = 0;
    A2[{3,1,1}] = 0;
    
Sebastian Wolf's avatar
Sebastian Wolf committed
58
    Tensor b2({4});
59
60
61
62
    b2[0] = 73;
    b2[1] = -73;
    b2[2] = 0;
    b2[3] = 0;
Baum's avatar
Baum committed
63
    
Sebastian Wolf's avatar
Sebastian Wolf committed
64
    Tensor x2({2,2});
Baum's avatar
Baum committed
65
66
67
    
    x2(k,i) = (b2(j)/A2(j,k,i));
    
68
69
    TEST((b2[0] - x2[{0,0}]) < 1e-14);
    TEST((b2[1] - x2[{0,1}]) < 1e-14);
Baum's avatar
Baum committed
70
71
    TEST((x2[{1,0}]) < 1e-14);
    TEST((x2[{1,1}]) < 1e-14);
72
});
73

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

static misc::UnitTest solve_vs_lsqr("Tensor", "solve vs least squares", [](){
	const size_t N = 500;
	Tensor A({N, N});
	for (size_t i=0; i<N; ++i) {
		if (i>0) A[{i, i-1}] = -1;
		if (i<N-1) A[{i, i+1}] = -1;
		A[{i,i}] = 2;
	}
	A[0] = 1;
	
	Tensor B = Tensor::random({N});

	Index i,j,k;
	
	// sparse
// 	LOG(blabla, "sparse start");
	Tensor X;
	solve(X, A, B);
// 	LOG(blabla, "sparse end");
	MTEST(frob_norm(A(i,j)*X(j)-B(i)) < 1e-10, "1 " << frob_norm(A(i,j)*X(j)-B(i)));
	
	A.use_dense_representation();
	// dense (cholesky)
// 	LOG(blabla, "chol start");
	Tensor X2;
	solve(X2, A, B);
// 	LOG(blabla, "chol end");
	MTEST(frob_norm(A(i,j)*X(j)-B(i)) < 1e-10, "2 " << frob_norm(A(i,j)*X(j)-B(i)));
	
	MTEST(approx_equal(X, X2, 1e-10), X.frob_norm() << " " << X2.frob_norm() << " diff: " << frob_norm(X-X2));
	
	// dense (LDL)
	A[0] = -0.9;
// 	LOG(blabla, "LDL start");
	solve(X, A, B);
// 	LOG(blabla, "LDL end");
	MTEST(frob_norm(A(i,j)*X(j)-B(i)) < 1e-10, "3 " << frob_norm(A(i,j)*X(j)-B(i)));
	
	// dense (LU)
	A[1] = -0.9;
// 	LOG(blabla, "LU start");
	solve(X, A, B);
// 	LOG(blabla, "LU end");
	MTEST(frob_norm(A(i,j)*X(j)-B(i)) < 1e-10, "4 " << frob_norm(A(i,j)*X(j)-B(i)));
	
	// dense (SVD)
	A[1] = -0.9;
// 	LOG(blabla, "SVD start");
	solve_least_squares(X, A, B);
// 	LOG(blabla, "SVD end");
	MTEST(frob_norm(A(i,j)*X(j)-B(i)) < 1e-10, "5 " << frob_norm(A(i,j)*X(j)-B(i)));
	
});

129
static misc::UnitTest tensor_solve_sparse("Tensor", "solve_sparse", [](){
130
	std::mt19937_64 &rnd = xerus::misc::randomEngine;
131
	std::normal_distribution<double> dist(0.0, 1.0);
132
	const size_t N = 100;
Sebastian Wolf's avatar
Sebastian Wolf committed
133
	std::uniform_int_distribution<size_t> eDist(1, N*N-1);
134
135
136
137
138
139
140
141
142
143
144
	
	Index i,j,k;
	
	Tensor id = Tensor::identity({N,N});
	Tensor r = Tensor({N}, [](size_t _i)->value_t{return double(_i);});
	Tensor x;
	x(i) = r(j) / id(j,i);
	MTEST(frob_norm(x-r) < 1e-14, "d " << frob_norm(x-r));
	
	r.use_sparse_representation();
	x(i) = r(j) / id(j,i);
145
146
	MTEST(frob_norm(x-r) < 1e-14, "d " << frob_norm(x-r));
	r.use_dense_representation();
147
148
149
150
151
152
153
154
155
156
157
	
	// consistency with dense solve:
	for (size_t n=0; n<N*3; ++n) {
		id[eDist(rnd)] = dist(rnd);
	}
	id.use_sparse_representation();
	
	// test faithful reconstruction
	internal::CholmodSparse idt(id.get_unsanitized_sparse_data(), N, N, false);
	Tensor id2({N,N});
	id2.get_unsanitized_sparse_data() = idt.to_map();
Sebastian Wolf's avatar
Sebastian Wolf committed
158
	MTEST(frob_norm(id-id2) < 1e-12, frob_norm(id-id2)); 
159
160
161
162
163
164
165
166
	
	Tensor fid(id);
	fid.use_dense_representation();
	Tensor fx;
	TEST(id.is_sparse());
	
	fx(i) = r(j) / fid(j,i);
	x(i) = r(j) / id(j,i);
Sebastian Wolf's avatar
Sebastian Wolf committed
167
168
169
	MTEST(frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x)<1e-12, frob_norm(id(j,i)*x(i) - r(j))/frob_norm(x));
	MTEST(frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x)<1e-12, frob_norm(fid(j,i)*fx(i) - r(j))/frob_norm(x));
	MTEST(frob_norm(fx-x)/frob_norm(x)<1e-12, frob_norm(fx-x)/frob_norm(x));
170
});
171

172
static misc::UnitTest tensor_solve_trans("Tensor", "solve_transposed", [](){
173
	std::mt19937_64 &rnd = xerus::misc::randomEngine;
174
175
176
177
178
179
180
181
182
183
184
185
186
187
	std::normal_distribution<double> dist(0.0, 1.0);
	const size_t N = 100;
	std::uniform_int_distribution<size_t> eDist(1,N*N-1);
	
	Index i,j,k;
	
	Tensor A = Tensor::identity({N,N});
	for (size_t n=0; n<N*3; ++n) {
		A[eDist(rnd)] = dist(rnd);
	}
	A.use_sparse_representation();
	Tensor At;
	At(i,j) = A(j,i);
	At.use_sparse_representation();
188
	
189
190
191
192
	Tensor r = Tensor({N}, [](size_t _i)->value_t{return double(_i);});
	Tensor x1, x2;
	x1(i) = r(j) / A(i,j);
	x2(i) = r(j) / At(j,i);
Sebastian Wolf's avatar
Sebastian Wolf committed
193
	MTEST(frob_norm(x1-x2) < 1e-12, "s " << frob_norm(x1-x2));
194
195
196
197
	
	A.use_dense_representation();
	At.use_dense_representation();
	
198
	Tensor x3, x4, residual;
199
200
	x3(i) = r(j) / A(i,j);
	x4(i) = r(j) / At(j,i);
201
202
	
	residual(i) = A(i,j) * (x3(j) - x4(j));
Sebastian Wolf's avatar
Sebastian Wolf committed
203
	MTEST(frob_norm(x3-x4) < 1e-12, "d " << frob_norm(x3-x4) << " residual: " << frob_norm(residual));
204
205
	
	residual(i) = A(i,j) * (x1(j) - x3(j));
Sebastian Wolf's avatar
Sebastian Wolf committed
206
	MTEST(frob_norm(x1-x3)/frob_norm(x1) < 1e-12, "sd " << frob_norm(x1-x3)/frob_norm(x1) << " residual: " << frob_norm(residual));
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
});


static misc::UnitTest tensor_solve_matrix("Tensor", "solve_matrix", [](){
	std::mt19937_64 &rnd = xerus::misc::randomEngine;
	std::uniform_int_distribution<size_t> nDist(1, 100);
	std::uniform_int_distribution<size_t> n2Dist(1, 10);
	std::uniform_int_distribution<size_t> dDist(1, 3);
	std::uniform_int_distribution<size_t> d2Dist(0, 3);
	std::normal_distribution<double> realDist(0, 1);
	
	Index i,j,k;
	
	for(size_t run = 0; run < 10; ++run) {
		std::vector<size_t> mDims, nDims, pDims;
		const size_t degM = dDist(rnd);
		const size_t degN = dDist(rnd);
		const size_t degP = d2Dist(rnd);
		for(size_t xi = 0; xi < degM; ++xi) { mDims.push_back(n2Dist(rnd)); }
		for(size_t xi = 0; xi < degN; ++xi) { nDims.push_back(n2Dist(rnd)); }
		for(size_t xi = 0; xi < degP; ++xi) { pDims.push_back(n2Dist(rnd)); }
		
		
		auto A = Tensor::random(mDims | nDims);
		A *= realDist(rnd);
		Tensor B;
		Tensor realX = Tensor::random(nDims | pDims);
		
		B(i^degM, k^degP) = A(i^degM, j^degN)*realX(j^degN, k^degP);
		
		auto factor = realDist(rnd);
		B *= factor;
		realX *= factor;
		
		Tensor X;
		
		solve_least_squares(X, A, B, degP);
		
		Tensor residual;
		
		residual(i^degM, k^degP) = A(i^degM, j^degN)*X(j^degN, k^degP) - B(i^degM, k^degP);
Benjamin Huber's avatar
Benjamin Huber committed
248
249
		MTEST(frob_norm(residual) < 1e-10, frob_norm(residual));
		
250
		
Benjamin Huber's avatar
Benjamin Huber committed
251
252
253
		X(j^degN, k^degP) = B(i^degM, k^degP) / A(i^degM, j^degN);  //solve_least_squares(X, A, B, degP);
		
		residual(i^degM, k^degP) = A(i^degM, j^degN)*X(j^degN, k^degP) - B(i^degM, k^degP);
254
255
		MTEST(frob_norm(residual) < 1e-10, frob_norm(residual));
	}
256
});