Written the deflation in cpp

Author: Gaspare99

CMakeLists.txt

@@ -7,6 +7,14 @@ project(QR_binding VERSION 1.0
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/src/pyclassify)
 set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/src/pyclassify)
 
+include(FetchContent)
+# Fetch Eigen 3.4.0 (or whatever version you need)
+FetchContent_Declare(
+  eigen
+  URL  https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
+)
+FetchContent_MakeAvailable(eigen)
+
 set(PYBIND11_FINDPYTHON ON)
 
 # To find pybind11 on Ulysses:
@@ -36,6 +44,8 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 add_compile_options(-O3 -Wall -Werror -Wpedantic)
 
 pybind11_add_module(cxx_utils ${CMAKE_SOURCE_DIR}/src/pyclassify/cxx_utils.cpp)
+
+target_link_libraries(cxx_utils   PRIVATE    Eigen3::Eigen  )
 set_target_properties(cxx_utils PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${CMAKE_LIBRARY_OUTPUT_DIRECTORY})
 
 target_compile_features(cxx_utils PUBLIC cxx_std_17)

src/pyclassify/cxx_utils.cpp

@@ -14,6 +14,7 @@
 
 
 namespace py=pybind11;	
+#include <pybind11/eigen.h>
 
 
 std::pair<std::vector<double>, std::vector<std::vector<double>> > 
@@ -585,6 +586,183 @@ secular_solver(
 }
 
 
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <unordered_set>
+
+/**
+ * Applies the deflation step in a divide-and-conquer eigenvalue algorithm.
+ *
+ * @param D          Diagonal entries of the matrix (as Eigen::VectorXd).
+ * @param v          Rank-one update vector (modified in-place).
+ * @param beta       Scalar multiplier for the rank-one update.
+ * @param tol_factor Factor to scale the deflation tolerance (default 1e-12).
+ * @return A tuple containing:
+ *     - deflated_eigvals: Vector of trivial eigenvalues (Eigen::VectorXd).
+ *     - deflated_eigvecs: Matrix whose columns are the trivial eigenvectors (Eigen::MatrixXd).
+ *     - D_keep: Remaining diagonal entries after deflation (Eigen::VectorXd).
+ *     - v_keep: Remaining rank-one vector entries after deflation (Eigen::VectorXd).
+ *     - P_final:   Combined permutation & Givens rotation as an Eigen::SparseMatrix<double>.
+ */
+
+
+/**
+ * Applies the deflation step in a divide-and-conquer eigenvalue algorithm.
+ *
+ * @param D          Diagonal entries of the matrix (as Eigen::VectorXd).
+ * @param v          Rank-one update vector (modified in-place).
+ * @param beta       Scalar multiplier for the rank-one update.
+ * @param tol_factor Factor to scale the deflation tolerance (default 1e-12).
+ * @return A tuple containing:
+ *     - deflated_eigvals: Vector of trivial eigenvalues (Eigen::VectorXd).
+ *     - deflated_eigvecs: Matrix whose columns are the trivial eigenvectors (Eigen::MatrixXd).
+ *     - D_keep: Remaining diagonal entries after deflation (Eigen::VectorXd).
+ *     - v_keep: Remaining rank-one vector entries after deflation (Eigen::VectorXd).
+ *     - P_final: Combined permutation & rotation as an Eigen::SparseMatrix<double>.
+ */
+std::tuple<
+    Eigen::VectorXd,
+    Eigen::MatrixXd,
+    Eigen::VectorXd,
+    Eigen::VectorXd,
+    Eigen::SparseMatrix<double>>
+deflateEigenpairs(
+    const Eigen::VectorXd& D,
+    Eigen::VectorXd v,
+    double beta,
+    double tol_factor = 1e-12
+) {
+    int n = D.size();
+    // 1) Build full matrix M and compute norm for tolerance
+    Eigen::MatrixXd M = D.asDiagonal();          
+    M += beta * v * v.transpose();               
+    double norm_T = M.norm();                    
+    double tol = tol_factor * norm_T;            
+
+    // 2) Prepare containers for deflation
+    std::vector<int> keep_indices, deflated_indices;
+    Eigen::VectorXd deflated_eigvals = Eigen::VectorXd::Zero(n);
+    std::vector<Eigen::VectorXd> deflated_eigvecs_list;
+    int j = 0;  
+
+    // 3) Zero-component deflation
+    for (int i = 0; i < n; ++i) {
+        if (std::abs(v(i)) < tol) {
+            deflated_eigvals(j) = D(i);
+            Eigen::VectorXd e_vec = Eigen::VectorXd::Zero(n);
+            e_vec(i) = 1.0;
+            deflated_eigvecs_list.push_back(e_vec);
+            deflated_indices.push_back(i);
+            ++j;
+        } else {
+            keep_indices.push_back(i);
+        }
+    }
+
+    // 4) Build permutation P: [keep_indices, deflated_indices]
+    std::vector<int> new_order;
+    new_order.reserve(n);
+    new_order.insert(new_order.end(), keep_indices.begin(), keep_indices.end());
+    new_order.insert(new_order.end(), deflated_indices.begin(), deflated_indices.end());
+    Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> P(n);
+    P.indices() = Eigen::VectorXi::Map(new_order.data(), n);
+
+    // 5) Extract subproblem D_keep and v_keep
+    Eigen::VectorXd D_keep(static_cast<int>(keep_indices.size()));
+    Eigen::VectorXd v_keep(static_cast<int>(keep_indices.size()));
+    for (int idx = 0; idx < static_cast<int>(keep_indices.size()); ++idx) {
+        D_keep(idx) = D(keep_indices[idx]);
+        v_keep(idx) = v(keep_indices[idx]);
+    }
+
+    // 6) Givens rotations for near-duplicate entries
+    std::unordered_set<int> to_check;
+    to_check.reserve(keep_indices.size());
+    for (int i = 0; i < static_cast<int>(keep_indices.size()); ++i)
+        to_check.insert(i);
+    std::vector<std::tuple<int,int,double,double>> rotations;
+    std::vector<int> vec_idx_list;
+    std::vector<int> to_check_copy(to_check.begin(), to_check.end());
+
+    for (size_t idx_i = 0; idx_i + 1 < to_check_copy.size(); ++idx_i) {
+        int i = to_check_copy[idx_i];
+        if (to_check.find(i) == to_check.end()) continue;
+        for (int k = i + 1; k < static_cast<int>(D_keep.size()); ++k) {
+            if (std::abs(D_keep(k) - D_keep(i)) < tol) {
+                to_check.erase(k);
+                double r = std::hypot(v_keep(i), v_keep(k));
+                double c = v_keep(i) / r;
+                double s = -v_keep(k) / r;
+                v_keep(i) = r;
+                v_keep(k) = 0.0;
+                rotations.emplace_back(i, k, c, s);
+                deflated_eigvals(j) = D_keep(i);
+                ++j;
+                // local eigenvector in full basis
+                Eigen::VectorXd tmp = Eigen::VectorXd::Zero(n);
+                tmp(k) = c;
+                tmp(i) = s;
+                deflated_eigvecs_list.push_back(P.transpose() * tmp);
+                vec_idx_list.push_back(k);
+            }
+        }
+    }
+
+    // 7) Final ordering after rotations
+    std::vector<int> final_order(to_check.begin(), to_check.end());
+    final_order.insert(final_order.end(), vec_idx_list.begin(), vec_idx_list.end());
+
+    // 8) Resize deflated_eigvals to actual number found
+    deflated_eigvals.conservativeResize(j);
+
+    // 9) Build P2: accumulate sparse Givens rotations
+    Eigen::SparseMatrix<double> P2(n,n);
+    P2.setIdentity();
+    for (auto &rot: rotations) {
+        int i,k; double c,s; std::tie(i,k,c,s) = rot;
+        Eigen::SparseMatrix<double> G(n,n);
+        G.setIdentity();
+        G.coeffRef(i,i) = c;
+        G.coeffRef(k,k) = c;
+        G.coeffRef(i,k) = -s;
+        G.coeffRef(k,i) = s;
+        P2 = P2 * G;
+    }
+
+    // 10) Build permutation matrix P3 from final_order
+    Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> P3(n);
+    P3.indices() = Eigen::VectorXi::Map(final_order.data(), n);
+
+    // 11) Convert P and P3 to sparse<double> and combine all transforms
+    Eigen::SparseMatrix<double> P_sparse  = P.toDenseMatrix().cast<double>().sparseView();
+    Eigen::SparseMatrix<double> P3_sparse = P3.toDenseMatrix().cast<double>().sparseView();
+    Eigen::SparseMatrix<double> P_final   = P3_sparse * P2 * P_sparse;
+
+    // 12) Extract final D_keep and v_keep for reduced problem
+    std::vector<int> final_keep(to_check.begin(), to_check.end());
+    Eigen::VectorXd D_keep_final(static_cast<int>(final_keep.size()));
+    Eigen::VectorXd v_keep_final(static_cast<int>(final_keep.size()));
+    for (int idx = 0; idx < static_cast<int>(final_keep.size()); ++idx) {
+        D_keep_final(idx) = D_keep(final_keep[idx]);
+        v_keep_final(idx) = v_keep(final_keep[idx]);
+    }
+
+    // 13) Assemble deflated eigenvectors matrix
+    Eigen::MatrixXd deflated_eigvecs(n, static_cast<int>(deflated_eigvecs_list.size()));
+    for (int col = 0; col < static_cast<int>(deflated_eigvecs_list.size()); ++col) {
+        deflated_eigvecs.col(col) = deflated_eigvecs_list[col];
+    }
+
+    return std::make_tuple(
+        deflated_eigvals,
+        deflated_eigvecs.transpose(),
+        D_keep_final,
+        v_keep_final,
+        P_final
+    );
+}
+
+
 // PYTHON BINDINGS USING PYBIND11
 
 PYBIND11_MODULE(cxx_utils, m) {
@@ -593,4 +771,5 @@ PYBIND11_MODULE(cxx_utils, m) {
     m.def("QR_algorithm", &QR_algorithm, py::arg("diag"), py::arg("off_diag"), py::arg("tol")=1e-8, py::arg("max_iter")=5000);
     m.def("Eigen_value_calculator", &Eigen_value_calculator, py::arg("diag"), py::arg("off_diag"), py::arg("tol")=1e-8, py::arg("max_iter")=5000);
     m.def("secular_solver_cxx", &secular_solver, py::arg("rho"), py::arg("d"), py::arg("v"), py::arg("indices"));
+    m.def("deflate_eigenpairs_cxx", &deflateEigenpairs, py::arg("D"), py::arg("v"), py::arg("beta"), py::arg("tol_factor") = 1e-12);
 }

src/pyclassify/eigenvalues.py

@@ -158,7 +158,7 @@ def power_method_numba(A, max_iter=500, tol=1e-7, x=None):
     return x @ A @ x
 
 
-@jit(nopython=True)
+#@jit(nopython=True)
 def Lanczos_PRO(A, q, m=None, tol=np.sqrt(np.finfo(float).eps)):
     r"""
     Perform the Lanczos algorithm for symmetric matrices.

src/pyclassify/parallel_tridiag_eigen.py

@@ -1,12 +1,15 @@
 from mpi4py import MPI
 import numpy as np
 from time import time
-from .cxx_utils import QR_algorithm, secular_solver_cxx
+from pyclassify.cxx_utils import QR_algorithm, secular_solver_cxx, deflate_eigenpairs_cxx
+from zero_finder import secular_solver_python as secular_solver
 from line_profiler import profile, LineProfiler
 import scipy.sparse as sp
 from line_profiler import LineProfiler
 import scipy
 from pyclassify.utils import make_symmetric
+from pyclassify.eigenvalues import EigenSolver, Lanczos_PRO
+
 
 profile = LineProfiler()
 
@@ -255,18 +258,22 @@ def parallel_tridiag_eigen(
         n1 = len(eigvals_left)
         D = np.concatenate((eigvals_left, eigvals_right))
         v_vec = np.concatenate((eigvecs_left[-1, :], eigvecs_right[0, :]))
-        deflated_eigvals, deflated_eigvecs, D_keep, v_keep, P = deflate_eigenpairs(
-            D, v_vec, tol_factor
+        deflated_eigvals, deflated_eigvecs, D_keep, v_keep, P = deflate_eigenpairs_cxx(
+            D, v_vec, beta, tol_factor
         )
+        # deflated_eigvals, deflated_eigvecs, D_keep, v_keep, P = deflate_eigenpairs(
+        #     D, v_vec, beta, tol_factor
+        # )
         reduced_dim = len(D_keep)
         if D_keep.size > 0:
             # M = np.diag(D_keep) + beta * np.outer(v_keep, v_keep)
             # lam , _= np.linalg.eigh(M)
             idx = np.argsort(D_keep)
             idx_inv = np.arange(0, reduced_dim)
             idx_inv = idx_inv[idx]
+            #lam, changing_position, delta = secular_solver( beta, D_keep[idx], v_keep[idx] )
             lam, changing_position, delta = secular_solver_cxx(
-                beta, D_keep[idx], v_keep[idx]
+                beta, D_keep[idx], v_keep[idx], np.arange(reduced_dim)
             )
             lam = np.array(lam)
             # #diff=lam_s-lam
@@ -429,7 +436,7 @@ def parallel_tridiag_eigen(
     #     if colnorm > 1e-14:
     #         local_block[:, c] /= colnorm
 
-    # gathered_blocks = comm.gather(local_block, root=0)
+    # gathranered_blocks = comm.gather(local_block, root=0)
     # gathered_cols = comm.gather(my_cols, root=0)
 
     # if rank == 0:
@@ -481,24 +488,27 @@ def parallel_eigen(
 if __name__ == "__main__":
     comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
-    n = 2000
+    n = 1000
     if rank == 0:
+        
         # import debugpy
         # port = 5678 + rank  # 5678 for rank 0, 5679 for rank 1
         # debugpy.listen(("localhost", port))
         # print(f"Rank {rank} waiting for debugger attach on port {port}")
         # debugpy.wait_for_client()
         np.random.seed(42)
-        main_diag = np.ones(n, dtype=np.float64) * 2
-        off_diag = np.ones(n - 1, dtype=np.float64)
-        eig = np.arange(1, n + 1)
-        A = np.diag(eig)
-        U = scipy.stats.ortho_group.rvs(n)
-
-        A = U @ A @ U.T
-        A = make_symmetric(A)
-        Lanc = EigenSolver(A)
-        _, main_diag, off_diag = Lanc.Lanczos_PRO(q=np.ones_like(eig), tol=1e-12)
+        main_diag = np.ones(n, dtype=np.float64) * 2.0
+        off_diag = np.ones(n - 1, dtype=np.float64) *1.0
+        # eig = np.arange(1, n + 1)
+        # A = np.diag(eig)
+        # U = scipy.stats.ortho_group.rvs(n)
+
+        # A = U @ A @ U.T
+        # A = make_symmetric(A)
+
+        # #Lanc = EigenSolver(A)
+        # #_, main_diag, off_diag = Lanc.Lanczos_PRO(q=np.ones_like(eig), tol=1e-12)
+        # _, main_diag, off_diag = Lanczos_PRO(A, np.ones_like(eig)*1.0, tol=1e-12)
         T = np.diag(main_diag) + np.diag(off_diag, 1) + np.diag(off_diag, -1)
         eig_numpy, eig_vec_numpy = np.linalg.eigh(T)
         # print(eig_numpy)
@@ -535,13 +545,9 @@ def parallel_eigen(
 
         print("Norm difference eigenaval", np.linalg.norm(eig_numpy - eigvals, np.inf))
 
-        for count, i in enumerate(eigvecs[0, :]):
-            if i < 0:
-                eigvecs[:, count] = (-1) * eigvecs[:, count]
+  
+        check_column_directions(eigvecs, eig_vec_numpy)
 
-        for count, i in enumerate(eig_vec_numpy[0, :]):
-            if i < 0:
-                eig_vec_numpy[:, count] = (-1) * eig_vec_numpy[:, count]
 
         # import sys
         # np.set_printoptions(threshold=sys.maxsize)
@@ -552,13 +558,13 @@ def parallel_eigen(
 
         # # print("Eigenvector solver:\n", eigvecs)
         # # print("Eigenvector numpy:\n", eig_vec_numpy)
-        # print("\n\n\nDifference :\n", eig_vec_numpy - eigvecs)
-        # diff= eig_vec_numpy-eigvecs
-        # flat_idx = np.argmax(diff)          # → 5   (counting row-major: 0..8)
+        #print("\n\n\nDifference :\n", eig_vec_numpy - eigvecs)
+        diff= eig_vec_numpy-eigvecs
+        flat_idx = np.argmax(diff)          # → 5   (counting row-major: 0..8)
 
         # # If you want row/column coordinates instead of the flattened index:
-        # row, col = np.unravel_index(flat_idx, diff.shape)   # → (1, 2)
-        # print(np.max(np.abs(diff),LaEig_vec axis=0))
+        row, col = np.unravel_index(flat_idx, diff.shape)   # → (1, 2)
+        print(np.max(np.abs(diff)))
         # print("\n\n", eig_vec_numpy[:, col], eigvecs[:, col])
         # print("Norm difference eigenvec", np.linalg.norm(eig_vec_numpy-eigvecs, np.inf))