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 rule¹cf. 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 theorem²loc. 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 nature³cf. 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\) exists⁴cf. 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

top