A subexponential bound for linear programming.
We present a simple randomized algorithm which solves linear programs with n constraints and d variables in expected O(nde(d ln(n+1))1/4) time in the unit cost model (where we count the number of arithmetic operations on the numbers in the input). The expectation is over the internal randomizations performed by the algorithm, and holds for any input. The algorithm is presented in an abstract framework, which facilitates its application to several other related problems. The algorithm has been presented in a previous work by the authors [ShW], but its analysis and the subexponential complexity bound are new.
Introduction
Linear programming is one of the basic problems in combinatorial optimization, and as
such has received considerable attention in the last four decades. Many algorithms have
been proposed for its solution, starting with the simplex algorithm and its relatives [Dan],
proceeding through the polynomial-time solutions of Khachiyan [Kha] and Karmarkar
[Kar], and continuing with several more recent techniques (reviewed below). While some
of the proposed algorithms have proven out to be extremely efficient in practice, analysis
of their running time has not been fully satisfactory so far. For example, the simplex algorithm
was shown to be exponential in the worst case [KlM]. The algorithms of Khachiyan
and of Karmarkar have polynomial bit-complexity, but the number of arithmetic operations
they perform depends on the size of the coefficients defining the input and it cannot
be bounded solely in terms of n (the number of constraints) and d (the number of vari-
ables). In this paper we assume a different model of computation, namely the real RAM,
widely used in computational geometry. Here the input may contain arbitrary real num-
bers, and each arithmetic operation with real numbers is charged unit cost. To distinguish
complexity bounds in this model from bounds in the bit-complexity model, we call them
combinatorial.
Until recently, the best known combinatorial bounds were exponential in either n or d;
a subexponential (randomized) bound is given in a very recent paper of Kalai [Kal] and also
in this paper. One of the major open problems in the area is to find a strongly polynomial
algorithm (i.e. of combinatorial polynomial complexity) for linear programming.
In this paper we describe a simple randomized algorithm which solves linear programs
with n inequalities and d variables with expected
d ln(n=
d )+O(
d+ln n) g)
arithmetic operations. In conjunction with Clarkson's linear programming algorithm
[Cla3], this gives an expected bound of
O(d
d
The expectation of the running time is with respect to the internal randomizations performed
by the algorithm, and holds for any input. This complexity matches that of a
recent algorithm due to Kalai [Kal]. (Except for the constant in the exponent, the only
significant difference is that Kalai has a version of his algorithm which runs in e O(
d) as
long as n is linear in d. We can guarantee this bound only for
d).) Chronologically
speaking, our algorithm was published first in [ShW], but with a weaker analysis of its
running time; Kalai's analysis with a subexponential bound came next, and immediately
afterwards we realized that our algorithm also has subexponential running time.
The algorithm is presented in an abstract framework which facilitates the application of
the algorithm to a large class of problems, including the computation of smallest enclosing
balls (or ellipsoids) of finite point sets in d-space, computing largest balls (ellipsoids) in
convex polytopes in d-space, computing the distance between polytopes in d-space, etc.
(however, we can guarantee a subexponential running time for only some of these problems;
see below for details).
To compare the complexity of our algorithm with other recent techniques, here is a brief
survey of some relevant recent literature. Megiddo [Meg2] has given the first deterministic
algorithm whose running time is of the form O(C d n), and is thus linear in n for any fixed
d. However the factor C d in his algorithm is 2 2 d
; an improvement to C
can be found
in [Dye1] and [Cla1]. Recently, a number of randomized algorithms have been presented
for the problem, see [DyF], [Cla3], [Sei], with a better dependence on d, where the best
expected running time is given by Clarkson [Cla3]: O(d log n). Nevertheless,
the dependence on d is still exponential. The more recent algorithm of Seidel [Sei] has
worse expected complexity of O(d!n), but is an extremely simple randomized incremental
algorithm. In [Wel] this algorithm was enhanced with a "move-to-front" heuristic, which
in practice has drastically improved the performance of the algorithm, but was (and still
is) very difficult to analyze. Our algorithm is another variant, in-between the techniques
of [Sei] and [Wel]. (It is interesting, that there are examples of linear programs (with
few constraints), where adding the move-to-front heuristic to our new algorithm gives
a significantly worse performance [Mat1].) Our algorithm also seems to behave well in
practice, and its analysis, as given here, also provides a considerably improved worst-case
upper bound on its expected complexity. Recently, derandomization techniques have been
applied to Clarkson's algorithm to obtain deterministic O(C d n) time LP algorithms with
C d of the order 2 O(d log d) [ChM].
The abstract framework we present considers the set H of n constraints, and a function
w which maps every subset of H to its optimal solution, where w satisfies a few simple
conditions; as it turns out, this is all which is needed to prove the correctness of the
algorithm, and to analyze its expected running time in terms of two primitive operations
'violation tests' and `basis computations'. It turns out that these primitive operations
can easily be implemented in polynomial time for linear programming, but that is by no
means clear for other instances of problems in the framework. For example in the case
of computing the smallest enclosing ball, a basis computation amounts to computing the
smallest ball enclosing d points in d-space. Only recently, G-artner [G-ar2] was able to
show that this can be done with expected e O(
d) arithmetic operations.
Clarkson's algorithm [Cla3] can also be shown to work in this framework, while the
algorithm of Seidel [Sei] (and the generalization to other applications in [Wel]) needs to
make a more explicit use of the geometry of the problems. A different framework has been
recently developed by Dyer [Dye2], which yields deterministic algorithms linear in n (with
larger constants in d).
The paper is organized as follows. In the next section we introduce some notations and
review basic observations on linear programming that are required for the presentation and
analysis of our algorithm in Section 3. The analysis culminates in a recurrence relationship,
whose solution is a nontrivial and interesting task in itself; Section 4 is devoted to that
solution. Finally, in Section 5 we describe the abstract framework, and mention a few
examples.
Notation and Basics
In the following we will prepare the notation required for our presentation. We model the
case of linear programming where the objective function is defined by the vertical vector
and the problem is to maximize c \Delta x over an intersection of given n
halfspaces.
We let '!' be the lexicographical ordering on R d , (i.e. x
1 , or x
and
2 , or so on). This ordering is extended to R d [ f+1; \Gamma1g ( +1 and \Gamma1
being special symbols) by the convention that
Let H be a set of closed halfspaces in d-space (which we call also constraints in d-space).
By v H we denote the lexicographically smallest point in the feasible region FH :=
v H is called the value of H . (It seems more standard to look for the backwards lexicographically
smallest point, i.e., the case when the reader accustomed
to this may prefer to think 'backwards' instead.) If FH is empty, then we let v
If FH is nonempty, but contains no minimum, then we let v
Intuitively, we can view +1 as a point which lies in all halfspaces and which dominates
all points. \Gamma1 may be seen as a point in minus infinity, or simply as a symbol for
'undefined'.
A basis B is a set of constraints with
of B. If a basis of H is a minimal subset B of H with
Constraint h is violated by H , if then this is equivalent to
Constraint h is extreme in H , if v H \Gammafhg ! v H . Note that h is extreme in H if
and only if h 2 H and h is violated by H \Gamma fhg.
The following properties are either trivial (like (i) and (ii)), or standard in linear
programming theory (like (iii)). We cite them explicitly here, because they constitute all
that is needed for the correctness and time analysis of the algorithms to be described.
H be sets of constraints in d-space, G ' H, and let h be a constraint
in d-space.
violated by G if and only if h is violated by H.
(iii) If v H ! +1, then any basis B ' H has exactly d constraints, and any F ' H has at
most d extreme constraints.
3 The Algorithm
Let H+ be the set of constraints x Given a set H of constraints, the
following algorithm will compute a basis of H [ H+ ; such a basis always exists, because
(The algorithm works for any initial basis B 0 instead of
in particular, we could take any basis B 0 ' H , if available, or we can take a set of
constraints d, for a symbolic value \Gamma! representing an arbitrarily small
We will always assume that v H[H+ ! +1 for the constraint sets we consider. It
is well-known that this condition can be obtained at the cost of an extra dimension and an
extra constraint. However, note that we have not made any general position assumptions;
so, e.g., we do not require that bases are unique.
Given a set H of n constraints, we might first consider the following trivial approach.
Remove a random constraint h, and compute a basis B for H \Gamma fhg recursively. If h is
not violated by B (or equivalently, if h is not violated by H \Gamma fhg), then B is a basis for
H and we are done. If h is violated, then we try again by removing a (hopefully different)
random constraint. Note that the probability that h is violated by H \Gamma fhg is at most d
since there are at most d extreme constraints in H .
So much for the basic idea, which we can already find in Seidel's LP algorithm [Sei];
however, in order to guarantee efficiency, we need some additional ingredients. First,
the procedure SUBEX lp actually solving the problem has two parameters: the set H of
constraints, and a basis C ' H ; in general, C is not a basis of H , and we call it a candidate
basis. It can be viewed as some auxiliary information we get for the computation of the
solution. Note that C has no influence on the output of the procedure (but it influences
its efficiency).
The problem of computing a basis of H [ H+ can now be solved by:
function procedure subex lp(H); /* H: set of n constraints in d-space;
returns a basis of H [ H+ .
return
The following pseudocode for SUBEX lp uses two primitive operations: The first is a test
'h is violated by B', for a constraint h and a basis B (violation test); this operation
can be implemented in time O(d), if we keep v B with the basis B. Second, SUBEX lp
assumes the availability of an operation, 'basis(B; h)', which computes a basis of B [fhg
for a d-element basis B and a constraint h (basis computation). This step corresponds to
a pivot step, and with an appropriate representation of B, it can be performed with O(d 2 )
operations. Note that if
all constraint sets H considered in an execution.
function procedure SUBEX lp(H; C); /* H: set of n constraints in d-space;
return C /* returns a basis of H.
else
choose a random h
if h is violated by B then /* , v B 62 h
return SUBEX lp(H; basis(B; h))
else
return
A simple inductive argument shows that the procedure returns the required answer.
This happens after a finite number of steps, since the first recursive call decreases the
number of constraints, while the second call increases the value of the candidate basis
(and there are only finitely many different bases).
For the analysis of the expected behavior of the algorithm, let us take a closer look at
the probability that the algorithm makes the second recursive call with candidate basis
'basis(B; h)'. As noted above, this happens exactly when h is extreme in h. Since we now
choose h from H \Gamma C (and C always has d elements), it follows that the probability for
being extreme is at most d
n\Gammad . Moreover, if d \Gamma k extreme constraints in H are already
in C, then the bound improves to k
n\Gammad ; in fact, there are never more 'bad' choices than
there are choices at all, and so the bound can be lowered to minfn\Gammad;kg
n\Gammad . We want to show
that the numerator decreases rapidly as we go down the recursion, and this will entail the
subexponential time bound.
The key notion in our analysis will be that of hidden dimension. For given G ' H
and a basis C ' G, the hidden dimension of the pair (G; C) measures how "close" is C
to the solution (basis) of G. The hidden dimension of (G; C) is d minus the number of
constraints h 2 C which must be contained in any basis B ' G with value greater than
or equal to the value of C.
Lemma 2 The hidden dimension of (G; C) is d \Gamma jfh 2 G; v G\Gammafhg
Proof. We need to show that a constraint h 2 G satisfies v G\Gammafhg appears in
G is a basis with
h is not in B, then B Conversely, if
be a basis for G \Gamma fhg. We have h is not
in B. This completes the proof of the lemma.
Let us remark that the hidden dimension of (G; C) depends on C only via the value
of C. An intuitive interpretation is that all the "local optima" defined by constraints of
G lying above the "level" v C are contained in a k-flat defined by some d \Gamma k bounding
hyperplanes of constraints of C, k being the hidden dimension. Hidden dimension zero
implies that C is a basis of G.
We enumerate the constraints in H in such a way that
ng
The ordering is not unique, but the parameter k emerging from this ordering is unique
and, by definition, it equals the hidden dimension of (H; C).
Note that for h are in C. Hence
the only choices for h which possibly entail a second call are h
is chosen, and
the first call (with candidate basis C) returns with basis B, . Then we
compute the hidden dimension of
the pair (H; C 0 ) for the second call is at most k \Gamma i.
The hidden dimension is monotone, i.e., if C ' G ' H , then the hidden dimension of
(G; C) does not exceed the hidden dimension of (H; C). This holds, because the constraints
are in C (and so in G), and v G\Gammafh i g - v H \Gammafh i g because G ' H .
We denote by t k (n) and b k (n) the maximum expected number of violation tests and
basis computations, respectively, entailed by a call SUBEX lp(H; C) with n constraints and
hidden dimension - k. The discussion above yields the following recurrence relationships:
Simple proofs by induction show that b k (n) - 2 k (n \Gamma d) and t k (n) - 2 k (n \Gamma d). It turns
out that for n not much larger than d, this is a gross overestimate, and in the next section
we will show that
d ln(n=
d )+O(
We have also every constraint is tested at most
once for violation by B, for each basis B ever appearing in the computation; we have tha
d) since we do not test the elements in B for violation by B, and we add 1 to
k (n) to account for the initial basis H+ .
Recall that we account O(d) arithmetic operations for a violation test, and O(d 2 ) arithmetic
operations for a basis computation. Note also that for the computation of v H[H+ ,
we have to add the d nonnegativity constraints before we initiate SUBEX lp. Anticipating
the solution of the recurrences, given in the next section, we thus conclude:
Theorem 3 Let H be a set of n constraints in d-space. The lexicographically smallest
nonnegative point v H[H+ in the intersection of the halfspaces in H can be computed by a
randomized algorithm with an expected number of
O(nd
d ln(n=
d )+O(
d+ln n)
operations.
Clarkson [Cla3] shows that a linear program with n constraints in d-space can be solved
with expected
O(d
stands for the complexity of solving a linear program
with 9d 2 constraints in d-space. We replace T d (9d 2 ) by our bound, and observe that, after
this replacement, the middle term in (3.2) is always dominated by the two other terms;
moreover, we can even omit the 'log n' factor in the last term in (3.3) without changing
the asymptotics of the whole expression. Thus, we obtain:
Corollary 4 The lexicographically smallest nonnegative point v H[H+ in the intersection
of the n halfspaces H in d-space can be computed by a randomized algorithm with an
expected number of
O(d
d log d) ) (3.3)
operations.
4 Solving the Recurrence
In this section we prove the promised bounds for b k (n). We first put - b k
and we let f(k; n) be the function which satisfies the recurrence for - b k (n) with equality
(thus majorizing - b k (n)).
Thus f(k; n) is the solution to the recurrence
and
We prove the following upper bound:
Lemma 5 For all n; k, f(k; n) holds. If n - e k=4
k, then
/s
O
ip
For a large range of values of n; k, this bound is essentially tight (at least the leading
term). In order to avoid various technicalities in the lower bound proof, we restrict the
range of n and we will not try to get as tight a bound as we could (concerning the order
of magnitude of the lower order terms). We show
Lemma 6 For k tending to infinity and n such that k - n - 2 o(k) , we have
s
We emphasize that this is a lower bound on f(k; n) and not on - b k (n). We do not have
such a lower bound on - b k (n), and it is conceivable that - b k (n) is much smaller. See [Mat1]
for recent related lower bounds.
Upper bound. For the proof of Lemma 5 we apply the technique of generating functions
(the easy inductive proof of the '1+2 k n'-bound has been mentioned in the previous section,
and will be omitted here). Define, for each n - 0,
If we multiply the recurrence (4.1) by z k and sum over k, we obtain
G
G
z jX
G
In other words, we have the recurrence
, it follows that
Y
Regarding z as a complex variable, we want to find the value of the coefficient f(k; n) of
z k in the Taylor expansion of G n (z). As is well known, this is equal to
f(k; n) =2-i
Z
z k+1 dz
where the integral is over a closed curve fl that contains the origin and does not contain
any other pole of G n (z).
We will choose for fl the circle 1. It is easily checked that none of the
denominators in the product defining G n (z) can vanish for this follows from the
inequality
This also implies that
which yields a bound of
Y
where
be an integer parameter to be determined later, and let us set
We estimate the terms in the product for q ? 2. We distinguish two cases. First, for j - q
we use the estimate
so, using Stirling's formula,
Y
q, the absolute values of the terms in the alternating sum in the parentheses
form a decreasing sequence. Hence if we stop the summation after some term, the error
will have the same sign as the first omitted term. So we have
(j
Then
q!
Finally we estimate
so we get
For most values of n; k, we set
which yields the bound
/s
O
ip
There are two extreme cases. If the above definition of q gives q ! 2 (this happens
when n is exponentially large compared to k, n ? e k=4
k), we use the easy bound of
comes close to
k or it is smaller, we set
and obtain the bound e O(
. (By playing with the estimates somewhat further, one can
get better bounds from (4.2) both for n !!
k (in this case we observe that that if n ! q,
the product in (4.3) only goes up to n, so that one gets C n instead of C q , and then a
better choice of q is possible), and for n very large, that is, n ?? e k=4
k.) This establishes
Lemma 5.
Lower bound. The proof of the lower bound is based on the following "explicit form"
for f(k; n):
Lemma 7
where the first sum is over all q-tuples (m
and the second sum is over all q-tuples (i
(Here q can also be 0;
this contributes one term equal to 1 to the sum. We can formally interpret this as a term
corresponding to the (unique) 0-tuple; this term, as the empty product, is equal to 1.)
Proof. As for the initial conditions, for only have the term with in the
sum, yielding the value 1. Similarly for only get a nonzero term for
again gives 1.
Consider the difference f(k; n) \Gamma f(k; n) is defined by (4.5). The
terms which appear in (4.5) for k; n and not for k; are precisely those with i
For the sum of such terms, we have
mq=1n
=n
The function defined by (4.5) thus satisfies both the initial conditions and the recurrence
relation.
For the proof of Lemma 6 let us first define several auxiliary parameters. Let
a function tending to 0 slowly enough, for instance log log k. Set
s
"q
Since we assume that k is large, we may neglect the effect of truncating q and m to
integers; similarly for other integer parameters in the sequel. The assumption log
guarantees that
Consider the inner sum in (4.5) for a given q-tuple (m
k. Then this sum is at least
One easily verifies (by induction, say) that 1
we use Stirling's formula to get the estimate so we can estimate the
above expression from below by
It is well-known (and easy to prove) that the number N(r; k) of r-tuples (m
with
r
. Let N   (r; k; m) denote the number of
r-tuples as above with the additional condition we claim that if q -
then
To see this it suffices to encode every r-tuple (m contributing to
a q-tuple contributing to N   (q; k; m). Such an encoding can be performed as follows: if
we leave this element as it is, otherwise we replace it by bm
equal to m and one element equal to m j +1\Gamma(m\Gamma1)bm
with the tuple (1; 5; 2; 8) will be transformed to (1; 3; 3; 2; 2; 3; 3; 3; 3; 1); note that
replaced by a block of t elements, their sum is t). The number of elements
of the resulting vector is
q, and we may
then pad the vector with ones to get exactly q elements. It is easy to see that this encoding
is reversible. By what has just been observed, the sum of elements in the resulting vector
will exceed the sum of the initial r-tuple by at most q. This proves (4.7).
With our value of m, we may choose use (4.7) to show that the
number of q-tuples (m k is at least
Combining this with (4.6) and observing that ln(n=m) -
k), we obtain
which proves Lemma 6.
5 An

Abstract

Framework
Let us consider optimization problems specified by pairs (H; w), where H is a finite set,
is a function with values in a linearly ordered set (W ; -); we assume
that W has a minimum value \Gamma1. The elements of H are called constraints, and for
is called the value of G. The goal is to compute a minimal subset BH of H
with the same value as H (from which, in general, the value of H is easy to determine),
assuming the availability of two basic operations to be specified below. It turns out that
the algorithm in Section 3 can be used to perform this computational task, as long as the
following axioms are
Axiom 1. (Monotonicity) For any F; G with F ' G ' H , we have
Axiom 2. (Locality) For any F ' G ' H with
implies that also w(F
If Axioms 1 and 2 hold, then we call (H; w) an LP-type problem. Obviously, linear programming
is an LP-type problem, if we set a constraint set G ' H . Then
the axioms coincide with Lemma 1 (i) and (ii). The notions needed in Section 3 carry over
in the obvious way: A basis B is a set of constraints with
for all proper subsets B 0 of B. For G ' H , if a basis of G is a minimal
subset B of G with Constraint h is violated by G, if w(G)
Constraint h is extreme in G, if w(G \Gamma fhg) ! w(G).
For the efficiency of the algorithm the following parameter is crucial: The maximum
cardinality of any basis is called the combinatorial dimension of (H; w), and is denoted by
dim(H; w).
We assume that the following primitive operations are available.
(Violation test) 'h is violated by B', for a constraint h and a basis B, tests
whether h is violated by B or not.
(Basis computation) 'basis(B; h)', for a constraint h and a basis B, computes
a basis of B [ fhg.
(Note that checking whether h violates B is equivalent to checking whether w(basis(B; h)) ?
w(B). This shows that the two primitive operations are closely related.) Now we have
all the ingredients necessary to apply our algorithm to an LP-type problem, provided we
have an initial basis B 0 for the first call of SUBEX lp. We can also show, using the simpler
inductive argument mentioned in Section 3, that the expected number of primitive operations
performed by the algorithm is O(2 ffi n), where is the
combinatorial dimension. However, in order to ensure the subexponential time bound we
need an extra condition:
Axiom 3. (Basis regularity) If B is a basis with and h is a
constraint, then every basis of B [ fhg has exactly dim(H; w) elements.
If Axioms 1-3 are satisfied, then we call (H; w) a basis-regular LP-type problem. We have
seen that linear programming in d-space is a basis-regular LP-type problem of combinatorial
dimension d, provided the program is feasible (Lemma 1 (iii) provides the stronger
property that every basis has the same cardinality). Actually, in the whole treatment
in Section 3 we were careful not to use other properties of linear programming except
for those formulated in Lemma 1. (We would be very glad to use any other properties
to obtain a faster algorithm, but we do not know how.) Of course, we have an extra
computational assumption in order to start the algorithm.
(Initial basis) An initial basis B 0 with exactly dim(H; w) elements is available.
(In the case of linear programming, the nonnegativity-constraints H+ can play the role of
the initial basis.) We may conclude
Theorem 8 Given a basis-regular LP-type problem (H; w), jH
dimension ffi and an initial basis B the algorithm SUBEX lp
computes a basis of H with an expected number of at most
ffi+ln n)
violation tests and basis computations.
As it turns out, Clarkson's algorithm can also be formulated and analysed within the
framework, and, if the basic cases (each involving O(ffi 2 ) constraints) are solved by our
algorithm, then the expected number of required violation tests is O(ffin
log ffi) ),
and the expected number of required basis computations is e O(
log ffi) log n.
Many problems can be shown to satisfy Axioms 1 and 2 (see below for a list of such
problems), but - except for linear programming - basis-regularity is not naturally satisfied
by any of them. However, we can artificially enforce Axiom 3 by the following
trick (due to Bernd G-artner, [G-ar1]): Let (H; w) be an LP-type problem of combinatorial
dimension ffi and with value set W . Then we define a pair (H; w 0 ), where
for G ' H , and the new value set is f\Gamma1g [
ordering. The straightforward proof of the following lemma is omitted.
Lemma 9 Given an LP-type problem (H; w), the pair (H; w 0 ) as defined above is a basis-
regular LP-type problem of combinatorial dimension dim(H; w). Any basis B complemented
by arbitrary dim(H; w) \Gamma jBj elements can serve as initial basis.
Hence, we can transform every LP-type problem into a basis-regular one, although we have
to be careful about the new interpretation of violation tests and basis computations. We
thus obtain an algorithm with an expected subexponential number of violation tests and
basis computations, but those primitive operations might be quite expensive. We exhibit
now two examples of LP-type problems where we can successfully apply our algorithm
(including efficient implementation of the primitive operations).
Smallest enclosing ball. The problem of computing the disk of smallest radius containing
a given set of n points in the plane goes back to J. J. Sylvester in 1857 [Syl].
The first linear time solution to this problem was provided by Megiddo [Meg1] (see this
reference for a short guide through previous O(n 3 ) and O(n log n) solutions). The general
problem of computing the smallest ball enclosing a set of n points in d-space can be solved
in linear time as long as the dimension is considered fixed, see [Meg1, Meg2] for a deterministic
algorithm, and [Wel] for a simple randomized algorithm; however both algorithms
are exponential in d.
Given a set P of n points in d-space, we define r(Q) as the smallest radius of a closed
ball containing Q ' P . It is well-known that this smallest radius exists, and that the ball
realizing this radius is unique. Moreover, there is always a subset B of Q containing at
most e.g. [Jun]). With these basic facts in hand,
it is easy to show that (P; r) is an LP-type problem of combinatorial dimension
Clearly, adding points (constraints) to a set cannot decrease the radius of the smallest
enclosing ball (and so monotonicity holds), and p is violated by Q ' P if and only of p lies
outside the unique smallest ball containing Q (this easily implies locality). The problem is
not basis-regular, so we have to apply the above transformation to validate our analysis.
While the violation test is an easy 'point in ball'-test, the basis computation amounts to
a non-trivial task: basically we have to compute the smallest ball enclosing d+ 2 points in
d-space. Recently, G-artner [G-ar2] was able to show that this problem can be solved with
expected e O(
d) arithmetic operations. Hence, we obtain
Corollary 10 The smallest enclosing ball of n points in d-space can be computed with an
expected number of
d) n; e 2
d ln(n=
d )+O(
operations.
Combining this with Clarkson's algorithm, the complexity reduces to the bound in (3.3)
for linear programming.
Distance between polytopes. Given two closed polytopes P 1 and P 2 , we want to
compute their distance, i.e.
with dist(a; b) denoting the Euclidean distance for points a and b. If the polytopes in-
tersect, then this distance is 0. If they do not intersect, then this distance equals the
maximum distance between two parallel hyperplanes separating the polytopes; such a
of hyperplanes is unique, and they are orthogonal to the segment connecting
two points a 2 P 1 and b 2 P 2 with It is now an easy exercise
to prove that there are always sets B 1 and B 2 of vertices of P 1 and P 2 , respectively, such
that denotes the convex hull of a
the distance is positive, a bound of d
We formulate the problem as an LP-type problem as follows: Let
i is the vertex set of P i , we assume that so every subset U of V has
a unique representation as 2. With this convention,
we define for empty, then we
define 1. The pair (V; w) constitutes now an LP-type problem, except that the
inequalities go the other way round: U ' W ' V implies that w(U) - w(W ), and 1
plays the role of \Gamma1. In order to see locality, observe that p is violated by U , if and only
lies between the unique pair of parallel hyperplanes which separate U 1 and U 2 and
have distance w(U ); this also shows how we can perform a violation test. From what we
mentioned above, the combinatorial dimension of this problem is at most d+2 (or d+ 1, if
the polytopes do not intersect). Hence, for a basis computation, we have to compute the
distance between two polytopes in d space with altogether d+3 vertices. Again, we invoke
[G-ar2] to ensure that this can be performed with expected e O(
d) arithmetic operations.
Corollary 11 The distance between two convex polytopes in d-space with altogether n
vertices can be computed with an expected number of
d) n; e 2
d ln(n=
d )+O(
operations.
The best previously published result was an expected O(n bd=2c ) randomized algorithm
(d considered fixed) in [Cla2]. The result in Corollary 11 can also be established if the
polytopes are given as intersections of n halfspaces, and again combining the result with
Clarkson's yields the bound in (3.3).
Other examples. There is quite a number of other examples which fit into the frame-
work, and thus can be solved in time linear in the number of constraints (the dimension
considered fixed). As we have mentioned before, the subexponential bound is a delicate
issue, which depends how efficiently we can solve small problems. We just provide a list
of examples, without giving further details.
(Smallest enclosing ellipsoid) Given n points in d-space, compute the ellipsoid
of smallest volume containing the points (combinatorial dimension d(d
see [DLL, Juh, Wel]).
The problem has been treated in a number of recent papers, [Pos, Wel, Dye2, ChM]. Here
the constraints are the points, and the value of a set of points is the volume of the smallest
enclosing ellipsoid (in its affine hull). Such an ellipsoid (also called L-owner-John ellipsoid)
is known to be unique, [DLL], from which locality easily follows; monotonicity is obviously
satisfied. The primitive operations can be treated by applying general methods for solving
systems of polynomial inequalities, but we cannot claim any subexponential time bounds
(of course, the bound linear in the number of points holds).
(Largest ellipsoid in polytope) Given a polytope in d-space as the intersection
of n halfspaces, compute the ellipsoid of largest volume contained in the
polytope (combinatorial dimension d(d
(Smallest intersecting ball) Given n closed convex objects in d-space, compute
the ball of smallest radius that intersects all of them (combinatorial dimension
In order to see the combinatorial dimension, we make the following observation. We
consider the Minkowski sum of each convex object with a closed ball of radius r (centered
at the origin). Then there is a ball of radius r intersecting all objects, if and only if
the intersection of all of these Minkowski sums is nonempty. By Helly's Theorem, this
intersection is nonempty, if and only if the intersection of any d + 1 of them is nonempty.
If r is the smallest radius which makes this intersection nonempty, then the interiors of
the Minkowski sums do not have a common point, and so there is a set B 0 of at most d+1
of them which do not have a common point. The corresponding set of objects contains
a basis (which is of cardinality at most d + 1), and the claimed combinatorial dimension
now follows easily.
(Angle-optimal placement of point in polygon) Let P be a star-shaped polygon
with n vertices in the plane (a polygon is star-shaped if there is a point inside
the polygon which sees all edges and vertices; the locus of all such points is
called the kernel). We want to compute a point p in the kernel of P , such that
after connecting p to all the vertices of P by straight edges, the minimal angle
between two adjacent edges is maximized (combinatorial dimension 3).
Unlike in the previous examples, it might not be obvious what the constraints are in this
problem. Let us assume that a polygon allows a placement of p with all angles occurring
at least ff. This restricts the locus of p to the intersection of the following convex regions:
(i) for every vertex v of the polygon, and its incident edges e 0 and e with inner (with
respect to P ) angle fi, there is a wedge of angle fi \Gamma 2ff with apex v, where p is forced
to lie. (ii) For every edge e of P with incident vertices v and v 0 , there is circular arc,
which contains all points which see vertices v and v 0 at an angle ff, and which lie on the
inner side of e; point p is forced to lie in the region bounded by e and this circular arc.
This suggests that we take the vertices (with incident edges) and the edges (with incident
vertices) as constraints (for the purpose of the algorithm ignoring that they stem from a
polygon). Thus the number of constraints is 2n and the combinatorial dimension is 3 by
a reasoning using Helly's theorem as before.
(Integer linear programming) Given n halfspaces and a vector c in d-space, compute
the point x with integer coordinates in the intersection of the halfspaces
which maximizes c \Delta x (combinatorial dimension 2 d , see [d'Oi, Bel, Sca]).
Although integer linear programming fits into the framework, it is a bad example in the
sense that the basis computation has no bound in the unit cost model. There are some
other problems which we did not mention so far, but may occur to the reader as natural
examples, e.g.,
(Largest ball in polytope) Given a polytope in d-space as the intersection of
halfspaces, compute the ball of largest radius contained in the polytope.
(Smallest volume annulus) Given n points in d-space, find two concentric
balls such that their symmetric difference contains the points and has minimal
volume.
Indeed, these problems are LP-type problems, but actually they can be directly formulated
as linear programs with d+ 1 and d+ 2, respectively, variables (the transformation for the
smallest volume annulus-problem can be found in [D-or]). Thus also the subexponential
time bounds hold.
Recently, Chazelle and Matou-sek, [ChM], gave a deterministic algorithm solving LP-
type problems in time O(ffi O(ffi) n), provided an additional axiom holds (together with an
additional computational assumption). Still these extra requirements are satisfied in many
natural LP-type problems. Matou-sek, [Mat2], investigates the problem of finding the best
solution satisfying all but k of the given constraints, for abstract LP-type problems as
defined here.
Nina Amenta, [Ame1], considers the following extension of the abstract framework:
Suppose we are given a family of LP-type problems (H; w - ), parameterized by a real
parameter -; the underlying ordered value set W has a maximum element +1 representing
infeasibility. The goal is to find the smallest - for which (H; w - ) is feasible, i.e. w -
+1. [Ame1] provides conditions under which the such a problem can be transformed into
a single LP-type problem, and she gives bounds on the resulting combinatorial dimension.
Related work can be found in [Ame2].
6 Discussion
We have presented a randomized subexponential algorithm which solves linear programs
and other related problems. Clearly, the challenging open problem is to improve on the
bounds provided in [Kal] and here, and to find a polynomial combinatorial algorithm for
linear programming.
In Section 4 we have shown that the bound we give is tight for the recursion (3.1)
we derived from the analysis. Even stronger, [Mat1] gives abstract examples of LP-type
problems of combinatorial dimension d and with 2d constraints, for which our algorithm
d) primitive operations. That is, in order to show a better bound of our
algorithm for linear programming, we have to use properties other than Axioms 1 through
3.
Rote, [Rot], and Megiddo, [Meg3], suggest dual 'one-permutation' variants of the al-
gorithm. It is interesting, that there are examples of linear programs, where a 'one-
permutation' variant of the algorithm suggested in [Mat1] seems to behave much worse
on certain linear programs than the original algorithm (this fact is substantiated only by
experimental results, [Mat1]); this has to be seen in contrast to the situation for Seidel's
linear programming algorithm, [Wel].

Acknowledgments

. The authors would like to thank Gil Kalai for providing a draft
copy of his paper, and Nina Amenta, Boris Aronov, Bernard Chazelle, Ken Clarkson,
Bernd G-artner, Nimrod Megiddo and G-unter Rote for several comments and useful discussions




--R


"Proc. 10th Annual ACM Symposium on Computational Geometry"
A theorem concerning the integer lattice
"Proc. 4th SIAM-ACM Symposium on Discrete Algo- rithms"
Linear programming in O(n
New applications of random sampling in computational geometry
"Proc. 29th IEEE Symposium on Foundations of Computer Science"

Linear Programming and Extensions
Algorithmen zur Ermittlung der Formab- weichung mit Koordinatenme-ger-aten
On a multidimensional search technique and its application to the Euclidean one-center problem
"Proc. 8th Annual ACM Symposium on Computational Geometry"
A randomized algorithm for fixed-dimensional linear pro- gramming

"Proc. 33rd IEEE Symposium on Foundations of Computer Science"


ACM Symposium on Theory of Computing"
A new polynomial-time algorithm for linear programming
Polynomial algorithm in linear programming
How good is the simplex algorithm
Lower bounds for a subexponential optimization algorithm
"Proc. 10th Annual ACM Symposium on Computational Geometry"
time algorithms for linear time programming in R 3 and related problems
Linear programming in linear time when the dimension is fixed
A note on subexponential simplex algorithms

"Proc. 16th Annual ACM Symposium on Theory of Computing"

"Proc. of the National Academy of Sciences of the United States of America"
Low dimensional linear programming and convex hulls made easy
"Proc. 9th Symposium on Theoretical Aspects of Computer Science"
A question in the geometry of situation
"New Results and New Trends in Computer Science"
--TR
A new polynomial-time algorithm for linear programming
Linear programming in O(n MYAMPERSANDtimes; 3<supscrpt>d<supscrpt>2</></>) time
On a multidimensional search technique and its application to the Euclidean one centre problem
Small-dimensional linear programming and convex hulls made easy
An optimal convex hull algorithm and new results on cuttings (extended abstract)
A subexponential randomized simplex algorithm (extended abstract)
A class of convex programs with applications to computational geometry
Linear Programming in Linear Time When the Dimension Is Fixed
A Combinatorial Bound for Linear Programming and Related Problems
Minimum spanning ellipsoids

--CTR
Bernd Grtner , Sven Schnherr, Exact primitives for smallest enclosing ellipses, Proceedings of the thirteenth annual symposium on Computational geometry, p.430-432, June 04-06, 1997, Nice, France
Kong , Qiuming Zhu, Incremental procedures for partitioning highly intermixed multi-class datasets into hyper-spherical and hyper-ellipsoidal clusters, Data & Knowledge Engineering, v.63 n.2, p.457-477, November, 2007
Michael Goldwasser, A survey of linear programming in randomized subexponential time, ACM SIGACT News, v.26 n.2, p.96-104, June 1995
Matthew J. Katz , Franck Nielsen, On piercing sets of objects, Proceedings of the twelfth annual symposium on Computational geometry, p.113-121, May 24-26, 1996, Philadelphia, Pennsylvania, United States
Nina Amenta, Bounded boxes, Hausdorff distance, and a new proof of an interesting Helly-type theorem, Proceedings of the tenth annual symposium on Computational geometry, p.340-347, June 06-08, 1994, Stony Brook, New York, United States
Nina Amenta, Helly theorems and generalized linear programming, Proceedings of the ninth annual symposium on Computational geometry, p.63-72, May 18-21, 1993, San Diego, California, United States
Bernard Chazelle , Ji Matouek, On linear-time deterministic algorithms for optimization problems in fixed dimension, Proceedings of the fourth annual ACM-SIAM Symposium on Discrete algorithms, p.281-290, January 25-27, 1993, Austin, Texas, United States
Piyush Kumar , Joseph S. B. Mitchell , E. Alper Yildirim, Approximate minimum enclosing balls in high dimensions using core-sets, Journal of Experimental Algorithmics (JEA), 8,
Gil Kalai, A subexponential randomized simplex algorithm (extended abstract), Proceedings of the twenty-fourth annual ACM symposium on Theory of computing, p.475-482, May 04-06, 1992, Victoria, British Columbia, Canada
Kenneth L. Clarkson, Las Vegas algorithms for linear and integer programming when the dimension is small, Journal of the ACM (JACM), v.42 n.2, p.488-499, March 1995
Sandeep Sen, Parallel multidimensional search using approximation algorithms: with applications to linear-programming and related problems, Proceedings of the eighth annual ACM symposium on Parallel algorithms and architectures, p.251-260, June 24-26, 1996, Padua, Italy
Henrik Bjrklund , Sergei Vorobyov, A combinatorial strongly subexponential strategy improvement algorithm for mean payoff games, Discrete Applied Mathematics, v.155 n.2, p.210-229, January, 2007
Ji Matouek, On geometric optimization with few violated constraints, Proceedings of the tenth annual symposium on Computational geometry, p.312-321, June 06-08, 1994, Stony Brook, New York, United States
Henrik Bjrklund , Sergei Vorobyov, Combinatorial structure and randomized subexponential algorithms for infinite games, Theoretical Computer Science, v.349 n.3, p.347-360, December 2005
Michael T. Goodrich, Efficient piecewise-linear function approximation using the uniform metric: (preliminary version), Proceedings of the tenth annual symposium on Computational geometry, p.322-331, June 06-08, 1994, Stony Brook, New York, United States
Michael T. Goodrich, Fixed-dimensional parallel linear programming via relative &egr;-approximations, Proceedings of the seventh annual ACM-SIAM symposium on Discrete algorithms, p.132-141, January 28-30, 1996, Atlanta, Georgia, United States
Michael Saks, Kleitman and combinatorics, Discrete Mathematics, v.257 n.2-3, p.225-247, 28 November
Bernard Chazelle, Computational geometry: a retrospective, Proceedings of the twenty-sixth annual ACM symposium on Theory of computing, p.75-94, May 23-25, 1994, Montreal, Quebec, Canada
randomized incremental algorithms;linear programming;combinatorial optimization;computational geometry