From cf813c5a64b213ffd46a7ae7a6add1a2f342b80c Mon Sep 17 00:00:00 2001 From: jingwuOUO Date: Mon, 26 Oct 2020 15:40:49 -0400 Subject: [PATCH] Revised as David's second review --- gtsam/linear/AcceleratedPowerMethod.h | 78 +++++++++++++++++++-------- gtsam/linear/PowerMethod.h | 23 ++++---- gtsam/sfm/ShonanAveraging.cpp | 10 ++-- 3 files changed, 75 insertions(+), 36 deletions(-) diff --git a/gtsam/linear/AcceleratedPowerMethod.h b/gtsam/linear/AcceleratedPowerMethod.h index a5148e1d4..eea007353 100644 --- a/gtsam/linear/AcceleratedPowerMethod.h +++ b/gtsam/linear/AcceleratedPowerMethod.h @@ -61,7 +61,7 @@ class AcceleratedPowerMethod : public PowerMethod { // vector as ritzVector explicit AcceleratedPowerMethod( const Operator &A, const boost::optional initial = boost::none, - const double initialBeta = 0.0) + double initialBeta = 0.0) : PowerMethod(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 { } } - // Update the ritzVector with beta and previous two ritzVector, and return - // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} || - Vector powerIteration(const Vector &x1, const Vector &x0, + // Run accelerated power iteration on the ritzVector with beta and previous + // two ritzVector, and return x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k + // - \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 - // x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k - \beta * x_{k-1} || - Vector powerIteration() const { - Vector y = powerIteration(this->ritzVector_, previousVector_, beta_); - previousVector_ = this->ritzVector_; + // Run accelerated power iteration on the ritzVector with beta and previous + // two ritzVector, and return x_{k+1} = A * x_k - \beta * x_{k-1} / || A * x_k + // - \beta * x_{k-1} || + Vector acceleratedPowerIteration () const { + 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_; - const double up = init_resid.dot( this->A_ * init_resid ); - const double down = init_resid.dot(init_resid); + 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 betas = {2 / 3 * maxBeta, 0.99 * maxBeta, maxBeta, - 1.01 * maxBeta, 1.5 * maxBeta}; + std::vector betas; Matrix R = Matrix::Zero(this->dim_, 10); - for (size_t i = 0; i < 10; i++) { - Vector x0 = Vector::Random(this->dim_); - x0.normalize(); - Vector x00 = Vector::Zero(this->dim_); + size_t T = 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) = powerIteration(x0, x00, betas[k]); - R.col(1) = powerIteration(R.col(0), x0, betas[k]); + R.col(0) = acceleratedPowerIteration(x0, x00, 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 diff --git a/gtsam/linear/PowerMethod.h b/gtsam/linear/PowerMethod.h index 8a577aa92..ae1d97cc7 100644 --- a/gtsam/linear/PowerMethod.h +++ b/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 - // than the tol then return true else false - bool converged(double tol) { + // After Perform power iteration on a single Ritz value, check if the Ritz + // residual for the current Ritz pair is less than the required convergence + // 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); - double error = (A_ * x - ritzValue_ * x).norm(); + const double ritzValue = x.dot(A_ * x); + 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; diff --git a/gtsam/sfm/ShonanAveraging.cpp b/gtsam/sfm/ShonanAveraging.cpp index 8d2bf4653..fab57c828 100644 --- a/gtsam/sfm/ShonanAveraging.cpp +++ b/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::computeMinEigenValueAP(const Values &values, const Matrix S = StiefelElementMatrix(values); auto A = computeA(S); - double minEigenValue; - bool success = PowerMinimumEigenValue(A, S, &minEigenValue, minEigenVector); + double minEigenValue = 0; + bool success = PowerMinimumEigenValue(A, S, minEigenValue, minEigenVector); if (!success) { throw std::runtime_error( "PowerMinimumEigenValue failed to compute minimum eigenvalue.");