Skip to content

Commit 83e882d

Browse files
dmbatesdmbates
committed
Use a macro for XtX temporarily. Clean up code.
1 parent b8da05b commit 83e882d

File tree

2 files changed

+42
-31
lines changed

2 files changed

+42
-31
lines changed

‎src/fastLm.cpp‎

Lines changed: 39 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -49,13 +49,19 @@ namespace lmsol {
4949
m_usePrescribedThreshold(false) {
5050
}
5151

52-
lm& lm::setThreshold(const RealScalar& threshold)
53-
{
52+
lm& lm::setThreshold(const RealScalar& threshold) {
5453
m_usePrescribedThreshold = true;
5554
m_prescribedThreshold = threshold;
5655
return *this;
5756
}
5857

58+
SelfAdjointView<MatrixXd,Lower> lm::XtX() const {
59+
return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().
60+
rankUpdate(m_X.adjoint());
61+
}
62+
// For some reason that function returning a SelfAdjointView encounters a bad_alloc error
63+
// Use a macro for the time being
64+
#define XtX MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().rankUpdate(m_X.adjoint())
5965
/** Returns the threshold that will be used by certain methods such as rank().
6066
*
6167
* The default value comes from experimenting (see "LU precision
@@ -64,13 +70,13 @@ namespace lmsol {
6470
*
6571
* @return The user-prescribed threshold or the default.
6672
*/
67-
RealScalar lm::threshold() const
68-
{
73+
RealScalar lm::threshold() const {
6974
return m_usePrescribedThreshold ? m_prescribedThreshold
7075
: numeric_limits<double>::epsilon() * m_p;
7176
}
7277

73-
ColPivQR::ColPivQR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
78+
ColPivQR::ColPivQR(const Map<MatrixXd> &X, const Map<VectorXd> &y)
79+
: lm(X, y) {
7480
ColPivHouseholderQR<MatrixXd> PQR(X); // decompose the model matrix
7581
Permutation Pmat(PQR.colsPermutation());
7682
m_r = PQR.rank();
@@ -89,7 +95,7 @@ namespace lmsol {
8995
m_coef.head(m_r) = Rinv * effects.head(m_r);
9096
m_coef = Pmat * m_coef;
9197
// create fitted values from effects
92-
// (can't use X * m_coef when X is rank-deficient)
98+
// (can't use X*m_coef if X is rank-deficient)
9399
effects.tail(m_df).setZero();
94100
m_fitted = PQR.householderQ() * effects;
95101
m_se.head(m_r) = Rinv.rowwise().norm();
@@ -98,36 +104,40 @@ namespace lmsol {
98104

99105
QR::QR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
100106
HouseholderQR<MatrixXd> QR(X);
101-
m_se = QR.matrixQR().topRows(m_p).
102-
triangularView<Upper>().solve(I_p()).rowwise().norm();
103107
m_coef = QR.solve(y);
104108
m_fitted = X * m_coef;
109+
m_se = QR.matrixQR().topRows(m_p).
110+
triangularView<Upper>().solve(I_p()).rowwise().norm();
105111
}
106112

107113

108114
Llt::Llt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
109-
LLT<MatrixXd> Ch(XtX());
115+
LLT<MatrixXd> Ch(XtX);
110116
m_coef = Ch.solve(X.adjoint() * y);
111117
m_fitted = X * m_coef;
112118
m_se = Ch.matrixL().solve(I_p()).colwise().norm();
113119
}
114120

115-
Ldlt::Ldlt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
116-
LDLT<MatrixXd> Ch(XtX());
117-
m_coef = Ch.solve(X.adjoint() * y);
118-
m_fitted = X * m_coef;
119-
m_se = Ch.solve(I_p()).diagonal().array().sqrt();
120-
}
121-
122-
inline DiagonalMatrix<double, Dynamic> Dplus(const ArrayXd& D, Index r, bool rev=false) {
121+
inline DiagonalMatrix<double, Dynamic> Dplus(const ArrayXd& D,
122+
Index r, bool rev=false) {
123123
VectorXd Di(VectorXd::Constant(D.size(), 0.));
124124
if (rev) Di.tail(r) = D.tail(r).inverse();
125125
else Di.head(r) = D.head(r).inverse();
126126
return DiagonalMatrix<double, Dynamic>(Di);
127127
}
128128

129+
Ldlt::Ldlt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
130+
LDLT<MatrixXd> Ch(XtX);
131+
ArrayXd D(Ch.vectorD());
132+
m_r = (D > D.maxCoeff() * threshold()).count();
133+
// FIXME: work out how to use Dplus with elements of D unsorted.
134+
m_coef = Ch.solve(X.adjoint() * y);
135+
m_fitted = X * m_coef;
136+
m_se = Ch.solve(I_p()).diagonal().array().sqrt();
137+
}
138+
129139
SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
130-
JacobiSVD<MatrixXd> UDV(X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV));
140+
JacobiSVD<MatrixXd> UDV(X.jacobiSvd(ComputeThinU|ComputeThinV));
131141
ArrayXd D(UDV.singularValues());
132142
m_r = (D > D[0] * threshold()).count();
133143
m_df = m_n - m_r;
@@ -137,16 +147,17 @@ namespace lmsol {
137147
m_se = VDi.rowwise().norm();
138148
}
139149

140-
SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
141-
SelfAdjointEigenSolver<MatrixXd> eig(XtX());
142-
ArrayXd D(eig.eigenvalues());
143-
m_r = (D > D[m_p - 1] * threshold()).count();
144-
D = D.sqrt();
145-
m_df = m_n - m_r;
146-
MatrixXd VDi(eig.eigenvectors() * Dplus(D, m_r, true));
147-
m_coef = VDi * VDi.adjoint() * X.adjoint() * y;
148-
m_fitted = X * m_coef;
149-
m_se = VDi.rowwise().norm();
150+
SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y)
151+
: lm(X, y) {
152+
SelfAdjointEigenSolver<MatrixXd> eig(XtX);
153+
ArrayXd D(eig.eigenvalues());
154+
m_r = (D > D[m_p - 1] * threshold()).count();
155+
D = D.sqrt();
156+
m_df = m_n - m_r;
157+
MatrixXd VDi(eig.eigenvectors() * Dplus(D, m_r, true));
158+
m_coef = VDi * VDi.adjoint() * X.adjoint() * y;
159+
m_fitted = X * m_coef;
160+
m_se = VDi.rowwise().norm();
150161
}
151162

152163
enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t};

‎src/fastLm.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,8 @@
2727
namespace lmsol {
2828
using Eigen::ArrayXd;
2929
using Eigen::ColPivHouseholderQR;
30+
using Eigen::ComputeThinU;
31+
using Eigen::ComputeThinV;
3032
using Eigen::DiagonalMatrix;
3133
using Eigen::Dynamic;
3234
using Eigen::HouseholderQR;
@@ -69,9 +71,7 @@ namespace lmsol {
6971
return MatrixXd::Identity(m_p, m_p);
7072
}
7173
RealScalar threshold() const;
72-
SelfAdjointView<MatrixXd,Lower> XtX() const {
73-
return MatrixXd(m_p, m_p).setZero().selfadjointView<Lower>().rankUpdate(m_X.adjoint());
74-
};
74+
SelfAdjointView<MatrixXd,Lower> XtX() const;
7575
const VectorXd& se() const {return m_se;}
7676
const VectorXd& coef() const {return m_coef;}
7777
const VectorXd& fitted() const {return m_fitted;}

0 commit comments

Comments
 (0)