Квазиньютоновские методы, или когда вторых производных для Атоса слишком много

При первом знакомстве с квазиньютоновскими методами можно удивиться дважды. Во-первых, после беглого взгляда на формулы охватывают сомнения, что это вообще может работать. Однако же они работают. Дальше кажется сомнительным, что они будут работать хорошо. И тем удивительнее видеть то, насколько они превосходят по скорости разнообразные вариации градиентного спуска, причем не на специально построенных задачах, а на самых настоящих, взятых из практики. И если после этого еще остаются сомнения вперемешку с интересом — то нужно разбираться в том, почему вообще работает это нечто.
Ранее уже рассматривалось происхождение и основные идеи, которые приводят в движение градиентные методы, включая метод Ньютона. А именно, мы опирались на ту информацию о поведении функции в окрестности текущего положения, которую дает нам незамысловатый математический анализ. Как минимум предполагалось, что информация о первых производных нам доступна. Что делать, если это все, что нам доступно? Градиентный спуск — наш приговор? Конечно да, если только вдруг не вспомнить, что мы имеем дело с процессом, в ходе которого целевая функция как следует прошаривается. И раз так, что почему бы нам не использовать накопленную информацию о поведении функции для того, чтобы сделать наше блуждание по ее поверхности чуть менее слепым?

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

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

Методы секущих


Напомню, что метод Ньютона для решения системы уравнений 0a0ec780406efe57ca6444290ccfde09.gif, основан на замене в окрестности некоторой близкой к решению точки 7790dd0efb4a03a4c876741804d9b559.gif функции 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif ее линейной аппроксимацией d15479f235f0d60ce8837c9043a0d2cc.gif, где 206f349991c0724c2fdce788124abe1c.gif — линейный оператор, который в случае, когда 7790dd0efb4a03a4c876741804d9b559.gif является вектором и 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif имеет частные производные по каждой переменной, совпадает с матрицей Якоби 4d2826ff6ba22f9f67cab70bfbe17a16.gif. Далее решается уравнение c9d8be3f2d70054db890ea34e3409544.gif и точка 2c4a7be5582848bfbcdd9ee141e7d764.gif принимается в качестве нового приближения к искомому решению. Это просто и это работает.

Но что делать, если матрицу Якоби мы по каким-либо причинам вычислить не можем? Первое, что в данном случае приходит в голову — если мы не можем вычислить частные производные аналитически, то вполне можем получить для них численное приближение. Самым простым (хотя и далеко не единственным) вариантом такого приближения может служить формула правых конечных разностей: 149708f5b8bab4374023295557622e82.gif, где 459e61aa08f7fe807167a596e7ebd8a9.gif — j-й базисный вектор. Матрицу, составленную из таких приближений, будем обозначать a5153399048e881eb8661304792b8c81.gif. Анализу того, насколько замена 206f349991c0724c2fdce788124abe1c.gif на a5153399048e881eb8661304792b8c81.gif в методе Ньютона влияет его сходимость, посвящено достаточно большое количество работ, но нас в данном случае интересует другой аспект. А именно, такое приближение требует вычисления функции в N дополнительных точках, и, кроме того, функция 6c4afd1002ddcfa43d07afbc9f103a9d.gif в этих точках интерполирует функцию 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif, т.е.

c565c94b4a37b9cd5f42fc1be92b2e15.gif

Не каждая аппроксимация матрицы Якоби обладает таким свойством, но каждая матрица аффинной функции, обладающей таким свойством, является аппроксимацией матрицы Якоби. Действительно, если 88f5c8dd7e9876a2d0e0980882f261da.gif и 5ade1ead2804a3bfaa8ffdf9122a179a.gif, то при 803ca4351b87edf1a13a2a2947772fa7.gif. Данное свойство, а именно — свойство интерполирования, дает нам конструктивный способ обобщить метода Ньютона.

Пусть 194ad1d42aa4320679b9498748ceb78d.gif — функция, удовлетворяющая требованию 06e34e3a7a0d0058ef351da74258a637.gif для некоторой системы линейно независимых векторов cf2deb64e8b0e4d34902a32a5fd93b7b.gif. Тогда такая функция называется секущей функции 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif, а определяющее ее уравнение — уравнением секущей. Если при этом система векторов cf2deb64e8b0e4d34902a32a5fd93b7b.gif полна (то есть их ровно N и они все еще линейно независимы), и, кроме того, система векторов 94fcf55798902795ffb670e35359d2af.gif линейно независима, то 77feecdd2ae9a4795d2f81f3eec18b1b.gif определяется однозначно.

Любой метод, построенный на основе локальной замены уравнения 0a0ec780406efe57ca6444290ccfde09.gif уравнением вида ccb55780f0c8e65187b0f4c9126be81c.gif, где 77feecdd2ae9a4795d2f81f3eec18b1b.gif удовлетворяет уравнению секущих, называют методом секущих.

Возникает справедливый вопрос о том, как наиболее рациональным способом построить секущую для функции 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif. Кажется очевидным следующий ход рассуждений: пусть в точке x уже построена аффинная модель, которая интерполирует заданную функцию в точках 462bcff32469c0ec5f8ccfc80534c05c.gif. Решение уравнения ccb55780f0c8e65187b0f4c9126be81c.gif дает нам новую точку 2c4a7be5582848bfbcdd9ee141e7d764.gif. Тогда для построения аффинной модели в точке 787cf7c3a3d374114b3a07305b7fa446.gif разумнее всего выбрать точки интерполяции так, что в них значение 01aa158fc8bc3d7f7f3b2807df8b4a5e.gif уже известно — то есть взять их из множества 33329722533f0b608b0994d2a5ba83fa.gif. Возможны разные варианты того, какие именно точки выбирать из множества ранее использованных. Например, можно взять в качестве точек интерполяции такие, в которых bb3e3acd5043b859fe89006d4cabe5a0.gif имеет наименьшее значение, либо просто первые 0558e93d918ff32e873b6a71703e9969.gif точек. В любом случае, кажется очевидным, что 95f75692ba0aeefcef24ae42714dbc1b.gif следует включать в множество точек интерполяции для новой аффинной модели. Таким образом, за f248e891effc6650d9d31fbefc54cbe4.gif шагов итерационного процесса в нашем множестве может оказаться до f248e891effc6650d9d31fbefc54cbe4.gif смещений, построенных по ранее пройденным точкам. Если процесс построен таким образом, что новая аффинная модель использует не более 4b28c13d5f5d658adb7478fbc9efc923.gif предыдущих значений, то такой процесс называют p-точечным методом секущих.

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

Вообще, самым устойчивым оказывается 2-х точечный метод секущих. То есть метод, в котором на каждой итерации нам придется вычислять дополнительно N-1 значений функции. Это явно не подходит для наших практических целей.

Тогда вопрос —, а к чему все это было?

Квазиньютоновские методы решения уравнений

Выход из положения прост, хотя и не очевиден. Если у нас нет технической возможности на основании уже вычисленных значений однозначно определить аффинную модель, удовлетворяющую уравнению секущих — то и не надо. Возьмем уравнение секущих за основу, но будем требовать, чтобы оно выполнялось только для некоторой неполной системы векторов e5ef2b43292735aa2a68afffb80bf520.gif. Иными словами, мы будем требовать выполнение условия интерполяции только для достаточно малого количества известных значений. Разумеется, в этом случае мы уже не можем гарантировать, что используемая в такой модели матрица будет стремиться к матрице Якоби, однако это нам и не потребуется. Добавив к этому, что аффинная модель должна интерполировать функцию в текущей точке, то есть 3b99d17fa378aeaf36097faef3830bd5.gif, мы получаем следующую формулировку метода секущих:

78cf304219a4bc90d4f900062bf2d027.gif

Бройден был первым, кто рассмотрел методы такого вида при m=1, назвав их квазиньютоновскими. Понятно, что условие секущих в данном случае позволяет однозначно идентифицировать матрицу c9d999d9a4e8bd3d6f8e50519d1dfaa8.gif только если наложить на нее дополнительные условия, и каждое такое дополнительное условие порождает отдельный метод. Сам Бройден рассуждал следующим образом:

поскольку движение в направлении 4b28c13d5f5d658adb7478fbc9efc923.gif от точки 46082f7d6471c3fabb832d8f94075758.gif к точке 1d056f3016bc715aacc23418d8629173.gif не дает нам никакой дополнительной информации о том, как изменяется функция в отличных от 4b28c13d5f5d658adb7478fbc9efc923.gif направлениях, то воздействие новой аффинной функции на вектор 9fcc76a21130891ea5d5b10efa979bff.gif должно отличаться от воздействия старой функции на тот же вектор тем меньше, чем больше отличается 9fcc76a21130891ea5d5b10efa979bff.gif от 4b28c13d5f5d658adb7478fbc9efc923.gif. В крайнем случае, когда 9fcc76a21130891ea5d5b10efa979bff.gif ортогонален 4b28c13d5f5d658adb7478fbc9efc923.gif, поведение новой функции вообще не должно отличаться от поведения старой.

Идея Бройдена гениальна в своей простоте. Действительно, если мы не имеем новой информации о поведении функции, то лучшее, что мы можем сделать — постараться не профукать старую. Тогда дополнительное условие

0b3d43d207b144b926274d0c81abccbf.gif для всех 9fcc76a21130891ea5d5b10efa979bff.gif, таких, что c169f6315171249a34b50b26a2975c6e.gif
позволяет однозначно определить матрицу нового преобразования — она получается добавлением к старой матрице поправки ранга 1.

522f361f94a7b5a2e9da68094983b21d.gif

Однако, не смотря на простоту и логичность сделанных Бройденом умозаключений, они не дают той точки опоры, которая могла бы послужить основой для построения других аналогичных методов. К счастью, существует более формальное выражение его идеи. А именно, построенная таким образом матрица 1477e7ca06155c3e43fd4a640e0f7f98.gif оказывается решением следующей задачи:

02896f3facb898d70f26abad02fe90a9.gif

Ограничение задачи представляет собой ни что иное, как уравнение секущей, а условие минимизации отражает наше желание сохранить как можно больше информации, содержащейся в матрице 107a45803b226180325815eaa7be8706.gif. Мерой расхождения между матрицами в данном случае выступает норма Фробениуса, в которой поставленная задачи имеет однозначное решение. Эта формулировка уже вполне может послужить отправной точкой для построения других методов. А именно, мы можем изменить как меру, посредством которой оцениваем вносимые изменения, так и ужесточить условия, накладываемые на матрицу. В общем, с такой формулировкой метода уже можно работать.

Квазиньютоновские методы оптимизации

Поняв основную идею, мы наконец можем вернуться к задачам оптимизации и заметить, что применение формулы Бройдена для пересчета аффинной модели, не слишком удачно ложится на нашу задачу. В самом деле, первая производная от функции-градиента 6b882ebe727121dcb5fc21b091044b5a.gif есть ни что иное, как матрица Гессе, которая по построению является симметричной. В то же время, обновление по правилу Бройдена приводит к несимметричной матрице 1477e7ca06155c3e43fd4a640e0f7f98.gif даже если 107a45803b226180325815eaa7be8706.gif была симметричной. Это не значит, что метод Бройдена не может быть применен для решения уравнения стационарной точки, но на основании такого правила обновления мы вряд ли сможем построить хорошие методы оптимизации. Вообще достаточно очевидно, что квазиньютоновский метод должен работать тем лучше, чем точнее система условий задачи описывает специфику конкретной матрицы Якоби.

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

03e167aa25e0f4aa6b8df8546552e79a.gif

Решением данной задачи является

df835674b94bab190bca3c18efed98ce.gif

Здесь 89cf84f292ed72d1b20755677688a054.gif, а формула пересчета матрицы носит имя своих создателей — Пауэлла, Шанно и Бройдена (PSB). Получаемая в результате матрица симметрична, но явно не является положительно определенной, если только вдруг 6c704047d3148fd7a8b563aaf79dd7f4.gif не окажется коллинеарным 4b28c13d5f5d658adb7478fbc9efc923.gif. А мы видели, что положительная определенность является крайне желательным в методах оптимизации.

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

f37eb10f4eb10d0c54acc9adab962f10.gif

Происхождение такой постановки вопроса — отдельная большая тема, однако интересно, что если матрица T такова, что 673872131fa6cb0f44e6839be0e448e7.gif (то есть G также является матрицей аффинного преобразования, удовлетворяющего уравнению секущих для направления p), то решение данной задачи оказывается независимым от выбора T и приводит к формуле обновления

135ea6c14ea8f63f961e83576f1be5d5.gif

известной как формула Давидона-Флетчера-Пауэлла. Данный способ обновления неплохо зарекомендовал себя на практике, поскольку обладает следующим свойством:

если e3ec441c17e1b43df108a7d8e15d3dd6.gif и 107a45803b226180325815eaa7be8706.gif положительно определена, то 1477e7ca06155c3e43fd4a640e0f7f98.gif также положительно определена.

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

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

3370b0af216ab9695789eeb586cf3604.gif

Ее решение — широко известная формула, открытая почти одновременно Бройденом, Флетчером, Гольдфарбом и Шанно (BFGS).

3fa840c7b6ec3de81eb02bb0e9240722.gif

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

Все описанные методы обновления матрицы предполагают внесение поправки ранга 2. Это позволяет легко и непринужденно обратить матрицу 1477e7ca06155c3e43fd4a640e0f7f98.gif, используя формулу Шермана-Моррисона и значение 5f63ac2d91f47a730fee01b5db38f3bd.gif.

ed09f80027e56f58a3502cc943758509.gif

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

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

Заключение

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

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

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

© Habrahabr.ru