CMA-ES

CMA-ES означає Коваріаційна матриця стратегії еволюції адаптації.

Еволюційна стратегія (ЕС) є стохастичною, похідною від методів чисельної оптимізації нелінійних або неопуклих проблем неперервної оптимізації. Вони належать до класу еволюційних алгоритмів і еволюційних обчислень. Еволюційний алгоритм в цілому засновано на принципі біологічної еволюції, а саме повторній взаємодії варіації (через мутації і рекомбінації) і відбору: в кожному поколінні (ітерації) нові особи (розв'язки) породжують зміни, а потім деякі генерації вибирають для наступного покоління на основі їх придатності або на основі цільової функції значення. Подібно до цього, в послідовності поколінь створюються все кращі генерації.

В еволюційній стратегії, нові рішення відбирають відповідно до багатовимірного нормального розподілу. Парні залежності між змінними в цьому розподілі представлені у коваріаційній матриці. Адаптація коваріаційної матриці (CMA) є методом для поновлення коваріаційної матриці цього розподілу. Це особливо корисно, якщо функція є погано обумовленою.

Адаптація коваріаційної матриці становить вивчення другого порядку моделі базової цільової функції, схоже на обчислення зворотної матриці Гессе в класичній оптимізації. На відміну від більшості класичних методів, виконується менше припущень про природу основної цільової функції. Тільки рейтинг між кандидатами на найкраще рішення використані для вивчення розподілу вибірки, і ні похідних, ні навіть самі значення функції не використовують.

Принципи

Концепція адаптації коваріаційної матриці. У поколіннях, розподіл форм може адаптуватися до еліпсоїдальної або форми хребта.

В CMA-ES алгоритмі використовують два основні принципи для адаптації параметрів розподілу пошуку.

По-перше, принцип максимальної правдоподібності, що заснована на ідеї, підвищення ймовірності успішного вирішення кандидатів і пошуку кроків. Середній розподіл оновлюється, та ймовірність вибрати минулі успішні рішення є максимальним. Коваріаційна матриця розподілу оновлюється (поступово), що збільшує ймовірність того що будуть повторені попередні успішні кроки пошуку. Обидва оновлення можна розглядати як природний градієнт спуску. Крім того, в результаті, CMA проводить повторний аналіз головних компонентів успішного кроку пошуку, зберігаючи при цьому всі головні осі. Оцінка розподілу алгоритмів і методу крос-ентропії ґрунтується на дуже схожих ідеях, але оцінка коваріаційної матриці на максимальній ймовірності успішного вирішення вказується замість успішного пошуку кроку.

По-друге, існують два шляхи еволюції розподілу — пошук і розвиток шляхів. Ці шляхи містять важливу інформацію про кореляцію між послідовними кроками. Зокрема, якщо послідовні кроки йдуть в тому ж напрямку, еволюція шляхів стають довгою. Еволюція шляху експлуатується у двох напрямках. Один шлях використовується для адаптації процедури коваріаційної матриці замість одного успішного кроку пошуку і полегшує набагато швидше, дисперсією збільшення сприятливих напрямків. Інший шлях використовується для проведення додаткових змін розміру кроку управління. Розміру кроку управління прагне зробити послідовні рухи розподілу. Розмір кроку контролю ефективно запобігає передчасному зближенню та призводить до оптимального розв'язку.

Алгоритм

Відповідно до найбільш часто використовуваних (μ / μ W, λ)-CMA-ES створено алгоритм, де в кожній ітерації зважене поєднання μ з найкращим з рішень λ нового кандидата використовується для оновлення параметрів розподілу.

Основний цикл складається з трьох основних частин:

1. відбір нових рішень,
2. зміна порядку з обраних рішень, заснованих на їх придатність,
3. оновлення внутрішніх змінних стану на основі повторного замовлення зразків.

Псевдокод алгоритму виглядає наступним чином.

 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}$  

Приклад коду в 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.5;          % 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 = 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);

Практичне використання

На відміну від більшості інших еволюційних алгоритмів, CMA-ES майже не містить параметрів. Однак розмір популяції λ може бути змінений користувачем, щоб змінити поведінку характеру пошуку. CMA-ES була емпірично успішною в сотнях додатків і вважається корисною, зокрема, у випадку неопуклої, не сепарабельної, погано обумовленою, мультимодальної ​​цільової функції. Розмірність простору пошуку коливається зазвичай від двох до кількох сотень. Якщо припустити, що оптимізація сценаріїв — це чорний ящик, де градієнтів немає і функція оцінки виражає вартість пошуку,тоді CMA-ES, ймовірно, поступається іншим методам в наступних умовах:

• у випадку маломірних функцій, наприклад, симплекс-метод або сурогатні методи, засновані на сепарабельних функцій без залежностей між змінними проєктування, зокрема, у випадку мультимодальності;
• у випадку опукло-квадратичної функції з низьким або середнім числом обумовленості в матриці Гессе. (Тут BFGS або NEWUOA, як правило, в десять разів швидше);
• у випадку функції, яка може бути вирішена за допомогою порівняно невеликого числа оцінок функції, наприклад, NEWUOA або багаторівневий пошук координат (MCS).

Посилання

• Короткий вступ в CMA-ES Н. Хансен
• Стратегія CMA еволюції: Навчальний посібник
• CMA-ES вихідний код сторінки