Download Linear Programming Algorithms: Geometric Approach and more Study notes Algorithms and Programming in PDF only on Docsity!
18 Linear Programming Algorithms (December 6)
In this lecture, we’ll see a few algorithms for actually solving linear programming problems. The most famous of these, the simplex method, was proposed by George Dantzig in 1947. Although most variants of the simplex algorithm performs well in practice, there is no sub-exponential upper bound on the running time of any simplex variant. However, if the dimension of the problem is considered a constant, there are several linear programming algorithms that run in linear time. I’ll describe a particularly simple randomized algorithm due to Raimund Seidel. My approach to describing these algorithms will rely much more heavily on geometric intuition than the usual linear-algebraic formalism. This works better for me, but your mileage may vary. For a more traditional description of the simplex algorithm, see Robert Vanderbei’s excellent textbook Linear Programming: Foundations and Extensions [Springer, 2001], which can be freely downloaded (but not legally printed) from the author’s website.
18.1 Bases, Feasibility, and Local Optimality
Consider the linear program max{c · x | Ax ≥ b, x ≥ 0 }, where A is an n × d constraint matrix, b is an n-dimensional coefficient vector, and c is a d-dimensional objective vector. We will interpret this linear program geometrically as looking for the lowest point in a convex polyhedron in IRd, described as the intersection of n + d halfspaces. As in the last lecture, we will consider only non-degenerate linear programs: At most d constraint hyperplanes pass through any point, and no constraint hyperplane is normal to the objective vector. A basis is a subset of d constraints, which by our non-degeneracy assumption must be linearly independent. The location of a basis is the unique point x that satisfies all d constraints with equality; geometrically, x is the unique intersection point of the d hyperplanes. The value of a basis is c · x, where x is the location of the basis. There are precisely
(n+d d
bases. Geometrically, the set of constraint hyperplanes defines a decomposition of IRd^ into convex polyhedra; this cell decomposition is called the arrangement of the hyperplanes. Every d-tuple of hyperplanes (i.e., every basis) defines a vertex of this arrangement (the location of the basis). I will use the words ‘vertex’ and ‘basis’ interchangeably. A basis is feasible if its location x satisfies all the linear constraints, or geometrically, if the point x is a vertex of the polyhedron. A basis is locally optimal if its location x is the optimal solution to the linear program with the same objective function and only the constraints in the basis. Geometrically, a basis is locally optimal if its location x is the lowest point in the intersection of those d halfspaces. A careful reading of the proof of the Strong Duality Theorem reveals that local optimality is the dual equivalent of feasibility; a basis is locally feasible for a linear program Π if and only if the same basis is feasible for the dual linear program q. For this reason, locally optimal bases are sometimes also called dual feasible. Every (d − 1)-tuple of hyperplanes defines a line in IRd^ that is broken into edges by the other hyperplanes. Most of these edges are line segments joining two vertices, but some are infinite rays, which we think of as segments joined to an artificial vertex called ∞. This collection of vertices and edges is called the 1-skeleton of the graph of the hyperplane arrangement. The graph of vertices and edges on the boundary of the feasible region is a subgraph of the arrangement graph. Two bases are said to be adjacent if they have d − 1 constraints in common; geometrically, two vertices are adjacent if they are joined by an edge in the arrangement graph. A pivot changes a basis B into some basis adjacent to B, or equivalently, moves a point from one vertex in the arrangement to an adjacent vertex.
The Weak Duality Theorem implies that the value of every feasible basis is less than or equal to the value of every locally optimal basis; every feasible vertex is higher than every locally optimal vertex. The Strong Duality Theorem implies that (under our non-degeneracy assumption), if a linear program has an optimal solution, it is the unique vertex that is both feasible and locally optimal. Moreover, the optimal solution is both the lowest feasible vertex and the highest locally optimal vertex.
18.2 The Simplex Algorithm: Primal View
From a geometric standpoint, Dantzig’s simplex algorithm is very simple. The input is a set of halfspaces H; we want the lowest vertex in the intersection of these halfspaces.
Simplex(H): if ∩H = ∅ return Infeasible x ← any feasible vertex while x is not locally optimal 〈〈pivot downward, maintaining feasibility 〉〉 x ← any feasible neighbor of x that is lower than x if x = ∞ return Unbounded return x.
Let’s ignore the first three lines for the moment. The algorithm maintains a feasible vertex x. At each pivot operation, it moves to a lower vertex, so the algorithm never visits the same vertex more than once. Thus, after at most
(n+d d
pivots, the algorithm either finds the optimal solution, or it discovers that the polyhedron is unbounded. Notice that we have left open the method for choosing which neighbor to choose at each pivot. There are several natural pivoting rules, but for most rules, there are input polyhedra that require an exponential number of pivots. No pivoting rule is known that guarantees a polynomial number of pivots in the worst case.
18.3 The Simplex Algorithm: Dual View
We can also geometrically interpret the execution of the simplex algorithm on the dual linear program q. Algebraically, there is no difference between these two algorithms; the only change is in how we interpret the linear algebra geometrically. Again, the input is a set of halfspaces H, and we want the lowest vertex in the intersection of these halfspaces. By the Strong Duality Theorem, this is the same as the highest locally-optimal vertex in the hyperplane arrangement.
Simplex(H): if there is no locally optimal vertex return Unbounded x ← any locally optimal vertex while x is not feasbile 〈〈pivot upward, maintaining local optimality 〉〉 x ← any locally-optimal neighbor of x that is higher than x if x = ∞ return Infeasible return x.
Simplex(H): x ← any vertex H˜ ← translation of H that makes x locally optimal while x is not locally optimal 〈〈pivot downward, maintaining feasibility 〉〉 x ← any feasible neighbor of x in Argh( H˜) that is lower than x if x = ∞ return Unbounded while x is not feasible 〈〈pivot upward maintaining local optimality 〉〉 x ← any locally optimal neighbor of x in Argh(H) that is higher than x if x = ∞ return Infeasible return x
18.5 Linear Expected Time for Fixed Dimensions
In most geometric applications of linear programming, the number of variables is a small constant, but the number of constraints may still be very large. The input to the following algorithm is a set H of n halfspaces and a set B of b hyperplanes. (B stands for basis.) The algorithm returns the lowest point in the intersection of the halfspaces in H and the hyperplanes B. At the top level of recursion, B is empty. I implicitly assume here that the linear program is bounded; if necessary, we can guarantee boundedness by adding a single halfspace to H. A point x violates a constraint h if it is not contained in the corresponding halfspace.
SeidelLP(H, B) : if |B| = d x ← ⋂ B if x violates any constraint in H return Infeasible else return x if |H ∪ B| = d return ⋂ (H ∪ B) h ← random element of H x ← SeidelLP(H \ h, B) (∗) if x = Infeasible return x else if x violates h return SeidelLP(H \ h, B ∪ ∂h) else return x
The point x recursively computed in line (∗) is not the optimal solution precisely when the random halfspace h is one of the d halfspaces that define the optimal solution. In other words, the probability of calling SeidelLP(H, B ∪ h) is exactly (d − b)/n. Thus, we have the following recurrence for the expected number of recursive calls for this algorithm:
T (n, b) =
1 if b = d or n + b = d T (n − 1 , b) +
d − b n · T (n − 1 , b + 1) otherwise
The recurrence is somewhat simpler if we let δ = d − b:
T (n, δ) =
1 if δ = 0 or n = δ T (n − 1 , δ) +
δ n · T (n − 1 , δ − 1) otherwise
It’s easy to prove by induction that T (n, δ) = O(δ! n):
T (n, δ) = T (n − 1 , δ) + δ n
· T (n − 1 , δ − 1)
≤ δ! (n − 1) +
δ n (δ − 1)! · (n − 1) [ind. hyp.]
= δ! (n − 1) + δ! n − 1 n ≤ δ! n
At the top level of recursion, we perform one violation test in O(d) time. In each of the base cases, we spend O(d^3 ) time computing the intersection point of d hyperplanes, and in the first base case, we spend O(dn) additional time testing for violations. More careful analysis implies that the algorithm runs in O(d! · n) expected time.
