Исходную задачу можно преобразовать к задаче безусловной оптимизации
$$ \min_x Q(x, \mu), $$где
$$ Q(x, \mu) = f(x) + \frac{\mu}{2}\sum\limits_{i=1}^mg^2_i(x), \quad \mu > 0 $$def QudraticPenaltyEquality(Q, gradQ, x0, get_mu, get_tau, **kwargs):
while True:
# Stop when norm of gradient of Q is less than current tau
x = MinimizeQ(Q, gradQ, x0, get_tau)
if global_convergence(x, Q, **kwargs):
break
mu = get_mu()
Q, gradQ = UpdateQ(mu)
x0 = UpdateStartPoint(x, Q)
return x
Теорема 1. Пусть для каждого $\mu$ безусловная задача имеет конечное глобальное решение. Тогда предельная точка последовательности решений безусловных задач при $\mu \to \infty$ является глобальным решением исходной задачи.
Теорема 2. Пусть $\tau_k \to 0$ и $\mu_k \to \infty$ и $\| Q'(x^*_k, \mu_k) \| \leq \tau_k$. Тогда
Для любой подпоследовательности $x^*_k \to x^*, \; k \in \mathcal{C}$ выполнено
$$ \lim_{k \in \mathcal{C}} \mu_k g_i(x_k) = \lambda_i^* $$для всех $i = 1,\ldots,m$, где $\lambda_i^*$ множители Лагранжа, удовлетсворяющие ККТ.
Пример. \begin{equation*} \begin{split} & \min -5x_1^2 + x_2^2\\ \text{s.t. }& x_1 = 1 \end{split} \end{equation*} Штрафная функция примет вид
$$ Q(x, \mu) = -5x_1^2 + x_2^2 + \frac{\mu}{2}(x_1 - 1)^2 $$Наблюдение При $\mu < 10$ функция $Q(x, \mu)$ неограничена снизу по $x$, и подзадача не имеет конечного решения
где $g'(x)$ - якобиан вектор-функции ограничений-равенств.
Около точки минимума в силу теоремы 2 справедливо следующее приближение
$$ Q''(x, \mu) \approx L''(x, \lambda^*) + \mu(g'(x))^{\top} g'(x) $$Итог: некоторые собственные значения гессиана $Q''(x, \mu)$ имеют порядок $\mu$, а некоторые от $\mu$ не зависят, что приводит к вырожденности при увеличении $\mu$
Следствие: поиск направления в методе Ньютона очень неточный
Введём новую переменную $\xi = \mu g'(x) p$ и запишем эту систему в виде $$\begin{bmatrix} f''(x) + \mu \sum\limits_{i=1}^m g_i(x)g_i''(x) & (g'(x))^{\top}\\ (g'(x)) & -\frac{1}{\mu} I \end{bmatrix} = \begin{bmatrix} -Q'(x, \mu) \\ 0 \end{bmatrix}$$
где $\mu > 0$
In [1]:
import cvxpy
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)
num_iters = 30
n = 20
m = 10
A = np.random.randn(m, n)
b = np.random.randn(m, 1)
# Initialize problem
x = cvxpy.Variable(shape=(n, 1))
f = cvxpy.norm(x, 2)
# Solve with CVXPY.
cvxpy.Problem(cvxpy.Minimize(f), [A*x == b]).solve(solver=cvxpy.SCS)
cvxpy_f = f.value
print ("Optimal value from CVXPY =", f.value)
In [2]:
# Solve with method of augmented Lagrangian
resid = A*x - b
mu = cvxpy.Parameter(value=1, nonneg=True)
penal_f = f + (mu/2.)*cvxpy.sum_squares(resid)
res_conv = []
f_conv = []
for t in range(num_iters):
cvxpy.Problem(cvxpy.Minimize(penal_f)).solve(solver=cvxpy.SCS)
mu.value = mu.value * 2
res_conv.append(np.linalg.norm(resid.value))
f_conv.append(f.value)
print("Optimal value from method of augmented Lagrangian =", f.value)
In [3]:
plt.figure(figsize=(8, 6))
fontsize=20
plt.plot(f_conv, label="Penalty")
plt.plot(np.ones(num_iters) * cvxpy_f, label="CVXPy")
plt.legend(fontsize=fontsize)
plt.xlabel("Number of iterations, $k$", fontsize=fontsize)
plt.ylabel("Objective, $f(x_k)$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
plt.grid(True)
_ = plt.yticks(fontsize=fontsize)
In [4]:
plt.figure(figsize=(8, 6))
fontsize=20
plt.semilogy(res_conv)
plt.xlabel("Number of iterations, $k$", fontsize=fontsize)
plt.ylabel("Norm of residuals, $\|Ax_k - b\|_2$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
plt.grid(True)
_ = plt.yticks(fontsize=fontsize)
Pro
Contra
где $g(\lambda) = \inf_x L(x, \lambda)$
где $\hat{x} = \arg\min_x L(x, \lambda_k)$
In [5]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)
import cvxpy as cvx
def dual_ascent(update_x, A, b, alpha, x0, lambda0, max_iter):
x = x0.copy()
lam = lambda0.copy()
conv_x = [x]
conv_lam = [lam]
for i in range(max_iter):
x = update_x(x, lam, A, b)
lam = lam + alpha * (A @ x - b)
conv_x.append(x.copy())
conv_lam.append(lam.copy())
return x, lam, conv_x, conv_lam
In [6]:
m, n = 10, 20
A = np.random.randn(m, n)
b = np.random.randn(m)
P = np.random.randn(n, n)
P = P.T @ P
c = np.random.randn(n)
spec = np.linalg.eigvalsh(P)
mu = spec.min()
print(mu)
x = cvx.Variable(n)
obj = 0.5 * cvx.quad_form(x, P) - c @ x
problem = cvx.Problem(cvx.Minimize(obj), [A @ x == b])
problem.solve(verbose=True)
print(np.linalg.norm(A @ x.value - b))
print(problem.value)
In [7]:
x0 = np.random.randn(n)
lam0 = np.random.randn(m)
max_iter = 100000
alpha = 1e-4
def f(x):
return 0.5 * x @ P @ x - c @ x
def L(x, lam):
return f(x) + lam @ (A @ x - b)
def update_x(x, lam, A, b):
return np.linalg.solve(P, c - A.T @ lam)
x_da, lam_da, conv_x_da, conv_lam_da = dual_ascent(update_x, A, b, alpha, x0, lam0, max_iter)
print(np.linalg.norm(A @ x_da - b))
print(0.5 * x_da @ P @ x_da - c @ x_da)
In [8]:
plt.figure(figsize=(10, 8))
plt.plot([f(x) for x in conv_x_da], label="Objective")
plt.plot(problem.value * np.ones(len(conv_x_da)), label="Traget value")
# plt.yscale("log")
plt.xscale("log")
plt.legend(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
In [9]:
plt.plot([L(x, lam) for x, lam in zip(conv_x_da, conv_lam_da)],
label="Lagrangian")
plt.legend(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
Out[9]:
In [10]:
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_da], label="$\|Ax - b\|_2$")
plt.legend(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.grid(True)
Мотивация: для метода штрафных функций решения подзадач могут нарушать ограничения, и известно только, что
$$ g_i(x^*_k) \approx \frac{\lambda^*}{\mu_k} \to 0, \quad \mu_k \to \infty $$Можно ли изменить $Q(x, \mu)$ так, чтобы избежать этого нарушения ограничений?
Идея: добавлять штраф не к целевой функции, а к функции Лагранжа. Аналог прямо-двойственного метода, так как за одну итерацию происходит обновление как прямых, так и двойственных переменных
Необходимое условие минимума $M(x_k, \lambda^k, \mu_k)$
$$ f'(x_k) + \sum\limits_{i=1}^m (\lambda^k_i + \mu_k g_i(x_k) ) g'_i(x_k) \approx 0 $$Отсюда получаем выражение для $\lambda^{k+1}$
$$ \lambda^{k+1}_i = \lambda^k_i + \mu_k g_i(x_k) $$где $\lambda$ - двойственные переменные для ограничений-равенств, $\nu$ - двойственная переменная для ограничений-неравенств.
In [11]:
def augmented_lagrangian(update_x, A, b, rho0, x0, lambda0, max_iter):
x = x0.copy()
lam = lambda0.copy()
conv_x = [x]
conv_lam = [lam]
rho = rho0
for i in range(max_iter):
x = update_x(x, lam, A, b)
lam = lam + rho * (A @ x - b)
conv_x.append(x.copy())
conv_lam.append(lam.copy())
return x, lam, conv_x, conv_lam
In [12]:
def update_x_al(x, lam, A, b):
return np.linalg.solve(P + rho * A.T @ A, c - A.T @ lam + A.T @ b)
rho = 10
max_iter = 1000
x_al, lam_al, conv_x_al, conv_lam_al = augmented_lagrangian(update_x_al, A, b, rho, x0, lam0, max_iter)
print(np.linalg.norm(A @ x_al - b))
print(0.5 * x_al @ P @ x_al - c @ x_al)
In [13]:
plt.plot([f(x) for x in conv_x_da], label="DA")
plt.plot([f(x) for x in conv_x_al], label="AL")
# plt.yscale("log")
plt.xscale("log")
plt.legend(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("Objective", fontsize=20)
Out[13]:
In [14]:
plt.plot([L(x, lam) for x, lam in zip(conv_x_da, conv_lam_da)],
label="DA")
plt.plot([L(x, lam) for x, lam in zip(conv_x_al, conv_lam_al)],
label="AL")
plt.legend(fontsize=20)
plt.xscale("log")
plt.xlabel("\# iterations", fontsize=20)
plt.xlabel("Lagrangian", fontsize=20)
Out[14]:
In [15]:
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_da], label="DA")
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_al], label="AL")
plt.legend(fontsize=20)
plt.xscale("log")
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("$\|Ax - b\|_2$", fontsize=20)
plt.grid(True)
plt.yticks(fontsize=20)
plt.xticks(fontsize=20)
Out[15]:
Задача станет такой
\begin{align*} & \min f(x) + I_{Ax = b} (z)\\ \text{s.t. } & x = z \end{align*}Для неё модифицированный лагранжиан примет вид
$$ L_{\rho}(x, z, \lambda) = f(x) + I_{Ax = b} (z) + \lambda^{\top}(x - z) + \frac{\rho}{2}\|x - z\|_2^2 $$
In [16]:
def admm(update_x, update_z, rho0, x0, z0, lambda0, max_iter):
x = x0.copy()
z = z0.copy()
lam = lambda0.copy()
conv_x = [x]
conv_z = [z]
conv_lam = [lam]
rho = rho0
for i in range(max_iter):
x = update_x(x, z, lam, A, b)
z = update_z(x, z, lam, A, b)
lam = lam + rho * (x - z)
conv_x.append(x.copy())
conv_z.append(z.copy())
conv_lam.append(lam.copy())
return x, z, lam, conv_x, conv_z, conv_lam
In [17]:
def update_x_admm(x, z, lam, A, b):
n = x.shape[0]
return np.linalg.solve(P + rho*np.eye(n), -lam + c + rho * z)
def update_z_admm(x, z, lam, A, b):
x_hat = lam / rho + x
return x_hat - A.T @ np.linalg.solve(A @ A.T, A @ x_hat - b)
In [18]:
z0 = np.random.randn(n)
lam0 = np.random.randn(n)
rho = 1000
x_admm, z_admm, lam_admm, conv_x_admm, conv_z_admm, conv_lam_admm = admm(update_x_admm,
update_z_admm,
rho, x0, z0, lam0,
max_iter=10)
print(f(x_admm))
In [19]:
plt.figure(figsize=(10, 8))
plt.plot([f(x) for x in conv_x_da], label="DA")
plt.plot([f(x) for x in conv_x_al], label="AL")
plt.plot([f(x) for x in conv_x_admm], label="ADMM x")
plt.plot([f(z) for z in conv_z_admm], label="ADMM z")
# plt.yscale("log")
plt.xscale("log")
plt.legend(fontsize=20)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("Objective", fontsize=20)
plt.grid(True)
plt.yticks(fontsize=20)
plt.xticks(fontsize=20)
Out[19]:
In [20]:
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_da], label="DA")
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_al], label="AL")
plt.semilogy([np.linalg.norm(A @ x - b) for x in conv_x_admm], label="ADMM")
plt.legend(fontsize=20)
plt.xscale("log")
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("$\|Ax - b\|_2$", fontsize=20)
plt.grid(True)
plt.yticks(fontsize=20)
plt.xticks(fontsize=20)
plt.show()
In [21]:
plt.semilogy([np.linalg.norm(x - z) for x, z in zip(conv_x_admm, conv_z_admm)])
plt.grid(True)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("$\|x_k - z_k\|_2$", fontsize=20)
plt.yticks(fontsize=20)
plt.show()
где $u_k = \lambda_k / \rho$
где $c^{\top}x$ определена на множестве $Ax = b$.
In [33]:
import scipy.optimize as scopt
m, n = 10, 200
A = np.random.rand(m, n)
b = np.random.rand(m)
c = np.random.rand(n)
scipy_linprog_conv = []
def callback_splin(cur_res):
scipy_linprog_conv.append(cur_res)
res = scopt.linprog(c, A_eq=A, b_eq=b,
bounds=[(0, None) for i in range(n)],
callback=callback_splin, method="simplex")
print(res)
In [34]:
def update_x_admm(x, z, lam, A, b):
n = x.shape[0]
m = A.shape[0]
C = np.block([[rho * np.eye(n), A.T], [A, np.zeros((m, m))]])
rhs = np.block([-lam - c + rho * z, b])
return np.linalg.solve(C, rhs)[:n]
def update_z_admm(x, z, lam, A, b):
x_hat = lam / rho + x
return np.clip(x_hat, 0, np.max(x_hat))
In [42]:
x0 = np.random.randn(n)
z0 = np.random.randn(n)
lam0 = np.random.randn(n)
rho = 1
x_admm, z_admm, lam_admm, conv_x_admm, conv_z_admm, conv_lam_admm = admm(update_x_admm,
update_z_admm,
rho, x0, z0, lam0, max_iter=10000)
print(c @ x_admm, res.fun)
In [43]:
print(c @ x_admm - res.fun, np.linalg.norm(x_admm - res.x))
In [44]:
plt.figure(figsize=(10, 8))
plt.plot([c @ x for x in conv_x_admm[:100]], label="ADMM")
plt.plot([c @ res.x for res in scipy_linprog_conv], label="Scipy")
plt.legend(fontsize=20)
plt.grid(True)
plt.xlabel("\# iterations", fontsize=20)
plt.ylabel("$c^{\\top}x_k$", fontsize=20)
plt.yticks(fontsize=20)
plt.xticks(fontsize=20)
Out[44]: