1 #include <Eigen/Core>
2 #include <Eigen/SparseCore>
3 #include <Eigen/SVD>
4 #include <iostream>
5 #include <random> // Requires C++ 11
6
7 #include <Spectra/contrib/PartialSVDSolver.h>
8
9 using namespace Spectra;
10
11 #include "catch.hpp"
12
13 using Matrix = Eigen::MatrixXd;
14 using Vector = Eigen::VectorXd;
15 using SpMatrix = Eigen::SparseMatrix<double>;
16
17 // Generate random sparse matrix
gen_sparse_data(int m,int n,double prob=0.5)18 SpMatrix gen_sparse_data(int m, int n, double prob = 0.5)
19 {
20 SpMatrix mat(m, n);
21 std::default_random_engine gen;
22 gen.seed(0);
23 std::uniform_real_distribution<double> distr(0.0, 1.0);
24 for (int i = 0; i < m; i++)
25 {
26 for (int j = 0; j < n; j++)
27 {
28 if (distr(gen) < prob)
29 mat.insert(i, j) = distr(gen) - 0.5;
30 }
31 }
32 return mat;
33 }
34
35 template <typename MatType>
run_test(const MatType & mat,int k,int m)36 void run_test(const MatType& mat, int k, int m)
37 {
38 PartialSVDSolver<MatType> svds(mat, k, m);
39 int nconv = svds.compute();
40
41 INFO("nconv = " << nconv);
42 REQUIRE(nconv == k);
43
44 Vector svals = svds.singular_values();
45 Matrix U = svds.matrix_U(k);
46 Matrix V = svds.matrix_V(k);
47
48 // SVD solver from Eigen
49 // Requires dense matrices
50 Matrix mat_dense = Matrix(mat);
51 Eigen::JacobiSVD<Matrix> svd(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
52 Vector svals_eigen = svd.singularValues();
53 Matrix U_eigen = svd.matrixU();
54 Matrix V_eigen = svd.matrixV();
55
56 double err = (svals - svals_eigen.head(k)).array().abs().maxCoeff();
57 INFO("Residual of singular values = " << err);
58 REQUIRE(err == Approx(0.0).margin(1e-9));
59
60 err = (U.array().abs() - U_eigen.leftCols(k).array().abs()).abs().maxCoeff();
61 INFO("Residual of left singular vectors = " << err);
62 REQUIRE(err == Approx(0.0).margin(1e-9));
63
64 err = (V.array().abs() - V_eigen.leftCols(k).array().abs()).abs().maxCoeff();
65 INFO("Residual of right singular vectors = " << err);
66 REQUIRE(err == Approx(0.0).margin(1e-9));
67 }
68
69 TEST_CASE("Partial SVD of tall dense matrix [1000x100]", "[svds_dense_tall]")
70 {
71 std::srand(123);
72
73 const Matrix A = Matrix::Random(1000, 100);
74 int k = 5;
75 int m = 10;
76
77 run_test<Matrix>(A, k, m);
78 }
79
80 TEST_CASE("Partial SVD of wide dense matrix [1000x100]", "[svds_dense_wide]")
81 {
82 std::srand(123);
83
84 const Matrix A = Matrix::Random(100, 1000);
85 int k = 5;
86 int m = 10;
87
88 run_test<Matrix>(A, k, m);
89 }
90
91 TEST_CASE("Partial SVD of tall sparse matrix [1000x100]", "[svds_sparse_tall]")
92 {
93 std::srand(123);
94
95 const SpMatrix A = gen_sparse_data(1000, 100, 0.1);
96 int k = 5;
97 int m = 10;
98
99 run_test<SpMatrix>(A, k, m);
100 }
101
102 TEST_CASE("Partial SVD of wide sparse matrix [1000x100]", "[svds_sparse_wide]")
103 {
104 std::srand(123);
105
106 const SpMatrix A = gen_sparse_data(100, 1000, 0.1);
107 int k = 5;
108 int m = 10;
109
110 run_test<SpMatrix>(A, k, m);
111 }
112