El algoritmo Nelder-Mead, también conocido como el método símplex de optimización, es una técnica heurística de búsqueda directa que no requiere el cálculo de derivadas. Es particularmente robusto para optimizar funciones no lineales donde el gradiente es desconocido o ruidoso. En este artículo, presentamos una implementación en C++ utilizando la biblioteca Eigen, siguienod la lógica de la función scipy.optimize.fmin de SciPy.

Características de la implementación inspirada en SciPy

A diferencia de la versión estándar del algoritmo, la implementación de SciPy (y esta versión en C++) incluye varias mejoras críticas para el uso en producción:

• Restricciones de caja (Bounds): Soporte para limitar el espacio de búsqueda mediante el recorte (clipping) de los vértices del símplex.
• Parámetros adaptativos: Ajuste automático de los coeficientes de reflexión, expansión y contracción según la dimensión del problema (basado en la propuesta de Gao y Han, 2012).
• Criterio de convergencia dual: El algoritmo solo se detiene cuando tanto la diferencia en los valores de la función como la distancia entre los puntos del símplex son menores a tolerancias específicas (fatol y xatol).

Implementación Técnica

Configuración del Optimizador (SolverConfig.hpp)

Definimos las estructuras necesarias para gestionar las opciones de búsqueda y los resultados obtenidos.

#ifndef SOLVER_CONFIG_HPP
#define SOLVER_CONFIG_HPP

#include <Eigen/Dense>
#include <vector>
#include <string>
#include <functional>

struct NelderMeadOptions {
    double tol_x = 1e-4;       // Tolerancia absoluta para x
    double tol_f = 1e-4;       // Tolerancia absoluta para f(x)
    int max_iter = -1;         // Límite de iteraciones
    int max_eval = -1;         // Límite de evaluaciones de función
    bool use_adaptive = true;  // Parámetros basados en la dimensión
    bool verbose = false;      // Mostrar progreso en consola
    
    // Pares (min, max) para cada variable
    std::vector<std::pair<double, double>> constraints;
};

struct OptimizationResult {
    Eigen::VectorXd best_x;
    double best_f;
    int iterations = 0;
    int eval_count = 0;
    bool success = false;
    std::string status_msg;
};

#endif

El Núcleo del Algoritmo (NelderMeadSolver.ipp)

La lógica principal utiliza operaciones vectorizadas de Eigen para calcular el centroide y las transformaciones del símplex.

#include <algorithm>
#include <iostream>

template <typename T = double>
OptimizationResult run_nelder_mead(
    std::function<T(const Eigen::Matrix<T, Eigen::Dynamic, 1>&)> objective,
    const Eigen::Matrix<T, Eigen::Dynamic, 1>& start_point,
    const NelderMeadOptions& opts = NelderMeadOptions()) 
{
    using Vector = Eigen::Matrix<T, Eigen::Dynamic, 1>;
    using Matrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;

    const int dim = start_point.size();
    
    // Parámetros de transformación
    double rho, chi, psi, sigma;
    if (opts.use_adaptive) {
        rho = 1.0;
        chi = 1.0 + 2.0 / dim;
        psi = 0.75 - 1.0 / (2.0 * dim);
        sigma = 1.0 - 1.0 / dim;
    } else {
        rho = 1.0; chi = 2.0; psi = 0.5; sigma = 0.5;
    }

    // Inicialización del símplex
    Matrix simplex(dim + 1, dim);
    Vector scores(dim + 1);
    
    auto apply_bounds = [&](const Vector& v) -> Vector {
        if (opts.constraints.empty()) return v;
        Vector limited = v;
        for (int i = 0; i < dim; ++i) {
            limited(i) = std::max(opts.constraints[i].first, 
                         std::min(opts.constraints[i].second, v(i)));
        }
        return limited;
    };

    simplex.row(0) = apply_bounds(start_point);
    for (int i = 0; i < dim; ++i) {
        Vector next_point = start_point;
        double step = (std::abs(next_point(i)) < 1e-6) ? 0.00025 : 0.05;
        next_point(i) += step * next_point(i);
        simplex.row(i + 1) = apply_bounds(next_point);
    }

    int fevals = 0;
    for (int i = 0; i <= dim; ++i) {
        scores(i) = objective(simplex.row(i).transpose());
        fevals++;
    }

    auto sort_simplex = [&]() {
        std::vector<int> indices(dim + 1);
        std::iota(indices.begin(), indices.end(), 0);
        std::sort(indices.begin(), indices.end(), [&](int a, int b) {
            return scores(a) < scores(b);
        });

        Matrix old_simplex = simplex;
        Vector old_scores = scores;
        for (int i = 0; i <= dim; ++i) {
            simplex.row(i) = old_simplex.row(indices[i]);
            scores(i) = old_scores(indices[i]);
        }
    };

    sort_simplex();
    int iterations = 0;
    int max_it = (opts.max_iter > 0) ? opts.max_iter : dim * 200;
    int max_fe = (opts.max_eval > 0) ? opts.max_eval : dim * 200;

    while (iterations < max_it && fevals < max_fe) {
        // Evaluación de convergencia
        double diff_f = (scores.array() - scores(0)).abs().maxCoeff();
        double diff_x = 0;
        for(int i=1; i<=dim; ++i) {
            diff_x = std::max(diff_x, (simplex.row(i) - simplex.row(0)).norm());
        }

        if (diff_f <= opts.tol_f && diff_x <= opts.tol_x) break;

        // Cálculo del centroide (excluyendo el peor punto)
        Vector centroid = simplex.topRows(dim).colwise().mean().transpose();
        
        // Reflexión
        Vector xr = apply_bounds((1 + rho) * centroid - rho * simplex.row(dim).transpose());
        double fxr = objective(xr); fevals++;

        if (fxr < scores(0)) {
            // Expansión
            Vector xe = apply_bounds((1 + rho * chi) * centroid - rho * chi * simplex.row(dim).transpose());
            double fxe = objective(xe); fevals++;
            if (fxe < fxr) {
                simplex.row(dim) = xe.transpose(); scores(dim) = fxe;
            } else {
                simplex.row(dim) = xr.transpose(); scores(dim) = fxr;
            }
        } else if (fxr < scores(dim - 1)) {
            simplex.row(dim) = xr.transpose(); scores(dim) = fxr;
        } else {
            // Contracción
            bool shrink = false;
            if (fxr < scores(dim)) {
                // Contracción externa
                Vector xc = apply_bounds((1 + psi * rho) * centroid - psi * rho * simplex.row(dim).transpose());
                double fxc = objective(xc); fevals++;
                if (fxc <= fxr) {
                    simplex.row(dim) = xc.transpose(); scores(dim) = fxc;
                } else shrink = true;
            } else {
                // Contracción interna
                Vector xcc = apply_bounds((1 - psi) * centroid + psi * simplex.row(dim).transpose());
                double fxcc = objective(xcc); fevals++;
                if (fxcc < scores(dim)) {
                    simplex.row(dim) = xcc.transpose(); scores(dim) = fxcc;
                } else shrink = true;
            }

            if (shrink) {
                for (int i = 1; i <= dim; ++i) {
                    simplex.row(i) = apply_bounds(simplex.row(0) + sigma * (simplex.row(i) - simplex.row(0)));
                    scores(i) = objective(simplex.row(i).transpose()); fevals++;
                }
            }
        }
        
        iterations++;
        sort_simplex();
    }

    OptimizationResult res;
    res.best_x = simplex.row(0);
    res.best_f = scores(0);
    res.iterations = iterations;
    res.eval_count = fevals;
    res.success = (iterations < max_it && fevals < max_fe);
    res.status_msg = res.success ? "Convergencia alcanzada." : "Límite excedido.";
    return res;
}

Ejemplo de uso: Función de Rosenbrock

La función de Rosenbrock es un benchmark clásico para algoritmos de optimización debido a su valle estrecho y curvado.

#include "SolverConfig.hpp"
#include "NelderMeadSolver.ipp"

int main() {
    auto rosenbrock = [](const Eigen::VectorXd& v) -> double {
        double x = v(0), y = v(1);
        return std::pow(1 - x, 2) + 100 * std::pow(y - x * x, 2);
    };

    Eigen::VectorXd initial_guess(2);
    initial_guess << -1.0, 2.0;

    NelderMeadOptions config;
    config.verbose = true;

    auto result = run_nelder_mead<double>(rosenbrock, initial_guess, config);

    std::cout << "Resultado: " << result.best_x.transpose() << std::endl;
    std::cout << "Valor f(x): " << result.best_f << std::endl;
    std::cout << "Estado: " << result.status_msg << std::endl;

    return 0;
}

Análisis de Convergencia

La robustez de esta versión radica en el manejo de los parámetros adaptativos. En dimensiones bajas (N=2), los coeficientes clásicos son suficientes, pero a medida que el número de variables aumenta, el símplex tiende a colapsar o degenerar. El esquema de Gao-Han ralentiza la contracción y acelera la expansión en espacios de alta dimensionalidad, permitiendo que el algoritmo explore de manera más eficiente antes de converger hacia el mínimo local.

Además, el uso de Eigen permite que las operaciones de álgebra lineal necesarias para el cálculo del centroide y las actualizcaiones de posición se ejecuten con una eficiencia óptima, aprovechando instrucciones vectoriales de la CPU (como AVX o SSE) si están disponibles.