Следующий шаг: произвольная достаточно гладкая функция на достаточно простом множестве - не обязательно полиэдральном.
Определение. Множество будем называть простым, если проекцию на него можно найти существенно быстрее (чаще всего аналитически) по сравнению с решением исходной задачи минимизации.
Идея: делать шаг градиентного спуска и проецировать полученную точку на допустимое множество $P$.
Теорема. Пусть $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 = 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(b)
grad = grad + A.T.dot(A.dot(x))
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)
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.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
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.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
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()
x_cvx = np.array(x.value).ravel()
print("CVX optimal value =", value)
In [11]:
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.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
cg = cs.FrankWolfe(f, grad, linsolver, ss.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
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 [12]:
myplot({"CG":eps, "PG": eps}, {"CG": time_cg, "PG": time_pg}, r"Accuracy, $\varepsilon$", "Time, s", xscale="log")
In [13]:
myplot({"CG":eps, "PG": eps}, {"CG": iter_cg, "PG": iter_pg}, r"Accuracy, $\varepsilon$", "Number of iterations", xscale="log")
In [14]:
def linsolver(gradient):
x = np.zeros(gradient.shape[0])
idx_min = np.argmin(gradient)
if gradient[idx_min] > 0:
x[idx_min] = 0
else:
x[idx_min] = 1
return x
In [15]:
def projection(y):
x = y.copy()
if np.all(x >= 0) and np.sum(x) <= 1:
return x
x = np.clip(x, 0, np.max(x))
if np.sum(x) <= 1:
return x
n = x.shape[0]
bget = False
x.sort()
x = x[::-1]
temp_sum = 0
t_hat = 0
for i in range(n - 1):
temp_sum += x[i]
t_hat = (temp_sum - 1.0) / (i + 1)
if t_hat >= x[i + 1]:
bget = True
break
if not bget:
t_hat = (temp_sum + x[n - 1] - 1.0) / n
return np.maximum(y - t_hat, 0)
In [30]:
x0 = np.random.rand(n) * 10
x0 = x0 / x0.sum()
cg = cs.FrankWolfe(f, grad, linsolver, ss.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
x_cg = cg.solve(x0=x0, max_iter=200, tol=1e-10)
print("Optimal value CG =", f(x_cg))
In [31]:
pg = cs.ProjectedGD(f, grad, projection, ss.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
x_pg = pg.solve(x0=x0, max_iter=200, tol=1e-10)
print("Optimal value PG =", f(x_pg))
In [32]:
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 [22]:
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.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
cg = cs.FrankWolfe(f, grad, linsolver, ss.Backtracking(rule_type="Armijo", rho=0.5, beta=0.1, init_alpha=1.))
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 [23]:
myplot({"CG":eps, "PG": eps}, {"CG": time_cg, "PG": time_pg},
r"Accuracy, $\varepsilon$", "Time, s", xscale="log")
In [24]:
myplot({"CG":eps, "PG": eps}, {"CG": iter_cg, "PG": iter_pg},
r"Accuracy, $\varepsilon$", "Number of iterations", xscale="log")
In [25]:
x = cvx.Variable(n)
obj = cvx.Minimize(0.5 * cvx.norm2(A * x - b)**2)
constr = [cvx.norm(x, 1) <= 1, x >= 0]
problem = cvx.Problem(objective=obj, constraints=constr)
value = problem.solve()
x_cvx = np.array(x.value).ravel()
print("CVX optimal value =", value)
Pro
Contra