<filename>continuous-optimization_17-18_3course/optimization_07.tex<gh_stars>10-100
\documentclass[a4paper, 12pt]{article}
%\input{header.sty}
%\usepackage{algorithm}
%\usepackage[noend]{algpseudocode}

\begin{document}
%\maketitle
\section{Лекция 7}
На прошлых лекциях мы рассмотрели метод Ньютона и его <<безгессианную>> модификацию HFN (Hessian-Free Newton), которая позволяла не вычислять гессиан, а требовала лишь процедуры умножения гессиана на вектор. HFN, хоть и позволяет экономить память, не является серебряной пулей, поскольку асимптотически сложен по времени. Хотелось бы создать метод, который бы работал быстрее классического метода Ньютона, но сохрнял бы в себе сверхлинейную локальную скорость сходимости. Заметьте, мы не стремимся к квадратичной скорости сходимости, но стремимся к сверхлинейной, поскольку уже её бывает достаточно. Кульминацией лекции станет метод L-BFGS, который, хоть и сходится линейно, обладает линейными затратами по памяти и сходится быстрее метода градиентного спуска.

\subsection{Квазиньютоновские методы}
Начинаем мы так же, как начинали при выводе метода Ньютона. Мы будем минимизировать не саму функцию $f(x_k)$, а её квадратичную модель $m_k(d) = f(x_k) + \nabla f(x_k)^T d + \frac{1}{2} d^T B_k d$ и надеяться, что мы таким образом будем косвенно минимизировать $f(x)$. Как обычно, будем строить последовательность решений $x_{k+1} = x_k + \alpha_k d_k$, где $d_k = B_k^{-1}\nabla f(x_k)$.

Основной идеей квазиньютоновских методов является то, что они стараются не вычислять гессиан $B_k$ каждый раз напрямую, а пытаются эффективно обновлять матрицу $B_k$ с предыдущей итерации. Естественно, мы хотим сохранить некоторые свойства гессиана, которые позволяют методу Ньютона добиться локальной квадратичной сходимости. В качестве такого условия выступает \textit{условие секущей}:
$$\nabla m_{k+1}(-\alpha_kd_k) = \nabla f(x_k) \Leftrightarrow \nabla f(x_{k+1}) - \nabla f(x_k) = B_{k+1}(x_{k+1} - x_k)$$
Чтобы разобраться в условии секущей, запишем формулу Ньютона-Лейбница:
\begin{align}\nabla f(x_{k+1}) - \nabla f(x_k) &= \int_{0}^{1} \nabla^2 f(x_k + \tau (x_{k+1} - x_k)) (x_{k+1} - x_k) d\tau \\  &= \left(\int_{0}^{1} \nabla^2 f(x_k + \tau (x_{k+1} - x_k))  d\tau\right)(x_{k+1} - x_k)\end{align}
То, что стоит в скобках, в действительности имеет смысл <<среднего гессиана>>, который, как видно, удовлетворяет условию секущей. В требовании условия секцщей мы выражаем своё желание, чтобы матрица $B_k$, которая некоторым образом обновляется после каждой итерации, была похожа на этот <<средний гессиан>>.

Введем теперь обозначения, к которым будем часто обращаться:
$$
s_k := x_{k+1} - x_k \qquad y_k := \nabla f(x_{k+1}) - \nabla f(x_k)
$$
В терминах $s_k$ и $y_k$ условие секущей запишется следующим образом: $B_{k+1}s_k = y_k$. Вспомним теперь, что для того, чтобы $d_k$ было направлением спуска, необходимо, чтобы матрица $B_k$ была положительно определена. Умножив обе части уравнения секущей слева на $s_k^T$, получим
$$s_k^T B_{k+1} s_k = s_k^Ty_k > 0$$
в силу положительно определенности $B_{k+1}$. Полученное условие является необходимым для выполнения условия секущей. Для доказательства того, что для какого-то метода спуска на любой итерации будет выполнено $s_k^Ty_k > 0$, пользуются следующими достаточными условиями.

\begin{Statement}
Если выполнено хотя бы одно из двух условий:
\begin{enumerate}
    \item $f$ --- сильно выпуклая функция;
    \item в работе метода на каждой итерации $d_k$ --- направление спуска, а $\alpha_k$ удовлетворяет условию Вульфа,
\end{enumerate}
то на любом шаге $s_k^Ty_k > 0$.
\end{Statement}
\begin{proof}
\\~

\begin{enumerate}
    \item Первое условие доказывается очень просто. Запишем условия $\mu$-сильной выпуклости в точках $x_1$ и $x_2$:
$$
\begin{cases}
f(x_2) \geq f(x_1) + \nabla f(x_1)^T (x_2 - x_1) + \frac{\mu}{2}\|x_1 - x_2\|^2 \\
f(x_1) \geq f(x_2) + \nabla f(x_2)^T (x_1 - x_2) + \frac{\mu}{2}\|x_1 - x_2\|^2
\end{cases}
$$
Сложив два неравенства и приведя подобные, получим:
$$\left(\nabla f(x_1) - \nabla f(x_2)\right)^T(x_2 - x_1) + \mu \|x_2 - x_1\|^2 \leq 0$$
откуда
$$\left(\nabla f(x_1) - \nabla f(x_2)\right)^T(x_1 - x_2) \geq \mu \|x_2 - x_1\|^2 > 0$$ для любых $x_1 \neq x_2$. А это значит, что тем более $s_k^Ty_k > 0$.

\item Для доказательства второго утверждения, запишем условие Вульфа:
$$\nabla f(x_{k+1})^T d_k \geq c_2 \nabla f(x_k)^Td_k$$
где $0 < c_2 < 1$. Вычтем теперь из обеих частей выражение $\nabla f(x_k)^T d_k$ и получим
$$\left(\underbrace{\nabla f(x_{k+1}) - \nabla f(x_k)}_{=y_k}\right)^Td_k \geq (c_2 - 1)\nabla f(x_k)^Td_k$$
Заметим теперь, что $c_2 - 1 < 0$, т.к. $c_2 < 1$, а $\nabla f(x_k)^Td_k < 0$ в силу того, что $d_k$ --- направление спуска по условию. Отсюда заключаем, что $y_k^Td_k > 0$. Если теперь умножить обе части на $\alpha_k > 0$ и вспомнить, что $\alpha_kd_k = x_{k+1} - x_k = s_k$, получим требуемое: $s_k^Ty_k > 0$.

\end{enumerate}
\end{proof}

Осталось теперь понять, как мы будем из матрицы $B_k$ получать матрицу $B_{k+1}$, т.е. необходимо задать процесс обновления. Заметим, что требований положительной определенности матрицы и уравнения секущей недостаточно, чтобы однозначно определеить матрицу $B_k$ (положительная определенность матрицы задает $n$ ограничений по критерию Сильвестра, уравнение секущей --- тоже $n$; но число элементов в матрице составляет $n^2$, и ни о какой однозначности говорить не приходится). Поэтому нам необходимо еще сильнее сузить зону поиска $B_{k+1}$. В квазиньютоновских методах обновление происходит по правилу $$B_{k+1} = B_k + \mathrm{LowRankUpdate}(B_k, s_k, y_k)$$
где $\mathrm{LowRankUpdate}(B_k, s_k, y_k)$ --- низкоранговая матрица (мы рассмотрим варианты с 1- и 2-ранговыми коррекциями). Оказывается, это ограничение уже позволяет определить матрицу $B_{k+1}$ однозначно.

Итак, небольшое резюме:
\begin{itemize}
    \item Квазиньютоновские методы не вычисляют матрицу-гессиан на каждой итерации, а строят некоторую последовательность матриц $B_k$, которые обновляются по правилу $B_{k+1} = B_k + \mathrm{LowRankUpdate}(B_k, s_k, y_k)$.
    \item Чтобы сохранить свойства гессиана, которые позволяют методу Ньютона сходиться локально квадратично, мы накладываем ограничения на матрицу $B_k$, а именно она должна удовлетворять условию секущей: $B_k s_k = y_k$.
    \item Необходимым условием для выполнения уравнения секущей является неравенство $s_k^Ty_k > 0$, убедиться в выполнении которого можно с помощью одного из достаточных условий, сформулированных выше.
    \item То, какую именно низкоранговую коррекцию мы будем проводить, определяет конкретный квазиньютоновский метод.
\end{itemize}

\subsection{Symmetric Rank 1 (SR1)}
Название метода говорит нам, что коррекцию $B_k$ мы будем производить с помощью симметричной одноранговой матрицы, т.е.
$$B_{k+1} = B_k + \sigma_kv_kv_k^T$$
где $v_k \in \mathbb{R}^n$ --- вектор, $\sigma_k \in \mathbb{R}$ --- число, которые мы хотим отыскать.

Запишем уравнение секущей: $B_{k+1} = (B_k + \sigma_kv_kv_k^T)s_k = B_ks_k + \sigma_k v_k (v_k^Ts_k)= y_k$. Заметим, что $v_k^Ts_k \in \mathbb{R}$ --- число. Отсюда получаем, что $$\sigma_k(v_k^Ts_k)v_k = y_k - B_ks_k$$
т.е. $v_k$ пропорционален $y_k - B_ks_k$. Обозначив за $\gamma_k$ коэффициент пропорциональности, запишем $v_k = \gamma_k(y_k - B_ks_k)$. Подставляя это выражение для $v_k$ обратно, мы получаем:
$$\left[\sigma_k\gamma_k^2(y_k - B_ks_k)^Ts_k\right](y_k - B_ks_k) = y_k - B_ks_k$$
откуда $\sigma_k\gamma_k^2(y_k - B_ks_k)^Ts_k = 1$ и
$$\sigma_k\gamma_k^2 = \frac{1}{(y_k - B_ks_k)^Ts_k}$$

Запишем, наконец, итоговую формулу для пересчета $B_{k+1}$:
$$B_{k+1}^{\mathrm{SR1}} = B_k + \sigma_kv_kv_k^T = B_k + \sigma_k\gamma_k^2 (y_k - B_ks_k)(y_k - B_ks_k)^T = B_k + \frac{(y_k - B_ks_k)(y_k - B_ks_k)^T}{(y_k - B_ks_k)^Ts_k}$$

Здесь и далее обозначим $H_k := B_k^{-1}$, т.к. для рассчета направления спуска нам нужна не сама матрица $B_k$, а обратная к ней. Из тождества Вудбери можно вывести формулу для обновления $H_{k+1}$:
$$H_{k+1}^{\mathrm{SR1}} = H_k + \frac{(H_ky_k - s_k)(H_ky_k - s_k)^T}{(H_ky_k - s_k)^Ts_k}$$

Заметим, что мы свели задачу к квадратичному времени. Может показаться, что затраты по времени, как и в методе Ньютона, кубические. Однако раскрыв скобки в выражении для $H_k$, мы увидим, что в действительности нам ни разу не придется перемножать две матрицы. Затраты по памяти, очевидно, квадратичные: нам придется хранить $H_k$ в памяти.

Однако необходимо помнить, что мы также требуем положительной определенности матрицы $B_{k+1}$. Проверим, не портит ли SR1-обновление положительную определенноть. Для этого предположим, что матрица $B_k > 0$, и рассмотрим произвольный вектор $x$:
$$x^TB_{k+1}^{\mathrm{SR1}}x = \underbrace{x^TB_kx}_{>0} + \frac{\overbrace{[(y_k - B_ks_k)^Tx]^2}^{\geq 0}}{\underbrace{y_k^Ts_k}_{>0} - \underbrace{s_k^TB_ks_k}_{>0}}$$
Мы видим, что в знаменателе дроби стоит разность двух положительных чисел, а значит знаменатель может быть отрицательным (следовательно, и $x^TB_{k+1}^{\mathrm{SR1}}x$ может быть меньше нуля). Такие случаи действительно возникают, поэтому на практике, чтобы поддерживать положительную определенность, матрицу обновляют по следующему правилу:
$$
B_{k+1}^{\mathrm{SR1}} = \begin{cases}
B_k, & \text{если } (y_k - B_ks_k)^Ts_k \leq \rho \|y_k - B_ks_k\|\|s_k\| \\
B_k + \frac{(y_k - B_ks_k)(y_k - B_ks_k)^T}{(y_k - B_ks_k)^Ts_k}, & \text{иначе}
\end{cases}
$$
где $0 < \rho < 1$.

Осталось лишь понять, как инициализировать $H_0$. Будем искать <<наилучшую>> матрицу среди диагональных матриц вида $\gamma_0I_n$. Для этого первый шаг сделаем согласно методу градиентного спуска:
$$x_1 = x_0 - \alpha_0\nabla f(x_0) \qquad s_0 = x_1 - x_0 \qquad y_0 = \nabla f(x_1) - \nabla f(x_0)$$
Очевидно, что в терминах обратной матрицы уравнение секущей можно переписать так: $B_ks_k = y_k \Leftrightarrow H_ky_k = s_k$. Для нахождения $H_0$ будем минимизировать невязку:
$$\|H_0y_0 - s_0\|^2 = \|\gamma_0y_0-s_0\|^2 \to \min_{\gamma_0}$$
Решением этой задачи, как нетрудно видеть, является
$$\gamma_0 = \frac{s_0^Ty_0}{y_0^Ty_0} \Rightarrow H_0 = \frac{s_0^Ty_0}{y_0^Ty_0}I_n$$

\subsection{BFGS}
В методе BFGS (первые буквы фамилий ученых: Broyden--Fletcher--Goldfarb--Shanno) коррекция производится матрицей ранга 2. Мы сначала сформулируем идею обновления матрицы, а потом покажем, что ранг корректировочной матрицы не больше 2.

Теперь давайте сразу требовать, чтобы $B_{k+1}$ была положительно определенной и симметричной. Каждая положительно определенная симметричная матрица может быть представлена в виде произведения $J_{k+1}J_{k+1}^T$ (такое разложение всегда существует, например, фактор Холецкого). Поэтому будем искать $B_{k+1}$ среди матриц вида $J_{k+1}J_{k+1}^T$, где $J_{k+1}$ --- некоторая невырожденная матрица. Запишем уравнение секущей:
$$B_{k+1}s_k = y_k \Leftrightarrow J_{k+1}\underbrace{J_{k+1}^Ts_k}_{u_k} = y_k$$
Напомним, что в добавок к тому, чтобы выполнялось условие секущей, мы хотим, чтобы $J_{k+1} - J_k$ было низкоранговой матрицей. Давайте выразим это желание в виде стремления $\|J_{k+1} - J_k\|^2$ к минимуму:
$$
\begin{cases}
    \|J_{k+1} - J_k\|_F^2 \to \min_{J_{k+1}} \\
    J_{k+1}u_k = y_k
\end{cases}
$$
Мы надеемся, что минимизировав $\|J_{k+1} - J_k\|^2_F$, мы получим разность $J_{k+1} - J_k$ низкого ранга.

Решив эту задачу с помощью методов множителей Лагранжа, мы приходим к решению
$$J_{k+1} = J_k - \frac{(J_ku_k - y_k)u_k^T}{u_k^Tu_k}$$
Учитывая то, что $u_k = J_{k+1}u_$

\end{document}
