Nelder-Mead optimization in C++

Question

I've implemented the Nelder-Mead optimization algorithm in C++. I found this document to be a very good explanation of how the algorithm works, but I'll try my best to explain it (with an example) here so it's clear for everyone.

The problem I'm trying to optimize is very simple : I'm trying to find the best a,b,c coefficients for a quadratic function to fit points on a quadratic function of form f(x) = 1.5*x^2 + 3.2*x + 7. So, of course, the optimal parameters I'm looking for are a=1.5, b=3.2, c=7.

In order to compute the error between the coefficients that I have and the expected coefficients, I use the following function :

#include "Eigen/Dense"
static double quadraticEquationError(Array<double, 1, 3> parameters)
{
    double searchedParamters[] = {1.5, 3.2, 7};

    auto x = Array<double, 100, 1>().setLinSpaced(100, 0, 90);
    auto expected = searchedParamters[0] * x*x + searchedParamters[1] * x + searchedParamters[2];
    auto result = parameters(0, 0) * x * x + parameters(0, 1) * x + parameters(0, 2);

    auto diff = result.matrix() - expected.matrix();

    return diff.norm();
}

It's very simple, but it's also a way for me to check if my algorithm works. The goal of this error function is to return 0, where I would have found my optimal parameters.

The Nelder-Mead algorithm takes as input an error function and some initial parameters (in my case, a=0, b=0, c=0). Then, it generates a simplex of N+1 dimensions (here, N = 3 since we are looking for 3 parameters) where each vertex of the simplex represents a new parameters set. After that, depending on the values of the worst, second worst and best sets, it replaces the worst point by a new value.

There are five operations that can be used : reflection, expansion, inside contraction, outside contraction and shrinking. The following figure shows the result of these operations for a 2D simplex (of 3 points).

In order to figure out which operation to use, we can look at the following diagram :

Progressively, the error will go down until a convergence criterion is reached. In my case, I decided to stop when the best error is under a certain value.

So, I hope this is pretty clear as to how the algorithm works, here's my code :

main.cpp

#include "library.h"
#include <iostream>
#include "Eigen/Dense"
#include "optimization/NelderMead.h"

 using namespace Optimization;

static double quadraticEquationError(Array<double, 1, 3> parameters)
{
    double searchedParamters[] = {1.5, 3.2, 7};

    auto x = Array<double, 100, 1>().setLinSpaced(100, 0, 90);
    auto expected = searchedParamters[0] * x*x + searchedParamters[1] * x + searchedParamters[2];
    auto result = parameters(0, 0) * x * x + parameters(0, 1) * x + parameters(0, 2);

    auto diff = result.matrix() - expected.matrix();

    return diff.norm();
}

 int main() {
     Array<double, 1, 3> initialParameters;
     initialParameters.row(0) << 0.0, 0.0, 0.0;

     auto nelderMead = NelderMead<3>(quadraticEquationError, initialParameters, 1, 2);

     nelderMead.optimize();

     return 0;
}

optimization/NelderMead.h

#include "Eigen/Dense"
#include <vector>
#include <limits>
#include <cstring>
#include <iostream>

using namespace Eigen;

namespace Optimization
{
    template<int nDims>
    class NelderMead {
    public:
        typedef Array<double, 1, nDims> FunctionParameters;
        typedef std::function<double(FunctionParameters)> ErrorFunction;

        NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
                   double minError, double initialEdgeLength,
                   double shrinkCoeff = 1, double contractionCoeff = 0.5,
                   double reflectionCoeff = 1, double expansionCoeff = 1)
                : errorFunction(errorFunction), minError(minError),
                shrinkCoeff(shrinkCoeff), contractionCoeff(contractionCoeff),
                reflectionCoeff(reflectionCoeff), expansionCoeff(expansionCoeff),
                worstValueId(-1), secondWorstValueId(-1), bestValueId(-1)
        {
            this->errors = std::vector(nDims + 1, std::numeric_limits<double>::max());

            const double b = initialEdgeLength / (nDims * SQRT2) * (sqrt(nDims + 1) - 1);
            const double a = initialEdgeLength / SQRT2;

            this->values = initial.replicate(nDims + 1, 1);

            for (int i = 0; i < nDims; i++)
            {
                FunctionParameters simplexRow;
                simplexRow.setConstant(b);
                simplexRow(0, i) = a;
                simplexRow += initial;

                this->values.row(i+1) = simplexRow;
            }
        }

        void optimize()
        {
            for (int i = 0; i < nDims+1; i++)
            {
                this->errors.at(i) = this->errorFunction(this->values.row(i));
            }

            this->invalidateIdsCache();

            while (this->errors.at(this->bestValueId) > this->minError)
            {
                step();
                auto bestError = this->errorFunction(this->best());
                auto worstError = this->errorFunction(this->worst());

                std::cout << "Best error " << std::to_string(bestError) << " with : " << this->best();
                std::cout << " Worst error " << std::to_string(worstError) << " with : " << this->worst();
                std::cout << '\n';
            }
        }

        void step()
        {
            auto meanWithoutWorst = this->getMeanWithoutWorst();

            auto reflectionOfWorst = this->getReflectionOfWorst(meanWithoutWorst);
            auto reflectionError = this->errorFunction(reflectionOfWorst);

            FunctionParameters newValue = reflectionOfWorst;
            double newError = reflectionError;

            bool shrink = false;

            if (reflectionError < this->errors.at(this->bestValueId))
            {
                auto expansionValue = this->expansion(meanWithoutWorst, reflectionOfWorst);
                double expansionError = this->errorFunction(expansionValue);

                if (expansionError < this->errors.at(this->bestValueId))
                {
                    newValue = expansionValue;
                    newError = expansionError;
                }
            }
            else if (reflectionError > this->errors.at(this->worstValueId))
            {
                newValue = this->insideContraction(meanWithoutWorst);
                newError = this->errorFunction(newValue);

                if (newError > this->errors.at(this->worstValueId)) { shrink = true; }
            }
            else if (reflectionError > this->errors.at(this->secondWorstValueId))
            {
                newValue = this->outsideContraction(meanWithoutWorst);
                newError = this->errorFunction(newValue);

                if (newError > reflectionError) { shrink = true; }
            }
            else
            {
                newValue = reflectionOfWorst;
                newError = reflectionError;
            }

            if (shrink)
            {
                this->shrink();
                this->invalidateIdsCache();
                return;
            }

            this->values.row(this->worstValueId) = newValue;
            this->errors.at(this->worstValueId) = newError;
            this->invalidateIdsCache();
        }

        inline FunctionParameters worst() { return this->values.row(this->worstValueId); }
        inline FunctionParameters best() { return this->values.row(this->bestValueId); }

    private:

        void shrink()
        {
            auto bestVertex = this->values.row(this->bestValueId);

            for (int i = 0; i < nDims + 1; i++)
            {
                if (i == this->bestValueId) { continue; }

                this->values.row(i) = bestVertex + this->shrinkCoeff * (this->values.row(i) - bestVertex);
                this->errors.at(i) = this->errorFunction(this->values.row(i));
            }
        }

        inline FunctionParameters expansion(FunctionParameters meanWithoutWorst, FunctionParameters reflection)
        {
            return reflection + this->expansionCoeff * (reflection - meanWithoutWorst);
        }

        inline FunctionParameters insideContraction(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst - this->contractionCoeff * (meanWithoutWorst - this->worst());
        }

        inline FunctionParameters outsideContraction(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst + this->contractionCoeff * (meanWithoutWorst - this->worst());
        }

        FunctionParameters getReflectionOfWorst(FunctionParameters meanWithoutWorst)
        {
            return meanWithoutWorst + this->reflectionCoeff * (meanWithoutWorst - this->worst());
        }

        FunctionParameters getMeanWithoutWorst()
        {
            FunctionParameters mean(0);
            for (int i = 0; i < nDims + 1; i++)
            {
                if (i == this->worstValueId) { continue; }

                mean += this->values.row(i);
            }

            // Not divided by nDims+1 because there's one ignored value.
            mean /= nDims;

            return mean;
        }

        void invalidateIdsCache() {
            double worstError = std::numeric_limits<double>::min();
            int worstId = -1;
            double secondWorstError = std::numeric_limits<double>::max();
            int secondWorstId = -1;
            double bestError = std::numeric_limits<double>::max();
            int bestId = -1;

            for (int i = 0; i < nDims + 1; i++)
            {
                auto error = this->errors.at(i);

                if (error > worstError)
                {
                    secondWorstError = worstError;
                    secondWorstId = worstId;
                    worstError = error;
                    worstId = i;
                }
                else if (error > secondWorstError)
                {
                    secondWorstError = error;
                    secondWorstId = i;
                }

                if (error < bestError)
                {
                    bestError = error;
                    bestId = i;
                }
            }

            // If we deal with a problem in 1D, it won't be set.
            if (secondWorstId == -1)
            {
                secondWorstId = worstId;
            }

            this->bestValueId = bestId;
            this->worstValueId = worstId;
            this->secondWorstValueId = secondWorstId;
        }

        ErrorFunction errorFunction;
        Array<double, nDims + 1, nDims> values;
        std::vector<double> errors;

        int worstValueId;
        int secondWorstValueId;
        int bestValueId;

        double minError;

        double shrinkCoeff;
        double expansionCoeff;
        double contractionCoeff;
        double reflectionCoeff;

        const double SQRT2 = sqrt(2);
    };
}

Now, the thing is I don't know much about C++, but I know a lot about C# so I think my code structure is fine. I tried to use templates and typedefs to make code more readable/generic, but I'd like my C++ skills to improve a lot since I'm starting a new job soon that will require me to write C++ (my employer knows I'm not super strong in C++ but I'd like to get up to speed).

I'm looking for parts of my code that aren't using the full potential of C++, although I'm open to any part of the code reviewed.

Answer 1

Usage

Documentation

I'll consider the text contents of the post as documentation and it looks good.

Consumers

It is important to have a straightforward consumption process. Usually it is having a CMakeLists script that searches for required dependencies and sets up required (!) compilation flags. The stuff like optimization levels and warnings can go into toolchain file.

Usage flow

At the moment it looks like a toy program. What is the usecase for the current implementation? I would expect something like the following:

template <typename T, std::size_t Dimensions>
using coefficients_t = Array<T, 1, Dimensions>;

template <typename T, std::size_t Dimensions>
struct optimization_state_t {
    coefficients_t<T, Dimensions> coefficients;
    std::function<double(coefficients_t<T, Dimensions>)> error_function;
    /*errors and what other internal state is needed */
}; 

template <typename T, std::size_t Dimensions>
optimization_state_t make_starting_state(coefficients_t<T, Dimensions> start = {}, std::function<double(coefficients_t<T, Dimensions>)> error_function = default_error_function);

template <typename T, std::size_t Dimensions>
void perform_step(optimization_state_t<T, Dimensions>& state);

template <typename T, std::size_t Dimensions>
optimization_state_t<T, Dimensions> optimize_until(coefficients_t<T, Dimensions> coefficients);

This way, my flow will be

decide on initial parameters -> optimize one step and see if the convergence speed is good -> try again or abort -> feed results somewhere else

or

decide on initial parameters -> optimize until some condition -> feed results somewhere else

I believe interfaces should support desired usage flow instead of just providing bits and pieces that together form a crooked road.

Note that error function is part of the optimization state, as the error values will lose meaning in incompatible error functions. As Scott Meyers said, interfaces should be "easy to use correctly and hard to use incorrectly".

---

Details

Proper types

I'm not an experienced user of Eigen, I believe Vector3d is a better fit due to the conversion into matrix in the error function and generally solving a system of equations. There is also a library called Blaze, which from what I remember makes better use of intel MKL (needless to say that on AMD CPUs performance will be subpar).

Use standard constants

There are well defined constants since C++20. Computing it from scratch using library features might change the value from release to release.

Consider the cost of accepting by value

To accept by value, a new instance of the type has to be made. If the type does dynamic memory allocation (new, malloc, etc) and there is no real need to copy/move/construct a new object, then it is better to avoid paying the performance cost of additional memory allocation by accepting by reference. On the other hand, if Array<double, 1, 3> does not do dynamic memory allocation, it is better to accept it by value since it is small and will not cause stackoverflow

Do not use checking version of subscripts

Subscript operators that do bounds checking prevent SIMD instruction generation for compilers. It is better to just do plain assert and call non-checking version.

Do not import names into current scope without the need

using namespace something is usually a bad idea inside a header. Aside from obvious name clashes, there is ADL. The bottom line is that if Eigen decides to add similar function to yours and the call is not qualified with a namespace, two erroneous outcomes will happen. First and the best one is that ambigious call error will be issued because both functions are at equal rank in overload resolution. The worst case is that the compiler will silently select unintended function because the programmer didn't know about this (and most people don't know about ADL).

---

Style

Usage of this

this is usually frowned upon unless there are some clever mixin tricks need to be done. I do not really know a strong reason why it is frowned upon, but arguments thrown around are that it makes it easy to confuse argument with member variable, that this can be a bit weird with lambdas and so on.

Naming

The standard uses snake_case for almost everything except template parameter names which are Pascal case.

There is also Google style guide, MISRA and the others. I believe it is important to just pick a well-known one and stick with it, because clang-format configs are already written for them.

Answer 2

double searchedParamters[] = {1.5, 3.2, 7};

You don't change the values within the function, so it ought to be const and there's no reason to create a local (stack-based) variable and copy the three elements into it on every invocation. Make this constexpr.

---

Now, the thing is I don't know much about C++, but I know a lot about C#

So, I can guess that you'll be completely overlooking the existence of const, and will probably not consider that C++ has value semantics.

---

static double quadraticEquationError(Array<double, 1, 3> parameters)

I don't know what Array is (your code posted is incomplete), but I thing you should not be passing it by value.

---

typedef std::function<double(FunctionParameters)> ErrorFunction;
Use using ErrorFunction = std::function<double(FunctionParameters)>; rather than the old typedef syntax.

        NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
                   double minError, double initialEdgeLength,
                   double shrinkCoeff = 1, double contractionCoeff = 0.5,
                   double reflectionCoeff = 1, double expansionCoeff = 1)
                : errorFunction(errorFunction), minError(minError),
                shrinkCoeff(shrinkCoeff), contractionCoeff(contractionCoeff),
                reflectionCoeff(reflectionCoeff), expansionCoeff(expansionCoeff),
                worstValueId(-1), secondWorstValueId(-1), bestValueId(-1)

The final three members that have constant initializers should use inline direct initializers in the class, and then you don't need to name them here.
Uniform Initialization is to use curly braces, but parens still work.

this->errors = std::vector(nDims + 1, std::numeric_limits<double>::max()); etc.
Don't use this-> to access members.

this->errors.at(i) = this->errorFunction(this->values.row(i));
really, don't do that!

---

The class member: const double SQRT2 = sqrt(2);
should (at least) be static, as it's always the same so there is no reason for each instance to have its own. Having a const member will disable the autogenerated assignment operator.

---

inline FunctionParameters expansion(FunctionParameters meanWithoutWorst, FunctionParameters reflection)
Don't use the inline keyword for functions defined inside the class. I notice you only use that for some of them, anyway. Such functions are automatically inline. (inline means they can be in the .h file, BTW. Nothing more.)

---

Your sample main calls optimize but doesn't read out any results. As a synopsis of how the class is to be used, that's not very useful since that essentially doesn't do anything.

Answer 3

Naming things

FunctionParameters is a bad name. There are many functions and many parameters. A much better name is Coordinates, as the type represents the coordinates of a vertex. You could even consider naming it Vertex, but consider that not every coordinate is an actual vertex (like meanWithoutWorst).

Use verbs for function names, and prefer nouns for variable names. Also be consistent with naming. For example, instead of expansion(), write getExpansion(), so it's a verb and matches getMeanWithoutWorst().

Also avoid having both variables and functions with the same name, like shrink. I would recommend keeping shrink() as the function name, but renaming the variable to doShrink. It's hard to follow the noun/verb split here.

Rethink how you are templating your code

Instead of templating your algorithm on the number of dimensions, and then forcing coordinates of vertices to be Array<double, 1, Ndims>, consider that the Nelder-Mead algorithm only needs some basic vector algebra and doesn't really care about the type of things, apart from the dimension. You could template it on the vertex coordinate type instead, and have it derive the dimensions of the coordinates from that:

template <typename Coordinates>
class NelderMead {
    static constexpr std::size_t nDims = Coordinates::ColsAtCompileTime + 1;
    ...

This still assumes Coordinates is a type that has a ColsAtCompileTime member value, but with a bit of work you could make this more generic. You can now use arrays of float or even std::complex for the vertex coordinates.

Also consider that you might not want to constrain the return type of the error function. Consider making the type of the error function a template parameter as well. You can use concepts to restrict it to only function types that take Coordinates as arguments:

template <typename Coordinates, typename ErrorFunction>
requires std::invocable<ErrorFunction, Coordinates>
class NelderMead {
    static constexpr std::size_t nVertices = Coordinates::ColsAtCompileTime + 1;
    using Error = std::invoke_result_t<ErrorFunction, Coordinates>;
    ...
public:
    NelderMead(ErrorFunction errorFunction, Coordinates initial, ...)
    ...
private:
    ErrorFunction errorFunction;
    std::array<Coordinate, nVertices> values;
    std::array<Error, nVertices> errors;
    ...
};

With this, class template argument deduction (CTAD) works and you no longer have to specify any template arguments when creating an instance of the NelderMead class. It also shows you how to avoid hardcoding double in your code.

Consider keeping values sorted by error

Instead of keeping track of bestId, worstId and secondWorstId, and introducing indirection in a lot of code, it might be better to keep the coordinates of the vertices sorted by their error. This is rather easy to do if you create a new struct to hold both a vertex's coordinates and error, and use that instead of having separate arrays/vectors for values and errors:

struct Vertex {
    Coordinates coordinates;
    ErrorType error;
};

std::array<Vertex, nVertices> vertices;

Sorting it on the error is especially easy with C++20's std::ranges::sort():

std::ranges::sort(vertices, {}, &Vertex::error);

This way, the best vertex is always vertices[0], second best vertices[1], and so on. Basically, you should sort whenever you called invalidateIdsCache(). Sorting is \$O(N \log N)\$, whereas your invalidateIdsCache() is \$O(N)\$, so when your coordinates have a large number of dimensions, it might not be worth it, but at least up to a few tens of dimensions it should be fine.

Potential infinite loop

The Nelder-Mead algorithm doesn't guarantee that the error of the best vertex of the simplex will decrease until zero. Like many other optimization algorithms, it can end up in a local minimum, which might be higher than the minError you provide it. Consider adding a maxSteps parameter to limit the number of steps done, or some other heuristic to detect that the error is not improving significantly anymore.

Answer 4

Include the appropriate headers

main.cpp includes "library.h" but doesn't seem to use anything from it (just as well, as you've not presented it for review). It also includes <iostream> and doesn't use it (though see Don't mix I/O with calculation below).

optimization/NelderMead.h includes <cstring> and doesn't appear to use it. It does however need <string> instead, for std::to_string() - though that is unnecessary, as shown below. It's also missing <cmath>, for std::sqrt().

Don't using namespace in headers

The header file has

using namespace Eigen;

That causes every translation unit that includes the header to have its global namespace filled with Eigen identifiers. That's harmful, and in the worst cases can cause bugs to appear when a newer version of the library is used. So don't import namespaces at global scope.

There's also using namespace Optimization; in the implementation file. That's less harmful, though I still recommend against it. Move the import inside main(), or (better) drop it altogether. After all, it adds more text and cognitive load than it saves!

Choose template type wisely

We have

template<int nDims>
class NelderMead

What happens if nDims is negative? Consider using an unsigned type for the dimension count.

Write initialiser lists in the order of execution

Object members are initialised in the order they are declared in the class. It helps readers if you use the same order in the constructor's initialiser list; this is particularly important when members depend on earlier members or the initialisers have side-effects.

We should be initialising all members (including errors and values); we can initialise members that don't depend on constructor arguments in their definitions.

That gives us:

    std::vector<double> errors{std::vector(nDims + 1, std::numeric_limits<double>::max())};
    int worstValueId = -1;
    int secondWorstValueId = -1;
    int bestValueId = -1;

    NelderMead(ErrorFunction errorFunction, FunctionParameters initial,
               double minError, double initialEdgeLength,
               double shrinkCoeff = 1, double contractionCoeff = 0.5,
               double reflectionCoeff = 1, double expansionCoeff = 1)
        : errorFunction{errorFunction},
          values{initial.replicate(nDims + 1, 1)},
          minError{minError},
          shrinkCoeff{shrinkCoeff},
          expansionCoeff{expansionCoeff},
          contractionCoeff{contractionCoeff},
          reflectionCoeff{reflectionCoeff}
    {

We don't need √2 in every object

We have a constant SQRT2 which doesn't need to be an instance member. Make it a static const and change its name so that it doesn't look like a preprocessor macro (we conventionally use all-caps for macro definitions, to identify hard-to-debug text substitution).

Given that it's only used in the constructor, we can reduce its scope significantly:

    static const auto root2 = std::sqrt(2);
    auto const a = initialEdgeLength / root2;
    auto const b = initialEdgeLength / (nDims * root2) * (std::sqrt(nDims + 1) - 1);

At that point, it's probably better to do away with the named constant entirely, and write std::sqrt(2) in both those places - any decent compiler will use constant values there, and it's just as clear to the reader.

Don't scatter this-> all over the code

In C++, when an object member is referenced, it means the one belonging to this, unless shadowed by a local variable (just don't do that). So cluttered statements such as

            this->errors.at(i) = this->errorFunction(this->values.row(i));

can be much simpler and therefore clearer:

            errors.at(i) = errorFunction(values.row(i));

Know your standard algorithms

There's no need to hand-code this loop:

        for (int i = 0; i < nDims+1; i++)
        {
            this->errors.at(i) = this->errorFunction(this->values.row(i));
        }

Instead, include <algorithm> and <ranges> (C++20) and use the standard-library transform():

        std::ranges::transform(values.rowwise(), errors.begin(), errorFunction);

Don't mix I/O with calculation

NelderMead::optimize() runs the algorithm and prints the result. This makes it inflexible, because calling code can't do anything different with the result. Separate the actions, so that the function either returns all the results (perhaps as an aggregate type) or sets appropriate member variables so that the optimum can be accessed via the object.

Removing printing from this function then removes the need to include <iostream> from the header (and likely transfers it to main.cpp, if we want to print the results from there).

While looking at that output operation, I see you're unnecessarily constructing strings for types that have << operator defined. There's no need to do that:

            std::cout << "Best error " << bestError << " with : " << best()
                      << " Worst error " << worstError << " with : " << worst()
                      << '\n';

Use const appropriately

worst() and best() don't need to change the state of the object, so should be declared const:

    FunctionParameters worst() const { return values.row(worstValueId); }
    FunctionParameters best() const { return values.row(bestValueId); }

I removed inline from these declarations, because functions declared in-class are automatically inline (i.e. can be defined in multiple translation units).

The same applies to a number of other functions.

Algorithms again

Instead of the logic in getMeanWithoutWorst() to sum all but one value, it is simpler to add all the values using std::accumulate (in <numeric>), then subtract the one we want to ignore:

    FunctionParameters getMeanWithoutWorst() const
    {
        auto const& rows = values.rowwise();
        auto const total = std::accumulate(rows.begin(), rows.end(), FunctionParameters{})
            - values.row(worstValueId);

        return total / nDims;
    }

A C++23 standard library should get std::ranges::accumulate() which would simplify further:

    FunctionParameters getMeanWithoutWorst() const
    {
        auto const total = std::ranges::accumulate(values.rowwise(), FunctionParameters{})
            - values.row(worstValueId);

        return total / nDims;
    }

Spelling

It sounds trivial, but spelling words correctly can really help when searching codebases!

searchedParamters should be searchedParameters.