fullTensor_factorisations.cxx 15.2 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
using namespace xerus;
Baum's avatar
Baum committed
25

26
static misc::UnitTest tensor_svd_id("Tensor", "SVD_Identity", [](){
Sebastian Wolf's avatar
Sebastian Wolf committed
27
28
29
30
    Tensor A({2,2,2,2});
    Tensor res1({2,2,4});
    Tensor res2({4,4});
    Tensor res3({4,2,2});
Baum's avatar
Baum committed
31
32
33
34
35
36
37
38
39
    
    A[{0,0,0,0}] = 1;
    A[{0,1,0,1}] = 1;
    A[{1,0,1,0}] = 1;
    A[{1,1,1,1}] = 1;    
    
    Index i, j, k, l, m, n;
    
    (res1(i,j,m), res2(m,n), res3(n,k,l)) = SVD(A(i,j,k,l));
Sebastian Wolf's avatar
Sebastian Wolf committed
40
41
42
	res2.reinterpret_dimensions({2,2,2,2});
	TEST(approx_entrywise_equal(res2, A));
	res2.reinterpret_dimensions({4,4});
43
44
	MTEST(frob_norm(res1(i^2, m)*res1(i^2, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(m, i^2)*res3(n, i^2) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
45
	
Sebastian Wolf's avatar
Sebastian Wolf committed
46
47
48
49
	(res1(m,i,j), res2(n,m), res3(k,n,l)) = SVD(A(i,j,k,l));
	res2.reinterpret_dimensions({2,2,2,2});
	TEST(approx_entrywise_equal(res2, A));
	res2.reinterpret_dimensions({4,4});
50
51
	MTEST(frob_norm(res1(m, i^2)*res1(n, i^2) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(k,m,l)*res3(k,n,l) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
52
});
Baum's avatar
Baum committed
53

54
static misc::UnitTest tensor_svd_zero("Tensor", "SVD_zero", [](){
Benjamin Huber's avatar
Benjamin Huber committed
55
56
57
58
59
60
61
62
63
64
65
66
    Tensor A({2,2,2,2});
    Tensor res1({2,2,4});
    Tensor res2({4,4});
    Tensor res3({4,2,2});
    
    Index i, j, k, l, m, n;
    
    (res1(i,j,m), res2(m,n), res3(n,k,l)) = SVD(A(i,j,k,l));
	MTEST(res2.dimensions[0] == 1, res2.dimensions[0]);
	MTEST(std::abs(res2[0]) < 1e-15, res2[0]);
	MTEST(frob_norm(res1(i^2, m)*res1(i^2, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(m, i^2)*res3(n, i^2) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
67
});
Benjamin Huber's avatar
Benjamin Huber committed
68

69
static misc::UnitTest tensor_svd_rnd512("Tensor", "SVD_Random_512x512", [](){
70
    Tensor A = Tensor::random({8,8,8,8,8,8});
Sebastian Wolf's avatar
Sebastian Wolf committed
71
72
73
74
    Tensor res1;
    Tensor res2;
    Tensor res3;
    Tensor res4;
Baum's avatar
Baum committed
75
76
77
78
79
80
     
    
    Index i, j, k, l, m, n, o, p, r, s;
    
    (res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(i,j,k,l,m,n));
    res4(i,j,k,l,m,n) = res1(i,j,k,o)*res2(o,p)*res3(p,l,m,n);
81
    TEST(approx_equal(res4, A, 1e-14));
82
83
	MTEST(frob_norm(res1(i^3, m)*res1(i^3, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(m, i^3)*res3(n, i^3) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
Baum's avatar
Baum committed
84
85
86
    
    (res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(l,k,m,i,j,n));
    res4(l,k,m,i,j,n) =  res1(i,j,k,o)*res2(o,p)*res3(p,l,m,n);
87
    TEST(approx_equal(res4, A, 1e-14));
88
89
	MTEST(frob_norm(res1(i^3, m)*res1(i^3, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(m, i^3)*res3(n, i^3) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
Baum's avatar
Baum committed
90
91
92
    
    (res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(l,i,m,k,j,n));
    res4(k,i,m,l,j,n) =  res1(i,j,l,o)*res2(o,p)*res3(p,k,m,n);
93
    TEST(approx_equal(res4, A, 1e-14));
94
95
	MTEST(frob_norm(res1(i^3, m)*res1(i^3, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(m, i^3)*res3(n, i^3) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
96
97
	
	(res1(i,o,k,j), res2(p,o), res3(l,n,m,p)) = SVD(A(l,i,m,k,j,n));
98
99
100
101
102
103
104
105
106
	res4(l,k,m,i,j,n) =  res1(k,o,i,j)*res2(p,o)*res3(l,n,m,p);
	TEST(approx_equal(res4, A, 1e-14));
	MTEST(frob_norm(res1(k,o,i,j)*res1(k,p,i,j) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(l,n,m,o)*res3(l,n,m,p) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " Vt not orthogonal");
	
	
	(res1(i,o,k,j), res2(p,o), res3(l,n,m,p)) = SVD(-1*A(l,i,m,k,j,n));
	res4(l,k,m,i,j,n) =  res1(k,o,i,j)*res2(p,o)*res3(l,n,m,p);
	TEST(approx_equal(-1*res4, A, 1e-14));
107
108
	MTEST(frob_norm(res1(k,o,i,j)*res1(k,p,i,j) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(l,n,m,o)*res3(l,n,m,p) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " Vt not orthogonal");
109
});
Baum's avatar
Baum committed
110

111
static misc::UnitTest tensor_svd_soft("Tensor", "SVD_soft_thresholding", [](){
112
    Tensor A = 10*Tensor::random({3,5,2,7,3,12});
Sebastian Wolf's avatar
Sebastian Wolf committed
113
    Tensor Ax, U, V, Us, Vs;
Sebastian Wolf's avatar
Sebastian Wolf committed
114
115
    Tensor S(Tensor::Tensor::Representation::Sparse);
	Tensor Ss(Tensor::Tensor::Representation::Sparse);
116
117
118
119
    
    Index i, j, k, l, m, n, o, p, r, s;
    
    (U(i,j,k,o), S(o,p), V(p,l,m,n)) = SVD(A(i,j,k,l,m,n));
120
    (Us(i,j,k,o), Ss(o,p), Vs(p,l,m,n)) = SVD(A(i,j,k,l,m,n), 7.3);
121
	
122
123
	U.resize_mode(U.degree()-1, Ss.dimensions[0]);
	V.resize_mode(0, Ss.dimensions[0]);
124
125
126
127
128
129
130
131
132
133
134
135
136
137
	
    TEST(approx_equal(U, Us, 1e-12));
    TEST(approx_equal(V, Vs, 1e-12));
	
	for(size_t x = 0; x < S.dimensions[0]; ++x) {
		if(x < Ss.dimensions[0]) {
			TEST(misc::approx_equal(Ss[{x,x}], std::max(0.0, S[{x, x}]-7.3), 3e-13));
		} else {
			TEST(S[{x,x}] <= 7.3);
		}
	}
	
	Ax(i,j,k,l,m,n) = U(i,j,k,o)* S(o,p)* V(p,l,m,n);
	TEST(approx_equal(A, Ax, 1e-12));
138
139
	MTEST(frob_norm(U(i,j,k,o)*U(i,j,k,p) - Tensor::identity(S.dimensions)(o, p)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(V(o,l,m,n)*V(p,l,m,n) - Tensor::identity(S.dimensions)(o, p)) < 1e-12, " Vt not orthogonal");
140
});
141

142
static misc::UnitTest tensor_svd_order_6("Tensor", "SVD_Random_Order_Six", [](){
143
    Tensor A = Tensor::random({9,7,5,5,9,7});
Sebastian Wolf's avatar
Sebastian Wolf committed
144
145
146
147
    Tensor res1;
    Tensor res2;
    Tensor res3;
    Tensor res4;
Baum's avatar
Baum committed
148
149
150
151
152
153
     
    
    Index i, j, k, l, m, n, o, p, r, s;
    
    (res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(i,j,k,l,m,n));
    res4(i,j,k,l,m,n) = res1(i,j,k,o)*res2(o,p)*res3(p,l,m,n);
154
    TEST(approx_equal(res4, A, 1e-14));
155
156
	MTEST(frob_norm(res1(i,j,k,o)*res1(i,j,k,p) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(o,l,m,n)*res3(p,l,m,n) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " Vt not orthogonal");
Baum's avatar
Baum committed
157
158
159
    
    (res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(m,j,k,l,i,n));
    res4(m,j,k,l,i,n) = res1(i,j,k,o)*res2(o,p)*res3(p,l,m,n);
160
    TEST(approx_equal(res4, A, 1e-14));
161
162
	MTEST(frob_norm(res1(i,j,k,o)*res1(i,j,k,p) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " U not orthogonal");
	MTEST(frob_norm(res3(o,l,m,n)*res3(p,l,m,n) - Tensor::identity(res2.dimensions)(o, p)) < 1e-12, " Vt not orthogonal");
163
});
Baum's avatar
Baum committed
164

165
static misc::UnitTest tensor_qr_rq_rnd6("Tensor", "QR_AND_RQ_Random_Order_Six", [](){
166
    Tensor A = Tensor::random({7,5,9,7,5,9});
Sebastian Wolf's avatar
Sebastian Wolf committed
167
168
169
170
171
172
173
174
    Tensor Q;
    Tensor R;
    Tensor Q2;
    Tensor R2;
    Tensor Q3;
    Tensor R3;
    Tensor Q4;
    Tensor res4;
Baum's avatar
Baum committed
175
176
177
178
179
180
    
    Index i, j, k, l, m, n, o, p, q, r;

    
    (Q(i,j,k,l), R(l,m,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q(i,j,k,o)*R(o,m,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
181
    MTEST(approx_equal(res4, A, 2e-15), "1 " << frob_norm(res4-A) << " / " << frob_norm(A));
182
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
183
	
Benjamin Huber's avatar
Benjamin Huber committed
184
	res4(l,m) = Q(i,j,k,l) * Q(i,j,k,m);
185
	res4.modify_diagonal_entries([](value_t &entry){entry -= 1;});
186
	TEST(frob_norm(res4) < 1e-12);
Baum's avatar
Baum committed
187
188
189
    
    (Q(i,j,k,l), R(l,m,n,r)) = QR(A(i,n,k,m,j,r));
    res4(i,n,k,m,j,r) = Q(i,j,k,o)*R(o,m,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
190
    MTEST(approx_equal(res4, A, 2e-15), "2 " << frob_norm(res4-A) << " / " << frob_norm(A));
191
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
Baum's avatar
Baum committed
192
193
194
    
    (Q2(i,k,l), R2(l,m,j,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q2(i,k,o)*R2(o,m,j,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
195
    MTEST(approx_equal(res4, A, 2e-15), "3 " << frob_norm(res4-A) << " / " << frob_norm(A));
Baum's avatar
Baum committed
196
197
198
    
    (Q3(i,m,j,k,l), R3(l,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q3(i,m,j,k,o)*R3(o,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
199
    MTEST(approx_equal(res4, A, 2e-15), "4 " << frob_norm(res4-A) << " / " << frob_norm(A));
200
201
202
	
	(Q(i,l,j,k,m), R(l,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q(i,o,j,k,m)*R(o,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
203
    MTEST(approx_equal(res4, A, 2e-15), "5 " << frob_norm(res4-A) << " / " << frob_norm(A));
204
	MTEST(frob_norm(Q(i,l,j,k,m)*Q(i,q,j,k,m) - Tensor::identity({Q.dimensions[1], Q.dimensions[1]})(l, q)) < 1e-12, " Q not orthogonal");
Baum's avatar
Baum committed
205
206
207
208
    
    
    (R(i,j,k,l), Q(l,m,n,r)) = RQ(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = R(i,j,k,o)*Q(o,m,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
209
    MTEST(approx_equal(res4, A, 5e-15), "6 " << frob_norm(res4-A) << " / " << frob_norm(A));
210
	MTEST(frob_norm(Q(p,m,n,r)*Q(q,m,n,r) - Tensor::identity({Q.dimensions[0], Q.dimensions[0]})(p, q)) < 1e-12, " Q not orthogonal");
Baum's avatar
Baum committed
211
212
213
    
    (R(i,j,k,l), Q(l,m,n,r)) = RQ(A(i,n,k,m,j,r));
    res4(i,n,k,m,j,r) = R(i,j,k,o)*Q(o,m,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
214
    MTEST(approx_equal(res4, A, 5e-15), "7 " << frob_norm(res4-A) << " / " << frob_norm(A));
215
	MTEST(frob_norm(Q(p,m,n,r)*Q(q,m,n,r) - Tensor::identity({Q.dimensions[0], Q.dimensions[0]})(p, q)) < 1e-12, " Q not orthogonal");
Baum's avatar
Baum committed
216
217
218
    
    (R2(i,m,j,k,l), Q2(l,n,r)) = RQ(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = R2(i,m,j,k,l)*Q2(l,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
219
    MTEST(approx_equal(res4, A, 2e-15), "8 " << frob_norm(res4-A) << " / " << frob_norm(A));
Baum's avatar
Baum committed
220
221
222
    
    (R3(i,k,l), Q3(l,m,j,n,r)) = RQ(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = R3(i,k,o)*Q3(o,m,j,n,r);
Benjamin Huber's avatar
Benjamin Huber committed
223
    MTEST(approx_equal(res4, A, 2e-15), "9 " << frob_norm(res4-A) << " / " << frob_norm(A));
224
225
226
	
	(R(l,i,k), Q(n,m,j,l,r)) = RQ(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = R(o,i,k)*Q(n,m,j,o,r);
Benjamin Huber's avatar
Benjamin Huber committed
227
    MTEST(approx_equal(res4, A, 2e-15), "10 " << frob_norm(res4-A) << " / " << frob_norm(A));
228
	MTEST(frob_norm(Q(n,m,j,q,r)*Q(n,m,j,p,r) - Tensor::identity({Q.dimensions[3], Q.dimensions[3]})(p, q)) < 1e-12, " Q not orthogonal");
229
230
231
232
233
    
    
    (R3(i,k,l), Q4(l,j,n,r)) = RQ(A(i,j,k,3,n,r));
    res4(i,j,k,n,r) = R3(i,k,o)*Q4(o,j,n,r);
    TEST(frob_norm(A(i,j,k,3,n,r) - res4(i,j,k,n,r)) < 1e-12);
234
	MTEST(frob_norm(Q4(q,j,n,r)*Q4(p,j,n,r) - Tensor::identity({Q4.dimensions[0], Q4.dimensions[0]})(p, q)) < 1e-12, " Q not orthogonal");
235
});
Baum's avatar
Baum committed
236

Benjamin Huber's avatar
Benjamin Huber committed
237

238
static misc::UnitTest tensor_qc("Tensor", "QC", [](){
239
    Tensor A = Tensor::random({2,2,2,2,2,2});
Sebastian Wolf's avatar
Sebastian Wolf committed
240
241
242
243
244
245
246
247
248
	Tensor B({2,3}, [](size_t i){return double(i);});
    Tensor Q;
    Tensor R;
    Tensor Q2;
    Tensor R2;
    Tensor Q3;
    Tensor R3;
    Tensor Q4;
    Tensor res4;
Benjamin Huber's avatar
Benjamin Huber committed
249
250
    
    Index i, j, k, l, m, n, o, p, q, r;
251
	
252
	(Q(i,j), R(j,k)) = QC(B(i,k));
Benjamin Huber's avatar
Benjamin Huber committed
253
    
254
    (Q(i,j,k,l), R(l,m,n,r)) = QC(A(i,j,k,m,n,r));
Benjamin Huber's avatar
Benjamin Huber committed
255
    res4(i,j,k,m,n,r) = Q(i,j,k,o)*R(o,m,n,r);
256
    TEST(approx_equal(res4, A, 1e-15));
257
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
Benjamin Huber's avatar
Benjamin Huber committed
258
    
259
    (Q(i,j,k,l), R(l,m,n,r)) = QC(A(i,n,k,m,j,r));
Benjamin Huber's avatar
Benjamin Huber committed
260
    res4(i,n,k,m,j,r) = Q(i,j,k,o)*R(o,m,n,r);
261
    TEST(approx_equal(res4, A, 1e-15));
262
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
Benjamin Huber's avatar
Benjamin Huber committed
263
    
264
    (Q2(i,k,l), R2(l,m,j,n,r)) = QC(A(i,j,k,m,n,r));
Benjamin Huber's avatar
Benjamin Huber committed
265
266
267
    res4(i,j,k,m,n,r) = Q2(i,k,o)*R2(o,m,j,n,r);
    TEST(approx_equal(res4, A, 1e-12));
    
268
    (Q3(i,m,j,k,l), R3(l,n,r)) = QC(A(i,j,k,m,n,r));
Benjamin Huber's avatar
Benjamin Huber committed
269
    res4(i,j,k,m,n,r) = Q3(i,m,j,k,o)*R3(o,n,r);
270
    TEST(approx_equal(res4, A, 1e-15));
271
272
273
	
	(Q(i,l,j,k,m), R(l,n,r)) = QC(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q(i,o,j,k,m)*R(o,n,r);
274
    TEST(approx_equal(res4, A, 1e-15));
275
	MTEST(frob_norm(Q(i,l,j,k,m)*Q(i,q,j,k,m) - Tensor::identity({Q.dimensions[1], Q.dimensions[1]})(l, q)) < 1e-12, " Q not orthogonal");
276
});
277

278
static misc::UnitTest tensor_sqr("Tensor", "Sparse_QR", [](){
279
    Tensor A =  Tensor::random({2,2,2,2,2,2}, 16);
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
	A.use_sparse_representation();
	Tensor B({2,3}, [](size_t i){return double(i);});
    Tensor Q;
    Tensor R;
    Tensor Q2;
    Tensor R2;
    Tensor Q3;
    Tensor R3;
    Tensor Q4;
    Tensor res4;
    
    Index i, j, k, l, m, n, o, p, q, r;
	
    (Q(i,j,k,l), R(l,m,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q(i,j,k,o)*R(o,m,n,r);
    TEST(approx_equal(res4, A, 1e-15));
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
    
    (Q(i,j,k,l), R(l,m,n,r)) = QR(A(i,n,k,m,j,r));
    res4(i,n,k,m,j,r) = Q(i,j,k,o)*R(o,m,n,r);
    TEST(approx_equal(res4, A, 1e-15));
	MTEST(frob_norm(Q(i,j,k,l)*Q(i,j,k,m) - Tensor::identity({R.dimensions[0], R.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
    
    (Q2(i,k,l), R2(l,m,j,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q2(i,k,o)*R2(o,m,j,n,r);
    TEST(approx_equal(res4, A, 1e-12));
    
    (Q3(i,m,j,k,l), R3(l,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q3(i,m,j,k,o)*R3(o,n,r);
    TEST(approx_equal(res4, A, 1e-15));
	
	(Q(i,l,j,k,m), R(l,n,r)) = QR(A(i,j,k,m,n,r));
    res4(i,j,k,m,n,r) = Q(i,o,j,k,m)*R(o,n,r);
    TEST(approx_equal(res4, A, 1e-15));
	MTEST(frob_norm(Q(i,l,j,k,m)*Q(i,q,j,k,m) - Tensor::identity({Q.dimensions[1], Q.dimensions[1]})(l, q)) < 1e-12, " Q not orthogonal");
315
});
316
317


318
static misc::UnitTest tensor_scq("Tensor", "Sparse_CQ", [](){
319
    Tensor A = Tensor::random({3,3,3,3,3,3}, 16);
320
321
322
323
324
325
326
327
328
329
330
	A.use_sparse_representation();
	Tensor Af(A);
	Af.use_dense_representation();
	A *= 2; Af *= 2;
    Tensor Qs, Qf, Cs, Cf, res;
    Index i, j, k, l, m, n, o, p, q, r;
	
    (Cs(i,j,k,l), Qs(l,m,n,r)) = CQ(A(i,j,k,m,n,r));
	(Cf(i,j,k,l), Qf(l,m,n,r)) = CQ(Af(i,j,k,m,n,r));
    res(i,j,k,m,n,r) = Cs(i,j,k,o)*Qs(o,m,n,r);
    TEST(approx_equal(res, A, 1e-15));
Ben Huber's avatar
Ben Huber committed
331
332
// 	TEST(approx_equal(Qs, Qf, 1e-15)); // NOTE apparently not true when there are small singular values
// 	TEST(approx_equal(Cs, Cf, 1e-15));// NOTE apparently not true when there are small singular values
333
	MTEST(frob_norm(Qs(l,i,j,k)*Qs(m,i,j,k) - Tensor::identity({Qs.dimensions[0], Qs.dimensions[0]})(l, m)) < 1e-12, " Q not orthogonal");
334
});
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362

static misc::UnitTest sparse_svd("Tensor", "SparseSVD", [](){
	Index i,j,k,l;
	size_t m1=400, m2=400;
	for (size_t n=10; n<=1000; n+=100) {
		Tensor A = Tensor::random({m1, m2}, 1*n);
		
		Tensor U, S, Vt;
// 		uint64 start = xerus::misc::uTime();
		(U(i,j), S(j,k), Vt(k,l)) = SVD(A(i,l));
// 		uint64 tSparse = xerus::misc::uTime() - start;
// 		XERUS_LOG(sparsity, U.is_sparse() << ' ' << S.is_sparse() << ' ' << Vt.is_sparse());
		Tensor T;
		T(i,l) = U(i,j)*S(j,k)*Vt(k,l);
		MTEST(approx_equal(T, A, 1e-14), n << ' ' << frob_norm(T-A) << '/' << frob_norm(A));
		U(i,j)=U(k,i)*U(k,j) - Tensor::identity(S.dimensions)(i,j);
		MTEST(frob_norm(U)<1e-10, n << ' ' << frob_norm(U));
		Vt(i,j)=Vt(i,k)*Vt(j,k) - Tensor::identity(S.dimensions)(i,j);
		MTEST(frob_norm(Vt)<1e-10, n << ' ' << frob_norm(Vt));
		
// 		A.use_dense_representation();
// 		start = xerus::misc::uTime();
// 		(U(i,j), S(j,k), Vt(k,l)) = SVD(A(i,l));
// 		uint64 tDense = xerus::misc::uTime() - start;
// 		XERUS_LOG(times, n << ' ' << tSparse << ' ' << tDense);
	}
});