Covariance Matrix Adaptation Evolution Strategy (CMA-ES)
1 491
0
Contents
Introduction
There are many algorithms in the world of optimization, but we are looking for the most powerful ones to solve the optimization problems of our trading robots. CMA-ES (Covariance Matrix Adaptation Evolution Strategy) is one of those rare examples where mathematical rigor is combined with biological intuition, creating an algorithm that not only solves optimization problems, but also learns to understand their structure.
The history of CMA-ES began in the late 1990s in research laboratories in Germany, where Nikolaus Hansen and Andreas Ostermeier posed a fundamental question: Is it possible to create an optimization algorithm that does not simply search for a solution, but adapts to the geometry of the problem? Traditional evolutionary algorithms generated offspring in spherical regions around parent individuals, which worked well for simple functions but proved ineffective for complex, ill-conditioned problems. Let's take a look at this interesting algorithm.
Implementation of the algorithm
Imagine searching for treasure on an irregularly shaped island. The usual approach is to search in all directions equally, as if the island were round. CMA-ES gradually learns the shape of the island and gradually shifts its search toward directions where the likelihood of finding treasure is higher. Moreover, it remembers successful routes and uses this memory to plan future searches.
CMA-ES is based on a deceptively simple equation: x_k ~ N(m, σ²C). But behind this simplicity lies a deep mathematical structure. Each symbol here carries important information about the state of the search: m is the current best guess about the location of the optimum, σ is a measure of how far we are willing to risk moving away from the known, while C is a covariance matrix that encodes our understanding of the function geometry. The only change we can justify is replacing the normal distribution with a power distribution, which means that the implementation will follow a modified equation: x_k ~ PowerDist(m, σ²C). This modification changes the nature of the space exploration (wider "jumps"), but preserves the fundamental adaptive nature of the algorithm.
The covariance matrix C is the true heart of the algorithm. It begins its life as a humble identity matrix representing a spherical distribution. But with each iteration it evolves, stretching along directions of rapid improvement and shrinking where progress is slow. Gradually, the sphere turns into an ellipse, then into an elongated ellipsoid, ideally oriented along the contours of the function being optimized.
The main innovation of CMA-ES is the concept of evolutionary paths. This is a kind of genetic memory of the algorithm, which remembers not only where the successful points were, but also how the algorithm arrived at them. The evolutionary path accumulates information about successive successful steps, creating a directional vector that points to the most promising search areas. The second evolutionary path performs a more subtle function: it controls step size. If the successive steps of the algorithm are correlated, that is, each subsequent step continues the direction of the previous one, then the step size increases - the algorithm "feels" that it is moving in the right direction. If the steps are random and uncorrelated, the step size decreases - perhaps the algorithm is already close to the optimum and needs to be searched more carefully.
Behind the biological metaphor of evolution in CMA-ES lies a rigorous mathematical principle - maximum likelihood estimation. The algorithm constantly asks itself: What distribution parameters make the observed successful points most probable? This transforms evolutionary optimization from a heuristic into a statistically based method.
Each update of the C covariance matrix consists of two components: rank-one update based on the evolutionary path and rank-μ update using information from the entire selected population. Rank-one update provides stability and the ability to capture long-term trends, while rank-μ update allows rapid adaptation to new information from the current generation.
In the context of CMA-ES, a modified Heaviside function is used, which plays an important role in the stagnation detection mechanism of the algorithm. The function compares the length of the evolutionary path with the expected length; if the path is too short, it is a sign of "staggering". Deactivation conditions (hsig = 0): the algorithm "wanders" randomly and the steps cancel each other out, then the step size is probably too large. During stagnation, the influence of rank-one update decreases and we rely more on information from the entire population. Activation conditions (hsig = 1): the algorithm makes progress in a certain direction, successive steps are correlated, then the step size is adequate to the current situation.
One of the deepest properties of CMA-ES is its invariance to transformations of the search space. The algorithm equally efficiently solves the function f(x) and the function f(Ax + b), where A is any non-singular matrix and b is any shift vector. This means that CMA-ES is independent of the choice of coordinate system. This invariance is not accidental; it is a direct consequence of the use of the maximum likelihood principle and the adaptation of the covariance matrix. The algorithm automatically detects a natural coordinate system for the problem, where the axes coincide with the principal directions of the function variation.
The beauty of theory must be combined with practical applicability. CMA-ES requires O(n²) of memory for storing the covariance matrix and O(n³) computations for eigendecomposition, which makes the algorithm applicable to problems with dimensions of up to several hundred variables. For large dimensions, specialized modifications have been developed: sep-CMA-ES is limited to diagonal covariance matrices, VkD-CMA uses variable dimensions, and LM-CMA applies the principles of limited memory. The implementation of these methods is currently beyond the scope of our article, but it may be possible to return to them in subsequent articles.
Figure 1. Illustration of the CMA-ES algorithm
The illustration shows the evolution through generations - three snapshots of the state of the algorithm (generations 1, 5 and 15), showing how the population gradually converges to the optimum,
Adaptation of the covariance matrix - the blue ellipses show how the shape of the distribution changes when the generation is: 1- round (identity matrix), 5 - elongated and rotated, 15 - finely tuned to local geometry.
Algorithm components: red dots are the entire population (λ descendants), green dots are the best solutions (μ parents), black dot m is the population mean, dashed circles are the isolines of the objective function.
The key formula is below, with a note about the power law modification. The illustration clearly demonstrates how the algorithm adapts its search strategy to the landscape of the function being optimized, gradually narrowing the search space and approaching the optimum.
Let's prepare the pseudocode.
Initialization
- Set algorithm parameters:
- Population size (lambda) = 50
- Number of parents (mu) = 25
- Learning rate for rank-1 update (c1) = 0.01
- Learning rate for rank-μ update (cμ) = 0.8
- Step size damping = 0.6
- Initial step size (sigma) = 0.3
- Calculate recombination weights:
- For each of the μ parents, calculate the weight as log(μ + 0.5) - log(i + 1)
- Normalize the weights so that their sum equals 1
- Calculate the effective selection mass of μ_eff
- Initialize strategy parameters:
- Calculate learning rates for evolution path (cs, cc)
- Compute the expected norm of a standard normal vector (chiN)
- Determine the interval for eigenvalue decomposition
- Initialize data structures:
- C covariance matrix = identity matrix
- B matrix of eigenvectors = identity matrix
- D vector of eigenvalues = unit vector
- pc and ps evolution paths = zero vectors
- Population mean m = random point in search space
- Create an initial population:
- Generate lambda random points in the search space
Basic cycle of evolution
Repeat until stop criterion is reached:
Step 1: Generating offspring
- For each of lambda descendants:
- Generate the z random vector from the standard normal distribution
- Apply the transformation: y = B × D × z
- Create a descendant: x_k = m + σ × y
- Apply search space constraints to x_k
Step 2: Evaluate and update
- Assess the fitness of all descendants
- Sort the population:
- Sort descendants in descending order of fitness
- Update the best solution found if necessary
- Update population mean:
- Store the previous mean m_old
- Calculate the new mean as the weighted sum of the μ best offspring: m_new = Σ(w_i × x_i) for i from 1 to μ
- Update evolution paths:
- Calculate the mean shift: y_w = (m_new - m_old) / σ
- Calculate C^(-1/2) × y_w via eigenvalue decomposition
- Update evolution path for step size: ps = (1 - cs) × ps + √(cs × (2 - cs) × μ_eff) × C^(-1/2) × y_w
- Check the stagnation condition (Heaviside function)
- Update the evolution path for the covariance matrix: pc = (1 - cc) × pc + stagnation_indicator × √(cc × (2 - cc) × μ_eff) × y_w
- Update the covariance matrix:
- Prepare rank-1 update from the outer product of pc
- Prepare the rank-μ update from the weighted sum of the best descendant's deviations
- Update C: C = (1 - c1 - cμ) × C + c1 × (pc × pc^T) + cμ × Σ(w_i × y_i × y_i^T)
- Ensure the matrix symmetry
- Update the step size:
- Calculate the adaptation coefficient based on the ps path length
- Update σ: σ = σ × exp((cs/damps) × (||ps||/chiN - 1))
- Limit σ within reasonable limits for numerical stability
- Eigenvalue decomposition (if necessary):
- If enough iterations have passed since the last decomposition:
- Perform the Jacobi decomposition of C = B × D² × B^T
- Sort eigenvalues and eigenvectors in descending order
- Ensure positive definiteness of the matrix
- If enough iterations have passed since the last decomposition:
- Check and correct the covariance matrix:
- Periodically check the C positive definiteness
- If necessary, add a small value to the diagonal
- Ensure the matrix symmetry
Now let's write the C_AO_CMAES derived from the C_AO class, which will implement the CMA-ES optimization algorithm.
Public fields, method:
- SetParams () - set CMA-ES parameter values from the params array.
- Init () - initializes the algorithm. Accepts minimum, maximum values and step size for each variable, as well as the number of epochs.
- Moving () - implements the main loop of the algorithm, responsible for generating new solutions.
- Revision () - evaluation of new solutions and updating of the algorithm state.
Private fields:
- CMA-ES parameters - the variables containing specific parameters of the CMA-ES algorithm: number of parents (mu), step size (sigma), learning rates (learningRateC1, learningRateCMu), damping factor (stepSizeDamping ) and others necessary for the algorithm to work.
- CMA-ES data structures - arrays used to store recombination weights, covMatrix covariance matrix, B eigenvectors, D eigenvalues, evolution paths (pc, ps), as well as other auxiliary data. mu_eff, counteval, eigeneval, chiN, eigenInterval values are used to control and manage the algorithm progress.
- Variables for performance optimization are designed to speed up calculations during the algorithm's operation, such as learning rates and damping.
- Auxiliary arrays used as temporary storage for vectors and matrices when performing various operations.
- Auxiliary methods, which perform individual steps of the CMA-ES algorithm, such as initializing the distribution, updating the distribution, computing the eigendecomposition, sorting the population, computing the weights, updating the mean, checking for positive definiteness, and ensuring positive definiteness.
Overall, the C_AO_CMAES class encapsulates the logic of the CMA-ES algorithm, providing an interface for initialization, parameter tuning, and performing optimization. It contains the parameters, data arrays, and methods needed to implement the algorithm. The separation of public and private fields provides access control to the internal components of the class.
//———————————————————————————————————————————————————————————————————— 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/en/articles/18227"; // Default parameters popSize = 50; // Default population size (lambda) mu = 25; // Number of parents (half of the population) learningRateC1 = 0.01; // Learning rate for rank-1 update learningRateCMu = 0.8; // Learning rate for rank-μ update stepSizeDamping = 0.6; // Damping for step size // Create and initialize the parameters array 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 [], // minimum values const double &rangeMaxP [], // maximum values const double &rangeStepP [], // step size const int epochsP = 0); void Moving (); void Revision (); //------------------------------------------------------------------ private: //--------------------------------------------------------- // CMA-ES specific parameters int mu; // Number of parents (selected points) double sigma; // Step size double learningRateC1; // Learning rate for rank-1 update double learningRateCMu; // Learning rate for rank-μ update double stepSizeDamping; // Damping factor for step size update // CMA-ES specific data structures double weights []; // Recombination weights double covMatrix []; // Covariance matrix (stored as a one-dimensional array) double B []; // C eigenvectors double D []; // C eigenvalues (square roots) double pc []; // Evolution path for C double ps []; // Evolution path for step-size control double mu_eff; // Effective selection mass int counteval; // Counter of function evaluations since the last decomposition int eigeneval; // Generation counter when decomposition was performed double chiN; // Expected norm N(0,I) int eigenInterval; // Interval for eigenvalue decomposition // Variables for performance optimization double cs; // Learning rate for the sigma path double cc; // Learning rate for rank-1 path double damps; // Damping for sigma double hsig_threshold; // Threshold for the Heaviside function // Auxiliary arrays double y_vec []; // Mutation vector double arindex []; // Array of indices for sorting double arfitness []; // Fitness array for sorting double temp_vec []; // Temporary vector for matrix operations double invsqrtC_times_yw []; // Temporary storage for C^(-1/2) * y_w // Caching variables for Box-Muller double cached_normal; bool has_cached; // Auxiliary methods void InitDistribution (); void UpdateDistribution (); void ComputeEigendecomposition (); double GetChiN (); void SortPopulation (); void ComputeWeights (); void UpdateMean (); bool CheckPositiveDefinite (); void EnforcePositiveDefinite (); }; //————————————————————————————————————————————————————————————————————
The Init method of the C_AO_CMAES class initializes the CMA-ES algorithm. It takes the minimum and maximum parameter values, step size and number of epochs as input. First, the StandardInit method is called to perform standard initialization. The has_cached flag is initialized to 'false' to indicate that there is no cached value for Box-Muller random number generation. Next, the sigma step size is initialized to the initial value of 0.3 (30% of the search range). The ComputeWeights method is called to compute the recombination weights, and the GetChiN method is called to compute the expected norm N(0, I).
cs, cc, damps and hsig_threshold parameters are calculated based on mu_eff (effective mass) and the number of coordinates (coords). These parameters are used to adapt the CMA-ES strategy during operation. The eigenInterval is calculated and set to determine how often the eigenvalue decomposition is performed. This parameter is configured to optimize performance.
Memory is allocated for CMA-ES-specific arrays such as the covMatrix covariance matrix, B eigenvectors, D eigenvalues, pc and ps evolution paths, y_vec, arindex, arfitness, temp_vec and invsqrtC_times_yw time vectors. Arrays are recreated only if the problem dimension needs to be changed.
The pc and ps evolution path arrays are initialized to zeros, the covMatrix covariance matrix and the B eigenvector matrix are initialized to the identity matrix, and the D eigenvalue array is initialized to ones. The InitDistribution method is called to initialize the initial distribution. Then the evaluation counters of the counteval and eigeneval functions are reset. The "revision" flag is set to 'true' to indicate that the fitness values should be recalculated. 'true' is returned on successful initialization.
In conclusion, the Init method performs all necessary initialization for the CMA-ES algorithm to work, including setting parameters, allocating memory for arrays, and initializing initial values.
//———————————————————————————————————————————————————————————————————— bool C_AO_CMAES::Init (const double &rangeMinP [], // minimum values const double &rangeMaxP [], // maximum values const double &rangeStepP [], // step size const int epochsP = 0) // number of epochs { if (!StandardInit (rangeMinP, rangeMaxP, rangeStepP)) return false; //------------------------------------------------------------------ // Initialize Box-Muller caching has_cached = false; // Initialize CMA-ES specific variables sigma = 0.3; // Initial step size (30% of search range) // Calculate the effective mass of the variance selection ComputeWeights (); // Expected norm N(0,I) chiN = GetChiN (); // Calculate and save strategy parameters 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); // Set the eigenvalue decomposition interval - tuning for performance eigenInterval = (int)(coords / (10.0 * MathSqrt (learningRateC1 + learningRateCMu))); eigenInterval = MathMax (1, eigenInterval); // Allocate arrays only once 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); // Initialize evolution paths with zeros ArrayInitialize (pc, 0); ArrayInitialize (ps, 0); // Fast initialization of the identity covariance matrix and decomposition 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; } // Initialize initial distribution InitDistribution (); // Reset calculation counters counteval = 0; eigeneval = 0; // Forced fitness recalculation revision = true; return true; } //————————————————————————————————————————————————————————————————————
The Moving method in the C_AO_CMAES class is responsible for generating new descendants (solutions) in the population. Transforming y = B * D * z: inside the loop over the population there is a nested loop that iterates over all the coordinates of each individual. At each i coordinate the y_vec [i] value is calculated. This is done by matrix multiplication. In fact, it is a multiplication of a vector of z random numbers by the matrix "B * D" matrix, where B are the eigenvectors of the covariance matrix, D are the eigenvalues of the covariance matrix, and z are random numbers generated using the u.PowerDistribution function.
Descendant generation: the (a [k].c [i]) value is calculated for each i coordinate of each k descendant. This is done by adding the scaled vector y_vec[i] to the current mean cB[i] (the center of the current distribution). Scaling is done by multiplying y_vec[i] by 'sigma' step size. Thus, a [k].c [i] = cB [i] + sigma * y_vec [i].
The u.SeInDiSp function ensures that the coordinate values of the offspring remain within the valid ranges defined by rangeMin[i], rangeMax[i] and rangeStep[i]. It can clip values, bounce them off boundaries, and quantize them to the nearest valid step value. After all descendants are generated, the "revision" flag is set to 'true', and the fitness function (quality) values for the generated descendants should be recalculated.
As a result, the method generates a new population of solutions based on the current distribution (cB average and the covariance matrix, represented by B and D) and the sigma step size. It uses random numbers to create small variations around the current mean to explore the search space. The SeInDiSp function ensures that new solutions remain within the admissible bounds, and the "revision" flag ensures that their fitness will then be evaluated.
//———————————————————————————————————————————————————————————————————— void C_AO_CMAES::Moving () { // Generate new lambda descendants for (int k = 0; k < popSize; k++) { // Apply the transformation 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); } // Generate a descendant: 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]); } } // Fitness recalculation mark revision = true; } //————————————————————————————————————————————————————————————————————
The Revision method of the C_AO_CMAES class is responsible for updating the algorithm parameters after evaluating the fitness function for a new population. First, the method checks the value of the "revision" flag. Then, the SortPopulation() method is called, which sorts the current population so that the best individuals are placed at the top. The UpdateDistribution() method updates the parameters of the CMA-ES distribution (the cB mean of the strategy vector, the covariance matrix, the sigma step size, the evolutionary paths) based on information about the best individuals in the population. Distribution update is a key step of CMA-ES, as it allows the algorithm to adapt its search strategy based on the obtained results.
As a result, the Revision method represents the main adaptation cycle of CMA-ES. It sorts the population, updates the distribution parameters based on the best individuals, and increments the computation counter.
//———————————————————————————————————————————————————————————————————— void C_AO_CMAES::Revision () { if (!revision) return; revision = false; // Sort the population by fitness SortPopulation (); // Update distribution parameters based on selected individuals UpdateDistribution (); // Update the calculation counter counteval++; } //————————————————————————————————————————————————————————————————————
The InitDistribution method of the C_AO_CMAES class is intended to initialize the initial distribution of the solution search. The method determines an initial mean value (cB) for each coordinate in the search space. For each i coordinate, the value of cB[i] is set randomly in the range between rangeMin[i] and rangeMax[i] using the u.RNDfromCI() function, with the initial mean value placed in the center of the valid range for each variable. The outer loop iterates over the entire population (popSize individuals). Inside this loop is an inner loop that iterates over all the coordinates of each individual. Coordinates are generated for each descendant. First, they are randomly selected from the given range RNDfromCI and then rounded to the sampling step using SeInDiSp.
The u.RNDfromCI function generates a uniformly distributed random number in the given interval. The u.SeInDiSp function ensures that the generated coordinates for each individual lie within the acceptable limits defined by rangeMin [i], rangeMax [i] and rangeStep [i].
The method initializes the initial population using random solutions uniformly distributed in the search space. The center of distribution (cB) is also chosen randomly. This provides a starting point for subsequent iterations of CMA-ES, in which the algorithm will adapt the search strategy to find the optimal solution.
//+------------------------------------------------------------------+ //| Initialize search distribution | //+------------------------------------------------------------------+ void C_AO_CMAES::InitDistribution () { // Set the initial mean to the center of the search space 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++) { // Generate a uniformly distributed point 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]); } } } //————————————————————————————————————————————————————————————————————
The GetChiN method computes the expected norm of a random vector drawn from N (0,I) for a given number of coordinates (coords). The result is approximated based on the dimension of the search space.
//+------------------------------------------------------------------+ //| Calculate the expected norm 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)); } //————————————————————————————————————————————————————————————————————
The SortPopulation is designed to sort a population of solutions by their objective function values (fitness) in order to determine the best solution. First, the method copies the fitness values of each individual from the population into the arfitness auxiliary array. The arindex array is also created, which initially contains the indices of the individuals. This is necessary so that after sorting by fitness it is possible to restore the correspondence between the sorted fitness and the original individual in the population.
The insertion sort algorithm is used. It iterates through the arfitness array, and for each element, inserts it into the correct place in the sorted part of the array, shifting larger elements to the right. It is important to note that the sorting is performed in descending order ( arfitness [j] < tempFitness in the 'while' loop), the goal is to maximize the target function (fitness). Along with arfitness, arindex is also sorted synchronously to track the position of the original indices of individuals.
After sorting, the method checks whether the fitness of the best individual is better than the current best solution, and if so, fB is updated and the best solution is replaced by the coordinates of the best individual from the population. Thus, the method not only sorts the population by fitness, but also keeps track of the best solution found at the moment.
//+------------------------------------------------------------------+ //| Sort the population by fitness | //+------------------------------------------------------------------+ void C_AO_CMAES::SortPopulation () { // Copy fitness values and indices 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; // Sort in descending order (for maximization) 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; } // Update the best solution if necessary if (arfitness [0] > fB) { fB = arfitness [0]; int bestIdx = (int)arindex [0]; ArrayCopy (cB, a [bestIdx].c, 0, 0, coords); } } //————————————————————————————————————————————————————————————————————
The UpdateMean method updates the population mean (cB) by weighted recombination of the best individuals. For each coordinate, a new average is calculated as the sum of the weighted coordinate values from the best individuals, where the weights are specified in the weights array.
//+------------------------------------------------------------------+ //| Update the mean using weighted recombination | //+------------------------------------------------------------------+ void C_AO_CMAES::UpdateMean () { // Weighted recombination: 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; } } //————————————————————————————————————————————————————————————————————
The ComputeWeights method calculates weights for weighted recombination, allocates the weights array by size mu, calculates log (mu + 0.5) for use in the weights formula. For each i from 0 to mu-1, it assigns a weight as the difference between log(mu + 0.5) and log(i+1), summing these weights, the method divides all the weights by the sum so that the sum is equal to 1, and calculates the sum of the squares of the weights, calculates mu_eff - the efficiency of using the weights, which is important for tuning the speed and variance in the algorithm.
This method provides optimal weights for recombination, taking into account the contribution of the best solutions.
//+------------------------------------------------------------------+ //| Calculate weighted recombination weights | //+------------------------------------------------------------------+ void C_AO_CMAES::ComputeWeights () { // Allocate the array of weights ArrayResize (weights, mu); // Pre-calculate log(mu + 0.5) double log_mu_plus_half = MathLog (mu + 0.5); // Calculate positive weights double sum = 0.0; for (int i = 0; i < mu; i++) { weights [i] = log_mu_plus_half - MathLog (i + 1); sum += weights [i]; } // Normalize weights 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]; } // Calculate the effective mass of the variance selection mu_eff = sum_weights * sum_weights / sum_squares; } //————————————————————————————————————————————————————————————————————
The UpdateDistribution method updates the distribution parameters, namely the covariance matrix (covMatrix) and the step size (sigma) in the CMA-ES algorithm. If enough generations have passed (counteval - eigeneval > eigenInterval), then the eigenvalue decomposition is calculated and the current mean (cB) is stored in the oldMean temporary array. The method calls UpdateMean to calculate a new mean, calculates the difference between the new and old mean divided by sigma, and also performs the multiplication of y_w by C^(-1/2), broken into several steps using the B and D decomposition, and stores the result in invsqrtC_times_yw.
The method then updates the evolution path for ps step size using invsqrtC_times_yw, determines whether progress is considered "stuck" based on the ps norm and the expected length, then updates the evolution path for the pc covariance matrix using y_w and the hsig progress indicator. The method first calculates the c1a rank-1 update matrix, then the cmu rank-μ update matrix based on the difference between individuals and the old mean divided by sigma using 'weights', updates the covMatrix covariance matrix using the rank-1 and rank-μ updates, and adjusts c1 when progress stalls. The EnforcePositiveDefinite is periodically called to ensure that the covariance matrix is positive definite. The sigma step size is updated based on the ps norm and is constrained between min_sigma and max_sigma for numerical stability.
Generally, UpdateDistribution adjusts the distribution parameters (covariance matrix and step size) based on the search history (evolutionary path) and population data to adapt to the landscape of the objective function.
//+------------------------------------------------------------------+ //| Update distribution parameters | //+------------------------------------------------------------------+ void C_AO_CMAES::UpdateDistribution () { // Check the necessity of eigenvalue decomposition if (counteval - eigeneval > eigenInterval) { ComputeEigendecomposition (); eigeneval = counteval; } // Preserve the old average double oldMean []; ArrayResize (oldMean, coords); ArrayCopy (oldMean, cB, 0, 0, coords); // Update average UpdateMean (); // Calculate the mean shift double y_w []; ArrayResize (y_w, coords); for (int j = 0; j < coords; j++) { y_w [j] = (cB [j] - oldMean [j]) / sigma; } // Calculate C^(-1/2) * y_w // Step 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]; } } // Step 2: D^(-1) * (B^T * y_w) for (int i = 0; i < coords; i++) { temp_vec [i] /= D [i]; } // Step 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]; } } // Update the evolution path for 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]; } // Heaviside function 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; // Update the evolution path for 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]; } // Prepare rank-1 update 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]; } } // Prepare rank-μ update double cmu []; ArrayResize (cmu, coords * coords); ArrayInitialize (cmu, 0.0); for (int k = 0; k < mu; k++) { int idx = (int)arindex [k]; // Calculate y_i = (x_i - m_old) / sigma for (int i = 0; i < coords; i++) { temp_vec [i] = (a [idx].c [i] - oldMean [i]) / sigma; } // Add the weighted outer product 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; } } } // Update C the covariance matrix double c1 = learningRateC1; double cmu_rate = learningRateCMu; // Adjust c1 if hsig is 'false' (progress stalled) if (!hsig) { c1 *= (1.0 - (1.0 - delta_hsig) * cc * (2.0 - cc)); } double one_minus_c1_cmu = 1.0 - c1 - cmu_rate; // Update C with rank-1 and rank-μ updates 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]; // Maintain symmetry if (i != j) { covMatrix [j * coords + i] = covMatrix [i * coords + j]; } } } // Ensure positive definiteness if (counteval % (10 * eigenInterval) == 0) { EnforcePositiveDefinite (); } // Update the sigma step size double exponent = (cs / damps) * (norm_ps / chiN - 1.0); sigma *= MathExp (exponent); // Limit sigma for numerical stability double min_sigma = 1e-16; double max_eigenvalue = D [0]; // D sorted in descending order 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; } //————————————————————————————————————————————————————————————————————
The ComputeEigendecomposition " performs the decomposition of the covMatrix symmetric covariance into eigenvalues and eigenvectors using the improved Jacobi method, the method creates a copy of the covMatrix so as not to modify the original data. Initially, B is set as the identity matrix representing the initial basis. It performs maximum iterations or until non-zero elements on the off-diagonal are less than "tolerance".
Finding the maximum off-diagonal element: Find the element with the largest absolute value off the diagonal to select for rotation. Next, if the maximum is less than "tolerance", the cycle is interrupted. The method calculates phi to perform the Jacobian rotation, modifies the elements of C_copy by applying the rotation, and updates the columns of B matrix (eigenvectors) according to the rotation. After the iterations, the method finds the square root of the diagonal elements (ensuring minimum positivity of 1e-14) and orders the eigenvalues and corresponding eigenvectors in descending order for further use. Thus, the method allows for accurate calculation of the eigenvalues and eigenvectors of the covariance matrix.
//+------------------------------------------------------------------+ //| Calculate the eigenvalue decomposition using the Jacobi method | //+------------------------------------------------------------------+ void C_AO_CMAES::ComputeEigendecomposition () { // Create a copy of the covariance matrix for decomposition double C_copy []; ArrayResize (C_copy, coords * coords); ArrayCopy (C_copy, covMatrix); // Initialize B as the identity matrix for (int i = 0; i < coords; i++) { for (int j = 0; j < coords; j++) { B [i * coords + j] = (i == j) ? 1.0 : 0.0; } } // Improved Jacobi eigenvalue decomposition int max_iterations = 10; //50 * coords; double tolerance = 0.01; //1e-14 * coords * coords; for (int iter = 0; iter < max_iterations; iter++) { // Find the largest off-diagonal element 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; } } } // Check convergence if (max_val < tolerance) break; // Calculate the rotation angle 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); // Update the matrix elements 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; // Update other elements in p and q rows/columns 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; } } // Update eigenvectors 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; } } // Extract eigenvalues and ensure positivity double min_eigenvalue = 1e-14; for (int i = 0; i < coords; i++) { D [i] = MathSqrt (MathMax (min_eigenvalue, C_copy [i * coords + i])); } // Sort eigenvalues and vectors in descending order 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) { // Exchange eigenvalues double temp = D [i]; D [i] = D [max_idx]; D [max_idx] = temp; // Exchange eigenvectors 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; } } } } //————————————————————————————————————————————————————————————————————
The CheckPositiveDefinite method checks whether the covariance matrix is positive definite. It quickly checks the positivity of diagonal elements and compares the minimum eigenvalue (from the D array sorted in descending order) with a small positive number 1e-14. If both checks pass, the method returns 'true'.
//+------------------------------------------------------------------+ //| Check the positive definiteness of the covariance matrix | //+------------------------------------------------------------------+ bool C_AO_CMAES::CheckPositiveDefinite () { // Quick check: all diagonal elements should be positive for (int i = 0; i < coords; i++) { if (covMatrix [i * coords + i] <= 0) return false; } // Check if the minimum eigenvalue is positive double min_eigenvalue = D [coords - 1]; // D is sorted in descending order return min_eigenvalue > 1e-14; } //————————————————————————————————————————————————————————————————————
The EnforcePositiveDefinite method ensures positive definiteness of the covMatrix covariance matrix . It performs several steps: it finds the minimum element on the diagonal and adds the missing "correction" value to the diagonal elements so that the minimum is equal to 1e-10, it maximizes the symmetry of the matrix by averaging each off-diagonal element with its transpose.
If the matrix is still not positive definite, an eigenvalue decomposition is performed using ComputeEigendecomposition, replacing all eigenvalues less than √1e-10 with √1e-10 to ensure positivity. The method then recalculates covMatrix using the B eigenvectors and the corrected D eigenvalues to obtain a corrected positive definite matrix. This approach ensures that the covariance matrix becomes positive definite and suitable for further calculations in the CMA-ES algorithm.
//+------------------------------------------------------------------+ //| Ensuring positive definiteness of the covariance matrix | //+------------------------------------------------------------------+ void C_AO_CMAES::EnforcePositiveDefinite () { // Method 1: Adding a small value to the diagonal double min_diag = 1e308; // A very large number 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; } } // Method 2: Ensuring symmetry 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 still not positive definite, perform the factorization and fix 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); } } // Reconstruct 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; } } } } //————————————————————————————————————————————————————————————————————
Test results
Now we can finally look at the results. The algorithm, as can be immediately noted, copes well and quickly with medium-sized problems, and with certain settings, with low-dimensional problems, but the most difficult, multidimensional problems require a disproportionate amount of time to execute, so they were effectively excluded from the results. Why does this happen? Traditionally, matrix calculations in MQL5 present a serious computational problem. When implementing matrix operations manually, the algorithm demonstrates unsatisfactory performance on high-dimensional problems, which significantly limits its practical applicability. However, with the advent of built-in classes for working with matrices, the situation changes dramatically. It is now critical to use the platform's native implementation of matrix operations for all computationally intensive tasks.
=============================
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%)
In the visualization, we can see that the algorithm works both with "large jumps" and focuses on the extrema of the function, exploring them thoroughly.
CMA-ES on the Hilly test function
CMA-ES on the Forest test function
CMA-ES on the Megacity test function
Based on the results, the CMA-ES algorithm ranks 38th in the ranking of the strongest population-based optimization algorithms.
# | AO | Description | Hilly | Hilly Final | Forest | Forest Final | Megacity (discrete) | Megacity Final | Final Result | % of 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 optimization | 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 | |
Summary
CMA-ES's ability to adapt the search strategy to the local geometry of the objective function makes it virtually indispensable for complex, ill-conditioned problems with unknown structure. The heavy tails of the power-law distribution allow the algorithm to make "long jumps," which is critical for escaping local optima. Its computational complexity of O(n²) memory and O(n³) time severely limits the algorithm's applicability. For high-dimensional problems (n > 100), the resource intensity becomes disproportionate to the benefits obtained. The running time on multidimensional functions grows so quickly that it makes the algorithm practically inapplicable for large n.
CMA-ES works on noisy functions, copes with discontinuous objective functions, is effective on multi-extreme landscapes, and the secret to this versatility lies in the algorithm's fundamental philosophy: instead of making assumptions about the function's structure, it carefully explores it, adapting its search strategy based on accumulated experience. CMA-ES changes the landscape of evolutionary computation by showing that biologically inspired algorithms can have rigorous mathematical foundations and that algorithm adaptation is not simply a parameter tweak, but a fundamental principle of learning the structure of a problem.
The algorithm is great for medium-sized problems, it is a specialized tool that looks like a worthy alternative in its niche. The algorithm combines mathematical elegance with practical efficiency, but requires careful application given its computational limitations. In such cases, fast MQL5 matrix operations have recently been introduced by the MQL5 language developers to calculate matrices; in this implementation of the algorithm, conventional sequential calculations were used.
Figure 2. Color gradation of algorithms according to the corresponding tests
Figure 3. Histogram of algorithm testing results (scale from 0 to 100, the higher the better, where 100 is the maximum possible theoretical result, in the archive there is a script for calculating the rating table)
CMAES pros and cons:
Pros:
- Good convergence on functions of medium dimension.
Cons:
- Scatter of results on low-dimensional functions.
- Critical resource intensity on large-scale functions.
The article is accompanied by an archive with the current versions of the algorithm codes. The author of the article is not responsible for the absolute accuracy in the description of canonical algorithms. Changes have been made to many of them to improve search capabilities. The conclusions and judgments presented in the articles are based on the results of the experiments.
Programs used in the article
| # | Name | Type | Description |
|---|---|---|---|
| 1 | #C_AO.mqh | Include | Parent class of population optimization algorithms |
| 2 | #C_AO_enum.mqh | Include | Enumeration of population optimization algorithms |
| 3 | TestFunctions.mqh | Include | Library of test functions |
| 4 | TestStandFunctions.mqh | Include | Test stand function library |
| 5 | Utilities.mqh | Include | Library of auxiliary functions |
| 6 | CalculationTestResults.mqh | Include | Script for calculating results in the comparison table |
| 7 | Testing AOs.mq5 | Script | The unified test stand for all population optimization algorithms |
| 8 | Simple use of population optimization algorithms.mq5 | Script | A simple example of using population optimization algorithms without visualization |
| 9 | Test_AO_CMAES.mq5 | Script | CMAES test bench |
Translated from Russian by MetaQuotes Ltd.
Original article: https://www.mql5.com/ru/articles/18227
Warning: All rights to these materials are reserved by MetaQuotes Ltd. Copying or reprinting of these materials in whole or in part is prohibited.
This article was written by a user of the site and reflects their personal views. MetaQuotes Ltd is not responsible for the accuracy of the information presented, nor for any consequences resulting from the use of the solutions, strategies or recommendations described.
Meta-Labeling the Classics (Part 1): Filtering and Sizing RSI Trades
Building a Dynamic STF Liquidity Sweep Indicator in MQL5
Joint Recurrence Quantification Analysis (JRQA) in MQL5: Detecting Simultaneous Recurrence in Two Series
Integrating AI into 3 Smart Money Concepts (SMC): OB, BOS, and FVG
- Free trading apps
- Over 8,000 signals for copying
- Economic news for exploring financial markets
You agree to website policy and terms of use