From bd02d57d82c25b967dc4c9c2a1b8df8d002bdf03 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 9 May 2018 05:40:08 +0200
Subject: [PATCH] Use Newton instead of Trust-Region for the Geodesic FE min
problems
In theory this is dangerous, because a pure Newton method will
only converge locally. However, the trust-region part had been
partially disabled anyway, and I have not noticed any problems.
Also, I think that from my well-posedness result for higher-order
GeoFEs follows that a Newton method will always converge if the
problem is well-posed (that is how the proof works, after all).
On the plus side, this patch brings roughly 5% speed increas.
---
dune/gfe/localgeodesicfefunction.hh | 5 +-
dune/gfe/targetspacertrsolver.cc | 162 +++++-----------------------
dune/gfe/targetspacertrsolver.hh | 58 ++--------
3 files changed, 33 insertions(+), 192 deletions(-)
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 1552f595..39602131 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -195,10 +195,7 @@ evaluate(const Dune::FieldVector<ctype, dim>& local) const
solver.setup(&assembler,
initialIterate,
1e-14, // tolerance
- 100, // maxTrustRegionSteps
- 2, // initial trust region radius
- 100, // inner iterations
- 1e-14 // inner tolerance
+ 100 // maxNewtonSteps
);
solver.solve();
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index 2e11e9ea..b00b59e1 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -13,39 +13,14 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::
setup(const AverageDistanceAssembler<TargetSpace>* assembler,
const TargetSpace& x,
double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int innerIterations,
- double innerTolerance)
+ int maxNewtonSteps)
{
assembler_ = assembler;
x_ = x;
tolerance_ = tolerance;
- maxTrustRegionSteps_ = maxTrustRegionSteps;
- initialTrustRegionRadius_ = initialTrustRegionRadius;
- innerIterations_ = innerIterations;
- innerTolerance_ = innerTolerance;
+ maxNewtonSteps_ = maxNewtonSteps;
this->verbosity_ = NumProc::QUIET;
minNumberOfIterations_ = 1;
-
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- // ////////////////////////////////
- // Create a projected gauss-seidel solver
- // ////////////////////////////////
-
- // First create a Gauss-seidel base solver
- innerSolverStep_ = std::auto_ptr<TrustRegionGSStep<MatrixType, CorrectionType> >(new TrustRegionGSStep<MatrixType, CorrectionType>);
-
- energyNorm_ = std::auto_ptr<EnergyNorm<MatrixType, CorrectionType> >(new EnergyNorm<MatrixType, CorrectionType>(*innerSolverStep_.get()));
-
- innerSolver_ = std::auto_ptr< ::LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(innerSolverStep_.get(),
- innerIterations,
- innerTolerance,
- energyNorm_.get(),
- Solver::QUIET));
-
- innerSolver_->useRelativeError_ = false;
-#endif
}
@@ -54,68 +29,34 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
{
assert(minNumberOfIterations_ > 0);
- MaxNormTrustRegion<blocksize,field_type> trustRegion(1, // we have only one block
- initialTrustRegionRadius_);
-
- field_type energy = assembler_->value(x_);
-
// /////////////////////////////////////////////////////
- // Trust-Region Solver
+ // Newton Solver
// /////////////////////////////////////////////////////
- for (size_t i=0; i<maxTrustRegionSteps_; i++) {
+ for (size_t i=0; i<maxNewtonSteps_; i++) {
if (this->verbosity_ == Solver::FULL) {
std::cout << "----------------------------------------------------" << std::endl;
- std::cout << " Trust-Region Step Number: " << i
- << ", radius: " << trustRegion.radius()
- << ", energy: " << energy << std::endl;
+ std::cout << " Newton Step Number: " << i
+ << ", energy: " << assembler_->value(x_) << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
}
- CorrectionType rhs(1); // length is 1 _block_
- CorrectionType corr(1); // length is 1 _block_
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- corr = 0;
-#endif
+ CorrectionType corr;
- MatrixType hesseMatrix(1,1);
+ MatrixType hesseMatrix;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- assembler_->assembleGradient(x_, rhs[0]);
- assembler_->assembleHessian(x_, hesseMatrix[0][0]);
-#else
- /** \todo Fix this sense copying */
- typename TargetSpace::EmbeddedTangentVector foo;
- assembler_->assembleEmbeddedGradient(x_, foo);
- rhs[0] = foo;
- assembler_->assembleEmbeddedHessian(x_, hesseMatrix[0][0]);
-#endif
+ typename TargetSpace::EmbeddedTangentVector rhs;
+ assembler_->assembleEmbeddedGradient(x_, rhs);
+ assembler_->assembleEmbeddedHessian(x_, hesseMatrix);
// The right hand side is the _negative_ gradient
rhs *= -1;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- dynamic_cast<LinearIterationStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->setProblem(hesseMatrix, corr, rhs);
-
- dynamic_cast<TrustRegionGSStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->obstacles_ = &trustRegion.obstacles();
-
- innerSolver_->preprocess();
-#endif
// /////////////////////////////
// Solve !
// /////////////////////////////
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- innerSolver_->solve();
-#else
Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
- GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
-#endif
-
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- corr = innerSolver_->iterationStep_->getSol();
-#endif
-
- //std::cout << "Corr: " << corr << std::endl;
+ GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix, corr, rhs, basis);
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
@@ -125,84 +66,31 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ std::cout << i+1 << " Newton steps were taken." << std::endl;
break;
}
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
-
+#if 0
TargetSpace newIterate = x_;
- newIterate = TargetSpace::exp(newIterate, corr[0]);
-
- field_type oldEnergy = energy;
- field_type energy = assembler_->value(newIterate);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- CorrectionType tmp(corr.size());
- tmp = 0;
- hesseMatrix.umv(corr, tmp);
- field_type modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-#else
- field_type modelDecrease = (rhs*corr) - 0.5 * hesseMatrix[0][0].energyScalarProduct(corr[0],corr[0]);
-#endif
+ newIterate = TargetSpace::exp(newIterate, corr);
- if (this->verbosity_ == NumProc::FULL) {
- std::cout << "Absolute model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "Relative model decrease: " << modelDecrease / energy
- << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
- }
-
- assert(modelDecrease >= -1e-15);
+ if (this->verbosity_ == NumProc::FULL)
+ {
+ field_type oldEnergy = assembler_->value(x_);
+ field_type energy = assembler_->value(newIterate);
- if (energy >= oldEnergy) {
- if (this->verbosity_ == NumProc::FULL)
+ if (energy >= oldEnergy)
printf("Richtung ist keine Abstiegsrichtung!\n");
}
- if (energy >= oldEnergy &&
- i>minNumberOfIterations_-1 &&
- (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Suspecting rounding problems" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
-
- x_ = newIterate;
- break;
- }
-
- // //////////////////////////////////////////////
- // Check for acceptance of the step
- // //////////////////////////////////////////////
- if ( (oldEnergy-energy) / modelDecrease > 0.9) {
- // very successful iteration
-
- x_ = newIterate;
- trustRegion.scale(2);
-
- } else if ( (oldEnergy-energy) / modelDecrease > 0.01
- || std::abs(oldEnergy-energy) < 1e-12) {
- // successful iteration
- x_ = newIterate;
-
- } else {
- // unsuccessful iteration
- trustRegion.scale(0.5);
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Unsuccessful iteration!" << std::endl;
- }
-
- // Write current energy
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "--- Current energy: " << energy << " ---" << std::endl;
-
+ // Actually take the step
+ x_ = newIterate;
+#else
+ x_ = TargetSpace::exp(x_, corr);
+#endif
}
}
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 10e2d87f..568d8314 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -1,17 +1,7 @@
#ifndef TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
#define TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
-
-//#define USE_GAUSS_SEIDEL_SOLVER
-
-#include <dune/istl/matrix.hh>
-
-#include <dune/solvers/common/boxconstraint.hh>
-#ifdef USE_GAUSS_SEIDEL_SOLVER
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
-#endif
-#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/common/numproc.hh>
#include <dune/gfe/symmetricmatrix.hh>
@@ -28,31 +18,16 @@ class TargetSpaceRiemannianTRSolver
// Centralize the field type here
typedef typename TargetSpace::ctype field_type;
- // Some types that I need
- // The types have the dynamic outer type because the dune-solvers solvers expect
- // this sort of type.
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- typedef Dune::Matrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
- typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> > CorrectionType;
-#else
- typedef Dune::Matrix<Dune::SymmetricMatrix<field_type, embeddedBlocksize> > MatrixType;
- typedef Dune::BlockVector<Dune::FieldVector<field_type, embeddedBlocksize> > CorrectionType;
-#endif
+ typedef Dune::SymmetricMatrix<field_type, embeddedBlocksize> MatrixType;
+ typedef Dune::FieldVector<field_type, embeddedBlocksize> CorrectionType;
public:
- TargetSpaceRiemannianTRSolver()
- // : IterativeSolver<std::vector<TargetSpace>, Dune::BitSetVector<blocksize> >(0,100,NumProc::FULL)
- {}
-
/** \brief Set up the solver using a monotone multigrid method as the inner solver */
void setup(const AverageDistanceAssembler<TargetSpace>* assembler,
const TargetSpace& x,
double tolerance,
- int maxTrustRegionSteps,
- double initialTrustRegionRadius,
- int innerIterations,
- double innerTolerance);
+ int maxNewtonSteps);
void solve();
@@ -70,32 +45,13 @@ protected:
/** \brief Tolerance of the solver */
double tolerance_;
- /** \brief The initial trust-region radius in the maximum-norm */
- double initialTrustRegionRadius_;
-
- /** \brief Maximum number of trust-region steps */
- size_t maxTrustRegionSteps_;
-
- /** \brief Maximum number of multigrid iterations */
- int innerIterations_;
-
- /** \brief Error tolerance of the multigrid QP solver */
- double innerTolerance_;
+ /** \brief Maximum number of Newton steps */
+ size_t maxNewtonSteps_;
/** \brief The assembler for the average-distance functional */
const AverageDistanceAssembler<TargetSpace>* assembler_;
-#ifdef USE_GAUSS_SEIDEL_SOLVER
- /** \brief The solver for the quadratic inner problems */
- std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;
-
- /** \brief The iteration step for the quadratic inner problems */
- std::auto_ptr<TrustRegionGSStep<MatrixType, CorrectionType> > innerSolverStep_;
-
- /** \brief Norm for the quadratic inner problems */
- std::auto_ptr<EnergyNorm<MatrixType, CorrectionType> > energyNorm_;
-#endif
- /** \brief Specify a minimal number of iterations the trust-region solver has to do
+ /** \brief Specify a minimal number of iterations the Newton solver has to do
*
* This is needed when working with automatic differentiation. While a very low
* number of iterations may be enough to precisely compute the value of a
--
GitLab