5

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.

FairyOnIce
  • 2,526
  • 7
  • 27
  • 48

2 Answers2

4

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).

Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725
1

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.

Pavan Yalamanchili
  • 12,021
  • 2
  • 35
  • 55
  • 2
    Pavan: All true, all irrelevant as ThePricess wants to work from R, not the command-line. – Dirk Eddelbuettel Jan 04 '13 at 21:39
  • @DirkEddelbuettel ThePirncess had this in the question "I learned there are many R executable/DLL that can be called *from* C code." So I assumed, she eventually wanted to have the decompositions in C rather than R. – Pavan Yalamanchili Jan 04 '13 at 21:43
  • 1
    In which case ThePrincess would be wrong/confused as that header is NOT for external use by C programs. Also not the "C-within-R" requirement though, so I think my read is closer than yours. – Dirk Eddelbuettel Jan 04 '13 at 21:44
  • I am not working on R-command-line. Sorry for the confusion. I want the decomposition result in C. I wanted to avoid using Rcpp because that would require me to change significant part of my C codes. – FairyOnIce Jan 04 '13 at 21:47
  • Nobody forces you to give up C, or C conventions. But you get data in and out with WAY less code using Rcpp, and your add-ons can still be C-style. Your choice. – Dirk Eddelbuettel Jan 04 '13 at 21:59
  • 2
    @ThePrincess just use Rcpp and RcppEigen or RcppGSL only for the part you need the svd, make it a separate function/header – pyCthon Jan 04 '13 at 22:02