Eigenvalues calculations in C-within-R codes

Question

I am writing R code which uses compiled C code.

From "Writing R Extensions" document, I learned there are many R executable/DLL that can be called from C code. The header file ‘Rmath.h’ lists many functions that are available and its source codes are listed on the this website: http://svn.r-project.org/R/trunk/src/nmath/

I need to calculate singular value decomposition of many matrices, however I do not find subroutines which does this on the above website. (So I am assuming that Rmath.h does not contain SVD subroutines) Is there simple way to do eigenvalue calculations in C-within-R code?

Thank you very much.

Answer 1

If you're open to using Rcpp and its related packages, the existing examples for the fastLm() show you how to do this with

• Eigen via RcppEigen
• Armadillo via RcppArmadillo
• the GSL via RcppGSL

where the latter two will use the same BLAS as R, and Eigen has something internal that can often be faster. All the packages implement lm() using the decompositions offered by the language (often using solve or related, but switching to SVD is straightforward once you have the toolchain working for you).

Edit: Here is an explicit example. Use the following C++ code:

#include <RcppArmadillo.h>   

// [[Rcpp::depends(RcppArmadillo)]] 

// [[Rcpp::export]]   
arma::vec getEigen(arma::mat M) { 
    return arma::eig_sym(M);      
}    

save in a file "eigenEx.cpp". Then all it takes is this R code:

library(Rcpp)               ## recent version for sourceCpp()
sourceCpp("eigenEx.cpp")    ## converts source file into getEigen() we can call

so that we can run this:

set.seed(42)        
X <- matrix(rnorm(3*3), 3, 3)
Z <- X %*% t(X)              

getEigen(Z)                 ## calls function created above

and I get the exact same eigen vector as from R. It really doesn't get much easier.

It also lets Armadillo chose what method to use for the Eigen decomposition, and as David hinted, this is something quicker than a full-blown SVD (see the Armadillo docs).

Answer 2

You can use one of many available Linear Algebra (lapack) libraries. Here is a link explaining how to get lapack libraries for windows. GOTOBLAS, and ACML are free. MKL is also free for non-commercial use. Once you have the libraries installed, the function you are looking for is sgesvd (for float) or dgesvd (for double).

Here are a couple of examples from Intel on how to use gesvd.

1. Row Major
2. Col Major

In case you are using Linux, Check out GNU SL and Eigen. These libraries can usually be installed from the package manager of the distribution.