| Метод | Скорость сходимости | Сложность | Аффинная инвариантность | Требования к $f(x)$ |
|---|---|---|---|---|
| Градиентный спуск | Глобально линейная | $O(n) + $ определение шага | Нет | Дифференцируема; градиент липшицев |
| Метод Ньютона | Локально квадратичная | $O(n^3) + $ определение шага | Да | Дважды диференцируема; гессиан липшицев, положительно определён |
Сложность вычисления можно уменьшить с помощью
Сложность вычисления и хранения можно уменьшить
причём при $\alpha \in (0, 1/L], f(x) \leq f_G(x)$, то есть $f_G$ - глобальная оценка $f(x)$
Идея: использовать промежуточную аппроксимацию вида
$$ 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) $$где первое условие даёт
$$ f'(x_{k+1}) - \alpha_k B_{k+1}h_k = f'(x_k), $$а второе выполняется автоматически.
Вопрос: всегда ли выполнено такое соотношение
между разностью градиентов и точек?
Hint: вспомините условие Вольфа
Вопрос: единственным ли образом определено $B_{k+1}$?
где $\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:
$$ 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}$
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
Идея
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
или
$$ h_{k+1} = -\left(I - \frac{s_k y_k^{\top}}{y_k^{\top}s_k}\right)f'(x_{k+1}) = -\hat{H}_{k+1} f'(x_{k+1}) $$удовлетворяет всем требованиям к матрице в методе BFGS и совпадает с формулой для обновления $H_k$, если $H_k = I$, то есть $k=1$ в методе LBFGS и $H_0 = I$
Первый способ
Задача
$$ \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}} $$
In [1]:
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)
In [2]:
n = 3000
m = 100
x0 = np.zeros(n)
max_iter = 100
tol = 1e-5
A = np.random.rand(m, n) * 10
In [3]:
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))
In [4]:
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
In [5]:
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),
}
In [6]:
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)))
In [10]:
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)
In [11]:
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})
In [8]:
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)))
In [9]:
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)
In [21]:
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),
}
In [22]:
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)))
In [23]:
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.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)
Pro:
Contra: