From 080437968bb2dc115c2f4d55ddf92adac0417096 Mon Sep 17 00:00:00 2001
From: Ricardo Ortiz <ricardo.ortiz@kitware.com>
Date: Fri, 15 Jan 2016 10:08:26 -0500
Subject: [PATCH] ENH: Rename the nonlinear solver
to just NewtonMethod. Update the solve routine
so it always goes into the armijo loop, subsequent
function norms are computed there. If the user wants
to disable the armijo updated then it should set the
armijoMax in the base class. The armijo function now
takes an extra parameter, the old function evaluation
norm.
---
Solvers/NewtonMethod.cpp | 155 +++++++++++++++++++++
Solvers/NewtonMethod.h | 260 ++++++++++++++++++++++++++++++++++++
Solvers/NonLinearSolver.cpp | 16 ++-
Solvers/NonLinearSolver.h | 2 +-
4 files changed, 429 insertions(+), 4 deletions(-)
create mode 100644 Solvers/NewtonMethod.cpp
create mode 100644 Solvers/NewtonMethod.h
diff --git a/Solvers/NewtonMethod.cpp b/Solvers/NewtonMethod.cpp
new file mode 100644
index 000000000..ff659b50a
--- /dev/null
+++ b/Solvers/NewtonMethod.cpp
@@ -0,0 +1,155 @@
+// This file is part of the SimMedTK project.
+// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
+// Rensselaer Polytechnic Institute
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+//
+//---------------------------------------------------------------------------
+//
+// Authors:
+//
+// Contact:
+//---------------------------------------------------------------------------
+
+#include "Solvers/NewtonMethod.h"
+
+#include "Solvers/IterativeLinearSolver.h"
+#include "Solvers/ConjugateGradient.h"
+#include <iostream>
+
+NewtonMethod::NewtonMethod():
+ linearSolver(std::make_shared<ConjugateGradient>()),
+ forcingTerm(0.9),
+ absoluteTolerance(1e-3),
+ relativeTolerance(1e-6),
+ gamma(0.9),
+ etaMax(0.9),
+ maxIterations(50),
+ useArmijo(true) {}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::solve(core::Vectord &x)
+{
+ if(!this->linearSolver)
+ {
+ return;
+ }
+
+ // Compute norms, set tolerances and other temporaries
+ double fnorm = this->nonLinearSystem->eval(x).norm();
+ double stopTolerance = this->absoluteTolerance + this->relativeTolerance * fnorm;
+ this->linearSolver->setTolerance(stopTolerance);
+ core::Vectord dx = x;
+
+ for(size_t i = 0; i < this->maxIterations; ++i)
+ {
+ if(fnorm < stopTolerance)
+ {
+ return;
+ }
+ this->updateJacobian(x);
+ this->linearSolver->solve(dx);
+ this->updateIterate(-dx,x);
+
+ double newNorm = fnorm;
+
+ newNorm = this->armijo(dx, x, fnorm);
+
+ if(this->forcingTerm > 0.0 && newNorm > stopTolerance)
+ {
+ double ratio = newNorm / fnorm; // Ratio of successive residual norms
+ this->updateForcingTerm(ratio, stopTolerance, fnorm);
+
+ // Reset tolerance in the linear solver according to the new forcing term
+ // to avoid over solving of the system.
+ this->linearSolver->setTolerance(this->forcingTerm);
+ }
+
+ fnorm = newNorm;
+ }
+}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::updateForcingTerm(const double ratio,
+ const double stopTolerance,
+ const double fnorm)
+{
+ double eta = this->gamma * ratio * ratio;
+ double forcingTermSqr = this->forcingTerm * this->forcingTerm;
+
+ // Save guard to prevent the forcing term to become too small for far away iterates
+ if(this->gamma * forcingTermSqr > .1)
+ {
+ eta = std::max(eta, this->gamma * forcingTermSqr);
+ }
+
+ this->forcingTerm = std::max(std::min(eta, this->etaMax), 0.5 * stopTolerance / fnorm);
+}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::setLinearSolver(std::shared_ptr< NewtonMethod::LinearSolverType > newLinearSolver)
+{
+ this->linearSolver = newLinearSolver;
+}
+
+//---------------------------------------------------------------------------
+std::shared_ptr< NewtonMethod::LinearSolverType > NewtonMethod::getLinearSolver() const
+{
+ return this->linearSolver;
+}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::updateJacobian(const core::Vectord &x)
+{
+ // Evaluate the Jacobian and sets the matrix
+ if(!this->jacobian)
+ {// TODO: Print message or log error.
+ return;
+ }
+
+ const core::SparseMatrixd &jacobianMatrix = this->jacobian(x);
+ const core::Vectord &f = this->nonLinearSystem->eval(x);
+
+ if(!(jacobianMatrix.innerSize() > 0))
+ {// TODO: Print message and/or log error.
+ return;
+ }
+
+ auto linearSystem = std::make_shared<LinearSolverType::LinearSystemType>(
+ jacobianMatrix, f);
+ this->linearSolver->setSystem(linearSystem);
+}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::setJacobian(const NonLinearSolver::JacobianType &newJacobian)
+{
+ this->jacobian = newJacobian;
+}
+
+//---------------------------------------------------------------------------
+const NonLinearSolver::JacobianType &NewtonMethod::getJacobian() const
+{
+ return this->jacobian;
+}
+
+//---------------------------------------------------------------------------
+void NewtonMethod::setAbsoluteTolerance(const double aTolerance)
+{
+ this->absoluteTolerance = aTolerance;
+}
+
+//---------------------------------------------------------------------------
+double NewtonMethod::getAbsoluteTolerance() const
+{
+ return this->absoluteTolerance;
+}
diff --git a/Solvers/NewtonMethod.h b/Solvers/NewtonMethod.h
new file mode 100644
index 000000000..8093a8782
--- /dev/null
+++ b/Solvers/NewtonMethod.h
@@ -0,0 +1,260 @@
+// This file is part of the SimMedTK project.
+// Copyright (c) Center for Modeling, Simulation, and Imaging in Medicine,
+// Rensselaer Polytechnic Institute
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+//
+//---------------------------------------------------------------------------
+//
+// Authors:
+//
+// Contact:
+//---------------------------------------------------------------------------
+
+#ifndef INEXACT_NEWTON_H
+#define INEXACT_NEWTON_H
+
+#include "Solvers/NonLinearSolver.h"
+#include "Solvers/LinearSolver.h"
+
+///
+/// @brief Newton method. This version of the newton method is based on the work
+/// by Tim Kelly and others at NC State University. The algorithm is globally convergent
+/// in the sense that for any initial iterate the iteration either converge to a
+/// root of F or fails. Global convergence is achieved using a line search sub-process
+/// and the Armijo rule.
+///
+class NewtonMethod : public NonLinearSolver
+{
+public:
+ using LinearSolverType = LinearSolver<core::SparseMatrixd>;
+
+public:
+ ///
+ /// \brief Constructors/Destructors.
+ ///
+ NewtonMethod();
+ ~NewtonMethod() = default;
+ NewtonMethod(const NewtonMethod &other) = delete;
+
+ NewtonMethod &operator=(const NewtonMethod &other) = delete;
+
+ ///
+ /// \brief Solve the non linear system of equations G(x)=0 using Newton's method.
+ ///
+ /// \param x Current iterate.
+ ///
+ void solve(core::Vectord &x) override;
+
+ ///
+ /// \brief Update forcing term according to Eisenstat-Walker criteria.
+ /// \see http://softlib.rice.edu/pub/CRPC-TRs/reports/CRPC-TR94463.pdf.
+ ///
+ /// \param ratio Ratio of consecutive function norms.
+ /// \param stopTolerance Linear combinarion of relative and absolute tolerances.
+ /// \param fnorm Current function norm.
+ ///
+ void updateForcingTerm(const double ratio,
+ const double stopTolerance,
+ const double fnorm);
+
+ ///
+ /// \brief Set LinearSolver.
+ ///
+ /// \param newLinearSolver Linear solver pointer.
+ ///
+ void setLinearSolver(std::shared_ptr<LinearSolverType> newLinearSolver);
+
+ ///
+ /// \brief Get LinearSolver.
+ ///
+ std::shared_ptr<LinearSolverType> getLinearSolver() const;
+
+ ///
+ /// \brief Update jacobians.
+ ///
+ /// \param x Current iterate.
+ ///
+ void updateJacobian(const core::Vectord &x);
+
+ ///
+ /// \brief Set Jacobian. This sets the jacobian function, used to compute the jacobian matrix.
+ ///
+ /// \param newJacobian Jacobian function.
+ ///
+ void setJacobian(const JacobianType &newJacobian);
+
+ ///
+ /// \brief Get Jacobian. Returns the Jacobian function.
+ ///
+ const JacobianType &getJacobian() const;
+
+ ///
+ /// \brief Set JacobianMatrix.
+ ///
+ /// \param newJacobianMatrix Replaces the stored jacobian matrix.
+ ///
+ void setJacobianMatrix(const core::SparseMatrixd &newJacobianMatrix);
+
+ ///
+ /// \brief Get JacobianMatrix. Returns jacobian matrix.
+ ///
+ const core::SparseMatrixd &getJacobianMatrix() const;
+
+ ///
+ /// \brief Get JacobianMatrix. Returns jacobian matrix.
+ ///
+ core::SparseMatrixd &getJacobianMatrix();
+
+ ///
+ /// \brief Set AbsoluteTolerance.
+ ///
+ /// \param aTolerance New absolute tolerance value.
+ ///
+ void setAbsoluteTolerance(const double aTolerance);
+
+ ///
+ /// \brief Get AbsoluteTolerance. Returns current tolerance value.
+ ///
+ double getAbsoluteTolerance() const;
+
+ ///
+ /// \brief Set RelativeTolerance.
+ ///
+ /// \param newRelativeTolerance New relative tolerance value
+ ///
+ void setRelativeTolerance(const double newRelativeTolerance)
+ {
+ this->relativeTolerance = newRelativeTolerance;
+ }
+
+ ///
+ /// \brief Get RelativeTolerance. Returns current relative tolerance value.
+ ///
+ double getRelativeTolerance() const
+ {
+ return this->relativeTolerance;
+ }
+
+ ///
+ /// \brief Set Gamma. Sets the gamma parameter used to update the forcing term.
+ ///
+ /// \param newGamma New gamma value.
+ ///
+ void setGamma(const double newGamma)
+ {
+ this->gamma = newGamma;
+ }
+
+ ///
+ /// \brief Get Gamma. Returns current gamma value.
+ ///
+ double getGamma() const
+ {
+ return this->gamma;
+ }
+
+ ///
+ /// \brief Set EtaMax. Maximum error tolerance for residual in inner
+ /// iteration. The inner iteration terminates when the relative linear residual is
+ /// smaller than eta*| F(x_c) |. eta is determined by the modified
+ /// Eisenstat-Walker formula if etamax > 0. If etamax < 0, then eta = |etamax|
+ /// for the entire iteration.
+ /// \see updateForcingTerm()
+ ///
+ void setEtaMax(const size_t newEtaMax)
+ {
+ this->etaMax = newEtaMax;
+ }
+
+ ///
+ /// \brief Get EtaMax. Returns current etaMax value.
+ ///
+ size_t getEtaMax() const
+ {
+ return this->etaMax;
+ }
+
+ ///
+ /// \brief Set MaxIterations. The maximum number of nonlinear iterations.
+ ///
+ void setMaxIterations(const size_t newMaxIterations)
+ {
+ this->maxIterations = newMaxIterations;
+ }
+
+ ///
+ /// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
+ ///
+ size_t getMaxIterations() const
+ {
+ return this->maxIterations;
+ }
+
+ ///
+ /// \brief Set the useArmijo flag. If useArmijo is true a line search is performed
+ /// using the Armijo-Goldstein condition.
+ ///
+ void setUseArmijo(const bool value)
+ {
+ this->useArmijo = value;
+ if(!value)
+ {
+ this->armijoMax = 0;
+ }
+ else
+ {
+ this->armijoMax = 30;
+ }
+ }
+
+ ///
+ /// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
+ ///
+ bool getUseArmijo() const
+ {
+ return this->useArmijo;
+ }
+
+ ///
+ /// \brief Set the useArmijo flag. If useArmijo is true a line search is performed
+ /// using the Armijo-Goldstein condition.
+ ///
+ void setForcingTerm(const double value)
+ {
+ this->forcingTerm = value;
+ }
+
+ ///
+ /// \brief Get MaxIterations. Returns current maximum nonlinear iterations.
+ ///
+ double getForcingTerm() const
+ {
+ return this->forcingTerm;
+ }
+
+
+private:
+ std::shared_ptr<LinearSolverType> linearSolver; ///> Linear solver to use. Default: Conjugate gradient.
+ JacobianType jacobian; ///> Jacobian matrix function.
+ double forcingTerm; ///> Method's forcing term (Default: 0.9).
+ double absoluteTolerance; ///> Tolerance for the method (Default: 1.0e-3).
+ double relativeTolerance; ///> Relative (to the rhs) tolerance (Default: 1.0e-6).
+ double gamma; ///> Internal parameter used to update the forcing term (Default: 0.9).
+ double etaMax; ///> Maximum tolerance for the linear solver (Default: 0.9).
+ size_t maxIterations; ///> Maximum number of nonlinear iterations (Default: 50).
+ bool useArmijo; ///> True if Armijo liner search is desired (Default: true).
+ std::vector<double> fnorms; ///> Consecutive function norms
+};
+
+#endif // INEXACT_NEWTON_H
diff --git a/Solvers/NonLinearSolver.cpp b/Solvers/NonLinearSolver.cpp
index 2550fb39e..c2006607b 100644
--- a/Solvers/NonLinearSolver.cpp
+++ b/Solvers/NonLinearSolver.cpp
@@ -35,10 +35,9 @@ NonLinearSolver::NonLinearSolver():
}
//---------------------------------------------------------------------------
-double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x)
+double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x, const double previousFnorm)
{
/// Temporaries used in the line search
- double previousFnorm = this->nonLinearSystem->getFunctionValue().norm();
std::array<double, 3> fnormSqr = {previousFnorm*previousFnorm, 0.0, 0.0};
std::array<double, 3> lambda = {this->sigma[0]*this->sigma[1], 1.0, 1.0};
@@ -51,8 +50,12 @@ double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x)
return currentFnorm;
}
+ // Save iterate in case this fails
+ auto x_old = x;
+
// Starts Armijo line search loop
- for(size_t i = 0; i < this->armijoMax; ++i)
+ size_t i;
+ for(i = 0; i < this->armijoMax; ++i)
{
/// Update x and keep books on lambda
this->updateIterate(-lambda[0]*dx,x);
@@ -75,6 +78,13 @@ double NonLinearSolver::armijo(const core::Vectord &dx, core::Vectord &x)
this->parabolicModel(fnormSqr, lambda);
}
+ if(i == this->armijoMax)
+ {
+ // TODO: Add to logger
+ std::cout << "Maximum number of Armijo iterations reached." << std::endl;
+ }
+
+
return currentFnorm;
}
diff --git a/Solvers/NonLinearSolver.h b/Solvers/NonLinearSolver.h
index 7fe77375b..9bfe111e0 100644
--- a/Solvers/NonLinearSolver.h
+++ b/Solvers/NonLinearSolver.h
@@ -60,7 +60,7 @@ public:
/// \param dx Computed direction.
/// \param x Current iterate.
///
- double armijo(const core::Vectord &dx, core::Vectord &x);
+ double armijo(const core::Vectord &dx, core::Vectord &x, const double previousFnorm);
///
/// \brief Three-point safeguarded parabolic model for a line search. Upon return
--
GitLab