Particle swarm CMA-ES Evolution strategy

Table of Contents

• Black box optimization
  • Introduction
  • Inclusions
  • Parameters

Black box optimization

In this example we show how to code PS-CMA-ES. This is just a simple variation to the CMA-ES, where you have multiple CMA-ES running. The the best solution across them is used to produce a drift velocity toward that point.

Simulation video 1

Introduction

In this example we try to find the global optimum of a function. In particular we are using the function F15 from the CEC 2005 benchmark test, to validate that PS-CMA-ES work. This example contain multiple files:

• f15_cec_const.hpp definitions of constants for the F15 function
• f15_cec_fun.hpp the function itself

The function is quite complicated and for reference please refere to the function F15 "Hybrid Composition" in the CEC 2005 test. The function can be called with hybrid_composition<dim>(x) where dim is the dimensionality and x is the point where is evaluated the function. The dimensionality can go from 1 to 50.

Considering to have a function \( f \) from \( \mathbb{R}^{dim} \) to \( \mathbb{R} \), the algorithm use a set of particles to find in parallel the global optimum of a function. The algorithm rather than try to find the global optimum sampling point randomly in the space, it uses a set of particles each of them having a gaussian sampling distribution \( e^{\sigma \cdot x^tCx} \) with C a \( dim \cdot dim \) matrix. At each step for each particle p lambda points are sampled using the sampling distribution centered on the particle position. The covariant matrix and sigma are is subsequently adjusted to favor sampling around the best sampled points. In order to do this the algorithm need the eigen-value decomposition of \( C = B^{t}DB \) where \( D \) is a diagonal Matrix and \( B \) is the Matrix of the Eigen-vector. In order to reduce or increase the sampling area the sigma is instead used. The algorithm use the vector path_s to detect stagnation of the particle movement, and use path_c(a transfomed version of path_s) to refine the sampling covariant matrix from the fact that the particle is "moving" toward that direction. PS-CMA-ES is just a variation in which every N_pos CMA-Es steps the CMA-ES is sampling distribution and position is biased toward the best founded point across all independent CMA-ES.

Explain the CMA-ES algorithm in detail is out of the purpose of this tutorial example. We will briefly go across the main step. For a full reference of the CMA-ES algoritm please refers to this paper. While for PS-CMA-ES refers to this paper.

Inclusions

In this example we use a set of particles so we will use vector_dist, we will use Eigen dense matrix. Because the standard dense matrix are not compatible with the vector_dist we will use EMatrix that are simple wrapper to Eigen::Matrix but compatible with vector_dist. Because EMatrix are compatible with all the Eigen value functions we can use them in all Eigen functions. For CMA-ES algorithm we also need Eigen-value eigen-vector decomposition and Jacobi-Rotation for the particle-swarm part.

 
// This example need EIGEN, if we do not have NOP the full example
 
#include "config.h"
 
#ifndef HAVE_EIGEN
 
int main(int argc, char* argv[])
{
    return 0;
}
 
#else
 
 
#define EIGEN_USE_LAPACKE
#include "Vector/vector_dist.hpp"
#include "DMatrix/EMatrix.hpp"
#include <Eigen/Eigenvalues>
#include <Eigen/Jacobi>
#include <limits>
#include "Vector/vector_dist.hpp"
#include <f15_cec_fun.hpp>
#include <boost/math/special_functions/sign.hpp>
 

PS-CMA-ES require several properties to be stored on the particles, some has been already explained. Here we explain the others.

• Zeta contain the lambda sampled points (before apply the covariant matrix tranformation) so it contain points samples on a gaussian of sigma=1 centered in zero
• ord Contain the sequrnce if we want the lambda generated points in order from the best to the worst
• stop If a flag that indicate that the CMA-ES reached some stop criteria
• fithist It contain historical information about the particles to penalize them in case the go out of boundary. It is 1:1 taken from cmaes.m (production version) this paper (or Google it)
• weight Same concept of fithist other information to penalize particles going out of the boundary
• validfit Same concept of fithist other information to penalize particles going out of the boundary
• xold It contain the previous position of the particles used in several calculations
• last_restart CMA-ES The CMA-ES sigma become very small the CMA-ES converged. At this point we can do two things, one is to stop the CMA-ES, the other is to restart-it to explore better the space. In case it restart. this parameter indicate at which iteration happen the last restart.
• iniphase Same concept of fithist other information to penalize particles going out of the boundary
• xmean_st This contain the new position of the particle it will be stored as particle position at the end of the CMA-ES step
• xmean_st This contain the new position of the particle in a space where we do not apply the covariant transformation. (In practice is a weighted sum of the Zeta samples points)

 
 
constexpr int sigma = 0;
constexpr int Cov_m = 1;
constexpr int B = 2;
constexpr int D = 3;
constexpr int Zeta = 4;
constexpr int path_s = 5;
constexpr int path_c = 6;
constexpr int ord = 7;
constexpr int stop = 8;
constexpr int fithist = 9;
constexpr int weight = 10;
constexpr int validfit = 11;
constexpr int xold = 12;
constexpr int last_restart = 13;
constexpr int iniphase = 14;
constexpr int xmean_st = 15;
constexpr int meanz_st = 16;
 
typedef vector_dist<dim,double, aggregate<double,
                                         EMatrixXd,
                                         EMatrixXd,
                                         Eigen::DiagonalMatrix<double,Eigen::Dynamic>,
                                         EVectorXd[lambda],
                                         EVectorXd,
                                         EVectorXd,
                                         int[lambda],
                                         int,
                                         double [hist_size],
                                         double [dim],
                                         double,
                                         EVectorXd,
                                         int,
                                         bool,
                                         EVectorXd,
                                         EVectorXd> > particle_type;
 

Parameters

CMA-ES and further PS-CMA-ES require some parameters in order to work. refers to the papers above to have a full explanation, but here is a short one

• dim Dimensionality of the test function
• lambda number of sample points taken at each iteration by CMA-ES suggested to use \( 4+floor(3*log(dim)) \)
• mu only mu best points are considered to adapt the Covariant matrix
• psoWeight How much the pso step bias the particle positions
• N_pso Number of CMA-ES step before do a PSO step (200 give the possibility to the CMA-ES to explore the neighborhood and narrow down at least a funnel)
• stopTolX stop criteria for CMA-ES. When the the sampling area is small enough stop
• StopToUpX stop criteria is the sampling area become too big
• restart_cma If the CMA-ES reach a stop criteria reinitialize and restart
• hist_size size of the array fit_hist (default should be mainly fine)