Двойственная функция
\begin{equation*} \begin{split} g(\mu) & = -b^{\top}\mu + \inf_x(f(x) + \mu^{\top}Ax) \\ & = -b^{\top}\mu - \sup_x((-A^{\top}\mu)^{\top}x -f(x)) \\ & = -b^{\top}\mu - f^*(-A^{\top}\mu) \end{split} \end{equation*}Двойственная задача
$$ \max_\mu -b^{\top}\mu - f^*(-A^{\top}\mu) $$Подход 1: найти сопряжённую функцию и решить безусловную задачу оптимизации
Трудности
или
$$ \begin{bmatrix} f' & A^{\top} \\ A & 0 \end{bmatrix} \begin{bmatrix} x^{\\*} \\ \mu^{\\*} \end{bmatrix} = \begin{bmatrix} 0 \\ b \end{bmatrix} $$Подход 2: решить нелинейную в общем случае систему методом Ньютона.
Вопрос: в каком случае система окажется линейной?
Из условий оптимальности имеем
$$ \begin{bmatrix} f''(x) & A^{\top} \\ A & 0 \end{bmatrix} \begin{bmatrix} v \\ w \end{bmatrix} = \begin{bmatrix} -f'(x) \\ 0 \end{bmatrix} $$Шаг метода Ньютона определён только для невырожденной матрицы!
Упражнение. Посчитайте за сколько итераций метод Ньютона сойдётся для квадратичной функции с ограничениями типа равенств.
Важно: начальная точка должна лежать в допустимом множестве!
def NewtonEqualityFeasible(f, gradf, hessf, A, b, stop_crit, line_search, x0, tol):
x = x0
n = x.shape[0]
while True:
newton_matrix = [[hessf(x), A.T], [A, 0]]
rhs = [-gradf(x), 0]
w = solve_lin_sys(newton_matrix, rhs)
h = w[:n]
if stop_crit(x, h, gradf(x), **kwargs) < tol:
break
alpha = line_search(x, h, f, gradf(x), **kwargs)
x = x + alpha * h
return x
Получим выражение для значения
$$ f(x) - \inf_v(\hat{f}(x + v) \; | \; A(x+v) = b), $$где $\hat{f}$ - квадратичная аппроксимация функции $f$.
Для этого
$$ \langle h^{\top} \rvert \cdot \quad f''(x)h + A^{\top}w = -f'(x) $$с учётом $Ah = 0$ получаем
$$ h^{\top}f''(x)h = -f'(x)^{\top}h $$Тогда
$$ \inf_v(\hat{f}(x + v) \; | \; A(x+v) = b) = f(x) - \frac{1}{2}h^{\top}f''(x)h $$Вывод: величина $h^{\top}f''(x)h$ является наиболее адекватным критерием остановки метода Ньютона.
Сходимость метода аналогична сходимости метода Ньютона для задачи безусловной оптимизации.
Теорема Пусть выполнены следующие условия
Тогда, метод Ньютона сходится к паре $(x^*, \mu^*)$ линейно, а при достижении достаточной близости к решению - квадратично.
где $r_p(x, \mu) = Ax - b$ и $r_d(x, \mu) = f'(x) + A^{\top}\mu$
или более подробно
$$ \begin{bmatrix} f''(x) & A^{\top}\\ A & 0 \end{bmatrix} \begin{bmatrix} z_p\\ z_d \end{bmatrix} = - \begin{bmatrix} r_p(x, \mu)\\ r_d(x, \mu) \end{bmatrix} = - \begin{bmatrix} f'(x) + A^{\top}\mu\\ Ax - b \end{bmatrix} $$def NewtonEqualityInfeasible(f, gradf, hessf, A, b, stop_crit, line_search, x0, mu0, tol):
x = x0
mu = mu0
n = x.shape[0]
while True:
z_p, z_d = ComputeNewtonStep(hessf(x), A, b)
if stop_crit(x, z_p, z_d, gradf(x), **kwargs) < tol:
break
alpha = line_search(x, z_p, z_d, f, gradf(x), **kwargs)
x = x + alpha * z_p
mu = z_d
return x
def linesearch(r, x, mu, z_p, z_d, c, beta):
alpha = 1
while np.linalg.norm(r(x + alpha * z_p, mu + alpha * z_d)) >= (1 - c * alpha) * np.linalg.norm(r(x, mu)):
alpha *= beta
return alpha
Результат аналогичен результаты для допустимой начальной точки
Теорема. Пусть
Тогда сходимость метода линейна при удалении от решении и квадратичная при достаточном приближении к решению.
где $f_i$ - выпуклые и дважды непрерывно дифференцируемы, $A \in \mathbb{R}^{p \times n}$ и $\mathrm{rank} \; A = p < n$.
Предполагаем, что задача строго разрешима, то есть выполняется условие Слейтера.
где $I_-$ - индикаторная функция
$$ I_-(u) = \begin{cases} 0, & u \leq 0\\ \infty, & u > 0 \end{cases} $$Проблема. Теперь целевая функция - недифференцируема.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
x = np.linspace(-2, 0, 100000, endpoint=False)
plt.figure(figsize=(10, 6))
for t in [0.1, 0.5, 1, 1.5, 2]:
plt.plot(x, -t * np.log(-x), label=r"$t = " + str(t) + "$")
plt.legend(fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel("u", fontsize=20)
Out[1]:
называется логарифмическим барьером. Её область определения - множество точек, для котороых ограничения типа неравенств выполняются строго.
Упражнение. Найдите градиент и гессиан $\phi(x)$
для $\lambda = \lambda^*(t)$ и $\mu = \mu^*$.
def BarrierMethod(f, x0, t0, tol, alpha, **kwargs):
x = x0
t = t0
while True:
x = SolveBarrierProblem(f, t, x, **kwargs)
if m * t < tol:
break
t *= alpha
return x
Простой поиск допустимой точки
\begin{equation*} \begin{split} & \min s\\ \text{s.t. } & f_i(x) \leq s\\ & Ax = b \end{split} \end{equation*}Похож на барьерный метод, но
При некоторых предположениях о начальной точке и начальном значении $\mu$, можно показать, что для достижения $\mu_k \leq \varepsilon$ потребуется
$$ \mathcal{O}\left(\sqrt{n}\log \left( \frac{1}{\varepsilon}\right)\right) $$итераций
Доказательство и все детали можно найти тут или тут
Исходная задача \begin{align*} &\min_x c^{\top}x \\ \text{s.t. } & Ax = b\\ & x_i \geq 0, \; i = 1,\dots, n \end{align*}
Аппроксимированная задача \begin{align*} &\min_x c^{\top}x {\color{red}{- \mu \sum\limits_{i=1}^n \ln x_i}} \\ \text{s.t. } & Ax = b\\ \end{align*} для некоторого $\mu > 0$
Пусть $s = \mu X^{-1}e$, тогда условия оптимальности можно переписать так:
Также $x > 0 \Rightarrow s > 0$
Из первой строки:
$$ -\mu X^{-2}p + A^{\top}\lambda = c - \mu X^{-1}e $$$$ -\mu Ap + AX^{2}A^{\top}\lambda = AX^2c - \mu AXe $$$$ AX^{2}A^{\top}\lambda = AX^2c - \mu AXe $$Так как $X \in \mathbb{S}^n_{++}$ и $A$ полного ранга, то уравнение имеет единственное решение $\lambda^*$.
Последовательное исключение переменных:
$Ap = -AH^{-1}(g + A^{\top}\lambda) = -AH^{-1}A^{\top}\lambda - AH^{-1}g = 0$
Здесь матрица $-AH^{-1}A^{\top}$ есть дополнение по Шуру матрицы $H$.
def PrimalBarrierLP(c, A, b, x0, mu0, rho, tol):
x = x0
mu = mu0
e = np.ones(x.shape[0])
while True:
p, lam = ComputeNewtonDirection(c, x, A, mu)
alpha = line_search(p, mu, c, x)
x = x + alpha * p
mu = rho * mu
if mu < tol and np.linalg.norm(x.dot(c - A.T.dot(lam)) - mu * e) < tol:
break
return x
In [2]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import scipy.optimize as scopt
import scipy.linalg as sclin
USE_COLAB = False
if not USE_COLAB:
plt.rc("text", usetex=True)
In [3]:
def NewtonLinConstraintsFeasible(f, gradf, hessf, A, x0, line_search, linsys_solver, args=(),
disp=False, disp_conv=False, callback=None, tol=1e-6, max_iter=100, **kwargs):
x = x0.copy()
n = x0.shape[0]
iteration = 0
lam = np.random.randn(A.shape[0])
while True:
gradient, hess = gradf(x, *args), hessf(x, *args)
h = linsys_solver(hess, A, gradient)
descent_dir = h[:n]
decrement = descent_dir.dot(hessf(x, *args).dot(descent_dir))
if decrement < tol:
if disp_conv:
print("Tolerance achieved! Decrement = {}".format(decrement))
break
alpha = line_search(x, descent_dir, f, gradf, args, **kwargs)
if alpha < 1e-16:
if disp_conv:
print("Step is too small!")
x = x + alpha * descent_dir
if callback is not None:
callback((descent_dir, x))
iteration += 1
if disp:
print("Current function val = {}".format(f(x, *args)))
print("Newton decrement = {}".format(decrement))
if iteration >= max_iter:
if disp_conv:
print("Maxiter exceeds!")
break
res = {"x": x, "num_iter": iteration, "tol": decrement}
return res
In [4]:
def simple_solver(hess, A, gradient):
n = hess.shape[0]
n_lin_row, n_lin_col = A.shape
modified_hess = np.zeros((n + n_lin_row, n + n_lin_row))
modified_hess[:n, :n] = hess
modified_hess[n:n + n_lin_row, :n_lin_col] = A
modified_hess[:n_lin_col, n:n + n_lin_row] = A.T
rhs = np.zeros(n + n_lin_row)
rhs[:n] = -gradient
h = np.linalg.solve(modified_hess, rhs)
return h
def elimination_solver(hess, A, gradient):
inv_hess_diag = np.divide(1.0, np.diag(hess))
inv_hess_grad = np.multiply(-inv_hess_diag, gradient)
rhs = A.dot(inv_hess_grad)
L_inv_hess = np.sqrt(inv_hess_diag)
AL_inv_hess = A * L_inv_hess
# print(AL_inv_hess.shape)
S = AL_inv_hess.dot(AL_inv_hess.T)
# cho_S = sclin.cho_factor(S)
# w = sclin.cho_solve(cho_S, rhs)
w = np.linalg.solve(S, rhs)
v = np.subtract(inv_hess_grad, np.multiply(inv_hess_diag, A.T.dot(w)))
# h = np.zeros(hess.shape[1] + A.shape[0])
# h[:hess.shape[1]] = v
# h[hess.shape[1]:hess.shape[1] + A.shape[0]] = w
return v
In [5]:
def backtracking(x, descent_dir, f, grad_f, args, **kwargs):
beta1 = kwargs["beta1"]
rho = kwargs["rho"]
alpha = 1
while f(x + alpha * descent_dir, *args) >= f(x, *args) + beta1 * alpha * grad_f(x, *args).dot(descent_dir) \
or np.isnan(f(x + alpha * descent_dir, *args)):
alpha *= rho
if alpha < 1e-16:
break
return alpha
In [6]:
def generate_KleeMinty_test_problem(n):
c = np.array([2**i for i in range(n)])
c = -c[::-1]
bounds = [(0, None) for i in range(n)]
b = np.array([5**(i+1) for i in range(n)])
a = np.array([1] + [2**(i+1) for i in range(1, n)])
A = np.zeros((n, n))
for i in range(n):
A[i:, i] = a[:n-i]
return c, A, b, bounds
In [7]:
n = 7
c, A, b, _ = generate_KleeMinty_test_problem(n)
eps = 1e-10
def f(x, c, mu):
n = c.shape[0]
return c.dot(x[:n]) - mu * np.sum(np.log(eps + x))
def gradf(x, c, mu):
grad = np.zeros(len(x))
n = c.shape[0]
grad[:n] = c - mu / (eps + x[:n])
grad[n:] = -mu / (eps + x[n:])
return grad
def hessf(x, c, mu):
return mu * np.diag(1. / (eps + x)**2)
A_lin = np.zeros((n, n + A.shape[0]))
A_lin[:n, :n] = A
A_lin[:n, n:n + A.shape[0]] = np.eye(A.shape[0])
mu = 0.1
In [8]:
scopt.check_grad(f, gradf, np.random.rand(n), c, mu)
Out[8]:
In [9]:
x0 = np.zeros(2*n)
x0[:n] = np.random.rand(n)
x0[n:2*n] = b - A.dot(x0[:n])
print(np.linalg.norm(A_lin.dot(x0) - b))
print(np.sum(x0 <= 1e-6))
In [10]:
hist_conv = []
def cl(x):
hist_conv.append(x)
res = NewtonLinConstraintsFeasible(f, gradf, hessf, A_lin, x0, backtracking, elimination_solver, (c, mu), callback=cl,
max_iter=2000, beta1=0.1, rho=0.7)
print("Decrement value = {}".format(res["tol"]))
fstar = f(res["x"], c, mu)
hist_conv_f = [np.abs(fstar - f(descdir_x[1], c, mu)) for descdir_x in hist_conv]
plt.figure(figsize=(12, 5))
plt.subplot(1,2,1)
plt.semilogy(hist_conv_f)
plt.xlabel("Number of iteration, $k$", fontsize=18)
plt.ylabel("$f^* - f_k$", fontsize=18)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
hist_conv_x = [np.linalg.norm(res["x"] - x[1]) for x in hist_conv]
plt.subplot(1,2,2)
plt.semilogy(hist_conv_x)
plt.xlabel("Number of iteration, $k$", fontsize=18)
plt.ylabel("$\| x_k - x^*\|_2$", fontsize=18)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)
plt.tight_layout()
In [11]:
def BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, linesearch, linsys_solver,
tol=1e-8, max_iter=500, disp_conv=False, **kwargs):
x = x0.copy()
n = x0.shape[0]
mu = mu0
while True:
res = NewtonLinConstraintsFeasible(f, gradf, hessf, A, x, linesearch, linsys_solver, (c, mu),
disp_conv=disp_conv, max_iter=max_iter, beta1=0.01, rho=0.5)
x = res["x"].copy()
if n * mu < tol:
break
mu *= rho_mu
return x
In [12]:
mu0 = 5
rho_mu = 0.5
x = BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, elimination_solver, max_iter=100)
%timeit BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, elimination_solver, max_iter=100)
%timeit BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, simple_solver, max_iter=100)
print(x[:n])
In [13]:
mu0 = 2
rho_mu = 0.5
n_list = range(3, 10)
n_iters = np.zeros(len(n_list))
times_simplex = np.zeros(len(n_list))
times_barrier_simple = np.zeros(len(n_list))
for i, n in enumerate(n_list):
print("Current dimension = {}".format(n))
c, A, b, bounds = generate_KleeMinty_test_problem(n)
time = %timeit -o -q scopt.linprog(c, A, b, bounds=bounds, options={"maxiter": 2**max(n_list) + 1}, method="simplex")
times_simplex[i] = time.best
A_lin = np.zeros((n, n + A.shape[0]))
A_lin[:n, :n] = A
A_lin[:n, n:n + A.shape[0]] = np.eye(A.shape[0])
x0 = np.zeros(2*n)
x0[:n] = np.random.rand(n)
x0[n:2*n] = b - A.dot(x0[:n])
time = %timeit -o -q BarrierPrimalLinConstr(f, gradf, hessf, A_lin, c, x0, mu0, rho_mu, backtracking, simple_solver)
times_barrier_simple[i] = time.best
In [14]:
plt.figure(figsize=(8, 5))
plt.semilogy(n_list, times_simplex, label="Simplex")
plt.semilogy(n_list, times_barrier_simple, label="Primal barrier")
plt.legend(fontsize=18)
plt.xlabel("Dimension, $n$", fontsize=18)
plt.ylabel("Computation time, sec.", fontsize=18)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)