RcppbdSVD

C++ Function Reference

1 Signature

svdeig BigDataStatMeth::RcppbdSVD(Eigen::MatrixXd &X, int k, int ncv, bool bcenter, bool bscale)

2 Description

Compute SVD decomposition using Spectra eigenvalue solver.

3 Parameters

• X (Eigen::MatrixXd &): Input matrix for SVD computation
• k (int): Number of singular values/vectors to compute (0 = auto-select)
• ncv (int): Number of Arnoldi vectors to use (0 = auto-select)
• bcenter (bool): Whether to center the data before computation
• bscale (bool): Whether to scale the data before computation

4 Returns

svdeig Structure containing U, S, V matrices and computation status

5 Details

Performs Singular Value Decomposition of a matrix using the Spectra library for eigenvalue computation. This function computes SVD by solving the eigenvalue problems X^{TX and XX}T.

6 Call Graph

7 Source Code

NoteImplementation

File: inst/include/hdf5Algebra/matrixSvd.hpp • Lines 141-204

inline svdeig RcppbdSVD( Eigen::MatrixXd& X, int k, int ncv, bool bcenter, bool bscale )
    {
        svdeig retsvd;
        
        Eigen::MatrixXd nX;
        const Eigen::MatrixXd& Xref = (bcenter || bscale)
            ? (nX = RcppNormalize_Data(X, bcenter, bscale, false), nX)
            : X;
        
        const int m   = static_cast<int>(Xref.rows());
        const int n   = static_cast<int>(Xref.cols());
        const int kmax = std::min(m, n);
        
        // ── Full SVD: delegate to LAPACK dgesdd (fastest for full decomposition) ──
        if (k <= 0 || k >= kmax) {
            // RcppbdSVD_lapack accepts a const ref-compatible T; pass a copy if needed
            Eigen::MatrixXd Xcopy = Xref;   // dgesdd_ overwrites the input
            return RcppbdSVD_lapack(Xcopy, false, false, false);
        }
        
        // ── Truncated SVD: one Spectra solve on the smaller Gram matrix ──
        // Choose the side that produces the smaller Gram matrix.
        // Then recover the other factor analytically: V = Xref^T * U * D^{-1}
        // (or U = Xref * V * D^{-1} for the m<n case).
        if (ncv <= 0)  ncv = std::min(k + std::max(k, 10), kmax);
        if (ncv <= k)  ncv = k + 1;
        ncv = std::min(ncv, kmax);
        
        if (m >= n) {
            // Smaller Gram: X^T * X  (n × n). Solve for V, recover U.
            Eigen::MatrixXd XTX = Xref.transpose() * Xref;
            Spectra::DenseSymMatProd<double> op(XTX);
            Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> eigs(op, k, ncv);
            eigs.init();
            eigs.compute(Spectra::SortRule::LargestAlge);
            if (eigs.info() != Spectra::CompInfo::Successful) {
                retsvd.bokuv = false; retsvd.bokd = false; return retsvd;
            }
            retsvd.d = eigs.eigenvalues().cwiseSqrt();   // singular values
            retsvd.v = eigs.eigenvectors();               // n × k
            // U = X * V * D^{-1}
            retsvd.u = (Xref * retsvd.v).array().rowwise()
                / retsvd.d.transpose().array();    // m × k
        } else {
            // Smaller Gram: X * X^T  (m × m). Solve for U, recover V.
            Eigen::MatrixXd XXT = Xref * Xref.transpose();
            Spectra::DenseSymMatProd<double> op(XXT);
            Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> eigs(op, k, ncv);
            eigs.init();
            eigs.compute(Spectra::SortRule::LargestAlge);
            if (eigs.info() != Spectra::CompInfo::Successful) {
                retsvd.bokuv = false; retsvd.bokd = false; return retsvd;
            }
            retsvd.d = eigs.eigenvalues().cwiseSqrt();
            retsvd.u = eigs.eigenvectors();               // m × k
            // V = X^T * U * D^{-1}
            retsvd.v = (Xref.transpose() * retsvd.u).array().rowwise()
                / retsvd.d.transpose().array();    // n × k
        }
        
        retsvd.bokuv = true;
        retsvd.bokd  = true;
        return retsvd;
    }

8 Usage Example

#include "BigDataStatMeth.hpp"

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