Также на странице курса.
где $S \subseteq \mathbb{R}^n$, $f_j: S \rightarrow \mathbb{R}, \; j = 0,\ldots,m$, $g_k: S \rightarrow \mathbb{R}, \; k=1,\ldots,p$
Все функции как минимум непрерывны.
Важный факт</span>: задачи нелинейной оптимизации
в их самой общей форме являются численно неразрешимыми!
если $x^*$ точка локального минимума дифференцируемой функции $f(x)$, тогда
$$ f'(x^*) = 0 $$если $x^*$ точка локального минимума дважды дифференцируемой функции $f(x)$, тогда
$$ f'(x^*) = 0 \quad \text{и} \quad f''(x^*) \succeq 0 $$пусть $f(x)$ дважды дифференцируемая функция, и пусть точка $x^*$ удовлетворяет условиям
$$ f'(x^*) = 0 \quad f''(x^*) \succ 0, $$тогда $x^*$ является точкой строго локального минимума функции $f(x)$.
Замечание: убедитесь, что Вы понимаете, как доказывать эти
результаты!
Но ведь $x^*$ неизвестна!
Тогда
\begin{align*} & \|x_{k+1} - x_k \| = \|x_{k+1} - x_k + x^* - x^* \| \leq \\ & \|x_{k+1} - x^* \| + \| x_k - x^* \| \leq 2\varepsilon \end{align*}Аналогично для сходимости по функции,
однако иногда можно оценить $f^*$!
Замечание: лучше использовать относительные изменения
этих величин!
Например $\dfrac{\|x_{k+1} - x_k \|_2}{\| x_k \|_2}$
Определение: оракулом называют некоторое абстрактное
устройство, которое отвечает на последовательные вопросы
метода
Аналогия из ООП:
Концепция чёрного ящика
Метод доверительных областей: фиксируется допустимый размер области по некоторой норме $\| \cdot \| \leq \alpha$ и модель целевой функции, которая хорошо её аппроксимирует в выбранной области.
Далее производится поиск направления $h_k$, минимизирующего модель целевой функции и не выводящего точку $x_k + h_k$ за пределы доверительной области
Для заданного класса задач сравнивают следующие величины:
1. Сублинейная
$$ \| x_{k+1} - x^* \|_2 \leq C k^{\alpha}, $$где $\alpha < 0$ и $ 0 < C < \infty$
2. Линейная (геометрическая прогрессия)
$$ \| x_{k+1} - x^* \|_2 \leq Cq^k, $$где $q \in (0, 1)$ и $ 0 < C < \infty$
3. Сверхлинейная
$$ \| x_{k+1} - x^* \|_2 \leq Cq^{k^p}, $$где $q \in (0, 1)$, $ 0 < C < \infty$ и $p > 1$
4. Квадратичная
$$ \| x_{k+1} - x^* \|_2 \leq C\| x_k - x^* \|^2_2, \qquad \text{или} \qquad \| x_{k+1} - x^* \|_2 \leq C q^{2^k} $$где $q \in (0, 1)$ и $ 0 < C < \infty$
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
USE_COLAB = False
if not USE_COLAB:
plt.rc("text", usetex=True)
import numpy as np
C = 10
alpha = -0.5
q = 0.9
num_iter = 50
sublinear = np.array([C * k**alpha for k in range(1, num_iter + 1)])
linear = np.array([C * q**k for k in range(1, num_iter + 1)])
superlinear = np.array([C * q**(k**2) for k in range(1, num_iter + 1)])
quadratic = np.array([C * q**(2**k) for k in range(1, num_iter + 1)])
plt.figure(figsize=(12,8))
plt.semilogy(np.arange(1, num_iter+1), sublinear,
label=r"Sublinear, $\alpha = -0.5$", linewidth=5)
# plt.semilogy(np.arange(1, num_iter+1), superlinear, linewidth=5,
# label=r"Superlinear, $q = 0.5, p=2$")
plt.semilogy(np.arange(1, num_iter+1), linear,
label=r"Linear, $q = 0.5$", linewidth=5)
# plt.semilogy(np.arange(1, num_iter+1), quadratic,
# label=r"Quadratic, $q = 0.5$", linewidth=5)
plt.xlabel("Number of iterations, $k$", fontsize=28)
plt.ylabel("Error rate upper bound", fontsize=28)
plt.legend(loc="best", fontsize=26)
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)
Что дают теоремы сходимости
Что НЕ дают теоремы сходимости
Мораль: нужно проявлять разумную осторожность
и здравый смысл!
Методы нулевого порядка: оракул возвращает только значение функции $f(x)$
Методы первого порядка: оракул возвращает значение функции $f(x)$ и её градиент $f'(x)$
Методы второго порядка: оракул возвращает значение функции $f(x)$, её градиент $f'(x)$ и гессиан $f''(x)$.
Вопрос: существуют ли методы более высокого порядка?
Идея из информатики первого семестра:
делим отрезок $[a,b]$ на две равные части
пока не найдём минимум унимодальной функции.
Тогда
$$ |x_{K+1} - x^*| \leq \frac{b_{K+1} - a_{K+1}}{2} = \left( \frac{1}{2} \right)^{\frac{N-1}{2}} (b - a) \approx 0.5^{K} (b - a) $$
In [55]:
def binary_search(f, a, b, epsilon, callback=None):
c = (a + b) / 2.0
while abs(b - a) > epsilon:
# Check left subsegment
y = (a + c) / 2.0
if f(y) <= f(c):
b = c
c = y
else:
# Check right subsegment
z = (b + c) / 2.0
if f(c) <= f(z):
a = y
b = z
else:
a = c
c = z
if callback is not None:
callback(a, b)
return c
In [56]:
def my_callback(a, b, left_bound, right_bound, approximation):
left_bound.append(a)
right_bound.append(b)
approximation.append((a + b) / 2.0)
In [57]:
import numpy as np
left_boud_bs = []
right_bound_bs = []
approximation_bs = []
callback_bs = lambda a, b: my_callback(a, b,
left_boud_bs, right_bound_bs, approximation_bs)
# Target unimodal function on given segment
f = lambda x: (x - 2) * x * (x + 2)**2 # np.power(x+2, 2)
# f = lambda x: -np.sin(x)
x_true = -2
# x_true = np.pi / 2.0
a = -3
b = -1.5
epsilon = 1e-8
x_opt = binary_search(f, a, b, epsilon, callback_bs)
print(np.abs(x_opt - x_true))
plt.figure(figsize=(10,6))
plt.plot(np.linspace(a,b), f(np.linspace(a,b)))
plt.title("Objective function", fontsize=28)
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)
Идея:
делить отрезок $[a,b]$ не на две равные насти,
а в пропорции "золотого сечения".
Оценим скорость сходимости аналогично методу дихотомии:
$$ |x_{K+1} - x^*| \leq b_{K+1} - a_{K+1} = \left( \frac{1}{\tau} \right)^{N-1} (b - a) \approx 0.618^K(b-a), $$где $\tau = \frac{\sqrt{5} + 1}{2}$.
In [20]:
def golden_search(f, a, b, tol=1e-5, callback=None):
tau = (np.sqrt(5) + 1) / 2.0
y = a + (b - a) / tau**2
z = a + (b - a) / tau
while b - a > tol:
if f(y) <= f(z):
b = z
z = y
y = a + (b - a) / tau**2
else:
a = y
y = z
z = a + (b - a) / tau
if callback is not None:
callback(a, b)
return (a + b) / 2.0
In [58]:
left_boud_gs = []
right_bound_gs = []
approximation_gs = []
cb_gs = lambda a, b: my_callback(a, b, left_boud_gs, right_bound_gs, approximation_gs)
x_gs = golden_search(f, a, b, epsilon, cb_gs)
print(f(x_opt))
print(f(x_gs))
print(np.abs(x_opt - x_true))
In [22]:
plt.figure(figsize=(10,6))
plt.semilogy(np.arange(1, len(approximation_bs) + 1), np.abs(x_true - np.array(approximation_bs, dtype=np.float64)), label="Binary search")
plt.semilogy(np.arange(1, len(approximation_gs) + 1), np.abs(x_true - np.array(approximation_gs, dtype=np.float64)), label="Golden search")
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel("Error rate upper bound", fontsize=26)
plt.legend(loc="best", fontsize=26)
plt.xticks(fontsize = 26)
_ = plt.yticks(fontsize = 26)
In [59]:
%timeit binary_search(f, a, b, epsilon)
%timeit golden_search(f, a, b, epsilon)
In [60]:
f = lambda x: np.sin(np.sin(np.sin(np.sqrt(x))))
x_true = (3 * np.pi / 2)**2
a = 2
b = 60
epsilon = 1e-8
plt.plot(np.linspace(a,b), f(np.linspace(a,b)))
plt.xticks(fontsize = 28)
_ = plt.yticks(fontsize = 28)
In [61]:
left_boud_bs = []
right_bound_bs = []
approximation_bs = []
callback_bs = lambda a, b: my_callback(a, b,
left_boud_bs, right_bound_bs, approximation_bs)
x_opt = binary_search(f, a, b, epsilon, callback_bs)
print(np.abs(x_opt - x_true))
In [62]:
left_boud_gs = []
right_bound_gs = []
approximation_gs = []
cb_gs = lambda a, b: my_callback(a, b, left_boud_gs, right_bound_gs, approximation_gs)
x_gs = golden_search(f, a, b, epsilon, cb_gs)
print(np.abs(x_opt - x_true))
In [27]:
plt.figure(figsize=(8,6))
plt.semilogy(np.abs(x_true - np.array(approximation_bs, dtype=np.float64)), label="Binary")
plt.semilogy(np.abs(x_true - np.array(approximation_gs, dtype=np.float64)), label="Golden")
plt.legend(fontsize=28)
plt.xticks(fontsize=28)
_ = plt.yticks(fontsize=28)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel("Error rate upper bound", fontsize=26)
Out[27]:
In [64]:
%timeit binary_search(f, a, b, epsilon)
%timeit golden_search(f, a, b, epsilon)
def DescentMethod(f, x0, epsilon, **kwargs):
x = x0
while StopCriterion(x, f, **kwargs) > epsilon:
h = ComputeDescentDirection(x, f, **kwargs)
alpha = SelectStepSize(x, h, f, **kwargs)
x = x + alpha * h
return x
Рассмотрим линейную аппроксимацию дифференцируемой функции $f$ вдоль некоторого направления убывания $h, \|h\|_2 = 1$:
$$ f(x + \alpha h) = f(x) + \alpha \langle f'(x), h \rangle + o(\alpha) $$Из условия убывания
$$ f(x) + \alpha \langle f'(x), h \rangle + o(\alpha) < f(x) $$и переходя к пределу при $\alpha \rightarrow 0$:
$$ \langle f'(x), h \rangle \leq 0 $$Также из неравенства Коши-Буняковского-Шварца
$$ \langle f'(x), h \rangle \geq -\| f'(x) \|_2 \| h \|_2 = -\| f'(x) \|_2 $$Таким образом, направление антиградиента
$$ h = -\dfrac{f'(x)}{\|f'(x)\|_2} $$даёт направление наискорейшего локального убывания функции$~f$.
В итоге метод имеет вид
$$ x_{k+1} = x_k - \alpha f'(x_k) $$Рассмотрим обыкновенное диференциальное уравнение вида:
$$ \frac{dx}{dt} = -f'(x(t)) $$и дискретизуем его на равномерной сетке с шагом $\alpha$:
$$ \frac{x_{k+1} - x_k}{\alpha} = -f'(x_k), $$где $x_k \equiv x(t_k)$ и $\alpha = t_{k+1} - t_k$ - шаг сетки.
Отсюда получаем выражение для $x_{k+1}$
$$ x_{k+1} = x_k - \alpha f'(x_k), $$которое в точности совпадает с выражением для градиентного спуска.
Такая схема называется явной или прямой схемой Эйлера.
Вопрос: какая схема называется неявной или обратной?
Глобальная оценка сверху на функцию $f$ в точке $x_k$:
$$ f(y) \leq f(x_k) + \langle f'(x_k), y - x_k \rangle + \frac{L}{2} \|y - x_k \|_2^2 = g(y), $$где $\lambda_{\max}(f''(x)) \leq L$ для всех допустимых $x$.
Справа — квадратичная форма, точка минимума которой имеет аналитическое выражение:
\begin{align*} & g'(y^*) = 0 \\ & f'(x_k) + L (y^* - x_k) = 0 \\ & y^* = x_k - \frac{1}{L}f'(x_k) = x_{k+1} \end{align*}Этот способ позволяет оценить значение шага как $\frac{1}{L}$. Однако часто константа $L$ неизвестна.
def GradientDescentMethod(f, x0, epsilon, **kwargs):
x = x0
while StopCriterion(x, f, **kwargs) > epsilon:
h = ComputeGradient(x, f, **kwargs)
alpha = SelectStepSize(x, h, f, **kwargs)
x = x - alpha * h
return x
Список подходов:
Требование достаточного убывания, требование существенного убывания и условие кривизны: для некоторых $\beta_1, \beta_2$, таких что $0 < \beta_1 < \beta_2 < 1$ найти $x_{k+1}$ такую что
Обычно коэффициенты выбирают так: $\beta_1 \in (0, 0.3)$, а $\beta_2 \in (0.9, 1)$
Требование существенного убывания и условие кривизны обеспечивают убывание функции по выбранному направлению $h_k$. Обычно выбирают одно из них.
In [1]:
%matplotlib notebook
import matplotlib.pyplot as plt
plt.rc("text", usetex=True)
import ipywidgets as ipywidg
import numpy as np
import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss
from tqdm import tqdm
In [2]:
f = lambda x: np.power(x, 2)
gradf = lambda x: 2 * x
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
def update(x0, step):
gd = methods.fo.GradientDescent(f, gradf, ss.ConstantStepSize(step))
_ = gd.solve(np.array([x0]), max_iter=10)
x_hist = gd.get_convergence()
x = np.linspace(-5, 5)
ax.clear()
ax.plot(x, f(x), color="r", label="$f(x) = x^2$")
y_hist = np.array([f(x) for x in x_hist])
x_hist = np.array(x_hist)
plt.quiver(x_hist[:-1], y_hist[:-1], x_hist[1:]-x_hist[:-1], y_hist[1:]-y_hist[:-1],
scale_units='xy', angles='xy', scale=1, width=0.005, color="green", label="Descent path")
ax.legend()
fig.canvas.draw()
step_slider = ipywidg.FloatSlider(value=0.8, min=0, max=1.2, step=0.1, description="Step")
x0_slider = ipywidg.FloatSlider(value=1.5, min=-4, max=4, step=0.1, description="Initial point")
_ = ipywidg.interact(update, x0=x0_slider, step=step_slider)
In [3]:
def plot_alpha(f, grad, x, h, alphas, beta1, beta2):
df = np.zeros_like(alphas)
for i, alpha in enumerate(alphas):
df[i] = f(x + alpha * h)
upper_bound = f(x) + beta1 * alphas * grad(x) * h
lower_bound = f(x) + beta2 * alphas * grad(x) * h
plt.plot(alphas, df, label=r"$f(x + \alpha h)$")
plt.plot(alphas, upper_bound, label="Upper bound")
plt.plot(alphas, lower_bound, label="Lower bound")
plt.xlabel(r"$\alpha$", fontsize=18)
plt.legend(loc="best", fontsize=18)
In [6]:
f = lambda x: x**2
grad = lambda x: 2 * x
beta1 = 0.1
beta2 = 0.9
x0 = 0.5
plot_alpha(f, grad, x0, -grad(x0), np.linspace(1e-3, 1.01, 10), beta1, beta2)
In [7]:
x_range = np.linspace(1e-10, 4)
plt.plot(x_range, x_range * np.log(x_range))
Out[7]:
In [8]:
x0 = 1
f = lambda x: x * np.log(x)
grad = lambda x: np.log(x) + 1
beta1 = 0.3
beta2 = 0.7
plot_alpha(f, grad, x0, -grad(x0), np.linspace(1e-3, 0.9, 10), beta1, beta2)
def SelectStepSize(x, f, h, rho, alpha0, beta1, beta2):
# 0 < rho < 1
# alpha0 - initial guess of step size
# beta1 and beta2 - constants from conditions
alpha = alpha0
# Check violating sufficient decrease and curvature conditions
while (f(x - alpha * h) >= f(x) + beta1 * alpha grad_f(x_k).dot(h)) and
(grad_f(x - alpha * h).dot(h) <= beta2 * grad_f(x_k).dot(h)):
alpha *= rho
return alpha
Теорема 1. Пусть
Тогда для градиентного метода выполнено:
$$ \lim\limits_{k \to \infty} f'(x_k) = 0, $$а функция монотонно убывает $f(x_{k+1}) < f(x_k)$.
Теорема 2. Пусть
Тогда
$$ f(x_k) - f^* \leq \dfrac{2L \| x_0 - x^*\|^2_2}{k+4} $$Теорема 3. Пусть
Тогда
$$ \| x_k - x^*\|_2 \leq \|x_0 - x^*\|_2 q^k, \qquad q = \max(|1 - \alpha l|, |1 - \alpha L|) < 1 $$и минимальное $q^* = \dfrac{L - \mu}{L + \mu}$ при $\alpha^* = \dfrac{2}{L + \mu}$
Из Теоремы 3 имеем
$$ q^* = \dfrac{L - \mu}{L + \mu} = \dfrac{L/\mu - 1}{L/\mu + 1} = \dfrac{M - 1}{M + 1}, $$где $M$ - оценка числа обусловленности $f''(x)$.
Вопрос: что такое число обусловленности матрицы?
Вопрос: какая геометрия у этого требования?
Мораль: необходимо сделать оценку $M$ как можно ближе к 1!
О том, как это сделать, Вам будет предложено подумать в домашнем задании :)
Поиск $\alpha_k$:
In [5]:
def GradientDescent(f, gradf, 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)
alpha = line_search(x, -gradient, **opt_arg)
x = x - alpha * gradient
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 [25]:
def my_f(x, A):
return 0.5 * x.dot(A.dot(x))
def my_gradf(x, A):
return A.dot(x)
In [26]:
plt.rc("text", usetex=True)
gammas = [0.1, 0.5, 1, 2, 3, 4, 5, 10, 20, 50, 100, 1000, 5000, 10000]
# gammas = [1]
num_iter_converg = []
for g in gammas:
A = np.array([[1, 0],
[0, g]], dtype=np.float64)
f = lambda x: my_f(x, A)
gradf = lambda x: my_gradf(x, A)
# x0 = np.random.rand(A.shape[0])
# x0 = np.sort(x0)
# x0 = x0[::-1]
x0 = np.array([g, 1], dtype=np.float64)
# print x0[1] / x0[0]
gd = methods.fo.GradientDescent(f, gradf, ss.ExactLineSearch4Quad(A))
x = gd.solve(x0, tol=1e-7, max_iter=100)
num_iter_converg.append(len(gd.get_convergence()))
plt.figure(figsize=(8, 6))
plt.loglog(gammas, num_iter_converg)
plt.xticks(fontsize = 20)
plt.yticks(fontsize = 20)
plt.xlabel(r"$\gamma$", fontsize=20)
plt.ylabel(r"Number of iterations with $\varepsilon = 10^{-7}$", fontsize=20)
Out[26]:
Пусть $A \in \mathbb{R}^{m \times n}$. Рассмотрим систему линейных неравенств: $Ax \leq 1$ при условии $|x_i| \leq 1$ для всех $i$.
Определение. Аналитическим центром системы неравенств $Ax \leq 1$ при условии $|x_i| \leq 1$ является решение задачи $$ f(x) = - \sum_{i=1}^m \log(1 - a_i^{\top}x) - \sum_{i = 1}^n \log (1 - x^2_i) \to \min_x $$ $$ f'(x) - ? $$
In [3]:
n = 100
m = 200
x0 = np.zeros(n)
A = np.random.rand(n, m)
In [4]:
import cvxpy as cvx
print(cvx.installed_solvers())
x = cvx.Variable(n)
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)
x = x.value
print("Optimal value =", prob.value)
In [5]:
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))
In [6]:
gd = methods.fo.GradientDescent(f, grad_f, ss.Backtracking("Armijo", rho=0.7, beta=0.1, init_alpha=1.))
x = gd.solve(x0, tol=1e-5, max_iter=100, disp=True)
x_conv = gd.get_convergence()
grad_conv = [np.linalg.norm(grad_f(x)) for x in x_conv]
plt.figure(figsize=(8,6))
plt.semilogy(grad_conv, label=r"$\| f'(x_k) \|_2$")
plt.semilogy([np.linalg.norm(x - np.array(x_k)) for x_k in x_conv], label=r"$\|x_k - x^*\|_2$")
plt.semilogy([np.linalg.norm(prob.value - f(np.array(x_k))) for x_k in x_conv], label=r"$\|f(x_k) - f^*\|_2$")
plt.semilogy([np.linalg.norm(np.array(x_conv[i]) - np.array(x_conv[i+1])) for i in range(len(x_conv) - 1)], label=r"$\|x_k - x_{k+1}\|_2$")
plt.semilogy([np.linalg.norm(f(np.array(x_conv[i])) - f(np.array(x_conv[i+1]))) for i in range(len(x_conv) - 1)], label=r"$\|f(x_k) - f(x_{k+1})\|_2$")
plt.xlabel(r"Number of iteration, $k$", fontsize=20)
plt.ylabel(r"Convergence rate", fontsize=20)
plt.xticks(fontsize = 20)
plt.yticks(fontsize = 20)
plt.legend(loc="best", fontsize=20)
plt.tight_layout()
Pro
Contra