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>

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 boost::optional<Vector> initial = boost::none)
      : A_(A), dim_(A.rows()), nrIterations_(0) {
    Vector x0;
    x0 = initial ? initial.get() : 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;

    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);
    }

    return isConverged;
  }

  double eigenvalue() const { return ritzValue_; }

  Vector eigenvector() const { return ritzVector_; }
};

}  // namespace gtsam