[Rcpp-commits] r3097 - pkg/RcppEigen/src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Jun 22 00:17:36 CEST 2011


Author: dmbates
Date: 2011-06-22 00:17:36 +0200 (Wed, 22 Jun 2011)
New Revision: 3097

Modified:
   pkg/RcppEigen/src/fastLm.cpp
   pkg/RcppEigen/src/fastLm.h
Log:
Add an SVD method for fastLm

Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp	2011-06-21 02:33:41 UTC (rev 3096)
+++ pkg/RcppEigen/src/fastLm.cpp	2011-06-21 22:17:36 UTC (rev 3097)
@@ -85,6 +85,19 @@
     m_unsc      = Ch.solve(MatrixXd::Identity(m_p, m_p));
 }
 
+SVD::SVD(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
+    double  prec = std::sqrt(std::numeric_limits<double>::epsilon());
+    SVDType  UDV = X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
+    VectorXd   D = UDV.singularValues();
+    m_r          = std::count_if(D.data(), D.data() + m_p,
+				 std::bind2nd(std::greater<double>(), D[0] * prec));
+    m_coef       = UDV.solve(y);
+    m_fitted     = X * m_coef;
+    MatrixXd VDi = UDV.matrixV() * DiagType(UDV.singularValues().array().inverse().matrix());
+    m_unsc       = m_unsc.selfadjointView<Eigen::Lower>().rankUpdate(VDi);
+    for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
+}
+
 enum {ColPivQR_t = 0, QR_t, LDLT_t, LLT_t, SVD_t, SymmEigen_t};
 
 static inline lm do_lm(const MMatrixXd &X, const MVectorXd &y, int type)
@@ -98,8 +111,8 @@
 	return LLT(X, y);
     case LDLT_t:
 	return LDLT(X, y);
-    // case SVD_t:
-    // 	return SVD(X, y);
+    case SVD_t:
+     	return SVD(X, y);
     // case SymmEigen_t:
     // 	return SymmEigen(X, y);
     }
@@ -150,85 +163,51 @@
 }
 
 #if 0
-extern "C" SEXP fastLmBench(SEXP Xs, SEXP ys) {
-    using namespace Eigen;
-    using namespace Rcpp;
-    try {
-	NumericMatrix X(Xs);
-	NumericVector y(ys);
-	Indices       n = X.nrow(), p = X.ncol();
-	int          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) {
+extern "C" SEXP fastLm(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();
+	size_t          n = X.nrow(), p = X.ncol(), r = 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);
+	const Map<MatrixXd>   A(X.begin(), n, p); // shares storage
+	const Map<VectorXd>   b(y.begin(), n);
+	
+	VectorXd     coef, fitted, unsc;
+	
+	LDLT<MatrixXd> Ch(SelfAdjointView<MatrixXd, Upper>(MatrixXd::Zero(p, p))
+			  .rankUpdate(A.adjoint()));
+	ArrayXd      diag = Ch.vectorD().array();            
+	double       dmax = *std::max_element(diag.data(), diag.data() + p),
+	    dmin = *std::min_element(diag.data(), diag.data() + p);
+	if (dmin <= 0. || dmin/dmax < 100 * std::numeric_limits<double>.epsilon()) {
+	    ColPivHouseholderQR<MatrixXd> mqr(A);      // decompose the model matrix
+	    r                  = mqr.rank();
+	    df                 = n - r;
+ 	    MatrixXd    Atrunc = (A * mqr.colsPermutation()).leftCols(r);
+ 	    HouseholderQR<MatrixXd> QR(Atrunc);
+ 	    VectorXd coefTrunc = QR.solve(b);
+ 	    fitted             = Atrunc * coefTrunc;
+ 	    coef.topRows(r)    = coefTrunc;
+ 	    std::fill(coef.data() + r, coef.data() + p, std::numeric_limits<double>::quiet_NaN());
+ 	    coef = mqr.colsPermutation() * coef;
+	    MatrixXd         R = QR.matrixQR().topRows(r);
+ 	    MatrixXd      Rinv = TriangularView<MatrixXd, Upper>(R).solve(MatrixXd::Identity(r, r));
+	    unsc.topRows(r)    = Rinv.rowwise.norm();
+ 	    std::fill(unsc.data() + r, unsc.data() + p, std::numeric_limits<double>::quiet_NaN());
+	} else {
+	    coef               = Ch.solve(A.adjoint() * b);
+	    fitted             = A * coef;
+	    unsc               = (Ch.solve(MatrixXd::Identity(p, p)).diagonal().array()).sqrt();
+	}
+	VectorXd    resid = b - fitted;
+	double          s = std::sqrt(resid.squaredNorm()/df);
+	VectorXf       se = unsc.array() * s;
 
 	LDLT           Ch(MatrixXd::Zero(p, p).selfAdjointView(Eigen::Lower).rankUpdate(A.adjoint()));
 	VectorXd     coef = Ch.solve(A.adjoint() * b);
@@ -239,7 +218,12 @@
 	std::copy(se.data(), se.data() + p, Rse.begin());
 			    
 	return List::create(_["coefficients"] = coef,
-			    _["se"]           = Rse,
+			    _["stderr"]       = se,
+                            _["rank"]         = r,
+                            _["s"]            = s,
+                            _["perm"]         = perm,
+			    _["residuals"]    = resid,
+			    _["fitted"]       = fitted,
 			    _["df.residual"]  = df);
     } catch( std::exception &ex ) {
 	forward_exception_to_r( ex );
@@ -249,36 +233,36 @@
     return R_NilValue; // -Wall
 }
 
-extern "C" SEXP fastLmSVD(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");
+// extern "C" SEXP fastLmSVD(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);
+// 	MatrixXd        A = Map<MatrixXd>(X.begin(), n, p); // shares storage
+// 	VectorXd        b = Map<VectorXd>(y.begin(), n);
 
-	VectorXd     coef = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);
-//	double         s2 = (b - A*coef).squaredNorm()/df;
+// 	VectorXd     coef = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);
+// //	double         s2 = (b - A*coef).squaredNorm()/df;
 
-	// ArrayXd        se = (Ch.solve(MatrixXd::Identity(p, p)).diagonal().array() * s2).sqrt();
-	// NumericVector Rse(p);	// should define a wrap method for ArrayXd, ArrayXXd, etc.
-	// std::copy(se.data(), se.data() + p, Rse.begin());
+// 	// ArrayXd        se = (Ch.solve(MatrixXd::Identity(p, p)).diagonal().array() * s2).sqrt();
+// 	// NumericVector Rse(p);	// should define a wrap method for ArrayXd, ArrayXXd, etc.
+// 	// 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
-}
+// 	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
+// }
 
 #endif

Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h	2011-06-21 02:33:41 UTC (rev 3096)
+++ pkg/RcppEigen/src/fastLm.h	2011-06-21 22:17:36 UTC (rev 3097)
@@ -26,6 +26,7 @@
 typedef Eigen::ColPivHouseholderQR<MatrixXd>           PivQRType;
 typedef Eigen::HouseholderQR<MatrixXd>                 QRType;
 typedef Eigen::JacobiSVD<MatrixXd>                     SVDType;
+typedef Eigen::DiagonalMatrix<double,Eigen::Dynamic>   DiagType;
 				// Types derived from decompositions
 typedef typename PivQRType::PermutationType            PermutationType;
 typedef typename PermutationType::IndicesType          IndicesType;
@@ -70,6 +71,11 @@
     QR(const MMatrixXd&, const MVectorXd&);
 };
 
+class SVD : public lm {
+public:
+    SVD(const MMatrixXd&, const MVectorXd&);
+};
+
 extern "C" SEXP fastLm(SEXP Xs, SEXP ys, SEXP types);
 
 #endif



More information about the Rcpp-commits mailing list