NelderMead.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco, Mattia Piazza and Enrico Bertolazzi.                       *
 *                                                                                               *
 * The Optimist project is distributed under the BSD 2-Clause License.                           *
 *                                                                                               *
 * Davide Stocco                          Mattia Piazza                        Enrico Bertolazzi *
 * University of Trento               University of Trento                  University of Trento *
 * davide.stocco@unitn.it            mattia.piazza@unitn.it           enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef OPTIMIST_OPTIMIZER_NELDER_MEAD_HH
#define OPTIMIST_OPTIMIZER_NELDER_MEAD_HH
 
#include "Optimist/Optimizer.hh"
 
namespace Optimist
{
  namespace Optimizer
  {
 
    /*\
     |   _   _      _     _           __  __                _
     |  | \ | | ___| | __| | ___ _ __|  \/  | ___  __ _  __| |
     |  |  \| |/ _ \ |/ _` |/ _ \ '__| |\/| |/ _ \/ _` |/ _` |
     |  | |\  |  __/ | (_| |  __/ |  | |  | |  __/ (_| | (_| |
     |  |_| \_|\___|_|\__,_|\___|_|  |_|  |_|\___|\__,_|\__,_|
     |
    \*/
 

    template <typename Real, Integer N>
    class NelderMead : public Optimizer<Real, N, NelderMead<Real, N>>
    {
    public:
      static constexpr bool requires_function{true};
      static constexpr bool requires_first_derivative{false};
      static constexpr bool requires_second_derivative{false};
 
      using typename Optimizer<Real, N, NelderMead<Real, N>>::Vector;
      using typename Optimizer<Real, N, NelderMead<Real, N>>::Matrix;
 
      using Simplex = std::vector<Vector>;
 
    private:
      Simplex m_s;          
      Vector  m_c;          
      Real    m_alpha{1.0}; 
      Real    m_gamma{2.0}; 
      Real    m_rho{0.5};   
      Real    m_sigma{0.5}; 
 
      public:
      NelderMead() {
        this->m_s.reserve(N+1);
      }

      std::string name_impl() const {return "NelderMead";}

      void initialize(Vector const & x_ini, Real delta)
      {
        this->m_s.clear();
        this->m_s.push_back(x_ini);
        Vector x;
        for (Integer i{0}; i < N; ++i) {
          x = x_ini;
          x(i) += delta;
          this->m_s.push_back(x);
        }
        this->m_c = this->centroid();
      }

      Vector centroid() const
      {
        return std::accumulate(this->m_s.begin(), this->m_s.end(), this->m_c) / N;
      }

      template <typename FunctionLambda>
      void sort(FunctionLambda && function)
      {
        std::sort(this->m_s.begin(), this->m_s.end(),
          [function, this] (Vector const & a, Vector const & b) {
            Real function_a, function_b;
            this->evaluate_function(std::forward<FunctionLambda>(function), a, function_a);
            this->evaluate_function(std::forward<FunctionLambda>(function), b, function_b);
            return function_a < function_b;
          }
        );
        this->m_c = this->centroid();
      }

      Vector reflection() const
      {
        return this->m_c + this->m_alpha * (this->m_c - this->m_s.back());
      }

      Vector expansion(Vector const & x) const
      {
        return this->m_c + this->m_gamma * (x - this->m_c);
      }

      Vector outward_contraction(Vector const & x) const
      {
        return this->m_c + this->m_rho * (x - this->m_c);
      }

      Vector inward_contraction() const
      {
        return this->m_c + this->m_rho * (this->m_s.back() - this->m_c);
      }

      void shrink()
      {
        for (Integer i{1}; i < static_cast<Integer>(this->m_s.size()); ++i) {
          this->m_s.at(i) = this->m_s.at(1) + this->m_sigma * (this->m_s.at(i) - this->m_s.at(1));
        }
      }

      bool solve_impl(FunctionLambda && function, Vector const & x_ini, Vector & x_sol)
      {
        // Setup internal variables
        this->reset();
 
        // Print header
        if (this->m_verbose) {this->header();}
 
        // Initialize variables
        Vector x_r, x_c, x_e, simplex_size;
        Real function_x_r, function_x_1, function_x_last, function_x_e, function_x_c;
 
        // Set initial iteration
        this->initialize(x_ini, 100.0);
 
        // Algorithm iterations
        for (this->m_iterations = 1; this->m_iterations < this->m_max_iterations; ++this->m_iterations)
        {
          // Sort the simplex points based on the function values
          this->sort(function);
 
          // Store trace
          this->store_trace(this->m_s.back());
 
          // Compute the reflection point
          x_r = this->reflection();
          this->evaluate_function(std::forward<FunctionLambda>(function), x_r, function_x_r);
          this->evaluate_function(std::forward<FunctionLambda>(function), this->m_s.at(1), function_x_1);
          this->evaluate_function(std::forward<FunctionLambda>(function), this->m_s.back(), function_x_last);
 
          if (function_x_1 <= function_x_r && function_x_r < function_x_last)
          {
            this->m_s.back() = x_r;
            if (this->m_verbose) {this->info(function_x_1, "Reflection");}
          }
          else if (function_x_r < function_x_1)
          {
            // Expansion
            x_e = this->expansion(x_r);
            this->evaluate_function(std::forward<FunctionLambda>(function), x_e, function_x_e);
            if (function_x_e < function_x_r) {
              this->m_s.back() = x_e;
            } else {
              this->m_s.back() = x_r;
            }
            if (this->m_verbose) {this->info(function_x_1, "Expansion");}
          }
          else if (function_x_r < function_x_last)
          {
            // Outward contraction
            x_c = this->outward_contraction(x_r);
            this->evaluate_function(std::forward<FunctionLambda>(function), x_c, function_x_c);
            if (function_x_c < function_x_r) {
              this->m_s.back() = x_c;
            } else {
              this->shrink();
            }
            if (this->m_verbose) {this->info(function_x_1, "Outward contraction");}
          }
          else
          {
            // Inward contraction
            x_c = this->inward_contraction();
            this->evaluate_function(std::forward<FunctionLambda>(function), x_c, function_x_c);
            if (function_x_c < function_x_last) {
              this->m_s.back() = x_c;
            } else {
              this->shrink();
            }
            if (this->m_verbose) {this->info(function_x_1, "Inward contraction");}
          }
          simplex_size = this->m_s.back() - this->m_s.front();
          if (simplex_size.norm() < this->m_tolerance) {
            this->m_converged = true;
            break;
          }
        }
 
        // Print bottom
        if (this->m_verbose) {this->bottom();}
 
        // Convergence data
        x_sol = this->m_s.at(0);
        return this->m_converged;
      }
 
    }; // class NelderMead
 
  } // namespace Optimizer
 
} // namespace Optimist
 
#endif // OPTIMIST_OPTIMIZER_NELDER_MEAD_HH