import numpy as np

import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss

import jax
import jax.numpy as jnp
from jax.config import config
config.update("jax_enable_x64", True)
import sklearn.datasets as skldata
n = 1000
m = 2000


X, y = skldata.make_classification(n_classes=2, n_features=n, n_samples=m, n_informative=n//3, random_state=42)
C = 1

"""
Explanation: Метод Ньютона
На прошлых семинарах...

Метод сопряжённых градиентов
Метод тяжёлого шарика
Ускоренный метод Нестерова

Идея метода Ньютона
Рассмотрим задачу
$$
\min\limits_{x\ \in \mathbb{R}^n} f(x).
$$

Градиентный спуск $\equiv$ линейная аппроксимация $f$
Метод Ньютона $\equiv$ квадратичная аппроксимация $f$:

$$
f(x + h) \approx f(x) + \langle f'(x), h \rangle + \frac{1}{2}h^{\top}f''(x)h \to \min_{h}
$$
Из необходимого условия минимума:
$$
f'(x) + f''(x) h = 0, \qquad h^* = -(f''(x))^{-1} f'(x)
$$
Является ли найденное направление направлением убывания?
Проверим знак скалярного произведения $\langle f'(x), h^* \rangle$.
$$
\langle f'(x), h^* \rangle = -(f')^{\top}(x) (f''(x))^{-1} f'(x) < 0 \Leftarrow f''(x) \succ 0
$$
Вопрос: а что если при некотором $k^*$ гессиан станет неопределён?
Метод Ньютона

Классический метод Ньютона: $\alpha_k \equiv 1$
Демпфированный метод Ньютона: $\alpha_k$ выбирается на каждой итерации по заданному правилу

```python
def NewtonMethod(f, x0, epsilon, **kwargs):
x = x0

while True:

    h = ComputeNewtonStep(x, f, **kwargs)

    if StopCriterion(x, f, h, **kwargs) &lt; epsilon:

        break

    alpha = SelectStepSize(x, h, f, **kwargs)

    x = x + alpha * h

return x

```
Теорема сходимости (Ю. Е. Нестеров Введение в выпуклую оптимизацию, $\S$ 1.2)
Теорема. Пусть функция $f(x)$
- дважды дифференцируема и её гессиан удовлетворяет условию Липшица с константой $M$
- существует точка локального минимума с положительно определённым гессианом
$$
f''(x^*) \succeq l\mathbf{I}, \; l > 0
$$

начальная точка $x_0$ достаточно близка к точке минимума, в частности

$$
\|x_0 - x^*\|_2 \leq \frac{2l}{3M}
$$
Тогда метод Ньютона сходится квадратично:
$$
\|x_{k+1} - x^ \|_2 \leq \dfrac{M\|x_k - x^\|^2_2}{2 (l - M\|x_k - x^*\|_2)}
$$
Упражнение
Примените метод Ньютона для поиска корня следующей функции
$$
\varphi(t) = \dfrac{t}{\sqrt{1+t^2}}
$$
и определите область сходимости.
Аффинная инвариантность
Рассмотрим функцию $f(x)$ и невырожденное преобразование с матрицей $A$. 
Выясним, как изменится шаг метода Ньютона после преобразования $A$.
Пусть $x = Ay$ и $g(y) = f(Ay)$. Тогда
$$
g(y + u) \approx g(y) + \langle g'(y), u \rangle + \frac{1}{2} u^{\top} g''(y) u \to \min_{u}
$$
и
$$
u^* = -(g''(y))^{-1} g'(y) \qquad y_{k+1} = y_k - (g''(y_k))^{-1} g'(y_k)
$$
или
\begin{align}
y_{k+1} & = y_k - (A^{\top}f''(Ay_k)A)^{-1} A^{\top}f'(Ay_k)\
& = y_k - A^{-1}(f''(Ay_k))^{-1}f'(Ay_k)
\end{align}
Таким образом,
$$
Ay_{k+1} = Ay_k - (f''(Ay_k))^{-1}f'(Ay_k) \quad x_{k+1} = x_k - (f''(x_k))^{-1}f'(x_k)
$$
Следовательно, направление метода Ньютона преобразуется при 
линейном преобразовани так же, как и координаты!
Использование метода сопряжённых градиентов в методе Ньютона

В методе Ньютона надо решать линейную систему $H(x_k) h_k = -f'(x_k)$ для поиска направлению движения
Если функция сильно выпуклая, то $H(x_k) \in \mathbb{S}^n_{++}$ и для решения системы можно использовать метод сопряжённых градиентов. В этом случае метод называют неточный метод Ньютона.
Что нового это даёт
Явно хранить гессиан не нужно, достаточно умножать его на вектор
Можно регулировать точность решения системы и не решать её очень точно вдали от решения. Важно: неточное решение может не быть направлением убывания! 
Сходимость будет только локально сверхлинейной при начале выбора шага с $\alpha_0 = 1$, как и в методе Ньютона



Вычислительная сложность и эксперименты
Узкие места метода Ньютона:

формирование и хранение гессиана
решение систем линейных уравнений

$$
f''(x_k)h = -f'(x_k)
$$
Сравнение с градиентным спуском
$$
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
"""


import cvxpy as cvx
x = cvx.Variable(n)

obj = cvx.Minimize(cvx.norm(x, 2)**2 / 2 + C * cvx.sum(cvx.logistic(cvx.multiply(-y, X * x))) / m)
prob = cvx.Problem(obj)
prob.solve(solver="SCS", verbose=True, max_iters=1000)
print("Optimal value =", prob.value)

"""
Explanation: Точное решение с помощью CVXPy
End of explanation
"""


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

grad_f = jax.jit(jax.grad(f))
hess_f = jax.jit(jax.hessian(f))
x0 = jnp.ones(n)
%timeit grad_f(x0).block_until_ready()
%timeit hess_f(x0).block_until_ready()

"""
Explanation: Вспомогательные функции
End of explanation
"""


def Newton(f, gradf, hessf, x0, epsilon, num_iter, line_search, 
                    disp=False, callback=None, **kwargs):
    x = x0.copy()
    iteration = 0
    opt_arg = {"f": f, "grad_f": gradf}
    for key in kwargs:
        opt_arg[key] = kwargs[key]
    while True:
        gradient = gradf(x)
        hess = hessf(x)
        h = np.linalg.solve(hess, -gradient)
        alpha = line_search(x, h, **opt_arg)
        x = x + alpha * h
        if callback is not None:
            callback(x)
        iteration += 1
        if disp:
            print("Current function val =", f(x))
            print("Current gradient norm = ", np.linalg.norm(gradf(x)))
        if np.linalg.norm(gradf(x)) < epsilon:
            break
        if iteration >= num_iter:
            break
    res = {"x": x, "num_iter": iteration, "tol": np.linalg.norm(gradf(x))}
    return res

"""
Explanation: Реализация метода Ньютона
End of explanation
"""


newton = methods.so.NewtonMethod(f, grad_f, hess_f,
                                 ss.Backtracking("Armijo", rho=0.1, beta=0.1, init_alpha=1.),
#                                  ss.ConstantStepSize(1.),
                                linsolver=jnp.linalg.solve)
# newton = methods.so.NewtonMethod(f, grad_f, hess_f, ss.ConstantStepSize(5e-1))
x_newton = newton.solve(x0, tol=1e-8, max_iter=100, disp=True)

gd = methods.fo.GradientDescent(f, grad_f, ss.Backtracking("Armijo", rho=0.1, beta=0.1, init_alpha=1.))
# agd = methods.fo.AcceleratedGD(f, grad_f, ss.ConstantStepSize(1e-2))
# gd = methods.fo.GradientDescent(f, grad_f, ss.ConstantStepSize(1e-5))
x_gd = gd.solve(x0, tol=1e-3, max_iter=100, disp=True)
# x_agd = agd.solve(x0, tol=1e-3, max_iter=3000, disp=True)

%matplotlib inline
import matplotlib.pyplot as plt
plt.rc("text", usetex=True)
    
plt.figure(figsize=(12, 8))
# Newton
plt.semilogy([np.linalg.norm(grad_f(x)) for x in newton.get_convergence()], label="$\| f'(x_k) \|^{N}_2$")
# Gradient
plt.semilogy([np.linalg.norm(grad_f(x)) for x in gd.get_convergence()], label="$\| f'(x_k) \|^{G}_2$")
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel(r"Convergence rate", fontsize=26)
plt.xticks(fontsize = 24)
plt.yticks(fontsize = 24)
plt.grid(True)
plt.legend(loc="best", fontsize=24)

"""
Explanation: Сравнение с градиентным спуском
End of explanation
"""


%timeit newton.solve(x0, tol=1e-3, max_iter=5000).block_until_ready()
%timeit agd.solve(x0, tol=1e-3, max_iter=5000).block_until_ready()

"""
Explanation: Сравнение времени работы
End of explanation
"""


import numpy as np
import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss
%matplotlib inline
import matplotlib.pyplot as plt
import scipy.optimize as scopt
plt.rc("text", usetex=True)

n = 3000
m = 100
x0 = np.zeros(n)
max_iter = 100
tol = 1e-5
A = np.random.rand(m, n) * 10

f = lambda x: -np.sum(np.log(1 - A.dot(x))) - np.sum(np.log(1 - x*x))
grad_f = lambda x: np.sum(A.T / (1 - A.dot(x)), axis=1) + 2 * x / (1 - np.power(x, 2))

def bb_method(f, gradf, x0, tol=1e-6, maxiter=100, callback=None, alpha_type=1):
    it = 0
    x_prev = x0.copy()
    current_tol = np.linalg.norm(gradf(x_prev))
    alpha = 1e-4
    while current_tol > tol and it < maxiter:
        it += 1
        current_grad = gradf(x_prev)
        if it != 1:
            g = current_grad - prev_grad
            if alpha_type == 1:
                alpha = g.dot(s) / g.dot(g)
            elif alpha_type == 2:
                alpha = s.dot(s) / g.dot(s)
        if callback:
            callback(x_prev)
        x_next = x_prev - alpha * current_grad
        current_tol = np.linalg.norm(gradf(x_next))
        prev_grad = current_grad
        s = x_next - x_prev
        x_prev = x_next
    if callback:
        callback(x_prev)
    return x_next

method = {
    "BB 1": methods.fo.BarzilaiBorweinMethod(f, grad_f, init_alpha=1e-4, type=1),
    "BFGS": methods.fo.BFGS(f, grad_f),
    "DFP": methods.fo.DFP(f, grad_f),
    "LBFGS": methods.fo.LBFGS(f, grad_f),
}

for m in method:
    print("\t Method {}".format(m))
    _ = method[m].solve(x0=x0, tol=tol, max_iter=max_iter, disp=True)

print("\t Method BFGS Scipy")
scopt_conv = []
scopt_res = scopt.minimize(f, x0, method="BFGS", jac=grad_f, callback=lambda x: scopt_conv.append(x), 
                           tol=tol, options={"maxiter": max_iter})
print("Result: {}".format(scopt_res.message))
if scopt_res.success:
    print("Convergence in {} iterations".format(scopt_res.nit))
print("Function value = {}".format(f(scopt_res.x)))

plt.figure(figsize=(8, 6))

for m in method:
    plt.semilogy([np.linalg.norm(grad_f(x)) for x in method[m].get_convergence()], label=m)

plt.semilogy([np.linalg.norm(grad_f(x)) for x in [x0] + scopt_conv], label="BFGS SciPy")
plt.ylabel("$\|f'(x_k)\|_2$", fontsize=18)
plt.xlabel("Number of iterations, $k$", fontsize=18)
plt.legend(fontsize=18)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)

for m in method:
    print("\t Method {}".format(m))
    %timeit method[m].solve(x0=x0, tol=tol, max_iter=max_iter)

%timeit scopt.minimize(f, x0, method="BFGS", jac=grad_f, tol=tol, options={"maxiter": max_iter})

"""
Explanation: Метод Ньютона даёт очень точное решение за длительное время
Градиентный спуска даёт не очень точное решение, но гораздо быстрее
Часто бывает, что очень точное решение не требуется, поэтому градиентный спуск может быть предпочтительнее

Pro & Contra
Pro

квадратичная сходимость вблизи решения
высокая точность полученного решения
аффинная инвариантность
параметры слабо влияют на скорость сходимости

Contra

необходимо хранить гессиан на каждой итерации: $O(n^2)$ памяти
необходимо решать линейные системы: $O(n^3)$ операций
гессиан может оказаться вырожден
гессиан может не быть положительно определён $\to$ направление $-(f''(x))^{-1}f'(x)$ может не быть направлением убывания

Сравнение с градиентным методом (Б.Т. Поляк Введение в оптимизацию, гл. 3, $\S$ 1 )
Метод | Скорость сходимости | Сложность | Аффинная инвариантность | Требования к $f(x)$
:---: | :---: | :---: | :---: | :---
Градиентный спуск | Глобально линейная | $O(n) + $ определение шага | Нет | Дифференцируема; градиент липшицев
Метод Ньютона | Локально квадратичная | $O(n^3) + $ определение шага | Да | Дважды диференцируема; гессиан липшицев, положительно определён
Что дальше?

Сложность: как избавиться от решения систем линейных уравнений и хранения гессиана?
Сходимость: как совместить локально квадратичную и глобально линейную скорости? Желательно получить глобально квадратичную сходимость!
Требования к $f(x)$ необходимо минимизировать
Квазиньютоновские методы частично решают эти проблемы

Как уменьшить сложность хранения и вычисления?


Сложность вычисления можно уменьшить с помощью

Квазиньютоновские методы, они же методы переменной метрики
Требуется хранение матрицы $n \times n$



Сложность вычисления и хранения можно уменьшить 

квазиньютоновские методы с ограниченной памятью, например L-BFGS (Limited Broyden-Fletcher-Goldfarb-Shanno)
НЕ требуется хранить матрицу
вместо этого хранятся $k \ll n$ векторов из $\mathbb{R}^n$



Единообразный способ получения метода Ньютона и градиентного спуска

градиентный метод получен из аппроксимации первого порядка:

$$
f_G(x) \approx f(y) + \langle f'(y), x - y \rangle + \frac{1}{2}(x-y)^{\top} \frac{1}{\alpha}I(x - y)
$$
причём при $\alpha \in (0, 1/L], f(x) \leq f_G(x)$, то есть $f_G$ - глобальная оценка $f(x)$
- метод Ньютона получен из аппроксимации второго порядка
$$
f_N(x) \approx f(y) + \langle f'(y), x - y \rangle + \frac{1}{2} (x-y)^{\top}f''(y)(x-y)
$$
Идея: использовать промежуточную аппроксимацию вида
$$
f_q(x) \approx f(y) + \langle f'(y), x - y \rangle + \frac{1}{2} (x-y)^{\top}{\color{red}{B(y)}}(x-y),
$$
которая даёт переход к следующей точке:
$$
x_{k+1} = x_k - \alpha_k B^{-1}_k f'(x_k) = x_k - \alpha_k H_k f'(x_k)
$$
Немного истории...

Первый квазиньютоновский метод придумал физик William Davidon в середине 1950-х для ускорения своих вычислений на ненадёжных компьютерах
Его статью с описанием предложенного метода не приняли к публикации, и она оставалась техническим отчётом <br></br> более 30 лет
Опубликована в 1991 году в первом выпуске SIAM Journal on Optimization 

Общая схема квазиньютоновских методов
```python
def QuasiNewtonMethod(f, x0, epsilon, **kwargs):
x = x0

H = I

while True:

    h = -H.dot(grad_f(x))

    if StopCriterion(x, f, h, **kwargs) &lt; epsilon:

        break

    alpha = SelectStepSize(x, h, f, **kwargs)

    x = x + alpha * h

    H = UpdateH(H, f(x), grad_f(x))

return x

```
Как искать $B_{k+1}$?
В точке $x_{k+1}$ имеем следующую аппрокисмацию:
$$
f_q(h) \approx f(x_{k+1}) + \langle f'(x_{k+1}), h \rangle + \frac{1}{2}h^{\top}B_{k+1}h
$$
Из определения, очевидно, что $B_{k+1} \in \mathbb{S}^n_{++}$.
Какие требования естественно наложить на $f_q(h)$?
$$
f_q'(-\alpha_k h_k) = f'(x_k) \qquad f'q(0) = f'(x{k+1}),
$$
где первое условие даёт
$$
f'(x_{k+1}) - \alpha_k B_{k+1}h_k = f'(x_k),
$$
а второе выполняется автоматически.
Квазиньютоновское уравнение (Secant equation)
Из первого условия получаем
$$
B_{k+1}s_k = y_k,
$$
где $s_k = x_{k+1} - x_k$ и $y_k = f'(x_{k+1}) - f'(x_k)$.
Это уравнение будет иметь решение только при $s^{\top}_k y_k > 0$. Почему?
Вопрос: единственным ли образом определено $B_{k+1}$?
Как однозначно определить $B_{k+1}$?
\begin{align}
& \min_B \| B_k - B \| \
\text{s.t. } & B = B^{\top}\
& Bs_k = y_k
\end{align}
DFP (Davidon-Fletcher-Powell)
$$
B_{k+1} = (I - \rho_k y_k s^{\top}_k)B_k(I - \rho_k s_ky^{\top}_k) + \rho_k y_k y^{\top}_k,
$$
где $\rho_k = \dfrac{1}{y^{\top}_k s_k}$,
или с помощью формулы Шермана-Морисона-Вудбери
$$
B^{-1}{k+1} = H{k+1} = H_k - \dfrac{H_ky_k y_k^{\top}H_k}{y^{\top}_kH_ky_k} + \dfrac{s_ks^{\top}_k}{y^{\top}_ks_k}
$$
Вопрос: какой ранг у разности матриц $B_{k+1} (H_{k+1})$ и $B_{k} (H_{k})$?
Вывод
Общая идея квазиньютоновских методов: 
вместо полного пересчёта гессиана на каждой итерации обновлять 
текущую его аппроксимацию с помощью легко вычислимого 
преобразования
BFGS
<img src="./bfgs.png" width=500>
Вопрос: какая естественная модификация метода DFP?
\begin{align}
& \min_H \| H_k - H \| \
\text{s.t. } & H = H^{\top}\
& Hy_k = s_k
\end{align}
Формула пересчёта для метода BFGS:
$$
H_{k+1} = (I - \rho_k s_ky^{\top}_k)H_k(I - \rho_k y_k s^{\top}_k) + \rho_k s_k s^{\top}_k,
$$
где $\rho_k = \dfrac{1}{y^{\top}_k s_k}$
Детали реализации

Не должно быть операций сложностью $O(n^3)$, то есть никаких матричных умножений и решений линейных систем (cf. реализацию в SciPy v. 0.18.1)
Только правило Вольфа гарантирует соблюдения условия кривизны $y_k^{\top}s_k > 0$
Параметры в правиле Вольфа обычно следующие
$\alpha_0 = 1$ необходим для сверхлинейной скорости
$\beta_1 = 10^{-4}$, $\beta_2 = 0.9$


Способы инициализации $H_0$
единичная матрица
$H_0 = \frac{y_0^{\top}s_0}{y_0^{\top}y_0}I$ после первого шага, но до вычисления $H_1$.При вычислении $x_1$  используется $H_0 = I$
$H_0 = \delta \|g_0\|^{-1}_2 I$, параметр $\delta$ необходимо заранее задать



```python
def update_H(x_next, x_current):
current_grad = grad(x_next)

s = x_next - x_current

y = current_grad - grad_mem[-1]

rho = 1. / y.dot(s)

if H is None:

    H = np.eye(x_current.shape[0]) / y.dot(y) / rho

Hy = H.dot(y)

Hys = np.outer(Hy, s)

ss = np.outer(s, s)

H = rho * ss + H - rho * Hys - rho * Hys.T + rho**2 * y.dot(Hy) * ss

x_current = x_next

return H

```
Сходимость
Теорема
Пусть $f$ дважды непрерывно дифференцируема и её гессиан липшицев, также пусть последовательность генерируемая методом BFGS сходится к точке $x^$ так что $\sum_{k=1}^{\infty} \|x_k - x^\| < \infty$. Тогда $x_k \to x^*$ сверхлинейно.
Самокоррекция

Если BFGS на некоторой итерации даёт плохую оценку обратного гессиана, то через несколько итераций это недоразумение будет автоматически исправлено, то есть метод сам скорректирует грубую оценку гессиана
Это свойство появляется только при правильном способе выбора шага, например при использовании правила Вольфа
Метод DFP существенно хуже корректирует неточные оценки обратного гессиана
Всё это будет ниже проиллюстрировано на примерах

BFGS с ограниченной памятью (L-BFGS)

В методе BFGS нужна не сама матрица $H$, а только функция умножения её на вектор 
Поскольку требуется локальная оценка гессиана, старые значения векторов $s$ и $y$ могут портить текущую оценку

Идея

Хранить $k \ll n$ последних векторов $s$ и $y$ - снижение требуемой памяти с $n^2$ до $kn$
Выполнение умножения на вектор рекурсивно, без явного формирования матрицы $H$

<img src="./l-bfgs-cite.png" width=600>
```python
def get_lbfgs_direction(x):
if H is None:

    current_grad = grad(x)

    return -current_grad

else:

    q = current_grad

    alpha = np.zeros(len(s_hist))

    rho = np.zeros(len(s_hist))

    for i in range(len(s_hist) - 1, -1, -1):

        rho[i] = 1. / s_hist[i].dot(y_hist[i])

        alpha[i] = s_hist[i].dot(q) * rho[i]

        q = q - alpha[i] * y_hist[i]

    r = q * H

    for i in range(len(s_hist)):

        beta = rho[i] * y_hist[i].dot(r)

        r = r + s_hist[i] * (alpha[i] - beta)

return -r

```
Barzilai-Borwein method

Первая статья об этом методе опубликована в 1988, в журнале IMA Journal of Numerical Analysis
Статья на NIPS 2016 о модификации этого метода в случае использования стохастической оценки градиента
Идея: комбинация идеи наискорейшего спуска и квазиньютоновского метода

Идея метода

Наискорейший спуск: $x_{k+1} = x_k - \alpha_k f'(x_k)$, $\alpha_k = \arg \min\limits_{\alpha > 0} f(x_{k+1})$
Метод Ньютона: $x_{k+1} = x_k - (f''(x_k))^{-1} f'(x_k)$
Аппроксимация гессиана диагональной матрицей:

$$
\alpha_k f'(x_k) = \alpha_k I f'(x_k) = \left( \frac{1}{\alpha_k} I \right)^{-1} f'(x_k) \approx f''(x_k))^{-1} f'(x_k)
$$

Как найти $\alpha_k$?

Снова квазиньютоновское уравнение (Secant equation)

Для точного гессиана
$$
f''(x_{k})(x_{k} - x_{k-1}) = f'(x_{k}) - f'(x_{k-1})
$$
Для приближения

$$
\alpha_k^{-1} s_{k-1} \approx y_{k-1}
$$

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

Три способа найти $\alpha_k$


Первый способ

Задача

$$
\min_{\beta} \|\beta s_{k-1} - y_{k-1} \|^2_2
$$

Решение

$$
\alpha = \frac{1}{\beta} = \frac{s^{\top}{k-1} s{k-1}}{s^{\top}{k-1} y{k-1}}
$$


Второй способ

Задача

$$
\min_{\alpha} \| s_{k-1} - \alpha y_{k-1} \|^2_2
$$

Решение

$$
\alpha = \frac{s^{\top}{k-1} y{k-1}}{y^{\top}{k-1} y{k-1}}
$$
- Третий способ называется немонотонный линейный поиск: специальная модификация правил Армихо, учитывающая историю изменений значения функции, статья 2004 г. в SIAM Journal on Optimization


Эксперименты
Поиск аналитического центра системы неравенств - пример самокоррекции метода BFGS
$$
f(x) = - \sum_{i=1}^m \log(1 - a_i^{\top}x) - \sum\limits_{i = 1}^n \log (1 - x^2_i) \to \min_x
$$
End of explanation
"""


n = 50
D = np.arange(1, n+1)
U = np.random.randn(n, n)
U, _ = np.linalg.qr(U)
A = U.dot(np.diag(D)).dot(U.T)
b = np.random.randn(n)
eig_vals = np.linalg.eigvals(A)
print("Condition number = {}".format(np.max(eig_vals) / np.min(eig_vals)))

f = lambda x: 0.5 * x.T.dot(A.dot(x)) - b.dot(x)
gradf = lambda x: A.dot(x) - b
x0 = np.random.randn(n)

method = {
    "BB 1": methods.fo.BarzilaiBorweinMethod(f, gradf, init_alpha=1e-4, type=1),
    "BB 2": methods.fo.BarzilaiBorweinMethod(f, gradf, init_alpha=1e-4, type=2),
    "BFGS": methods.fo.BFGS(f, gradf),
    "DFP": methods.fo.DFP(f, gradf),
    "GD": methods.fo.GradientDescent(f, gradf, ss.ExactLineSearch4Quad(A, b)),
    "LBFGS": methods.fo.LBFGS(f, gradf, hist_size=10),
}

for m in method:
    print("\t Method {}".format(m))
    _ = method[m].solve(x0=x0, tol=tol, max_iter=max_iter, disp=True)

print("\t Method BFGS Scipy")

scopt_conv = []
scopt_res = scopt.minimize(f, x0, method="BFGS", jac=gradf, callback=lambda x: scopt_conv.append(x), 
                           tol=tol, options={"maxiter": max_iter})
print("Result: {}".format(scopt_res.message))
if scopt_res.success:
    print("Convergence in {} iterations".format(scopt_res.nit))
print("Function value = {}".format(f(scopt_res.x)))

print("\t Method L-BFGS Scipy")

scopt_lbfgs_conv = []
scopt_res = scopt.minimize(f, x0, method="L-BFGS-B", jac=gradf, tol=tol, 
                           options={"maxiter": max_iter, 'maxcor': 10, "ftol": 1e-10, "gtol": 1e-6},
                           callback=lambda x: scopt_lbfgs_conv.append(x),
                           )
print("Result: {}".format(scopt_res.message))
if scopt_res.success:
    print("Convergence in {} iterations".format(scopt_res.nit))
print("Function value = {}".format(f(scopt_res.x)))

plt.figure(figsize=(12, 8))
fontsize = 26
for m in method:   
    plt.semilogy([np.linalg.norm(gradf(x)) for x in method[m].get_convergence()], label=m)

plt.semilogy([np.linalg.norm(gradf(x)) for x in [x0] + scopt_conv], label='BFGS SciPy')

plt.semilogy([np.linalg.norm(gradf(x)) for x in [x0] + scopt_lbfgs_conv], label='LBFGS SciPy')
plt.legend(fontsize=fontsize)
plt.ylabel("$\|f'(x_k)\|_2$", fontsize=fontsize)
plt.xlabel("Number of iterations, $k$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)

for m in method:
    print("\t Method {}".format(m))
    %timeit method[m].solve(x0=x0, tol=tol, max_iter=max_iter)

%timeit scopt.minimize(f, x0, method="BFGS", jac=gradf, tol=tol, options={"maxiter": max_iter})
%timeit scopt.minimize(f, x0, method="L-BFGS-B", jac=gradf, tol=tol, options={"maxiter": max_iter, 'maxcor': 10, "ftol": 1e-10, "gtol": 1e-6})

"""
Explanation: Плохо обусловленная квадратичная задача
End of explanation
"""


