[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