Следующий шаг: произвольная достаточно гладкая функция на достаточно простом множестве - не обязательно полиэдральном.
Определение. Множество будем называть простым, если проекцию на него можно найти существенно быстрее (чаще всего аналитически) по сравнению с решением исходной задачи минимизации.
где $A^+$ - псевдообратная матрица. Если $A$ полного ранга и столбцы линейно-независимы, тогда $A^+ = (A^{\top}A)^{-1}A^{\top}$.
$x_k(\alpha) = \pi_P (x_k - \alpha f'(x_k))$
Теорема. Пусть $f$ выпуклая дифференцируемая функция и её градиент липшицев на $P$ с константой $L$. Пусть $P$ выпуклое и замкнутое множество и $0 < \alpha < 2 / L$.
Тогда
где $I(x|P) = \begin{cases} 0, & x \in P \\ \infty, & \text{иначе} \end{cases}$ - индикаторная функция множества.
Pro
Contra
Определение. Множество $D$ будем называть простым, если можно найти решение следующей задачи
$$ \min_{x \in D} c^{\top}x $$существенно быстрее (чаще всего аналитически) по сравнению с решением исходной задачи минимизации.
Замечание 1: отличие этого определения от предыдущего в линейности целевой функции (была квадратичная), поэтому простых множеств для этого определения больше.
Замечание 2: иногда на допустимое множество легко найти проекцию, но задача линейного программирования является неограниченной. Например, для множества
$$ D = \{ x \in \mathbb{R}^n \; | \; x_i \geq 0 \}, $$проекция на которое очевидна, решение задачи линейного программирования равно $-\infty$, если есть хотя бы одна отрицательная компонента вектора $c$. Теорема с объяснением будет ниже.
Идея: делать шаг не по градиенту, а по направлению, которое точно не выведет из допустимого множества.
Аналогия с градиентным спуском: линейная аппроксимация на допустимом множестве:
$$ f(x_k + s_k) = f(x_k) + \langle f'(x_k), s_k \rangle \to \min_{{\color{red}{s_k \in D}}} $$Начинать поиск нужно с $\alpha_k = 1$
или
$$ f(x) - f(x^*) \leq -\min_{s \in D} \langle f'(x), s - x\rangle = \max_{s \in D} \langle f'(x), x - s\rangle =\max_{s \in D} -\langle f'(x), s \rangle = - \min_{s \in D} \langle f'(x), s \rangle = g(x) $$Получили аналог зазора двойственности для контроля точности и устойчивости решения.
Теорема 4.2.1. Пусть $X$ - выпуклый компакт и $f(x)$ - дифференцируемая функция на $X$ с Липшицевым градиентом. Шаг выбирается по правилу Армихо. Тогда для любого ${\color{red}{x_0 \in X}}$
Теорема (прямая).(Convex Optimization: Algorithms and Complexity, Th 3.8.) Пусть $f$ выпуклая и дифференцируемая функция и её градиент Липшицев с константой $L$. Множество $X$ - выпуклый компакт диаметра $d > 0$. Тогда метод условного градиента с шагом $\alpha_k = \frac{2}{k + 1}$ сходится как
$$ f(x^*) - f(x_k) \leq \dfrac{2d^2L}{k + 2}, \quad k \geq 1 $$Теорема (двойственная) см. эту статью. После выполнения $K$ итераций метода условного градиента для выпуклой и непрерывно дифференцируемой функции для функции $g$ и любого $k \leq K$ выполнено
$$ g(x_k) \leq \frac{2\beta C_f}{K+2} (1 + \delta), $$где $\beta \approx 3$, $\delta$ - точность решения промежуточных задач, $C_f$ - оценка кривизны $f$ на множестве $D$
$$ C_f = \sup_{x, s \in D; \gamma \in [0,1]} \frac{2}{\gamma^2}\left(f(x + \gamma(s - x)) - f(x) - \langle \gamma(s - x), f'(x)\rangle\right) $$Аргумент супремума так же известен как дивергенция Брегмана.
Определение. Atomic norm называется следующая функция
$$ \|x\|_{\mathcal{D}} = \inf_{t \geq 0} \{ t | x \in t\mathcal{D} \} $$Она является нормой, если симметрична и $0 \in \mathrm{int}(\mathcal{D})$
Таблица взята из статьи
In [1]:
def func(x, A, b):
return 0.5 * np.linalg.norm(A.dot(x) - b)**2
f = lambda x: func(x, A, b)
def grad_f(x, A, b):
grad = A.T.dot(A.dot(x) - b)
return grad
grad = lambda x: grad_f(x, A, b)
In [2]:
def linsolver(gradient):
x = np.zeros(gradient.shape[0])
pos_grad = gradient > 0
neg_grad = gradient < 0
x[pos_grad] = np.zeros(np.sum(pos_grad == True))
x[neg_grad] = np.ones(np.sum(neg_grad == True))
return x
In [3]:
def projection(y):
return np.clip(y, 0, 1)
In [4]:
import liboptpy.constr_solvers as cs
import liboptpy.step_size as ss
import numpy as np
from tqdm import tqdm
n = 200
m = 100
A = np.random.randn(m, n)
x_true = np.random.rand(n)
b = A.dot(x_true) + 0.01 * np.random.randn(m)
eigvals = np.linalg.eigvalsh(A.T @ A)
L = np.max(eigvals)
In [5]:
def myplot(x, y, xlab, ylab, xscale="linear", yscale="log"):
plt.figure(figsize=(10, 8))
plt.xscale(xscale)
plt.yscale(yscale)
for key in y:
plt.plot(x[key], y[key], label=key)
plt.xticks(fontsize=24)
plt.yticks(fontsize=24)
plt.legend(loc="best", fontsize=24)
plt.xlabel(xlab, fontsize=24)
plt.ylabel(ylab, fontsize=24)
In [6]:
x0 = np.random.rand(n)
cg = cs.FrankWolfe(f, grad, linsolver, ss.ConstantStepSize(1e-2))
x_cg = cg.solve(x0=x0, max_iter=200, tol=1e-10, disp=1)
print("Optimal value CG =", f(x_cg))
In [7]:
pg = cs.ProjectedGD(f, grad, projection, ss.ConstantStepSize(1 / L))
x_pg = pg.solve(x0=x0, max_iter=200, tol=1e-10, disp=1)
print("Optimal value PG =", f(x_pg))
In [8]:
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)
y_hist_f_cg = [f(x) for x in cg.get_convergence()]
y_hist_f_pg = [f(x) for x in pg.get_convergence()]
myplot({"CG": range(1, len(y_hist_f_cg) + 1), "PG": range(1, len(y_hist_f_pg) + 1)},
{"CG": y_hist_f_cg, "PG": y_hist_f_pg}, "Number of iteration",
r"Objective function, $\frac{1}{2}\|Ax - b\|^2_2$")
In [9]:
import cvxpy as cvx
x = cvx.Variable(n)
obj = cvx.Minimize(0.5 * cvx.norm(A * x - b, 2)**2)
constr = [x >= 0, x <= 1]
problem = cvx.Problem(objective=obj, constraints=constr)
value = problem.solve(solver=cvx.SCS)
x_cvx = np.array(x.value).ravel()
print("CVX optimal value =", value)
In [20]:
eps = [10**(-i) for i in range(8)]
time_pg = np.zeros(len(eps))
time_cg = np.zeros(len(eps))
iter_pg = np.zeros(len(eps))
iter_cg = np.zeros(len(eps))
pg = cs.ProjectedGD(f, grad, projection, ss.ConstantStepSize(1 / L))
cg = cs.FrankWolfe(f, grad, linsolver, ss.ConstantStepSize(1 / L))
for i, tol in tqdm(enumerate(eps)):
res = %timeit -o -q pg.solve(x0=x0, tol=tol, max_iter=100000)
time_pg[i] = res.average
iter_pg[i] = len(pg.get_convergence())
res = %timeit -o -q cg.solve(x0=x0, tol=tol, max_iter=100000)
time_cg[i] = res.average
iter_cg[i] = len(cg.get_convergence())
In [21]:
myplot({"CG":eps, "PG": eps}, {"CG": time_cg, "PG": time_pg}, r"Accuracy, $\varepsilon$", "Time, s", xscale="log")
In [22]:
myplot({"CG":eps, "PG": eps}, {"CG": iter_cg, "PG": iter_pg}, r"Accuracy, $\varepsilon$", "Number of iterations", xscale="log")
Pro
Contra