PowerMethod.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/base/Matrix.h>
#include <gtsam/base/Vector.h>
 
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <random>
#include <vector>
#include <optional>
 
namespace gtsam {
 
using Sparse = Eigen::SparseMatrix<double>;

template <class Operator>
class PowerMethod {
 protected:
  const Operator &A_;
 
  const int dim_;  // dimension of Matrix A
 
  size_t nrIterations_;  // number of iterations
 
  double ritzValue_;   // Ritz eigenvalue
  Vector ritzVector_;  // Ritz eigenvector
 
 public:

  explicit PowerMethod(const Operator &A,
                       const std::optional<Vector> initial = {})
      : A_(A), dim_(A.rows()), nrIterations_(0) {
    Vector x0;
    x0 = initial ? *initial : Vector::Random(dim_);
    x0.normalize();
 
    // initialize Ritz eigen value
    ritzValue_ = 0.0;
 
    // initialize Ritz eigen vector
    ritzVector_ = powerIteration(x0);
  }

  Vector powerIteration(const Vector &x) const {
    Vector y = A_ * x;
    y.normalize();
    return y;
  }
  
  Vector powerIteration() const { return powerIteration(ritzVector_); }

  bool converged(double tol) const {
    const Vector x = ritzVector_;
    // store the Ritz eigen value
    const double ritzValue = x.dot(A_ * x);
    const double error = (A_ * x - ritzValue * x).norm();
    return error < tol;
  }

  size_t nrIterations() const { return nrIterations_; }

  bool compute(size_t maxIterations, double tol) {
    // Starting
    bool isConverged = false;
    if (maxIterations == 0) return isConverged;
 
    Vector product = A_ * ritzVector_;
 
    for (size_t i = 0; i < maxIterations && !isConverged; i++) {
      ++nrIterations_;
      // update the ritzVector after power iteration
      ritzVector_ = product;
      ritzVector_.normalize();
 
      // Reuse this product for the Ritz value, residual, and next iteration.
      product = A_ * ritzVector_;
      ritzValue_ = ritzVector_.dot(product);
      const double error = (product - ritzValue_ * ritzVector_).norm();
      isConverged = error < tol;
    }
 
    return isConverged;
  }

  double eigenvalue() const { return ritzValue_; }

  Vector eigenvector() const { return ritzVector_; }
};
 
}  // namespace gtsam