I had to port a simplex/Nelder-Mead optimizer that I already have in Python in C++. As for the Python version, I tried to be as generic as possible but as efficient as possible, so the state is no longer a dictionary, but a simple structure.

I could have used the Numerical Recipes version, but the licence cost is not worth it, and the code is not generic enough, not explained enough. And also there are some design decisions that are questionable (one method = one responsibility).

Design

The main choice was to use Eigen for the inner computation in the code. This library is perfectly fitted for numerical computation in C++. The parameter set for the function could be any matrix class, but Eigen’s ones are ensured to be compatible inside the code. The simplex is stored as an Eigen dynamic array to cover any dimension possible cases.

There are several tools that are pretty much generic in a simplex algorithm: the state (current parameters, iteration…) and the stopping criterion (a conjunction of a maximum iteration and a relative value difference). So these concepts will be implemented as separate generic classes.

The main class, Simplex, will encompass the actual simplex algorithm. There are three inner central methods: the one that creates the first state (initialize_polytope()), the one that will detect the best parameter set, the worst one and the worst second (find_best_worst_near_worst()) and the one that create a new candidate based on the current simplex and the selected worst parameter set (create_new_parameters()).

The main method is optimize(), building on the class attributes that were set before and the main methods that were exposed before. The algorithm follows strictly the implementation on the Wikipedia page. The sizes of expansion, contraction could have been made template, or even modified at compile time, it is an exercise left for those that want a fully generic class.

A static builder was written to ease the creation of a new instance. It doesn’t do anything else than instantiating the class with the constructor, its only goal is to remove the need for knowing the exact template parameter and to fully specify them. This can be done because of the auto keyword in C++11.

Using it

Let’s start with a function to optimize. I’ve chosen the Rosenbrock function as it is not complex to compute but still difficult to optimize. For the parameter kind, I’ve also used Eigen.

typedef Eigen::Vector2f ParameterTypeRosenbrock;
 
// Rosenbrock cost function
struct Rosenbrock
{
  typedef float DataType;
  typedef ParameterTypeRosenbrock ParameterType;
  float operator()(const ParameterTypeRosenbrock& parameters) const
  {
    return (parameters(1, 0) - parameters(0, 0) * parameters(0, 0)) * (parameters(1, 0) - parameters(0, 0) * parameters(0, 0))
        + (1 - parameters(0, 0)) * (1 - parameters(0, 0));
  }
};

Now, I can create the optimizer with all the different parameters.

int main(int argc, char** argv)
{
  ParameterTypeRosenbrock start_point;
  start_point << 10, 10;
 
  Rosenbrock fun;
 
  long max_iterations = 1000;
  float ftol = 0.00001;
 
  auto optimizer = Optimization::Local::build_simplex( // Builder to generate the correct optimizer
      fun,                                             // Used to infer function type
      Optimization::Local::make_and_criteria(Optimization::Local::IterationCriterion(max_iterations),
          Optimization::Local::RelativeValueCriterion<float>(ftol))); // Composite stoping criterion
 
  optimizer.set_start_point(start_point);  // Starting parameters
  optimizer.set_delta(1);                  // Simplex size
  optimizer.optimize(fun);                 // Optimization start
 
  std::cout << "Starting point: " << start_point << std::endl;
  std::cout << "Starting value: " << fun(start_point) << std::endl;
  std::cout << "Best point: " << optimizer.get_best_parameters() << std::endl; // Retrieve the best parameters
  std::cout << "Best value: " << optimizer.get_best_value() << std::endl;      // Retrieve the best value
}

With C++11, it is easy to find the type of the optimizer, thanks to the builder. It is still cumbersome as the builder function needs the function instance to infer the inner data type as well as the parameter set type.

Starting point: 10
10
Starting value: 8181
Best point: 1
1
Best value: 5.68434e-014

The simplex does not need such a complex generix framework. The implementation usually uses the criterion, so I didn’t have to implement it the way I did. Still, it is interesting for a more complex optimizer like a quasi-Newton or conjugate-gradient.

Conclusion

This code is not world shaking but it explains how I consider optimization. The other optimizers from scikits.optimization could be written in the same way and they could be better understood by someone versed in object oriented languages than usual procedural code.

As usual, it can be found on github: https://github.com/mbrucher/CppOptimization