sparse SVD via two QC decompositions -> SVD of core (closes #32)

src/unitTests/fullTensor_factorisations.cxx

@@ -332,3 +332,31 @@ static misc::UnitTest tensor_scq("Tensor", "Sparse_CQ", [](){
 // 	TEST(approx_equal(Cs, Cf, 1e-15));// NOTE apparently not true when there are small singular values
 	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");
 });
+
+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);
+	}
+});
+

src/xerus/blasLapackWrapper.cpp

@@ -318,8 +318,8 @@ namespace xerus {
 			REQUIRE(_n <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for BLAS/Lapack");
 			REQUIRE(_m <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for BLAS/Lapack");
 			
-			REQUIRE(_m > 0, "Dimension n must be larger than zero");
-			REQUIRE(_n > 0, "Dimension m must be larger than zero");
+			REQUIRE(_m > 0, "Dimension m must be larger than zero");
+			REQUIRE(_n > 0, "Dimension n must be larger than zero");
 			
 			XERUS_PA_START;
 			

src/xerus/cholmod_wrapper.cpp

@@ -220,10 +220,12 @@ namespace xerus { namespace internal {
 		REQUIRE(_m < std::numeric_limits<long>::max() && _n < std::numeric_limits<long>::max() && _A.size() < std::numeric_limits<long>::max(),
 			"sparse matrix given to qc decomposition too large for suitesparse"
 		);
+		REQUIRE(_m>0 && _n>0, "invalid matrix of dimensions " << _m << 'x' << _n);
 		CholmodSparse A(_A, _m, _n, _transposeA);
 		cholmod_sparse *Q, *R;
 		SuiteSparse_long *E;
 		long rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?long(std::min(_m,_n)):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
+		rank = std::max(rank, 1l);
 		CholmodSparse Qs(Q);
 		CholmodSparse Rs(R);
 		INTERNAL_CHECK(E == nullptr, "IE: sparse QR returned a permutation despite fixed ordering?!");
@@ -247,6 +249,7 @@ namespace xerus { namespace internal {
 		SuiteSparse_long *E;
 		// decompose A^T = q^T*r^T
 		long rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?long(std::min(_m,_n)):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
+		rank = std::max(rank, 1l);
 		CholmodSparse Qs(Q);
 		CholmodSparse Rs(R);
 		INTERNAL_CHECK(E == nullptr, "IE: sparse QR returned a permutation despite fixed ordering?!");

src/xerus/tensor.cpp

@@ -1426,16 +1426,37 @@ namespace xerus {
 		
 		size_t lhsSize, rhsSize, rank;
 		std::tie(lhsSize, rhsSize, rank) = calculate_factorization_sizes(_input, _splitPos);
-		prepare_factorization_output(_U, _Vt, _input, _splitPos, rank, Tensor::Representation::Dense);
 		
 		std::unique_ptr<value_t[]> tmpS(new value_t[rank]);
 		
+		// sparse SVD becomes inefficient when the matrix is not sparse enough
+		// sparse SVD is about equally fast to dense SVD when there are about N = 1.55*(min(m,n)+(max-min)/5) entries set
+		// will calculate with 2 instead of 1.55 to make sure that we will certainly be slower with sparse
+		// (note that the algorithm is quadratic in the sparsity)
+		size_t min = std::min(lhsSize, rhsSize);
+		size_t max = std::max(lhsSize, rhsSize);
+		if (_input.is_sparse() && _input.sparsity() > 2*(min+(max-min)/5)) {
+			_input.use_dense_representation();
+		}
+		
 		// Calculate the actual SVD
 		if(_input.is_sparse()) {
-			LOG_ONCE(warning, "Sparse SVD not yet implemented. falling back to the dense SVD");
-			_input.use_dense_representation();
-			blasWrapper::svd(_U.override_dense_data(), tmpS.get(), _Vt.override_dense_data(), _input.get_unsanitized_dense_data(), lhsSize, rhsSize);
+			// calculate QC in both directions
+			calculate_qc(_U, _input, _input, _splitPos);
+			calculate_cq(_input, _Vt, _input, 1);
+			_input.use_dense_representation(); // in the very unlikely case that is it not already dense
+			
+			// then calculate SVD only of remaining (hopefully small) core
+			Tensor UPrime, VPrime;
+			std::tie(lhsSize, rhsSize, rank) = calculate_factorization_sizes(_input, 1);
+			prepare_factorization_output(UPrime, VPrime, _input, 1, rank, Tensor::Representation::Dense);
+			blasWrapper::svd(UPrime.override_dense_data(), tmpS.get(), VPrime.override_dense_data(), _input.get_unsanitized_dense_data(), lhsSize, rhsSize);
+			
+			// contract U*UPrime and VPrime*V to obtain SVD (UU', S, V'V) from orthogonal U and V as wel as the SVD (U', S, V')
+			contract(_U, _U, UPrime, 1);
+			contract(_Vt, VPrime, _Vt, 1);
 		} else {
+			prepare_factorization_output(_U, _Vt, _input, _splitPos, rank, Tensor::Representation::Dense);
 			blasWrapper::svd(_U.override_dense_data(), tmpS.get(), _Vt.override_dense_data(), _input.get_unsanitized_dense_data(), lhsSize, rhsSize);
 		}