AcceleratedPowerMethod.h
Go to the documentation of this file.
1 /* ----------------------------------------------------------------------------
2 
3  * GTSAM Copyright 2010-2019, Georgia Tech Research Corporation,
4  * Atlanta, Georgia 30332-0415
5  * All Rights Reserved
6  * Authors: Frank Dellaert, et al. (see THANKS for the full author list)
7 
8  * See LICENSE for the license information
9 
10  * -------------------------------------------------------------------------- */
11 
20 #pragma once
21 
23 
24 namespace gtsam {
25 
27 
50 template <class Operator>
51 class AcceleratedPowerMethod : public PowerMethod<Operator> {
52 
53  double beta_ = 0; // a Polyak momentum term
54 
55  Vector previousVector_; // store previous vector
56 
57  public:
63  const Operator &A, const std::optional<Vector> initial = {},
64  double initialBeta = 0.0)
66  // initialize Ritz eigen vector and previous vector
67  this->ritzVector_ = initial ? *initial : Vector::Random(this->dim_);
68  this->ritzVector_.normalize();
69  previousVector_ = Vector::Zero(this->dim_);
70 
71  // initialize beta_
72  beta_ = initialBeta;
73  }
74 
81  const double beta) const {
82  Vector y = this->A_ * x1 - beta * x0;
83  y.normalize();
84  return y;
85  }
86 
93  Vector y = acceleratedPowerIteration(this->ritzVector_, previousVector_, beta_);
94  return y;
95  }
96 
101  double estimateBeta(const size_t T = 10) const {
102  // set initial estimation of maxBeta
103  Vector initVector = this->ritzVector_;
104  const double up = initVector.dot( this->A_ * initVector );
105  const double down = initVector.dot(initVector);
106  const double mu = up / down;
107  double maxBeta = mu * mu / 4;
108  size_t maxIndex;
109  std::vector<double> betas;
110 
111  Matrix R = Matrix::Zero(this->dim_, 10);
112  // run T times of iteration to find the beta that has the largest Rayleigh quotient
113  for (size_t t = 0; t < T; t++) {
114  // after each t iteration, reset the betas with the current maxBeta
115  betas = {2 / 3 * maxBeta, 0.99 * maxBeta, maxBeta, 1.01 * maxBeta,
116  1.5 * maxBeta};
117  // iterate through every beta value
118  for (size_t k = 0; k < betas.size(); ++k) {
119  // initialize x0 and x00 in each iteration of each beta
120  Vector x0 = initVector;
121  Vector x00 = Vector::Zero(this->dim_);
122  // run 10 steps of accelerated power iteration with this beta
123  for (size_t j = 1; j < 10; j++) {
124  if (j < 2) {
125  R.col(0) = acceleratedPowerIteration(x0, x00, betas[k]);
126  R.col(1) = acceleratedPowerIteration(R.col(0), x0, betas[k]);
127  } else {
128  R.col(j) = acceleratedPowerIteration(R.col(j - 1), R.col(j - 2),
129  betas[k]);
130  }
131  }
132  // compute the Rayleigh quotient for the randomly sampled vector after
133  // 10 steps of power accelerated iteration
134  const Vector x = R.col(9);
135  const double up = x.dot(this->A_ * x);
136  const double down = x.dot(x);
137  const double mu = up / down;
138  // store the momentum with largest Rayleigh quotient and its according index of beta_
139  if (mu * mu / 4 > maxBeta) {
140  // save the max beta index
141  maxIndex = k;
142  maxBeta = mu * mu / 4;
143  }
144  }
145  }
146  // set beta_ to momentum with largest Rayleigh quotient
147  return betas[maxIndex];
148  }
149 
156  bool compute(size_t maxIterations, double tol) {
157  // Starting
158  bool isConverged = false;
159 
160  for (size_t i = 0; i < maxIterations && !isConverged; i++) {
161  ++(this->nrIterations_);
162  Vector tmp = this->ritzVector_;
163  // update the ritzVector after accelerated power iteration
165  // update the previousVector with ritzVector
166  previousVector_ = tmp;
167  // update the ritzValue
168  this->ritzValue_ = this->ritzVector_.dot(this->A_ * this->ritzVector_);
169  isConverged = this->converged(tol);
170  }
171 
172  return isConverged;
173  }
174 };
175 
176 } // namespace gtsam
Vector acceleratedPowerIteration(const Vector &x1, const Vector &x0, const double beta) const
const Operator & A_
Definition: PowerMethod.h:64
AcceleratedPowerMethod(const Operator &A, const std::optional< Vector > initial={}, double initialBeta=0.0)
Scalar * y
Compute maximum Eigenpair with accelerated power method.
double mu
Eigen::MatrixXd Matrix
Definition: base/Matrix.h:39
Rot2 R(Rot2::fromAngle(0.1))
Values initial
Eigen::VectorXd Vector
Definition: Vector.h:38
Power method for fast eigenvalue and eigenvector computation.
Compute maximum Eigenpair with power method.
Definition: PowerMethod.h:58
double estimateBeta(const size_t T=10) const
Eigen::Triplet< double > T
bool compute(size_t maxIterations, double tol)
traits
Definition: chartTesting.h:28
static Symbol x0('x', 0)
Pose3 x1
Definition: testPose3.cpp:663
const G double tol
Definition: Group.h:86
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
bool converged(double tol) const
Definition: PowerMethod.h:113
std::ptrdiff_t j
Point2 t(10, 10)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:33:53