<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Beyond gradient descent, vol. 2
Что такое оптимальные методы и ускоренный метод Нестерова
Схема получения оценок снизу на сложность методов и задач
Фиксируем класс функций $\mathcal{F}$
Фиксируем класс методов оптимизации $\mathcal{M}$
Ищем настолько плохую функцию из класса $\mathcal{F}$, что любой метод из класса $\mathcal{M}$ сходится не лучше, чем некоторая оценка
Такая оценка называется оценкой снизу
Пример
Фиксируем класс методов
Рассмотрим такие методы, что
$$
x_{k+1} = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_k)}
$$
Далее в рамках этого семинара для краткости только такие методы будем называть методами первого порядка
Весь последующий анализ НЕ применим, если 
$$
x_{k+1} = x_0 + G(f'(x_0), \ldots, f'(x_k)),
$$
где $G$ - некоторая нелинейная функция
С такими методами мы познакомимся на одном из ближайших занятий
Фиксируем класс функций
Выпуклые функции с липшицевым градиентом
Теорема. Существует выпуклая функция с Липшицевым градиентом, такая что
$$
f(x_t) - f^ \geq \frac{3L\|x_0 - x^\|_2^2}{32(t+1)^2},
$$
где $x_k = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_{k-1})}$, $1 \leq k \leq t$
Привести пример такой функции и доказать эту теорему Вам надо в домашнем задании.
Сильно выпуклые функции с липшицевым градиентом
Теорема. Существует сильно выпуклая функция с Липшицевым градиентом, такая что
$$
f(x_t) - f^ \geq \frac{\mu}{2}\left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^{2t}\|x_0 - x^\|_2^2,
$$
где $\kappa = \frac{L}{\mu}$ и $x_k = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_{k-1})}$, $1 \leq k \leq t$
Привести пример такой функции и доказать эту теорему Вам надо в домашнем задании.
Оценки сходимости известных методов
Оценки сходимости для градиентного спуска
Step1:  Главное про ускоренный метод Нестерова
Теоретически оптимальная скорость сходимости для выпуклых и сильно выпуклых функций
Необходимо подбирать шаг
Волнообразное поведение может быть подавлено с помощью рестартов
Нелинейные композиции градиентов могут дать более быстрые на практике методы (см самое начало семинара)
Оригинальное доказательство сложно понять на интуитивном уровне, есть гораздо более простые подходы, будут рассказаны позднее
Адаптивный выбор константы $L$
В методе Нестерова шаг у градиента равен $\frac{1}{L}$, где $L$ - константа Липшица градиента
Однако она обычно неизвестна
НО! У нас есть условие на эту константу, а именно
$$
f(y) \leq f(x) + \langle f'(x), y - x \rangle + \frac{L}{2}\|x - y\|_2^2
$$
Аналогия с правилами Армихо и прочими такого типа
Более подробно про этот метод смотрите в книге лектора 
```python
def backtracking_L(f, h, grad, x, L0, rho)
Step2:  Эксперимент на неквадратичной задаче
$$
f(w) = \frac12 \|w\|2^2 + C \frac1m \sum{i=1}^m \log (1 + \exp(- y_i \langle x_i, w \rangle)) \to \min_w
$$
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import liboptpy.unconstr_solvers.fo as fo
import liboptpy.step_size as ss
import liboptpy.base_optimizer as base
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(42)
n = 100
A = np.random.randn(n, n)
A = A.T.dot(A)
A_eigvals = np.linalg.eigvalsh(A)
mu = np.min(A_eigvals)
A = A - mu * np.eye(n)
x_true = np.random.randn(n)
b = A.dot(x_true)
f = lambda x: 0.5 * x.dot(A.dot(x)) - b.dot(x)
grad = lambda x: A.dot(x) - b
A_eigvals = np.linalg.eigvalsh(A)
L = np.max(A_eigvals)
mu = np.min(A_eigvals)
print(L, mu)
class HeavyBall(base.LineSearchOptimizer):
    def __init__(self, f, grad, step_size, beta, **kwargs):
        super().__init__(f, grad, step_size, **kwargs)
        self._beta = beta
    
    def get_direction(self, x):
        self._current_grad = self._grad(x)
        return -self._current_grad
    def _f_update_x_next(self, x, alpha, h):
        if len(self.convergence) < 2:
            return x + alpha * h
        else:
            return x + alpha * h + self._beta * (x - self.convergence[-2])
    
    def get_stepsize(self):
        return self._step_size.get_stepsize(self._grad_mem[-1], self.convergence[-1], len(self.convergence))
beta_test = 0.8
methods = {
    "GD fixed": fo.GradientDescent(f, grad, ss.ConstantStepSize(1 / L)),
    "GD Armijo": fo.GradientDescent(f, grad, 
            ss.Backtracking("Armijo", rho=0.5, beta=0.1, init_alpha=1.)),
    r"HB, $\beta = {}$".format(beta_test): HeavyBall(f, grad, ss.ConstantStepSize(1 / L), beta=beta_test),
    "Nesterov": fo.AcceleratedGD(f, grad, ss.ConstantStepSize(1 / L)),
}
x0 = np.random.randn(n)
max_iter = 4000
tol = 1e-6
for m in methods:
    _ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
figsize = (10, 8)
fontsize = 26
plt.figure(figsize=figsize)
for m in methods:
    plt.semilogy([np.linalg.norm(grad(x)) for x in methods[m].get_convergence()], label=m)
plt.legend(fontsize=fontsize, loc="best")
plt.xlabel("Number of iteration, $k$", fontsize=fontsize)
plt.ylabel(r"$\| f'(x_k)\|_2$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)
for m in methods:
    print(m)
    %timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
Explanation: Beyond gradient descent, vol. 2
Что такое оптимальные методы и ускоренный метод Нестерова
Схема получения оценок снизу на сложность методов и задач
Фиксируем класс функций $\mathcal{F}$
Фиксируем класс методов оптимизации $\mathcal{M}$
Ищем настолько плохую функцию из класса $\mathcal{F}$, что любой метод из класса $\mathcal{M}$ сходится не лучше, чем некоторая оценка
Такая оценка называется оценкой снизу
Пример
Фиксируем класс методов
Рассмотрим такие методы, что
$$
x_{k+1} = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_k)}
$$
Далее в рамках этого семинара для краткости только такие методы будем называть методами первого порядка
Весь последующий анализ НЕ применим, если 
$$
x_{k+1} = x_0 + G(f'(x_0), \ldots, f'(x_k)),
$$
где $G$ - некоторая нелинейная функция
С такими методами мы познакомимся на одном из ближайших занятий
Фиксируем класс функций
Выпуклые функции с липшицевым градиентом
Теорема. Существует выпуклая функция с Липшицевым градиентом, такая что
$$
f(x_t) - f^ \geq \frac{3L\|x_0 - x^\|_2^2}{32(t+1)^2},
$$
где $x_k = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_{k-1})}$, $1 \leq k \leq t$
Привести пример такой функции и доказать эту теорему Вам надо в домашнем задании.
Сильно выпуклые функции с липшицевым градиентом
Теорема. Существует сильно выпуклая функция с Липшицевым градиентом, такая что
$$
f(x_t) - f^ \geq \frac{\mu}{2}\left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^{2t}\|x_0 - x^\|_2^2,
$$
где $\kappa = \frac{L}{\mu}$ и $x_k = x_0 + \texttt{span}{f'(x_0), \ldots, f'(x_{k-1})}$, $1 \leq k \leq t$
Привести пример такой функции и доказать эту теорему Вам надо в домашнем задании.
Оценки сходимости известных методов
Оценки сходимости для градиентного спуска: напоминание
Пусть
$f(x)$ дифференцируема на $\mathbb{R}^n$
$f(x)$ выпукла 
$f'(x)$ удовлетворяет условию Липшица с константой $L$
$\alpha = \dfrac{1}{L}$
Тогда 
$$
f(x_k) - f^ \leq \dfrac{2L \| x_0 - x^\|^2_2}{k+4}
$$
Пусть
$f(x)$ дифференцируема на $\mathbb{R}^n$, 
градиент $f(x)$ удовлетворяет условию Липшица с константой $L$
$f(x)$ является сильно выпуклой с константой $\mu$
$\alpha = \dfrac{2}{\mu + L}$
Тогда для градиентного спуска выполнено:
$$
\| x_k - x^ \|_2 \leq \left( \dfrac{\kappa - 1}{\kappa + 1} \right)^k \|x_0 - x^\|_2
$$ 
$$f(x_k) - f^ \leq \dfrac{L}{2} \left( \dfrac{\kappa - 1}{\kappa + 1} \right)^{2k} \| x_0 - x^\|^2_2,
$$
где $\kappa = \frac{L}{\mu}$
Оценки сходимости для метода сопряжённых градиентов: напоминание
Для сильно выпуклой квадратичной функции
$$
f(x) = \frac{1}{2}x^{\top}Ax - b^{\top}x
$$ 
и метода сопряжённых градиентов справедлива следующая оценка сходимости 
$$
2 (f_k - f^) = \| x_{k} - x^ \|_A \leq 2\left( \dfrac{\sqrt{\kappa(A)} - 1}{\sqrt{\kappa(A)} + 1} \right)^k \|x_0 - x^*\|_A,
$$
где $\kappa(A) = \frac{\lambda_1(A)}{\lambda_n(A)} = \frac{L}{\mu}$ - число обусловленности матрицы $A$, $\lambda_1(A) \geq ... \geq \lambda_n(A) > 0$ - собственные значения матрицы $A$
Can we do better?
Существует ли метод, который сходится в соответствии с нижними оценками для
- произвольной сильно выпуклой функции с липшицевым градиентом (не только квадратичной)?
- произвольной выпуклой функции с липшицевым градиентом?
Ускоренный метод Нестерова
Предложен в 1983 г. Ю.Е. Нестеровым
<img src="./nesterov.jpeg">
Одна из возможных форм записи 
\begin{equation}
\begin{split}
& y_0 = x_0 \
& x_{k+1} = y_k - \alpha_k f'(y_k)\
& y_{k+1} = x_{k+1} + \frac{k}{k + 3} (x_{k+1} - x_k)
\end{split}
\end{equation}
Сравните с методом тяжёлого шарика
Tакже не обязательно монотонен
Для любителей геометрии есть альтернативный метод под названием geometric descent с такой же скоростью сходимости 
Геометрическая интерпретация ускоренного метода Нестерова
<img src="nesterov_plot.png" width=600>
Теорема сходимости
Пусть $f$ выпукла с Липшицевым градиентом, а шаг $\alpha_k = \frac{1}{L}$. Тогда ускоренный метод Нестерова сходится как 
$$
f(x_k) - f^ \leq \frac{2L \|x_0 - x^\|_2^2}{(k+1)^2}
$$ 
Пусть $f$ сильно выпукла с липшицевым градиентом. Тогда ускоренный метод Нестерова при шаге $\alpha_k = \frac{1}{L}$ сходится как
$$
f(x_k) - f^* \leq L\|x_k - x_0\|_2^2 \left(1 - \frac{1}{\sqrt{\kappa}} \right)^k
$$
Тестовая задача
$$
f(x) = \frac{1}{2}x^{\top}Ax - b^{\top}x \to \min_x,
$$
где матрица $A$ положительно полуопределённая, то есть функция НЕ является сильно выпуклой
End of explanation
beta_test = 0.9
methods = {
    "GD Armijo": fo.GradientDescent(f, grad, 
            ss.Backtracking("Armijo", rho=0.5, beta=0.1, init_alpha=1.)),
    r"HB, $\beta = {}$".format(beta_test): HeavyBall(f, grad, ss.ConstantStepSize(1 / L), beta=beta_test),
    "Nesterov": fo.AcceleratedGD(f, grad, ss.ConstantStepSize(1 / L)),
    "Nesterov adaptive": fo.AcceleratedGD(f, grad, ss.Backtracking(rule_type="Lipschitz", rho=0.5, init_alpha=1)),
}
x0 = np.zeros(n)
max_iter = 2000
tol = 1e-6
for m in methods:
    _ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
figsize = (10, 8)
fontsize = 26
plt.figure(figsize=figsize)
for m in methods:
    plt.semilogy([np.linalg.norm(grad(x)) for x in methods[m].get_convergence()], label=m)
plt.legend(fontsize=fontsize, loc="best")
plt.xlabel("Number of iteration, $k$", fontsize=fontsize)
plt.ylabel(r"$\| f'(x_k)\|_2$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)
for m in methods:
    print(m)
    %timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
Explanation: Главное про ускоренный метод Нестерова
Теоретически оптимальная скорость сходимости для выпуклых и сильно выпуклых функций
Необходимо подбирать шаг
Волнообразное поведение может быть подавлено с помощью рестартов
Нелинейные композиции градиентов могут дать более быстрые на практике методы (см самое начало семинара)
Оригинальное доказательство сложно понять на интуитивном уровне, есть гораздо более простые подходы, будут рассказаны позднее
Адаптивный выбор константы $L$
В методе Нестерова шаг у градиента равен $\frac{1}{L}$, где $L$ - константа Липшица градиента
Однако она обычно неизвестна
НО! У нас есть условие на эту константу, а именно
$$
f(y) \leq f(x) + \langle f'(x), y - x \rangle + \frac{L}{2}\|x - y\|_2^2
$$
Аналогия с правилами Армихо и прочими такого типа
Более подробно про этот метод смотрите в книге лектора 
```python
def backtracking_L(f, h, grad, x, L0, rho):
L = L0
fx = f(x)
gradx = grad(x)
while True:
    y = x - 1 / L * h
    if f(y) &lt;= fx - 1 / L * gradx.dot(h) + 1 / (2 * L) * h.dot(h):
        break
    else:
        L = L * rho
return L
```
Эксперименты
End of explanation
import jax
import jax.numpy as jnp
from jax.config import config
config.update("jax_enable_x64", True)
import sklearn.datasets as skldata
n = 300
m = 1000
X, y = skldata.make_classification(n_classes=2, n_features=n, n_samples=m, n_informative=n//3, random_state=42)
C = 1
@jax.jit
def f(w):
    return jnp.linalg.norm(w)**2 / 2 +  C * jnp.mean(jnp.logaddexp(jnp.zeros(X.shape[0]), -y * (X @ w)))
# def grad(w):
#     denom = scspec.expit(-y * X.dot(w))
#     return w - C * X.T.dot(y * denom) / X.shape[0]
autograd_f = jax.jit(jax.grad(f))
x0 = jnp.ones(n)
print("Initial function value = {}".format(f(x0)))
print("Initial gradient norm = {}".format(jnp.linalg.norm(autograd_f(x0))))
beta_test = 0.9
L_trial = (1 + C*100)
methods = {
    "GD Armijo": fo.GradientDescent(f, autograd_f, 
            ss.Backtracking("Armijo", rho=0.5, beta=0.01, init_alpha=1.)),
    r"HB, $\beta = {}$".format(beta_test): HeavyBall(f, autograd_f, ss.ConstantStepSize(1 / L_trial), beta=beta_test),
    "Nesterov": fo.AcceleratedGD(f, autograd_f, ss.ConstantStepSize(1 / L_trial)),
    "Nesterov adaptive": fo.AcceleratedGD(f, autograd_f, ss.Backtracking(rule_type="Lipschitz", rho=0.5, init_alpha=1)),
}
max_iter = 2000
tol = 1e-6
for m in methods:
    _ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
figsize = (10, 8)
fontsize = 26
plt.figure(figsize=figsize)
for m in methods:
    plt.semilogy([np.linalg.norm(autograd_f(x)) for x in methods[m].get_convergence()], label=m)
plt.legend(fontsize=fontsize, loc="best")
plt.xlabel("Number of iteration, $k$", fontsize=fontsize)
plt.ylabel(r"$\| f'(x_k)\|_2$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)
for m in methods:
    print(m)
    %timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
Explanation: Эксперимент на неквадратичной задаче
$$
f(w) = \frac12 \|w\|2^2 + C \frac1m \sum{i=1}^m \log (1 + \exp(- y_i \langle x_i, w \rangle)) \to \min_w
$$
End of explanation
<END_TASK>