Исходную задачу можно преобразовать к задаче безусловной оптимизации
$$ \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 = 10
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()
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()
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.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.yticks(fontsize=fontsize)
Pro
Contra
Мотивация: для метода штрафных функций решения подзадач могут нарушать ограничения, и известно только, что
$$ 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 [5]:
import cvxpy
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rc("text", usetex=True)
num_iters = 10
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()
print ("Optimal value from CVXPY =", f.value)
In [14]:
# Solve with method of augmented Lagrangian
resid = A*x - b
y = cvxpy.Parameter(shape=(m, 1))
y.value = np.zeros((m, 1))
aug_lagr = f + y.T*resid + (mu/2.)*cvxpy.sum_squares(resid)
for t in range(num_iters):
cvxpy.Problem(cvxpy.Minimize(aug_lagr)).solve()
print(y.value, resid.value)
y = cvxpy.Parameter(shape=(m, 1), value=y.value + mu * resid.value)
# temp = y.value
# y.value = temp + mu * resid.value
print("Optimal value from method of augmented Lagrangian =", f.value)
In [3]:
mus = [0.1, 1, 10, 100]
conv_res = {}
conv_obj = {}
for mu in mus:
conv_res[mu] = np.zeros(num_iters)
conv_obj[mu] = np.zeros(num_iters)
x = cvxpy.Variable(shape=(n, 1))
f = cvxpy.norm(x, 2)
resid = A*x - b
y = cvxpy.Parameter(shape=(m, 1))
y.value = np.zeros((m, 1))
aug_lagr = f + y.T*resid + (mu/2.)*cvxpy.sum_squares(resid)
for t in range(num_iters):
cvxpy.Problem(cvxpy.Minimize(aug_lagr)).solve()
y.value += mu*resid.value
conv_res[mu][t] = np.linalg.norm(resid.value)
conv_obj[mu][t] = aug_lagr.value
In [4]:
plt.figure(figsize=(8, 6))
fontsize=20
for mu in mus:
plt.plot(conv_obj[mu], label=r"$\mu = {}$".format(mu))
plt.legend(fontsize=fontsize)
plt.xlabel("Number of iterations, $k$", fontsize=fontsize)
plt.ylabel("Augmented Lagrangian, $M_k$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)
In [5]:
plt.figure(figsize=(8, 6))
fontsize=20
for mu in mus:
plt.semilogy(conv_res[mu], label=r"$\mu = {}$".format(mu))
plt.legend(fontsize=fontsize)
plt.xlabel("Number of iterations, $k$", fontsize=fontsize)
plt.ylabel("Norm of residuals, $\|Ax_k - b\|_2$", fontsize=fontsize)
plt.xticks(fontsize=fontsize)
_ = plt.yticks(fontsize=fontsize)