[Rcpp-devel] RcppEigen now on CRAN

[Douglas Bates]

On Mon, Jun 27, 2011 at 3:49 PM, baptiste auguie
<baptiste.auguie at googlemail.com> wrote:
> Hi,
>
> Another august piece work, thanks.
> As a happy user of Armadillo (guess what, this week my packages even
> passed win-builder! Yeah!) I have a general question: between (Rcpp),
> RcppGSL, RcppArmadillo, RcppEigen, we now have so much choice that it
> may be difficult to decide which framework to use for a particular
> project. Are there comparative charts somewhere, in terms of function
> coverage, speed? (those I found did not cover these three libraries).

There is a benchmark comparing different methods of obtaining the
coefficients and the standard errors for linear regression models.
You can run the benchmark as

library(RcppEigen)
source(file.path(system.file("examples", package="RcppEigen"),
"lmBenchmark.R"), echo=TRUE)

Be aware that it takes a while to run.  On my laptop the results,
without RcppGSL which performs terribly on this benchmark, are

lm benchmark for n = 100000 and p = 40: nrep = 20
     test  relative elapsed user.self sys.self
3    LDLt  1.000000   3.381      3.31     0.00
7     LLt  1.007394   3.406      3.31     0.00
5 SymmEig  2.159716   7.302      7.08     0.08
8    arma  5.733215  19.384     19.13     0.01
6      QR  6.087548  20.582     20.24     0.09
2   PivQR  6.506655  21.999     21.65     0.07
1  lm.fit 10.447797  35.324     32.44     2.68
4     SVD 43.231588 146.166    142.41     1.77

I probably should skip the SVD method on large examples.  Also, I
forgot a print() in the source code, which is why you need to use
echo=TRUE.

Sourcing that file defines a function called do_bench with three
arguments: n, p and nrep.  The script requires the package rbenchmark
and will use RcppArmadillo and/or RcppGSL if they are available.

On this example Armadillo is slightly faster than the QR and
Column-pivoted QR decomposition as implemented in RcppEigen.  However,
the RcppArmadillo and the RcppGSL fastLm functions do not use a
rank-revealing decomposition whereas R's lm.fit and the LDLt, SymmEig,
PivQR and SVD methods from RcppEigen do.  The problem with using a
non-rank-revealing decomposition is that you can produce nonsense
answers and not know it.

> It could also be very useful to have a kind of Rosetta stone, akin to
> http://arma.sourceforge.net/docs.html#syntax but comparing (R) | Rcpp
> |  RcppGSL | RcppArmadillo | RcppEigen, ideally with wise comments
> from advanced users. Where could we build such meta-documentation?

Sounds like a good idea.  I haven't tried any of the Wikis associated
with R but would be willing to contribute if someone had a
recommendation on which one to use.
> Just a thought,
>
> Best regards,
>
> baptiste
>
>
>
>
> On 28 June 2011 07:55, Douglas Bates <bates at stat.wisc.edu> wrote:
>> A new package, RcppEigen, which provides linkage from R to the Eigen
>> C++ template library for linear algebra (http://eigen.tuxfamily.org)
>> via Rcpp, is now available on CRAN.    The Eigen library is a pure
>> template library that does not depend on the BLAS or Lapack, thus
>> avoiding some difficulties such as Baptiste Auguste encountered
>> regarding inconsistent versions of the Lapack routines being available
>> on different platforms.
>>
>> A great deal of care has gone into the development of the Eigen
>> library to ensure good performance through the use of vectorization
>> instructions such as SSE2, SSE3 and ARM_NEON, when available, and by
>> using block structured algorithms.  The files RcppEigen/src/fastLm.h
>> and RcppEigen/src/fastLm.cpp in the source package provide examples of
>> the use of some of the decompositions available in Eigen.  Eigen
>> supports both dense and sparse matrices, with various types of
>> triangular, symmetric (or selfAdjoint) and diagonal adaptors.  To
>> learn about the facilities available in Eigen I would recommend
>> starting with the tutorial at
>> http://eigen.tuxfamily.org/dox/index.html
>>
>> Like the Rcpp and RcppArmadillo packages, RcppEigen has a plugin for
>> cxxfunction from the inline package, to allow for rapid prototyping.
>> There is also a function RcppEigen.package.skeleton to aid in creating
>> a package that uses the RcppEigen headers.
>>
>> Because Eigen essentially creates its own BLAS functions, the size of
>> the libs directory for packages using RcppEigen can be very large and
>> compiling source files that include the RcppEigen.h file can be
>> relatively slow.  As with most template libraries, the idea is to
>> off-load the complexity of the code onto the compiler to sort out,
>> with the penalty that instantiating templates may slow down the
>> compiler.
>>
>> I enclose an example of a short C++ function using RcppEigen and a
>> slightly longer and more complex C++ function, compiled with inline.
>> Some things to notice:
>>
>> - in keeping with the functional programming semantics of R, it is a
>> good idea to declare Rcpp objects created from the input SEXP
>> arguments as const.
>>
>> - a Eigen matrix of double precision values has type Eigen::MatrixXd,
>> similarly Eigen::MatrixXi for integer matrices and Eigen::MatrixXcd
>> for matrices of std::complex<double> (see the tutorial for the
>> explanation of the X in MatrixXd)
>>
>> - the Eigen::VectorXd class is a column vector of doubles, similarly
>> Eigen::VectorXi and Eigen::VectorXcd.  A row vector is
>> Eigen::RowVectorXd, etc.  These are all just specializations of the
>> Matrix template (Vectors have 1 column, RowVectors have 1 row, fixed
>> at compile time).
>>
>> - to obtain an Eigen matrix or vector from an Rcpp object, without
>> copying the contents, use the Eigen::Map adaptor class.
>>
>> - the Eigen classes ArrayXd and ArrayXXd are similar representations
>> to VectorXd and MatrixXd but allow for element-wise operations.  The
>> .array() method creates a view of a Vector or Matrix object as an
>> Array.  The .matrix() method maps the other way.
>>
>> - the Eigen::Matrix and Eigen::Array classes have methods for
>> Rcpp::wrap defined in the RcppEigen headers.
>>
>> - Adaptor methods like .triangularView<type>(),
>> .selfadjointView<type>() and .diagonalView() (applied to a Vector) are
>> the common ways of restricting to specialized forms of matrices.
>>
>> - most methods for Eigen classes return *this as a reference so that
>> methods can be chained.
>>
>> - like Armadillo, Eigen performs lazy evaluation of expressions to
>> reduce the number of intermediate results that must be stored.
>>
>> As with any new libraries, Eigen takes some getting used to, but
>> overall I have found it very good to work with and recommend its use.
>>
>> _______________________________________________
>> 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
>>
>>
>