Linear Programming

In the following section, Set Theory, Topology and Nonstandard Analysis are presupposed.

Diameter theorem for polytopes: Every 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 all remaining dimensions. Overcoming every minimal distance requires at most two edges and yields the factor 2.\(\square\)

Definition: Let \(\vartheta := \omega\,{}_{\epsilon}\omega\) and \({|| \cdot ||}_{1}\) the 1-norm. If \(x \in X \subseteq {}^{\omega}\mathbb{R}^{\omega}\) passes to the next step, it is followed by \({}^{*}\) where \(\Delta x := {x}^{*} – x\). The unit vector of \(x \in X^{*}\) is \({_1}{x} := x/||x||\). 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}({\epsilon}^{|\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¹Dantzig, George B.: Lineare Programmierung und Erweiterungen; 1. Aufl.; 1966; Springer; Berlin. may solve an LP almost always in \(\mathcal{O}({\omega}^3\vartheta)\).

Proof and algorithm: Let 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}\). A (big) \(x^+ \in {}^{\omega}\mathbb{R}_{\geq0}\) achieves \(x^* \in {}^{\omega}\mathbb{R}_{\geq0}^n\) for \(x := x^* – {x^+}{\upharpoonright}_n\). Solving max \(\{-h \in {}^{\omega}\mathbb{R}_{\le 0} : Ax – b \le h{\upharpoonright}_m\}\) 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 \({b_k}/{a}_{kj} :=\) min \(\{{b_i}/{a}_{ij} : {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}\) for \(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{\upharpoonright}^{T}_n}\;{_1}{c}^{*}\).

If \(c^Tx^* = c^Tx\) is repeated, relax all 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 may require to solve an LP with \(c > 0\) and \(b = 0\). Then the relaxation (perturbation) is reverted. 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 as the lexicographical method does. 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.. Bland’s rule is not optimal, too. 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\)

Bit-square theorem (bit\({}^2\) theorem for short): Repetitions allow more than \(\mathcal{O}(n^2)\) multiplications and additions to be reduced to \(\mathcal{O}(n^2)\) for an algorithm of numbers of bit length \(m \le n\) with \(m, n \in {}^{\omega}\mathbb{N}^{*}.\square\)

Examples: The matrix product \((c_{ij}) := \left({\LARGE{\textbf{+}}}_{k=1}^n{a_{ik}b_{kj}}\right)\) of \((a_{ij}), (b_{ij}) \in {}^{\omega}\mathbb{C}^{n \times n}\) over products \({\LARGE{\textbf{+}}}_{k=1}^n{b_{jk}4^{km}}\) with \(a_{ij}\) and \(i, j \in {}^{\omega}\mathbb{N}_{\le n}^{*}\), determining an \(n\)-dimensional Taylor series (see Nonstandard Analysis), the Gauss-Jordan algorithm, the simplex and the Gram-Schmidt orthogonalisation method.

Remark: If \(m\) is too large, it will be appropriately decreased computationally. Coding quaternions enables concerted matrix multiplication in the cloud by comparing simultaneously.

Theorem: Intex methods solve almost every solvable LP with inter-/extrapolations in \(\mathcal{O}({\vartheta}^3)\).

Proof and algorithm: Let \(z := 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, \rho], Ax – b \le r{\upharpoonright}_m,\) \(c – {A}^{T}y \le r{\upharpoonright}_n\}\) have the radius \(\rho := s|\min \; \{b_1, …, b_m, -c_1, …, -c_n\}|\) and the scaling factor \(s \in [1, 2]\).

It follows \(0{\upharpoonright}_{z} \in \partial P_{\rho}\). By the strong duality theorem³loc. cit., p. 60 – 65., LP min \(\{ r \in [0, \rho] : (x, y)^T \in P_r\}\) solves 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}\,\check{p}_k + \text{max}\,\check{p}_k\) 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, \rho] : t \in {}^{\omega}\mathbb{R}_{> 0}, t(x_0, y_0)^T \in P_r\}\) approximates \(g\) better and achieves \(r \le \check{\rho}\).

Repeat this for \(t^o(x^o, y^o)^T\) until \(g \in P_0\) is computed in \(\mathcal{O}({}_2\rho^2dmn)\) if it exists.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_{\grave{z}}\) for \(p^* := p + wq\) and \(w \in {}^{\omega}\mathbb{R}_{\ge 0}\) may be solved, too. If \(\text{min}_k \Delta r_k r = 0\) follows, stop computing, otherwise repeat until min \(r = 0\) or min \(r > 0\) is sure.\(\square\)

Remark: Details are only published when no misuse for non-transparent or bad decisions must be feared. Simplex method and face algorithm⁴Pan, Ping-Qi: Linear Programming Computation; 1st Ed.; 2014; Springer; New York., p. 580 f. solve LPs faster for small \(m\) and \(n\). Intex methods are faster than all known (worst-case) LP-algorithms in \(\mathcal{O}(\omega^{37/18}\vartheta)\). Rounding errors provoked in \(b\) and \(c\) make any solution unique.

Remark: If necessary, the same small modulus relaxes constraints temporarily. The runtime for \(m, n \in {}^{\nu}\mathbb{N}^*\) is \(\mathcal{O}(dmn)\). Optimising it for distributed computing in \({}^{\nu}\mathbb{R}^{\nu}\) reduces it to \(\mathcal{O}(1)\). Integer solutions \(v_j\) satisfying \(\acute{v}_0 \le \lfloor v_j \rfloor – v_j\) prove them 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.).

Remark: Intex methods make use of the fact that the calculated geometric centre and the centre of the sphere coincide in \(P_0\). The former is secured by rotating the coordinate axes (in pairs) like \(c\) with the aid of a simple rotation matrix in a comparable runtime. Numbers of length \(\mathcal{O}({\omega})\) can only be processed in \(\mathcal{O}(\vartheta)\) as is generally known. The following results can also be easily transferred to complex numbers:

Conclusion: 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}\) any second one in \(\mathcal{O}({\vartheta}^3)\) where \({y}^{o}\) may be treated analogously.\(\square\)

Conclusion: 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 matrix \(A \in {}^{\omega}\mathbb{R}^{n \times n}\) in \(\mathcal{O}({\vartheta}^3).\square\)

Conclusion: LPs max \(\{{x}_{j} : Ax = 0\}\) yield all solutions to \(Ax = b\) in \(\mathcal{O}({\vartheta}^3)\). Matrix \(A\) is regular if and only if LP max \(\{{||x||}_{1} : Ax = 0\} = 0.\square\)

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

Conclusion: The polynomial \(c_{\acute{n}}x^{\acute{n}} + … + c_1x + c_0\) can be determined in \(\mathcal{O}({\vartheta}^3)\) by measurement \((x_j, y_j) \in {}^{\omega}\mathbb{R}^2\) and residues \(r_j \in {}^{\omega}\mathbb{R}\) for \(j = 1, …, m\) and \((c, x)^T \in {}^{\omega}\mathbb{R}^{\grave{n}}\) from \(y = r + Xc\) and \(X^Tr = 0\) where \(X = \left(x_j^k\right)\) for \(k = 0, …, \acute{n}\) (also under additional constraints⁶cf. loc. cit., p. 196 ff..\(\square\)

Conclusion: The DFT (see Scientific Computing) can determine the Padé approximant \({a(x)}/{b(x)}={{\LARGE{\textbf{+}}}_{m=0}^n {a_m x^m}}/{{\LARGE{\textbf{+}}}_{m=0}^n {b_m x^m}}\) of \(f(x)\) for \(x \in {}^{\omega}\mathbb{R}\) in \(\mathcal{O}({\vartheta}^3).\square\)

Corollary: Intex methods and two-dimensional bisection or Newton’s methods may solve most solvable convex programmes min \(\{{f}_{1}(x) : x \in {}^{\omega}\mathbb{R}^{n}, {({f}_{2}(x), …, {f}_{m}(x))}^{T} \le 0\}\) where all \({f}_{k} \in {}^{\omega}\mathbb{R}\) are convex functions for \(k = 1, …, m\) in polynomial runtime, if the number of operands \({x}_{j}\) of all \({f}_{k}\) is \(\le {\omega}^{\nu-3}\) and if an \(x\) exists such that \({f}_{k}(x) < 0\) for all \(k > 1\)⁷cf. Geiger, Carl; Kanzow, Christian: Theorie und Numerik restringierter Optimierungsaufgaben; 1. Aufl.; 2002; Springer; Berlin, p. 56 ff..\(\square\)

Corollary: LP max \(\{x : y_{\acute{m}} – xy_m = b_{\acute{m}}, y_n = y_0 = 0, x \le x_0 \in {}^{\omega}\mathbb{R}\}\) can determine for \(m = 1, …, n\) by Horner scheme an \(x \in {}^{\omega}\mathbb{R}\) solving \(n\)-polynomial \(b_{\acute{n}}x^{\acute{n}} + … + b_1x + b_0 = 0\) in \(\mathcal{O}({\vartheta}^2)\). Aborting TS can analogously solve every continuous inequation system with convex solution set.\(\square\)

© 2008-2021 by Boris Haase

top