Given vectors $c \in \mathbb{R}^n$, $b \in \mathbb{R}^m$ and matrix $A \in \mathbb{R}^{m \times n}$ such that $m < n$ and $\mathrm{rank}(A) = m$
Producing optimal quantity of product with resources constraints: \begin{align*} &\max_x c^{\top}x \\ \text{s.t. } & Ax \leq b\\ & x_i \geq 0, \; i = 1,\dots, n, \end{align*} where $x_i$ - quantity of the $i$-th product, $c_i$ is a revenue from the $i$-th product, $b_i$ - available quantity of the $i$-th resource, $a_{ij}$ is a quantity of the $i$-th resource, which is required to produce unit of the $j$-th product.
Flows in networks: transport problem, max flow problem, minimal cost path in communication network to pass message
Regression problem in $\ell_1$ and $\ell_{\infty}$ norms can be formulated as linear programming
Questions:
Theorem
All extreme points of polytope correspond to vertices of the polytope.
Q: how formalize and perform these steps?
Given extreme point $x$, corresponding basis matrix $B$ and set of indices $\mathcal{B}$.
Compute $u = B^{-1}a_{j^*}$
if there are positive entries, then compute
$$ \theta^* = \min_{\{i | u_i > 0\}} \frac{x_{\mathcal{B}(i)}}{u_i} $$
Select such index $\ell$ that
$$ \theta^* = \frac{x_{\mathcal{B}(\ell)}}{u_{\ell}}. $$
Compose new basis matrix $\hat{B}$ through replacing column $a_{\mathcal{B}(\ell)}$ with column $a_{j^*}$. New extreme point $\hat{x}$, corresponding to the basis matrix $\hat{B}$, is defined as
$$ \begin{align*} & \hat{x}_{j^*} = \theta^*\\ & \hat{x}_{\mathcal{B}(k)} = x_{\mathcal{B}(k)} - \theta^*u_k, \text{if } k \neq \ell \end{align*} $$
Matrix $U$ has rank-one
The best total complexity is $O(m^2)$, if reduced costs are computed with pivoting, and the worst total complexity is $O(mn)$, if all reduced costs are computed.
Assume
Then simplex method stops after finite number of iterations and
gives one of the following result
Definition. Extreme point is called degenerate, if more than $n - m$ of its entries are zero
Q: what is geometric interpretation of degeneracy?
R. Bland is American methematician,
one of the author of the oriented matroids theory.
To find initial extreme point compose the following auxilliary problem assuming that $b_i \geq 0, \; i =1, \dots,m$. This assumption is easy to satisfied with multiplication corresponding rows of $A$ and elements of $b$ by $-1$. \begin{align*} & \min_{z, y} y_1 + \ldots + y_m \\ \text{s.t. } & Az + y = b\\ & z \geq 0, \; y \geq 0 \end{align*}
Idea: unite two phases of two-phase simplex method
into the single phase
\begin{align*} & \min_{z, y} c^{\top}z + M(y_1 + \ldots + y_m) \\ \text{s.t. } & Az + y = b\\ & z \geq 0, \; y \geq 0, \end{align*}where $M$ is arbitrary large positive real number.
Usually it is unknown in advance,
therefore you can use it as free parameter that can be done large enough if necessary
Examples are available here
In [2]:
import scipy.optimize as scopt
import numpy as np
n = 1000
m = 10
c = 10 * np.random.rand(n)
b = np.random.rand(m)
A = np.random.randn(m, n)
res = scopt.linprog(c, A, b, bounds=[(-1, None) for i in range(n)])
print(res)
In the following problem \begin{align*} & \max_{x \in \mathbb{R}^n} 2^{n-1}x_1 + 2^{n-2}x_2 + \dots + 2x_{n-1} + x_n\\ \text{s.t. } & x_1 \leq 5\\ & 4x_1 + x_2 \leq 25\\ & 8x_1 + 4x_2 + x_3 \leq 125\\ & \ldots\\ & 2^n x_1 + 2^{n-1}x_2 + 2^{n-2}x_3 + \ldots + x_n \leq 5^n\\ & x \geq 0 \end{align*} starting from the point $x_0 = 0$ and and following to trajectory of simplex method one has to visit $2^n - 1$ vertices.
Exercise: solve this problem for $n = 2$ and $n = 3$, and generalize solution for arbitrary $n$.
In [3]:
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 [4]:
n = 5
c, A, b, bounds = generate_KleeMinty_test_problem(n)
print(c)
print(A)
print(b)
print(bounds)
In [5]:
res = scopt.linprog(c, A, b, bounds=bounds)
print(res)
In [6]:
n_list = range(3, 16)
n_iters = np.zeros(len(n_list))
times = np.zeros(len(n_list))
for i, n in enumerate(n_list):
c, A, b, bounds = generate_KleeMinty_test_problem(n)
res = scopt.linprog(c, A, b, bounds=bounds, options={"maxiter": 2**max(n_list)})
time = %timeit -o scopt.linprog(c, A, b, bounds=bounds, options={"maxiter": 2**max(n_list) + 1})
n_iters[i] = res.nit
times[i] = time.best
In [7]:
USE_COLAB = False
%matplotlib inline
import matplotlib.pyplot as plt
if not USE_COLAB:
plt.rc("text", usetex=True)
plt.figure(figsize=(20,5))
plt.subplot(1, 2, 1)
plt.plot(n_list, n_iters - np.array([2**n - 1 for n in n_list]), label="$K_t - K_{exp}$")
# plt.semilogy(n_list, [2**n - 1 for n in n_list], label="Theory")
plt.xlabel("Dimension, $n$", fontsize=24)
plt.ylabel("Number of iterations, $K$", fontsize=24)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.legend(fontsize=18)
plt.subplot(1, 2, 2)
plt.semilogy(n_list, times)
plt.xlabel("Dimension, $n$", fontsize=24)
plt.ylabel("Computation time", fontsize=24)
plt.xticks(fontsize=18)
_ = plt.yticks(fontsize=18)