RcppPseudoinv

C++ Function Reference

1 Signature

Eigen::MatrixXd BigDataStatMeth::RcppPseudoinv(Eigen::MatrixXd *A, Rcpp::Nullable< int > threads=R_NilValue)

2 Description

Compute pseudoinverse of in-memory matrix.

3 Parameters

  • A (Eigen::MatrixXd *): Input matrix to compute pseudoinverse of
  • threads (Rcpp::Nullable< int >): Number of threads for parallel processing (optional)

4 Returns

Pseudoinverse matrix

5 Details

Computes the Moore-Penrose pseudoinverse using SVD decomposition with parallel processing support.

6 Call Graph

Function dependencies

7 Source Code

File: inst/include/hdf5Algebra/matrixPseudoinverse.hpp • Lines 81-129

inline Eigen::MatrixXd RcppPseudoinv(Eigen::MatrixXd* A, 
                                            Rcpp::Nullable<int> threads = R_NilValue)
{
    
    char Schar='S';
    char Cchar='C';
    int ione = 1; 
    double done = 1.0;
    double dzero = 0.0;
    
    int m = A->rows();
    int n = A->cols();
    int lda = m;
    int ldu = std::max(1,m);
    int ldvt =  std::max(1, std::min(m, n));
    int k = std::min(m,n);
    int lwork = std::max( 1, 4*std::min(m,n)* std::min(m,n) + 7*std::min(m, n) );
    int info = 0;
    
    Eigen::VectorXd s = Eigen::VectorXd::Zero(k);
    Eigen::VectorXd work = Eigen::VectorXd::Zero(lwork);
    Eigen::MatrixXd u = Eigen::MatrixXd::Zero(ldu,k);
    Eigen::MatrixXd vt = Eigen::MatrixXd::Zero(ldvt,n);
    
    
    // dgesvd_( char JOBU, char JOBVT, int M, int N, double* A, int LDA, double* S, double* U, int LDU, double* VT, int LDVT, double WORK, int LWORK, int INFO  );
    dgesvd_( &Schar, &Schar, &m, &n, A->data(), &lda, s.data(), u.data(), &ldu, vt.data(), &ldvt, work.data(), &lwork, &info);
    Eigen::MatrixXd pinv = Eigen::MatrixXd::Zero(n,m);
    
    //.OpenMP.// omp_set_dynamic(1);
    
#pragma omp parallel for num_threads(get_number_threads(threads, R_NilValue))
    for (int i = 0; i < k; i++){
        double tempS;
        if(s[i] > 1.0e-9)
            tempS = 1.0/s[i];
        else
            tempS = s[i];
        
        // zscal_ (int* N, double* DA, double* DX, int* INCX )
        dscal_( &m, &tempS, &(u(i*ldu)), &ione );
    }
    // dgemm_( char TRANSA, char TRANSB, int M, int N, int K, double ALPHA, double* A, int LDA, double* B, int LDB, double BETA, double* C, int LDC )   
    dgemm_( &Cchar, &Cchar, &n, &m, &k, &done, vt.data(), &k, u.data(), &m, &dzero, pinv.data(), &n );
    
    
    return(pinv);
    
}

8 Usage Example

#include "BigDataStatMeth.hpp"

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