A first prototype of an standard trust-region solver using the infinity-norm...

A first prototype of an standard trust-region solver using the infinity-norm for the trust-reegion constraints.
It can be used to solve smooth nonlinear (nonconvex) unconstrained minimization problems

dune/solvers/solvers/CMakeLists.txt

@@ -9,4 +9,6 @@ install(FILES
     solver.hh
     tcgsolver.cc
     tcgsolver.hh
+    trustregionsolver.cc
+    trustregionsolver.hh
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers/solvers)

dune/solvers/solvers/Makefile.am

@@ -3,7 +3,7 @@ SUBDIRS =
 solversdir = $(includedir)/dune/solvers/solvers
 solvers_HEADERS = cgsolver.cc cgsolver.hh iterativesolver.cc iterativesolver.hh \
                   loopsolver.cc loopsolver.hh quadraticipopt.hh solver.hh \
-                  tcgsolver.cc tcgsolver.hh
+                  tcgsolver.cc tcgsolver.hh trustregionsolver.cc trustregionsolver.hh
 
 EXTRA_DIST = CMakeLists.txt
 

dune/solvers/solvers/trustregionsolver.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/io.hh>
+
+// Using a monotone multigrid as the inner solver
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+template <class ProblemType, class VectorType, class MatrixType>
+void TrustRegionSolver<ProblemType,VectorType,MatrixType>::solve()
+{
+    MonotoneMGStep<MatrixType,VectorType>* mgStep = nullptr;
+
+    // if the inner solver is a monotone multigrid set up a max-norm trust-region
+    if (dynamic_cast<::LoopSolver<VectorType>*>(innerSolver_.get())) {
+        mgStep = dynamic_cast<MonotoneMGStep<MatrixType,VectorType>*>(dynamic_cast<::LoopSolver<VectorType>*>(innerSolver_.get())->iterationStep_);
+
+    }
+
+    // setup the trust-region in the maximums norm
+    std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles(x_.size());
+    for (size_t i=0; i<trustRegionObstacles.size(); i++)
+        for (int j=0; j<blocksize; j++) {
+            trustRegionObstacles[i].lower(j) = -initialRadius_;
+            trustRegionObstacles[i].upper(j) = initialRadius_;
+        }
+
+    mgStep->obstacles_ = &trustRegionObstacles;
+
+    // /////////////////////////////////////////////////////
+    //   Trust-Region Solver
+    // /////////////////////////////////////////////////////
+
+    double oldEnergy = problem_->energy(x_);
+
+    for (int i=0; i<this->maxIterations_; i++) {
+
+        Dune::Timer totalTimer;
+        if (this->verbosity_ == Solver::FULL) {
+            std::cout << "----------------------------------------------------" << std::endl;
+            std::cout << "      Trust-Region Step Number: " << i
+                      << ",     radius: " << initialRadius_
+                      << ",     energy: " << oldEnergy << std::endl;
+            std::cout << "----------------------------------------------------" << std::endl;
+        }
+
+        // Assemble the quadratic problem
+        Dune::Timer gradientTimer;
+        problem_->assembleQP(x_);
+
+        // Compute gradient norm to monitor convergence
+        VectorType gradient = problem_->f_;
+        for (size_t j=0; j<gradient.size(); j++)
+            for (int k=0; k<gradient[j].size(); k++)
+                if ((*mgStep->ignoreNodes_)[j][k])
+                    gradient[j][k] = 0;
+
+        if (this->verbosity_ == Solver::FULL)
+            std::cout << "Gradient norm: " << this->errorNorm_->operator()(gradient) << std::endl;
+
+        if (this->verbosity_ == Solver::FULL)
+            std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
+        VectorType corr(x_.size());
+        corr = 0;
+
+        mgStep->setProblem(problem_->A_, corr, problem_->f_);
+
+        innerSolver_->preprocess();
+
+        ///////////////////////////////
+        //    Solve !
+        ///////////////////////////////
+
+        // Solve the subproblem until we find an acceptable iterate
+        if (this->verbosity_ == NumProc::FULL)
+            std::cout << "Solve quadratic problem..." << std::endl;
+        while(true) {
+
+            Dune::Timer solutionTimer;
+            innerSolver_->solve();
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+
+            corr = mgStep->getSol();
+
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+            if (corr.infinity_norm() < this->tolerance_) {
+                if (this->verbosity_ != NumProc::QUIET)
+                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
+                return;
+            }
+
+            // ////////////////////////////////////////////////////
+            //   Check whether trust-region step can be accepted
+            // ////////////////////////////////////////////////////
+
+            VectorType newIterate = x_;
+            newIterate += corr;
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Infinity norm of new iterate : " << newIterate.infinity_norm() << std::endl;
+
+            field_type energy    = problem_->energy(newIterate);
+
+            // compute the model decrease
+            field_type modelDecrease = problem_->modelDecrease(corr);
+
+            field_type relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+            if (this->verbosity_ == NumProc::FULL) {
+                std::cout << "Absolute model decrease: " << modelDecrease
+                    << ",  functional decrease: " << oldEnergy - energy << std::endl;
+                std::cout << "Relative model decrease: " << relativeModelDecrease
+                    << ",  functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+            }
+
+            // assure that the inner solver works correctly
+            assert(modelDecrease >= 0);
+
+            if (energy >= oldEnergy) {
+                if (this->verbosity_ == NumProc::FULL)
+                    printf("No energy descent!\n");
+            }
+
+            if (energy >= oldEnergy &&
+                    (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
+                if (this->verbosity_ == NumProc::FULL)
+                    std::cout << "Suspecting rounding problems, stopping here!" << std::endl;
+
+                if (this->verbosity_ != NumProc::QUIET)
+                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+                x_ = newIterate;
+                return;
+            }
+
+            // //////////////////////////////////////////////
+            //   Check for acceptance of the step
+            // //////////////////////////////////////////////
+            if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+                // very successful iteration
+
+                x_ = newIterate;
+                scaleTrustRegionRadius(trustRegionObstacles,enlargeRadius_);
+                std::cout << "Very Successful iteration, enlarging radius! " <<initialRadius_<<std::endl;
+
+                // current energy becomes 'oldEnergy' for the next iteration
+                oldEnergy = energy;
+                break;
+
+            } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+                    || std::abs(oldEnergy-energy) < 1e-12) {
+                // successful iteration
+                x_ = newIterate;
+
+                // current energy becomes 'oldEnergy' for the next iteration
+                oldEnergy = energy;
+                break;
+
+            } else {
+
+                // unsuccessful iteration
+
+                // Decrease the trust-region radius
+                scaleTrustRegionRadius(trustRegionObstacles,shrinkRadius_);
+
+                if (this->verbosity_ == NumProc::FULL)
+                    std::cout << "Unsuccessful iteration, shrinking radius! " <<initialRadius_<<std::endl;
+            }
+        }
+        std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+    }
+
+}

dune/solvers/solvers/trustregionsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_TRUST_REGION_SOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_TRUST_REGION_SOLVER_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+/** \brief A standard Trust-region solver with infinity-norm for the trust-region constraints, to solve an unconstrained nonlinear (possibly nonconvex) minimization problem. 
+ *
+ *  The energy functional is provided by the ProblemType class. It has to provide the following methods:
+ *
+ *      void assembleQP(const VectorType& iterate)
+ *      field_type energy(const VectorType& iterate)
+ *      field_type modelDecrease(const VectorType& iterate)
+ *
+ *      and member: 
+ *          MatrixType A_
+ *          VectorType f_
+ *
+ *      The method 'assembleQP' should assemble the local quadratic approximation at 'iterate'
+ *      The local model is to be given by the quadratic part A_ and the linear part f_
+ *
+ *      The method 'energy' should evaluate the nonlinear energy at 'iterate'.
+ *
+ *      The method 'modelDecrease' should compute the decrease in the model-energy. 
+ *
+ *  The inner solver that has to be provided should be able to solve non-convex box-constrained quadratic problems,
+ *  if the nonlinear functional is non-convex itself.
+ *
+ *  \tparam ProblemType A class representing the energy functional
+ *
+ * */
+template <class ProblemType, class VectorType, class MatrixType>
+class TrustRegionSolver
+    : public IterativeSolver<VectorType>
+{
+    const static int blocksize = VectorType::value_type::dimension;
+    typedef IterativeSolver<VectorType> Base;
+    typedef typename VectorType::field_type field_type;
+
+public:
+
+    TrustRegionSolver()
+        : Base(0,100,NumProc::FULL), enlargeRadius_(2), shrinkRadius_(0.5),
+          goodIteration_(0.9), badIteration_(0.01)
+    {}
+
+    TrustRegionSolver(field_type tolerance, int maxIterations, NumProc::VerbosityMode verbosity=NumProc::FULL)
+        : Base(tolerance, maxIterations, verbosity), enlargeRadius_(2), shrinkRadius_(0.5),
+          goodIteration_(0.9), badIteration_(0.01)
+    {}
+
+    TrustRegionSolver(const Dune::ParameterTree& trConfig)
+        : Base(0,100,NumProc::FULL)
+    {
+        setupTrParameter(trConfig);
+    }
+
+    //! Setup the standard trust-region parameter
+    void setupTrParameter(const Dune::ParameterTree& trConfig) {
+
+        this->tolerance_     = trConfig.get<field_type>("tolerance");
+        this->maxIterations_ = trConfig.get<int>("maxIterations");
+        this->verbosity_     = trConfig.get("verbosity",NumProc::FULL);
+        initialRadius_       = trConfig.get<field_type>("initialRadius");
+        maxRadius_           = trConfig.get<field_type>("maxRadius");
+        enlargeRadius_       = trConfig.get("enlargeRadius",2.0);
+        shrinkRadius_        = trConfig.get("shrinkRadius",0.5);
+        goodIteration_       = trConfig.get("goodIteration",0.9);
+        badIteration_        = trConfig.get("badIteration",0.01);
+
+    }
+
+
+    /** \brief Set up the solver. */
+    void setup(const Dune::ParameterTree& trConfig,
+               Norm<VectorType>& errorNorm,
+               ProblemType& problem, Solver& innerSolver) {
+
+        setupTrParameter(trConfig);
+        this->errorNorm_ = &errorNorm;
+        problem_ = Dune::stackobject_to_shared_ptr(problem);
+        innerSolver_ = Dune::stackobject_to_shared_ptr(innerSolver);
+    }
+
+    //! Solve the problem.
+    void solve();
+
+    //! Set the initial iterate
+    void setInitialIterate(const VectorType& x) {x_ = x;}
+
+    VectorType getSol() const {return x_;}
+
+protected:
+
+    //! Adjust the trust-region obstacles
+    void scaleTrustRegionRadius(std::vector<BoxConstraint<field_type,blocksize> >& obs, const field_type& scaling)
+    {
+        if (scaling*initialRadius_>maxRadius_)
+            initialRadius_ = maxRadius_;
+        else
+            initialRadius_ *= scaling;
+
+        for (size_t i=0; i<obs.size(); i++)
+            for (int j=0; j<blocksize; j++) {
+                obs[i].lower(j) = -initialRadius_;
+                obs[i].upper(j) = initialRadius_;
+        }
+    }
+
+    /** \brief The problem type. */
+    std::shared_ptr<ProblemType> problem_;
+
+    /** \brief The solution vector */
+    VectorType x_;
+
+    /** \brief The solver for the quadratic inner problems */
+    std::shared_ptr<Solver> innerSolver_;
+
+    /** \brief The initial trust-region radius in the maximum-norm */
+    field_type initialRadius_;
+
+    /** \brief The maximal trust-region radius in the maximum-norm */
+    field_type maxRadius_;
+
+    /** \brief If the iteration is very successfull we enlarge the radius by this factor.*/
+    field_type enlargeRadius_;
+
+    /** \brief If the iteration is unsuccessfull we shrink the radius by this factor.*/
+    field_type shrinkRadius_;
+
+    /** \brief If the ratio between predicted and achieved decrease is above this number, the iteration is very successful.*/
+    field_type goodIteration_;
+
+    /** \brief If the ratio between predicted and achieved decrease is lower than this number, the iteration is unsuccessful.*/
+    field_type badIteration_;
+};
+
+#include "trustregionsolver.cc"
+
+#endif