Алгоритм Бройдена–Флетчера–Гольдфарба–Шанно

Метод оптимизации

В числовой оптимизации алгоритм Бройдена–Флетчера–Гольдфарба–Шанно ( BFGS ) является итеративным методом решения задач нелинейной оптимизации без ограничений. ^[1] Как и связанный с ним метод Дэвидона–Флетчера–Пауэлла , BFGS определяет направление спуска , предварительно обуславливая градиент информацией о кривизне. Он делает это путем постепенного улучшения приближения к матрице Гессе функции потерь , полученной только из оценок градиента (или приближенных оценок градиента) с помощью обобщенного метода секущих . ^[2]

Поскольку обновления матрицы кривизны BFGS не требуют обращения матрицы , ее вычислительная сложность составляет всего , по сравнению с методом Ньютона . Также широко используется L-BFGS , который является версией BFGS с ограниченной памятью, которая особенно подходит для задач с очень большим количеством переменных (например, >1000). Вариант BFGS-B обрабатывает простые ограничения ящика. ^[3] . Матрица BFGS также допускает компактное представление , что делает ее более подходящей для больших задач с ограничениями. ${\displaystyle {\mathcal {O}}(n^{2})}$ ${\displaystyle {\mathcal {O}}(n^{3})}$

Алгоритм назван в честь Чарльза Джорджа Бройдена , Роджера Флетчера , Дональда Голдфарба и Дэвида Шенно . ^{[4] [5] [6] [7]}

Обоснование

Задача оптимизации заключается в минимизации , где — вектор в , а — дифференцируемая скалярная функция. Нет ограничений на значения, которые могут принимать. ${\displaystyle f(\mathbf {x} )}$ ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbb {R} ^{n}}$ ${\displaystyle f}$ ${\displaystyle \mathbf {x} }$

Алгоритм начинается с начальной оценки оптимального значения и продолжается итеративно для получения более точной оценки на каждом этапе. ${\displaystyle \mathbf {x} _{0}}$

Направление поиска p _k на этапе k задается решением аналога уравнения Ньютона:

${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k}),}$

где — приближение к матрице Гессе в , которая обновляется итеративно на каждом этапе, а — градиент функции, вычисленной в x _k . Затем используется линейный поиск в направлении p _k для нахождения следующей точки x _{k +1} путем минимизации по скаляру ${\displaystyle B_{k}}$ ${\displaystyle \mathbf {x} _{k}}$ ${\displaystyle \nabla f(\mathbf {x} _{k})}$ ${\displaystyle f(\mathbf {x} _{k}+\gamma \mathbf {p} _{k})}$ ${\displaystyle \gamma >0.}$

Квази-Ньютоновское условие, налагаемое на обновление, равно ${\displaystyle B_{k}}$

${\displaystyle B_{k+1}(\mathbf {x} _{k+1}-\mathbf {x} _{k})=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k}).}$

Пусть и , тогда удовлетворяет ${\displaystyle \mathbf {y} _{k}=\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}$ ${\displaystyle \mathbf {s} _{k}=\mathbf {x} _{k+1}-\mathbf {x} _{k}}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$,

что является уравнением секанса.

Условие кривизны должно быть выполнено для того, чтобы быть положительно определенным, что можно проверить, предварительно умножив уравнение секущей на . Если функция не является сильно выпуклой , то условие должно быть выполнено явно, например, путем нахождения точки x _{k +1 ,} удовлетворяющей условиям Вульфа , которые влекут за собой условие кривизны, с помощью линейного поиска. ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}>0}$ ${\displaystyle B_{k+1}}$ ${\displaystyle \mathbf {s} _{k}^{T}}$

Вместо того чтобы требовать вычисления полной матрицы Гессе в точке как , приближенный Гессиан на этапе k обновляется путем добавления двух матриц: ${\displaystyle \mathbf {x} _{k+1}}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}=B_{k}+U_{k}+V_{k}.}$

Оба и являются симметричными матрицами ранга один, но их сумма является матрицей обновления ранга два. Матрицы обновления BFGS и DFP отличаются от своего предшественника матрицей ранга два. Другой более простой метод ранга один известен как симметричный метод ранга один, который не гарантирует положительной определенности . Чтобы сохранить симметрию и положительную определенность , форму обновления можно выбрать как . Налагая условие секущей, . Выбирая и , мы можем получить: ^[8] ${\displaystyle U_{k}}$ ${\displaystyle V_{k}}$ ${\displaystyle B_{k+1}}$ ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$ ${\displaystyle B_{k+1}\mathbf {s} _{k}=\mathbf {y} _{k}}$ ${\displaystyle \mathbf {u} =\mathbf {y} _{k}}$ ${\displaystyle \mathbf {v} =B_{k}\mathbf {s} _{k}}$

${\displaystyle \alpha ={\frac {1}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}},}$

${\displaystyle \beta =-{\frac {1}{\mathbf {s} _{k}^{T}B_{k}\mathbf {s} _{k}}}.}$

Наконец, подставляем и и получаем уравнение обновления : ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle B_{k+1}=B_{k}+\alpha \mathbf {u} \mathbf {u} ^{\top }+\beta \mathbf {v} \mathbf {v} ^{\top }}$ ${\displaystyle B_{k+1}}$

${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}.}$

Алгоритм

Рассмотрим следующую задачу безусловной оптимизации, где — нелинейная целевая функция. $${\displaystyle {\begin{aligned}{\underset {\mathbf {x} \in \mathbb {R} ^{n}}{\text{minimize}}}\quad &f(\mathbf {x} ),\end{aligned}}}$$ ${\displaystyle f:\mathbb {R} ^{n}\to \mathbb {R} }$

Из начального предположения и начального предположения матрицы Гессе следующие шаги повторяются по мере сходимости к решению: ${\displaystyle \mathbf {x} _{0}\in \mathbb {R} ^{n}}$ ${\displaystyle B_{0}\in \mathbb {R} ^{n\times n}}$ ${\displaystyle \mathbf {x} _{k}}$

1. Получите направление, решив . ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle B_{k}\mathbf {p} _{k}=-\nabla f(\mathbf {x} _{k})}$
2. Выполнить одномерную оптимизацию ( поиск по линии ) для нахождения приемлемого размера шага в направлении, найденном на первом шаге. Если выполняется точный поиск по линии, то . На практике обычно достаточно неточного поиска по линии с приемлемым удовлетворяющим условиям Вульфа . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Установить и обновить . ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$.
5. ${\displaystyle B_{k+1}=B_{k}+{\frac {\mathbf {y} _{k}\mathbf {y} _{k}^{\mathrm {T} }}{\mathbf {y} _{k}^{\mathrm {T} }\mathbf {s} _{k}}}-{\frac {B_{k}\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} }B_{k}^{\mathrm {T} }}{\mathbf {s} _{k}^{\mathrm {T} }B_{k}\mathbf {s} _{k}}}}$.

Сходимость можно определить, наблюдая за нормой градиента; при наличии некоторого можно остановить алгоритм, когда Если инициализировано с помощью , первый шаг будет эквивалентен градиентному спуску , но дальнейшие шаги все больше и больше уточняются с помощью , приближения к гессиану. ${\displaystyle \epsilon >0}$ ${\displaystyle ||\nabla f(\mathbf {x} _{k})||\leq \epsilon .}$ ${\displaystyle B_{0}}$ ${\displaystyle B_{0}=I}$ ${\displaystyle B_{k}}$

Первый шаг алгоритма выполняется с использованием обратной матрицы , которую можно эффективно получить, применив формулу Шермана–Моррисона к шагу 5 алгоритма, что дает ${\displaystyle B_{k}}$

${\displaystyle B_{k+1}^{-1}=\left(I-{\frac {\mathbf {s} _{k}\mathbf {y} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)B_{k}^{-1}\left(I-{\frac {\mathbf {y} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}\right)+{\frac {\mathbf {s} _{k}\mathbf {s} _{k}^{T}}{\mathbf {y} _{k}^{T}\mathbf {s} _{k}}}.}$

Это можно эффективно вычислить без временных матриц, распознавая, что является симметричным, и что и являются скалярами, используя такое расширение, как ${\displaystyle B_{k}^{-1}}$ ${\displaystyle \mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k}}$ ${\displaystyle \mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}$

${\displaystyle B_{k+1}^{-1}=B_{k}^{-1}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {B_{k}^{-1}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }B_{k}^{-1}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}.}$

Поэтому, чтобы избежать инверсии матрицы, можно аппроксимировать обратную матрицу Гессе вместо самой матрицы Гессе: ^[9] ${\displaystyle H_{k}{\overset {\operatorname {def} }{=}}B_{k}^{-1}.}$

Из начального предположения и приблизительной обратной матрицы Гессе следующие шаги повторяются по мере сходимости к решению: ${\displaystyle \mathbf {x} _{0}}$ ${\displaystyle H_{0}}$ ${\displaystyle \mathbf {x} _{k}}$

1. Получите направление, решив . ${\displaystyle \mathbf {p} _{k}}$ ${\displaystyle \mathbf {p} _{k}=-H_{k}\nabla f(\mathbf {x} _{k})}$
2. Выполнить одномерную оптимизацию ( поиск по линии ) для нахождения приемлемого размера шага в направлении, найденном на первом шаге. Если выполняется точный поиск по линии, то . На практике обычно достаточно неточного поиска по линии с приемлемым удовлетворяющим условиям Вульфа . ${\displaystyle \alpha _{k}}$ ${\displaystyle \alpha _{k}=\arg \min f(\mathbf {x} _{k}+\alpha \mathbf {p} _{k})}$ ${\displaystyle \alpha _{k}}$
3. Установить и обновить . ${\displaystyle \mathbf {s} _{k}=\alpha _{k}\mathbf {p} _{k}}$ ${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+\mathbf {s} _{k}}$
4. ${\displaystyle \mathbf {y} _{k}={\nabla f(\mathbf {x} _{k+1})-\nabla f(\mathbf {x} _{k})}}$.
5. ${\displaystyle H_{k+1}=H_{k}+{\frac {(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}+\mathbf {y} _{k}^{\mathrm {T} }H_{k}\mathbf {y} _{k})(\mathbf {s} _{k}\mathbf {s} _{k}^{\mathrm {T} })}{(\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k})^{2}}}-{\frac {H_{k}\mathbf {y} _{k}\mathbf {s} _{k}^{\mathrm {T} }+\mathbf {s} _{k}\mathbf {y} _{k}^{\mathrm {T} }H_{k}}{\mathbf {s} _{k}^{\mathrm {T} }\mathbf {y} _{k}}}}$.

В задачах статистической оценки (таких как максимальное правдоподобие или байесовский вывод) достоверные интервалы или доверительные интервалы для решения могут быть оценены из обратной конечной матрицы Гессе ^{[ требуется ссылка ]} . Однако эти величины технически определяются истинной матрицей Гессе, и приближение BFGS может не сходиться к истинной матрице Гессе. ^[10]

Дальнейшее развитие событий

Формула обновления BFGS в значительной степени опирается на то, что кривизна строго положительна и отделена от нуля. Это условие выполняется, когда мы выполняем линейный поиск с условиями Вульфа на выпуклой цели. Однако некоторые реальные приложения (например, методы последовательного квадратичного программирования) регулярно выдают отрицательные или почти нулевые кривизны. Это может произойти при оптимизации невыпуклой цели или при использовании подхода доверительной области вместо линейного поиска. Также возможно получение ложных значений из-за шума в цели. ${\displaystyle \mathbf {s} _{k}^{\top }\mathbf {y} _{k}}$

В таких случаях можно использовать одно из так называемых затухающих обновлений BFGS (см. ^[11] ), которые модифицируют и/или для получения более надежного обновления. ${\displaystyle \mathbf {s} _{k}}$ ${\displaystyle \mathbf {y} _{k}}$

Известные реализации

Известные реализации с открытым исходным кодом:

• ALGLIB реализует BFGS и его версию с ограниченным объемом памяти на C++ и C#
• GNU Octavefsolve использует в своей работе разновидность BFGS с расширениями доверенных областей .
• GSL реализует BFGS как gsl_multimin_fdfminimizer_vector_bfgs2. ^[12]
• В R алгоритм BFGS (и версия L-BFGS-B, которая допускает ограничения по ящикам) реализован как опция базовой функции optim(). ^[13]
• В SciPy функция scipy.optimize.fmin_bfgs реализует BFGS. ^[14] Также возможно запустить BFGS с использованием любого из алгоритмов L-BFGS , установив параметр L на очень большое число.
• В Julia пакет Optim.jl реализует BFGS и L-BFGS в качестве опции решателя для функции optimize() (среди прочих опций). ^[15]

Известные фирменные реализации включают в себя:

• Программное обеспечение для крупномасштабной нелинейной оптимизации Artelys Knitro реализует, среди прочего, алгоритмы BFGS и L-BFGS.
• В MATLAB Optimization Toolbox функция fminunc использует BFGS с кубическим линейным поиском, когда размер задачи установлен на «средний масштаб».
• Mathematica включает BFGS.
• LS-DYNA также использует BFGS для решения неявных проблем.

Смотрите также

• Алгоритм BHHH
• Формула Дэвидона–Флетчера–Пауэлла
• Градиентный спуск
• L-BFGS
• Алгоритм Левенберга–Марквардта
• Метод Нелдера–Мида
• Поиск по шаблону (оптимизация)
• Квази-ньютоновские методы
• Симметричный ранг один
• Компактное квазиньютоновское представление

Ссылки

1. Флетчер, Роджер (1987), Практические методы оптимизации (2-е изд.), Нью-Йорк: John Wiley & Sons , ISBN 978-0-471-91547-8
2. Деннис, Дж. Э. младший ; Шнабель, Роберт Б. (1983), «Методы секущих для безусловной минимизации», Численные методы безусловной оптимизации и нелинейные уравнения , Энглвуд Клиффс, Нью-Джерси: Prentice-Hall, стр. 194–215, ISBN 0-13-627216-9
3. Берд, Ричард Х.; Лу, Пэйхуан; Нокедаль, Хорхе; Чжу, Цию (1995), «Алгоритм с ограниченной памятью для оптимизации с ограниченными ограничениями», SIAM Journal on Scientific Computing , 16 (5): 1190–1208, CiteSeerX 10.1.1.645.5814 , doi : 10.1137/0916069 
4. Бройден, К. Г. (1970), «Сходимость класса алгоритмов минимизации двойного ранга», Журнал Института математики и ее приложений , 6 : 76–90, doi : 10.1093/imamat/6.1.76
5. Флетчер, Р. (1970), «Новый подход к алгоритмам с переменной метрикой», Computer Journal , 13 (3): 317–322, doi : 10.1093/comjnl/13.3.317
6. Голдфарб, Д. (1970), «Семейство переменных метрических обновлений, полученных с помощью вариационных средних», Математика вычислений , 24 (109): 23–26, doi : 10.1090/S0025-5718-1970-0258249-6
7. Шэнно, Дэвид Ф. (июль 1970 г.), «Обусловленность квазиньютоновских методов минимизации функций», Mathematics of Computation , 24 (111): 647–656, doi : 10.1090/S0025-5718-1970-0274029-X , MR  0274029
8. Флетчер, Роджер (1987), Практические методы оптимизации (2-е изд.), Нью-Йорк: John Wiley & Sons , ISBN 978-0-471-91547-8
9. Нокедаль, Хорхе; Райт, Стивен Дж. (2006), Численная оптимизация (2-е изд.), Берлин, Нью-Йорк: Springer-Verlag , ISBN 978-0-387-30303-1
10. Ge, Ren-pu; Powell, MJD (1983). "Сходимость матриц переменной метрики в неограниченной оптимизации". Математическое программирование . 27 (2). 123. doi :10.1007/BF02591941. S2CID  8113073.
11. Хорхе Нокедаль; Стивен Дж. Райт (2006), Численная оптимизация
12. "GNU Scientific Library — документация GSL 2.6". www.gnu.org . Получено 22.11.2020 .
13. "R: Универсальная оптимизация". stat.ethz.ch . Получено 22.11.2020 .
14. "scipy.optimize.fmin_bfgs — Справочное руководство SciPy v1.5.4". docs.scipy.org . Получено 22.11.2020 .
15. "Optim.jl Настраиваемые параметры". julianlsolvers .

Дальнейшее чтение

• Авриэль, Мордехай (2003), Нелинейное программирование: анализ и методы , Dover Publishing, ISBN 978-0-486-43227-4
• Боннан, Ж. Фредерик; Жильбер, Ж. Шарль; Лемарешаль, Клод ; Сагастисабаль, Клаудия А. (2006), «Ньютоновские методы», Численная оптимизация: теоретические и практические аспекты (второе изд.), Берлин: Springer, стр. 51–66, ISBN 3-540-35445-X
• Флетчер, Роджер (1987), Практические методы оптимизации (2-е изд.), Нью-Йорк: John Wiley & Sons , ISBN 978-0-471-91547-8
• Луенбергер, Дэвид Г .; Йе, Инью (2008), Линейное и нелинейное программирование , Международная серия по исследованию операций и науке управления, т. 116 (третье изд.), Нью-Йорк: Springer, стр. xiv+546, ISBN 978-0-387-74502-2, г-н  2423726
• Келли, CT (1999), Итерационные методы оптимизации , Филадельфия: Общество промышленной и прикладной математики, стр. 71–86, ISBN 0-89871-433-8