human
I am working on the file "Spring2021/acc_grad.ipynb".
The first task is:
Резюме

Сопряжённые направления
Метод сопряжённых градиентов
Сходимость
Эксперименты

Метод тяжёлого шарика (heavy-ball method)

Предложен в 1964 г. Б.Т. Поляком
<img src="./polyak.jpeg">
Для квадратичной целевой функции зигзагообразное поведение градиентного спуска обусловлено неоднородностью направлений
Давайте учитывать предыдущие направления для поиска новой точки
Метод тяжёлого шарика

$$
x_{k+1} = x_k - \alpha_k f'(x_k) + {\color{red}{\beta_k(x_k - x_{k-1})}}
$$

Помимо параметра шага вдоль антиградиента $\alpha_k$ появился ещё один параметр $\beta_k$

Геометрическая интерпретация метода тяжёлого шарика
Картинка отсюда 
<img src="./heavy_ball.png" width=600 align="center">
Теорема сходимости
Пусть $f$ сильно выпукла с Липшицевым градиентом. Тогда для 
$$
\alpha_k = \frac{4}{(\sqrt{L} + \sqrt{\mu})^2}
$$ 
и 
$$
\beta_k = \max(|1 - \sqrt{\alpha_k L}|^2, |1 - \sqrt{\alpha_k \mu}|^2)
$$ 
справедлива следующая оценка сходимости
$$ 
\left\| \begin{bmatrix} x_{k+1} - x^ \ x_k - x^ \end{bmatrix} \right\|_2
\leq 
\left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k \left \|
\begin{bmatrix} x_1 - x^ \ x_0 - x^ \end{bmatrix}
\right \|_2
$$

Совпадает с оценкой снизу для методов первого порядка!
Оптимальные параметры $\alpha_k$ и $\beta_k$ определяются через неизвестные константы $L$ и $\mu$

Схема доказательства

Перепишем метод как 

\begin{align}
\begin{bmatrix} x_{k+1}\ x_k \end{bmatrix} = \begin{bmatrix} (1 + \beta_k)I & -\beta_k I\ I & 0 \end{bmatrix} \begin{bmatrix} x_k\ x_{k-1} \end{bmatrix} + \begin{bmatrix} -\alpha_k f'(x_k)\ 0 \end{bmatrix}
\end{align}

Используем теорему из анализа

\begin{align}
\begin{bmatrix} x_{k+1} - x^\ x_k - x^ \end{bmatrix} = \underbrace{ \begin{bmatrix} (1 + \beta_k)I - \alpha_k \int_0^1 f''(x(\tau))d\tau & -\beta_k I\ I & 0 \end{bmatrix}}_{=A_k}\begin{bmatrix} x_k - x^\ x_{k-1} - x^\end{bmatrix},
\end{align}
где $x(\tau) = x_k + \tau(x^* - x_k) $
- В силу интегральной теоремы о среднем $A_k(x) = \int_0^1 f''(x(\tau))d\tau = f''(z)$, поэтому $L$ и $\mu$ ограничивают...
Can you write Python program for it?
 ... Can you provide a the Python program you requested?
assistant
#Here is the Python program you requested.

import liboptpy.base_optimizer as base
import numpy as np
import liboptpy.unconstr_solvers.fo as fo
import liboptpy.step_size as ss
import matplotlib.pyplot as plt
%matplotlib inline

plt.rc("text", usetex=True)

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))

np.random.seed(42)
n = 100
A = np.random.randn(n, n)
A = A.T.dot(A)
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)
print("Condition number = {}".format(L / mu))
alpha_opt = 4 / (np.sqrt(L) + np.sqrt(mu))**2 
beta_opt = np.maximum((1 - np.sqrt(alpha_opt * L))**2, 
                      (1 - np.sqrt(alpha_opt * mu))**2)
print(alpha_opt, beta_opt)
beta_test = 0.95

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),
    "HB optimal": HeavyBall(f, grad, ss.ConstantStepSize(alpha_opt), beta = beta_opt),
    "CG": fo.ConjugateGradientQuad(A, b)
}
x0 = np.random.randn(n)
max_iter = 5000
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...