Linear Programming

In the following section, Set Theory and Nonstandard Analysis are presupposed. The simplex and the centre method can solve a LP mostly in cubic and quadratic runtime respectively.

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: If the finite stint is dropped, the theorem can analogously be transferred to polyhedra.

Definition: Let \({|| \cdot ||}_{1}\) the 1-norm. Passing \(x \in X \subseteq {}^{\omega}\mathbb{R}^{\omega}\) to the next step, lets it follow by \({}^{*}\) where \(\Delta x := {x}^{*} – x\). The unit vector of \(x \in X^{*}\) is \({_1}{x} := x/||x||\). 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.\(\triangle\)

Theorem: The simplex method¹Dantzig, George B.: Lineare Programmierung und Erweiterungen; 1. Aufl.; 1966; Springer; Berlin. is exponential in the worst case.

Proof and algorithm: Let the LP max \(\{{c}^{T}x : c \in {}^{\omega}\mathbb{R}^{n}, x \in P\}\) have the feasible domain \(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}^{*}\}\). Its dual 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}\) if \(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. Dividing 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\) maybe repeated and reversed later. 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)}}\underline{1}^{T}_n}{_1}{c}^{*}\).

If \(c^Tx^* = c^Tx\) is repeated, perturb, which means relax, the constraints with \({b}_{i} = 0\) by \(d \in {}^{\omega}\mathbb{R}_{>0}\). For \({b}_{i} := ||{a}_{i}||\), drop \(d\) in the tableau. If another multiple vertex is encountered for once, simply add \(||{a}_{i}||\) to every \({b}_{i}\). Leaving it, after which the relaxation is reverted, may require to solve an LP with \(c > 0\) and \(b = 0\).

This task can also finally arise, if further solutions of the LPs are to be determined for at least two \(c_j = 0\). This avoids having to compute several minimal relaxations that differ strongly but are not selected according to optimality points of view, as with the lexicographical method. Bland’s rule is not optimal, too.

The rectangle rule simply allows to compute \({c}_{j}^{*}, {a}_{ij}^{*}\) and \({b}_{i}^{*}\)²cf. Vanderbei, Robert J.: Linear Programming; 3rd Ed.; 2008; Springer; New York, p. 63.. Despite the diameter theorem for polytopes, no pivoting rule can prevent the simplex method from being exponential as the case may be: An exponential “drift” for e. g. Klee-Minty or Jeroslow polytopes may always provide an unfavourable edge. The result follows in accordance with the state of research.\(\square\)

Remarks: The following method has been developed in 35 years of intensive research since the meaning and scope of LPs were early clear. It could turn out to be one of the most important algorithms ever, because despite various (numerical) difficulties that it presented, it is speed that counts for great problems. The simplex method, on the other hand, allows simple and precise calculations with rational numbers if the initial problems are appropriately formulated.

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

Proof and algorithm: Let \(z := \grave{m} + n\) and \(d \in [0, 1]\) the density of \(A\). First, normalise and scale \({b}^{T}y – {c}^{T}x \le 0, Ax \le b\) as well as \({A}^{T}y \ge c\). Let \(P_r := \{(x, y)^T \in {}^{\omega}\mathbb{R}_{\ge 0}^{z} : {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\}\) have the radius \(\check{r} := s|\min \; \{b_1, …, b_m, -c_1, …, -c_n\}|\) and the scaling factor \(s \in [1, 2]\). It follows \(\underline{0}_{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_0\). For \(p_k^* := (\text{min}\,p_k + \text{max}\,p_k)/2\) and \(k = 1, …, \grave{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^o, y^o, r^o)^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{\grave{z}}\). Repeat this for \(t^o(x^o, y^o)^T\) until \(g \in P_0\) is computed in \(\mathcal{O}({}_z\check{r} {}_e\check{r}dmn)\) if it exists. Numbers of length \(\mathcal{O}({\omega})\) can only be processed in \(\mathcal{O}(\vartheta)\) as is generally known.

Solving all two-dimensional LPs \(\text{min}_k r_k\) by bisection methods for \(r_k \in {}^{\omega}\mathbb{R}_{\ge 0}\) and \(k = 1, …, z\) in \(\mathcal{O}({\vartheta}^2)\) each time determines \(q \in {}^{\omega}\mathbb{R}^k\) where \(q_k := \Delta p_k \Delta r_k/r\) and \(r := \text{min}_k \Delta r_k\). Let simplified \(|\Delta p_1| = … = |\Delta p_{z}|\). Here min \(r_z\) for \(p^* := p + wq\) and \(w \in {}^{\omega}\mathbb{R}_{\ge 0}\) would be also to solve. If \(\text{min}_k \Delta r_k r = 0\) follows, stop, otherwise repeat until min \(r = 0\) or min \(r > 0\) is sure. If necessary, constraints are temporarily relaxed by the same small modulus.\(\square\)

Conclusion: If neither a primally feasible \(x\) nor a dual solution \(y\) needs to be computed, the runtime of the LP max \(\{{c}^{T}x : c \in {}^{\omega}\mathbb{R}^{n}, x \in {P}_{\ge 0}\}\) can roughly be halved by setting \(r := {c}^{T}x\).

Remarks: Simplex method and face algorithm⁴Pan, Ping-Qi: Linear Programming Computation; 1st Ed.; 2014; Springer; New York., p. 580 f. may solve the LP faster for small \(m\) and \(n\). The current stock of constraints or variables can easily be changed because the centre method does not durably transform constraints and is faster than all known (worst-case) LP-algorithms in \(\mathcal{O}(\omega^{37/18}\vartheta)\). Details are first published when no misuse for non-transparent or bad decisions must be feared.

Remarks: The centre method makes use of the fact that the calculated geometric centre and the centre of the sphere coincide in \(P_0\). The former can be secured by rotating the coordinate axes (in pairs) with the aid of a simple rotation matrix in a comparable running time. For this purpose, their number may have to be increased to an even number. It is advantageous that the calculations can of course also be carried out independently of an optimisation question.

Remarks: If the centre method is optimised for distributed computing in \({}^{\nu}\mathbb{R}^{\nu}\), its runtime only amounts to \(\mathcal{O}(1)\). If some \(p_k\) are replaced by \(\lfloor p_k \rfloor\) or \(\lceil p_k \rceil\) after solving the continuous problem, it proves successful for (mixed) integer problems and (non-) convex (Pareto) optimisation (according to nature⁵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. 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)\) as well as matrix multiplication to be executed in \(\mathcal{O}(\omega{\vartheta})\) by parallelising.\(\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 so that \({f}_{i}(x) < 0\) for all \(i > 1\)⁶cf. Bertsekas, Dimitri P.: Nonlinear Programming; 3rd Ed.; 2016; Athena Scientific; Belmont., p. 589 ff..\(\square\)

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

© 2008-2021 by Boris Haase

top