В градиентном спуске направления убывания - анти-градиенты, но для функций с плохо обусловленным гессианом сходимость медленная.
Идея: двигаться вдоль направлений, которые гарантируют сходимость за $n$ шагов.
Определение. Множество ненулевых векторов $\{p_0, \ldots, p_l\}$ называется сопряжённым относительно матрицы $A \in \mathbb{S}^n_{++}$, если
$$ p^{\top}_iAp_j = 0, \qquad i \neq j $$Утверждение. Для любой $x_0 \in \mathbb{R}^n$ последовательность $\{x_k\}$, генерируемая методом сопряжённых направлений, сходится к решению системы $Ax = b$ максимум за $n$ шагов.
def ConjugateDirections(x0, A, b, p):
x = x0
r = A.dot(x) - b
for i in range(len(p)):
alpha = - (r.dot(p[i])) / (p[i].dot(A.dot(p[i])))
x = x + alpha * p[i]
r = A.dot(x) - b
return x
Идея: новое направление $p_k$ ищется в виде $p_k = -r_k + \beta_k p_{k-1}$, где $\beta_k$ выбирается, исходя из требования сопряжённости $p_k$ и $p_{k-1}$:
$$ \beta_k = \dfrac{p^{\top}_{k-1}Ar_k}{p^{\top}_{k-1}Ap^{\top}_{k-1}} $$Таким образом, для получения следующего сопряжённого направления $p_k$ необходимо хранить только сопряжённое направление $p_{k-1}$ и остаток $r_k$ с предыдущей итерации.
Вопрос: как находить размер шага $\alpha_k$?
Теорема Пусть после $k$ итераций $x_k \neq x^*$. Тогда
Теорема 1. Если матрица $A$ имеет только $r$ различных собственных значений, то метод сопряжённых градиентов cойдётся за $r$ итераций.
Теорема 2. Имеет место следующая оценка сходимости
$$ \| x_{k} - x^* \|_A \leq 2\left( \dfrac{\sqrt{\kappa(A)} - 1}{\sqrt{\kappa(A)} + 1} \right)^k \|x_0 - x^*\|_A, $$где $\|x\|_A = x^{\top}Ax$ и $\kappa(A) = \frac{\lambda_1(A)}{\lambda_n(A)}$ - число обусловленности матрицы $A$, $\lambda_1(A) \geq ... \geq \lambda_n(A)$ - собственные значения матрицы $A$
Замечание: сравните коэффициент геометрической прогрессии с аналогом в градиентном спуске.
Упражнение Проверьте численно, насколько быстро растёт обусловленность матрицы из векторов $\{b, Ab, ... \}$
Доказательство
def ConjugateGradientQuadratic(x0, A, b, eps):
r = A.dot(x0) - b
p = -r
while np.linalg.norm(r) > eps:
alpha = r.dot(r) / p.dot(A.dot(p))
x = x + alpha * p
r_next = r + alpha * A.dot(p)
beta = r_next.dot(r_next) / r.dot(r)
p = -r_next + beta * p
r = r_next
return x
Идея: использовать градиенты $f'(x_k)$ неквадратичной функции вместо остатков $r_k$ и линейный поиск шага $\alpha_k$ вместо аналитического вычисления. Получим метод Флетчера-Ривса.
def ConjugateGradientFR(f, gradf, x0, eps):
x = x0
grad = gradf(x)
p = -grad
while np.linalg.norm(gradf(x)) > eps:
alpha = StepSearch(x, f, gradf, **kwargs)
x = x + alpha * p
grad_next = gradf(x)
beta = grad_next.dot(grad_next) / grad.dot(grad)
p = -grad_next + beta * p
grad = grad_next
if restart_condition:
p = -gradf(x)
return x
In [4]:
import numpy as np
n = 100
# Random
# A = np.random.randn(n, n)
# A = A.T.dot(A)
# Clustered eigenvalues
A = np.diagflat([np.ones(n//4), 10 * np.ones(n//4), 100*np.ones(n//4), 1000* np.ones(n//4)])
U = np.random.rand(n, n)
Q, _ = np.linalg.qr(U)
A = Q.dot(A).dot(Q.T)
A = (A + A.T) * 0.5
print("A is normal matrix: ||AA* - A*A|| =", np.linalg.norm(A.dot(A.T) - A.T.dot(A)))
b = np.random.randn(n)
# Hilbert matrix
# A = np.array([[1.0 / (i+j - 1) for i in range(1, n+1)] for j in range(1, n+1)])
# b = np.ones(n)
f = lambda x: 0.5 * x.dot(A.dot(x)) - b.dot(x)
grad_f = lambda x: A.dot(x) - b
x0 = np.zeros(n)
In [5]:
USE_COLAB = False
%matplotlib inline
import matplotlib.pyplot as plt
if not USE_COLAB:
plt.rc("text", usetex=True)
plt.rc("font", family='serif')
if USE_COLAB:
!pip install git+https://github.com/amkatrutsa/liboptpy
import seaborn as sns
sns.set_context("talk")
eigs = np.linalg.eigvalsh(A)
plt.semilogy(np.unique(eigs))
plt.ylabel("Eigenvalues", fontsize=20)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
In [6]:
import scipy.optimize as scopt
def callback(x, array):
array.append(x)
In [7]:
scopt_cg_array = []
scopt_cg_callback = lambda x: callback(x, scopt_cg_array)
x = scopt.minimize(f, x0, method="CG", jac=grad_f, callback=scopt_cg_callback)
x = x.x
print("||f'(x*)|| =", np.linalg.norm(A.dot(x) - b))
print("f* =", f(x))
In [8]:
def ConjugateGradientQuadratic(x0, A, b, tol=1e-8, callback=None):
x = x0
r = A.dot(x0) - b
p = -r
while np.linalg.norm(r) > tol:
alpha = r.dot(r) / p.dot(A.dot(p))
x = x + alpha * p
if callback is not None:
callback(x)
r_next = r + alpha * A.dot(p)
beta = r_next.dot(r_next) / r.dot(r)
p = -r_next + beta * p
r = r_next
return x
In [11]:
import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss
max_iter = 70
print("\t CG quadratic")
cg_quad = methods.fo.ConjugateGradientQuad(A, b)
x_cg = cg_quad.solve(x0, tol=1e-7, max_iter=max_iter, disp=True)
print("\t Gradient Descent")
gd = methods.fo.GradientDescent(f, grad_f, ss.ExactLineSearch4Quad(A, b))
x_gd = gd.solve(x0, tol=1e-7, max_iter=max_iter, disp=True)
print("Condition number of A =", abs(max(eigs)) / abs(min(eigs)))
In [12]:
plt.figure(figsize=(8,6))
plt.semilogy([np.linalg.norm(grad_f(x)) for x in cg_quad.get_convergence()], label=r"$\|f'(x_k)\|^{CG}_2$", linewidth=2)
plt.semilogy([np.linalg.norm(grad_f(x)) for x in scopt_cg_array[:max_iter]], label=r"$\|f'(x_k)\|^{CG_{PR}}_2$", linewidth=2)
plt.semilogy([np.linalg.norm(grad_f(x)) for x in gd.get_convergence()], label=r"$\|f'(x_k)\|^{G}_2$", linewidth=2)
plt.legend(loc="best", fontsize=20)
plt.xlabel(r"Iteration number, $k$", fontsize=20)
plt.ylabel("Convergence rate", fontsize=20)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
In [13]:
print([np.linalg.norm(grad_f(x)) for x in cg_quad.get_convergence()])
In [14]:
plt.figure(figsize=(8,6))
plt.plot([f(x) for x in cg_quad.get_convergence()], label=r"$f(x^{CG}_k)$", linewidth=2)
plt.plot([f(x) for x in scopt_cg_array], label=r"$f(x^{CG_{PR}}_k)$", linewidth=2)
plt.plot([f(x) for x in gd.get_convergence()], label=r"$f(x^{G}_k)$", linewidth=2)
plt.legend(loc="best", fontsize=20)
plt.xlabel(r"Iteration number, $k$", fontsize=20)
plt.ylabel("Function value", fontsize=20)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
In [15]:
import numpy as np
import sklearn.datasets as skldata
import scipy.special as scspec
n = 300
m = 1000
X, y = skldata.make_classification(n_classes=2, n_features=n, n_samples=m, n_informative=n//3)
C = 1
def f(w):
return np.linalg.norm(w)**2 / 2 + C * np.mean(np.logaddexp(np.zeros(X.shape[0]), -y * X.dot(w)))
def grad_f(w):
denom = scspec.expit(-y * X.dot(w))
return w - C * X.T.dot(y * denom) / X.shape[0]
# f = lambda x: -np.sum(np.log(1 - A.T.dot(x))) - np.sum(np.log(1 - x*x))
# grad_f = lambda x: np.sum(A.dot(np.diagflat(1 / (1 - A.T.dot(x)))), axis=1) + 2 * x / (1 - np.power(x, 2))
x0 = np.zeros(n)
print("Initial function value = {}".format(f(x0)))
print("Initial gradient norm = {}".format(np.linalg.norm(grad_f(x0))))
In [16]:
def ConjugateGradientFR(f, gradf, x0, num_iter=100, tol=1e-8, callback=None, restart=False):
x = x0
grad = gradf(x)
p = -grad
it = 0
while np.linalg.norm(gradf(x)) > tol and it < num_iter:
alpha = utils.backtracking(x, p, method="Wolfe", beta1=0.1, beta2=0.4, rho=0.5, f=f, grad_f=gradf)
if alpha < 1e-18:
break
x = x + alpha * p
if callback is not None:
callback(x)
grad_next = gradf(x)
beta = grad_next.dot(grad_next) / grad.dot(grad)
p = -grad_next + beta * p
grad = grad_next.copy()
it += 1
if restart and it % restart == 0:
grad = gradf(x)
p = -grad
return x
In [17]:
import scipy.optimize as scopt
import liboptpy.restarts as restarts
n_restart = 60
tol = 1e-5
max_iter = 600
scopt_cg_array = []
scopt_cg_callback = lambda x: callback(x, scopt_cg_array)
x = scopt.minimize(f, x0, tol=tol, method="CG", jac=grad_f, callback=scopt_cg_callback, options={"maxiter": max_iter})
x = x.x
print("\t CG by Polak-Rebiere")
print("Norm of garient = {}".format(np.linalg.norm(grad_f(x))))
print("Function value = {}".format(f(x)))
print("\t CG by Fletcher-Reeves")
cg_fr = methods.fo.ConjugateGradientFR(f, grad_f, ss.Backtracking("Wolfe", rho=0.9, beta1=0.1, beta2=0.4, init_alpha=1.))
x = cg_fr.solve(x0, tol=tol, max_iter=max_iter, disp=True)
print("\t CG by Fletcher-Reeves with restart n")
cg_fr_rest = methods.fo.ConjugateGradientFR(f, grad_f, ss.Backtracking("Wolfe", rho=0.9, beta1=0.1, beta2=0.4,
init_alpha=1.), restarts.Restart(n // n_restart))
x = cg_fr_rest.solve(x0, tol=tol, max_iter=max_iter, disp=True)
print("\t Gradient Descent")
gd = methods.fo.GradientDescent(f, grad_f, ss.Backtracking("Wolfe", rho=0.9, beta1=0.1, beta2=0.4, init_alpha=1.))
x = gd.solve(x0, max_iter=max_iter, tol=tol, disp=True)
In [18]:
plt.figure(figsize=(8, 6))
plt.semilogy([np.linalg.norm(grad_f(x)) for x in cg_fr.get_convergence()], label=r"$\|f'(x_k)\|^{CG_{FR}}_2$ no restart", linewidth=2)
plt.semilogy([np.linalg.norm(grad_f(x)) for x in cg_fr_rest.get_convergence()], label=r"$\|f'(x_k)\|^{CG_{FR}}_2$ restart", linewidth=2)
plt.semilogy([np.linalg.norm(grad_f(x)) for x in scopt_cg_array], label=r"$\|f'(x_k)\|^{CG_{PR}}_2$", linewidth=2)
plt.semilogy([np.linalg.norm(grad_f(x)) for x in gd.get_convergence()], label=r"$\|f'(x_k)\|^{G}_2$", linewidth=2)
plt.legend(loc="best", fontsize=16)
plt.xlabel(r"Iteration number, $k$", fontsize=20)
plt.ylabel("Convergence rate", fontsize=20)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
In [19]:
%timeit scopt.minimize(f, x0, method="CG", tol=tol, jac=grad_f, options={"maxiter": max_iter})
%timeit cg_fr.solve(x0, tol=tol, max_iter=max_iter)
%timeit cg_fr_rest.solve(x0, tol=tol, max_iter=max_iter)
%timeit gd.solve(x0, tol=tol, max_iter=max_iter)
Рассмотрим такие методы, что
$$ x_{k+1} = x_0 + \texttt{span}\{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$
Привести пример такой функции и доказать эту теорему Вам надо в домашнем задании.
Пусть
$\alpha = \dfrac{1}{L}$
Тогда
$$ f(x_k) - f^* \leq \dfrac{2L \| x_0 - x^*\|^2_2}{k+4} $$
Пусть
$\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$
Картинка отсюда
Пусть $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 $$где $x(\tau) = x_k + \tau(x^* - x_k) $
где $T_i = \begin{bmatrix} 1 + \beta_k - \alpha_k \lambda_i & -\beta_k \\ 1 & 0 \end{bmatrix}$ и $\simeq$ обозначает равенство спектральных радиусов поскольку матрица перестановки является ортогональной
In [117]:
USE_COLAB = False
if USE_COLAB:
!pip install git+https://github.com/amkatrutsa/liboptpy
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
if not USE_COLAB:
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))
In [122]:
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
In [126]:
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
In [127]:
for m in methods:
_ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
In [128]:
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)
In [20]:
for m in methods:
print(m)
%timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
In [3]:
n = 300
m = 1000
X, y = skldata.make_classification(n_classes=2, n_features=n, n_samples=m, n_informative=n//3)
C = 1
def f(w):
return np.linalg.norm(w)**2 / 2 + C * np.mean(np.logaddexp(np.zeros(X.shape[0]), -y * X.dot(w)))
def grad(w):
denom = scspec.expit(-y * X.dot(w))
return w - C * X.T.dot(y * denom) / X.shape[0]
# f = lambda x: -np.sum(np.log(1 - A.T.dot(x))) - np.sum(np.log(1 - x*x))
# grad_f = lambda x: np.sum(A.dot(np.diagflat(1 / (1 - A.T.dot(x)))), axis=1) + 2 * x / (1 - np.power(x, 2))
x0 = np.zeros(n)
print("Initial function value = {}".format(f(x0)))
print("Initial gradient norm = {}".format(np.linalg.norm(grad_f(x0))))
In [67]:
alpha_test = 5e-3
beta_test = 0.9
methods = {
r"GD, $\alpha_k = {}$".format(alpha_test): fo.GradientDescent(f, grad, ss.ConstantStepSize(alpha_test)),
"GD Armijo": fo.GradientDescent(f, grad,
ss.Backtracking("Armijo", rho=0.7, beta=0.1, init_alpha=1.)),
r"HB, $\beta = {}$".format(beta_test): HeavyBall(f, grad, ss.ConstantStepSize(alpha_test), beta=beta_test),
}
# x0 = np.random.rand(n)
x0 = np.zeros(n)
max_iter = 400
tol = 1e-5
In [68]:
for m in methods:
print(m)
_ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
scopt_cg_array = []
scopt_cg_callback = lambda x: callback(x, scopt_cg_array)
x = scopt.minimize(f, x0, tol=tol, method="CG", jac=grad, callback=scopt_cg_callback, options={"maxiter": max_iter})
x = x.x
In [69]:
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.semilogy([np.linalg.norm(grad(x)) for x in scopt_cg_array], label="CG PR")
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)
In [43]:
for m in methods:
print(m)
%timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
In [130]:
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)
In [131]:
beta_test = 0.9
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 = 2000
tol = 1e-6
In [132]:
for m in methods:
_ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
In [133]:
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)
In [48]:
for m in methods:
print(m)
%timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
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) <= fx - 1 / L gradx.dot(h) + 1 / (2 * L) h.dot(h):
break
else:
L = L * rho
return L
In [30]:
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
In [31]:
for m in methods:
_ = methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)
In [32]:
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)
In [33]:
for m in methods:
print(m)
%timeit methods[m].solve(x0=x0, max_iter=max_iter, tol=tol)