diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp index 53d1f07..1e2d0c7 100644 --- a/include/sparseir/_linalg.hpp +++ b/include/sparseir/_linalg.hpp @@ -266,6 +266,9 @@ void reflectorApply(Eigen::VectorBlock, -1, 1, tr } } +/* +A will be modified in-place. +*/ template std::pair, int> rrqr(MatrixX& A, T rtol = std::numeric_limits::epsilon()) { using std::abs; @@ -491,11 +494,12 @@ truncateQRResult(const Eigen::MatrixX& Q, const Eigen::MatrixX& R, int k) // Truncated SVD (TSVD) template std::tuple, Eigen::VectorX, Eigen::MatrixX> -tsvd(Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { +tsvd(const Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { // Step 1: Apply RRQR to A QRPivoted A_qr; int k; - std::tie(A_qr, k) = rrqr(A, rtol); + Eigen::MatrixX A_ = A; // create a copy of A + std::tie(A_qr, k) = rrqr(A_, rtol); // Step 2: Truncate QR Result to rank k auto tqr = truncate_qr_result(A_qr, k); auto p = A_qr.jpvt; @@ -511,16 +515,22 @@ tsvd(Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { // Do not use the svd_jacobi function directly. // Better to write a wrrapper function for the SVD. Eigen::JacobiSVD svd; + + // The following comment is taken from Julia's implementation + // # RRQR is an excellent preconditioner for Jacobi. One should then perform + // # Jacobi on RT svd.compute(R_trunc.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); + // Reconstruct A from QR factorization Eigen::PermutationMatrix perm(p.size()); - perm.indices() = invperm(p); + perm.indices() = p; Eigen::MatrixX U = Q_trunc * svd.matrixV(); // implement invperm Eigen::MatrixX V = (perm * svd.matrixU()); Eigen::VectorX s = svd.singularValues(); + // TODO: Create a return type for truncated SVD return std::make_tuple(U, s, V); } diff --git a/test/_linalg.cxx b/test/_linalg.cxx index d7b1eab..70fd63d 100644 --- a/test/_linalg.cxx +++ b/test/_linalg.cxx @@ -256,7 +256,10 @@ TEST_CASE("RRQR Trunc", "[linalg]") { auto QR = truncate_qr_result(A_qr, k); auto Q = QR.first; auto R = QR.second; - + REQUIRE(Q.rows() == m); + REQUIRE(Q.cols() == k); + REQUIRE(R.rows() == k); + REQUIRE(R.cols() == n); MatrixX A_rec = Q * R * getPropertyP(A_qr).transpose(); REQUIRE(A_rec.isApprox(Aorig, 1e-5 * A.norm())); } @@ -266,43 +269,45 @@ TEST_CASE("TSVD", "[linalg]") { using std::pow; // double { - for (auto tol : {1e-14, 1e-13}) { - VectorX x = VectorX::LinSpaced(201, -1, 1); - MatrixX Aorig(201, 51); + for (auto tol : {1e-14}) { + int N1 = 201; + int N2 = 51; + VectorX x = VectorX::LinSpaced(N1, -1, 1); + //MatrixX Aorig(201, 51); + MatrixX Aorig(N1, N2); for (int i = 0; i < Aorig.cols(); i++) { - //Aorig.col(i) = x.array().pow(i); - for (int j = 0; j < Aorig.rows(); j++) { - Aorig(j, i) = pow(x(j), i); - } + Aorig.col(i) = x.array().pow(i); + //for (int j = 0; j < Aorig.rows(); j++) { + //Aorig(j, i) = pow(x(j), i); + //} } - MatrixX A = Aorig; + MatrixX A = Aorig; // create a copy of Aorig + auto tsvd_result = tsvd(A, double(tol)); tsvd(Aorig, double(tol)); auto U = std::get<0>(tsvd_result); auto s = std::get<1>(tsvd_result); auto V = std::get<2>(tsvd_result); int k = s.size(); + auto S_diag = s.asDiagonal(); - // U * S_diag * V.transpose(); - // std::cout << "U: " << U.rows() << "," << U.cols() << std::endl; - // std::cout << "S: " << S_diag.rows() << "," << S_diag.cols() << std::endl; - // std::cout << "V: " << V.rows() << "," << V.cols() << std::endl; - auto B = U * S_diag * V.transpose() - Aorig; - std::cout << Aorig.norm() << std::endl; // Oh...? - std::cout << B.norm() << std::endl; // Oh...? - bool test_completed = false; - REQUIRE(!test_completed); - //REQUIRE((U * S_diag * V).isApprox(Aorig, tol * A.norm())); - /* + auto Areconst = U * S_diag * V.transpose(); + auto diff = (A - Areconst).norm() / A.norm(); + // std::cout << "diff " << diff << std::endl; + // std::cout << "Areconst " << Areconst.norm() << std::endl; + // std::cout << "Aorig " << Aorig.norm() << std::endl; + // std::cout << "norm diff" << Aorig.norm() - Areconst.norm() << std::endl; + + REQUIRE(Areconst.isApprox(Aorig, tol * Aorig.norm())); REQUIRE((U.transpose() * U).isIdentity()); REQUIRE((V.transpose() * V).isIdentity()); - REQUIRE(std::is_sorted(S.data(), S.data() + S.size(), std::greater())); + REQUIRE(std::is_sorted(s.data(), s.data() + s.size(), std::greater())); REQUIRE(k < std::min(A.rows(), A.cols())); - Eigen::JacobiSVD> svd(A.cast()); - REQUIRE(S.isApprox(svd.singularValues().head(k).cast())); - */ + Eigen::JacobiSVD> svd(Aorig.cast()); + REQUIRE(s.isApprox(svd.singularValues().head(k))); + REQUIRE(S_diag.toDenseMatrix().isApprox(svd.singularValues().head(k).asDiagonal().toDenseMatrix())); } } } @@ -382,7 +387,7 @@ TEST_CASE("SVD of VERY triangular 2x2", "[linalg]") { REQUIRE(cu.hi() == 1.0); REQUIRE(su.hi() == 1e-100); REQUIRE(smax.hi() == 1.0); - REQUIRE(smin.hi() == 1e-200); // so cloe + REQUIRE(smin.hi() == 1e-200); REQUIRE(cv.hi() == 1e-100); REQUIRE(sv.hi() == 1.0); U << cu, -su, su, cu; diff --git a/test/freq.cxx b/test/freq.cxx index 5b01aaf..8113880 100644 --- a/test/freq.cxx +++ b/test/freq.cxx @@ -52,7 +52,7 @@ TEST_CASE("freq", "[Imaginary value calculation for MatsubaraFreq]") double beta = 3.0; BosonicFreq bf(2); std::complex expected_value_im(0, 2 * M_PI / beta); - std::cout << bf.value_im(beta) << std::endl; + // std::cout << bf.value_im(beta) << std::endl; REQUIRE(bf.value_im(beta) == expected_value_im); }