fullTensor_solve.cxx 7.28 KB
Newer Older
Baum's avatar
Baum committed
1
// Xerus - A General Purpose Tensor Library
2
// Copyright (C) 2014-2017 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
});