Homepage of Boris Haase

Previous | Next

#70: Improvement Linear Programming on 27.11.2017

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

Proof and algorithm: Let P := {x ∈ cn : Ax ≤ b, b ∈ cm, A ∈ cm×n, m, n ∈ cℕ*} be the feasible domain of the LP max {dTx : d ∈ cn, x ∈ P}. By taking the dual or setting x := x+ - x- with x+, x- ≥ 0, we obtain x ≥ 0. We first solve max {-z : Ax - b ≤ (z, , z)Tcm, z ≥ 0} to obtain a feasible x when b ≥ 0 does not hold. Initial and target value are z := |min {b1, ..., bm}| resp. z = 0 and we begin with x := 0 as in the first case. Pivoting if necessary, we may assume that b ≥ 0.

Let i, j, k ∈ cℕ* and let aiT the i-th row vector of A. If dj ≤ 0 for all j, the LP is solved. If for some dj > 0 (dj = 0), aij ≤ 0 for all i, the LP is positively unbounded (for now, we may drop dj and A.j as well as bi and ai, but only when aij < 0 holds). Otherwise, we divide all aiTx ≤ bi by ||ai|| and all dj and aij by the minimum of |aij| such that aij ≠ 0 for each j to improve the runtime performance, if possible. This will be reversed later. If necessary, renormalise by ||ai||.

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 dj > 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 dT(x* - x)/||x* - x|| resp. dj2/(1 + ||A.j||2) in the k-th constraint. If there are multiple maxima, select max djbk/akj according to the rule of best pivot value or alternatively (less well) the smallest angle min Σ dj*/||d*||.

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 dj are 0. Along the chosen path, the objective function increases (effectively) strictly monotonically. We can then simply calculate dj*, 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 O(mn), if it is solvable.

Proof and algorithm: We compute the LP min {h ∈ c≥0 : bTy - dTx ≤ h, 0 ≤ x ∈ cn, 0 ≤ y ∈ cm, Ax - b ≤ (h, , h)Tcm, d - ATy ≤ (h, , h)Tcn} via the (dual) programme min {bTy : 0 ≤ y ∈ cm, ATy ≥ d} for the (primal) programme max {dTx : d ∈ cn, x ∈ P≥0}. First, we normalise bTy - dTx ≤ 0, Ax ≤ b and ATy ≥ d. We compute all scalar products dTai resp. -dj and sort them. Let v := (xT, yT)Tcm+n. We start with v := 0.

We compute x and analogously y from some equations aiTx = bi resp. xj = 0 with the largest corresponding dTai resp. -dj. Initial value for the height h is the minimum needed for the inequalities to be valid. If x and y are optimal, we stop. Otherwise, we successively define a few times all vk* := (max vk + min vk)/2. We solve the two-dimensional LPs min hk in hk, vkc≥0. We determine a favourable direction u ∈ cm+n with uk := Δvk Δhk/min Δhk.

Thereafter, we solve the two-dimensional LP min h in h, t ∈ c≥0 with v* := v + tu. We always continue with the v* that belongs to the smallest h resp. hk reached so far. We solve the two-dimensional LPs via bisection method. We may basically extrapolate h und v in O(m + n). When we cannot leave v, we relax all constraints except v ≥ 0 a bit and undo that at the end of one iteration step in O(mn).

The diagonal of a cube with dimension O(max (m, n)) representing h, which has the target value 0, does not matter in the calculation. Since we process in each iteration step a fix amount of h, the claim follows by the strong duality theorem ([775], p. 60 - 65). The source code is not published here, since we can solve (average) LPs faster by the simplex- or the face-algorithm ([916], p. 580 f.).⃞

Corollary: Every linear system Ax = b for x ∈ cn can be solved in O(mn), if we write it as LP min {h ∈ c≥0 : ±(Ax - b) ≤ (h, , h)Tcm} and if it is solvable.⃞

Corollary: An eigenvector x ∈ cn \ {0} of the matrix A for the eigenvalue λ ∈ cℝ can be determined in O(mn) with the LP min {h ∈ c≥0 : ±(Ax - λx) ≤ (h, , h)Tcm}, if this is solvable.⃞

Corollary: Every convex programme min {f1(x) : x ∈ cn, (f2(x), , fm(x))T ≤ 0} where the ficℝ are convex functions is strongly polynomial and may be solved by the normal method and two-dimensional Newton's methods, which handle the fi, in O(p), where p ∈ cℕ* denotes the number of operands xj of the fi, assuming that it is solvable and that an x exists so that fi(x) < 0 for all i > 1 (see [939], p. 589 ff.).⃞

Remarks: The normal method is numerically very stable, since the initial data are barely altered, and currently the fastest known worst-case LP-solving algorithm. It may also be applied to (mixed) integer problems and branch-and-bound, in particular for nonconvex optimisation. It can be better paralleled than the simplex method. It can easily be extended to convex programmes with a vector-valued objective function f1.

Diameter theorem for polytopes: The diameter of an n-dimensional polytope defined by m constraints with m, n ∈ cℕ* 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 d are integers wlog and that m and n are fixed. By dualising, a full-dimensional LP may be obtained as described before within the normal method. This shows that the problem of (mixed) integer linear programming is not NP-complete:

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

code of simplex method

© 27.11.2017 by Boris Haase

Valid XHTML 1.0 • privacy policy • disclaimer • pdf-version • bibliography • subjects • definitions • statistics • php-code • rss-feed • mwiki • top