Исходная задача \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 [26]:
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 [27]:
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 [28]:
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 [29]:
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 [30]:
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 [31]:
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 [32]:
scopt.check_grad(f, gradf, np.random.rand(n), c, mu)
Out[32]:
In [33]:
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 [34]:
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 [35]:
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 [36]:
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 [37]:
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)
res = scopt.linprog(c, A, b, bounds=bounds, options={"maxiter": 2**max(n_list) + 1})
time = %timeit -o -q scopt.linprog(c, A, b, bounds=bounds, options={"maxiter": 2**max(n_list) + 1})
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 [38]:
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)
In [39]:
def generate_linprog_problem(n, m=10):
x = np.random.rand(n)
A = np.random.randn(m, n)
b = A.dot(x)
c = np.random.randn(n)
c[c < 0] = 1
return c, A, b, x
In [40]:
n = 20
c, A, b, x0 = generate_linprog_problem(n)
res = scopt.linprog(c, A_eq=A, b_eq=b, bounds=[(0, None) for i in range(n)])
m = A.shape[0]
# x0 = np.random.rand(n)
# resid = b - A[:, m:].dot(x0[m:])
# x0[:m] = np.linalg.solve(A[:m, :m], resid)
print(np.linalg.norm(A.dot(x0) - b))
mu0 = 1
rho_mu = 0.5
In [41]:
x_bar_simple = BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking,
simple_solver, max_iter=100, tol=1e-8)
x_bar_elim = BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking,
elimination_solver, tol=1e-8)
print(np.linalg.norm(res["x"] - x_bar_simple))
print(np.linalg.norm(res["x"] - x_bar_elim))
print(c.dot(res["x"]))
print(c.dot(x_bar_elim))
In [42]:
n_list = [10*i for i in range(2, 63, 10)]
times_simple = np.zeros(len(n_list))
times_elim = np.zeros(len(n_list))
mu0 = 5
rho_mu = 0.5
for i, n in enumerate(n_list):
print("Current dimension = {}".format(n))
c, A, b, x0 = generate_linprog_problem(n)
time = %timeit -o BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, simple_solver)
times_simple[i] = time.average
time = %timeit -o BarrierPrimalLinConstr(f, gradf, hessf, A, c, x0, mu0, rho_mu, backtracking, elimination_solver)
times_elim[i] = time.average
In [43]:
plt.semilogy(n_list, times_elim, label="Elimination solver")
plt.semilogy(n_list, times_simple, label="Simple solver")
# plt.semilogy(n_list, 10**(-6)*np.array(n_list)**3 / 3)
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)