RcppbdEigen_spectra

C++ Function Reference

1 Signature

eigdecomp BigDataStatMeth::RcppbdEigen_spectra(const Eigen::MatrixXd &X, int k, const std::string &which="LM", int ncv=0, bool bcenter=false, bool bscale=false, double tol=1e-10, int max_iter=1000)

2 Description

Eigenvalue decomposition using Spectra (compatible with BigDataStatMeth SVD version)

3 Parameters

• X (const Eigen::MatrixXd &): Input matrix
• k (int): Number of eigenvalues to compute
• which (const std::string &): Which eigenvalues to compute (LM, SM, LR, SR, LI, SI, LA, SA)
• ncv (int): Number of Arnoldi vectors (if 0, uses auto-selection)
• bcenter (bool): Whether to center the data
• bscale (bool): Whether to scale the data
• tol (double): Convergence tolerance for Spectra
• max_iter (int): Maximum iterations for Spectra

4 Returns

eigdecomp structure containing results

5 Details

Uses the same Spectra version and patterns as the existing SVD implementation in BigDataStatMeth. Based on the RcppbdSVD function in matrixSvd.hpp.

6 Call Graph

7 Source Code

NoteImplementation

File: inst/include/hdf5Algebra/matrixEigenDecomposition.hpp • Lines 179-282

inline eigdecomp RcppbdEigen_spectra(const Eigen::MatrixXd& X, int k, const std::string& which = "LM", 
                                         int ncv = 0, bool bcenter = false, bool bscale = false,
                                         double tol = 1e-10, int max_iter = 1000) {
        
        eigdecomp reteig;
        Eigen::MatrixXd nX;
        // int nconv = 0;
        
        try {
            
            int n = X.rows();
            if (X.rows() != X.cols()) {
                Rf_error("Matrix must be square for eigendecomposition");
                return reteig;
            }
            
            // Better parameter selection following RSpectra defaults
            std::tie(k, ncv) = validateSpectraParams(n, k, ncv);
            
            // Improved symmetry detection
            reteig.is_symmetric = isMatrixSymmetric(X);
            
            if (reteig.is_symmetric) {
                // Use symmetric solver - following SVD pattern
                Eigen::MatrixXd Xcp;
                if (bcenter == true || bscale == true) {
                    nX = RcppNormalize_Data(X, bcenter, bscale, false);
                    Xcp = nX;  // For symmetric case, use matrix directly
                } else {
                    Xcp = X;
                }
                
                Spectra::DenseSymMatProd<double> op(Xcp);
                Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> eigs(op, k, ncv);
                
                // // Set tolerance and max iterations
                // eigs.set_max_iter(max_iter);
                
                // Initialize and compute with appropriate sort rule
                eigs.init();
                Spectra::SortRule sort_rule = getSymmetricSortRule(which);
                //..// nconv = eigs.compute(sort_rule, max_iter, tol);
                (void)eigs.compute(sort_rule, max_iter, tol);
                
                // Retrieve results
                if (eigs.info() == Spectra::CompInfo::Successful) {
                    reteig.eigenvalues_real = eigs.eigenvalues();
                    reteig.eigenvalues_imag = Eigen::VectorXd::Zero(k);
                    reteig.eigenvectors_real = eigs.eigenvectors();
                    reteig.eigenvectors_imag = Eigen::MatrixXd::Zero(n, k);
                    reteig.bconv = true;
                } else {
                    reteig.bconv = false;
                }
                
            } else {
                // For non-symmetric matrices, use general solver
                Eigen::MatrixXd Xcp;
                if (bcenter == true || bscale == true) {
                    nX = RcppNormalize_Data(X, bcenter, bscale, false);
                    Xcp = nX;
                } else {
                    Xcp = X;
                }
                
                Spectra::DenseGenMatProd<double> op(Xcp);
                Spectra::GenEigsSolver<Spectra::DenseGenMatProd<double>> eigs(op, k, ncv);
                
                // // Set tolerance and max iterations
                // eigs.set_max_iter(max_iter);
                
                // Initialize and compute with appropriate sort rule
                eigs.init();
                Spectra::SortRule sort_rule = getGeneralSortRule(which);
                //..// nconv = eigs.compute(sort_rule, max_iter, tol);
                (void)eigs.compute(sort_rule, max_iter, tol);
                
                // Retrieve results
                if (eigs.info() == Spectra::CompInfo::Successful) {
                    Eigen::VectorXcd eigenvals = eigs.eigenvalues();
                    Eigen::MatrixXcd eigenvecs = eigs.eigenvectors();
                    
                    reteig.eigenvalues_real = eigenvals.real();
                    reteig.eigenvalues_imag = eigenvals.imag();
                    reteig.eigenvectors_real = eigenvecs.real();
                    reteig.eigenvectors_imag = eigenvecs.imag();
                    reteig.bconv = true;
                } else {
                    reteig.bconv = false;
                }
            }
            
            reteig.bcomputevectors = true;
            
        } catch(std::exception &ex) {
            Rcpp::Rcout << "C++ exception RcppbdEigen_spectra: " << ex.what();
            reteig.bconv = false;
        } catch (...) {
            Rf_error("C++ exception RcppbdEigen_spectra (unknown reason)");
            reteig.bconv = false;
        }
        
        return reteig;
    }

8 Usage Example

#include "BigDataStatMeth.hpp"

// Example usage
auto result = RcppbdEigen_spectra(...);