CMA-ES

CMA-ES stands for Covariance Matrix Adaptation Evolution Strategy. Evolution strategies (ES) are stochastic, derivative-free methods for numerical optimization of non-linear or non-convex continuous optimization problems. They belong to the class of evolutionary algorithms and evolutionary computation. An evolutionary algorithm is broadly based on the principle of biological evolution, namely the repeated interplay of variation (via recombination and mutation) and selection: in each generation (iteration) new individuals (candidate solutions, denoted as x) are generated by variation, usually in a stochastic way, of the current parental individuals. Then, some individuals are selected to become the parents in the next generation based on their fitness or objective function value f(x). Like this, over the generation sequence, individuals with better and better f-values are generated.

In an evolution strategy, new candidate solutions are sampled according to a multivariate normal distribution in the \mathbb{R}^n. Recombination amounts to selecting a new mean value for the distribution. Mutation amounts to adding a random vector, a perturbation with zero mean. Pairwise dependencies between the variables in the distribution are represented by a covariance matrix. The covariance matrix adaptation (CMA) is a method to update the covariance matrix of this distribution. This is particularly useful, if the function f is ill-conditioned.

Adaptation of the covariance matrix amounts to learning a second order model of the underlying objective function similar to the approximation of the inverse Hessian matrix in the Quasi-Newton method in classical optimization. In contrast to most classical methods, fewer assumptions on the nature of the underlying objective function are made. Only the ranking between candidate solutions is exploited for learning the sample distribution and neither derivatives nor even the function values themselves are required by the method.

Principles

Illustration of an actual optimization run with covariance matrix adaptation on a simple two-dimensional problem. The spherical optimization landscape is depicted with solid lines of equal f-values. The population (dots) is much larger than necessary, but clearly shows how the distribution of the population (dotted line) changes during the optimization. On this simple problem, the population concentrates over the global optimum within a few generations.

Two main principles for the adaptation of parameters of the search distribution are exploited in the CMA-ES algorithm.

First, a maximum-likelihood principle, based on the idea to increase the probability of successful candidate solutions and search steps. The mean of the distribution is updated such that the likelihood of previously successful candidate solutions is maximized. The covariance matrix of the distribution is updated (incrementally) such that the likelihood of previously successful search steps is increased. Both updates can be interpreted as a natural gradient descent. Also, in consequence, the CMA conducts an iterated principal components analysis of successful search steps while retaining all principal axes. Estimation of distribution algorithms and the Cross-Entropy Method are based on very similar ideas, but estimate (non-incrementally) the covariance matrix by maximizing the likelihood of successful solution points instead of successful search steps.

Second, two paths of the time evolution of the distribution mean of the strategy are recorded, called search or evolution paths. These paths contain significant information about the correlation between consecutive steps. Specifically, if consecutive steps are taken in a similar direction, the evolution paths become long. The evolution paths are exploited in two ways. One path is used for the covariance matrix adaptation procedure in place of single successful search steps and facilitates a possibly much faster variance increase of favorable directions. The other path is used to conduct an additional step-size control. This step-size control aims to make consecutive movements of the distribution mean orthogonal in expectation. The step-size control effectively prevents premature convergence yet allowing fast convergence to an optimum.

Algorithm

In the following the most commonly used (μ/μw, λ)-CMA-ES is outlined, where in each iteration step a weighted combination of the μ best out of λ new candidate solutions is used to update the distribution parameters. The main loop consists of three main parts: 1) sampling of new solutions, 2) re-ordering of the sampled solutions based on their fitness, 3) update of the internal state variables based on the re-ordered samples. A pseudocode of the algorithm looks as follows.

 set \lambda  // number of samples per iteration, at least two, generally > 4
 initialize m, \sigma, C=I, p_\sigma=0, p_c=0  // initialize state variables
 while not terminate  // iterate
    for i in \{1...\lambda\}  // sample \lambda new solutions and evaluate them
       x_i = sample_multivariate_normal(mean=m, covariance_matrix=\sigma^2 C )
       f_i = fitness(x_i)
    x_{1...\lambda}  x_{s(1)...s(\lambda)} with s(i) = argsort(f_{1...\lambda}, i)  // sort solutions
    m' = m  // we need later m - m' and x_i - m'       
    m  update_m(x_1, ... , x_\lambda)  // move mean to better solutions 
    p_\sigma  update_ps(p_\sigma, \sigma^{-1} C^{-1/2} (m - m'))  // update isotropic evolution path
    p_c  update_pc(p_c, \sigma^{-1}(m - m'), ||p_\sigma||)  // update anisotropic evolution path
    C  update_C(C, p_c, {(x_1 - m')}/{\sigma},... , {(x_\lambda - m')}/{\sigma})  // update covariance matrix
    \sigma  update_sigma(\sigma, ||p_\sigma||)  // update step-size using isotropic path length
 return m or x_1

The order of the five update assignments is relevant. In the following, the update equations for the five state variables are specified.

Given are the search space dimension n and the iteration step k. The five state variables are

m_k\in\mathbb{R}^n, the distribution mean and current favorite solution to the optimization problem,
\sigma_k>0, the step-size,
 C_k, a symmetric and positive definite n\times n covariance matrix with  C_0 = I and
 p_\sigma\in\mathbb{R}^n, p_c\in\mathbb{R}^n, two evolution paths, initially set to the zero vector.

The iteration starts with sampling \lambda>1 candidate solutions x_i\in\mathbb{R}^n from a multivariate normal distribution \textstyle \mathcal{N}(m_k,\sigma_k^2 C_k), i.e. for i=1,...,\lambda


\begin{align}
  x_i \ &\sim\ \mathcal{N}(m_k,\sigma_k^2 C_k)      
     \\&\sim\ m_k + \sigma_k\times\mathcal{N}(0,C_k) 
\end{align}

The second line suggests the interpretation as perturbation (mutation) of the current favorite solution vector m_k (the distribution mean vector). The candidate solutions  x_i are evaluated on the objective function f:\mathbb{R}^n\to\mathbb{R} to be minimized. Denoting the f-sorted candidate solutions as


    \{x_{i:\lambda}\;|\;i=1\dots\lambda\} = \{x_i\;|\;i=1\dots\lambda\} \;\;\text{and}\;\; 
     f(x_{1:\lambda})\le\dots\le f(x_{\mu:\lambda})\le f(x_{\mu+1:\lambda}) \dots,

the new mean value is computed as


\begin{align}
  m_{k+1} &= \sum_{i=1}^{\mu} w_i\, x_{i:\lambda} 
  \\ &= m_k + \sum_{i=1}^{\mu} w_i\, (x_{i:\lambda} - m_k) 
\end{align}

where the positive (recombination) weights  w_1 \ge w_2 \ge \dots \ge w_\mu > 0 sum to one. Typically, \mu \le \lambda/2 and the weights are chosen such that \textstyle \mu_w := 1 / \sum_{i=1}^\mu w_i^2 \approx \lambda/4. The only feedback used from the objective function here and in the following is an ordering of the sampled candidate solutions due to the indices i:\lambda.

The step-size \sigma_k is updated using cumulative step-size adaptation (CSA), sometimes also denoted as path length control. The evolution path (or search path) p_\sigma is updated first.


  p_\sigma \gets \underbrace{(1-c_\sigma)}_{\!\!\!\!\!\text{discount factor}\!\!\!\!\!}\, p_\sigma 
    + \overbrace{\sqrt{1 - (1-c_\sigma)^2}}^{
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\text{complements for discounted variance}
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!} \underbrace{\sqrt{\mu_w} 
     \,C_k^{\;-1/2} \, \frac{\overbrace{m_{k+1} - m_k}^{\!\!\!\text{displacement of}\; m\!\!\!}}{\sigma_k}}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
                      \text{distributed as}\; \mathcal{N}(0,I)\;\text{under neutral selection}
                      \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}

  \sigma_{k+1} = \sigma_k \times \exp\bigg(\frac{c_\sigma}{d_\sigma}
                          \underbrace{\left(\frac{\|p_\sigma\|}{E\|\mathcal{N}(0,I)\|} - 1\right)}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
    \text{unbiased about 0 under neutral selection}
    \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
}\bigg)

where

c_\sigma^{-1}\approx n/3 is the backward time horizon for the evolution path p_\sigma and larger than one,
\mu_w=\left(\sum_{i=1}^\mu w_i^2\right)^{-1} is the variance effective selection mass and 1 \le \mu_w \le \mu by definition of w_i,
C_k^{\;-1/2} = \sqrt{C_k}^{\;-1} = \sqrt{C_k^{\;-1}} is the unique symmetric square root of the inverse of C_k, and
d_\sigma is the damping parameter usually close to one. For d_\sigma=\infty or c_\sigma=0 the step-size remains unchanged.

The step-size \sigma_k is increased if and only if \|p_\sigma\| is larger than the expected value

\begin{align}E\|\mathcal{N}(0,I)\| &= \sqrt{2}\,\Gamma((n+1)/2)/\Gamma(n/2) 
  \\&\approx \sqrt{n}\,(1-1/(4\,n)+1/(21\,n^2)) \end{align}

and decreased if it is smaller. For this reason, the step-size update tends to make consecutive steps C_k^{-1}, in that after the adaptation has been successful \textstyle\left(\frac{m_{k+2}-m_{k+1}}{\sigma_{k+1}}\right)^T\! C_k^{-1} \frac{m_{k+1}-m_{k}}{\sigma_k} \approx 0.[1]

Finally, the covariance matrix is updated, where again the respective evolution path is updated first.


  p_c \gets \underbrace{(1-c_c)}_{\!\!\!\!\!\text{discount factor}\!\!\!\!\!}\, 
            p_c + 
     \underbrace{\mathbf{1}_{[0,\alpha\sqrt{n}]}(\|p_\sigma\|)}_{\text{indicator function}} 
     \overbrace{\sqrt{1 - (1-c_c)^2}}^{
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\text{complements for discounted variance}
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}
     \underbrace{\sqrt{\mu_w} 
     \, \frac{m_{k+1} - m_k}{\sigma_k}}_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
                      \text{distributed as}\; \mathcal{N}(0,C_k)\;\text{under neutral selection}
                      \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}

  C_{k+1} = \underbrace{(1 - c_1 - c_\mu + c_s)}_{\!\!\!\!\!\text{discount factor}\!\!\!\!\!}
               \, C_k + c_1 \underbrace{p_c p_c^T}_{
   \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
   \text{rank one matrix}
   \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!} 
         + \,c_\mu \underbrace{\sum_{i=1}^\mu w_i \frac{x_{i:\lambda} - m_k}{\sigma_k} 
             \left( \frac{x_{i:\lambda} - m_k}{\sigma_k} \right)^T}_{
                     \text{rank} \;\min(\mu,n)\; \text{matrix}}

where T denotes the transpose and

c_c^{-1}\approx n/4 is the backward time horizon for the evolution path p_c and larger than one,
\alpha\approx 1.5 and the indicator function \mathbf{1}_{[0,\alpha\sqrt{n}]}(\|p_\sigma\|) evaluates to one iff \|p_\sigma\|\in[0,\alpha\sqrt{n}] or, in other words, \|p_\sigma\|\le\alpha\sqrt{n}, which is usually the case,
c_s = (1 - \mathbf{1}_{[0,\alpha\sqrt{n}]}(\|p_\sigma\|)^2) \,c_1 c_c (2-c_c) makes partly up for the small variance loss in case the indicator is zero,
c_1 \approx 2 / n^2 is the learning rate for the rank-one update of the covariance matrix and
c_\mu \approx \mu_w / n^2 is the learning rate for the rank-\mu update of the covariance matrix and must not exceed 1 - c_1.

The covariance matrix update tends to increase the likelihood for p_c and for (x_{i:\lambda} - m_k)/\sigma_k to be sampled from \mathcal{N}(0,C_{k+1}). This completes the iteration step.

The number of candidate samples per iteration, \lambda, is not determined a priori and can vary in a wide range. Smaller values, for example \lambda=10, lead to more local search behavior. Larger values, for example \lambda=10n with default value \mu_w \approx \lambda/4, render the search more global. Sometimes the algorithm is repeatedly restarted with increasing \lambda by a factor of two for each restart.[2] Besides of setting \lambda (or possibly \mu instead, if for example \lambda is predetermined by the number of available processors), the above introduced parameters are not specific to the given objective function and therefore not meant to be modified by the user.

Example code in MATLAB/Octave

function xmin=purecmaes   % (mu/mu_w, lambda)-CMA-ES

  % --------------------  Initialization --------------------------------  
  % User defined input parameters (need to be edited)
  strfitnessfct = 'frosenbrock';  % name of objective/fitness function
  N = 20;               % number of objective variables/problem dimension
  xmean = rand(N,1);    % objective variables initial point
  sigma = 0.3;          % coordinate wise standard deviation (step size)
  stopfitness = 1e-10;  % stop if fitness < stopfitness (minimization)
  stopeval = 1e3*N^2;   % stop after stopeval number of function evaluations
  
  % Strategy parameter setting: Selection  
  lambda = 4+floor(3*log(N));  % population size, offspring number
  mu = lambda/2;               % number of parents/points for recombination
  weights = log(mu+1/2)-log(1:mu)'; % muXone array for weighted recombination
  mu = floor(mu);        
  weights = weights/sum(weights);     % normalize recombination weights array
  mueff=sum(weights)^2/sum(weights.^2); % variance-effectiveness of sum w_i x_i

  % Strategy parameter setting: Adaptation
  cc = (4+mueff/N) / (N+4 + 2*mueff/N);  % time constant for cumulation for C
  cs = (mueff+2) / (N+mueff+5);  % t-const for cumulation for sigma control
  c1 = 2 / ((N+1.3)^2+mueff);    % learning rate for rank-one update of C
  cmu = min(1-c1, 2 * (mueff-2+1/mueff) / ((N+2)^2+mueff));  % and for rank-mu update
  damps = 1 + 2*max(0, sqrt((mueff-1)/(N+1))-1) + cs; % damping for sigma 
                                                      % usually close to 1
  % Initialize dynamic (internal) strategy parameters and constants
  pc = zeros(N,1); ps = zeros(N,1);   % evolution paths for C and sigma
  B = eye(N,N);                       % B defines the coordinate system
  D = ones(N,1);                      % diagonal D defines the scaling
  C = B * diag(D.^2) * B';            % covariance matrix C
  invsqrtC = B * diag(D.^-1) * B';    % C^-1/2 
  eigeneval = 0;                      % track update of B and D
  chiN=N^0.5*(1-1/(4*N)+1/(21*N^2));  % expectation of 
                                      %   ||N(0,I)|| == norm(randn(N,1))
  
  % -------------------- Generation Loop --------------------------------
  counteval = 0;  % the next 40 lines contain the 20 lines of interesting code 
  while counteval < stopeval
    
      % Generate and evaluate lambda offspring
      for k=1:lambda,
          arx(:,k) = xmean + sigma * B * (D .* randn(N,1)); % m + sig * Normal(0,C) 
          arfitness(k) = feval(strfitnessfct, arx(:,k)); % objective function call
          counteval = counteval+1;
      end
    
      % Sort by fitness and compute weighted mean into xmean
      [arfitness, arindex] = sort(arfitness); % minimization
      xold = xmean;
      xmean = arx(:,arindex(1:mu))*weights;   % recombination, new mean value
    
      % Cumulation: Update evolution paths
      ps = (1-cs)*ps ... 
            + sqrt(cs*(2-cs)*mueff) * invsqrtC * (xmean-xold) / sigma; 
      hsig = norm(ps)/sqrt(1-(1-cs)^(2*counteval/lambda))/chiN < 1.4 + 2/(N+1);
      pc = (1-cc)*pc ...
            + hsig * sqrt(cc*(2-cc)*mueff) * (xmean-xold) / sigma;

      % Adapt covariance matrix C
      artmp = (1/sigma) * (arx(:,arindex(1:mu))-repmat(xold,1,mu));
      C = (1-c1-cmu) * C ...                  % regard old matrix  
           + c1 * (pc*pc' ...                 % plus rank one update
                   + (1-hsig) * cc*(2-cc) * C) ... % minor correction if hsig==0
           + cmu * artmp * diag(weights) * artmp'; % plus rank mu update

      % Adapt step size sigma
      sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1)); 
    
      % Decomposition of C into B*diag(D.^2)*B' (diagonalization)
      if counteval - eigeneval > lambda/(c1+cmu)/N/10  % to achieve O(N^2)
          eigeneval = counteval;
          C = triu(C) + triu(C,1)'; % enforce symmetry
          [B,D] = eig(C);           % eigen decomposition, B==normalized eigenvectors
          D = sqrt(diag(D));        % D is a vector of standard deviations now
          invsqrtC = B * diag(D.^-1) * B';
      end
    
      % Break, if fitness is good enough or condition exceeds 1e14, better termination methods are advisable 
      if arfitness(1) <= stopfitness || max(D) > 1e7 * min(D)
          break;
      end

  end % while, end generation loop

  xmin = arx(:, arindex(1)); % Return best point of last iteration.
                             % Notice that xmean is expected to be even
                             % better.
  
% ---------------------------------------------------------------  
function f=frosenbrock(x)
    if size(x,1) < 2 error('dimension must be greater one'); end
    f = 100*sum((x(1:end-1).^2 - x(2:end)).^2) + sum((x(1:end-1)-1).^2);

Theoretical Foundations

Given the distribution parametersmean, variances and covariancesthe normal probability distribution for sampling new candidate solutions is the maximum entropy probability distribution over \mathbb{R}^n, that is, the sample distribution with the minimal amount of prior information built into the distribution. More considerations on the update equations of CMA-ES are made in the following.

Variable Metric

The CMA-ES implements a stochastic variable-metric method. In the very particular case of a convex-quadratic objective function

 f(x) = {\textstyle\frac{1}{2}}(x-x^*)^T H (x-x^*)

the covariance matrix C_k adapts to the inverse of the Hessian matrix H, up to a scalar factor and small random fluctuations. More general, also on the function g \circ f, where  g is strictly increasing and therefore order preserving and f is convex-quadratic, the covariance matrix C_k adapts to H^{-1}, up to a scalar factor and small random fluctuations.

Maximum-Likelihood Updates

The update equations for mean and covariance matrix maximize a likelihood while resembling an expectation-maximization algorithm. The update of the mean vector m maximizes a log-likelihood, such that

 m_{k+1} = \arg\max_{m} \sum_{i=1}^\mu w_i \log p_\mathcal{N}(x_{i:\lambda} | m)

where

 \log p_\mathcal{N}(x) = 
   - \frac{1}{2} \log\det(2\pi C) - \frac{1}{2} (x-m)^T C^{-1} (x-m)

denotes the log-likelihood of x from a multivariate normal distribution with mean m and any positive definite covariance matrix C. To see that m_{k+1} is independent of C remark first that this is the case for any diagonal matrix C, because the coordinate-wise maximizer is independent of a scaling factor. Then, rotation of the data points or choosing C non-diagonal are equivalent.

The rank-\mu update of the covariance matrix, that is, the right most summand in the update equation of C_k, maximizes a log-likelihood in that

 \sum_{i=1}^\mu w_i \frac{x_{i:\lambda} - m_k}{\sigma_k} 
             \left( \frac{x_{i:\lambda} - m_k}{\sigma_k} \right)^T 
     = \arg\max_{C} \sum_{i=1}^\mu w_i \log p_\mathcal{N}\left(\left.\frac{x_{i:\lambda} - m_k}{\sigma_k} \right| C\right)

for \mu\ge n (otherwise C is singular, but substantially the same result holds for \mu < n). Here,  p_\mathcal{N}(x | C) denotes the likelihood of x from a multivariate normal distribution with zero mean and covariance matrix C. Therefore, for c_1=0 and c_\mu=1, C_{k+1} is the above maximum-likelihood estimator. See estimation of covariance matrices for details on the derivation.

Natural Gradient Descent in the Space of Sample Distributions

Akimoto et al.[3] and Glasmachers et al.[4] discovered independently that the update of the distribution parameters resembles the descend in direction of a sampled natural gradient of the expected objective function value E f (x) (to be minimized), where the expectation is taken under the sample distribution. With the parameter setting of c_\sigma=0 and c_1=0, i.e. without step-size control and rank-one update, CMA-ES can thus be viewed as an instantiation of Natural Evolution Strategies (NES).[3][4] The natural gradient is independent of the parameterization of the distribution. Taken with respect to the parameters θ of the sample distribution p, the gradient of E f (x) can be expressed as

 \begin{align} 
  {\nabla}_{\!\theta} E(f(x) | \theta) 
     &= \nabla_{\!\theta} \int_{\mathbb R^n}f(x) p(x) \mathrm{d}x
  \\ &= \int_{\mathbb R^n}f(x) \nabla_{\!\theta} p(x) \mathrm{d}x
  \\ &= \int_{\mathbb R^n}f(x) p(x) \nabla_{\!\theta} \ln p(x) \mathrm{d}x
  \\ &= E(f(x) \nabla_{\!\theta} \ln p(x|\theta))
\end{align}

where p(x)=p(x|\theta) depends on the parameter vector \theta, the so-called score function, \nabla_{\!\theta} \ln p(x|\theta) = \frac{\nabla_{\!\theta} p(x)}{p(x)} , indicates the relative sensitivity of p w.r.t. θ, and the expectation is taken with respect to the distribution p. The natural gradient of E f (x), complying with the Fisher information metric (an informational distance measure between probability distributions and the curvature of the relative entropy), now reads

 \begin{align} 
  \tilde{\nabla} E(f(x) | \theta) 
  &= F^{-1}_\theta \nabla_{\!\theta} E(f(x) | \theta) 
\end{align}

where the Fisher information matrix  F_{\theta} is the expectation of the Hessian of -lnp and renders the expression independent of the chosen parameterization. Combining the previous equalities we get

 \begin{align} 
  \tilde{\nabla} E(f(x) | \theta) 
  &= F^{-1}_\theta E(f(x) \nabla_{\!\theta} \ln p(x|\theta))
  \\ &= E(f(x) F^{-1}_\theta \nabla_{\!\theta} \ln p(x|\theta))
\end{align}

A Monte Carlo approximation of the latter expectation takes the average over λ samples from p

 \tilde{\nabla} \widehat{E}_\theta(f) := -\sum_{i=1}^\lambda \overbrace{w_i}^{\!\!\!\!\text{preference weight}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!} \underbrace{F^{-1}_\theta \nabla_{\!\theta} \ln p(x_{i:\lambda}|\theta)}_{\!\!\!\!\!\text{candidate direction from }x_{i:\lambda}\!\!\!\!\!}
\quad\mathrm{with~}w_i = -f(x_{i:\lambda})/\lambda

where the notation i:\lambda from above is used and therefore w_i are monotonously decreasing in i.

Ollivier et al.[5] finally found a rigorous formulation for the more robust weights, w_i, as they are defined in the CMA-ES (weights are zero for i > μ), formulated as consistent estimator for the CDF of f(X), X\sim p(.|\theta) at the point f(x_{i:\lambda}), composed with a fixed monotonous decreased transformation w, i.e.,

w_i = w\left(\frac{\mathsf{rank}(f(x_{i:\lambda})) - 1/2}{\lambda}\right)

Let

\theta = [m_k^T \mathrm{vec}(C_k)^T \sigma_k]^T \in \mathbb{R}^{n+n^2+1}

such that  p(.|\theta) is the density of the multivariate normal distribution \mathcal N(m_k,\sigma_k^2 C_k). Then, we have an explicit expression for the inverse of the Fisher information matrix where \sigma_k is fixed

F^{-1}_{\theta | \sigma_k} = \left[\begin{array}{cc}\sigma_k^2 C_k&0\\ 0&2 C_k\otimes C_k\end{array}\right]

and for

\ln p(x|\theta) = \ln p(x|m_k,\sigma_k^2 C_k) = -\frac{1}{2}(x-m_k)^T \sigma_k^{-2} C_k^{-1} (x-m_k) 
         \,-\, \frac{1}{2}\ln\det(2\pi\sigma_k^2 C_k)

and, after some calculations, the updates in the CMA-ES turn out as[3]

 \begin{align} 
   m_{k+1} 
     &= m_k - \underbrace{[\tilde{\nabla} \widehat{E}_\theta(f)]_{1,\dots, n}}_{
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
     \text{natural gradient for mean}
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
     } 
   \\
     &= m_k + \sum_{i=1}^\lambda w_i (x_{i:\lambda} - m_k) 
   \end{align}

and

 \begin{align} 
   C_{k+1} 
     &= C_k + c_1(p_c p_c^T - C_k)  
         - c_\mu\,\mathrm{mat}(\overbrace{[\tilde{\nabla} \widehat{E}_\theta(f)]_{n+1,\dots,n+n^2}}^{
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
        \text{natural gradient for covariance matrix} 
     \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
     })\\
     &= C_k + c_1(p_c p_c^T - C_k) 
        + c_\mu \sum_{i=1}^\lambda w_i \left(\frac{x_{i:\lambda} - m_k}{\sigma_k} \left(\frac{x_{i:\lambda} - m_k}{\sigma_k}\right)^T - C_k\right) 
   \end{align}

where mat forms the proper matrix from the respective natural gradient sub-vector. That means, setting c_1=c_\sigma=0, the CMA-ES updates descend in direction of the approximation  \tilde{\nabla} \widehat{E}_\theta(f) of the natural gradient while using different step-sizes (learning rates) for the orthogonal parameters m and C respectively.

Stationarity or Unbiasedness

It is comparatively easy to see that the update equations of CMA-ES satisfy some stationarity conditions, in that they are essentially unbiased. Under neutral selection, where x_{i:\lambda} \sim \mathcal N(m_k,\sigma_k^2 C_k), we find that

 E(m_{k+1}\,|\, m_k) = m_k

and under some mild additional assumptions on the initial conditions

 E(\log \sigma_{k+1} \,|\, \sigma_k) = \log \sigma_k

and with an additional minor correction in the covariance matrix update for the case where the indicator function evaluates to zero, we find

 E(C_{k+1} \,|\, C_k) = C_k

Invariance

Invariance properties imply uniform performance on a class of objective functions. They have been argued to be an advantage, because they allow to generalize and predict the behavior of the algorithm and therefore strengthen the meaning of empirical results obtained on single functions. The following invariance properties have been established for CMA-ES.

Any serious parameter optimization method should be translation invariant, but most methods do not exhibit all the above described invariance properties. A prominent example with the same invariance properties is the Nelder–Mead method, where the initial simplex must be chosen respectively.

Convergence

Conceptual considerations like the scale-invariance property of the algorithm, the analysis of simpler evolution strategies, and overwhelming empirical evidence suggest that the algorithm converges on a large class of functions fast to the global optimum, denoted as x^*. On some functions, convergence occurs independently of the initial conditions with probability one. On some functions the probability is smaller than one and typically depends on the initial m_0 and \sigma_0. Empirically, the fastest possible convergence rate in k for rank-based direct search methods can often be observed (depending on the context denoted as linear or log-linear or exponential convergence). Informally, we can write

\|m_k - x^*\| \;\approx\; \|m_0 - x^*\| \times e^{-ck}

for some c>0, and more rigorously

\frac{1}{k}\sum_{i=1}^k\log\frac{\|m_i - x^*\|}{\|m_{i-1} - x^*\|} 
   \;=\; \frac{1}{k}\log\frac{\|m_k - x^*\|}{\|m_{0} - x^*\|}
    \;\to\; -c < 0 \quad\text{for}\; k\to\infty\;,

or similarly,

E\log\frac{\|m_k - x^*\|}{\|m_{k-1} - x^*\|}
    \;\to\; -c < 0 \quad\text{for}\; k\to\infty\;.

This means that on average the distance to the optimum decreases in each iteration by a "constant" factor, namely by \exp(-c). The convergence rate c is roughly 0.1\lambda/n, given \lambda is not much larger than the dimension n. Even with optimal \sigma and C, the convergence rate c cannot largely exceed 0.25\lambda/n, given the above recombination weights w_i are all non-negative. The actual linear dependencies in \lambda and n are remarkable and they are in both cases the best one can hope for in this kind of algorithm. Yet, a rigorous proof of convergence is missing.

Interpretation as Coordinate System Transformation

Using a non-identity covariance matrix for the multivariate normal distribution in evolution strategies is equivalent to a coordinate system transformation of the solution vectors,[6] mainly because the sampling equation


\begin{align}
  x_i &\sim\ m_k + \sigma_k\times\mathcal{N}(0,C_k)
    \\
  &\sim\ m_k + \sigma_k \times C_k^{1/2}\mathcal{N}(0,I) 
\end{align}

can be equivalently expressed in an "encoded space" as


  \underbrace{C_k^{-1/2}x_i}_{\text{represented in the encode space}
    \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!} 
     \sim\ \underbrace{C_k^{-1/2} m_k} {} + \sigma_k \times\mathcal{N}(0,I)

The covariance matrix defines a bijective transformation (encoding) for all solution vectors into a space, where the sampling takes place with identity covariance matrix. Because the update equations in the CMA-ES are invariant under linear coordinate system transformations, the CMA-ES can be re-written as an adaptive encoding procedure applied to a simple evolution strategy with identity covariance matrix.[6] This adaptive encoding procedure is not confined to algorithms that sample from a multivariate normal distribution (like evolution strategies), but can in principle be applied to any iterative search method.

Performance in Practice

In contrast to most other evolutionary algorithms, the CMA-ES is, from the users perspective, quasi parameter-free. The user has to choose an initial solution point, m_0\in\mathbb{R}^n, and the initial step-size, \sigma_0>0. Optionally, the number of candidate samples λ (population size) can be modified by the user in order to change the characteristic search behavior (see above) and termination conditions can or should be adjusted to the problem at hand.

The CMA-ES has been empirically successful in hundreds of applications and is considered to be useful in particular on non-convex, non-separable, ill-conditioned, multi-modal or noisy objective functions. The search space dimension ranges typically between two and a few hundred. Assuming a black-box optimization scenario, where gradients are not available (or not useful) and function evaluations are the only considered cost of search, the CMA-ES method is likely to be outperformed by other methods in the following conditions:

On separable functions, the performance disadvantage is likely to be most significant in that CMA-ES might not be able to find at all comparable solutions. On the other hand, on non-separable functions that are ill-conditioned or rugged or can only be solved with more than 100 n function evaluations, the CMA-ES shows most often superior performance.

Variations and Extensions

The (1+1)-CMA-ES[7] generates only one candidate solution per iteration step which becomes the new distribution mean if it is better than the current mean. For c_c=1 the (1+1)-CMA-ES is a close variant of Gaussian adaptation. Some Natural Evolution Strategies are close variants of the CMA-ES with specific parameter settings. Natural Evolution Strategies do not utilize evolution paths (that means in CMA-ES setting c_c=c_\sigma=1) and they formalize the update of variances and covariances on a Cholesky factor instead of a covariance matrix. The CMA-ES has also been extended to multiobjective optimization as MO-CMA-ES.[8] Another remarkable extension has been the addition of a negative update of the covariance matrix with the so-called active CMA.[9]

With the advent of niching methods in evolutionary strategies, the question of an optimal niche radius arises. An "adaptive individual niche radius" is introduced in [10]

See also

References

  1. Hansen, N. (2006), "The CMA evolution strategy: a comparing review", Towards a new evolutionary computation. Advances on estimation of distribution algorithms, Springer, pp. 1769–1776
  2. Auger, A.; N. Hansen (2005). "A Restart CMA Evolution Strategy With Increasing Population Size" (PDF). 2005 IEEE Congress on Evolutionary Computation, Proceedings. IEEE. pp. 1769–1776.
  3. 1 2 3 Akimoto, Y.; Y. Nagata; I. Ono; S. Kobayashi (2010). "Bidirectional Relation between CMA Evolution Strategies and Natural Evolution Strategies". Parallel Problem Solving from Nature, PPSN XI. Springer. pp. 154–163.
  4. 1 2 Glasmachers, T.; T. Schaul; Y. Sun; D. Wierstra; J. Schmidhuber (2010). "Exponential Natural Evolution Strategies" (PDF). Genetic and Evolutionary Computation Conference GECCO. Portland, OR.
  5. Ollivier, Y.; Arnold, L.; Auger, A.; Hansen, N. (2013). "Information-Geometric Optimization Algorithms: A Unifying Picture via Invariance Principles". arXiv:1106.3708v2.
  6. 1 2 Hansen, N. (2008). "Adpative Encoding: How to Render Search Coordinate System Invariant". Parallel Problem Solving from Nature, PPSN X. Springer. pp. 205–214.
  7. Igel, C.; T. Suttorp; N. Hansen (2006). "A Computational Efficient Covariance Matrix Update and a (1+1)-CMA for Evolution Strategies" (PDF). Proceedings of the Genetic and Evolutionary Computation Conference (GECCO). ACM Press. pp. 453–460.
  8. Igel, C.; N. Hansen; S. Roth (2007). "Covariance Matrix Adaptation for Multi-objective Optimization". Evolutionary Computation (MIT press) 15 (1): 1–28. doi:10.1162/evco.2007.15.1.1. PMID 17388777.
  9. Jastrebski, G.A.; D.V. Arnold (2006). "Improving Evolution Strategies through Active Covariance Matrix Adaptation". 2006 IEEE World Congress on Computational Intelligence, Proceedings. IEEE. pp. 9719–9726. doi:10.1109/CEC.2006.1688662.
  10. Shir, Ofer M.; Bäck, Thomas (2006). Parallel Problem Solving from Nature-PPSN IX. Springer. pp. 142–151.

Bibliography

External links

This article is issued from Wikipedia - version of the Wednesday, March 23, 2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.