Commit 08043796 authored by Ricardo Ortiz's avatar Ricardo Ortiz
Browse files

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.
parent 41df08ca
Pipeline #4821 passed with stage
// 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;
}
// 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
......@@ -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;
}
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment