12345qiupeng
/
gtsam
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

in-place house vector (was about 10% of all mallocs in Urban)

release/4.3a0
Frank Dellaert 2010-03-11 15:04:31 +00:00
parent b726e8e5e2
commit 67fb7fa9ff
3 changed files with 19 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
cpp/Matrix.cpp
View File

@ -641,9 +641,9 @@ inline void householder_manual(Matrix &A, size_t k) {
for(size_t r = j ; r < m; r++) for(size_t r = j ; r < m; r++)
xjm(r-j) = A(r,j); xjm(r-j) = A(r,j);
// calculate the Householder vector // calculate the Householder vector, in place
double beta; Vector vjm; double beta = houseInPlace(xjm);
boost::tie(beta,vjm) = house(xjm); Vector& vjm = xjm;
// do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w' // do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
householder_update(A, j, beta, vjm); householder_update(A, j, beta, vjm);

20
cpp/Vector.cpp
View File

@ -327,18 +327,18 @@ namespace gtsam {
} }
/* ************************************************************************* */ /* ************************************************************************* */
pair<double, Vector > house(const Vector &x) { // imperative version, pass in x
const double x0 = x(0); double houseInPlace(Vector &v) {
const double x0 = v(0);
const double x02 = x0*x0; const double x02 = x0*x0;
// old code - GSL verison was actually a bit slower // old code - GSL verison was actually a bit slower
const double sigma = inner_prod(trans(x),x) - x02; const double sigma = inner_prod(v,v) - x02;
double beta = 0.0;
Vector v(x); v(0) = 1.0; v(0) = 1.0;
if( sigma == 0.0 ) if( sigma == 0.0 )
beta = 0.0; return 0.0;
else { else {
double mu = sqrt(x02 + sigma); double mu = sqrt(x02 + sigma);
if( x0 <= 0.0 ) if( x0 <= 0.0 )
@ -347,9 +347,15 @@ namespace gtsam {
v(0) = -sigma / (x0 + mu); v(0) = -sigma / (x0 + mu);
const double v02 = v(0)*v(0); const double v02 = v(0)*v(0);
beta = 2.0 * v02 / (sigma + v02);
v = v / v(0); v = v / v(0);
return 2.0 * v02 / (sigma + v02);
} }
}
/* ************************************************************************* */
pair<double, Vector > house(const Vector &x) {
Vector v(x);
double beta = houseInPlace(v);
return make_pair(beta, v); return make_pair(beta, v);
} }

3
cpp/Vector.h
View File

@ -250,6 +250,9 @@ Vector operator/(double s, const Vector& v);
*/ */
std::pair<double,Vector> house(const Vector &x); std::pair<double,Vector> house(const Vector &x);
/** beta = house(x) computes the HouseHolder vector in place */
double houseInPlace(Vector &x);
/** /**
* Weighted Householder solution vector, * Weighted Householder solution vector,
* a.k.a., the pseudoinverse of the column * a.k.a., the pseudoinverse of the column