Тесты

Оптимизация по алгоритму фейерверков

Джеймс Маккафри

Исходный код можно скачать по ссылке

James McCaffreyВо многих программных системах с машинным обучением используют числовую оптимизацию. Например, в классификации по логистической регрессии (logistic regression classification) обучение классификатора — это процесс поиска набора значений для весов, связанных с входными переменными так, чтобы для набора обучающих данных вычисленные выходные значения близко соответствовали известным значениям выходных переменных. Иначе говоря, цель заключается в оптимизации (минимизации) ошибки между вычисленными и нужными выходными значениями.

Существует много различных алгоритмов численной оптимизации. Обратное распространение ошибок (back-propagation) основано на классических численных методах и часто применяется в нейронных сетях. Оптимизация роя частиц (particle swarm optimization) имитирует групповое поведение стаи птиц. Эволюционный алгоритм оптимизации (разновидность генетического алгоритма оптимизации) моделирует биологические процессы репликации хромосом.

В этой статье описывается сравнительно новый (впервые опубликованный в 2010 году) метод оптимизации, который называют оптимизацией по алгоритму фейерверков (fireworks algorithm optimization, FAO). Этот метод не моделирует и не имитирует явным образом поведение фейерверков, но, когда алгоритм визуализируется, получаемое изображение напоминает геометрию взрывающихся фейерверков.

Лучший способ понять, куда я клоню в этой статье, и получить представление о том, что такое FAO, — взглянуть на демонстрационную программу на рис. 1. Цель этой программы — использовать FAO для нахождения минимального значения функции Экли (Ackley’s function) с десятью переменными. Функция Экли является стандартным эталонным тестом для оценки эффективности алгоритмов численной оптимизации. В демонстрационной версии есть минимальное значение 0.0, находящееся в (0, 0, 0, 0, 0, 0, 0, 0, 0, 0). Функция Экли весьма изощренная, поскольку существует множество локальных минимальных решений, в которые могут попасть алгоритмы поиска до нахождения одного глобального минимума. График функции Экли двух переменных показан на рис. 2.

Оптимизация по алгоритму фейерверков в действии
Рис. 1. Оптимизация по алгоритму фейерверков в действии

Функция Экли двух переменных
Рис. 2. Функция Экли двух переменных

Демонстрационная программа задает количество фейерверков равным 5. Чем больше фейерверков, тем выше шансы на нахождение истинно оптимального решения, но за счет производительности. FAO обычно хорошо работает с количеством фейерверков от 3 до 10. FAO — процесс итеративный и требует задания максимального значения счетчика цикла. Переменную — счетчик цикла в оптимизации машинного обучения часто называют эпохой (epoch), и в демонстрации ей присваивается значение 1000 (т. е. максимум 1000 итераций). Это довольно мало и рассчитано лишь на демонстрацию; в реальных сценариях типичны значения от 10 000 до 100 000.

В демонстрационной программе, выполнение которой иллюстрирует рис. 1, наименьшая ошибка, связанная с лучшей позицией, найденной на данный момент, отображается через каждые 100 итераций (эпох). Заметьте, что FAO очень быстро начинает конвергенцию. После 1000 итераций лучшая из найденных позиций — (0.034 0.098 0.003 0.132 –0.054 0.181 –0.018 0.051 0.004 –0.023), а связанная с ней ошибка равна 0.41. Если максимальное количество эпох увеличить до 10 000, то FAO фактически найдет оптимальное решение. FAO, как и большинство алгоритмов численной оптимизации, не гарантирует нахождение оптимального решения в реальных сценариях, где таковое решение не известно.

В этой статье предполагается, что владеете навыками программирования хотя бы на среднем уровне, но ничего не знаете о численной оптимизации или о FAO. Демонстрационная программа написана на C#, но у вас не должно возникнуть особых трудностей в рефакторинге кода под другой язык, такой как JavaScript или Python.

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

Общая структура программы

Общая структура программы с небольшими правками для экономии места представлена на рис. 3. Чтобы создать демонстрационную программу, я запустил Visual Studio и создал новое консольное приложение на C# с названием FireworksAlgorithm. В этой программе нет значимых зависимостей от .NET, поэтому подойдет любая недавняя версия Visual Studio.

Рис. 3. Общая структура программы

using System;
using System.Collections.Generic;
namespace FireworksAlgorithm
{
  class FireworksProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin fireworks algorithm demo\n");
      Console.WriteLine("Goal is to solve Ackley's function");
      Console.WriteLine("Function has min value = 0.0 at (0, .. ))");
      int dim = 10;
      int n = 5;    // число фейерверков
      int maxEpochs = 10000;
      Console.WriteLine("\nSetting Ackley's dimension to " + dim);
      Console.WriteLine("Setting number fireworks to " + n);
      Console.WriteLine("Setting maxEpochs to " + maxEpochs);
      Console.WriteLine("\nBegin algorithm\n");
      double[] bestPosition = Solve(dim, n, maxEpochs);
      Console.WriteLine("\nAlgorithm complete");
      Console.WriteLine("\nBest solution found: ");
      for (int i = 0; i < dim; ++i)
        Console.Write(bestPosition[i].ToString("F3") + " ");
      Console.WriteLine();
      double error = Error(bestPosition);
      Console.WriteLine("\nError of best solution found = " +
        error.ToString("F5"));
      Console.WriteLine("\nEnd fireworks algorithm demo\n");
      Console.ReadLine();
    } // Main
    public static double Error(double[] position) { . . }
    public static double[] Solve(int dim, int n, int maxEpochs) { . . }
    private static int[] PickDimensions(int dim, int z, int seed) { . . }
    private static double YMax(Info[] fireworks) { . . }
    private static double YMin(Info[] fireworks) { . . }
    private static int[] NumberSparks(Info[] fireworks, int m,
      double a, double b) { . . }
    private static double[] Amplitudes(Info[] fireworks, int A,
      int epoch, int maxEpochs, double minX, double maxX) { . . }
    private static double MinAmplitude(int epoch, int maxEpochs,
      double minX, double maxX) { . . }
    private static void AddGaussianSparks(Info[] fireworks,
      List<Info>[] sparksList, int dim, int mHat, int epoch, double minX,
      double maxX, double[] bestPosition, ref double bestError, Random rnd)
  }
  public class Info
  {
    public double[] position;
    public double error;
  }
} // ns

После загрузки кода шаблона в редактор Visual Studio я переименовал в окне Solution Explorer файл Program.cs в более описательный FireworksProgram.cs, и Visual Studio автоматически переименовала класс Program за меня. В начале кода я удалил все выражения using, которые указывали на ненужные пространства имен, оставив только ссылки на System и Collections.Generic.

Я кодировал демонстрацию, используя подход со статическими методами, а не ООП. Вся управляющая логика находится в методе Main, который вызывает два открытых метода: Solve и Error. Важные выражения вызовов выглядят так:

int dim = 10;
int n = 5;
int maxEpochs = 1000;
double[] bestPosition = Solve(dim, n, maxEpochs)
double error = Error(bestPosition);

Переменная dim хранит количество измерений задачи. В случае машинного обучения она должна была бы содержать, как правило, количество определяемых весовых значений модели. Переменная n — это число фейерверков. Я использую по возможности предельно лаконичные имена переменных вроде n, которые применяются и в научных статьях по FAO, чтобы вам было легче читать эти статьи как дополнительный источник. Алгоритм FAO содержится в методе Solve. Метод Error принимает позицию (массив из десяти значений), вычисляет значение функции Экли для этих значений и возвращает среднее суммы квадратов разности между вычисленными выходными значениями и известным минимум 0.0.

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

Вникаем в алгоритм

Схема алгоритма фейерверков показана на рис. 4. Изображение представляет упрощенную (надуманную) задачу минимизации, а не функцию Экли из демонстрационной программы. Цель надуманной задачи — минимизировать некую произвольную функцию, имеющую две входные переменные, X и Y, где функция имеет минимум 0.0 при X = 0 и Y = 0.

Алгоритм фейерверков
Рис. 4. Алгоритм фейерверков

Fireworks algorithm optimization Оптимизация по алгоритму фейерверков
(worst spark) (худший огонек)
firework[0] firework[0]
firework[1] firework[1]
firework[2] firework[2]
min value at (0,0) мин. значение при (0,0)
(best spark) (лучший огонек)
(Gaussian spark) (гауссов огонек)

На рис. 4 используется три фейерверка. В каждом фейерверке есть так называемые огоньки (sparks). Существует два вида огоньков: обычные и гауссовы. На рис. 5 показан алгоритм фейерверков в виде очень высокоуровневого псевдокода.

Рис. 5. Алгоритм фейерверков в виде высокоуровневого псевдокода

Инициализируем три фейерверка в случайных позициях
Цикл until done
  for each фейерверк
    Вычисляем радиус разлета (amplitude) фейерверка
    Вычисляем число обычных огоньков
    Генерируем обычные огоньки
  end for
  Генерируем особые гауссовы огоньки
  Оцениваем каждый огонек
  Из списка огоньков выбираем три,
    которые будут позициями новых фейерверков
  Создаем три новых фейерверка
Конец цикла
Возвращаем позицию лучшего найденного огонька

В случае изображения на рис. 4 позиции трех фейерверков обозначаются точечными маркерами в кружках по координатам (–6, 2), (3, 4) и (3, –3). Заметьте, что первый фейерверк является худшим, так как он дальше всех от целевой позиции (0, 0), а третий фейерверк — лучшим, потому он ближе всего. Радиус разлета (amplitudes) обозначается пунктирными окружностями вокруг каждого фейерверка. Хорошие фейерверки имеют меньшие радиусы разлета, а плохие фейерверки — большие.

Каждый фейерверк будет генерировать разное количество обычных (не гауссовых) огоньков. Плохие фейерверки будут генерировать меньше огоньков, чем хорошие. На рис. 4 firework[0] генерирует три огонька, firework[1] — четыре, а firework[2] — пять. Каждый обычный огонек располагается внутри области разлета его родительского фейерверка. Поскольку у хороших фейерверков малый радиус разлета и относительно много огоньков, алгоритм при поиске фокусируется рядом с хорошими фейерверками. Пренебрегать плохими фейерверками не следует, так как они могли бы привести к хорошему решению, которое находится вне диапазона текущих фейерверков. Гауссовы огоньки генерируются так, чтобы они располагались между фейерверком и лучшей известной позицией любого огонька, встреченного во время поиска.

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

Процесс генерации фейерверков, а затем огоньков идет в цикле до тех пор, пока не будет встречено некое условие остановки. Такое условие — это просто максимальное значение счетчика цикла в случае демонстрационной программы. При использовании FAO для машинного обучения условием остановки может быть порог ошибки: когда происходит выход за пороговое значение, зависящее от конкретной решаемой задачи, цикл обработки прекращается.

Реализация алгоритма фейерверков

Определение метода Solve начинается так:

public static double[] Solve(int dim, int n, int maxEpochs)
{
  int m = n * 10;
  int mHat = 5;
  double a = 0.04;
  double b = 0.8;
  int A = 40;
  double minX = -10.0;
  double maxX = 10.0;
...

Переменная m — это общее количество обычных огоньков, создаваемых для n фейерверков. Переменная mHat (названа в соответствии с научными статьями, где она обозначается строчной буквой m с символом ^) определяет число генерируемых особых гауссовых огоньков. Переменные a и b контролируют минимальное и максимальное количества огоньков для конкретного фейерверка. Переменная A является максимальным радиусом разлета для фейерверка. Переменные minX и maxX содержат наименьшее и наибольшее значения для любого значения в массиве position объекта Info.

Метод Solve создает n начальных фейерверков:

Random rnd = new Random(3);
Info[] fireworks = new Info[n];
for (int i = 0; i < n; ++i)
{
  fireworks[i] = new Info();
  fireworks[i].position = new double[dim];
  for (int j = 0; j < dim; ++j)
    fireworks[i].position[j] = (maxX - minX) * rnd.NextDouble() + minX;
  fireworks[i].error = Error(fireworks[i].position);
}

Random-объект rnd изначально использует зародышевое значение 3 только потому, что это значение дает репрезентативный прогон демонстрационной программы. Каждый из n фейерверков хранится как объект Info в массиве. Код инициализации выбирает случайные значения между minX и maxX для каждой ячейки массива position. В некоторых специфических сценариях машинного обучения предпочтительнее инициализировать исходные фейерверки так, чтобы они далеко отстояли друг от друга.

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

double[] bestPosition = new double[dim];
for (int k = 0; k < dim; ++k)
  bestPosition[k] = fireworks[0].position[k];
double bestError = fireworks[0].error; // произвольно
List<Info>[] sparksList = new List<Info>[n];
for (int i = 0; i < n; ++i)
  sparksList[i] = new List<Info>();

Основной цикл обработки начинается следующим образом:

int epoch = 0;
while (epoch < maxEpochs)
{
  if (epoch % 100 == 0) // Отображаем прогресс через каждые 100 итераций
  {
    Console.Write("epoch = " + epoch);
    Console.WriteLine(" error at best known position = " +
      bestError.ToString("F4"));
  }
...

Здесь прогресс отображается через каждые 100 итераций. Возможно, вы захотите передавать булев флаг — переменную с именем progress для управления тем, будет ли вообще отображаться прогресс:

int[] numberSparks = NumberSparks(fireworks, m, a, b);
double[] amplitudes = Amplitudes(fireworks, A, epoch, 
  maxEpochs, minX, maxX);
for (int i = 0; i < n; ++i)
  sparksList[i].Clear(); // число измененных огоньков

Далее с помощью вспомогательных методов NumberSparks и Amplitudes вычисляются количество огоньков в каждом фейерверке и максимальный радиус разлета каждого фейерверка. Число огоньков в фейерверке пропорционально ошибке фейерверка, чтобы хорошие фейерверки большую долю общего количества обычных огоньков — m. Аналогично каждый радиус разлета тоже пропорционален ошибке, чтобы у хороших фейерверков был меньший радиус разлета.

Потом создается экземпляр каждого огонька:

for (int i = 0; i < n; ++i)
{
  double amp = amplitudes[i]; // Радиус разлета текущего фейерверка
  int ns = numberSparks[i];   // Количество огоньков для текущего фейерверка
  for (int j = 0; j < ns; ++j) // Цикл по каждому огоньку текущего фейерверка
  {
    Info spark = new Info(); // у огонька есть позиция и ошибка
    spark.position = new double[dim]; // Выделение пространства (конструктор этого не делает)
    for (int k = 0; k < dim; ++k) // Позиция огонька в границах его родительского фейерверка
      spark.position[k] = fireworks[i].position[k];

Позиция огонька основана на позиции его фейерверка. Затем случайным образом выбирается несколько (z) измерений из массива position, а вместо них в этот массив добавляется случайно формируемая замена:

int z = (int)Math.Round(dim * rnd.NextDouble());
int[] dimensions = PickDimensions(dim, z, epoch);
for (int ii = 0; ii < dimensions.Length; ++ii)
{
  double h = amp * 2 * rnd.NextDouble() - 1;
  int k = dimensions[ii]; // для удобства
    spark.position[k] += h; // замена из родителя
  if (spark.position[k] < minX || spark.position[k] > maxX)
    spark.position[k] = (maxX - minX) * rnd.NextDouble() + minX;
}
spark.error = Error(spar.position);
sparksList[i].Add(spark);

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

if (spark.error < bestError)
    {
      bestError = spark.error;
      for (int k = 0; k < dim; ++k)
        bestPosition[k] = spark.position[k];
    }
  } // каждый новый обычный огонек
} // каждый фейерверк

Теперь генерируются гауссовы огоньки:

AddGaussianSparks(fireworks, sparksList, dim, mHat,
  epoch, minX, maxX, bestPosition, ref bestError, rnd);

Вспомогательный метод AddGaussianSparks генерирует особые огоньки так, чтобы их позиции имели тенденцию располагаться между случайно выбранным фейерверком и лучшей из известных позиций. Метод Solve завершает работу нахождением лучшего и худшего из сгенерированных огоньков и использованием их позиций для двух новых фейерверков на следующей итерации. Остальные n–2 фейерверка создаются, используя случайно выбираемые огоньки:

...
    // Находим лучший и худший огоньки
    // Создаем два новых фейерверка
    // Создаем остальные n-2 фейерверка

    ++epoch;
  } // основной цикл
  return bestPosition;
} // метод Solve

Детали реализации см. в пакете исходного кода, сопутствующего этой статье.

Несколько замечаний

Оригинальная статья, в которой был представлен алгоритм фейерверков, — «Fireworks Algorithm for Optimization» — написана Y. Tan и Y. Zhu в 2010 г. В этой статье было несколько ошибок и факторов, которые сделали этот алгоритм непригодным для реальной оптимизации. Эти недочеты были быстро замечены другими исследователями. Моя статья основана на научной статье S. Zheng, A. Janecek и Y. Tan «Enhanced Fireworks Algorithm» за 2013 г. Вы можете найти обе статьи в нескольких местах Интернета. Я внес довольно много упрощений и мелких изменений в алгоритмы, представленные в обеих научных статьях. По моему мнению, пока недостаточно научно обоснованных фактов для ответа на вопрос, лучше или хуже оптимизация на основе фейерверков, чем другие алгоритмы оптимизации. Но, безусловно, сам по себе этот алгоритм весьма интересен.


Джеймс Маккафри (Dr. James McCaffrey) работает на Microsoft Research в Редмонде (штат Вашингтон). Принимал участие в создании нескольких продуктов Microsoft, в том числе Internet Explorer и Bing. С ним можно связаться по адресу jammc@microsoft.com.

Выражаю благодарность за рецензирование статьи экспертам Microsoft Тодду Беллоу (Todd Bello), Кеннету Гриффину (Kenneth Griffin), Йонг Лю (Yong Liu), Брайену Меллинджеру (Brian Mellinger) и Есину Саке (Esin Saka).