Skip to content

[Rcpp-devel] Principal Component in Armadillo

4 messages · João Daniel Nunes Duarte, Dirk Eddelbuettel

Original
João Daniel Nunes Duarte
# 
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?

Cheers,

Jo?o Daniel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130311/3d351363/attachment.html>
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
João Daniel Nunes Duarte
# 
Thanks for your explanation, Dirk! Indeed, I looked at Armadillo
documentation, but it was not clear what was the kind of return.

Regarding the linking error, I had created the package using
Rcpp.package.skeleton (because I didn't know I was going to use
RcppArmadillo), and then I made the necessary changes according to
documentation (
http://romainfrancois.blog.free.fr/index.php?post/2010/05/19/RcppArmadillo-0.2.1).
So I created the package again using "RcppArmadillo.package.skeleton" and
everything worked fine!

Actually, what "RcppArmadillo.package.skeleton" creates is slightly
different from what is described on documentation (maybe it is system
dependent):

- Documentation

# Makevars
PKG_LIBS = $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()" )
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

# Makevars.win
PKG_LIBS = $(shell Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS)
$(BLAS_LIBS) $(FLIBS)

- Documentation

# Makevars
PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS)
$(BLAS_LIBS) $(FLIBS)

# Makevars.win
PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()")
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)


So, it seems that the R libraries was not being imported during the package
compilation; problem solved!

Regards,

Jo?o Daniel


2013/3/11 Dirk Eddelbuettel <edd at debian.org>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130312/c52fc95d/attachment.html>
Dirk Eddelbuettel
#
On 12 March 2013 at 14:10, Jo?o Daniel Nunes Duarte wrote:
| Thanks for your explanation, Dirk! Indeed, I looked at Armadillo documentation,
| but it was not clear what was the kind of return.
| 
| Regarding the linking error, I had created the package using
| Rcpp.package.skeleton (because I didn't know I was going to use RcppArmadillo),
| and then I made the necessary changes according to documentation (http://
| romainfrancois.blog.free.fr/index.php?post/2010/05/19/RcppArmadillo-0.2.1). So

Eek. That is a post that is almost three years old.  Why don't you look at
the documentation in the package you use?

| I created the package again using "RcppArmadillo.package.skeleton" and
| everything worked fine!
| 
| Actually, what "RcppArmadillo.package.skeleton" creates is slightly different
| from what is described on documentation (maybe it is system dependent):
| 
| - Documentation
| 
| # Makevars
| PKG_LIBS = $(shell $(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()" ) $(LAPACK_LIBS)
| $(BLAS_LIBS) $(FLIBS)
| 
| # Makevars.win
| PKG_LIBS = $(shell Rscript.exe -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $
| (BLAS_LIBS) $(FLIBS)
| 
| - Documentation
| 
| # Makevars
| PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS) $
| (BLAS_LIBS) $(FLIBS)
| 
| # Makevars.win
| PKG_LIBS = $(shell $(R_HOME)/bin/Rscript.exe -e "Rcpp:::LdFlags()") $
| (LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

You labeled both sides 'Documentation' ...   

Anyway, these are all identical apart from how Rscript is invoked.  They key
is the return of LdFlags combined with three values given.

Dirk