Skip to content

Commit eab75d4

Browse files
dmbatesdmbates
committed
Added another fastLm method using the Lapack subroutine dgesdd for the SVD
1 parent 78f5e14 commit eab75d4

File tree

5 files changed

+46
-18
lines changed

5 files changed

+46
-18
lines changed

‎inst/examples/lmBenchmark.R‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,9 @@ exprs$lm.fit <- expression(stats::lm.fit(mm, y))
2121
exprs$PivQR <- expression(.Call("fastLm", mm, y, 0L, PACKAGE="RcppEigen"))
2222
## LDLt Cholesky decomposition with rank detection
2323
exprs$LDLt <- expression(.Call("fastLm", mm, y, 2L, PACKAGE="RcppEigen"))
24-
## SVD
24+
## SVD using the Lapack subroutine dgesdd and Eigen support
25+
exprs$GESDD <- expression(.Call("fastLm", mm, y, 6L, PACKAGE="RcppEigen"))
26+
## SVD (the JacobiSVD class from Eigen)
2527
exprs$SVD <- expression(.Call("fastLm", mm, y, 4L, PACKAGE="RcppEigen"))
2628
## eigenvalues and eigenvectors of X'X
2729
exprs$SymmEig <- expression(.Call("fastLm", mm, y, 5L, PACKAGE="RcppEigen"))

‎src/Makevars‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,5 +7,4 @@ PKG_CXXFLAGS=-I../inst/include -DNDEBUG
77
## to enable error checking. However, checking for errors is flagged
88
## in R CMD check as an undesirable practice.
99

10-
PKG_LIBS=`$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"`
11-
10+
PKG_LIBS=`$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)

‎src/Makevars.win‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
## This assumes that we can call Rscript to ask Rcpp about its locations
22
## Use the R_HOME indirection to support installations of multiple R version
3-
PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
3+
PKG_LIBS = $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()") $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
44

55
PKG_CPPFLAGS = -I../inst/include -I. -DNDEBUG
66

‎src/fastLm.cpp‎

Lines changed: 35 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
// <http://www.gnu.org/licenses/>.
2222

2323
#include "fastLm.h"
24+
#include <R_ext/Lapack.h>
2425

2526
namespace lmsol {
2627
using Rcpp::_;
@@ -54,18 +55,14 @@ namespace lmsol {
5455
return *this;
5556
}
5657

57-
ArrayXd lm::Dplus(const ArrayXd& D) {
58-
ArrayXd Di(m_p);
59-
double comp(D.maxCoeff() * threshold());
60-
for (int j = 0; j < m_p; ++j) Di[j] = (D[j] < comp) ? 0. : 1./D[j];
61-
m_r = (Di != 0.).count();
62-
return Di;
58+
inline ArrayXd lm::Dplus(const ArrayXd& d) {
59+
ArrayXd di(d.size());
60+
double comp(d.maxCoeff() * threshold());
61+
for (int j = 0; j < d.size(); ++j) di[j] = (d[j] < comp) ? 0. : 1./d[j];
62+
m_r = (di != 0.).count();
63+
return di;
6364
}
6465

65-
MatrixXd lm::I_p() const {
66-
return MatrixXd::Identity(m_p, m_p);
67-
}
68-
6966
MatrixXd lm::XtX() const {
7067
return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
7168
rankUpdate(m_X.adjoint());
@@ -138,6 +135,31 @@ namespace lmsol {
138135
m_se = Ch.solve(I_p()).diagonal().array().sqrt();
139136
}
140137

138+
int gesdd(MatrixXd& A, ArrayXd& S, MatrixXd& Vt) {
139+
int info, mone = -1, m = A.rows(), n = A.cols();
140+
std::vector<int> iwork(8 * n);
141+
double wrk;
142+
if (m < n || S.size() != n || Vt.rows() != n || Vt.cols() != n)
143+
throw std::invalid_argument("dimension mismatch in gesvd");
144+
F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
145+
&m, Vt.data(), &n, &wrk, &mone, &iwork[0], &info);
146+
int lwork(wrk);
147+
std::vector<double> work(lwork);
148+
F77_CALL(dgesdd)("O", &m, &n, A.data(), &m, S.data(), A.data(),
149+
&m, Vt.data(), &n, &work[0], &lwork, &iwork[0], &info);
150+
return info;
151+
}
152+
153+
GESDD::GESDD(const Map<MatrixXd>& X, const Map<VectorXd> &y) : lm(X, y) {
154+
MatrixXd U(X), Vt(m_p, m_p);
155+
ArrayXd S(m_p);
156+
if (gesdd(U, S, Vt)) throw std::runtime_error("error in gesdd");
157+
MatrixXd VDi(Vt.adjoint() * Dplus(S).matrix().asDiagonal());
158+
m_coef = VDi * U.adjoint() * y;
159+
m_fitted = X * m_coef;
160+
m_se = VDi.rowwise().norm();
161+
}
162+
141163
SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
142164
JacobiSVD<MatrixXd> UDV(X.jacobiSvd(ComputeThinU|ComputeThinV));
143165
MatrixXd VDi(UDV.matrixV() *
@@ -157,7 +179,7 @@ namespace lmsol {
157179
m_se = VDi.rowwise().norm();
158180
}
159181

160-
enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t};
182+
enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t, GESDD_t};
161183

162184
static inline lm do_lm(const Map<MatrixXd> &X, const Map<VectorXd> &y, int type) {
163185
switch(type) {
@@ -173,6 +195,8 @@ namespace lmsol {
173195
return SVD(X, y);
174196
case SymmEigen_t:
175197
return SymmEigen(X, y);
198+
case GESDD_t:
199+
return GESDD(X, y);
176200
}
177201
throw invalid_argument("invalid type");
178202
return ColPivQR(X, y); // -Wall

‎src/fastLm.h‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -29,8 +29,6 @@ namespace lmsol {
2929
using Eigen::ColPivHouseholderQR;
3030
using Eigen::ComputeThinU;
3131
using Eigen::ComputeThinV;
32-
using Eigen::DiagonalMatrix;
33-
using Eigen::Dynamic;
3432
using Eigen::HouseholderQR;
3533
using Eigen::JacobiSVD;
3634
using Eigen::LDLT;
@@ -66,7 +64,7 @@ namespace lmsol {
6664
lm(const Map<MatrixXd>&, const Map<VectorXd>&);
6765

6866
ArrayXd Dplus(const ArrayXd& D);
69-
MatrixXd I_p() const;
67+
MatrixXd I_p() const {return MatrixXd::Identity(m_p, m_p);}
7068
MatrixXd XtX() const;
7169

7270
// setThreshold and threshold are based on ColPivHouseholderQR methods
@@ -98,6 +96,11 @@ namespace lmsol {
9896
QR(const Map<MatrixXd>&, const Map<VectorXd>&);
9997
};
10098

99+
class GESDD : public lm {
100+
public:
101+
GESDD(const Map<MatrixXd>&, const Map<VectorXd>&);
102+
};
103+
101104
class SVD : public lm {
102105
public:
103106
SVD(const Map<MatrixXd>&, const Map<VectorXd>&);

0 commit comments

Comments
 (0)