Основано на статье Smooth minimization of non-smooth functions by Y. Nesterov
Функция негладкая, но имеет структуру $f = h + g$, где $h$ - гладкая выпуклая, $g$ - негладкая выпуклая
Проксимальный градиентный метод
$$ \epsilon \sim O\left(\frac{L}{k}\right), \quad k \sim O\left(\frac{L}{\epsilon}\right) $$
Ускоренный проксимальный градиентный метод
$$ \epsilon \sim O\left(\frac{L}{k^2}\right), \quad k \sim O\left(\sqrt{\frac{L}{\epsilon}}\right) $$
ускоренными методами для гладких задач
зависит от $\frac{L_{\mu}}{\epsilon_{\mu}}$, где $L_{\mu}$ - константа Липшица градиента функции $f_{\mu}$.
In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
plt.rc("text", usetex=True)
x = np.linspace(-3, 3, 10000)
f = np.abs(x)
def huber(x_range, mu):
h = np.zeros_like(x_range)
for i, x in enumerate(x_range):
if np.abs(x) <= mu:
h[i] = x**2 / (2 * mu)
else:
h[i] = np.abs(x) - mu / 2
return h
mu = 0.5
huber_f = huber(x, mu)
plt.figure(figsize=(12, 8))
plt.plot(x, f, label="$|x|$")
plt.plot(x, huber_f, label=r"Huber $(\mu = {})$".format(mu))
plt.grid(True)
plt.xticks(fontsize=28)
plt.yticks(fontsize=28)
plt.legend(fontsize=28)
Out[1]:
откуда получаем оценку на сходимость по функции
$$ f(x) - f^* \leq f_{\mu}(x) - f_{\mu}^* + \frac{m\mu}{2} $$аналогично одномерному случаю + правила дифференцирования сложной функции (выведите это выражение!)
тогда
$$ O\left(\sqrt{\frac{L_{\mu}}{\epsilon_{\mu}}}\right) = O\left(\frac{1}{\epsilon}\right) $$для ускоренного градиентного метода
In [2]:
import liboptpy.unconstr_solvers as methods
import liboptpy.step_size as ss
m, n = 500, 100
A = np.random.randn(m, n)
x_true = np.random.randn(n)
b = A.dot(x_true)
In [3]:
f = lambda x: np.linalg.norm(A.dot(x) - b, 1)
subgrad = lambda x: A.T.dot(np.sign(A.dot(x) - b))
alpha = 1e-3
s = 1e-1
sg_methods = {
"SM 1/k": methods.fo.SubgradientMethod(f, subgrad, ss.InvIterStepSize()),
"SM fixed={}".format(alpha): methods.fo.SubgradientMethod(f, subgrad, ss.ConstantStepSize(alpha)),
"SM scaled fix, s={}".format(s): methods.fo.SubgradientMethod(f, subgrad,
ss.ScaledConstantStepSize(s)),
}
In [4]:
x0 = np.random.randn(n)
max_iter = 50000
In [5]:
for m in sg_methods:
_ = sg_methods[m].solve(x0=x0, max_iter=max_iter)
In [6]:
plt.figure(figsize=(10, 8))
for m in sg_methods:
plt.semilogy([f(x) for x in sg_methods[m].get_convergence()], label=m)
plt.legend(fontsize=20)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel(r"Objective, $f(x_k)$", fontsize=26)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
In [11]:
def huber(x, mu):
res = np.zeros_like(x)
res[np.abs(x) <= mu] = x[np.abs(x) <= mu]**2 / (2 * mu)
res[np.abs(x) >= mu] = np.abs(x[np.abs(x) >= mu]) - mu / 2
return res
def grad_huber(x, mu):
res = np.zeros_like(x)
res[np.abs(x) <= mu] = x[np.abs(x) <= mu] / mu
res[np.abs(x) >= mu] = np.sign(x[np.abs(x) >= mu])
return res
In [7]:
mu = 1e-2
print(mu)
fmu = lambda x: np.sum(huber(A.dot(x) - b, mu))
grad_fmu = lambda x: A.T @ grad_huber(A.dot(x) - b, mu)
In [12]:
grad_methods = {
"GD Armijo": methods.fo.GradientDescent(fmu, grad_fmu, ss.Backtracking("Armijo", beta=0.1, rho=0.5, init_alpha=1)),
"GD": methods.fo.GradientDescent(fmu, grad_fmu, ss.ConstantStepSize(mu/np.linalg.norm(A, 2)**2)),
"Fast GD": methods.fo.AcceleratedGD(fmu, grad_fmu, ss.ConstantStepSize(mu/np.linalg.norm(A, 2)**2))
}
In [13]:
max_iter = 10000
for m in grad_methods:
_ = grad_methods[m].solve(x0=x0, max_iter=max_iter)
In [14]:
plt.figure(figsize=(10, 8))
for m in grad_methods:
plt.semilogy([fmu(x) for x in grad_methods[m].get_convergence()], label=m)
plt.legend(fontsize=20)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel(r"Objective, $f_{\mu}(x_k)$", fontsize=26)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
In [15]:
plt.figure(figsize=(10, 8))
for m in grad_methods:
plt.semilogy([np.linalg.norm(grad_fmu(x)) for x in grad_methods[m].get_convergence()], label=m)
plt.legend(fontsize=20)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel(r"$\|f_{\mu}'(x_k)\|_2$", fontsize=26)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
In [16]:
plt.figure(figsize=(10, 8))
for m in grad_methods:
plt.semilogy([f(x) for x in grad_methods[m].get_convergence()], label=m)
for m in sg_methods:
plt.semilogy([f(x) for x in sg_methods[m].get_convergence()], label=m)
plt.legend(fontsize=20)
plt.xlabel(r"Number of iterations, $k$", fontsize=26)
plt.ylabel(r"Objective, $f(x_k)$", fontsize=26)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
In [17]:
for m in grad_methods:
x_m = grad_methods[m].get_convergence()[-1]
print(m, "\n\t", np.linalg.norm(x_true - x_m), f(x_m) - f(x_true))
for m in sg_methods:
x_m = sg_methods[m].get_convergence()[-1]
print(m, "\n\t", np.linalg.norm(x_true - x_m), f(x_m) - f(x_true))
выпукла всегда
Теорема. Если $f$ - сильно выпуклая функция с константой $m$, сопряжённая к которой
$$ f^*(y) = \sup_{x} \; (y^{\top}x - f(x)) $$тогда
1. Так как $f(x)$ - сильно выпуклая функция, то у функции $y^{\top} x - f(x)$ единственная точка максимума для каждого $y$. Обозначим её $x_y$
2. В силу условия оптимальности выполнено
$$ y = f'(x_y) \quad f^*(y) = \langle y, x_y \rangle - f(x_y) $$3. Тогда для произвольного $u$ выполнено
$$ f^*(u) = \sup_v \; (u^{\top}v - f(v)) \geq u^{\top}x_y - f(x_y) = x_y^{\top}(u - y) + x_y^{\top}y - f(x_y) = x_y^{\top}(u - y) + f^*(y) $$4. То есть по определению субдифференциала $x_y \in \partial f^*(y)$, но так как $x_y$ единственная точка максимума функции $y^{\top} x - f(x)$, то $\partial f^*(y) = \{x_y\}$ и функция $f^*$ дифференцируема
5. В итоге $\nabla f^*(y) = x_y$
или
$$ \| \nabla f^* (v) - \nabla f^*(u) \|_2 \leq \frac{1}{m}\| u - v\|_2 $$Определение. Функция $d$ называется функцией близости на выпуклом замкнутом множестве $C$ если
Центром множества будем называть точку
Дополнительные предположения:
Тогда выполнено
где $h$ - выпукла, замкнута и имеет ограниченную область определения
- Оценка сверху
$$
f_{\mu}(x) \leq \sup_{y \in \mathrm{dom}(h)} \; ((Ax + b)^{\top}y - h(y)) = f(x)
$$
- Оценка снизу
\begin{align*}
f_{\mu}(x) & = \sup_{y \in \mathrm{dom}(h)} \; ((Ax + b)^{\top}y - h(y) - \mu d(y)) \\
& \geq \sup_{y \in \mathrm{dom}(h)} \; ((Ax + b)^{\top}y - h(y) - \mu D) = f(x) - \mu D
\end{align*}где $\square$ обозначает инфимальную конволюцию двух функций
Таким образом, получаем следующую интерпретацию проксимального метода
$$ x_{k+1} = prox_{\mu g}(x_k) = x_k - \mu g'_{\mu}(x_k) $$как градиентного спуска для сглаженной аппроксимации исходной функции.
Теорема. $x^*$ точка минимума функции $g$ iff $x^*$ точка минимума для $g_{\mu}$
Доказательство
1. $x^*$ точка минимума функции $g$ iff $x^* = prox_{g}(x^*)$
2. Если $x^*$ точка минимума $g$, тогда $x^* = prox_{g}(x^*) = x^* - \mu g'_{\mu}(x^*)$, то есть $g'_{\mu}(x^*) = 0$, следовательно в силу выпуклости $g_{\mu}(x)$ точка $x^*$ и для неё является точкой минимума.
3. Если $x^*$ точка минимума для функции $g_{\mu}(x)$ тогда $g'_{\mu}(x^*) = 0$, но тогда $prox_{\mu g}(x^*) = x^*$, то есть $x^*$ неподвижная точка проксимального оператора и тем самым показано, что $x^*$ точка минимума $g$
для $h$ такой что $f(x) = h^*(Ax + b)$
то как правильнее подходить к её решению?