L-BFGS алгоритм минимизации функции многих переменных

Метод Ньютона - классический метод, который в курсе по методам оптимизации обычно идет первым. Этот музейный экспонат давно уже не применяется на практике, но он послужил основой для создания целого семейства квази-Ньютоновских методов, представителем которых является L-BFGS алгоритм.

Квази-Ньютоновские методы: принцип работы

Классический метод Ньютона ипользует гессиан функции. Шаг метода определяется, как произведение матрицы, обратной к гессиану, на градиент функции. Если функция является положительно определенной квадратичной формой, то за один шаг данного метода мы окажемся в её минимуме. В случае знаконеопределенной квадратичной формы, у которой нет минимума, мы сойдемся к седловой точке или к максимуму. Одним словом, метод ищет стационарную точку квадратичной формы.

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

Квази-Ньютоновские методы решают эту проблему следующим образом: вместо гессиана используется его положительно определенная аппроксимация. Если гессиан положительно определен, то мы совершаем шаг по методу Ньютона. Если гессиан знаконеопределен, то перед совершением шага по методу Ньютона мы модифицируем гессиан так, чтобы он был положительно определен.

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

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

L-BFGS схема обновления гессиана

Гессиан функции доступен далеко не всегда, гораздо чаще мы можем вычислить только градиент функции. Поэтому используют следующую схему работы: на основе N последовательных вычислений градиента строится гессиан функции и совершается квази-Ньютоновский шаг. Существует специальная формула, позволяющая итеративно получать аппроксимацию гессиана, причем на каждом шаге аппроксимирующая матрица остается положительно определенной. В данном алгоритме используется BFGS-схема обновления, названная по первым буквам имен Broyden-Fletcher-Goldfarb-Shanno (если быть точным, то эта формула строит не сам гессиан, а обратную к нему матрицу; таким образом, не надо тратить время на её обращение).

Буква L в названии схемы происходит от слов "limited memory". В случае больших размерностей объем памяти порядка N², требуемый для хранения гессиана, оказывается слишком большой нагрузкой, также как и затраты машинного времени на его обработку. Поэтому вместо использования N значений градиента для построения гессиана можно обойтись меньшим числом значений, позволяющим использовать объем памяти порядка N·M. Обычно на практике M выбирают в промежутке от 3 до 7, в сложных случаях можно увеличить эту константу до 20. Разумеется, в результате такой экономии мы получим не сам гессиан, а лишь его аппроксимацию. С одной стороны, при этом замедляется сходимость. С другой, скорость работы может даже вырасти. На первый взгляд парадоксальное, это утверждение не содержит противоречий: сходимость измеряется числом итераций алгоритма, в то время, как скорость работы - числом тактов процессора, потраченных на вычисления.

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

Разностные схемы и аналитический градиент

Если известен градиент функции, то алгоритму требуется намного меньше итераций для сходимости, чем методам, не использующим информацию о градиенте. Одно значение градиента в плане информативности эквивалентно N значениям функции, так что такое различие в быстродействии вполне объяснимо. Вместе с тем, многое зависит от того, как именно вычисляется градиент. Если градиент вычисляется по разностной схеме, то уменьшение числа итераций будет компенсировано пропорциональным ростом их трудоемкости из-за использования разностной схемы. Если градиент известен в аналитической форме и эффективно вычисляется, то L-BFGS алгоритм будет значительно быстрее.

Замечание #1
Не вычисляйте градиент функции на основе двухточечной разностной формулы - она недостаточно точна. В ряде случаев алгоритм просто не сможет работать, и завершится с сообщение об ошибке. Используйте хотя бы четырехточечную формулу.

Использование алгоритма и обратная коммуникация

Алгоритм оптимизации в ходе своей работы должен получать значения функции/градиента в выбранных им точках. В большинстве программных пакетов эта проблема решается путем передачи указателя на функцию (C++, Delphi) или делегата (C#), который осуществляет эту операцию.

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

Таким образом, работа с алгоритмом оптимизации осуществляется в следующей последовательности:

Подготовка структуры данных MinLBFGSState при помощи подпрограммы инициализации алгоритма MinLBFGSCreate.
Установка дополнительных параметров алгоритма оптимизации при помощи подпрограмм MinLBFGSSetXXXXX.
Вызов подпрограммы MinLBFGSIteration.
Если подпрограмма вернула False, работа алгоритма завершена и минимум найден (сам минимум может быть получен при помощи подпрограммы MinLBFGSResults).
Если подпрограмма вернула True, подпрограмма требует информацию о функции и градиенте. Вычислите функцию/градиент и поместите в структуру MinLBFGSState.
После того, как вся требуемая информация загружена в структуру MinLBFGSState, требуется повторно вызвать подпрограмму MinLBFGSIteration.