СМА-ES

Стратегия эволюции адаптации ковариационной матрицы (CMA-ES) — это особый вид стратегии численной оптимизации . Стратегии эволюции (ES) — это методы для без производных численной оптимизации нелинейных стохастические или невыпуклых задач непрерывной оптимизации . Они относятся к классу эволюционных алгоритмов и эволюционных вычислений . Эволюционный алгоритм в целом основан на принципе биологической эволюции , а именно на повторяющемся взаимодействии вариаций (посредством рекомбинации и мутации) и отбора: в каждом поколении (итерации) появляются новые особи (кандидаты на решение, обозначаемые как $x$) генерируются вариациями текущих родительских особей, обычно стохастическим образом. Затем некоторые люди отбираются, чтобы стать родителями в следующем поколении, на основе их приспособленности или целевой функции. значения $f(x)$. Таким образом, люди со все лучшим и лучшим $f$-значения генерируются в ходе последовательности генерации.

В стратегии развития новые возможные решения обычно отбираются в соответствии с многомерным нормальным распределением . $\mathbb {R} ^{n}$. Рекомбинация сводится к выбору нового среднего значения для распределения. Мутация представляет собой добавление случайного вектора, возмущения с нулевым средним значением. Парные зависимости между переменными в распределении представлены ковариационной матрицей . Адаптация ковариационной матрицы (CMA) — это метод обновления ковариационной матрицы этого распределения. Это особенно полезно, если функция $f$ является плохо кондиционированным .

Адаптация ковариационной матрицы сводится к изучению модели второго порядка базовой целевой функции , аналогичной аппроксимации обратной матрицы Гессе в квазиньютоновском методе классической оптимизации . В отличие от большинства классических методов, делается меньше предположений относительно базовой целевой функции. Поскольку используется только ранжирование (или, что то же самое, сортировка) возможных решений, метод не требует ни производных, ни даже (явной) целевой функции. Например, рейтинг может быть получен в результате парных соревнований между кандидатами решений в турнире по швейцарской системе .

Принципы

Иллюстрация фактического запуска оптимизации с адаптацией ковариационной матрицы для простой двумерной задачи. Сферический ландшафт оптимизации изображен сплошными линиями равного размера. $f$-ценности. Популяция (точки) намного больше, чем необходимо, но ясно показывает, как меняется распределение совокупности (пунктирная линия) в ходе оптимизации. На этой простой задаче население концентрируется на достижении глобального оптимума в течение нескольких поколений.

В алгоритме 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 do  // iterate
    for $i$ in $\{1\ldots \lambda \}$ do  // sample $\lambda$ new solutions and evaluate them
        $x_{i}={}$sample_multivariate_normal(mean${}=m$, covariance_matrix${}=\sigma ^{2}C$)
        $f_{i}=\operatorname {fitness} (x_{i})$
    $x_{1\ldots \lambda }$ ← $x_{s(1)\ldots s(\lambda )}$ with $s(i)=\operatorname {argsort} (f_{1\ldots \lambda },i)$ // sort solutions
    $m'=m$  // we need later $m-m'$ and $x_{i}-m'$       
    $m$ ← update_m$(x_{1},\ldots ,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 ,\ldots ,(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}$

Порядок пяти назначений обновлений имеет значение: $m$ сначала необходимо обновить, $p_{\sigma }$ и $p_{c}$ необходимо обновить перед $C$, и $\sigma$ должен обновляться в последнюю очередь. Уравнения обновления для пяти переменных состояния указаны ниже.

Даны измерения пространства поиска. $n$ и шаг итерации $k$. Пять переменных состояния:

$m_{k}\in \mathbb {R} ^{n}$, среднее распределение и текущее любимое решение задачи оптимизации,

$\sigma _{k}>0$, размер шага,

$C_{k}$, симметричный и положительно определенный $n\times n$ ковариационная матрица с $C_{0}=I$ и

$p_{\sigma }\in \mathbb {R} ^{n},p_{c}\in \mathbb {R} ^{n}$, два пути эволюции, первоначально заданные для нулевого вектора.

Итерация начинается с выборки $\lambda >1$ возможные решения $x_{i}\in \mathbb {R} ^{n}$ из многомерного нормального распределения $\textstyle {\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})$, то есть для $i=1,\ldots ,\lambda$

${\begin{aligned}x_{i}\ &\sim \ {\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})\\&\sim \ m_{k}+\sigma _{k}\times {\mathcal {N}}(0,C_{k})\end{aligned}}$

Вторая строка предполагает интерпретацию как несмещенное возмущение (мутацию) текущего вектора избранного решения. $m_{k}$ (вектор среднего распределения). Возможные решения $x_{i}$ оцениваются по целевой функции $f:\mathbb {R} ^{n}\to \mathbb {R}$ быть сведено к минимуму. Обозначая $f$-отсортированные возможные решения как

$\{x_{i:\lambda }\mid i=1\dots \lambda \}=\{x_{i}\mid i=1\dots \lambda \}{\text{ and }}f(x_{1:\lambda })\leq \dots \leq f(x_{\mu :\lambda })\leq f(x_{\mu +1:\lambda })\leq \cdots ,$

новое среднее значение вычисляется как

${\begin{aligned}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{aligned}}$

где положительные (рекомбинационные) веса $w_{1}\geq w_{2}\geq \dots \geq w_{\mu }>0$ сумма до одного. Обычно $\mu \leq \lambda /2$ а веса выбираются так, что $\textstyle \mu _{w}:=1/\sum _{i=1}^{\mu }w_{i}^{2}\approx \lambda /4$. Единственная обратная связь, используемая от целевой функции здесь и далее, — это упорядочивание выборочных решений-кандидатов по индексам $i:\lambda$.

Размер шага $\sigma _{k}$ обновляется с использованием совокупной адаптации размера шага (CSA), иногда также называемой контролем длины пути . Путь эволюции (или путь поиска) $p_{\sigma }$ обновляется первым.

$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 }\|}{\operatorname {E} \|{\mathcal {N}}(0,I)\|}}-1\right)} _{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\text{unbiased about 0 under neutral selection}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}{\bigg )}$

где

$c_{\sigma }^{-1}\approx n/3$ это обратный временной горизонт для пути эволюции $p_{\sigma }$ и больше единицы ( $c_{\sigma }\ll 1$ напоминает константу экспоненциального затухания как $(1-c_{\sigma })^{k}\approx \exp(-c_{\sigma }k)$ где $c_{\sigma }^{-1}$ это связанный с ним срок службы и $c_{\sigma }^{-1}\ln(2)\approx 0.7c_{\sigma }^{-1}$ период полураспада),

$\mu _{w}=\left(\sum _{i=1}^{\mu }w_{i}^{2}\right)^{-1}$ - эффективная масса выбора дисперсии и $1\leq \mu _{w}\leq \mu$ по определению $w_{i}$,

$C_{k}^{\;-1/2}={\sqrt {C_{k}}}^{\;-1}={\sqrt {C_{k}^{\;-1}}}$ уникальным симметричным квадратным корнем обратного числа является $C_{k}$, и

$d_{\sigma }$ – параметр демпфирования обычно близок к единице. Для $d_{\sigma }=\infty$ или $c_{\sigma }=0$ размер шага остается неизменным.

Размер шага $\sigma _{k}$ увеличивается тогда и только тогда, когда $\|p_{\sigma }\|$ больше ожидаемого значения

${\begin{aligned}\operatorname {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{aligned}}$

и уменьшается, если оно меньше. По этой причине обновление размера шага имеет тенденцию делать последовательные шаги. $C_{k}^{-1}$-сопряженное , то есть после успешной адаптации $\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$.

Наконец, ковариационная матрица обновляется , причем сначала снова обновляется соответствующий путь эволюции.

$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}} _{\operatorname {rank} \min(\mu ,n){\text{ matrix}}}$

где $T$ обозначает транспонирование и

$c_{c}^{-1}\approx n/4$ это обратный временной горизонт для пути эволюции $p_{c}$ и больше единицы,

$\alpha \approx 1.5$ и индикаторная функция $\mathbf {1} _{[0,\alpha {\sqrt {n}}]}(\|p_{\sigma }\|)$ оценивается как один iff $\|p_{\sigma }\|\in [0,\alpha {\sqrt {n}}]$ или, другими словами, $\|p_{\sigma }\|\leq \alpha {\sqrt {n}}$, что обычно и бывает,

$c_{s}=(1-\mathbf {1} _{[0,\alpha {\sqrt {n}}]}(\|p_{\sigma }\|)^{2})\,c_{1}c_{c}(2-c_{c})$ частично компенсирует небольшую потерю дисперсии в случае, если индикатор равен нулю,

$c_{1}\approx 2/n^{2}$ - скорость обучения для обновления ковариационной матрицы первого ранга и

$c_{\mu }\approx \mu _{w}/n^{2}$ это скорость обучения для ранга $\mu$ обновление ковариационной матрицы и не должно превышать $1-c_{1}$.

Обновление ковариационной матрицы тенденцию увеличивать вероятность имеет $p_{c}$ и для $(x_{i:\lambda }-m_{k})/\sigma _{k}$ быть отобранным из ${\mathcal {N}}(0,C_{k+1})$. На этом этап итерации завершен.

Количество образцов-кандидатов на итерацию, $\lambda$, не определяется априори и может изменяться в широких пределах. Меньшие значения, например $\lambda =10$, приводят к более локальному поиску. Большие значения, например $\lambda =10n$ со значением по умолчанию $\mu _{w}\approx \lambda /4$, сделайте поиск более глобальным. Иногда алгоритм неоднократно перезапускается с увеличением $\lambda$ в два раза за каждый перезапуск. Помимо настройки $\lambda$ (или возможно $\mu$ вместо этого, если, например, $\lambda$ предопределено количеством доступных процессоров), введенные выше параметры не являются специфичными для данной целевой функции и, следовательно, не предназначены для изменения пользователем.

Пример кода в 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.
end
% ---------------------------------------------------------------  
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);
end

Теоретические основы

Учитывая параметры распределения — среднее значение, дисперсии и ковариации — нормальное распределение вероятностей для выборки новых решений-кандидатов представляет собой распределение вероятностей максимальной энтропии по $\mathbb {R} ^{n}$, то есть выборочное распределение с минимальным количеством априорной информации, встроенной в распределение. Дополнительные соображения по уравнениям обновления CMA-ES приведены ниже.

Переменная метрика

CMA-ES реализует стохастический метод переменной метрики . В весьма частном случае выпукло-квадратичной целевой функции

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

ковариационная матрица $C_{k}$ адаптируется к обратной матрице Гессе $H$, с точностью до скалярного множителя и малых случайных флуктуаций. В более общем плане, также о функции $g\circ f$, где $g$ строго возрастает и, следовательно, сохраняет порядок, ковариационная матрица $C_{k}$ адаптируется к $H^{-1}$, с точностью до скалярного множителя и малых случайных флуктуаций. Для соотношения выбора $\lambda /\mu \to \infty$ (и, следовательно, численность населения $\lambda \to \infty$), $\mu$ выбранные решения дают эмпирическую ковариационную матрицу, отражающую обратный гессиан, даже в стратегиях эволюции без адаптации ковариационной матрицы. Этот результат был доказан для $\mu =1$ на статической модели, опираясь на квадратичное приближение.

Обновления с максимальной вероятностью

Уравнения обновления для среднего и ковариационной матрицы максимизируют правдоподобие , напоминая алгоритм максимизации ожидания . Обновление среднего вектора $m$ максимизирует логарифмическое правдоподобие, такое, что

$m_{k+1}=\arg \max _{m}\sum _{i=1}^{\mu }w_{i}\log p_{\mathcal {N}}(x_{i:\lambda }\mid m)$

где

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

обозначает логарифмическую вероятность $x$ из многомерного нормального распределения со средним $m$ и любая положительно определенная ковариационная матрица $C$. Чтобы увидеть это $m_{k+1}$ не зависит от $C$ сначала отметим, что это справедливо для любой диагональной матрицы $C$, поскольку покоординатный максимизатор не зависит от масштабного коэффициента. Затем вращение точек данных или выбор $C$ недиагональные эквивалентны.

Ранг- $\mu$ обновление ковариационной матрицы, то есть самого правого слагаемого в уравнении обновления $C_{k}$, максимизирует логарифмическую вероятность в этом

$\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)$

для $\mu \geq n$ (в противном случае $C$ сингулярно, но по существу тот же результат справедлив и для $\mu <n$). Здесь, $p_{\mathcal {N}}(x|C)$ обозначает вероятность $x$ из многомерного нормального распределения с нулевым средним и ковариационной матрицей $C$. Поэтому для $c_{1}=0$ и $c_{\mu }=1$, $C_{k+1}$ — приведенная выше оценка максимального правдоподобия . см . в разделе «Оценка ковариационных матриц» Подробную информацию о выводе .

градиентный спуск в пространстве

Акимото и др. и Глазмахерс и др. независимо обнаружило, что обновление параметров распределения напоминает спуск в направлении выборочного естественного градиента ожидаемого значения целевой функции. $Ef(x)$ (подлежит минимизации), где математическое ожидание берется под выборочное распределение. При настройке параметра $c_{\sigma }=0$ и $c_{1}=0$, то есть без контроля размера шага и обновления первого ранга, CMA-ES, таким образом, можно рассматривать как реализацию стратегий естественной эволюции (NES). Естественный не градиент зависит от параметризации распределения. Принимая во внимание параметры θ выборочного распределения p , градиент $Ef(x)$ может быть выражено как

${\begin{aligned}{\nabla }_{\!\theta }E(f(x)\mid \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\\&=\operatorname {E} (f(x)\nabla _{\!\theta }\ln p(x\mid \theta ))\end{aligned}}$

где $p(x)=p(x\mid \theta )$ зависит от вектора параметров $\theta$. Так называемая функция оценки , $\nabla _{\!\theta }\ln p(x\mid \theta )={\frac {\nabla _{\!\theta }p(x)}{p(x)}}$, указывает относительную чувствительность p относительно θ , а математическое ожидание берется относительно распределения p . градиент Естественный ​ ​ $Ef(x)$, соответствующий информационной метрике Фишера (мера информационного расстояния между распределениями вероятностей и кривизной относительной энтропии ), теперь читается как

${\begin{aligned}{\tilde {\nabla }}\operatorname {E} (f(x)\mid \theta )&=F_{\theta }^{-1}\nabla _{\!\theta }\operatorname {E} (f(x)\mid \theta )\end{aligned}}$

где Фишера информационная матрица $F_{\theta }$ является ожиданием гессиана от −ln p и делает выражение независимым от выбранной параметризации. Объединяя предыдущие равенства, получаем

${\begin{aligned}{\tilde {\nabla }}\operatorname {E} (f(x)\mid \theta )&=F_{\theta }^{-1}\operatorname {E} (f(x)\nabla _{\!\theta }\ln p(x\mid \theta ))\\&=\operatorname {E} (f(x)F_{\theta }^{-1}\nabla _{\!\theta }\ln p(x\mid \theta ))\end{aligned}}$

Приближение Монте-Карло последнего ожидания берет среднее значение по λ выборкам из p

${\tilde {\nabla }}{\widehat {E}}_{\theta }(f):=-\sum _{i=1}^{\lambda }\overbrace {w_{i}} ^{\!\!\!\!{\text{preference weight}}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!}\underbrace {F_{\theta }^{-1}\nabla _{\!\theta }\ln p(x_{i:\lambda }\mid \theta )} _{\!\!\!\!\!{\text{candidate direction from }}x_{i:\lambda }\!\!\!\!\!}\quad {\text{with }}w_{i}=-f(x_{i:\lambda })/\lambda$

где обозначение $i:\lambda$ сверху используется и поэтому $w_{i}$ монотонно уменьшаются в $i$.

Оливье и др. наконец нашел строгий вывод для весов, $w_{i}$, как они определены в CMA-ES. Веса представляют собой непротиворечивую CDF асимптотически оценку $f(X)$ в точках г. $i$ порядка статистика $f(x_{i:\lambda })$, как определено выше, где $X\sim p(.|\theta )$, составленный с фиксированным монотонно убывающим преобразованием $w$, то есть,

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

Эти веса делают алгоритм нечувствительным к конкретным $f$-ценности. Более кратко, используя CDF оценщик $f$ вместо $f$ пусть алгоритм зависит только от ранга $f$-ценностей, но не их основного распределения. Это делает алгоритм инвариантным к строго возрастанию $f$-преобразования. Теперь мы определяем

$\theta =[m_{k}^{T}\operatorname {vec} (C_{k})^{T}\sigma _{k}]^{T}\in \mathbb {R} ^{n+n^{2}+1}$

такой, что $p(\cdot \mid \theta )$ - плотность многомерного нормального распределения ${\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})$. Тогда у нас есть явное выражение для обратной информационной матрицы Фишера, где $\sigma _{k}$ фиксировано

$F_{\theta \mid \sigma _{k}}^{-1}=\left[{\begin{array}{cc}\sigma _{k}^{2}C_{k}&0\\0&2C_{k}\otimes C_{k}\end{array}}\right]$

и для

$\ln p(x\mid \theta )=\ln p(x\mid 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})$

и после некоторых вычислений обновления в CMA-ES оказываются такими:

${\begin{aligned}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{aligned}}$

и

${\begin{aligned}C_{k+1}&=C_{k}+c_{1}(p_{c}p_{c}^{T}-C_{k})-c_{\mu }\operatorname {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{aligned}}$

где mat образует правильную матрицу из соответствующего подвектора естественного градиента. Это означает, что установка $c_{1}=c_{\sigma }=0$, обновления CMA-ES спускаются в направлении аппроксимации ${\tilde {\nabla }}{\widehat {E}}_{\theta }(f)$ естественного градиента при использовании разных размеров шага (скорость обучения 1 и $c_{\mu }$) для ортогональных параметров $m$ и $C$ соответственно. Более поздние версии допускают другую скорость обучения для среднего значения. $m$ также. Самая последняя версия CMA-ES также использует другую функцию. $w$ для $m$ и $C$ с отрицательными значениями только для последнего (так называемая активная CMA).

Стационарность или несмещенность

Сравнительно легко увидеть, что уравнения обновления CMA-ES удовлетворяют некоторым условиям стационарности, поскольку они по существу несмещены. При нейтральном отборе, где $x_{i:\lambda }\sim {\mathcal {N}}(m_{k},\sigma _{k}^{2}C_{k})$, мы находим это

$\operatorname {E} (m_{k+1}\mid m_{k})=m_{k}$

и при некоторых мягких дополнительных предположениях о начальных условиях

$\operatorname {E} (\log \sigma _{k+1}\mid \sigma _{k})=\log \sigma _{k}$

и с дополнительной незначительной поправкой в ​​обновлении ковариационной матрицы для случая, когда индикаторная функция равна нулю, мы находим

$\operatorname {E} (C_{k+1}\mid C_{k})=C_{k}$

Инвариантность

Свойства инвариантности подразумевают единообразную производительность для класса целевых функций. Утверждалось, что они являются преимуществом, поскольку позволяют обобщать и предсказывать поведение алгоритма и, следовательно, усиливать значение эмпирических результатов, полученных для отдельных функций. Для CMA-ES установлены следующие свойства инвариантности.

• Инвариантность относительно сохраняющих порядок преобразований значения целевой функции $f$, в том, что для любого $h:\mathbb {R} ^{n}\to \mathbb {R}$ поведение идентично на $f:x\mapsto g(h(x))$ для всех строго возрастающих $g:\mathbb {R} \to \mathbb {R}$. Эту инвариантность легко проверить, поскольку только $f$В алгоритме используется ранжирование, инвариантное относительно выбора $g$.
• Масштабная инвариантность , то есть для любого $h:\mathbb {R} ^{n}\to \mathbb {R}$ поведение не зависит от $\alpha >0$ для целевой функции $f:x\mapsto h(\alpha x)$ данный $\sigma _{0}\propto 1/\alpha$ и $m_{0}\propto 1/\alpha$.
• Инвариантность относительно вращения пространства поиска в том, что для любого $h:\mathbb {R} ^{n}\to \mathbb {R}$ и любой $z\in \mathbb {R} ^{n}$ поведение на $f:x\mapsto h(Rx)$ не зависит от ортогональной матрицы $R$, данный $m_{0}=R^{-1}z$. В более общем плане алгоритм также инвариантен относительно общих линейных преобразований. $R$ когда дополнительно исходная ковариационная матрица выбирается как $R^{-1}{R^{-1}}^{T}$.

Любой серьезный метод оптимизации параметров должен быть трансляционно-инвариантным, но большинство методов не обладают всеми описанными выше свойствами инвариантности. Ярким примером с теми же свойствами инвариантности является метод Нелдера–Мида , где исходный симплекс должен быть выбран соответственно.

Конвергенция

Концептуальные соображения, такие как свойство алгоритма масштабной инвариантности, анализ более простых стратегий эволюции и подавляющее количество эмпирических данных, позволяют предположить, что алгоритм быстро сходится на большом классе функций к глобальному оптимуму, обозначаемому как $x^{*}$. На некоторых функциях сходимость происходит независимо от начальных условий с вероятностью единица. Для некоторых функций вероятность меньше единицы и обычно зависит от начального значения. $m_{0}$ и $\sigma _{0}$. Эмпирически самая быстрая возможная скорость сходимости в $k$ для методов прямого поиска на основе рангов часто можно наблюдать (в зависимости от контекста это называется линейной сходимостью , лог-линейной или экспоненциальной сходимостью). Неформально мы можем написать

$\|m_{k}-x^{*}\|\;\approx \;\|m_{0}-x^{*}\|\times e^{-ck}$

для некоторых $c>0$и более строго

${\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 \;,$

или аналогично,

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

Это означает, что в среднем расстояние до оптимума уменьшается на каждой итерации в «постоянный» раз, а именно на $\exp(-c)$. Скорость сходимости $c$ примерно $0.1\lambda /n$, данный $\lambda$ ненамного больше размера $n$. Даже при оптимальном $\sigma$ и $C$, скорость сходимости $c$ не может существенно превышать $0.25\lambda /n$, учитывая указанные выше веса рекомбинации $w_{i}$ все неотрицательны. Реальные линейные зависимости в $\lambda$ и $n$ замечательны и в обоих случаях являются лучшими, на которые можно надеяться в этом алгоритме. Однако строгое доказательство сходимости отсутствует.

Интерпретация как преобразование системы координат

Использование нетождественной ковариационной матрицы для многомерного нормального распределения в стратегиях эволюции эквивалентно преобразованию системы координат векторов решения: главным образом потому, что уравнение выборки

${\begin{aligned}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{aligned}}$

может быть эквивалентно выражено в «закодированном пространстве» как

$\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)$

Матрица ковариации определяет биективное преобразование (кодирование) для всех векторов решений в пространство, где выборка происходит с единичной матрицей ковариации. Поскольку уравнения обновления в CMA-ES инвариантны относительно преобразований линейной системы координат, CMA-ES можно переписать как процедуру адаптивного кодирования, применяемую к простой стратегии эволюции с единичной ковариационной матрицей. Эта процедура адаптивного кодирования не ограничивается алгоритмами выборки из многомерного нормального распределения (например, стратегии эволюции), но в принципе может применяться к любому итерационному методу поиска.

на

В отличие от большинства других эволюционных алгоритмов , CMA-ES, с точки зрения пользователя, не содержит квазипараметров. Пользователь должен выбрать начальную точку решения, $m_{0}\in \mathbb {R} ^{n}$и начальный размер шага, $\sigma _{0}>0$. При желании количество выборок-кандидатов λ (размер популяции) может быть изменено пользователем, чтобы изменить характерное поведение поиска (см. Выше), а условия завершения могут или должны быть скорректированы в соответствии с рассматриваемой проблемой.

CMA-ES оказался эмпирически успешным в сотнях приложений и считается полезным, в частности, для невыпуклых, неразделимых, плохо обусловленных, мультимодальных или зашумленных целевых функций. Одно исследование оптимизации «черного ящика» показало, что он превзошел 31 другой алгоритм оптимизации, особенно хорошо действуя на «сложных функциях» или в пространствах поиска большей размерности.

Размерность пространства поиска обычно колеблется от двух до нескольких сотен. Предполагая сценарий оптимизации «черного ящика», где градиенты недоступны (или бесполезны), а оценки функций являются единственной рассматриваемой стоимостью поиска, метод CMA-ES, вероятно, будет превзойден другими методами в следующих условиях:

• на маломерных функциях, скажем $n<5$, например, с помощью нисходящего симплекс-метода или суррогатных методов (например, кригинга с ожидаемым улучшением);
• на разделимых функциях без или с незначительными зависимостями между переменными расчета, в частности, в случае мультимодальности или большой размерности, например, путем дифференциальной эволюции ;
• на (почти) выпукло -квадратичных функциях с низким или умеренным числом обусловленности матрицы Гессе , где BFGS , NEWUOA или SLSQP обычно работают как минимум в десять раз быстрее;
• о функциях, которые уже могут быть решены сравнительно небольшим числом вычислений, скажем, не более $10n$, где CMA-ES часто работает медленнее, чем, например, NEWUOA или многоуровневый поиск координат (MCS).

В случае с разделяемыми функциями недостаток производительности, вероятно, будет наиболее существенным, поскольку CMA-ES вообще не сможет найти сопоставимые решения. С другой стороны, для неразделимых функций, которые плохо обусловлены или трудноразрешимы или могут быть решены только с помощью более чем $100n$ При оценке функций CMA-ES чаще всего показывает превосходную производительность.

Вариации и расширения

(1+1)-CMA-ES генерирует только одно решение-кандидат на шаг итерации, которое становится новым средним значением распределения, если оно лучше текущего среднего значения. Для $c_{c}=1$ (1+1)-CMA-ES является близким вариантом гауссовой адаптации . Некоторые стратегии естественной эволюции являются близкими вариантами CMA-ES с особыми настройками параметров. Стратегии естественной эволюции не используют пути эволюции (это означает, что в условиях CMA-ES $c_{c}=c_{\sigma }=1$) и они формализуют обновление дисперсий и ковариаций по фактору Холецкого вместо ковариационной матрицы. CMA-ES также был расширен до многокритериальной оптимизации как MO-CMA-ES. Еще одним замечательным расширением стало добавление отрицательного обновления ковариационной матрицы с помощью так называемого активного CMA. Использование дополнительного активного обновления CMA в настоящее время считается вариантом по умолчанию.

См. также

• Глобальная оптимизация – Отделение математики
• Стохастическая оптимизация — метод оптимизации
• Оптимизация без производных – Математическая дисциплина
• Алгоритм оценки распределения - Семейство методов стохастической оптимизации

Библиография

• Хансен Н., Остермайер А. (2001). Полностью дерандомизированная самоадаптация в стратегиях эволюции. Эволюционные вычисления , 9 (2), стр. 159–195. [1]
• Хансен Н., Мюллер С.Д., Кумутсакос П. (2003). Снижение временной сложности стратегии дерандомизированной эволюции с адаптацией ковариационной матрицы (CMA-ES). Эволюционные вычисления , 11 (1), стр. 1–18. [2]
• Хансен Н., Керн С. (2004). Оценка стратегии развития CMA на мультимодальных тестовых функциях. В книге Синь Яо и др., редакторы, «Параллельное решение проблем из природы» – PPSN VIII , стр. 282–291, Springer. [3]
• Игель С., Хансен Н., Рот С. (2007). Адаптация ковариационной матрицы для многокритериальной оптимизации. Эволюционные вычисления , 15 (1), стр. 1–28. [4]

Внешние ссылки

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