Linear Programming

Linear Programming
Linear Programming

In the following section, Set Theory, Topology and Nonstandard Analysis are presupposed. The exponential simplex and the polynomial centre method solve linear programmes (LPs).

Diameter theorem for polytopes: The diameter of an \(n\)-dimensional polytope defined by \(m\) constraints for \(m, n \in {}^{\omega}\mathbb{N}_{\ge 2}\) is at most \(2(m + n – 3)\).

Proof: At most \(\acute{m}\) hyperplanes can be assembled into an incomplete cycle of dimension 2 and there are at most \(n – 2\) alternatives sidewards in the remaining dimensions. Overcoming every minimal distance requires at most two edges and yields the factor 2.\(\square\)

Remark: Dropping the requirement of finiteness, the theorem can be extended to polyhedra analogously.

Definition: Let \(\vartheta := \omega\,{}_{e}\omega, \underline{u}_n := (u, …, u)^T \in{}^{\omega}\mathbb{R}^{n}\) and \({|| \cdot ||}_{1}\) the 1-norm. If \(x \in X \subseteq {}^{\omega}\mathbb{R}^{\omega}\) is passed to the next step, it is followed by \({}^{*}\) where \(\Delta x := {x}^{*} – x\). The unit vector of \(x\) is \({_1}{x} := x/||x||\), where \(_1{0}\) is undefined. Let \(x, y \in X \) in \((x, y, …)^T\) be row vectors. If a method requires computation time in seconds and memory in bits of \(\mathcal{O}({\omega}^{\mathcal{O}(1)})\), it is polynomial. If one of both quantities is \(\mathcal{O}({e}^{|\mathcal{O}(\omega)|})\), the method is exponential. Let the eigenproduct (yet: determinant) of a square matrix be the product of its eigenvalues.\(\triangle\)

Theorem: The simplex method is exponential.

Proof and algorithm: Let \(P := \{x \in {}^{\omega}\mathbb{R}^{n} : Ax \le b, b \in {}^{\omega}\mathbb{R}^{m}, A \in {}^{\omega}\mathbb{R}^{m \times n}, m, n \in {}^{\omega}\mathbb{N}^{*}\}\) be the feasible domain of the LP max \(\{{c}^{T}x : c \in {}^{\omega}\mathbb{R}^{n}, x \in P\}\). The dual of the latter gives \({x}^{*} \in {}^{\omega}\mathbb{R}_{\ge 0}^{m}\). Putting \(x := {x}^{+} – {x}^{-}\) with \({x}^{+}, {x}^{-} \ge 0\) yields \({x}^{*} \in {}^{\omega}\mathbb{R}_{\ge 0}^{2n}\). Solving max \(\{-h \in {}^{\omega}\mathbb{R}_{\le 0} : Ax – b \le \underline{h}_n\}\) yields \(x \in P_{\ge 0}\) when \(b \ge 0\) does not hold. Let \(|\text{min } \{{b}_{1}, …, {b}_{m}\}|\) be the initial value for \(h\) and 0 its target value. Start with \(x := 0\). Only one pivot step causes that \({b}^{*} \ge 0\). Let \(i, j, k \in {}^{\omega}\mathbb{N}^{*}\) and let \({a}_{i}^{T}\) the \(i\)-th row vector of \(A\).

If \({c}_{j} \le 0\) for all \(j\), the LP is solved. If for some \({c}_{j} > 0\) also \({a}_{ij} \le 0\) for all \(i\), the LP is positively unbounded. If for some \({c}_{j} = 0\) also \({a}_{ij} \le 0\) for all \(i\), drop \({c}_{j}\) and \({A}_{.j}\) as well as \({b}_{i}\) and \({a}_{i}\), but only when \({a}_{ij} < 0\) holds. The inequality \({a}_{ij}{x}_{j} \ge 0 > {b}_{i}\) for all \(j\) has no solution, too. If necessary, divide all \({a}_{i}^{T}x \le {b}_{i}\) by \(||{a}_{i}||\) and all \({c}_{j}\) and \({a}_{ij}\) by the minimum of \(|{a}_{ij}|\) such that \({a}_{ij} \ne 0\) for each \(j\). This will be reversed later. If necessary, renormalise by \(||{a}_{i}||\). Redundant constraints (with \({a}_{i} \le 0\)) may always be removed.

For each \({c}_{j} > 0\) and non-base variable \({x}_{j}\) select min \(\{{b}_{k}/{a}_{kj} : {a}_{ij} > 0\}\) to obtain a next \({x}^{*} \in P_{\ge 0}\) where \({x}_{j}^{*} = {x}_{j} + {b}_{k}/{a}_{kj}\). To select the steepest edge, pick the pivot \({a}_{kj}\) corresponding to \({x}_{j}\) that maximises \({c}^{T}{_1}{\Delta x}\) or \({c}_{j}^{2}/\grave{v}_{j}^{2}\) with \(v_j := ||{A}_{.j}||\) in the \(k\)-th constraint. Multiple maxima allow to use the rule of best pivot value max\({}_{k,j} {c}_{j}{b}_{k}/{a}_{kj}\) or (slower) the smallest angle min \({{_{(1)}}(1, …, 1)}^{T}{_1}{c}^{*}\). If always \(c^Tx^* = c^Tx\) holds, perturb, which means relax, the constraints with \({b}_{i} = 0\) by the same, minimal modulus.

For \({b}_{i} := ||{a}_{i}||\), the modulus need not to be written into the tableau. If another multiple vertex is encountered, despite this being unlikely, simply increase the earlier \({b}_{i}\) by \(||{a}_{i}||\). Leaving a multiple vertex, after which the relaxation is reverted, may require to solve an LP with \(c > 0\) and \(b = 0\). Along the chosen path, the objective function increases otherwise strictly monotonically.

Eventually, \({c}_{j}^{*}, {a}_{ij}^{*}\) and \({b}_{i}^{*}\) can be simply computed using the rectangle rule1cf. Vanderbei, Robert J.: Linear Programming; 3rd Ed.; 2008; Springer; New York, p. 63. Despite the diameter theorem for polytopes, the simplex method is not polynomial under any given set of pivoting rules in the worst-case scenario: An exponential “drift” for e. g. Klee-Minty or Jeroslow polytopes may provide a large deviation from the shortest path by forcing the selection of an unfavourable edge for every step. The result follows in accordance with the state of research.\(\square\)

Theorem: The centre method solves every solvable LP in \(\mathcal{O}(\omega{\vartheta}^2)\).

Proof and algorithm: First, normalise and scale \({b}^{T}y – {c}^{T}x \le 0, Ax \le b\) as well as \({A}^{T}y \ge c\). Let \(d \in [0, 1]\) the density of \(A\) and \(P_r := \{(x, y)^T : x \in {}^{\omega}\mathbb{R}_{\ge 0}^{n}, y \in {}^{\omega}\mathbb{R}_{\ge 0}^{m},{b}^{T}y – {c}^{T}x \le r \in [0, \check{r}], Ax – b \le \underline{r}_m, c – {A}^{T}y \le \underline{r}_n\}\) for the radius \(\check{r} := s|\min \; \{b_1, …, b_m, -c_1, …, -c_n\}|\) and the scaling factor \(s \in [1, 2]\). Putting the number \(z := \grave{m} + n\) implies \(\underline{0}_{\acute{z}} \in \partial P_{\check{r}}\). By the strong duality theorem2loc. cit., p. 60 – 65 the LP min \(\{r \in [0, \check{r}] : (x, y)^T \in P_r\}\) solves the LPs max \(\{{c}^{T}x : c \in {}^{\omega}\mathbb{R}^{n}, x \in {P}_{\ge 0}\}\) and min \(\{{b}^{T}y : y \in {}^{\omega}\mathbb{R}_{\ge 0}^{m}, {A}^{T}y \ge c\}\).

Its solution is the geometric centre \(g\) of the polytope \(P_r\). For \(p_k^* := (\text{min}\,p_k + \text{max}\,p_k)/2\) and \(k = 1, …, z\) approximate \(g\) by \(p_0 := (x_0, y_0, r_0)^T\) until \(||\Delta p||_1\) is sufficiently small. The solution \(t^o(x_1, y_1, r_1)^T\) of the two-dimensional LP min \(\{r \in [0, \check{r}] : t \in {}^{\omega}\mathbb{R}_{> 0}, t(x_0, y_0)^T \in P_r\}\) approximates \(g\) better and achieves \(r \le \check{r}/\sqrt{z}\). Repeat this for \(t^o(x_1, y_1)^T\) until \(g \in P_0\) is computed in \(\mathcal{O}({}_z\check{r} {}_e\check{r}dmn)\) if it exists. Finally, numbers of length \(\mathcal{O}({\omega})\) can only be processed in \(\mathcal{O}(\vartheta).\square\)

Remarks: If the centre method is optimised for distributed computing in \({}^{\nu}\mathbb{R}^{\nu}\), its runtime only amounts to \(\mathcal{O}(1)\). It is also well-suited for (mixed) integer problems and (non-) convex (Pareto) optimisation (according to nature3cf. Vasuki, A: Nature-Inspired Optimization Algorithms; 1st Ed.; 2020; CRC Press; Boca Raton). Rounding errors can be kept small by using a modified Kahan-Babuška-Neumaier summation. Loops may be parallelised. The transfer to complex numbers is easy. All of this holds also for that what follows.

Conclusion: The LP max \(\{{||x – {x}^{o}||}_{1} : {c}^{T}x = {c}^{T}{x}^{o}, Ax \le b, x – {x}^{o} \in {[-1, 1]}^{n}, x \in {}^{\omega}\mathbb{R}_{\ge 0}^{n}\}\) can determine for the first solution \({x}^{o}\) a second one in \(\mathcal{O}(\omega{\vartheta}^2)\) if any, where \({y}^{o}\) may be treated analogously.\(\square\)

Conclusion: The LP max \(\{\nu|\lambda_j| + ||x_j||_1: Ax_j = \lambda_j x_j, |\lambda_j| \in [0, r_j], x_j \in {[-1, 1]}^{n*}\}\) can determine for every \(j = 1, …, n\) the eigenvalue \(\lambda_j \in {}^{\omega}\mathbb{R}\) and the eigenvector \(x_j \in {}^{\omega}\mathbb{R}^{n}\) of the matrix \(A \in {}^{\omega}\mathbb{R}^{n \times n}\) in \(\mathcal{O}(\omega{\vartheta}^2).\square\)

Conclusion: The LP min \(\{r \in [0, s \, \text{max } \{|{b}_{1}|, …, |{b}_{m}|\}] : \pm(Ax – b) \le \underline{r}_m\}\) can determine an \(x \in {}^{\omega}\mathbb{R}^{n}\) of every solvable linear system (LS) \(Ax = b\) in \(\mathcal{O}(\omega{\vartheta}^2)\). The LPs max \(\{{x}_{j} : Ax = 0\}\) yield all solutions to the LS. The matrix \(A\) is regular if and only if the LP max \(\{{||x||}_{1} : Ax = 0\} = 0.\square\)

Conclusion: Let \({\alpha }_{j} := {A}_{.j}^{-1}\) for \(j = 1, …, n\) concerning the matrix \({A}^{-1} \in {}^{\omega}\mathbb{R}^{n \times n}\) and let \({\delta}_{ij}\) the Kronecker delta. A regular \(A\) has an eigenproduct \(\ne 0\) and allows every LS \({A \alpha }_{j} = {({\delta}_{1j}, …, {\delta}_{nj})}^{T}\) to be solved in \(\mathcal{O}(\omega{\vartheta}^2).\square\)

Corollary: Every solvable convex programme min \(\{{f}_{1}(x) : x \in {}^{\omega}\mathbb{R}^{n}, {({f}_{2}(x), …, {f}_{m}(x))}^{T} \le 0\}\) where the \({f}_{i} \in {}^{\omega}\mathbb{R}\) are convex functions for \(i = 1, …, m\) may be solved by the centre method and two-dimensional bisection or Newton’s methods in polynomial runtime, if the number of operands \({x}_{j}\) of the \({f}_{i}\) is \(\le {\omega}^{\nu-3}\) and if an \(x\) exists4cf. Bertsekas, Dimitri P.: Nonlinear Programming; 3rd Ed.; 2016; Athena Scientific; Belmont, p. 589 ff. so that \({f}_{i}(x) < 0\) for all \(i > 1.\square\)

Please, input linear programme (separators spaces, last row
objective function, first column right-hand sides as in the example): not the number of the correct decimal places matters, but the time complexity and the number of steps of the simplex method!

© 2008-2020 by Boris Haase



1 cf. Vanderbei, Robert J.: Linear Programming; 3rd Ed.; 2008; Springer; New York, p. 63
2 loc. cit., p. 60 – 65
3 cf. Vasuki, A: Nature-Inspired Optimization Algorithms; 1st Ed.; 2020; CRC Press; Boca Raton
4 cf. Bertsekas, Dimitri P.: Nonlinear Programming; 3rd Ed.; 2016; Athena Scientific; Belmont, p. 589 ff.