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

Fixed not solve errors and detailed documentation

release/4.3a0
jingwuOUO 2020-10-22 18:05:09 -04:00
parent 32bf81efea
commit 56300ca23c
4 changed files with 33 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
gtsam/linear/AcceleratedPowerMethod.h
View File

@ -63,16 +63,10 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
const Operator &A, const boost::optional<Vector> initial = boost::none,
const double initialBeta = 0.0)
: PowerMethod<Operator>(A, initial) {
Vector x0;
// initialize ritz vector
x0 = initial ? initial.get() : Vector::Random(this->dim_);
Vector x00 = Vector::Random(this->dim_);
x0.normalize();
x00.normalize();
// initialize Ritz eigen vector and previous vector
previousVector_ = powerIteration(x0, x00, beta_);
this->ritzVector_ = powerIteration(previousVector_, x0, beta_);
this->ritzVector_ = initial ? initial.get() : Vector::Random(this->dim_);
this->ritzVector_.normalize();
previousVector_ = Vector::Zero(this->dim_);
// initialize beta_
if (!initialBeta) {
@ -80,9 +74,10 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
} else {
beta_ = initialBeta;
}
}
}
// Update the ritzVector with beta and previous two ritzVector
// Update the ritzVector with beta and previous two ritzVector, and return
// x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
Vector powerIteration(const Vector &x1, const Vector &x0,
const double beta) const {
Vector y = this->A_ * x1 - beta * x0;
@ -90,7 +85,8 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
return y;
}
// Update the ritzVector with beta and previous two ritzVector
// Update the ritzVector with beta and previous two ritzVector, and return
// x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
Vector powerIteration() const {
Vector y = powerIteration(this->ritzVector_, previousVector_, beta_);
previousVector_ = this->ritzVector_;
@ -101,8 +97,8 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
void estimateBeta() {
// set beta
Vector init_resid = this->ritzVector_;
const double up = init_resid.transpose() * this->A_ * init_resid;
const double down = init_resid.transpose().dot(init_resid);
const double up = init_resid.dot( this->A_ * init_resid );
const double down = init_resid.dot(init_resid);
const double mu = up / down;
double maxBeta = mu * mu / 4;
size_t maxIndex;
@ -111,11 +107,12 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
Matrix R = Matrix::Zero(this->dim_, 10);
for (size_t i = 0; i < 10; i++) {
Vector x0 = Vector::Random(this->dim_);
x0.normalize();
Vector x00 = Vector::Zero(this->dim_);
for (size_t k = 0; k < betas.size(); ++k) {
for (size_t j = 1; j < 10; j++) {
if (j < 2) {
Vector x0 = this->ritzVector_;
Vector x00 = previousVector_;
R.col(0) = powerIteration(x0, x00, betas[k]);
R.col(1) = powerIteration(R.col(0), x0, betas[k]);
} else {
@ -123,8 +120,8 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
}
}
const Vector x = R.col(9);
const double up = x.transpose() * this->A_ * x;
const double down = x.transpose().dot(x);
const double up = x.dot(this->A_ * x);
const double down = x.dot(x);
const double mu = up / down;
if (mu * mu / 4 > maxBeta) {
maxIndex = k;

6
gtsam/linear/PowerMethod.h
View File

@ -61,19 +61,19 @@ class PowerMethod {
// initialize Ritz eigen value
ritzValue_ = 0.0;
// initialize Ritz eigen vectors
// initialize Ritz eigen vector
ritzVector_ = Vector::Zero(dim_);
ritzVector_ = powerIteration(x0);
}
// Update the vector by dot product with A_
// Update the vector by dot product with A_, and return A * x / || A * x ||
Vector powerIteration(const Vector &x) const {
Vector y = A_ * x;
y.normalize();
return y;
}
// Update the vector by dot product with A_
// Update the vector by dot product with A_, and return A * x / || A * x ||
Vector powerIteration() const { return powerIteration(ritzVector_); }
// After Perform power iteration on a single Ritz value, if the error is less

2
gtsam/linear/tests/testAcceleratedPowerMethod.cpp
View File

@ -96,7 +96,7 @@ TEST(AcceleratedPowerMethod, useFactorGraph) {
AcceleratedPowerMethod<Matrix> apf(L.first, initial);
apf.compute(50, 1e-4);
EXPECT_DOUBLES_EQUAL(ev1, apf.eigenvalue(), 1e-8);
EXPECT(assert_equal(ev2, apf.eigenvector(), 3e-5));
EXPECT(assert_equal(ev2, apf.eigenvector(), 4e-5));
}
/* ************************************************************************* */

28
gtsam/sfm/ShonanAveraging.cpp
View File

@ -515,42 +515,42 @@ static bool PowerMinimumEigenValue(
double minEigenvalueNonnegativityTolerance = 10e-4) {
// a. Compute dominant eigenpair of S using power method
PowerMethod<Sparse> lmOperator(A);
PowerMethod<Sparse> pmOperator(A);
const bool lmConverged = lmOperator.compute(
const bool pmConverged = pmOperator.compute(
maxIterations, 1e-5);
// Check convergence and bail out if necessary
if (!lmConverged) return false;
if (!pmConverged) return false;
const double lmEigenValue = lmOperator.eigenvalue();
const double pmEigenValue = pmOperator.eigenvalue();
if (lmEigenValue < 0) {
if (pmEigenValue < 0) {
// The largest-magnitude eigenvalue is negative, and therefore also the
// minimum eigenvalue, so just return this solution
*minEigenValue = lmEigenValue;
*minEigenValue = pmEigenValue;
if (minEigenVector) {
*minEigenVector = lmOperator.eigenvector();
*minEigenVector = pmOperator.eigenvector();
minEigenVector->normalize(); // Ensure that this is a unit vector
}
return true;
}
const Sparse C = lmEigenValue * Matrix::Identity(A.rows(), A.cols()) - A;
const Sparse C = pmEigenValue * Matrix::Identity(A.rows(), A.cols()) - A;
const boost::optional<Vector> initial = perturb(S.row(0));
AcceleratedPowerMethod<Sparse> minShiftedOperator(C, initial);
AcceleratedPowerMethod<Sparse> apmShiftedOperator(C, initial);
const bool minConverged = minShiftedOperator.compute(
maxIterations, minEigenvalueNonnegativityTolerance / lmEigenValue);
const bool minConverged = apmShiftedOperator.compute(
maxIterations, minEigenvalueNonnegativityTolerance / pmEigenValue);
if (!minConverged) return false;
*minEigenValue = lmEigenValue - minShiftedOperator.eigenvalue();
*minEigenValue = pmEigenValue - apmShiftedOperator.eigenvalue();
if (minEigenVector) {
*minEigenVector = minShiftedOperator.eigenvector();
*minEigenVector = apmShiftedOperator.eigenvector();
minEigenVector->normalize(); // Ensure that this is a unit vector
}
if (numIterations) *numIterations = minShiftedOperator.nrIterations();
if (numIterations) *numIterations = apmShiftedOperator.nrIterations();
return true;
}