Метод коллинеарных градиентов

Метод коллинеарных градиентов (МКГ) — итерационный метод направленного поиска локального экстремума гладкой функции многих переменных $J(u)\colon \mathbb {R} ^{n}\to \mathbb {R}$ с движением к экстремуму вдоль вектора $d\in \mathbb {R} ^{n}$ такого, где градиенты $\nabla J(u)$ и $\nabla J(u+d)$ коллинеарные. Это метод перового порядка (использует только первые производные $\nabla J$) с квадратичной скоростью сходимости. Может применяться к функциям высокой размерности $n$ с несколькими локальными экстремумами. МКГ можно отнести к семейству методов Truncated Newton method

Коллинеарные векторы $\nabla J(u^{k})$ и $\nabla J(u^{k_{\ast }})$ с направлением минимизации $d^{k}$ для выпуклой функции, $n=2$

Идея метода

Для гладкой функции $J(u)$ в относительно большой окрестности точки $u^{k}$ найдётся точка $u^{k_{\ast }}$, где градиенты $\nabla J^{k}\,{\overset {\textrm {def}}{=}}\,\nabla J(u^{k})$ и $\nabla J^{k_{\ast }}\,{\overset {\textrm {def}}{=}}\,\nabla J(u^{k_{\ast }})$ коллинеарные. Направлением на экстремум $u_{\ast }$ из точки $u^{k}$ будет направление $d^{k}={(u^{k_{\ast }}-u}^{k})$. Вектор $d^{k}$ указывает на максимум или на минимум в зависимости от положения точки $u^{k_{\ast }}$. Она может быть спереди или сзади от $u^{k}$ относительно направления на $u_{\ast }$ (см. рисунок). Далее будем рассматривать минимизацию.

Очередная итерация МКГ:

(1) $\quad u^{k+1}=u^{k}+b^{k}d^{k},\quad k\in \left\{0,1\dots \right\},$

где оптимальное $b^{k}\in \mathbb {R}$ находится аналитически из предположения квадратичности одномерной функции $J(u^{k}+bd^{k})$:

(2) $\quad b^{k}=\left(1-{\frac {\langle \nabla J(u^{k_{\ast }},d^{k}\rangle }{\langle \nabla J(u^{k}),d^{k}\rangle }}\right)^{-1},\quad \forall u^{k_{\ast }}.$

Угловые скобки — это скалярное произведение в евклидовом пространстве $\mathbb {R} ^{n}$. Если $J(u)$ выпуклая функция в окрестности $u^{k}$, то для передней точки $u^{k_{\ast }}$ получаем число $b^{k}>0$, для задней $b^{k}<0$. Делаем шаг (1).

Для строго выпуклой квадратичной функции $J(u)$ шаг МКГ

 $b^{k}d^{k}=-H^{-1}\nabla J^{k},$

т.е. это шаг метода Ньютона (метод второго порядка с квадратичной скоростью сходимости), где $H$ — матрица Гессе. Такие шаги обеспечивают МКГ квадратичную скорость сходимости.

В общем случае, если $J(u)$ имеет переменную выпуклость и возможны седловые точки, то следует контролировать направление минимизации по углу $\gamma$ между векторами $\nabla J^{k}$ и $d^{k}$. Если $\cos(\gamma )={\frac {\langle \nabla J^{k},d^{k}\rangle }{||\nabla J(u^{k})||\;||d^{k}||}}\geq 0$, то $d^{k}$ — это направление максимизации и в (1) следует брать $b^{k}$ с обратным знаком.

Поиск коллинеарных градиентов

Коллинеарность градиентов оценивается невязкой их ортов, которая имеет вид системы $n$ уравнений для поиска корня $u=u^{k_{\ast }}$:

(3) $\quad r^{k}(u)={\frac {\nabla J(u)}{||\nabla J(u)||}}s-{\frac {\nabla J(u^{k})}{||\nabla J(u^{k})||}}=0\in \mathbb {R} ^{n},$

где знак $s=\operatorname {sgn} \langle \nabla J(u),\nabla J(u^{k})\rangle$ позволяет одинаково оценивать коллинеарность градиентов по одну или разные стороны от минимума $u_{\ast }$, $||r^{k}(u)||\leq {\sqrt {2}}$.

Система (3) решается итерационно (подитерации $l\,$) методом сопряжённых градиентов в предположении, что она линейна в окрестности $u^{k}$:

(4) $\quad u^{k_{l+1}}=u^{k_{l}}+\tau ^{l}p^{l},\quad l=1,2\ldots ,$

где вектор $p^{l}\,{\overset {\textrm {def}}{=}}\,p(u^{k_{l}})=-r^{l}+{\beta ^{l}p}^{l-1}$, $r^{l}\,{\overset {\textrm {def}}{=}}\,r(u^{k_{l}})$, $\beta ^{l}=||r^{l}||^{2}/||r^{l-1}||^{2},\beta ^{1,n,2n...}=0$, $\tau ^{l}=||r^{l}||^{2}/\langle p^{l},H^{l}p^{l}\rangle$, произведение матрицы Гессе $H^{l}$ на $p^{l}$ находится численным дифференцированием:

(5) $\quad H^{l}p^{l}\approx {\frac {r(u^{k_{h}})-r(u^{k_{l}})}{h/||p^{l}||}},$

где $u^{k_{h}}=u^{k_{l}}+hp^{l}/||p^{l}||$, $h$ — малое положительное число такое, что $\langle p^{l},H^{l}p^{l}\rangle \neq 0$.

Начальное приближение задаётся под 45° ко всем осям координат длинной $\delta ^{k}$:

(6) $\quad u_{i}^{k_{1}}=u_{i}^{k}+{\frac {\delta ^{k}}{\sqrt {n}}}\operatorname {sgn} {\ \nabla _{i}J}^{k},\quad i=1\ldots n.$

Начальный радиус $\delta ^{k}$-окрестности точки $u^{k}$ корректируется:

(7) $\quad \delta ^{k}=\max \left[\min \left(\delta ^{k-1}{\frac {||\nabla J(u^{k})||}{||\nabla J(u^{k-1})||}},\delta ^{0}\right),\delta _{m}\right],\quad k>0.$

Необходимо $||u^{k_{l}}-u^{k}||\geq \delta ^{m},\quad l\geq 1$. Здесь малое положительное число $\delta _{m}$ заметно больше машинного эпсилон.

Подитерации $l$ завершаются при выполнении хотя бы одного из условий:

1. $||r^{l}||\leq c_{1}{\sqrt {2}},\quad 0\leq c_{1}<1$ — достигнута точность;
2. $\left|{\frac {||r^{l}||-||r^{l-1}||}{||r^{l}||}}\right|\leq c_{1},\quad l>1$ — прекратилась сходимость;
3. $l\leq l_{max}=\operatorname {integer} \left|c_{2}\ln c_{1}\ln n\right|,\quad c_{2}\geq 1$ — избыточность подитераций.

Алгоритм выбора направления минимизации

• Параметры: $c_{1},c_{2},\delta ^{0},\delta _{m}=10^{-15}\delta ^{0},h=10^{-5}$.
• Входные данные: $n,k=0,u^{k},J(u^{k}),\nabla J(u^{k}),\nabla J^{k}$.

1. $l=1$. Если $k>0$ задаём $\delta ^{k}$ из (7).
2. Находим $u^{k_{l}}$ из (6).
3. Вычисляем $\nabla J^{l},||\nabla J^{l}||$ и находим $r^{l}$ из (3) при $u=u^{k_{l}}$.
4. Если $||r^{l}||\leq c_{1}{\sqrt {2}}\,$ или $l\geq l_{max}$, или $||u^{k_{l}}-u^{k}||<\delta _{m}$, или {$\,l>1$ и $\left|{\frac {||r^{l}||-||r^{l-1}||}{||r^{l}||}}\right|\leq c_{1}$}, то принимаем $u^{k_{\ast }}=u^{k_{l}}$, возвращаем $\nabla J\left(u^{k_{\ast }}\right)$, $d^{k}={(u^{k_{\ast }}-u}^{k})$, стоп $\color {red}\blacksquare$.
5. Если $l\neq 1,n,2n,3n\ldots$, задаём $\beta ^{l}=||r^{l}||^{2}/||r^{l-1}||^{2}$, иначе $\beta ^{l}=0$.
6. Вычисляем $p^{l}=-r^{l}+\beta ^{l}p^{l-1}$.
7. Находим шаговый множитель $\tau ^{l}$ для подитераций:
   1. запоминаем $u^{k_{l}}$, $\nabla J^{l}$, $||\nabla J^{l}||$, $r^{l}$, $||r^{l}||$;
   2. задаём $u^{k_{h}}=u^{k_{l}}+hp^{l}/||p^{l}||$, вычисляем $\nabla J(u^{k_{h}})$, $r\left(u^{k_{h}}\right)$ и находим $H^{l}p^{l}$ из (5), присваиваем $w\leftarrow \langle p^{l},H^{l}p^{l}\rangle$;
   3. если $w=0$, тогда $h\leftarrow 10h$, возвращаемся к шагу 7.2;
   4. восстанавливаем $u^{k_{l}}$, $\nabla J^{l}$, $||\nabla J^{l}||$, $r^{l}$, $||r^{l}||$;
   5. находим $\tau ^{l}=||r^{l}||^{2}/w$.
8. Делаем подитерацию $u^{k_{l+1}}$ из (4).
9. $l\leftarrow l+1$, переходим к шагу 3.

Параметр $c_{2}=3\div 5$. Для функций без седловых точек рекомендуется $c_{1}\approx 10^{-8}$, $\delta \approx {10}^{-5}$. Для «обхода» седловых точек рекомендуется $c_{1}\approx 0.1$, $\delta \approx 0.1$.

Описанный алгоритм позволяет приблизительно найти коллинеарные градиенты из системы уравнений (3). Полученное направление $b^{k}d^{k}$ для алгоритма МКГ (1) будет приблизительным направлением Ньютона (truncated Newton method).

Демонстрации

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

Тестовая функция «повёрнутый эллипсоид»

Строго выпуклая квадратичная функция:

Минимизация МКГ, $n=2$

$J(u)=\sum _{i=1}^{n}\left(\sum _{j=1}^{i}u_{j}\right)^{2},\quad u_{\ast }=(0...0).$

На рисунке для ${\color {red}n=2}$ заданы три чёрные стартовые точки $u^{0}$. Серые точки — подитерации $u^{0_{l}}$ с $\delta ^{0}=0.5$ (показано пунктиром, завышено для демонстрации). Параметры $c_{1}=10^{-8}$, $c_{2}=4$. Для всех $u^{0}$ потребовалась одна итерация и подитераций $l$ не более двух.

При ${\color {red}n=1000}$ (параметр $\delta ^{0}={10}^{-5}$) с начальной точкой $u^{0}=(-1...1)$ МКГ достиг $u_{\ast }$ с точностью 1 % за 3 итерации и 754 вычисления $J$ и $\nabla J$. Другие методы первого порядка: Квазиньютоновский BFGS (работа с матрицами) потребовал 66 итераций и 788 вычислений; сопряжённых градиентов (Fletcher-Reeves) — 274 итерации и 2236 вычислений; конечно-разностный метод Ньютона — 1 итерация и 1001 вычислений. Метод Ньютона второго порядка — 1 итерация.

С ростом размерности $\color {red}n$, вычислительные погрешности при реализации условия коллинеарности (3), могут заметно возрастать. Поэтому МКГ, по сравнению с методом Ньютона, в рассматриваемом примере потребовал более одной итерации.

Минимизация МКГ: 3 итерации и 16 вычислений $J$ и $\nabla J$

Тестовая функция Розенброка

$J(u)=100(u_{1}^{2}-u_{2})^{2}+(u_{1}-1)^{2},\quad u_{\ast }=(1,1).$

Параметры те же, кроме $\delta ^{0}={10}^{-5}$. Траектория спуска МКГ полностью совпадает с методом Ньютона. На рисунке синяя начальная точка $u^{0}=\left(-0.8;-1.2\right)$, красная — $u_{\ast }$. В каждой точке нарисованы орты градиентов.

Тестовая функция Химмельблау

$J(u)=(u_{1}^{2}+u_{2}-11)^{2}+(u_{1}+u_{2}^{2}-7)^{2}.$

Параметры $c_{1}=0.1$, $\delta =0.05$.

Минимизация МКГ: 7 итераций и 22 вычисления $J$ и $\nabla J$. Красные линии — $\cos {\gamma }\geq 0$.	Минимизация методом Ньютона: 9 итераций ($b^{k}=1$)
Метод сопряжённых градиентов (Fletcher-Reeves): 9 итерации и 62 вычисления $J$ и $\nabla J$	Квазиньютоновский BFGS: 6 итераций и 55 вычислений $J$ и $\nabla J$. Красная линия (нарушение условия кривизны) — метод наискорейшего спуска.

МКГ является очень экономичным по количеству вычислений $J$ и $\nabla J$. Благодаря формуле (2), он не требует затратных вычислений шагового множителя $b^{k}$ посредством линейного поиска (например, методом золотого сечения и т.п.).