[Rcpp-devel] RcppEigen now on CRAN
Douglas Bates
bates at stat.wisc.edu
Mon Jun 27 21:55:03 CEST 2011
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.
-------------- next part --------------
R version 2.13.0 (2011-04-13)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(inline)
> library(RcppEigen)
Loading required package: Rcpp
> src <- '
+ const NumericMatrix foo(foo_);
+ if (foo.rows() != foo.cols())
+ throw std::invalid_argument("foo_ must be a square, numeric matrix");
+
+ Eigen::MatrixXd fooinv = Eigen::Map<Eigen::MatrixXd>(foo.begin(), foo.rows(), foo.cols()).inverse();
+ return wrap(fooinv);
+ '
> ff <- cxxfunction(signature(foo_ = "matrix"), src, "RcppEigen")
>
> set.seed(1)
> (mm <- matrix(rnorm(16), nc = 4))
[,1] [,2] [,3] [,4]
[1,] -0.6264538 0.3295078 0.5757814 -0.62124058
[2,] 0.1836433 -0.8204684 -0.3053884 -2.21469989
[3,] -0.8356286 0.4874291 1.5117812 1.12493092
[4,] 1.5952808 0.7383247 0.3898432 -0.04493361
> (mmi <- solve(mm))
[,1] [,2] [,3] [,4]
[1,] -0.6213910 0.2532768 0.1732266 0.44441565
[2,] 1.5301385 -0.8532288 -0.8240908 0.26742439
[3,] -0.4197372 0.5574457 0.8748252 0.22924604
[4,] -0.5605101 -0.1913025 0.1990296 -0.09383139
> (ffi <- ff(mm))
[,1] [,2] [,3] [,4]
[1,] -0.6213910 0.2532768 0.1732266 0.44441565
[2,] 1.5301385 -0.8532288 -0.8240908 0.26742439
[3,] -0.4197372 0.5574457 0.8748252 0.22924604
[4,] -0.5605101 -0.1913025 0.1990296 -0.09383139
>
> src2 <- '
+ const NumericMatrix X(X_);
+ int n = X.rows(), p = X.cols();
+ const NumericVector y(y_);
+ if (y.size() != n)
+ throw std::invalid_argument("length(y_) != nrow(X_)");
+ if (n <= p)
+ throw std::invalid_argument("nrow(X_) <= ncol(X_)");
+
+ typedef Eigen::MatrixXd MatrixXd;
+ typedef Eigen::Map<MatrixXd> MMatrixXd; // shares storage
+ typedef Eigen::ColPivHouseholderQR<MatrixXd> ColPivQR;
+ typedef ColPivQR::PermutationType Permutation;
+ typedef Eigen::VectorXd VectorXd;
+ typedef Eigen::Map<VectorXd> MVectorXd;
+
+ const MMatrixXd A(X.begin(), X.rows(), X.cols());
+ const MVectorXd b(y.begin(), y.size());
+ const ColPivQR QR(A);
+ int rank = QR.rank();
+ VectorXd coef = QR.solve(b);
+ MatrixXd R = QR.matrixQR().topLeftCorner(p, p);
+ // strangely, there is no inverse method for triangular matrices
+ MatrixXd Rinv = R.triangularView<Eigen::Upper>().solve(MatrixXd::Identity(p, p));
+ VectorXd fitted = A * coef;
+ VectorXd resid = b - fitted;
+ double s2 = resid.squaredNorm()/(n - rank);
+ VectorXd se = QR.colsPermutation() *
+ (Rinv.rowwise().squaredNorm().array() * s2).sqrt().matrix();
+
+ return List::create(_["coefficients"] = coef,
+ _["df.residual"] = n - rank,
+ _["fitted.values"] = fitted,
+ _["rank"] = rank,
+ _["residuals"] = resid,
+ _["se"] = se,
+ _["s2"] = s2);
+ '
> fastLm <- cxxfunction(signature(X_ = "matrix", y_ = "numeric"), src2, "RcppEigen")
>
> lmfit <- with(datasets::trees, lm.fit(cbind(1, log(Girth)), log(Volume)))
> str(lmfit)
List of 8
$ coefficients : Named num [1:2] -2.35 2.2
..- attr(*, "names")= chr [1:2] "x1" "x2"
$ residuals : num [1:31] 0.0298 -0.0483 -0.1087 -0.0223 0.0727 ...
$ effects : Named num [1:31] -18.2218 2.8152 -0.1029 -0.0223 0.0721 ...
..- attr(*, "names")= chr [1:31] "x1" "x2" "" "" ...
$ rank : int 2
$ fitted.values: num [1:31] 2.3 2.38 2.43 2.82 2.86 ...
$ assign : NULL
$ qr :List of 5
..$ qr : num [1:31, 1:2] -5.57 0.18 0.18 0.18 0.18 ...
..$ qraux: num [1:2] 1.18 1.26
..$ pivot: int [1:2] 1 2
..$ tol : num 1e-07
..$ rank : int 2
..- attr(*, "class")= chr "qr"
$ df.residual : int 29
> flm <- with(datasets::trees, fastLm(cbind(1, log(Girth)), log(Volume)))
> str(flm)
List of 7
$ coefficients : num [1:2] -2.35 2.2
$ df.residual : int 29
$ fitted.values: num [1:31] 2.3 2.38 2.43 2.82 2.86 ...
$ rank : int 2
$ residuals : num [1:31] 0.0298 -0.0483 -0.1087 -0.0223 0.0727 ...
$ se : num [1:2] 0.2307 0.0898
$ s2 : num 0.0132
> all.equal(unname(lmfit$coefficients), flm$coefficients)
[1] TRUE
> all.equal(lmfit$residuals, flm$residuals)
[1] TRUE
> (csum <-coef(summary(lm(log(Volume) ~ log(Girth), datasets::trees))))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.353325 0.23066284 -10.20244 4.180317e-11
log(Girth) 2.199970 0.08983455 24.48913 6.364197e-21
> all.equal(flm$se, unname(csum[, "Std. Error"]))
[1] TRUE
>
> proc.time()
user system elapsed
21.11 1.25 22.02
More information about the Rcpp-devel
mailing list