Loading...
Searching...
No Matches
AcceleratedPowerMethod.h
Go to the documentation of this file.
/* ----------------------------------------------------------------------------
* GTSAM Copyright 2010-2019, Georgia Tech Research Corporation,
* Atlanta, Georgia 30332-0415
* All Rights Reserved
* Authors: Frank Dellaert, et al. (see THANKS for the full author list)
* See LICENSE for the license information
* -------------------------------------------------------------------------- */
#pragma once
#include <gtsam/linear/PowerMethod.h>
namespace gtsam {
using Sparse = Eigen::SparseMatrix<double>;
template <class Operator>
class AcceleratedPowerMethod : public PowerMethod<Operator> {
double beta_ = 0; // a Polyak momentum term
Vector previousVector_; // store previous vector
public:
explicit AcceleratedPowerMethod(
const Operator &A, const boost::optional<Vector> initial = boost::none,
double initialBeta = 0.0)
: PowerMethod<Operator>(A, initial) {
// initialize Ritz eigen vector and previous vector
this->ritzVector_ = initial ? initial.get() : Vector::Random(this->dim_);
this->ritzVector_.normalize();
previousVector_ = Vector::Zero(this->dim_);
// initialize beta_
beta_ = initialBeta;
}
Vector acceleratedPowerIteration (const Vector &x1, const Vector &x0,
const double beta) const {
Vector y = this->A_ * x1 - beta * x0;
y.normalize();
return y;
}
Vector acceleratedPowerIteration () const {
Vector y = acceleratedPowerIteration(this->ritzVector_, previousVector_, beta_);
return y;
}
double estimateBeta(const size_t T = 10) const {
// set initial estimation of maxBeta
Vector initVector = this->ritzVector_;
const double up = initVector.dot( this->A_ * initVector );
const double down = initVector.dot(initVector);
const double mu = up / down;
double maxBeta = mu * mu / 4;
size_t maxIndex;
std::vector<double> betas;
Matrix R = Matrix::Zero(this->dim_, 10);
// run T times of iteration to find the beta that has the largest Rayleigh quotient
for (size_t t = 0; t < T; t++) {
// after each t iteration, reset the betas with the current maxBeta
betas = {2 / 3 * maxBeta, 0.99 * maxBeta, maxBeta, 1.01 * maxBeta,
1.5 * maxBeta};
// iterate through every beta value
for (size_t k = 0; k < betas.size(); ++k) {
// initialize x0 and x00 in each iteration of each beta
Vector x0 = initVector;
Vector x00 = Vector::Zero(this->dim_);
// run 10 steps of accelerated power iteration with this beta
for (size_t j = 1; j < 10; j++) {
if (j < 2) {
R.col(0) = acceleratedPowerIteration(x0, x00, betas[k]);
R.col(1) = acceleratedPowerIteration(R.col(0), x0, betas[k]);
} else {
R.col(j) = acceleratedPowerIteration(R.col(j - 1), R.col(j - 2),
betas[k]);
}
}
// compute the Rayleigh quotient for the randomly sampled vector after
// 10 steps of power accelerated iteration
const Vector x = R.col(9);
const double up = x.dot(this->A_ * x);
const double down = x.dot(x);
const double mu = up / down;
// store the momentum with largest Rayleigh quotient and its according index of beta_
if (mu * mu / 4 > maxBeta) {
// save the max beta index
maxIndex = k;
maxBeta = mu * mu / 4;
}
}
}
// set beta_ to momentum with largest Rayleigh quotient
return betas[maxIndex];
}
bool compute(size_t maxIterations, double tol) {
// Starting
bool isConverged = false;
for (size_t i = 0; i < maxIterations && !isConverged; i++) {
++(this->nrIterations_);
Vector tmp = this->ritzVector_;
// update the ritzVector after accelerated power iteration
this->ritzVector_ = acceleratedPowerIteration();
// update the previousVector with ritzVector
previousVector_ = tmp;
// update the ritzValue
this->ritzValue_ = this->ritzVector_.dot(this->A_ * this->ritzVector_);
isConverged = this->converged(tol);
}
return isConverged;
}
};
} // namespace gtsam
PowerMethod.h
Power method for fast eigenvalue and eigenvector computation.
gtsam
Global functions in a separate testing namespace.
Definition chartTesting.h:28
gtsam::AcceleratedPowerMethod::acceleratedPowerIteration
Vector acceleratedPowerIteration(const Vector &x1, const Vector &x0, const double beta) const
Run accelerated power iteration to get ritzVector with beta and previous two ritzVector x0 and x00,...
Definition AcceleratedPowerMethod.h:80
gtsam::AcceleratedPowerMethod::AcceleratedPowerMethod
AcceleratedPowerMethod(const Operator &A, const boost::optional< Vector > initial=boost::none, double initialBeta=0.0)
Constructor from aim matrix A (given as Matrix or Sparse), optional intial vector as ritzVector.
Definition AcceleratedPowerMethod.h:62
gtsam::AcceleratedPowerMethod::estimateBeta
double estimateBeta(const size_t T=10) const
Tuning the momentum beta using the Best Heavy Ball algorithm in Ref(3), T is the iteration time to fi...
Definition AcceleratedPowerMethod.h:101
gtsam::AcceleratedPowerMethod::acceleratedPowerIteration
Vector acceleratedPowerIteration() const
Run accelerated power iteration to get ritzVector with beta and previous two ritzVector x0 and x00,...
Definition AcceleratedPowerMethod.h:92
gtsam::AcceleratedPowerMethod::compute
bool compute(size_t maxIterations, double tol)
Start the accelerated iteration, after performing the accelerated iteration, calculate the ritz error...
Definition AcceleratedPowerMethod.h:156
gtsam::PowerMethod::A_
const Operator & A_
Const reference to an externally-held matrix whose minimum-eigenvalue we want to compute.
Definition PowerMethod.h:63
gtsam::PowerMethod::PowerMethod
PowerMethod(const Operator &A, const boost::optional< Vector > initial=boost::none)
Construct from the aim matrix and intial ritz vector.
Definition PowerMethod.h:77
gtsam::PowerMethod::converged
bool converged(double tol) const
After Perform power iteration on a single Ritz value, check if the Ritz residual for the current Ritz...
Definition PowerMethod.h:112