Added more detailed documentation
parent
f86ad95582
commit
32bf81efea
|
|
@ -80,8 +80,7 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
|
||||||
} else {
|
} else {
|
||||||
beta_ = initialBeta;
|
beta_ = initialBeta;
|
||||||
}
|
}
|
||||||
|
}
|
||||||
}
|
|
||||||
|
|
||||||
// Update the ritzVector with beta and previous two ritzVector
|
// Update the ritzVector with beta and previous two ritzVector
|
||||||
Vector powerIteration(const Vector &x1, const Vector &x0,
|
Vector powerIteration(const Vector &x1, const Vector &x0,
|
||||||
|
|
|
||||||
|
|
@ -89,14 +89,16 @@ class PowerMethod {
|
||||||
// Return the number of iterations
|
// Return the number of iterations
|
||||||
size_t nrIterations() const { return nrIterations_; }
|
size_t nrIterations() const { return nrIterations_; }
|
||||||
|
|
||||||
// Start the power/accelerated iteration, after updated the ritz vector,
|
// Start the power/accelerated iteration, after performing the
|
||||||
// calculate the ritz error, repeat this operation until the ritz error converge
|
// power/accelerated power iteration, calculate the ritz error, repeat this
|
||||||
int compute(size_t maxIterations, double tol) {
|
// operation until the ritz error converge. If converged return true, else
|
||||||
|
// false.
|
||||||
|
bool compute(size_t maxIterations, double tol) {
|
||||||
// Starting
|
// Starting
|
||||||
bool isConverged = false;
|
bool isConverged = false;
|
||||||
|
|
||||||
for (size_t i = 0; i < maxIterations; i++) {
|
for (size_t i = 0; i < maxIterations; i++) {
|
||||||
++nrIterations_ ;
|
++nrIterations_;
|
||||||
ritzVector_ = powerIteration();
|
ritzVector_ = powerIteration();
|
||||||
isConverged = converged(tol);
|
isConverged = converged(tol);
|
||||||
if (isConverged) return isConverged;
|
if (isConverged) return isConverged;
|
||||||
|
|
|
||||||
|
|
@ -18,18 +18,16 @@
|
||||||
* @brief Check eigenvalue and eigenvector computed by accelerated power method
|
* @brief Check eigenvalue and eigenvector computed by accelerated power method
|
||||||
*/
|
*/
|
||||||
|
|
||||||
|
#include <CppUnitLite/TestHarness.h>
|
||||||
#include <gtsam/base/Matrix.h>
|
#include <gtsam/base/Matrix.h>
|
||||||
#include <gtsam/base/VectorSpace.h>
|
#include <gtsam/base/VectorSpace.h>
|
||||||
#include <gtsam/inference/Symbol.h>
|
#include <gtsam/inference/Symbol.h>
|
||||||
#include <gtsam/linear/GaussianFactorGraph.h>
|
|
||||||
#include <gtsam/linear/AcceleratedPowerMethod.h>
|
#include <gtsam/linear/AcceleratedPowerMethod.h>
|
||||||
|
#include <gtsam/linear/GaussianFactorGraph.h>
|
||||||
#include <CppUnitLite/TestHarness.h>
|
|
||||||
|
|
||||||
#include <Eigen/Core>
|
#include <Eigen/Core>
|
||||||
#include <Eigen/Dense>
|
#include <Eigen/Dense>
|
||||||
#include <Eigen/Eigenvalues>
|
#include <Eigen/Eigenvalues>
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <random>
|
#include <random>
|
||||||
|
|
||||||
|
|
@ -70,7 +68,7 @@ TEST(AcceleratedPowerMethod, useFactorGraph) {
|
||||||
for (size_t j = 0; j < 3; j++) {
|
for (size_t j = 0; j < 3; j++) {
|
||||||
fg.add(X(j), -I_1x1, X(j + 1), I_1x1, Vector1::Zero(), model);
|
fg.add(X(j), -I_1x1, X(j + 1), I_1x1, Vector1::Zero(), model);
|
||||||
}
|
}
|
||||||
fg.add(X(3), -I_1x1, X(0), I_1x1, Vector1::Zero(), model); // extra row
|
fg.add(X(3), -I_1x1, X(0), I_1x1, Vector1::Zero(), model); // extra row
|
||||||
|
|
||||||
// Get eigenvalues and eigenvectors with Eigen
|
// Get eigenvalues and eigenvectors with Eigen
|
||||||
auto L = fg.hessian();
|
auto L = fg.hessian();
|
||||||
|
|
@ -79,22 +77,26 @@ TEST(AcceleratedPowerMethod, useFactorGraph) {
|
||||||
// Check that we get zero eigenvalue and "constant" eigenvector
|
// Check that we get zero eigenvalue and "constant" eigenvector
|
||||||
EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
|
EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
|
||||||
auto v0 = solver.eigenvectors().col(0);
|
auto v0 = solver.eigenvectors().col(0);
|
||||||
for (size_t j = 0; j < 3; j++)
|
for (size_t j = 0; j < 3; j++) EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);
|
||||||
EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);
|
|
||||||
|
|
||||||
size_t maxIdx = 0;
|
size_t maxIdx = 0;
|
||||||
for (auto i =0; i<solver.eigenvalues().rows(); ++i) {
|
for (auto i = 0; i < solver.eigenvalues().rows(); ++i) {
|
||||||
if (solver.eigenvalues()(i).real() >= solver.eigenvalues()(maxIdx).real()) maxIdx = i;
|
if (solver.eigenvalues()(i).real() >= solver.eigenvalues()(maxIdx).real())
|
||||||
|
maxIdx = i;
|
||||||
}
|
}
|
||||||
// Store the max eigenvalue and its according eigenvector
|
// Store the max eigenvalue and its according eigenvector
|
||||||
const auto ev1 = solver.eigenvalues()(maxIdx).real();
|
const auto ev1 = solver.eigenvalues()(maxIdx).real();
|
||||||
auto ev2 = solver.eigenvectors().col(maxIdx).real();
|
auto ev2 = solver.eigenvectors().col(maxIdx).real();
|
||||||
|
|
||||||
Vector initial = Vector4::Random();
|
Vector disturb = Vector4::Random();
|
||||||
|
disturb.normalize();
|
||||||
|
Vector initial = L.first.row(0);
|
||||||
|
double magnitude = initial.norm();
|
||||||
|
initial += 0.03 * magnitude * disturb;
|
||||||
AcceleratedPowerMethod<Matrix> apf(L.first, initial);
|
AcceleratedPowerMethod<Matrix> apf(L.first, initial);
|
||||||
apf.compute(50, 1e-4);
|
apf.compute(50, 1e-4);
|
||||||
EXPECT_DOUBLES_EQUAL(ev1, apf.eigenvalue(), 1e-8);
|
EXPECT_DOUBLES_EQUAL(ev1, apf.eigenvalue(), 1e-8);
|
||||||
EXPECT(assert_equal(-ev2, apf.eigenvector(), 3e-5));
|
EXPECT(assert_equal(ev2, apf.eigenvector(), 3e-5));
|
||||||
}
|
}
|
||||||
|
|
||||||
/* ************************************************************************* */
|
/* ************************************************************************* */
|
||||||
|
|
|
||||||
|
|
@ -18,18 +18,16 @@
|
||||||
* @brief Check eigenvalue and eigenvector computed by power method
|
* @brief Check eigenvalue and eigenvector computed by power method
|
||||||
*/
|
*/
|
||||||
|
|
||||||
|
#include <CppUnitLite/TestHarness.h>
|
||||||
#include <gtsam/base/Matrix.h>
|
#include <gtsam/base/Matrix.h>
|
||||||
#include <gtsam/base/VectorSpace.h>
|
#include <gtsam/base/VectorSpace.h>
|
||||||
#include <gtsam/inference/Symbol.h>
|
#include <gtsam/inference/Symbol.h>
|
||||||
#include <gtsam/linear/GaussianFactorGraph.h>
|
#include <gtsam/linear/GaussianFactorGraph.h>
|
||||||
#include <gtsam/linear/PowerMethod.h>
|
#include <gtsam/linear/PowerMethod.h>
|
||||||
|
|
||||||
#include <CppUnitLite/TestHarness.h>
|
|
||||||
|
|
||||||
#include <Eigen/Core>
|
#include <Eigen/Core>
|
||||||
#include <Eigen/Dense>
|
#include <Eigen/Dense>
|
||||||
#include <Eigen/Eigenvalues>
|
#include <Eigen/Eigenvalues>
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <random>
|
#include <random>
|
||||||
|
|
||||||
|
|
@ -68,7 +66,7 @@ TEST(PowerMethod, useFactorGraph) {
|
||||||
for (size_t j = 0; j < 3; j++) {
|
for (size_t j = 0; j < 3; j++) {
|
||||||
fg.add(X(j), -I_1x1, X(j + 1), I_1x1, Vector1::Zero(), model);
|
fg.add(X(j), -I_1x1, X(j + 1), I_1x1, Vector1::Zero(), model);
|
||||||
}
|
}
|
||||||
fg.add(X(3), -I_1x1, X(0), I_1x1, Vector1::Zero(), model); // extra row
|
fg.add(X(3), -I_1x1, X(0), I_1x1, Vector1::Zero(), model); // extra row
|
||||||
|
|
||||||
// Get eigenvalues and eigenvectors with Eigen
|
// Get eigenvalues and eigenvectors with Eigen
|
||||||
auto L = fg.hessian();
|
auto L = fg.hessian();
|
||||||
|
|
@ -77,12 +75,12 @@ TEST(PowerMethod, useFactorGraph) {
|
||||||
// Check that we get zero eigenvalue and "constant" eigenvector
|
// Check that we get zero eigenvalue and "constant" eigenvector
|
||||||
EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
|
EXPECT_DOUBLES_EQUAL(0.0, solver.eigenvalues()[0].real(), 1e-9);
|
||||||
auto v0 = solver.eigenvectors().col(0);
|
auto v0 = solver.eigenvectors().col(0);
|
||||||
for (size_t j = 0; j < 3; j++)
|
for (size_t j = 0; j < 3; j++) EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);
|
||||||
EXPECT_DOUBLES_EQUAL(-0.5, v0[j].real(), 1e-9);
|
|
||||||
|
|
||||||
size_t maxIdx = 0;
|
size_t maxIdx = 0;
|
||||||
for (auto i =0; i<solver.eigenvalues().rows(); ++i) {
|
for (auto i = 0; i < solver.eigenvalues().rows(); ++i) {
|
||||||
if (solver.eigenvalues()(i).real() >= solver.eigenvalues()(maxIdx).real()) maxIdx = i;
|
if (solver.eigenvalues()(i).real() >= solver.eigenvalues()(maxIdx).real())
|
||||||
|
maxIdx = i;
|
||||||
}
|
}
|
||||||
// Store the max eigenvalue and its according eigenvector
|
// Store the max eigenvalue and its according eigenvector
|
||||||
const auto ev1 = solver.eigenvalues()(maxIdx).real();
|
const auto ev1 = solver.eigenvalues()(maxIdx).real();
|
||||||
|
|
@ -92,7 +90,7 @@ TEST(PowerMethod, useFactorGraph) {
|
||||||
PowerMethod<Matrix> pf(L.first, initial);
|
PowerMethod<Matrix> pf(L.first, initial);
|
||||||
pf.compute(50, 1e-4);
|
pf.compute(50, 1e-4);
|
||||||
EXPECT_DOUBLES_EQUAL(ev1, pf.eigenvalue(), 1e-8);
|
EXPECT_DOUBLES_EQUAL(ev1, pf.eigenvalue(), 1e-8);
|
||||||
auto actual2 = pf.eigenvector();
|
auto actual2 = pf.eigenvector();
|
||||||
EXPECT(assert_equal(-ev2, actual2, 3e-5));
|
EXPECT(assert_equal(-ev2, actual2, 3e-5));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -482,6 +482,9 @@ Sparse ShonanAveraging<d>::computeA(const Matrix &S) const {
|
||||||
// Perturb the initial initialVector by adding a spherically-uniformly
|
// Perturb the initial initialVector by adding a spherically-uniformly
|
||||||
// distributed random vector with 0.03*||initialVector||_2 magnitude to
|
// distributed random vector with 0.03*||initialVector||_2 magnitude to
|
||||||
// initialVector
|
// initialVector
|
||||||
|
// ref : Part III. C, Rosen, D. and Carlone, L., 2017, September. Computational
|
||||||
|
// enhancements for certifiably correct SLAM. In Proceedings of the
|
||||||
|
// International Conference on Intelligent Robots and Systems.
|
||||||
Vector perturb(const Vector &initialVector) {
|
Vector perturb(const Vector &initialVector) {
|
||||||
// generate a 0.03*||x_0||_2 as stated in David's paper
|
// generate a 0.03*||x_0||_2 as stated in David's paper
|
||||||
int n = initialVector.rows();
|
int n = initialVector.rows();
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue