Русский Português
preview
Estrategia de Evolución de Adaptación de la Matriz de Covarianza — Covariance Matrix Adaptation Evolution Strategy (CMA-ES)

Estrategia de Evolución de Adaptación de la Matriz de Covarianza — Covariance Matrix Adaptation Evolution Strategy (CMA-ES)

MetaTrader 5Trading |
20 0
Andrey Dik
Andrey Dik

Contenido

  1. Introducción
  2. Implementación del algoritmo
  3. Resultados de las pruebas


Introducción

En el mundo de la optimización existen muchos algoritmos, pero buscamos los más potentes para resolver los problemas de optimización de nuestros robots comerciales. El CMA-ES (Covariance Matrix Adaptation Evolution Strategy) es uno de esos raros ejemplos en los que el rigor matemático se combina con la intuición biológica, creando un algoritmo que no solo resuelve problemas de optimización, sino que también aprende a comprender su estructura.

La historia del CMA-ES comenzó a finales de la década de 1990 en laboratorios de investigación en Alemania, donde Nikolaus Hansen y Andreas Ostermeier plantearon una pregunta fundamental: ¿es posible crear un algoritmo de optimización que no solo busque una solución, sino que se adapte a la geometría del problema? Los algoritmos evolutivos tradicionales generaban descendencia en regiones esféricas alrededor de los individuos progenitores, lo cual funcionaba bien para funciones simples, pero resultaba ineficaz para problemas complejos y mal condicionados. Así pues, vamos a echar un vistazo a este interesante algoritmo.


Implementación del algoritmo

Imagine buscar un tesoro en una isla de forma irregular. El método habitual consiste en buscar en todas las direcciones por igual, como si la isla fuera redonda. El CMA-ES estudia gradualmente la forma de la isla y comienza a buscar principalmente en aquellas direcciones donde la probabilidad de encontrar tesoros es mayor. Además, recuerda las rutas exitosas y usa esta memoria para planificar búsquedas futuras.

El CMA-ES se basa en una ecuación engañosamente simple: x_k ~ N(m, σ²C), Pero tras esta simplicidad se esconde una profunda estructura matemática. Cada símbolo aquí contiene información importante sobre el estado de la búsqueda: m es la mejor estimación actual sobre la ubicación del óptimo, σ es una medida de cuánto estamos dispuestos a alejarnos de los conocido, y C es una matriz de covarianza que codifica nuestra comprensión de la geometría de la función. El único cambio que podemos justificar es reemplazar la distribución normal por una distribución de potencia, lo cual significa que la implementación seguirá una fórmula modificada: x_k ~ PowerDist(m, σ²C). Esta modificación cambia la naturaleza de la exploración espacial (saltos más amplios), pero conserva la naturaleza adaptativa fundamental del algoritmo.

La matriz de covarianza C es el verdadero corazón del algoritmo. Esta comienza su vida como una sencilla matriz identidad que representa una distribución esférica, pero con cada iteración evoluciona, expandiéndose en áreas de rápida mejora y reduciéndose donde el progreso es lento. Gradualmente, la esfera se transforma en una elipse, y luego en un elipsoide alargado, idealmente orientado a lo largo de los contornos de la función que se está optimizando.

La principal innovación de CMA-ES es el concepto de trayectorias evolutivas. Se trata de una especie de memoria genética del algoritmo que recuerda no solo dónde se produjeron los aciertos, sino también cómo llegó el algoritmo a ellos. La trayectoria evolutiva acumula información sobre los sucesivos pasos exitosos, creando un vector direccional que apunta a las áreas de búsqueda más prometedoras. La segunda vía evolutiva cumple una función más sutil, a saber, controla el tamaño del paso. Si los pasos sucesivos del algoritmo están correlacionados, es decir, cada paso subsiguiente continúa la dirección del anterior, entonces el tamaño del paso aumenta: el algoritmo "siente" que se está moviendo en la dirección correcta. Si los pasos son aleatorios y no están correlacionados, el tamaño del paso disminuye; quizá el algoritmo ya esté cerca del óptimo y necesite una búsqueda más cuidadosa.

Detrás de la metáfora biológica de la evolución en CMA-ES se esconde un principio matemático riguroso: la maximización de la verosimilitud. El algoritmo se pregunta constantemente: ¿qué parámetros de la distribución hacen que los puntos de éxito observados sean los más probables? Esto transforma la optimización evolutiva de un método heurístico en un método basado en estadísticas.

Cada actualización de la matriz de covarianza C consta de dos componentes: la actualización rank-one basada en la trayectoria evolutiva y la actualización rank-μ que utiliza información de toda la población muestreada. La actualización rank-one proporciona estabilidad y la capacidad de capturar tendencias a largo plazo, mientras que la actualización rank-μ permite una rápida adaptación a la nueva información de la generación actual.

En el contexto del CMA-ES, se utiliza una función de Heaviside modificada que desempeña un papel importante en el mecanismo de detección de estancamiento del algoritmo. La función compara la longitud de la trayectoria evolutiva con la longitud esperada; si la trayectoria es demasiado corta, es señal de "estancamiento". Condiciones de desactivación (hsig = 0): el algoritmo "vaga" aleatoriamente y los pasos se cancelan entre sí, entonces el tamaño del paso es probablemente demasiado grande. Durante el estancamiento, la influencia de las actualizaciones rank-one disminuye y dependemos más de la información de toda la población. Condiciones de activación (hsig = 1): el algoritmo avanza en una dirección determinada, los pasos sucesivos están correlacionados y, por lo tanto, el tamaño del paso es adecuado a la situación actual.

Una de las propiedades más importantes del CMA-ES es su invariancia ante las transformaciones del espacio de búsqueda. El algoritmo resuelve con igual eficiencia la función f(x) y la función f(Ax + b), donde A es cualquier matriz no singular, mientras que b es cualquier vector de desplazamiento. Esto significa que CMA-ES es independiente de la elección del sistema de coordenadas. Esta invariancia no es accidental; supone una consecuencia directa del uso del principio de máxima verosimilitud y de la adaptación de la matriz de covarianza. El algoritmo detecta automáticamente un sistema de coordenadas natural para el problema donde los ejes coinciden con las direcciones principales de variación de la función.

La belleza de la teoría debe combinarse con la aplicabilidad práctica. El CMA-ES requiere una memoria de O(n²) para almacenar la matriz de covarianza y un cálculo de O(n³) para su propia descomposición, lo cual hace que el algoritmo sea aplicable a problemas con dimensionalidades de hasta varios cientos de variables. Para dimensionalidades grandes, se han desarrollado modificaciones especializadas: sep-CMA-ES se limita a matrices de covarianza diagonales, VkD-CMA utiliza dimensiones variables y LM-CMA aplica los principios de memoria limitada. La implementación de estos métodos actualmente escapa al alcance de este artículo, pero es posible que se retomen en artículos posteriores.

cmaes-algorithm

Figura 1. Ilustración del algoritmo CMA-ES

La ilustración muestra la evolución a través de las generaciones: tres instantáneas del estado del algoritmo (generaciones 1, 5 y 15) que muestran cómo la población converge gradualmente al óptimo.

Adaptación de la matriz de covarianza: las elipsis azules muestran cómo cambia la forma de la distribución conforme avanza la generación: 1 - redonda (matriz identidad), 5 - alargada y rotada, 15 - ajustada con precisión a la geometría local.
Componentes del algoritmo: los puntos rojos representan la población total (λ descendientes), los puntos verdes son las mejores soluciones (μ padres), el punto negro m es la media de la población y los círculos punteados son las isolíneas de la función objetivo.

La fórmula clave se encuentra a continuación, con una nota sobre la modificación de la ley de potencias. La ilustración demuestra claramente cómo el algoritmo adapta su estrategia de búsqueda al entorno de la función optimizada, reduciendo gradualmente el espacio de búsqueda y acercándose al óptimo.

Pasemos ahora a escribir el pseudocódigo del algoritmo.

Inicialización

  1. Establecemos los parámetros del algoritmo:
    • Tamaño de la población (lambda) = 50
    • Número de padres (mu) = 25
    • Tasa de aprendizaje para la actualización de rango 1 (c1) = 0,01
    • Tasa de aprendizaje para la actualización de rango-μ (cμ) = 0,8
    • Amortiguación del tamaño del paso = 0,6
    • Tamaño del paso inicial (sigma) = 0,3
  2. Calcular los pesos de recombinación:
    • Para cada uno de los μ padres, calculamos el peso como log(μ + 0.5) - log(i + 1).
    • Normalizamos los pesos de modo que su suma sea igual a 1.
    • Calculamos la masa efectiva de elección μ_eff
  3. Inicializamos los parámetros de la estrategia:
    • Calculamos las tasas de aprendizaje para las trayectorias evolutivas (cs, cc)
    • Determinamos el valor esperado de la distribución normal estándar (chiN).
    • Determinamos el intervalo para la descomposición en valores propios.
  4. Inicializamos las estructuras de datos:
    • Matriz de covarianza C = matriz de identidad
    • Matriz de autovectores B = matriz de identidad
    • Vector de valores propios D = vector unitario
    • Trayectorias de evolución pc y ps = vectores cero
    • Media poblacional m = punto aleatorio en el espacio de búsqueda
  5. Creamos una población inicial:
    • Generamos puntos aleatorios lambda en el espacio de búsqueda

El ciclo básico de la evolución

Repetimos el procedimiento hasta que se cumpla el criterio de parada:

Etapa 1: Generación de descendientes 

  1. Para cada uno de los descendientes de lambda:
    • Generamos un vector aleatorio z a partir de la distribución normal estándar.
    • Aplicamos la transformación: y = B × D × z
    • Creamos un descendiente: x_k = m + σ × y
    • Aplicamos las restricciones de espacio de búsqueda a x_k

Etapa 2: Evaluación y actualización 

  1. Evaluación de la aptitud de todos los descendientes.
  2. Clasificación de la población:
    • Clasificamos los descendientes en orden descendente de aptitud
    • Actualizamos la mejor solución encontrada si es necesario.
  3. Actualización de la media poblacional:
    • Guardamos la media antigua m_old
    • Calculamos la nueva media como la suma ponderada de los μ mejores descendientes: m_new = Σ(w_i × x_i) para i desde 1 hasta μ
  4. Actualización de las rutas de evolución:
    • Calculamos el desplazamiento medio: y_w = (m_new - m_old) / σ
    • Calculamos C^(-1/2) × y_w mediante descomposición en valores propios.
    • Actualizamos la ruta de evolución para el tamaño del paso: ps = (1 - cs) × ps + √(cs × (2 - cs) × μ_eff) × C^(-1/2) × y_w
    • Comprobamos la condición de estancamiento (función de Heaviside).
    • Actualizamos la trayectoria de evolución para la matriz de covarianza: pc = (1 - cc) × pc + indicador_de_estancamiento × √(cc × (2 - cc) × μ_eff) × y_w
  5. Actualización de la matriz de covarianza:
    • Preparamos la actualización de rango 1 a partir del producto exterior de pc
    • Preparamos una actualización de rango-μ a partir de la suma ponderada de las desviaciones de los mejores descendientes.
    • Actualizamos C: C = (1 - c1 - cμ) × C + c1 × (pc × pc^T) + cμ × Σ(w_i × y_i × y_i^T)
    • Aseguramos la simetría de la matriz
  6. Actualización del tamaño del paso:
    • Calculamos el coeficiente de adaptación según la longitud de la ruta ps.
    • Actualizamos de σ: σ = σ × exp((cs/damps) × (||ps||/chiN - 1))
    • Limitamos σ dentro de límites razonables para la estabilidad numérica.
  7. Descomposición en valores propios (si es necesario):
    • Si han transcurrido suficientes iteraciones desde la última descomposición:
      • Realizamos la expansión de Jacobi de C = B × D² × B^T
      • Clasificamos los valores propios y los vectores propios en orden descendente.
      • Aseguramos la positividad definida de la matriz
  8. Comprobación y corrección de la matriz de covarianza:
    • Comprobamos periódicamente la positividad definida de C.
    • Si es necesario, añadimos un pequeño valor a la diagonal.
    • Aseguramos la simetría de la matriz

Ahora vamos a escribir una clase " C_AO_CMAES", derivada de la clase " C_AO", que implementará el algoritmo de optimización CMA-ES.

Campos públicos, método:

  • SetParams () - establece los valores de los parámetros CMA-ES del array "params".
  • Init () - inicializa el algoritmo. Admite valores mínimos y máximos, así como el tamaño del paso para cada variable y el número de épocas.
  • Moving () - implementa el ciclo principal del algoritmo, que es el responsable de generar nuevas soluciones.
  • Revision () - evalúa nuevas soluciones y actualiza el estado del algoritmo.

Campos privados:

  • Parámetros CMA-ES - variables que contiene parámetros específicos del algoritmo CMA-ES: número de padres (mu), tamaño del paso (sigma), tasas de aprendizaje (learningRateC1, learningRateCMu), factor de amortiguación (stepSizeDamping) y otros necesarios para que el algoritmo funcione.
  • Estructuras de datos CMA-ES - arrays utilizados para almacenar los pesos de recombinación (weights), la matriz de covarianza (covMatrix), los vectores propios (B), los valores propios (D), las rutas de evolución (pc, ps), así como otros datos auxiliares. Las variables  "mu_eff", "counteval", "eigeneval", "chiN", "eigenInterval" se utilizan para controlar y gestionar el progreso del algoritmo.
  • Variables de optimización del rendimiento: están diseñadas para acelerar los cálculos durante la ejecución del algoritmo, como las tasas de aprendizaje y la amortiguación.
  • Arrays auxiliaresse utilizan como almacenamiento temporal para vectores y matrices al realizar diversas operaciones.
  • Métodos auxiliares: realizan pasos individuales del algoritmo CMA-ES, como la inicialización de la distribución, la actualización de la distribución, el cálculo de la factorización de valores propios, la clasificación de la población, el cálculo de los pesos, la actualización de la media, la comprobación de la positividad definida y la garantía de la positividad definida.

En general, la clase " C_AO_CMAES" encapsula la lógica del algoritmo CMA-ES, ofreciendo una interfaz para la inicialización, el ajuste de parámetros y la realización de la optimización. Contiene los parámetros, los arrays de datos y los métodos necesarios para implementar el algoritmo. La separación de campos públicos y privados ofrece control de acceso a los componentes internos de la clase.

//————————————————————————————————————————————————————————————————————
class C_AO_CMAES : public C_AO
{
  public: //----------------------------------------------------------
  ~C_AO_CMAES () { }
  C_AO_CMAES ()
  {
    ao_name = "CMAES";
    ao_desc = "Covariance Matrix Adaptation Evolution Strategy";
    ao_link = "https://www.mql5.com/ru/articles/18227";

    // Параметры по умолчанию
    popSize         = 50;          // Размер популяции по умолчанию (lambda)
    mu              = 25;          // Количество родителей (половина популяции)
    learningRateC1  = 0.01;        // Скорость обучения для ранг-1 обновления
    learningRateCMu = 0.8;         // Скорость обучения для ранг-μ обновления
    stepSizeDamping = 0.6;         // Демпфирование для размера шага

    // Создание и инициализация массива параметров
    ArrayResize (params, 5);
    params [0].name = "popSize";         params [0].val = popSize;
    params [1].name = "mu";              params [1].val = mu;
    params [2].name = "learningRateC1";  params [2].val = learningRateC1;
    params [3].name = "learningRateCMu"; params [3].val = learningRateCMu;
    params [4].name = "stepSizeDamping"; params [4].val = stepSizeDamping;
  }

  void SetParams ()
  {
    popSize         = (int)params [0].val;
    mu              = (int)params [1].val;
    learningRateC1  = params      [2].val;
    learningRateCMu = params      [3].val;
    stepSizeDamping = params      [4].val;
  }

  bool Init (const double &rangeMinP  [],  // минимальные значения
             const double &rangeMaxP  [],  // максимальные значения
             const double &rangeStepP [],  // размер шага
             const int     epochsP = 0);

  void Moving   ();
  void Revision ();

  //------------------------------------------------------------------
  private: //---------------------------------------------------------
  // Специфичные параметры CMA-ES
  int mu;                  // Количество родителей (выбранных точек)
  double sigma;            // Размер шага
  double learningRateC1;   // Скорость обучения для ранг-1 обновления
  double learningRateCMu;  // Скорость обучения для ранг-μ обновления
  double stepSizeDamping;  // Фактор демпфирования для обновления размера шага

  // Специфичные структуры данных CMA-ES
  double weights   [];      // Веса рекомбинации
  double covMatrix [];      // Ковариационная матрица (хранится как одномерный массив)
  double B  [];             // Собственные векторы C
  double D  [];             // Собственные значения (квадратные корни) C
  double pc [];             // Путь эволюции для C
  double ps [];             // Сопряженный путь эволюции для sigma
  double mu_eff;            // Эффективная масса выбора дисперсии
  int    counteval;         // Счетчик вычислений функции с последнего разложения
  int    eigeneval;         // Счетчик генераций, когда выполнялось разложение
  double chiN;              // Ожидаемая норма N(0,I)
  int    eigenInterval;     // Интервал для разложения на собственные значения

  // Переменные для оптимизации производительности
  double cs;                // Скорость обучения для пути sigma
  double cc;                // Скорость обучения для пути ранг-1
  double damps;             // Демпфирование для sigma
  double hsig_threshold;    // Порог для функции Хевисайда

  // Вспомогательные массивы
  double y_vec     [];         // Вектор мутации
  double arindex   [];         // Массив индексов для сортировки
  double arfitness [];         // Массив фитнеса для сортировки
  double temp_vec  [];         // Временный вектор для матричных операций
  double invsqrtC_times_yw []; // Временное хранилище для C^(-1/2) * y_w

  // Переменные кэширования для Бокса-Мюллера
  double cached_normal;
  bool   has_cached;

  // Вспомогательные методы
  void   InitDistribution          ();
  void   UpdateDistribution        ();
  void   ComputeEigendecomposition ();
  double GetChiN                   ();
  void   SortPopulation            ();
  void   ComputeWeights            ();
  void   UpdateMean                ();
  bool   CheckPositiveDefinite     ();
  void   EnforcePositiveDefinite   ();
};
//————————————————————————————————————————————————————————————————————

El método "Init" de la clase "C_AO_CMAES" inicializa el algoritmo CMA-ES. Recibe como entrada los valores mínimo y máximo de los parámetros, el tamaño del paso y el número de épocas. En primer lugar, llamamos al método "StandardInit" para realizar la inicialización estándar. El indicador "has_cached" se inicializa en "false" para indicar que no hay ningún valor guardado en la caché para la generación de números aleatorios de Box-Muller. A continuación, el tamaño del paso "sigma" se inicializa con un valor inicial de 0,3 (30% del rango de búsqueda). El método ComputeWeights se llama para calcular los pesos de recombinación, mientras que el método GetChiN se llama para calcular la norma esperada N(0, I).

Los parámetros "cs", "cc", "damps" y "hsig_threshold" se calculan según "mu_eff" (masa efectiva) y el número de coordenadas (coords). Estos parámetros se usan para adaptar la estrategia CMA-ES durante el funcionamiento. El parámetro "eigenInterval" se calcula y se establece para determinar con qué frecuencia se realiza la descomposición de valores propios. Esta configuración está diseñada para optimizar el rendimiento.

Se asigna memoria para los arrays específicos de CMA-ES, como la matriz de covarianza "covMatrix", los vectores propios "B", los valores propios "D", las rutas de evolución "pc" y "ps", los vectores temporales "y_vec", "arindex", "arfitness", "temp_vec" e "invsqrtC_times_yw". Los arrays se recrean únicamente si es necesario cambiar la dimensión del problema.

Los arrays de trayectoria de evolución "pc" y "ps" se inicializan a cero, la matriz de covarianza "covMatrix" y la matriz de vectores propios "B" se inicializan a la matriz identidad, y el array de valores propios "D" se inicializa a unidades. El método "InitDistribution" se llama para inicializar la distribución inicial. A continuación, se reinician los contadores de evaluación de las funciones "counteval" y "eigeneval". El indicador "revision" se establece en "true" para indicar que los valores de aptitud deben recalcularse. Se retornará "true" si la inicialización es exitosa.

Finalmente, el método " Init" realiza toda la inicialización necesaria para que funcione el algoritmo CMA-ES, incluyendo la configuración de parámetros, la asignación de memoria para los arrays y la inicialización de los valores iniciales. 

//————————————————————————————————————————————————————————————————————
bool C_AO_CMAES::Init (const double &rangeMinP [],  // минимальные значения
                       const double &rangeMaxP [],  // максимальные значения
                       const double &rangeStepP [], // размер шага
                       const int     epochsP = 0)   // количество эпох
{
  if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false;

  //------------------------------------------------------------------
  // Инициализация кэширования Бокса-Мюллера
  has_cached = false;

  // Инициализация специфичных переменных CMA-ES
  sigma          = 0.3;  // Начальный размер шага (30% от диапазона поиска)

  // Вычисление эффективной массы выбора дисперсии
  ComputeWeights ();

  // Ожидаемая норма N(0,I)
  chiN = GetChiN ();

  // Вычисление и сохранение параметров стратегии
  cs = (mu_eff + 2.0) / (coords + mu_eff + 5.0);
  cc = (4.0 + mu_eff / coords) / (coords + 4.0 + 2.0 * mu_eff / coords);
  damps = 1.0 + 2.0 * MathMax (0.0, MathSqrt ((mu_eff - 1.0) / (coords + 1.0)) - 1.0) + cs;
  hsig_threshold = 1.4 + 2.0 / (coords + 1.0);

  // Установка интервала разложения на собственные значения - настройка для производительности
  eigenInterval = (int)(coords / (10.0 * MathSqrt (learningRateC1 + learningRateCMu)));
  eigenInterval = MathMax (1, eigenInterval);

  // Выделение массивов - выделяем только один раз
  ArrayResize (covMatrix, coords * coords);
  ArrayResize (B, coords * coords);
  ArrayResize (D, coords);
  ArrayResize (pc, coords);
  ArrayResize (ps, coords);
  ArrayResize (y_vec, coords);
  ArrayResize (arindex, popSize);
  ArrayResize (arfitness, popSize);
  ArrayResize (temp_vec, coords);
  ArrayResize (invsqrtC_times_yw, coords);

  // Инициализация путей эволюции нулями
  ArrayInitialize (pc, 0);
  ArrayInitialize (ps, 0);

  // Быстрая инициализация единичной ковариационной матрицы и разложения
  ArrayInitialize (covMatrix, 0.0);
  ArrayInitialize (B, 0.0);

  for (int i = 0; i < coords; i++)
  {
    covMatrix [i * coords + i] = 1.0;
    B [i * coords + i] = 1.0;
    D [i] = 1.0;
  }

  // Инициализация начального распределения
  InitDistribution ();

  // Сброс счетчиков вычислений
  counteval = 0;
  eigeneval = 0;

  // Принудительный пересчет фитнеса
  revision = true;

  return true;
}
//————————————————————————————————————————————————————————————————————

El método "Moving" de la clase "C_AO_CMAES" es responsable de generar nuevos descendientes (soluciones) en la población. Transformación y = B * D * z: dentro del ciclo sobre la población hay un ciclo anidado que itera sobre todas las "coordenadas" de cada individuo. En cada coordenada "i" se calcula el valor "y_vec[i]". Esto se realiza usando la multiplicación de matrices. De hecho, supone la multiplicación de un vector de números aleatorios "z" por una matriz "B * D", donde "B" son los vectores propios de la matriz de covarianza, "D" son los valores propios de la matriz de covarianza y "z" son números aleatorios generados mediante la función "u.PowerDistribution". 

Generación de descendientes: para cada coordenada "i" de cada descendiente "k", se calcula el valor (a[k].c[i]). Esto se hace sumando a la media actual "cB[i]" (el centro de la distribución actual) el vector escalado "y_vec[i]". El escalado se logra multiplicando "y_vec[i]" por el tamaño del paso "sigma". Por consiguiente, a[k].c[i] = cB[i] + sigma * y_vec[i].

La función "u.SeInDiSp" garantiza que los valores de coordenadas de los descendientes permanezcan dentro de los rangos válidos definidos por "rangeMin[i]", "rangeMax[i]" y "rangeStep[i]". Esta puede recortar valores, hacerlos rebotar en los límites y cuantificarlos al valor de paso válido más cercano. Una vez generados todos los descendientes, el indicador "revision" se establece en "true" y se deben recalcular los valores de la función de aptitud (calidad) para los descendientes generados.

Como resultado, el método genera una nueva población de soluciones basada en la distribución actual (media " cB" y la matriz de covarianza representada por " B". y " D") y el tamaño del paso " sigma". Este usa números aleatorios para crear pequeñas variaciones alrededor de la media actual y así explorar el espacio de búsqueda. La función " SeInDiSp" garantiza que las nuevas decisiones se mantengan dentro de límites aceptables y la bandera de "revision" garantiza que se evaluará la calidad de estas decisiones.

//————————————————————————————————————————————————————————————————————
void C_AO_CMAES::Moving ()
{
  // Генерация lambda новых потомков
  for (int k = 0; k < popSize; k++)
  {

    // Применение преобразования y = B*D*z
    for (int i = 0; i < coords; i++)
    {
      y_vec [i] = 0.0;
      for (int j = 0; j < coords; j++)
      {
        y_vec [i] += B [i * coords + j] * D [j] * u.PowerDistribution (0.0, -8, 8, 20);
      }

      // Генерация потомка: x_k = m + σ * y
      a [k].c [i] = cB [i] + sigma * y_vec [i];
      a [k].c [i] = u.SeInDiSp (a [k].c [i], rangeMin [i], rangeMax [i], rangeStep [i]);
    }
  }

  // Отметка для пересчета фитнеса
  revision = true;
}
//————————————————————————————————————————————————————————————————————

El método "revision" de la clase " C_AO_CMAES" es responsable de actualizar los parámetros del algoritmo después de evaluar la función de aptitud para una nueva población. En primer lugar, el método comprueba el valor del indicador "revision".  A continuación, se llama al método "SortPopulation()", que clasifica la población actual de manera que los mejores individuos se coloquen al principio. El método "UpdateDistribution()" actualiza los parámetros de la distribución CMA-ES (la media "cB" del vector de estrategia, la matriz de covarianza, el tamaño del paso "sigma" y las rutas evolutivas) basándose en información sobre los mejores individuos de la población. La actualización de la distribución supone un paso clave de CMA-ES, ya que permite que el algoritmo adapte su estrategia de búsqueda en función de los resultados obtenidos. 

Como resultado, el método " Revision" representa el ciclo principal de adaptación de CMA-ES. Este clasifica la población, actualiza los parámetros de distribución en función de los mejores individuos e incrementa el contador de cálculos. 

//————————————————————————————————————————————————————————————————————
void C_AO_CMAES::Revision ()
{
  if (!revision) return;

  revision = false;

  // Сортировка популяции по фитнесу
  SortPopulation ();

  // Обновление параметров распределения на основе выбранных особей
  UpdateDistribution ();

  // Обновление счетчика вычислений
  counteval++;
}
//————————————————————————————————————————————————————————————————————

El método "InitDistribution" de la clase "C_AO_CMAES" se encarga de inicializar la distribución inicial de la búsqueda de soluciones. El método determina un valor medio inicial (cB) para cada coordenada en el espacio de búsqueda. Para cada coordenada "i", el valor de "cB[i]" se establece aleatoriamente en el rango entre "rangeMin[i]" y "rangeMax[i]" usando la función "u.RNDfromCI()", con el valor medio inicial colocado en el centro del rango válido para cada variable. El ciclo exterior itera sobre toda la población (popSize individuos). Dentro de este ciclo hay un ciclo interno que itera sobre todas las coordenadas de cada individuo. Para cada descendiente, se generan sus coordenadas. Primero se eligen aleatoriamente del rango dado "RNDfromCI" y luego se redondean al paso de muestreo usando "SeInDiSp".

La función " u.RNDfromCI" genera un número aleatorio con distribución uniforme en un intervalo dado. La función " u.SeInDiSp" garantiza que las coordenadas generadas para cada individuo estén dentro de los límites aceptables definidos rangoMin[i], rangoMax[i] y rangeStep[i].

El método inicializa la población inicial con soluciones aleatorias distribuidas uniformemente en el espacio de búsqueda. El centro de distribución (cB) también se selecciona aleatoriamente. Esto ofrece un punto de partida para las iteraciones posteriores de CMA-ES, en las que el algoritmo adaptará la estrategia de búsqueda para encontrar la solución óptima.

//+------------------------------------------------------------------+
//| Инициализация распределения поиска                               |
//+------------------------------------------------------------------+
void C_AO_CMAES::InitDistribution ()
{
  // Установка начального среднего в центр пространства поиска
  for (int i = 0; i < coords; i++)
  {
    cB [i] = u.RNDfromCI (rangeMin [i], rangeMax [i]);
  }

  for (int k = 0; k < popSize; k++)
  {
    for (int i = 0; i < coords; i++)
    {
      // Генерация равномерно распределенной точки
      a [k].c [i] = u.RNDfromCI (rangeMin [i], rangeMax [i]);
      a [k].c [i] = u.SeInDiSp (a [k].c [i], rangeMin [i], rangeMax [i], rangeStep [i]);
    }
  }
}
//————————————————————————————————————————————————————————————————————

El método "GetChiN" calcula la esperanza matemática de la norma de la distribución normal estándar N (0,I) para un número dado de coordenadas (coords). El resultado se aproxima según la dimensionalidad del espacio de búsqueda.

//+------------------------------------------------------------------+
//| Вычисление ожидаемой нормы N(0,I)                                |
//+------------------------------------------------------------------+
double C_AO_CMAES::GetChiN ()
{
  double n = (double)coords;
  return MathSqrt (n) * (1.0 - 1.0 / (4.0 * n) + 1.0 / (21.0 * n * n));
}
//————————————————————————————————————————————————————————————————————

Método  "SortPopulation" está diseñado para clasificar una población de soluciones según los valores de su función objetivo (aptitud) con el fin de determinar la mejor solución. En primer lugar, el método copia los valores de aptitud de cada individuo de la población en un array auxiliar llamado "arfitness". También se crea un array llamado "arindex", que inicialmente contiene los índices de los individuos. Esto es necesario para que, después de clasificar según la aptitud, sea posible restablecer la correspondencia entre la aptitud clasificada y el individuo original en la población.

Se utiliza el algoritmo de clasificación por inserción (insertion sort). Este itera a través del array "arfitness" y, para cada elemento, lo inserta en el lugar correcto de la parte clasificada del array, desplazando los elementos más grandes hacia la derecha. Es importante considerar que la clasificación se realiza en orden descendente (arfitness[j] < tempFitness en el ciclo while), el objetivo es maximizar la función objetivo (fitness). Junto con "arfitness", "arindex" también se clasifica de forma sincrónica para rastrear la posición de los índices originales de los individuos.

Tras la clasificación, el método comprueba si la aptitud del mejor individuo es mejor que la de la mejor solución actual y, en caso afirmativo, se actualiza "fB" y la mejor solución se sustituye por las coordenadas del mejor individuo de la población. De este modo, el método no solo clasifica la población según su aptitud, sino que también realiza un seguimiento de la mejor solución encontrada hasta el momento.

//+------------------------------------------------------------------+
//| Сортировка популяции по фитнесу                                  |
//+------------------------------------------------------------------+
void C_AO_CMAES::SortPopulation ()
{
  // Копирование значений фитнеса и индексов
  for (int i = 0; i < popSize; i++)
  {
    arindex [i] = i;
    arfitness [i] = a [i].f;
  }

  for (int i = 1; i < popSize; i++)
  {
    double tempFitness = arfitness [i];
    double tempIndex = arindex [i];
    int j = i - 1;

    // Сортировка по убыванию (для максимизации)
    while (j >= 0 && arfitness [j] < tempFitness)
    {
      arfitness [j + 1] = arfitness [j];
      arindex [j + 1] = arindex [j];
      j--;
    }

    arfitness [j + 1] = tempFitness;
    arindex [j + 1] = tempIndex;
  }

  // Обновление лучшего решения при необходимости
  if (arfitness [0] > fB)
  {
    fB = arfitness [0];
    int bestIdx = (int)arindex [0];
    ArrayCopy (cB, a [bestIdx].c, 0, 0, coords);
  }
}
//————————————————————————————————————————————————————————————————————

El método  "UpdateMean" actualiza el valor medio de la población (cB) mediante la recombinación ponderada de los mejores individuos. Para cada coordenada, se calcula una nueva media como la suma de los valores de coordenadas ponderadas de los mejores individuos, donde los pesos se especifican en el array " weights" .

//+------------------------------------------------------------------+
//| Обновление среднего с использованием взвешенной рекомбинации     |
//+------------------------------------------------------------------+
void C_AO_CMAES::UpdateMean ()
{
  // Взвешенная рекомбинация: m^(g+1) = Σ w_i * x_{i:λ}^(g+1)
  for (int j = 0; j < coords; j++)
  {
    double meanSum = 0.0;
    for (int i = 0; i < mu; i++)
    {
      int idx = (int)arindex [i];
      meanSum += weights [i] * a [idx].c [j];
    }
    cB [j] = meanSum;
  }
}
//————————————————————————————————————————————————————————————————————

El método ComputeWeights calcula los pesos para la recombinación ponderada, asigna un array "weights" de tamaño "mu" y calcula "log(mu + 0.5)" para usarlo en la fórmula de pesos. Para cada "i" desde "0" hasta "mu-1", asigna un peso como la diferencia entre log(mu + 0.5) y log(i+1); sumando estos pesos, el método divide todos los pesos por la suma para que la suma sea igual a 1, y calcula la suma de los cuadrados de los pesos, calcula "mu_eff", la eficiencia del uso de los pesos, que es importante para ajustar la velocidad y la varianza en el algoritmo.

Este método ofrece pesos óptimos para la recombinación, considerando la contribución de las mejores soluciones.

//+------------------------------------------------------------------+
//| Вычисление весов взвешенной рекомбинации                         |
//+------------------------------------------------------------------+
void C_AO_CMAES::ComputeWeights ()
{
  // Выделение массива весов
  ArrayResize (weights, mu);

  // Предварительное вычисление log(mu + 0.5)
  double log_mu_plus_half = MathLog (mu + 0.5);

  // Вычисление положительных весов
  double sum = 0.0;
  for (int i = 0; i < mu; i++)
  {
    weights [i] = log_mu_plus_half - MathLog (i + 1);
    sum += weights [i];
  }

  // Нормализация весов
  double sum_weights = 0.0;
  double sum_squares = 0.0;
  for (int i = 0; i < mu; i++)
  {
    weights [i] /= sum;
    sum_weights += weights [i];
    sum_squares += weights [i] * weights [i];
  }

  // Вычисление эффективной массы выбора дисперсии
  mu_eff = sum_weights * sum_weights / sum_squares;
}
//————————————————————————————————————————————————————————————————————

El método "UpdateDistribution" actualiza los parámetros de distribución, a saber, la matriz de covarianza (covMatrix) y el tamaño del paso (sigma) en el algoritmo CMA-ES. Si han transcurrido suficientes generaciones (counteval - eigeneval > eigenInterval), entonces se calculará la descomposición de valores propios y la media actual (cB) se almacenará en el array temporal "oldMean". El método llama a "UpdateMean" para calcular una nueva media, calcula la diferencia entre la media nueva y la antigua dividida por "sigma", y también realiza la multiplicación de "y_w" por C^(-1/2), dividida en varios pasos usando la descomposición "B" y "D", y almacena el resultado en "invsqrtC_times_yw".

El método actualiza la ruta de evolución para el tamaño de paso "ps" usando "invsqrtC_times_yw", determina si el progreso se considera "estancado" en función de la norma "ps" y la longitud esperada, y luego actualiza la ruta de evolución para la matriz de covarianza "pc" usando "y_w" y el indicador de progreso "hsig". El método primero calcula una matriz de actualización de rango 1, c1a, luego una matriz de actualización de rango μ, cmu, basada en la diferencia entre los individuos y la media anterior dividida por sigma usando los pesos, actualiza la matriz de covarianza covMatrix usando las actualizaciones de rango 1 y rango μ, y ajusta c1 cuando el progreso se estanca. La función "EnforcePositiveDefinite" se llama periódicamente para garantizar que la matriz de covarianza sea definida de forma positiva. El tamaño del paso "sigma" se actualiza según la norma "ps" y está restringido entre "min_sigma" y "max_sigma" para garantizar la estabilidad numérica.

En general, " UpdateDistribution" ajusta los parámetros de distribución (matriz de covarianza y tamaño del paso) según la historia de búsqueda (trayectoria evolutiva) y los datos de la población para adaptarse al entorno de la función objetivo.

//+------------------------------------------------------------------+ 
//| Обновление параметров распределения                              |
//+------------------------------------------------------------------+
void C_AO_CMAES::UpdateDistribution ()
{
  // Проверка необходимости разложения на собственные значения
  if (counteval - eigeneval > eigenInterval)
  {
    ComputeEigendecomposition ();
    eigeneval = counteval;
  }

  // Сохранение старого среднего
  double oldMean [];
  ArrayResize (oldMean, coords);
  ArrayCopy (oldMean, cB, 0, 0, coords);

  // Обновление среднего
  UpdateMean ();

  // Вычисление смещения среднего
  double y_w [];
  ArrayResize (y_w, coords);
  for (int j = 0; j < coords; j++)
  {
    y_w [j] = (cB [j] - oldMean [j]) / sigma;
  }

  // Вычисление C^(-1/2) * y_w
  // Шаг 1: B^T * y_w
  ArrayInitialize (temp_vec, 0.0);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      temp_vec [i] += B [j * coords + i] * y_w [j];
    }
  }

  // Шаг 2: D^(-1) * (B^T * y_w)
  for (int i = 0; i < coords; i++)
  {
    temp_vec [i] /= D [i];
  }

  // Шаг 3: B * D^(-1) * B^T * y_w
  ArrayInitialize (invsqrtC_times_yw, 0.0);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      invsqrtC_times_yw [i] += B [i * coords + j] * temp_vec [j];
    }
  }

  // Обновление пути эволюции для sigma
  double norm_ps_squared = 0.0;
  for (int i = 0; i < coords; i++)
  {
    ps [i] = (1.0 - cs) * ps [i] + MathSqrt (cs * (2.0 - cs) * mu_eff) * invsqrtC_times_yw [i];
    norm_ps_squared += ps [i] * ps [i];
  }

  // Функция Хевисайда
  double norm_ps = MathSqrt (norm_ps_squared);
  double expected_length = MathSqrt (1.0 - MathPow (1.0 - cs, 2.0 * counteval)) * chiN;
  bool hsig = norm_ps / expected_length < hsig_threshold;

  // Обновление пути эволюции для C
  double delta_hsig = hsig ? 1.0 : 0.0;
  for (int i = 0; i < coords; i++)
  {
    pc [i] = (1.0 - cc) * pc [i] + delta_hsig * MathSqrt (cc * (2.0 - cc) * mu_eff) * y_w [i];
  }

  // Подготовка ранг-1 обновления
  double c1a [];
  ArrayResize (c1a, coords * coords);
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j <= i; j++)
    {
      c1a [i * coords + j] = c1a [j * coords + i] = pc [i] * pc [j];
    }
  }

  // Подготовка ранг-μ обновления
  double cmu [];
  ArrayResize (cmu, coords * coords);
  ArrayInitialize (cmu, 0.0);

  for (int k = 0; k < mu; k++)
  {
    int idx = (int)arindex [k];

    // Вычисление y_i = (x_i - m_old) / sigma
    for (int i = 0; i < coords; i++)
    {
      temp_vec [i] = (a [idx].c [i] - oldMean [i]) / sigma;
    }

    // Добавление взвешенного внешнего произведения
    for (int i = 0; i < coords; i++)
    {
      for (int j = 0; j <= i; j++)
      {
        double update = weights [k] * temp_vec [i] * temp_vec [j];
        cmu [i * coords + j] += update;
        if (i != j) cmu [j * coords + i] += update;
      }
    }
  }

  // Обновление ковариационной матрицы C
  double c1 = learningRateC1;
  double cmu_rate = learningRateCMu;

  // Корректировка c1 если hsig false (застой прогресса)
  if (!hsig)
  {
    c1 *= (1.0 - (1.0 - delta_hsig) * cc * (2.0 - cc));
  }

  double one_minus_c1_cmu = 1.0 - c1 - cmu_rate;

  // Обновление C с ранг-1 и ранг-μ обновлениями
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j <= i; j++)
    {
      covMatrix [i * coords + j] = one_minus_c1_cmu * covMatrix [i * coords + j] +
                                   c1 * c1a [i * coords + j] +
                                   cmu_rate * cmu [i * coords + j];

      // Сохранение симметрии
      if (i != j)
      {
        covMatrix [j * coords + i] = covMatrix [i * coords + j];
      }
    }
  }

  // Обеспечение положительной определенности
  if (counteval % (10 * eigenInterval) == 0)
  {
    EnforcePositiveDefinite ();
  }

  // Обновление размера шага sigma
  double exponent = (cs / damps) * (norm_ps / chiN - 1.0);
  sigma *= MathExp (exponent);

  // Ограничение sigma для численной стабильности
  double min_sigma = 1e-16;
  double max_eigenvalue = D [0]; // D отсортирован по убыванию
  double max_sigma = 1e4 * MathMax (1.0, MathSqrt (max_eigenvalue));

  if (sigma < min_sigma) sigma = min_sigma;
  else
    if (sigma > max_sigma) sigma = max_sigma;
}
//————————————————————————————————————————————————————————————————————

El método  "ComputeEigendecomposition" realiza la descomposición de la matriz de covarianza simétrica " covMatrix" en valores propios y vectores propios utilizando el método Jacobi mejorado; el método crea una copia de la "covMatrix" para no modificar los datos originales. Inicialmente, "B" se define como la matriz identidad que representa la base inicial. Realiza el número máximo de iteraciones o hasta que los elementos distintos de cero en la diagonal principal sean menores que "tolerance". 

Búsqueda del elemento máximo fuera de la diagonal: encuentra el elemento con el valor absoluto máximo fuera de la diagonal para seleccionarlo para la rotación. A continuación, si el valor máximo es inferior a la "tolerance", el ciclo se interrumpe. El método calcula "phi" para realizar la rotación jacobiana, modifica los elementos de "C_copy" aplicando la rotación y actualiza las columnas de la matriz "B" (vectores propios) de acuerdo con la rotación. Tras varias iteraciones, el método halla la raíz cuadrada de los elementos diagonales (garantizando una positividad mínima de 1e-14) y clasifica los valores propios y los vectores propios correspondientes en orden descendente para su posterior uso. De este modo, el método permite calcular con precisión los valores propios y los vectores propios de la matriz de covarianza.

//+------------------------------------------------------------------+
//| Вычисление разложения на собственные значения методом Якоби      |
//+------------------------------------------------------------------+
void C_AO_CMAES::ComputeEigendecomposition ()
{
  // Создание копии ковариационной матрицы для разложения
  double C_copy [];
  ArrayResize (C_copy, coords * coords);
  ArrayCopy (C_copy, covMatrix);

  // Инициализация B как единичной матрицы
  for (int i = 0; i < coords; i++)
  {
    for (int j = 0; j < coords; j++)
    {
      B [i * coords + j] = (i == j) ? 1.0 : 0.0;
    }
  }

  // Улучшенное разложение Якоби на собственные значения
  int max_iterations = 10; //50 * coords;
  double tolerance   = 0.01; //1e-14 * coords * coords;

  for (int iter = 0; iter < max_iterations; iter++)
  {
    // Поиск наибольшего вне диагонального элемента
    double max_val = 0.0;
    int p = 0, q = 1;

    for (int i = 0; i < coords - 1; i++)
    {
      for (int j = i + 1; j < coords; j++)
      {
        double val = MathAbs (C_copy [i * coords + j]);
        if (val > max_val)
        {
          max_val = val;
          p = i;
          q = j;
        }
      }
    }

    // Проверка сходимости
    if (max_val < tolerance) break;

    // Вычисление угла поворота
    double app = C_copy [p * coords + p];
    double aqq = C_copy [q * coords + q];
    double apq = C_copy [p * coords + q];

    double phi = 0.5 * MathArctan (2.0 * apq / (aqq - app + 1e-14));
    double c = MathCos (phi);
    double s = MathSin (phi);

    // Обновление элементов матрицы
    double app_new = c * c * app - 2 * c * s * apq + s * s * aqq;
    double aqq_new = s * s * app + 2 * c * s * apq + c * c * aqq;

    C_copy [p * coords + p] = app_new;
    C_copy [q * coords + q] = aqq_new;
    C_copy [p * coords + q] = C_copy [q * coords + p] = 0.0;

    // Обновление других элементов в строках/столбцах p и q
    for (int i = 0; i < coords; i++)
    {
      if (i != p && i != q)
      {
        double aip = C_copy [i * coords + p];
        double aiq = C_copy [i * coords + q];
        C_copy [i * coords + p] = C_copy [p * coords + i] = c * aip - s * aiq;
        C_copy [i * coords + q] = C_copy [q * coords + i] = s * aip + c * aiq;
      }
    }

    // Обновление собственных векторов
    for (int i = 0; i < coords; i++)
    {
      double bip = B [i * coords + p];
      double biq = B [i * coords + q];
      B [i * coords + p] = c * bip - s * biq;
      B [i * coords + q] = s * bip + c * biq;
    }
  }

  // Извлечение собственных значений и обеспечение положительности
  double min_eigenvalue = 1e-14;
  for (int i = 0; i < coords; i++)
  {
    D [i] = MathSqrt (MathMax (min_eigenvalue, C_copy [i * coords + i]));
  }

  // Сортировка собственных значений и векторов по убыванию
  for (int i = 0; i < coords - 1; i++)
  {
    int max_idx = i;
    for (int j = i + 1; j < coords; j++)
    {
      if (D [j] > D [max_idx]) max_idx = j;
    }

    if (max_idx != i)
    {
      // Обмен собственных значений
      double temp = D [i];
      D [i] = D [max_idx];
      D [max_idx] = temp;

      // Обмен собственных векторов
      for (int k = 0; k < coords; k++)
      {
        temp = B [k * coords + i];
        B [k * coords + i] = B [k * coords + max_idx];
        B [k * coords + max_idx] = temp;
      }
    }
  }
}
//————————————————————————————————————————————————————————————————————

El método " CheckPositiveDefinite" comprueba si la matriz de covarianza es definida positiva. Comprueba rápidamente la positividad de los elementos diagonales y compara el valor propio mínimo (del array " D" , clasificado en orden descendente) con un número positivo pequeño 1e-14 . Si ambas comprobaciones se superan, el método devuelve "true".

//+------------------------------------------------------------------+
//| Проверка положительной определенности ковариационной матрицы     |
//+------------------------------------------------------------------+
bool C_AO_CMAES::CheckPositiveDefinite ()
{
  // Быстрая проверка: все диагональные элементы должны быть положительными
  for (int i = 0; i < coords; i++)
  {
    if (covMatrix [i * coords + i] <= 0) return false;
  }

  // Проверка минимального собственного значения на положительность
  double min_eigenvalue = D [coords - 1]; // D отсортирован по убыванию
  return min_eigenvalue > 1e-14;
}
//————————————————————————————————————————————————————————————————————

El método " EnforcePositiveDefinite" garantiza la positividad definida de la matriz de covarianza " covMatrix". Realiza varios pasos: encuentra el elemento mínimo en la diagonal y agrega el valor de "correction" faltante a los elementos diagonales para que el mínimo sea igual a 1e-10, maximiza la simetría de la matriz promediando cada elemento fuera de la diagonal con su transpuesta.

Si la matriz todavía no es definida positiva, se realiza una descomposición en valores propios utilizando "ComputeEigendecomposition", reemplazando todos los valores propios menores que √1e-10 con √1e-10 para garantizar la positividad. El método luego recalcula "covMatrix" usando los vectores propios de "B" y los valores propios corregidos de "D" para obtener una matriz definida positiva corregida. Este enfoque garantiza que la matriz de covarianza sea definida positiva y adecuada para cálculos posteriores en el algoritmo CMA-ES.

//+------------------------------------------------------------------+
//| Обеспечение положительной определенности ковариационной матрицы  |
//+------------------------------------------------------------------+
void C_AO_CMAES::EnforcePositiveDefinite ()
{
  // Метод 1: Добавление малого значения к диагонали
  double min_diag = 1e308; // Очень большое число
  for (int i = 0; i < coords; i++)
  {
    if (covMatrix [i * coords + i] < min_diag)
    {
      min_diag = covMatrix [i * coords + i];
    }
  }

  if (min_diag < 1e-10)
  {
    double correction = 1e-10 - min_diag;
    for (int i = 0; i < coords; i++)
    {
      covMatrix [i * coords + i] += correction;
    }
  }

  // Метод 2: Обеспечение симметрии
  for (int i = 0; i < coords; i++)
  {
    for (int j = i + 1; j < coords; j++)
    {
      double avg = (covMatrix [i * coords + j] + covMatrix [j * coords + i]) * 0.5;
      covMatrix [i * coords + j] = covMatrix [j * coords + i] = avg;
    }
  }

  // Если все еще не положительно определена, выполнить разложение и исправить
  if (!CheckPositiveDefinite ())
  {
    ComputeEigendecomposition ();

    double min_eigenvalue = 1e-10;
    for (int i = 0; i < coords; i++)
    {
      if (D [i] < MathSqrt (min_eigenvalue))
      {
        D [i] = MathSqrt (min_eigenvalue);
      }
    }

    // Восстановление C = B * D^2 * B^T
    ArrayInitialize (covMatrix, 0.0);
    for (int i = 0; i < coords; i++)
    {
      for (int j = 0; j <= i; j++)
      {
        double sum = 0.0;
        for (int k = 0; k < coords; k++)
        {
          sum += B [i * coords + k] * D [k] * D [k] * B [j * coords + k];
        }
        covMatrix [i * coords + j] = covMatrix [j * coords + i] = sum;
      }
    }
  }
}
//————————————————————————————————————————————————————————————————————


Resultados de las pruebas

Ahora por fin podemos ver los resultados. Como podemos observar de inmediato, el algoritmo funciona bien y con rapidez con problemas de tamaño mediano y, en ciertas configuraciones, con problemas de baja dimensionalidad, pero los problemas multidimensionales más difíciles requieren una cantidad desproporcionada de tiempo para ejecutarse, por lo que fueron excluidos de los resultados. ¿Con qué está relacionado esto? Tradicionalmente, los cálculos matriciales en MQL5 presentan un grave problema computacional. Al implementar operaciones matriciales manualmente, el algoritmo muestra un rendimiento insatisfactorio en problemas de alta dimensionalidad, lo cual limita significativamente su aplicabilidad práctica. Sin embargo, con la llegada de las clases integradas para trabajar con matrices, la situación cambia drásticamente. Ahora es fundamental usar la implementación nativa de la plataforma para las operaciones matriciales en todas las tareas que requieren una gran capacidad de cálculo.

CMAES|Covariance Matrix Adaptation Evolution Strategy|50.0|25.0|0.01|0.8|1.0|
=============================
5 Hilly's; Func runs: 10000; result: 0,7625797883550075
25 Hilly's; Func runs: 10000; result: 0,7208874560706138
500 Hilly's; Func runs: 10000; result: 0.0
=============================
5 Forest's; Func runs: 10000; result: 0,8205636421348295
25 Forest's; Func runs: 10000; result: 0,7961602627346933
500 Forest's; Func runs: 10000; result: 0.0
=============================
5 Megacity's; Func runs: 10000; result: 0,7584615384615383
25 Megacity's; Func runs: 10000; result: 0,49076923076923074
500 Megacity's; Func runs: 10000; result: 0.0
=============================
All score: 4,34942 (48,33%)

En la visualización, podemos ver que el algoritmo funciona tanto con "grandes saltos" como centrándose en los extremos de la función, examinándolos cuidadosamente.

Hilly

CMA-ES en la función de prueba Hilly

Forest

CMA-ES en la función de prueba Forest

Megacity

CMA-ES en la función de prueba Megacity

Según los resultados del trabajo, el algoritmo CMA-ES ocupa el puesto 38 en la clasificación de los algoritmos de optimización basados en poblaciones más potentes.

AO
Description
Hilly
Hilly
Final
Forest
Forest
Final
Megacity (discrete)
Megacity
Final
Final
Result
% de
MAX
10 p (5 F) 50 p (25 F) 1000 p (500 F) 10 p (5 F) 50 p (25 F) 1000 p (500 F) 10 p (5 F) 50 p (25 F) 1000 p (500 F)
1 ANS across neighbourhood search 0,94948 0,84776 0,43857 2,23581 1,00000 0,92334 0,39988 2,32323 0,70923 0,63477 0,23091 1,57491 6,134 68,15
2 CLA code lock algorithm (joo) 0,95345 0,87107 0,37590 2,20042 0,98942 0,91709 0,31642 2,22294 0,79692 0,69385 0,19303 1,68380 6,107 67,86
3 AMOm animal migration ptimization M 0,90358 0,84317 0,46284 2,20959 0,99001 0,92436 0,46598 2,38034 0,56769 0,59132 0,23773 1,39675 5,987 66,52
4 (P+O)ES (P+O) evolution strategies 0,92256 0,88101 0,40021 2,20379 0,97750 0,87490 0,31945 2,17185 0,67385 0,62985 0,18634 1,49003 5,866 65,17
5 CTA comet tail algorithm (joo) 0,95346 0,86319 0,27770 2,09435 0,99794 0,85740 0,33949 2,19484 0,88769 0,56431 0,10512 1,55712 5,846 64,96
6 TETA time evolution travel algorithm (joo) 0,91362 0,82349 0,31990 2,05701 0,97096 0,89532 0,29324 2,15952 0,73462 0,68569 0,16021 1,58052 5,797 64,41
7 SDSm stochastic diffusion search M 0,93066 0,85445 0,39476 2,17988 0,99983 0,89244 0,19619 2,08846 0,72333 0,61100 0,10670 1,44103 5,709 63,44
8 BOAm billiards optimization algorithm M 0,95757 0,82599 0,25235 2,03590 1,00000 0,90036 0,30502 2,20538 0,73538 0,52523 0,09563 1,35625 5,598 62,19
9 AAm archery algorithm M 0,91744 0,70876 0,42160 2,04780 0,92527 0,75802 0,35328 2,03657 0,67385 0,55200 0,23738 1,46323 5,548 61,64
10 ESG evolution of social groups (joo) 0,99906 0,79654 0,35056 2,14616 1,00000 0,82863 0,13102 1,95965 0,82333 0,55300 0,04725 1,42358 5,529 61,44
11 SIA simulated isotropic annealing (joo) 0,95784 0,84264 0,41465 2,21513 0,98239 0,79586 0,20507 1,98332 0,68667 0,49300 0,09053 1,27020 5,469 60,76
12 BBO biogeography based optimization 0,94912 0,69456 0,35031 1.99399 0,93820 0,67365 0,25682 1,86867 0,74615 0,48277 0,17369 1.40261 5.265 58,50
13 ACS artificial cooperative search 0,75547 0,74744 0,30407 1,80698 1,00000 0,88861 0,22413 2,11274 0,69077 0,48185 0,13322 1,30583 5,226 58,06
14 DA dialectical algorithm 0,86183 0,70033 0,33724 1,89940 0,98163 0,72772 0,28718 1,99653 0,70308 0,45292 0,16367 1,31967 5,216 57,95
15 BHAm black hole algorithm M 0,75236 0,76675 0,34583 1,86493 0,93593 0,80152 0,27177 2,00923 0,65077 0,51646 0,15472 1,32195 5.196 57.73
16 ASO anarchy society optimization 0,84872 0,74646 0,31465 1,90983 0,96148 0,79150 0,23803 1,99101 0,57077 0,54062 0,16614 1,27752 5,178 57,54
17 RFO royal flush optimization (joo) 0,83361 0,73742 0,34629 1,91733 0,89424 0,73824 0,24098 1,87346 0,63154 0,50292 0,16421 1,29867 5,089 56,55
18 AOSm atomic orbital search M 0,80232 0,70449 0,31021 1,81702 0,85660 0,69451 0,21996 1,77107 0,74615 0,52862 0,14358 1,41835 5,006 55,63
19 TSEA turtle shell evolution algorithm (joo) 0,96798 0,64480 0,29672 1,90949 0,99449 0,61981 0,22708 1,84139 0,69077 0,42646 0,13598 1,25322 5,004 55,60
20 DE differential evolution 0,95044 0,61674 0,30308 1,87026 0,95317 0,78896 0,16652 1,90865 0,78667 0,36033 0,02953 1,17653 4,955 55,06
21 SRA successful restaurateur algorithm (joo) 0,96883 0,63455 0,29217 1,89555 0,94637 0,55506 0,19124 1,69267 0,74923 0,44031 0,12526 1,31480 4,903 54,48
22 CRO chemical reaction optimisation 0,94629 0,66112 0,29853 1,90593 0,87906 0,58422 0,21146 1,67473 0,75846 0,42646 0,12686 1,31178 4,892 54,36
23 BIO blood inheritance optimization (joo) 0,81568 0,65336 0,30877 1,77781 0,89937 0,65319 0,21760 1,77016 0,67846 0,47631 0,13902 1,29378 4,842 53,80
24 BSA bird swarm algorithm 0,89306 0,64900 0,26250 1,80455 0,92420 0,71121 0,24939 1,88479 0,69385 0,32615 0,10012 1,12012 4,809 53,44
25 HS harmony search 0,86509 0,68782 0,32527 1,87818 0,99999 0,68002 0,09590 1,77592 0,62000 0,42267 0,05458 1,09725 4,751 52,79
26 SSG saplings sowing and growing 0,77839 0,64925 0,39543 1,82308 0,85973 0,62467 0,17429 1,65869 0,64667 0,44133 0,10598 1,19398 4,676 51,95
27 BCOm bacterial chemotaxis optimization M 0,75953 0,62268 0,31483 1,69704 0,89378 0,61339 0,22542 1,73259 0,65385 0,42092 0,14435 1,21912 4,649 51,65
28 ABO african buffalo optimization 0,83337 0,62247 0,29964 1,75548 0,92170 0,58618 0,19723 1,70511 0,61000 0,43154 0,13225 1,17378 4,634 51,49
29 (PO)ES (PO) evolution strategies 0,79025 0,62647 0,42935 1,84606 0,87616 0,60943 0,19591 1,68151 0,59000 0,37933 0,11322 1,08255 4,610 51,22
30 FBA fractal-based Algorithm 0,79000 0,65134 0,28965 1.73099 0,87158 0,56823 0,18877 1.62858 0,61077 0,46062 0,12398 1,19537 4.555 50.61
31 TSm tabu search M 0,87795 0,61431 0,29104 1,78330 0,92885 0,51844 0,19054 1,63783 0,61077 0,38215 0,12157 1,11449 4,536 50,40
32 BSO brain storm optimization 0,93736 0,57616 0,29688 1,81041 0,93131 0,55866 0,23537 1,72534 0,55231 0,29077 0,11914 0,96222 4,498 49,98
33 WOAm wale optimization algorithm M 0,84521 0,56298 0,26263 1,67081 0,93100 0,52278 0,16365 1,61743 0,66308 0,41138 0,11357 1,18803 4,476 49,74
34 AEFA artificial electric field algorithm 0,87700 0,61753 0,25235 1,74688 0,92729 0,72698 0,18064 1,83490 0,66615 0,11631 0,09508 0,87754 4,459 49,55
35 AEO artificial ecosystem-based optimization algorithm 0,91380 0,46713 0,26470 1,64563 0,90223 0,43705 0,21400 1,55327 0,66154 0,30800 0,28563 1,25517 4,454 49,49
36 CAm camel algorithm M 0,78684 0,56042 0,35133 1.69859 0,82772 0,56041 0,24336 1.63149 0,64846 0,33092 0,13418 1,11356 4.444 49.37
37 ACOm ant colony optimization M 0,88190 0,66127 0,30377 1,84693 0,85873 0,58680 0,15051 1,59604 0,59667 0,37333 0,02472 0,99472 4,438 49,31
38 CMAES covariance_matrix_adaptation_evolution_strategy 0,76258 0,72089 0,00000 1,48347 0,82056 0,79616 0,00000 1.61672 0,75846 0,49077 0,00000 1.24923 4.349 48.33
39 BFO-GA bacterial foraging optimization - ga 0,89150 0,55111 0,31529 1,75790 0,96982 0,39612 0,06305 1,42899 0,72667 0,27500 0,03525 1,03692 4,224 46,93
40 SOA simple optimization algorithm 0,91520 0,46976 0,27089 1,65585 0,89675 0,37401 0,16984 1,44060 0,69538 0,28031 0,10852 1,08422 4.181 46,45
41 ABHA artificial bee hive algorithm 0,84131 0,54227 0,26304 1,64663 0,87858 0,47779 0,17181 1,52818 0,50923 0,33877 0,10397 0,95197 4.127 45,85
42 ACMO atmospheric cloud model optimization 0,90321 0,48546 0,30403 1,69270 0,80268 0,37857 0,19178 1,37303 0,62308 0,24400 0,10795 0,97503 4,041 44,90
43 ADAMm adaptive moment estimation M 0,88635 0,44766 0,26613 1.60014 0,84497 0,38493 0,16889 1,39880 0,66154 0,27046 0,10594 1,03794 4.037 44.85
44 CGO chaos game optimization 0,57256 0,37158 0,32018 1,26432 0,61176 0,61931 0,62161 1,85267 0,37538 0,21923 0,19028 0,78490 3,902 43,35
45 CROm coral reefs optimization M 0,78512 0,46032 0,25958 1,50502 0,86688 0,35297 0,16267 1,38252 0,63231 0,26738 0,10734 1,00703 3,895 43,27
RW random walk 0,48754 0,32159 0,25781 1,06694 0,37554 0,21944 0,15877 0,75375 0,27969 0,14917 0,09847 0,52734 2.348 26.09


Conclusiones

La capacidad del CMA-ES para adaptar la estrategia de búsqueda a la geometría local de la función objetivo la hace prácticamente indispensable para problemas complejos y mal condicionados con estructura desconocida. Las colas pesadas de la distribución de ley de potencias permiten que el algoritmo realice "saltos largos", lo que resulta fundamental para escapar de los óptimos locales. La complejidad computacional de O(n²) en la memoria y O(n³) en el tiempo limita severamente la aplicabilidad del algoritmo. Para problemas de alta dimensionalidad (n > 100), la intensidad de los recursos se vuelve desproporcionada con respecto a los beneficios obtenidos. El tiempo de ejecución en funciones multidimensionales crece tan rápidamente que hace que el algoritmo resulte prácticamente inaplicable para valores grandes de n.

El CMA-ES funciona con funciones ruidosas, maneja funciones objetivo discontinuas, es eficaz en entornos con múltiples extremos, y el secreto de esta versatilidad reside en la filosofía fundamental del algoritmo: en lugar de hacer suposiciones sobre la estructura de la función, la explora cuidadosamente, adaptando su estrategia de búsqueda según la experiencia acumulada. El CMA-ES transforma el panorama de la computación evolutiva al demostrar que los algoritmos inspirados en la biología pueden tener fundamentos matemáticos rigurosos y que la adaptación de algoritmos no es simplemente un ajuste de parámetros, sino un principio fundamental para aprender la estructura de un problema.

El algoritmo es idóneo para problemas de tamaño mediano, pues supone una herramienta especializada que parece una alternativa valiosa en su nicho. El algoritmo combina elegancia matemática con eficiencia práctica, pero requiere una aplicación cuidadosa debido a sus limitaciones computacionales. En estos casos, los desarrolladores del lenguaje MQL5 han introducido recientemente operaciones rápidas con matrices para calcularlas; en esta implementación del algoritmo, se usaron cálculos secuenciales convencionales.

tab

Figura 2. Clasificación por colores de los algoritmos según las pruebas respectivas

chart

Figura 3. Histograma de los resultados de las pruebas de algoritmos (en una escala de 0 a 100, cuanto más mejor, donde 100 es el máximo resultado teórico posible, el script para calcular la tabla de puntuación está en el archivo)

Ventajas y desventajas del algoritmo CMAES:

Ventajas:

  1. Buena convergencia en funciones de dimensionalidad media.

Desventajas:

  1. Dispersión de resultados en funciones de baja dimensionalidad.
  2. Consumo crítico de recursos en funciones de alta dimensionalidad.

Adjuntamos al artículo un archivo con las versiones actuales de los códigos de los algoritmos. El autor de este artículo no se responsabiliza de la exactitud absoluta de la descripción de los algoritmos canónicos, muchos de ellos han sido modificados para mejorar las capacidades de búsqueda. Las conclusiones y juicios presentados en los artículos se basan en los resultados de los experimentos realizados.


Programas usados en el artículo

# Nombre Tipo Descripción
1 #C_AO.mqh
Archivo de inclusión
Clase padre de algoritmos de optimización basados en la población
2 #C_AO_enum.mqh
Archivo de inclusión
Enumeración de los algoritmos de optimización basados en la población
3 TestFunctions.mqh
Archivo de inclusión
Biblioteca de funciones de prueba
4
TestStandFunctions.mqh
Archivo de inclusión
Biblioteca de funciones del banco de pruebas
5
Utilities.mqh
Archivo de inclusión
Biblioteca de funciones auxiliares
6
CalculationTestResults.mqh
Archivo de inclusión
Script para calcular los resultados en una tabla comparativa
7
Testing AOs.mq5
Script Banco de pruebas único para todos los algoritmos de optimización basados en la población
8
Simple use of population optimization algorithms.mq5
Script
Ejemplo sencillo de utilización de algoritmos de optimización basados en la población sin visualización
9
Test_AO_CMAES.mq5
Script Banco de pruebas para CMAES

Traducción del ruso hecha por MetaQuotes Ltd.
Artículo original: https://www.mql5.com/ru/articles/18227

Archivos adjuntos |
CMAES.ZIP (297.43 KB)
Utilizando redes neuronales en MetaTrader Utilizando redes neuronales en MetaTrader
En el artículo se muestra la aplicación de las redes neuronales en los programas de MQL, usando la biblioteca de libre difusión FANN. Usando como ejemplo una estrategia que utiliza el indicador MACD se ha construido un experto que usa el filtrado con red neuronal de las operaciones. Dicho filtrado ha mejorado las características del sistema comercial.
Características del Wizard MQL5 que debe conocer (Parte 65): Uso de los patrones FrAMA y Force Index Características del Wizard MQL5 que debe conocer (Parte 65): Uso de los patrones FrAMA y Force Index
La media móvil adaptativa fractal (FrAMA) y el oscilador Force Index son otro par de indicadores que podrían utilizarse conjuntamente en un asesor experto de MQL5. Estos dos indicadores se complementan en cierta medida, ya que FrAMA es un indicador de seguimiento de tendencias, mientras que el Force Index es un oscilador basado en el volumen. Como siempre, utilizamos el asistente MQL5 para explorar rápidamente el potencial que puede tener esta combinación.
Particularidades del trabajo con números del tipo double en MQL4 Particularidades del trabajo con números del tipo double en MQL4
En estos apuntes hemos reunido consejos para resolver los errores más frecuentes al trabajar con números del tipo double en los programas en MQL4.
Redes neuronales en el trading: Pipeline inteligente de previsiones (Mezcla dispersa de expertos) Redes neuronales en el trading: Pipeline inteligente de previsiones (Mezcla dispersa de expertos)
Hoy le proponemos familiarizarnos con la implementación práctica de un bloque de mezcla dispersa de expertos para series temporales en el entorno de computación OpenCL. Este artículo explica paso a paso el funcionamiento de la convolución multiventana enmascarada, así como la organización del aprendizaje de gradientes en condiciones de múltiples flujos de información.