## #67: Improvement Linear Programming on 04.10.2016

In the following section, we solve linear programmes (LPs). The well-known simplex method is described, together with the maximally quadratic and furthermore strongly polynomial normal method as a solution of Smale's 9th problem. We may generalise the normal method to (non-) convex programmes (with vector-valued functions). The diameter theorem for polytopes is proven. It is shown that (mixed) integer LPs are polynomial.

Proof and algorithm: Let M := {x ∈ κn : Ax ≤ b, b ∈ κm, A ∈ κm×n, m, n ∈ κℕ*} be the feasible domain of the LP max {cTx : c ∈ κn, x ∈ M}. By taking the dual or setting x := x+ - x- with x+, x- ≥ 0, we obtain x ≥ 0. We first solve max {-z : Ax - ze ≤ b, z ∈ κ≥0} with objective z = 0 and e = (1, ..., 1)Tκm to obtain a feasible x when b ≥ 0 does not hold. We begin with z := |min {b1, ..., bm}| and x := 0 as in the first case. Pivoting if necessary, we may assume that b ≥ 0.

Let i, j, k ∈ κℕ* and let aiT the i-th row vector of A. If cj ≤ 0 for all j, the LP is solved. If for some cj > 0, aij ≤ 0 for all i, the LP is positively unbounded. Otherwise, we divide all aiTx ≤ bi by ||ai|| and all cj and aij by the minimum of |aij| such that aij ≠ 0 for each j. This will be reversed later. If necessary, renormalise by ||ai||. This yields good runtime performance even on strongly deformed polytopes.

In each step, we can remove multiple constraints and such with ai ≤ 0, since they are redundant (avoidable by adding an extra slack variable in each case). The second case is analogous. If in both cases bi = 0 and ai ≥ 0 for some i, then the LP has maximum 0 and solution x = 0 if b ≥ 0, otherwise it has no solutions. In each step, for each cj > 0 and non-base variable xj, we select the minimum ratio bk/akj for aij > 0.

The variables with * are considered in the next step. The next potential vertex is given by xj* = xj + bk/akj for feasible x*. To select the steepest edge, select the pivot akj corresponding to xj that maximises cT(x* - x)/||x* - x|| i.e. cj2/(1 + Σ aij2) in the k-th constraint. If there are multiple maxima, select max cjbk/akj or alternatively the smallest angle min Σ cj*/||c*||, according to the rule of best pivot value.

If there are more than n values bi equal 0 and we cannot directly maximise the objective function, we relax (perturb) the constraints with bi = 0 by the same, minimal modulus. These do not need to be written into the tableau: We simply set bi = ||ai||. If another multiple vertex is encountered, despite this being unlikely, simply increase the earlier bi by ||ai||.

The cost of eliminating a multiple vertex, after which we revert the relaxation, corresponds to the cost of solving an LP with b = 0. The same task is potentially required at the end of the process when checking whether the LP admits any other solutions if at least two of the cj are 0. Along the chosen path, the objective function increases (effectively) strictly monotonically. We can then simply calculate cj*, aij* and bi* using the rectangle rule.

In the worst-case scenario, the simplex method is not strongly polynomial despite the diameter formula for polytopes (see below) under any given set of pivoting rules, since an exponential "drift" can be constructed with Klee-Minty or Jeroslow polytopes, or others, creating a large deviation from the shortest path by forcing the selection of the least favourable edge. This is consistent with existing proofs. The result follows.⃞

Theorem: The normal method solves the LP in at most O(mn) and is strongly polynomial.

Proof and algorithm: The normal method includes the two phases of the simplex method and solves the LP if possible by moving in a direction calculated from one of O(n) orthogonal directions inside M towards the potential maximum. We solve every two-dimensional LP in O(m) by using the bisection method. The method begins identically to the simplex method, the reduction to one phase in O(mn) via the dual programme, too.

We define the height h1 := cTx and wlog eliminate x1 with c1 > 0. If we cannot leave x, we set bi* := bi + d for all aiTx = bi with a small d ∈ κ>0. We begin each iteration step by solving the LPs max hj in hj, u ∈ κ≥0 for all j ≥ 1 with x* := (hj, x2, ..., xj-1, u, xj+1, ..., xn)T and constant xk for 1 ≠ k ≠ j. We then calculate a favourable direction r ∈ κn with rj := xjΔhj/max Δhj. We solve the LP max h1 in h1, v ∈ κ≥0 with x* := (h1, x2 + vr2, ..., xn + vrn)T.

We always continue with the x* belonging to the greatest value hj reached so far for j ≥ 1. We successively define xj* := (max xj + min xj)/2 for all j > 1, until we have changed these xj* sufficiently, if possible. Then we set bi* := bi again. We generally may extrapolate from preceding results in O(m + n). Thus, we have finished one iteration step in O(mn). Since we process all hj independently from m and n in O(1), the claim follows.⃞

Remarks: We terminate the procedure if we cannot increase h1 anymore or if the LP appears to be unlimited, since we reach a ceiling for h1. If b ≥ 0 does not hold, we first solve the LP max h1 subject to Ax + q ≤ b and x ≥ 0 with q ∈ κm, z := |min {b1, ..., bm}|, together with qi := h1 - z for bi < 0 and qi := 0 otherwise. The initial value is h1 = 0 and the target value is h1 = z. The dual programme to max {cTx : c ∈ κn, x ∈ M, x ≥ 0} reads min {bTy : y ∈ κm, y ≥ 0, ATy ≥ c}.

Corollary: Every linear system (LS) of equations Ax = b with x ∈ κn may be solved in at most O(mn), provided that a solution exists.

Proof: We write Ax = b as Ax ≤ b and -Ax ≤ -b where x = x+ - x- with x+, x- ≥ 0.⃞

Theorem: Every convex programme min {f1(x) : x ∈ κn, x ≥ 0, (f2(x), …, fm(x))T ≤ 0} where the fiκℝ are convex posynomials is strongly polynomial and may be solved by the normal method and Newton's method, which handles the fi, in at most O(p), assuming that it is solvable, where p ∈ κℕ* denotes the number of operands xj of the fi and the objective function f1 is linearised.

Proof: The claim follows from the existence of the normal method.⃞

Remarks: The normal method is currently the fastest known LS/LP-solving algorithm and is numerically very stable, since the initial data are barely altered. It can also be applied to branch and bound, in particular for nonconvex optimisation. The normal method is a better candidate for parallelisation than the simplex method. It can easily be extended to convex programmes with vector-valued or other convex fi.

Diameter theorem for polytopes: The diameter of a n-dimensional polytope defined by m constraints with m, n ∈ κℕ* is at most max (2(m - n), 0).

Proof: Each vertex of a (potentially deformed) hypercube is formed by at most n hyperplanes. If we complete the polytope by adding or removing hyperplanes, the claim follows for the chosen path, since each step adds at most two additional edges. This theorem can be extended to polyhedra analogously by dropping the requirement of finiteness.⃞

Further thoughts: Gomory or Lenstra cuts can find an integer solution of the original problem in polynomial time if we additionally assume that a, b, and c are integers wlog and that m and n are fixed. By finding any redundant constraints and hyperplanes, a full-dimensional LP may be obtained. This shows that the problem of (mixed) integer linear programming is not NP-complete:

Theorem: (Mixed) integer LPs may be solved in polynomial time.⃞