[Rcpp-commits] r3251 - pkg/RcppEigen/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 31 23:20:31 CET 2011
Author: dmbates
Date: 2011-10-31 23:20:31 +0100 (Mon, 31 Oct 2011)
New Revision: 3251
Modified:
pkg/RcppEigen/src/fastLm.cpp
pkg/RcppEigen/src/fastLm.h
Log:
Extend rank calculations to other decompositions
Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp 2011-10-31 19:48:01 UTC (rev 3250)
+++ pkg/RcppEigen/src/fastLm.cpp 2011-10-31 22:20:31 UTC (rev 3251)
@@ -43,7 +43,6 @@
m_p(X.cols()),
m_coef(VectorXd::Constant(m_p, ::NA_REAL)),
m_r(::NA_INTEGER),
- m_df(m_n - m_p),
m_fitted(m_n),
m_se(VectorXd::Constant(m_p, ::NA_REAL)),
m_usePrescribedThreshold(false) {
@@ -55,6 +54,18 @@
return *this;
}
+ ArrayXd lm::Dplus(const ArrayXd& D) {
+ ArrayXd Di(m_p);
+ double comp(D.maxCoeff() * threshold());
+ for (int j = 0; j < m_p; ++j) Di[j] = (D[j] < comp) ? 0. : 1./D[j];
+ m_r = (Di != 0.).count();
+ return Di;
+ }
+
+ MatrixXd lm::I_p() const {
+ return MatrixXd::Identity(m_p, m_p);
+ }
+
MatrixXd lm::XtX() const {
return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
rankUpdate(m_X.adjoint());
@@ -85,7 +96,6 @@
triangularView<Upper>().solve(I_p()).rowwise().norm();
return;
}
- m_df = m_n - m_r;
MatrixXd Rinv(PQR.matrixQR().topLeftCorner(m_r, m_r).
triangularView<Upper>().
solve(MatrixXd::Identity(m_r, m_r)));
@@ -94,7 +104,7 @@
m_coef = Pmat * m_coef;
// create fitted values from effects
// (can't use X*m_coef if X is rank-deficient)
- effects.tail(m_df).setZero();
+ effects.tail(m_n - m_r).setZero();
m_fitted = PQR.householderQ() * effects;
m_se.head(m_r) = Rinv.rowwise().norm();
m_se = Pmat * m_se;
@@ -116,19 +126,13 @@
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().selfadjointView<Lower>());
- ArrayXd D(Ch.vectorD());
- m_r = (D > D.maxCoeff() * threshold()).count();
- // FIXME: work out how to use Dplus with elements of D unsorted.
+ Dplus(Ch.vectorD()); // to set the rank
+//FIXME: Check on the permutation in the LDLT and incorporate it in
+//the coefficients and the standard error computation.
+// m_coef = Ch.matrixL().adjoint().
+// solve(Dplus(D) * Ch.matrixL().solve(X.adjoint() * y));
m_coef = Ch.solve(X.adjoint() * y);
m_fitted = X * m_coef;
m_se = Ch.solve(I_p()).diagonal().array().sqrt();
@@ -136,10 +140,8 @@
SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
JacobiSVD<MatrixXd> UDV(X.jacobiSvd(ComputeThinU|ComputeThinV));
- ArrayXd D(UDV.singularValues());
- m_r = (D > D[0] * threshold()).count();
- m_df = m_n - m_r;
- MatrixXd VDi(UDV.matrixV() * Dplus(D, m_r));
+ MatrixXd VDi(UDV.matrixV() *
+ DiagonalMatrix(Dplus(UDV.singularValues().array()).matrix()));
m_coef = VDi * UDV.matrixU().adjoint() * y;
m_fitted = X * m_coef;
m_se = VDi.rowwise().norm();
@@ -148,14 +150,11 @@
SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y)
: lm(X, y) {
SelfAdjointEigenSolver<MatrixXd> eig(XtX().selfadjointView<Lower>());
- 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();
+ MatrixXd VDi(eig.eigenvectors() *
+ DiagonalMatrix(Dplus(eig.eigenvalues().array()).sqrt().matrix());
+ 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};
@@ -199,14 +198,16 @@
}
VectorXd resid = y - ans.fitted();
- double s = resid.norm() / std::sqrt(double(ans.df()));
+ int rank = ans.rank();
+ int df = (rank == ::NA_INTEGER) ? n - X.cols() : n - rank;
+ double s = resid.norm() / std::sqrt(double(df));
// Create the standard errors
VectorXd se = s * ans.se();
return List::create(_["coefficients"] = coef,
_["se"] = se,
- _["rank"] = ans.rank(),
- _["df.residual"] = ans.df(),
+ _["rank"] = rank,
+ _["df.residual"] = df,
_["residuals"] = resid,
_["s"] = s,
_["fitted.values"] = ans.fitted());
Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h 2011-10-31 19:48:01 UTC (rev 3250)
+++ pkg/RcppEigen/src/fastLm.h 2011-10-31 22:20:31 UTC (rev 3251)
@@ -58,7 +58,6 @@
Index m_p; /**< number of columns of X */
VectorXd m_coef; /**< coefficient vector */
int m_r; /**< computed rank or NA_INTEGER */
- int m_df; /**< residual degrees of freedom */
VectorXd m_fitted; /**< vector of fitted values */
VectorXd m_se; /**< standard errors */
RealScalar m_prescribedThreshold; /**< user specified tolerance */
@@ -66,18 +65,17 @@
public:
lm(const Map<MatrixXd>&, const Map<VectorXd>&);
- // setThreshold and threshold are based on ColPivHouseholderQR methods
- MatrixXd I_p() const {
- return MatrixXd::Identity(m_p, m_p);
- }
- RealScalar threshold() const;
- MatrixXd XtX() const;
- const VectorXd& se() const {return m_se;}
- const VectorXd& coef() const {return m_coef;}
- const VectorXd& fitted() const {return m_fitted;}
- int df() const {return m_df;}
- int rank() const {return m_r;}
- lm& setThreshold(const RealScalar&);
+ ArrayXd Dplus(const ArrayXd& D);
+ MatrixXd I_p() const;
+ MatrixXd XtX() const;
+
+ // setThreshold and threshold are based on ColPivHouseholderQR methods
+ RealScalar threshold() const;
+ const VectorXd& se() const {return m_se;}
+ const VectorXd& coef() const {return m_coef;}
+ const VectorXd& fitted() const {return m_fitted;}
+ int rank() const {return m_r;}
+ lm& setThreshold(const RealScalar&);
};
class ColPivQR : public lm {
More information about the Rcpp-commits
mailing list