From c3c64dd164ece9c683247e5598dc39695518c863 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Wed, 6 Nov 2024 20:25:26 -0500 Subject: [PATCH 1/6] WIP --- include/sparseir/_linalg.hpp | 8 ++++++-- test/_linalg.cxx | 39 +++++++++++++++++++++--------------- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp index 53d1f07..028d89b 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; + 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; diff --git a/test/_linalg.cxx b/test/_linalg.cxx index 481605b..d6365c0 100644 --- a/test/_linalg.cxx +++ b/test/_linalg.cxx @@ -266,14 +266,19 @@ 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}) { + //auto N1 = 201; + //auto N2 = 51; + auto N1 = 20; + auto N2 = 20; + 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; @@ -284,15 +289,17 @@ TEST_CASE("TSVD", "[linalg]") { 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...? - REQUIRE(0==1); - //REQUIRE((U * S_diag * V).isApprox(Aorig, tol * A.norm())); + //std::cout << "U " << U.rows() << " " << U.cols() << std::endl; + //std::cout << "V " << V.rows() << " " << V.cols() << std::endl; + //std::cout << "S_diag " << S_diag.rows() << " " << S_diag.cols() << std::endl; + 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; + + REQUIRE(Areconst.isApprox(Aorig, tol * Aorig.norm())); + //a /* REQUIRE((U.transpose() * U).isIdentity()); REQUIRE((V.transpose() * V).isIdentity()); From 292c424838394d8b2146eb14bdddbfabc18a49e1 Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Wed, 6 Nov 2024 20:26:11 -0500 Subject: [PATCH 2/6] WIP --- test/_linalg.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/test/_linalg.cxx b/test/_linalg.cxx index d6365c0..9e3e04b 100644 --- a/test/_linalg.cxx +++ b/test/_linalg.cxx @@ -297,6 +297,7 @@ TEST_CASE("TSVD", "[linalg]") { 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())); //a From 8df37061ec99cf520049f0d9db78bfc45d3044f1 Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Thu, 7 Nov 2024 15:45:19 +0900 Subject: [PATCH 3/6] Fix bug --- include/sparseir/_linalg.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp index 028d89b..441606f 100644 --- a/include/sparseir/_linalg.hpp +++ b/include/sparseir/_linalg.hpp @@ -518,7 +518,7 @@ tsvd(const Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { svd.compute(R_trunc.transpose(), Eigen::ComputeThinU | Eigen::ComputeThinV); Eigen::PermutationMatrix perm(p.size()); - perm.indices() = invperm(p); + perm.indices() = p; Eigen::MatrixX U = Q_trunc * svd.matrixV(); // implement invperm From b4ff34e7c2d29dfa8fa7db98f1768b7bdb7eeffc Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Thu, 7 Nov 2024 15:45:27 +0900 Subject: [PATCH 4/6] Add comments --- include/sparseir/_linalg.hpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/include/sparseir/_linalg.hpp b/include/sparseir/_linalg.hpp index 441606f..1e2d0c7 100644 --- a/include/sparseir/_linalg.hpp +++ b/include/sparseir/_linalg.hpp @@ -498,7 +498,7 @@ tsvd(const Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { // Step 1: Apply RRQR to A QRPivoted A_qr; int k; - Eigen::MatrixX A_ = A; + 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); @@ -515,8 +515,13 @@ tsvd(const 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() = p; @@ -525,6 +530,7 @@ tsvd(const Eigen::MatrixX& A, T rtol = std::numeric_limits::epsilon()) { 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); } From 8c4a42902f04fc748665abe61cf53deaf906182c Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Thu, 7 Nov 2024 15:45:56 +0900 Subject: [PATCH 5/6] Complete tests :tada: --- test/_linalg.cxx | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/test/_linalg.cxx b/test/_linalg.cxx index 9e3e04b..680d939 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())); } @@ -267,10 +270,8 @@ TEST_CASE("TSVD", "[linalg]") { // double { for (auto tol : {1e-14}) { - //auto N1 = 201; - //auto N2 = 51; - auto N1 = 20; - auto N2 = 20; + int N1 = 201; + int N2 = 51; VectorX x = VectorX::LinSpaced(N1, -1, 1); //MatrixX Aorig(201, 51); MatrixX Aorig(N1, N2); @@ -281,17 +282,21 @@ TEST_CASE("TSVD", "[linalg]") { //} } - 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(); - //std::cout << "U " << U.rows() << " " << U.cols() << std::endl; - //std::cout << "V " << V.rows() << " " << V.cols() << std::endl; - //std::cout << "S_diag " << S_diag.rows() << " " << S_diag.cols() << std::endl; + + std::cout << "U " << U.rows() << " " << U.cols() << std::endl; + std::cout << "V " << V.rows() << " " << V.cols() << std::endl; + std::cout << "S_diag " << S_diag.rows() << " " << S_diag.cols() << std::endl; + auto Areconst = U * S_diag * V.transpose(); auto diff = (A - Areconst).norm() / A.norm(); std::cout << "diff " << diff << std::endl; @@ -300,16 +305,15 @@ TEST_CASE("TSVD", "[linalg]") { std::cout << "norm diff" << Aorig.norm() - Areconst.norm() << std::endl; REQUIRE(Areconst.isApprox(Aorig, tol * Aorig.norm())); - //a - /* + 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())); } } } @@ -389,7 +393,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; From 0a528bdf6842aebcf22e26162ad8fc716c63aa7a Mon Sep 17 00:00:00 2001 From: Satoshi Terasaki Date: Thu, 7 Nov 2024 15:47:54 +0900 Subject: [PATCH 6/6] Suppress output --- test/_linalg.cxx | 12 ++++-------- test/freq.cxx | 2 +- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/test/_linalg.cxx b/test/_linalg.cxx index 680d939..7b0a8e9 100644 --- a/test/_linalg.cxx +++ b/test/_linalg.cxx @@ -293,16 +293,12 @@ TEST_CASE("TSVD", "[linalg]") { auto S_diag = s.asDiagonal(); - std::cout << "U " << U.rows() << " " << U.cols() << std::endl; - std::cout << "V " << V.rows() << " " << V.cols() << std::endl; - std::cout << "S_diag " << S_diag.rows() << " " << S_diag.cols() << std::endl; - 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; + // 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())); 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); }