Fix bugs in tsvd

include/sparseir/_linalg.hpp

@@ -266,6 +266,9 @@ void reflectorApply(Eigen::VectorBlock<Eigen::Block<Eigen::MatrixX<T>, -1, 1, tr
     }
 }
 
+/*
+A will be modified in-place.
+*/
 template <typename T>
 std::pair<QRPivoted<T>, int> rrqr(MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
     using std::abs;
@@ -491,11 +494,12 @@ truncateQRResult(const Eigen::MatrixX<T>& Q, const Eigen::MatrixX<T>& R, int k)
 // Truncated SVD (TSVD)
 template <typename T>
 std::tuple<Eigen::MatrixX<T>, Eigen::VectorX<T>, Eigen::MatrixX<T>>
-tsvd(Eigen::MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
+tsvd(const Eigen::MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
     // Step 1: Apply RRQR to A
     QRPivoted<T> A_qr;
     int k;
-    std::tie(A_qr, k) = rrqr<T>(A, rtol);
+    Eigen::MatrixX<T> A_ = A; // create a copy of A
+    std::tie(A_qr, k) = rrqr<T>(A_, rtol);
     // Step 2: Truncate QR Result to rank k
     auto tqr = truncate_qr_result<T>(A_qr, k);
     auto p = A_qr.jpvt;
@@ -511,16 +515,22 @@ tsvd(Eigen::MatrixX<T>& A, T rtol = std::numeric_limits<T>::epsilon()) {
     // Do not use the svd_jacobi function directly.
     // Better to write a wrrapper function for the SVD.
     Eigen::JacobiSVD<decltype(R_trunc)> 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<Dynamic, Dynamic> perm(p.size());
-    perm.indices() = invperm(p);
+    perm.indices() = p;
 
     Eigen::MatrixX<T> U = Q_trunc * svd.matrixV();
     // implement invperm
     Eigen::MatrixX<T> V = (perm * svd.matrixU());
 
     Eigen::VectorX<T> s = svd.singularValues();
+    // TODO: Create a return type for truncated SVD
     return std::make_tuple(U, s, V);
 }
 

test/_linalg.cxx

@@ -256,7 +256,10 @@ TEST_CASE("RRQR Trunc", "[linalg]") {
                 auto QR = truncate_qr_result<DDouble>(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<DDouble> 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<double> x = VectorX<double>::LinSpaced(201, -1, 1);
-                MatrixX<double> Aorig(201, 51);
+            for (auto tol : {1e-14}) {
+                int N1 = 201;
+                int N2 = 51;
+                VectorX<double> x = VectorX<double>::LinSpaced(N1, -1, 1);
+                //MatrixX<double> Aorig(201, 51);
+                MatrixX<double> 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<double> A = Aorig;
+                MatrixX<double> A = Aorig; // create a copy of Aorig
+
                 auto tsvd_result = tsvd<double>(A, double(tol));
                 tsvd<double>(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<DDouble>()));
+                REQUIRE(std::is_sorted(s.data(), s.data() + s.size(), std::greater<DDouble>()));
                 REQUIRE(k < std::min(A.rows(), A.cols()));
 
-                Eigen::JacobiSVD<MatrixX<DDouble>> svd(A.cast<DDouble>());
-                REQUIRE(S.isApprox(svd.singularValues().head(k).cast<DDouble>()));
-                */
+                Eigen::JacobiSVD<MatrixX<double>> svd(Aorig.cast<double>());
+                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;

test/freq.cxx

@@ -52,7 +52,7 @@ TEST_CASE("freq", "[Imaginary value calculation for MatsubaraFreq]")
     double beta = 3.0;
     BosonicFreq bf(2);
     std::complex<double> 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);
 }