[Rcpp-commits] r3063 - in pkg/RcppEigen: R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 14 22:08:03 CEST 2011
Author: dmbates
Date: 2011-06-14 22:08:03 +0200 (Tue, 14 Jun 2011)
New Revision: 3063
Modified:
pkg/RcppEigen/R/fastLm.R
pkg/RcppEigen/man/fastLm.Rd
pkg/RcppEigen/src/fastLm.cpp
Log:
Change order of fastLmPure arguments and the various fastLm* C++ functions. Add fastLmChol1 and fastLmChol2.
Modified: pkg/RcppEigen/R/fastLm.R
===================================================================
--- pkg/RcppEigen/R/fastLm.R 2011-06-14 00:17:36 UTC (rev 3062)
+++ pkg/RcppEigen/R/fastLm.R 2011-06-14 20:08:03 UTC (rev 3063)
@@ -17,11 +17,11 @@
## You should have received a copy of the GNU General Public License
## along with RcppEigen. If not, see <http://www.gnu.org/licenses/>.
-fastLmPure <- function(y, X) {
+fastLmPure <- function(X, y) {
stopifnot(is.matrix(X), is.numeric(y), NROW(y)==nrow(X))
- .Call("fastLm", y, X, package="RcppEigen")
+ .Call("fastLm", X, y, package="RcppEigen")
}
fastLm <- function(X, ...) UseMethod("fastLm")
@@ -31,7 +31,7 @@
X <- as.matrix(X)
y <- as.numeric(y)
- res <- fastLmPure(y, X)
+ res <- fastLmPure(X, y)
names(res$coefficients) <- colnames(X)
res$call <- match.call()
@@ -87,12 +87,12 @@
fastLm.formula <- function(formula, data=list(), ...) {
mf <- model.frame(formula=formula, data=data)
- x <- model.matrix(attr(mf, "terms"), data=mf)
+ X <- model.matrix(attr(mf, "terms"), data=mf)
y <- model.response(mf)
- res <- fastLm.default(x, y, ...)
- names(res$coefficients) <- colnames(x)
+ res <- fastLm.default(X, y, ...)
res$call <- match.call()
+ # I think this is redundant. The formula is available as res$call$formula
res$formula <- formula
res
}
Modified: pkg/RcppEigen/man/fastLm.Rd
===================================================================
--- pkg/RcppEigen/man/fastLm.Rd 2011-06-14 00:17:36 UTC (rev 3062)
+++ pkg/RcppEigen/man/fastLm.Rd 2011-06-14 20:08:03 UTC (rev 3063)
@@ -10,7 +10,7 @@
decomposition from the \code{Eigen} linear algebra library.
}
\usage{
-fastLmPure(y, X)
+fastLmPure(X, y)
fastLm(X, \dots)
\method{fastLm}{default}(X, y, \dots)
\method{fastLm}{formula}(formula, data = list(), \dots)
@@ -59,7 +59,7 @@
\code{fastLmPure} returns a list with three components:
\item{coefficients}{a vector of coefficients}
\item{rank}{a scalar denoting the computed rank of the model matrix}
- \item{df}{a scalar denoting the degrees of freedom in the model}
+ \item{df.residual}{a scalar denoting the degrees of freedom in the model}
\item{stderr}{a vector of the (estimated) standard errors of the coefficient estimates}
\item{s}{a scalar - the root mean square for residuals}
\item{residuals}{the vector of residuals}
Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp 2011-06-14 00:17:36 UTC (rev 3062)
+++ pkg/RcppEigen/src/fastLm.cpp 2011-06-14 20:08:03 UTC (rev 3063)
@@ -21,7 +21,7 @@
#include <RcppEigen.h>
-extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
+extern "C" SEXP fastLm(SEXP Xs, SEXP ys) {
using namespace Eigen;
using namespace Rcpp;
try {
@@ -67,7 +67,7 @@
return List::create(_["coefficients"] = coef,
_["rank"] = r,
- _["df"] = df,
+ _["df.residual"] = df,
_["perm"] = mqr.colsPermutation().indices(),
_["stderr"] = se,
_["s"] = s,
@@ -83,7 +83,101 @@
return R_NilValue; // -Wall
}
+extern "C" SEXP fastLmBench(SEXP Xs, SEXP ys) {
+ using namespace Eigen;
+ using namespace Rcpp;
+ try {
+ NumericMatrix X(Xs);
+ NumericVector y(ys);
+ size_t n = X.nrow(), p = X.ncol();
+ ptrdiff_t df = n - p;
+ if ((size_t)y.size() != n)
+ throw std::invalid_argument("size mismatch");
+
+ MatrixXd A = Map<MatrixXd>(X.begin(), n, p); // shares storage
+ VectorXd b = Map<VectorXd>(y.begin(), n);
+ HouseholderQR<MatrixXd> Aqr(A);
+ VectorXd coef = Aqr.solve(b);
+ double s = std::sqrt((b - A*coef).squaredNorm()/df);
+ MatrixXd Rinv((TriangularView<MatrixXd, Upper>(Aqr.matrixQR().topRows(p)))
+ .solve(MatrixXd::Identity(p, p)));
+ VectorXd se = Rinv.rowwise().norm() * s;
+ return List::create(_["coefficients"] = coef,
+ _["se"] = se,
+ _["df.residual"] = df);
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
+ return R_NilValue; // -Wall
+}
+extern "C" SEXP fastLmChol1(SEXP Xs, SEXP ys) {
+ using namespace Eigen;
+ using namespace Rcpp;
+ try {
+ NumericMatrix X(Xs);
+ NumericVector y(ys);
+ size_t n = X.nrow(), p = X.ncol();
+ ptrdiff_t df = n - p;
+ if (df <= 0l)
+ throw std::invalid_argument("nrow(X) > ncol(X) not satisfied");
+ if ((size_t)y.size() != n)
+ throw std::invalid_argument("size mismatch");
+
+ MatrixXd A = Map<MatrixXd>(X.begin(), n, p); // shares storage
+ VectorXd b = Map<VectorXd>(y.begin(), n);
+
+ LLT<MatrixXd> Ch(SelfAdjointView<MatrixXd, Lower>(MatrixXd::Zero(p, p)).rankUpdate(A.adjoint()));
+ VectorXd coef = Ch.solve(A.adjoint() * b);
+ double s = std::sqrt((b - A*coef).squaredNorm()/df);
+
+ VectorXd se = (Ch.matrixL().solve(MatrixXd::Identity(p, p))).colwise().norm() * s;;
+
+ return List::create(_["coefficients"] = coef,
+ _["se"] = se,
+ _["df.residual"] = df);
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
+ return R_NilValue; // -Wall
+}
+
+extern "C" SEXP fastLmChol2(SEXP Xs, SEXP ys) {
+ using namespace Eigen;
+ using namespace Rcpp;
+ try {
+ NumericMatrix X(Xs);
+ NumericVector y(ys);
+ size_t n = X.nrow(), p = X.ncol();
+ ptrdiff_t df = n - p;
+ if ((size_t)y.size() != n)
+ throw std::invalid_argument("size mismatch");
+
+ MatrixXd A = Map<MatrixXd>(X.begin(), n, p); // shares storage
+ VectorXd b = Map<VectorXd>(y.begin(), n);
+
+ LDLT<MatrixXd> Ch(SelfAdjointView<MatrixXd, Lower>(MatrixXd::Zero(p, p)).rankUpdate(A.adjoint()));
+ VectorXd coef = Ch.solve(A.adjoint() * b);
+ double s2 = (b - A*coef).squaredNorm()/df;
+
+ ArrayXd se = (Ch.solve(MatrixXd::Identity(p, p)).diagonal().array() * s2).sqrt();
+ NumericVector Rse(p);
+ std::copy(se.data(), se.data() + p, Rse.begin());
+
+ return List::create(_["coefficients"] = coef,
+ _["se"] = Rse,
+ _["df.residual"] = df);
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
+ return R_NilValue; // -Wall
+}
More information about the Rcpp-commits
mailing list