Рассмотрим задачу
$$ \min\limits_{x\ \in \mathbb{R}^n} f(x). $$Из необходимого условия минимума:
$$ f'(x) + f''(x) h = 0, \qquad h^* = -(f''(x))^{-1} f'(x) $$def NewtonMethod(f, x0, epsilon, **kwargs):
x = x0
while True:
h = ComputeNewtonStep(x, f, **kwargs)
if StopCriterion(x, f, h, **kwargs) < epsilon:
break
alpha = SelectStepSize(x, h, f, **kwargs)
x = x + alpha * h
return x
Теорема. Пусть функция $f(x)$
Тогда метод Ньютона сходится квадратично:
$$ \|x_{k+1} - x^* \|_2 \leq \dfrac{M\|x_k - x^*\|^2_2}{2 (l - M\|x_k - x^*\|_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) $$или
Таким образом,
$$ 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) $$Следовательно, направление метода Ньютона преобразуется при
линейном преобразовани так же, как и координаты!
Узкие места метода Ньютона:
In [1]:
import numpy as np
USE_COLAB = False
if USE_COLAB:
!pip install git+https://github.com/amkatrutsa/liboptpy
import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss
n = 1000
m = 200
x0 = np.zeros((n,))
A = np.random.rand(n, m) * 10
In [2]:
import cvxpy as cvx
x = cvx.Variable((n, 1))
obj = cvx.Minimize(cvx.sum(-cvx.log(1 - A.T * x)) -
cvx.sum(cvx.log(1 - cvx.square(x))))
prob = cvx.Problem(obj)
prob.solve(solver="SCS", verbose=True, max_iters=1000)
print("Optimal value =", prob.value)
In [3]:
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))
hess_f = lambda x: (A.dot(np.diagflat(1 / (1 - A.T.dot(x))**2))).dot(A.T) + np.diagflat(2 * (1 + x**2) / (1 - x**2)**2)
In [4]:
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
In [5]:
newton = methods.so.NewtonMethod(f, grad_f, hess_f, ss.Backtracking("Armijo", rho=0.9, beta=0.1, init_alpha=1.))
x_newton = newton.solve(x0, tol=1e-6, max_iter=50, disp=True)
gd = methods.fo.GradientDescent(f, grad_f, ss.Backtracking("Armijo", rho=0.9, beta=0.1, init_alpha=1.))
x_gd = gd.solve(x0, tol=1e-6, max_iter=50, disp=True)
In [6]:
%matplotlib inline
import matplotlib.pyplot as plt
if not USE_COLAB:
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.legend(loc="best", fontsize=24)
Out[6]:
In [7]:
%timeit newton.solve(x0, tol=1e-6, max_iter=50)
%timeit gd.solve(x0, tol=1e-6, max_iter=50)
Pro
Contra
| Метод | Скорость сходимости | Сложность | Аффинная инвариантность | Требования к $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), $$а второе выполняется автоматически.
Вопрос: единственным ли образом определено $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
Первый способ
Задача
$$ \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 [7]:
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 [8]:
n = 3000
m = 100
x0 = np.zeros(n)
max_iter = 100
tol = 1e-5
A = np.random.rand(m, n) * 10
In [9]:
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 [10]:
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 [11]:
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 [12]:
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 [13]:
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 [14]:
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 [15]:
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 [16]:
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 [39]:
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 [42]:
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)))
In [43]:
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)
In [45]:
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})
Pro:
Contra: