[Rcpp-devel] Principal Component in Armadillo

[Dirk Eddelbuettel]

On 11 March 2013 at 22:42, João Daniel Nunes Duarte wrote:
| Hello,
| 
| I'm trying to calculate principal component using 'princomp' function from
| RcppArmadillo. Here's the cpp code:
| 
| #include <RcppArmadillo.h>
| 
| RcppExport SEXP pca(SEXP mats) { 
|   try {
|    
|     Rcpp::NumericMatrix matr(mats);
|     int n = matr.nrow(), k = matr.ncol();   
| 
|     arma::mat mat(matr.begin(), n, k, false);
|     arma::colvec pca;
|    
|     pca = arma::princomp(mat);
|  
|     return Rcpp::wrap(mat);
|    
|   } catch(...) {
|     ::Rf_error("c++ error");
|   }
| }
| 
| However, when I "R CMD check" the package, I get the following error:
| 
| ** testing if installed package can be loaded
| Error in dyn.load(file, DLLpath = DLLpath, ...) :
|   unable to load shared object '/home/tecto/cpp/Rcpp/amora.Rcheck/amora/libs/
| amora.so':
|   /home/tecto/cpp/Rcpp/amora.Rcheck/amora/libs/amora.so: undefined symbol:
| dgesvd_
| Error: loading failed
| Execution halted
| 
| I've read the Armadillo documentation for that function (http://
| arma.sourceforge.net/docs.html#princomp), however it was not clear to me how to
| use it correctly.
| 
| Does my code have any mistake?

Yes, several in fact.  princomp() returns a matrix, not a vector. You were
trying to return mat when you probably meant to return pca.

As for your linking error: should not happen.  On all systems, RcppArmadillo
should get these BLAS functions from R by linking to R.  Unless your R is
built in a weird way.  

Here is what I get in the slightly corrected (and rewritten to use
sourceCpp() instead) example below:


R> sourceCpp("/tmp/joao.cpp")

R> M <- matrix(1:9, 3, 3)

R> localpca(M)
   1.0000   4.0000   7.0000
   2.0000   5.0000   8.0000
   3.0000   6.0000   9.0000

   0.5774   0.8165        0
   0.5774  -0.4082  -0.7071
   0.5774  -0.4082   0.7071

        [,1]      [,2]      [,3]
[1,] 0.57735  0.816497  0.000000
[2,] 0.57735 -0.408248 -0.707107
[3,] 0.57735 -0.408248  0.707107
R> 


The program follows:

-----------------------------------------------------------------------------

#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
Rcpp::NumericMatrix localpca(Rcpp::NumericMatrix matr) {
  int n = matr.nrow(), k = matr.ncol();
  arma::mat mat(matr.begin(), n, k, false);
  Rcpp::Rcout << mat << std::endl;
  arma::mat pca = arma::princomp(mat);
  Rcpp::Rcout << pca << std::endl;
  return Rcpp::wrap(pca);
}

/*** R
M <- matrix(1:9, 3, 3)
localpca(M)
*/

-----------------------------------------------------------------------------

And for what it is worth, R returns the same (modulo ordering):

R> prcomp(M)
Standard deviations:
[1] 1.73205 0.00000 0.00000

Rotation:
         PC1       PC2       PC3
[1,] 0.57735  0.000000  0.816497
[2,] 0.57735 -0.707107 -0.408248
[3,] 0.57735  0.707107 -0.408248



Hth,  Dirk

| 
| Cheers,
| 
| João Daniel
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com