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

[noreply at r-forge.r-project.org]

Author: dmbates
Date: 2011-06-13 23:02:30 +0200 (Mon, 13 Jun 2011)
New Revision: 3057

Modified:
   pkg/RcppEigen/src/fastLm.cpp
Log:
Handle the rank-deficient case.


Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp	2011-06-13 21:01:50 UTC (rev 3056)
+++ pkg/RcppEigen/src/fastLm.cpp	2011-06-13 21:02:30 UTC (rev 3057)
@@ -25,7 +25,6 @@
     using namespace Eigen;
     using namespace Rcpp;
     try {
-//	static double NaN = std::numeric_limits<double>::quiet_NaN();
 	NumericMatrix X(Xs);
 	NumericVector y(ys);
 	size_t n = X.nrow(), p = X.ncol();
@@ -36,19 +35,36 @@
 	VectorXd b = Map<VectorXd>(y.begin(), n);
     
 	ColPivHouseholderQR<MatrixXd> mqr(A);      // decompose the model matrix
-	if (!mqr.isInjective())
-	    throw std::invalid_argument("X is rank deficient");
-	size_t               r = mqr.rank();
-	ptrdiff_t           df = n - r;
-	VectorXd          coef = mqr.solve(b); 
-	VectorXd        fitted = A * coef;
-	VectorXd           res = b - fitted;
-	double               s = std::sqrt(res.squaredNorm()/df);
+	size_t        r = mqr.rank();
+	ptrdiff_t    df = n - r;
+	VectorXd   coef(p);
+	VectorXd fitted(n);
+	VectorXd     se(p);
+	MatrixXd   Rinv;
+	if (r < p) {		// Handle the rank-deficient case
+	    MatrixXd Atrunc = (A * mqr.colsPermutation()).leftCols(r);
+	    HouseholderQR<MatrixXd> QR(Atrunc);
+	    VectorXd coefTrunc = QR.solve(b);
+	    fitted = Atrunc * coefTrunc;
+	    MatrixXd   R = QR.matrixQR().topLeftCorner(r, r);
+	    Rinv = TriangularView<MatrixXd, Upper>(R).solve(MatrixXd::Identity(r, r));
+	    coef.topRows(r) = coefTrunc;
+	    std::fill(coef.data() + r, coef.data() + p,
+		      std::numeric_limits<double>::quiet_NaN());
+	    coef = mqr.colsPermutation() * coef;
+	    std::fill(se.data() + r, se.data() + p,
+		      std::numeric_limits<double>::quiet_NaN());
+	} else {
+	    coef = mqr.solve(b); 
+	    fitted = A * coef;
+	    MatrixXd  R = mqr.matrixQR().topRows(p);
+	    Rinv = TriangularView<MatrixXd, Upper>(R).solve(MatrixXd::Identity(p, p));
+	}
+	VectorXd      res = b - fitted;
+	double          s = std::sqrt(res.squaredNorm()/df);
+	se.topRows(r) = Rinv.rowwise().norm() * s;
+	se = mqr.colsPermutation() * se;
 
-	MatrixXd  R = mqr.matrixQR().topLeftCorner(r, r); // need to create this explicitly
-	MatrixXd  Rinv(TriangularView<MatrixXd, Upper>(R).solve(MatrixXd::Identity(r, r)));
-	VectorXd se = mqr.colsPermutation() * Rinv.rowwise().norm() * s; // standard errors
-
 	return List::create(_["coefficients"]  = coef,
 			    _["rank"]          = r,
 			    _["df"]            = df,