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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Oct 26 21:01:02 CEST 2011


Author: dmbates
Date: 2011-10-26 21:01:02 +0200 (Wed, 26 Oct 2011)
New Revision: 3228

Modified:
   pkg/RcppEigen/src/fastLm.cpp
   pkg/RcppEigen/src/fastLm.h
Log:
Use a macro for XtX temporarily.  Clean up code.


Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp	2011-10-26 18:59:34 UTC (rev 3227)
+++ pkg/RcppEigen/src/fastLm.cpp	2011-10-26 19:01:02 UTC (rev 3228)
@@ -49,13 +49,19 @@
 	  m_usePrescribedThreshold(false) {
     }
 
-    lm& lm::setThreshold(const RealScalar& threshold)
-    {
+    lm& lm::setThreshold(const RealScalar& threshold) {
 	m_usePrescribedThreshold = true;
 	m_prescribedThreshold = threshold;
 	return *this;
     }
 
+    SelfAdjointView<MatrixXd,Lower> lm::XtX() const {
+	return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
+	    rankUpdate(m_X.adjoint());
+    }
+// For some reason that function returning a SelfAdjointView encounters a bad_alloc error
+// Use a macro for the time being
+#define XtX MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().rankUpdate(m_X.adjoint())
     /** Returns the threshold that will be used by certain methods such as rank().
      * 
      *  The default value comes from experimenting (see "LU precision
@@ -64,13 +70,13 @@
      *
      *  @return The user-prescribed threshold or the default.
      */
-    RealScalar lm::threshold() const
-    {
+    RealScalar lm::threshold() const {
 	return m_usePrescribedThreshold ? m_prescribedThreshold
 	    : numeric_limits<double>::epsilon() * m_p; 
     }
 
-    ColPivQR::ColPivQR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+    ColPivQR::ColPivQR(const Map<MatrixXd> &X, const Map<VectorXd> &y)
+	: lm(X, y) {
 	ColPivHouseholderQR<MatrixXd> PQR(X); // decompose the model matrix
 	Permutation                  Pmat(PQR.colsPermutation());
 	m_r                               = PQR.rank();
@@ -89,7 +95,7 @@
 	m_coef.head(m_r)                  = Rinv * effects.head(m_r);
 	m_coef                            = Pmat * m_coef;
 				// create fitted values from effects
-				// (can't use X * m_coef when X is rank-deficient)
+				// (can't use X*m_coef if X is rank-deficient)
 	effects.tail(m_df).setZero();
 	m_fitted                          = PQR.householderQ() * effects;
 	m_se.head(m_r)                    = Rinv.rowwise().norm();
@@ -98,36 +104,40 @@
     
     QR::QR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
 	HouseholderQR<MatrixXd> QR(X);
+	m_coef                     = QR.solve(y);
+	m_fitted                   = X * m_coef;
 	m_se                       = QR.matrixQR().topRows(m_p).
 	    triangularView<Upper>().solve(I_p()).rowwise().norm();
-	m_coef                     = QR.solve(y);
-	m_fitted                   = X * m_coef;
     }
     
     
     Llt::Llt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	LLT<MatrixXd>  Ch(XtX());
+	LLT<MatrixXd>  Ch(XtX);
 	m_coef            = Ch.solve(X.adjoint() * y);
 	m_fitted          = X * m_coef;
 	m_se              = Ch.matrixL().solve(I_p()).colwise().norm();
     }
     
+    inline DiagonalMatrix<double, Dynamic> Dplus(const ArrayXd& D,
+						 Index r, bool rev=false) {
+	VectorXd   Di(VectorXd::Constant(D.size(), 0.));
+	if (rev) Di.tail(r)  = D.tail(r).inverse();
+	else Di.head(r)      = D.head(r).inverse();
+	return DiagonalMatrix<double, Dynamic>(Di);
+    }
+
     Ldlt::Ldlt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	LDLT<MatrixXd> Ch(XtX());
+	LDLT<MatrixXd> Ch(XtX);
+	ArrayXd         D(Ch.vectorD());
+	m_r               = (D > D.maxCoeff() * threshold()).count();
+	// FIXME: work out how to use Dplus with elements of D unsorted.
 	m_coef            = Ch.solve(X.adjoint() * y);
 	m_fitted          = X * m_coef;
 	m_se              = Ch.solve(I_p()).diagonal().array().sqrt();
     }
     
-    inline DiagonalMatrix<double, Dynamic> Dplus(const ArrayXd& D, Index r, bool rev=false) {
-	VectorXd   Di(VectorXd::Constant(D.size(), 0.));
-	if (rev) Di.tail(r)  = D.tail(r).inverse();
-	else Di.head(r)      = D.head(r).inverse();
-	return DiagonalMatrix<double, Dynamic>(Di);
-    }
-
     SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	JacobiSVD<MatrixXd>  UDV(X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV));
+	JacobiSVD<MatrixXd>  UDV(X.jacobiSvd(ComputeThinU|ComputeThinV));
 	ArrayXd                D(UDV.singularValues());
 	m_r                      = (D > D[0] * threshold()).count();
 	m_df                     = m_n - m_r;
@@ -137,16 +147,17 @@
 	m_se                     = VDi.rowwise().norm();
     }
 
-    SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
-	SelfAdjointEigenSolver<MatrixXd> eig(XtX());
-	ArrayXd                            D(eig.eigenvalues());
-	m_r                                  = (D > D[m_p - 1] * threshold()).count();
-	D                                    = D.sqrt();
-	m_df                                 = m_n - m_r;
-	MatrixXd                         VDi(eig.eigenvectors() * Dplus(D, m_r, true));
-	m_coef                               = VDi * VDi.adjoint() * X.adjoint() * y;
-	m_fitted                             = X * m_coef;
-	m_se                                 = VDi.rowwise().norm();
+    SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y)
+	: lm(X, y) {
+	SelfAdjointEigenSolver<MatrixXd> eig(XtX);
+	ArrayXd                      D(eig.eigenvalues());
+	m_r                            = (D > D[m_p - 1] * threshold()).count();
+	D                              = D.sqrt();
+	m_df                           = m_n - m_r;
+	MatrixXd                   VDi(eig.eigenvectors() * Dplus(D, m_r, true));
+	m_coef                         = VDi * VDi.adjoint() * X.adjoint() * y;
+	m_fitted                       = X * m_coef;
+	m_se                           = VDi.rowwise().norm();
     }
 
     enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t};

Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h	2011-10-26 18:59:34 UTC (rev 3227)
+++ pkg/RcppEigen/src/fastLm.h	2011-10-26 19:01:02 UTC (rev 3228)
@@ -27,6 +27,8 @@
 namespace lmsol {
     using Eigen::ArrayXd;
     using Eigen::ColPivHouseholderQR;
+    using Eigen::ComputeThinU;
+    using Eigen::ComputeThinV;
     using Eigen::DiagonalMatrix;
     using Eigen::Dynamic;
     using Eigen::HouseholderQR;
@@ -69,9 +71,7 @@
 	    return MatrixXd::Identity(m_p, m_p);
 	}
 	RealScalar                threshold() const;
-	SelfAdjointView<MatrixXd,Lower> XtX() const {
-	    return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().rankUpdate(m_X.adjoint());
-	};
+	SelfAdjointView<MatrixXd,Lower> XtX() const;
 	const VectorXd&                  se() const {return m_se;}
 	const VectorXd&                coef() const {return m_coef;}
 	const VectorXd&              fitted() const {return m_fitted;}



More information about the Rcpp-commits mailing list