Revised as David's second review

gtsam/linear/AcceleratedPowerMethod.h

@@ -61,7 +61,7 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
   // vector as ritzVector
   explicit AcceleratedPowerMethod(
       const Operator &A, const boost::optional<Vector> initial = boost::none,
-      const double initialBeta = 0.0)
+      double initialBeta = 0.0)
       : PowerMethod<Operator>(A, initial) {
     // initialize Ritz eigen vector and previous vector
     this->ritzVector_ = initial ? initial.get() : Vector::Random(this->dim_);
@@ -76,61 +76,97 @@ class AcceleratedPowerMethod : public PowerMethod<Operator> {
     }
   }
 
-  // Update the ritzVector with beta and previous two ritzVector, and return
+  // Run accelerated power iteration on the ritzVector with beta and previous
-  // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
+  // two ritzVector, and return x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k
-  Vector powerIteration(const Vector &x1, const Vector &x0,
+  // - \beta * x_{k-1} ||
+  Vector acceleratedPowerIteration (const Vector &x1, const Vector &x0,
                         const double beta) const {
     Vector y = this->A_ * x1 - beta * x0;
     y.normalize();
     return y;
   }
 
-  // Update the ritzVector with beta and previous two ritzVector, and return
+  // Run accelerated power iteration on the ritzVector with beta and previous
-  // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} ||
+  // two ritzVector, and return x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k
-  Vector powerIteration() const {
+  // - \beta * x_{k-1} ||
-    Vector y = powerIteration(this->ritzVector_, previousVector_, beta_);
+  Vector acceleratedPowerIteration () const {
-    previousVector_ = this->ritzVector_;
+    Vector y = acceleratedPowerIteration(this->ritzVector_, previousVector_, beta_);
     return y;
   }
 
   // Tuning the momentum beta using the Best Heavy Ball algorithm in Ref(3)
   void estimateBeta() {
     // set beta
-    Vector init_resid = this->ritzVector_;
+    Vector initVector = this->ritzVector_;
-    const double up = init_resid.dot( this->A_ * init_resid );
+    const double up = initVector.dot( this->A_ * initVector );
-    const double down = init_resid.dot(init_resid);
+    const double down = initVector.dot(initVector);
     const double mu = up / down;
     double maxBeta = mu * mu / 4;
     size_t maxIndex;
-    std::vector<double> betas = {2 / 3 * maxBeta, 0.99 * maxBeta, maxBeta,
+    std::vector<double> betas;
-                                 1.01 * maxBeta, 1.5 * maxBeta};
 
     Matrix R = Matrix::Zero(this->dim_, 10);
-    for (size_t i = 0; i < 10; i++) {
+    size_t T = 10;
-      Vector x0 = Vector::Random(this->dim_);
+    // run T times of iteration to find the beta that has the largest Rayleigh quotient
-      x0.normalize();
+    for (size_t t = 0; t < T; t++) {
-      Vector x00 = Vector::Zero(this->dim_);
+      // 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) = powerIteration(x0, x00, betas[k]);
+            R.col(0) = acceleratedPowerIteration(x0, x00, betas[k]);
-            R.col(1) = powerIteration(R.col(0), x0, betas[k]);
+            R.col(1) = acceleratedPowerIteration(R.col(0), x0, betas[k]);
           } else {
-            R.col(j) = powerIteration(R.col(j - 1), R.col(j - 2), betas[k]);
+            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
     beta_ = betas[maxIndex];
   }
+
+  // Start the accelerated iteration, after performing the
+  // accelerated iteration, calculate the ritz error, repeat this
+  // operation until the ritz error converge. If converged return true, else
+  // false.
+  bool compute(size_t maxIterations, double tol) {
+    // Starting
+    bool isConverged = false;
+
+    for (size_t i = 0; i < maxIterations; 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);
+      if (isConverged) return isConverged;
+    }
+
+    return isConverged;
+  }
 };
 
 }  // namespace gtsam

gtsam/linear/PowerMethod.h

@@ -66,23 +66,24 @@ class PowerMethod {
     ritzVector_ = powerIteration(x0);
   }
 
-  // Update the vector by dot product with A_, and return A * x / || A * x ||
+  // Run power iteration on the vector, 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_, and return A * x / || A * x ||
+  // Run power iteration on the vector, 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
+  // After Perform power iteration on a single Ritz value, check if the Ritz
-  // than the tol then return true else false
+  // residual for the current Ritz pair is less than the required convergence
-  bool converged(double tol) {
+  // tol, return true if yes, else false
+  bool converged(double tol) const {
     const Vector x = ritzVector_;
     // store the Ritz eigen value
-    ritzValue_ = x.dot(A_ * x);
+    const double ritzValue = x.dot(A_ * x);
-    double error = (A_ * x - ritzValue_ * x).norm();
+    const double error = (A_ * x - ritzValue * x).norm();
     return error < tol;
   }
 
@@ -90,18 +91,20 @@ class PowerMethod {
   size_t nrIterations() const { return nrIterations_; }
 
   // Start the power/accelerated iteration, after performing the
-  // power/accelerated power iteration, calculate the ritz error, repeat this
+  // power/accelerated iteration, calculate the ritz error, repeat this
   // operation until the ritz error converge. If converged return true, else
   // false.
   bool compute(size_t maxIterations, double tol) {
     // Starting
     bool isConverged = false;
 
-    for (size_t i = 0; i < maxIterations; i++) {
+    for (size_t i = 0; i < maxIterations && !isConverged; i++) {
       ++nrIterations_;
+      // update the ritzVector after power iteration
       ritzVector_ = powerIteration();
+      // update the ritzValue 
+      ritzValue_ = ritzVector_.dot(A_ * ritzVector_);
       isConverged = converged(tol);
-      if (isConverged) return isConverged;
     }
 
     return isConverged;

gtsam/sfm/ShonanAveraging.cpp

@@ -509,7 +509,7 @@ Vector perturb(const Vector &initialVector) {
 // compute the maximum eigenpair (\theta, \vect{v}) of C = \lamda_dom * I - A by
 // accelerated power method. Then return (\lamda_dom - \theta, \vect{v}).
 static bool PowerMinimumEigenValue(
-    const Sparse &A, const Matrix &S, double *minEigenValue,
+    const Sparse &A, const Matrix &S, double &minEigenValue,
     Vector *minEigenVector = 0, size_t *numIterations = 0,
     size_t maxIterations = 1000,
     double minEigenvalueNonnegativityTolerance = 10e-4) {
@@ -528,7 +528,7 @@ static bool PowerMinimumEigenValue(
   if (pmEigenValue < 0) {
     // The largest-magnitude eigenvalue is negative, and therefore also the
     // minimum eigenvalue, so just return this solution
-    *minEigenValue = pmEigenValue;
+    minEigenValue = pmEigenValue;
     if (minEigenVector) {
       *minEigenVector = pmOperator.eigenvector();
       minEigenVector->normalize();  // Ensure that this is a unit vector
@@ -545,7 +545,7 @@ static bool PowerMinimumEigenValue(
 
   if (!minConverged) return false;
 
-  *minEigenValue = pmEigenValue - apmShiftedOperator.eigenvalue();
+  minEigenValue = pmEigenValue - apmShiftedOperator.eigenvalue();
   if (minEigenVector) {
     *minEigenVector = apmShiftedOperator.eigenvector();
     minEigenVector->normalize();  // Ensure that this is a unit vector
@@ -723,8 +723,8 @@ double ShonanAveraging<d>::computeMinEigenValueAP(const Values &values,
   const Matrix S = StiefelElementMatrix(values);
   auto A = computeA(S);
 
-  double minEigenValue;
+  double minEigenValue = 0;
-  bool success = PowerMinimumEigenValue(A, S, &minEigenValue, minEigenVector);
+  bool success = PowerMinimumEigenValue(A, S, minEigenValue, minEigenVector);
   if (!success) {
     throw std::runtime_error(
         "PowerMinimumEigenValue failed to compute minimum eigenvalue.");