[Rcpp-devel] RcppEigen now on CRAN
Douglas Bates
bates at stat.wisc.edu
Tue Jun 28 03:16:49 CEST 2011
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
>>
>>
>
More information about the Rcpp-devel
mailing list