full_book

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

следующая точка y1=y0–f0′/f0′′ равнялась бы (–y0). В силу нечетности первой производной и четности второй производной, обязательно на следующем шаге выполнится равенство y2=–y1=y0. Поэтому сходимости к точке минимума из точки y0 не будет.

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

Предварительно укажем на некоторые свойства выбранного класса функций.

Свойство 7.3. Любая функция f(y) из класса Φm, M имеет единственный минимум.

Это непосредственно следует из сильной выпуклости функции (смотрите пункт 1.4.1.1).

Свойство 7.4. Для функций f(y) из класса Φm,M существует взаимосвязь между ошибкой по координате и ошибкой по значению функции, выражаемая соотношением

0,5m||y–y*||2 ≤ f(y) – f(y*) ≤ 0,5M||y–y*||2.

ДОКАЗАТЕЛЬСЬВО сразу вытекает из (7.22) и разложения f(y)–f(y*) = (y–y*)ТΓ f(θ) (y–y*)/2.

Свойство 7.5. Для функций f(y) из класса Φm,M существует взаимосвязь между ошибкой по значению функции и нормой градиента

m(1+m/M)( f(y) – f(y*) ) ≤ || f(y)||2.

ДОКАЗАТЕЛЬСТВО. Раскладывая f(y*) в окрестности точки y, получим для приращения на классе Φm,M

f(y)–f(y*) = ( f(y), y–y* ) – (y–y*)ТΓ f(θ) (y–y*)/2 ≤ || f(y)|| ||y–y*||–0,5m||y–y*||2.

Усиливая это неравенство слева с использованием свойства 7.4, придем к оценке

||y–y*|| ≤ || f(y)|| / m.

Используя еще раз два предыдущих неравенства, а также правое неравенство из свойства 7.4, получим

f(y)–f(y*) ≤ || f(y)||2 / m – 0,5m||y–y*||2 ≤ || f(y)||2 / m – (m/M)( f(y)–f(y*)).

Отсюда непосредственно вытекает справедливость свойства.

Поведение метода Ньютона с регулировкой шага на классе Φm,M определяется следующей теоремой [39]. Ее доказательство опирается на приведенные выше свойства.

Теорема 7.4 Метод Ньютона с регулировкой шага на функциях f Φm,M для любой начальной точки y0 порождает последовательность точек y k, сходящуюся к точке минимума y* со сверхлинейной скоростью.

ДОКАЗАТЕЛЬСТВО распадается на две части. Во-первых, устанавливается сам факт сходимости (рассуждения здесь во многом аналогичны тем, которые использовались при обосновании сходимости градиентного метода). Во-вторых, оценивается скорость сходимости, исходя из того, что любая начальная точка гарантированно приводит в сколь угодно малую окрестность решения. Это

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

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

Для доказательства сходимости рассмотрим точки, лежащие на луче,

f( y(x)) − f ( yk ) = ( f ( yk ) , y(x) − yk ) + 0,5 ( y(x) − yk )T Γ f (θ) ( y(x) − yk ) ≤

≤x( f ( yk ) , d k ) + x2 0,5M d k 2 =ϕ(x).

= − ( f ( yk ) , d k )2 (2M d k 2 ) =

= − ( f ( yk ) , (Γ f ( yk ))−1 f ( yk ) )( Γ f ( yk )d k , d k ) (2M d k 2 ) ≤ ≤ −(m (2M 2 )) f ( yk ) 2 ≤ 0.

Таким образом, либо норма градиента равна нулю (для функций класса Фm,M это означает точное определение глобального минимума), либо изменение функции на шаге строго отрицательно, т.е. имеет место строгое монотонное убывание последовательности значений функции на траектории поиска. В силу ограниченности снизу значений функции из класса Фm,M, получаем, что

f ( yk+1 ) − f ( yk ) → 0 при k → ∞, но тогда из предыдущего неравенства

|| f (y k)||→ 0 при k→∞. В рассматриваемом классе это означает сходимость по координате (см. доказательство свойства 7.5). Сходимость к решению из любой начальной точки доказана.

Оценим скорость сходимости. Будем использовать сокращенные

место следующие оценки, первая из которых дважды использует свойство 7.4.

)Замечание. Если дополнительно потребовать от функции f Φm,M существования непрерывных третьих производных, можно доказать, что метод Ньютона–Рафсона будет сходиться квадратично.

Таким образом, модификация Ньютона–Рафсона расширяет область сходимости метода Ньютона.

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

7.4.2. Стратегии модификации матриц Гессе при нарушении их положительной определенности

Вторая модификация метода Ньютона связана с преодолением случаев отсутствия положительной определенности матрицы вторых производных. Следует обратить внимание на то, что при нарушении положительной определенности Γk (а значит и (Γ k)–1, если она существует) направление смещения d k = (Γk)–1(– f(y k)) может не быть направлением убывания. Действительно, производная функции f в точке y k по направлению d k оценивается следующим образом

∂f(y k)/∂ d k = ( f(y k),d k )= – ( f(y k),(Γk)–1 f(y k) ).

Врассматриваемом случае знак этого произведения не определен и может

оказаться положительным. В этом случае метод Ньютона–Рафсона применить нельзя, т.к. при его использовании смещение вдоль d k будет нулевым.

ПРОСТОЙ АЛГОРИТМ. Простейший выход из этого положения может состоять в проверке знака скалярного произведения ( f(y k), d k ) на каждом шаге. Если окажется, что ( f(y k), d k )< 0, то выполняется обычный шаг по методу Ньютона– Рафсона, если же ( f(y k), d k )≥ 0, то проводится замена направления поиска на антиградиентное направление d k = – f(y k) и выполняется шаг по методу наискорейшего градиентного спуска.

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

Таким образом, необходима более гибкая модификация метода. Необходимо как-то использовать матрицу Γk, хотя непосредственно ее применить нельзя. Основная идея, на основе которой выполняется модификация матриц, состоит в том, чтобы заменить Γk на достаточно близкую к ней (в смысле некоторой нормы) положительно определенную матрицу Γk и затем использовать ее в итерационном соотношении метода Ньютона–Рафсона

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

Наиболее естественным представляется использование спектрального разложения. Определим для Γk набор собственных чисел λl,…,λN и систему ортонормированных собственных векторов ul,…,uN. Тогда возможно следующее представление

Γk = λ1u 1 (u 1) T +…+ λ N u N(u N)T = ULUT,

где матрица U составлена из вектор–столбцов ul,…,uN, L – диагональная матрица с числами λl,…,λN по диагонали. Такое представление матрицы называется спектральным разложением. Если положительная определенность Γk нарушена, то существует λ i≤ 0. Матрица Γk строится так, что у нее сохраняются все собственные векторы ul,…,uN, а собственные числа заменяются на новые λ i так, что

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

где в диагональной матрице L на диагонали используются числа λl,…, λN.

При этом подпространство локальной положительной кривизны функции f(y) сохраняется, а подпространство отрицательной кривизны функции заменяется подпространством малой положительной кривизны. Если построить квадратичную аппроксимацию функции f(y) по результатам ее испытания в точке yk и затем заменить в ней матрицу вторых производных Γk на модифицированную матрицу (7.24)–(7.25), то произойдут качественные изменения в структуре аппроксимации.

Рис.7. 8. Изменение изолиний при замене матрицы

На рис.7.8 на примере пространства двух переменных показаны изменения в линиях равного уровня квадратичной аппроксимации P k(y) после замещения знаконеопределенной матрицы Γk положительно определенной Γk. На этом рисунке видно, как направление метода Ньютона d k, приводящее в стационарную точку поверхности Pk(x) заменяется новым направлением d k, приводящим в минимум модифицированной аппроксимации. Для создания наглядных представлений об изменении характера аппроксимирующей поверхности при замене Γk на положительно определенную Γk., полезно обратиться к иллюстрациям, представленным на рис. 7.8– 7.9.

Поведение функции f(x), представленное на рис.7.9, соответствует изолиниям, показанным на рис.7.8. На рисунке следует обратить внимание на то, что кривизна модифицированной аппроксимирующей поверхности P k(y), построенной по измененной матрице Γk, рассматриваемая в направлении u1 положительной локальной кривизны поверхности f(y), совпадает с локальной кривизной самой f(y) в этом направлении. В то же время, в направлении u2 отрицательная кривизна заменена малой положительной.

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

Рис. 7. 9. Изменение вида аппроксимирующей поверхности в случае знаконеопределенной матрицы Гессе

Несколько иное соответствие между первоначальной и модифицированной аппроксимациями возникает в том случае, когда матрица Γk отрицательно определена. В этом случае все собственные направления матрицы Γk трансформируются в направления малой положительной кривизны (рис.7.10). Это приводит к замене направления Ньютона d k на новое — d k, ориентированное на точку минимума измененной аппроксимации.

Рис. 7. 10. Влияние модификации матрицы при ее отрицательной определенности

Таким образом, метод коррекции матрицы на основе спектрального разложения весьма нагляден и прост для понимания. Однако этот подход имеет один существенный недостаток – большие затраты по вычислениям, связанные с поиском собственных векторов и чисел симметричной матрицы. Необходимый объем вычислений для построения спектрального разложения оценивается как 4N3.

При разработке вычислительных методов оптимизации для построения положительно определенных матриц Γk по исходным матрицам Γk вместо

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

спектрального разложения часто используется модифицированное разложение Холесского [4], вычислительная реализация которого проще и требует меньшего числа операций.

В теории матриц известно, что для любой симметричной положительно определенной матрицы Γ существует нижняя треугольная матрица L с единичной диагональю и диагональная матрица D с положительной диагональю, что справедливо разложение Холесского Γ = LDL Т, т.е.

для i < j ≤ N.

Из (7.27)–(7.29) легко получить формулы для определения коэффициентов разложения Холесского. Расчет их значений производится строка за строкой для матриц D и LT . Порядок вычисления этих коэффициентов удобно пояснить следующей диаграммой:

d11 l12, l13, l14, …; d22 l23, l24, …; d33 l34, l35, …

Получим

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

Это обеспечит положительность с «запасом» элементов модифицированной диагональной матрицы D. Отметим, что условие (7.32) не может являться единственным условием модификации. Действительно, близость к нулю некоторых модифицированных элементов d i j при их дальнейшем использовании в (7.31) может привести к лавинообразному росту элементов l i j при вычислениях. При обычном разложении Холесского это невозможно, т.к. из (7.28) вытекает, что

требование, чтобы для модифицированных элементов разложения выполнялось

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

Выполним необходимую модификацию за счет изменения только диагональных элементов матрицы Γ. Обозначим через ∆ i добавки к элементам γ i i. Тогда согласно (7.32)–(7.33) должно выполняться:

Элементы (7.36), (7.37) определят модифицированные матрицы LT, D. В

γ * = max { γ ii : i = 1,K, N}, ξ = max { | γ i j | :1 ≤ i ≤ N , 1 ≤ j ≤ N , i ≠ j }. (7.39)

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

АЛГОРИТМ МОДИФИЦИРОВАННОГО МЕТОДА НЬЮТОНА–РАФСОНА.

ШАГ 0. Задаются начальная точка y 0, параметры выбора коэффициента одномерного шага 0 <µ < η <<1, 0<σ<<1; параметр останова ε и параметр

модификации δ >0. Полагается k = 0 .

ШАГ 1. Вычисляются f k = f(y k), f k = f(y k), Γ k = Γ f(y k).

ШАГ 2. Вычисляется β 2 по формулам (7.38), (7.39). Строятся матрицы модифицированного разложения Холесского Lk, Dk для матрицы Γk.

ШАГ 3. Определяется модифицированное направление Ньютоновского шагаd k = ( Γk)–1(– f k) путем последовательного решения двух систем линейных уравнений с треугольными матрицами

L kv k = – f k

( D k ( L k)T) d k = v k

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

)Замечание 1. Известно, что выполнение модифицированного разложения

Холесского требует около (1/6)N 3 операций, а последующее определение d k требует числа операций порядка N 2.

)Замечание 2. Модифицированный метод Ньютона сохраняет поисковые

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

знаконеопределенность матриц Γk ), поскольку новое направление поиска d k, в силу гарантированной положительной определенности матриц Γk, обязательно является направлением строгого локального убывания и соответствует направлению антиградиента в некоторой новой метрике пространства. Кроме того, метод обладает сверхлинейной скоростью сходимости в окрестности решения, если функция в этой окрестности дважды непрерывно дифференцируема и принадлежит классу Фm,M из (7.22).

Справедливость последнего утверждения следует из того, что в указанной области модификация матрицы производиться не будет, т.е. Γk = Γk, а, следовательно, метод будет точно совпадать с методом Ньютона–Рафсона.

7.5.Методы первого порядка, явно изменяющие метрику пространства

Вэтом разделе будет продолжено изучение методов, основанных на квадратичной модели поведения минимизируемой функции. В отличие от предыдущего раздела, будет рассмотрена группа методов, которые хотя и используют предположение о гладкости функции и близости ее к квадратичной, но не измеряют матриц вторых производных. Будет рассмотрено несколько групп таких методов. Методы первой группы (квазиньютоновские методы) строят оценки матриц Гессе и используют их вместо истинных матриц вторых производных, точно также, как это делает метод Ньютона–Рафсона. Методы второй группы — методы растяжения также основаны на построении вспомогательных матриц, используемых для перемасштабирования пространства, но эти матрицы не являются оценками Гессиана и строятся на основе эвристических принципов. Эти методы будут рассмотрены в данном разделе.

Методы третьей группы явно никаких матриц не строят, хотя неявно метрику пространства изменяют (методы сопряженных направлений). Они рассматриваются в разделе 7.6.

7.5.1.Квазиньютоновские методы. Рекуррентные соотношения для оценок матриц Гессе по измерениям градиента в основных точках поиска

Методы этого класса относятся к методам первого порядка. Они используют результаты испытаний, состоящих в вычислении f(y k), f(y k). Предполагается, что функция f обладает свойствами, соответствующими квадратичной модели.

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

Следовательно, у функции f существует симметричная матрица вторых производных, недоступная непосредственному измерению.

Казалось бы, в этих условиях самым естественным являлось конечноразностное оценивание Гессиана в каждой точке y k по измерениям градиента на множестве узлов, размещенных с некоторым шагом h в окрестности данной точки. Получив оценку можно применить модифицированный метод Ньютона– Рафсона. Однако данный подход требует слишком большого объема вычислений и, кроме того, связан со значительными погрешностям оценивания. Это приводит к тому, что данный подход обычно не применяется. Оказывается оценки Гессиана можно строить без дополнительных вычислений градиента функции f(y). Именно такой подход используется в квазиньютоновских методах.

Идея, положенная в основу квазиньютоновских методов, состоит в том, чтобы по результатам измерения градиентов функции f в точках y k траектории поиска попытаться построить симметричную матрицу Gk, являющуюся оценкой кривизны поверхности f(y) на траектории поиска, т.е. оценивающую Γ f, или построить симметричную матрицу Hk, являющуюся оценкой обратной матрицы ( Γ f )–1. После выполнения k-го испытания, Gk рассматривается как оценка матрицы вторых производных Γ f(y k), а Hk — как оценка ( Γ f )–1. В зависимости от того, какая именно матричная оценка используется, точка очередного измерения будет выбираться по правилу

Методы вида (7.40), (7.41) выбирают направления перемещения, совпадающие с антиградиентными направлениями функции f, вычисленными в некотором измененном пространстве. Фактически, на каждом шаге выполняется изменение метрики пространства переменных (поворот осей и перемасштабирование) за счет домножения градиента на соответствующую матрицу. Поэтому второе название этих методов – методы переменной метрики.

При определенных условиях через N шагов матрица GN будет близка к матрице

Γ f(y N), а HN — к ( Γ f (y N))–1

Рассмотрим основные идеи, лежащие в основе построения оценок Gk и Hk. Рассмотрим случай, когда f(y) = (yТΓy)/2+(c,y)+b — квадратичная функция с симметричной положительно определенной матрицей Γ.

Пусть

Потребуем, чтобы таким же условиям удовлетворяли, соответственно, оценка Gk+1 матрицы Γ и оценка Hk матрицы Γ –1, построенные по (k+1)–му измерениям градиента. А именно, пусть выполняется требование

Условия (7.43) называются квазиньютоновскими условиями.

ННГУ. Учебный курс «Модели и методы конечномерной оптимизации». Часть 2. Нелинейное программирование и многоэкстремальная оптимизация.

Наложим дополнительные требования на матрицы оценок. Поскольку сама матрица Γ симметрична, потребуем выполнения свойства симметрии от матрицы Gk +1, ( Hk+1), положив

Будем определять ее новое значение путем коррекции предыдущей матрицы

Gk+1 = Gk +Uk (или Hk+1 = Hk +Uk).

где поправки U k строятся в виде матриц ранга 1 и находятся из условий (7.43), (7.44).

Эти условия определяют поправку неединственным образом. Рассмотрим возможные приемы ее построения на примере оценок Gk. Оценки вида Hk строятся аналогично.

Простейший способ определения поправочной матрицы, предложенный Бройденом, состоит в том, чтобы составить ее из вектор–столбцов невязок вида z k–Gk ∆ k, помноженных на специально подобранные числа. Нетрудно проверить, что при произвольных v k, для которых (v k)T∆ k ≠0, нужная оценка определяется следующей формулой

Если положить в ней v k = (z k–Gk ∆ k), то получим формулу Бройдена (B–формулу)

Существуют и другие способы оценивания. Например, можно несимметричную поправку в формуле (7.45) заменить похожей симметричной. Непосредственной проверкой можно убедиться (полезно в качестве упражнения выполнить эту проверку), что для любого вектора v k такого, что (v k)Т∆ k≠0, соотношение

дает оценочную матрицу, удовлетворяющую (7.43),(7.44).

Положив в (7.45) v k = z k, получим формулу Девидона–Флетчера–Пауэлла

(DFP–формулу), представимую в следующем виде

Gk +1 = Gk – Gk∆ k(∆ k)ТGk /( (∆k)TGk∆k) + z k(z k)Т/( (z k)T∆k) + (∆ k)TGk∆ kw k(w k)Т, (7.48)

где

Поскольку в (7.48),(7.49) w k(w k)Т – симметричная матрица и, как можно показать, (w k)Т∆ k = 0 (проверьте это в качестве упражнения), то, в силу (7.43), последнее слагаемое в (7.48) можно отбросить. В результате будет получена

сокращенная формула Девидона–Флетчера–Пауэлла. Ее называют также формулой Бройдена–Флетчера–Гольдфарба–Шанно (BFGHформула).

В приведенных итерационных соотношениях начальное значение выбирается в виде G0=E (единичная матрица).

Аналогичные формулы для оценивания обратной матрицы Гессе имеют