CIS677: Notes for lecture 8: Convexity and Algorithms

CIS677
Notes for lecture 8: Convexity and Algorithms


Outline

Before going back into the mathematics of faces and polarity, let's consider a couple simple algorithms from the point of view of duality.

A dual view of Jarvis march

In the algorithm of Jarvis for planar convex hulls, a line is seen to "rotate" about the sites; at each moment the line rotates about a site, until another site is hit by the line.

The rotating lines correspond to supporting halfspaces of the sites.

The polars of these halfspaces are therefore all in the polar set of the convex hull of the sites. The polars S* of the sites are halfspaces, and the polar of conv(S) is the intersection P(S*) of halfspaces.

In the polar setting, the polar of the "rotating line" is a point that moves along a line bounding the polar halfspace of some site. The motion begins at a vertex of the halfspace intersection.

The process of finding the site hit by the rotating line in "primal space" corresponds to finding the first line hit by the moving point in "dual space".

A dual view of incremental algorithms

In an incremental convex hull algorithm, the convex hull of a subset R of S is maintained, and sites are added one by one to R.

Dually, the intersection P(R*) is maintained, and halfspaces are added one by one from S*. Each additional halfspace s* clips away a portion of the previous intersection. The new vertices and edge in the dual space correspond to the new edges and vertex in the primal space.

Notice that the "clipped off" portion is itself a polygon, the intersection of the complement of s* with P(R*).

Vertices and faces

In Lecture 7, we found that a point v is a vertex iff v is not in conv(S\{v}). From this it follows that:

Lemma. If P=conv(S) is a polytope, then
vert(P) is a subset of S, and P=conv(vert(P)).
Proof. Suppose v is a vertex of P; then conv(S\{v}) is not equal to conv(S), and so v must be in S.

So vert(P) is a subset of S, and so P=conv(S) contains conv(vert(P)).

If x is in S but not in vert(P), then conv(S)=conv(S\{x}), and so by induction on the size of S\vert(P), conv(S)=conv(vert(P)). QED

Lemma. If P is a polyhedron, the intersection of two faces of P is a face of P.

Proof. Suppose f and g are faces, so there are supporting halfspaces

hf := {x | a·x <= b}, and
hg := {x | a'·x <= b'}
whose bounding planes intersect P at f and g. Clearly the halfspace
h := {x | (a+a')·x <= b+b'},
has P a subset of h, and for any x in P,
(a+a')·x = b+b'
only when both
a·x = b and a'·x = b';
hence the intersection of f and g is the intersection of P with the bounding plane of h, and so is a face. QED
Lemma. Let P be a polyhedron. If g is in faces(P), and f is in faces(g), then f is in faces(P).
Proof. Suppose
hg := {x | ag·x <= bg}
witnesses that g is a face of P: that is, hg is supporting and g is the intersection of hg with P. Suppose
hf := {x | af·x <= bf}
is a witness that f is a face of g.

Notice that hf does not necessarily contain P. Now for value z, consider

hz := {x | (ag + zaf)x <= bg + zbf}.
That is, h0=hg. Notice that f is in the bounding plane of hz, for any z.

Suppose P=conv(S1)+cone(S0). Now pick z>0, but small enough that all S1 and S0 are included in hz. Then P is in hz, but the intersection of P with the bounding plane of hz is f. Hence f is a face of P. QED

Corollary. Let P be a polyhedron. If g is in faces(P), then
faces(g)={f in faces(P) | f is a subset of g}.

Proof. The previous lemma implies that the faces of g are a subset of the faces of P.

Suppose f is a face of P that lies in g. If h witnesses that f is a face of P, then h also witnesses that f is a face of g. QED

In particular, the vertices of f in vert(P) are a subset of the vertices of P, and so if P=conv(S), then vert(f) is a subset S, and f=conv(f intersect S).

Faces and lattices

We've seen that for a polyhedron P and a,b,c in faces(P):

P is a face of P; if a is a face of b is a face of c, then a is a face of c.
and easily:
if a is a face of b and b is a face of a, then a=b.

That is, the faces can be partially ordered under the "is a face of" relation. The fact that the intersection of two faces is a face implies also that the faces form a mathematical structure called a lattice.

Faces of Cones

When v is a vertex of a polyhedron P, cone({v}) is an extreme ray of C:=cone(lift(P)); that is, cone({v}) is a 1-dimensional face of C. The set of extreme rays of a cone C is denoted extr(C), and it can be shown that C=cone(extr(C)), which is analogous to P=conv(vert(P)) for a polytope P. (Actually, this requires that C have a vertex, which not all cones have.)

Faces, Cones, and Polarity

We've seen that the polar of a cone is cone. There is a fundamental relation between the faces of a cone C and its polar C*.

Suppose f is a face of C, and h={x | h*·x <= 0} is a witness to this. Since h is supporting, h* is a member of C*. Let fo be the subset of C* comprising all h* for which h is a witness for f. That is,

fo := {y | y·x <0 for all x in C\f, y·x = 0 for all x in f}.
Moreover,
fo = {y | y·x <0 for all x in extr(C)\extr(f), y·x = 0 for all x in extr(f)},
and so fo is a face of C*: it is the intersection of C* with the a set of planes associated with the vertices of f.

From the description of fo, it's immediate that

f is a face of g iff go is a face of fo.

Is fo=f*?

What is fo if C is given as an intersection of halfspaces? We suppose for now that f is determined by the bounding planes of some subset Hf of H.

Lemma. If C=P(H), and f is a face of C given by
f = {x | h*·x <0 for h in H\Hf, h*·x = 0 for h in Hf},
then fo=cone(Hf*).

Proof. Any member y of Hf* has

y·x = 0 for all x in f, and
y·x < 0 for all x in C\f,
and any positive combinations of such y's preserves these conditions, so cone(Hf*) is a subset of fo.

On the other hand, if y is in fo, then y is in cone(H*), so y=sumh* in H* ch h*, with all ch>=0. If ch>0 for h* not in Hf*, then y·x < 0 for x in f, which cannot happen for y in fo. Hence y is in cone(H<fsup>*). QED

Finally:

Lemma. If C is a cone and f is a face of C, then foo=f.

Proof. If C=cone(S), and f= cone(Sf), then

fo = {y | y·h* <0 for all h in (S\Sf)*, y·h* = 0 for all h in Sf*},
and foo = cone(Sf**)= cone(Sf)=f, using the previous lemma. QED

Thus, the mapping f-->fo is a 1-1 correspondence between the faces of C and the faces of C*, and reverses the "is a face of" relation.

A similar mapping holds between the faces of a polyhedron and the faces of its polar.

From this relation, we can get some additional facts about polyhedra, for example, that a face is the intersection of the facets containing it.

Complexity of the convex hull

For now, we'll only state the following:

Upper Bound Theorem. A polytope with n vertices in d dimensions has Theta(nfloor(d/2)) facets.

That is, for d=2 or 3, floor(d/2)=1 so the hull has linear complexity. In fact, using Euler's relation, a polytope in three dimension whose n vertices are in general position has exactly 2n-4 facets.

In four or five dimensions, a polyhedron with n vertices can have as many as Theta(n2) facets.

Algorithms for convex hulls

If a set S of points is in general position, then the facets of conv(S) are simplices and so the faces of conv(S) are trivially obtainable if the vert(f) is known for each facet f of conv(S). This information is also sufficient to give the adjacencies between facets: Two facets are adjacent iff they share a common ridge: they have dim(P)-2 vertices in common. Such pairs of facets can be found relatively quickly using radix-sort-based methods. Before discussing a randomized incremental algorithm:

Delaunay triangulations and convex hulls

As I've mentioned before, there is a very close relation between Delaunay triangulations and convex hulls; the Delaunay graph of a set of sites in d dimensions can be found by computing the convex hull of a corresponding set of sites in d+1 dimensions.

How does this work for planar Delaunay graphs?

Define a mapping dlift on points in the plane as follows; for a point x in the plane, let

dlift(x) := [x x2] = [x x·x] = (x1, x2, x12+x22),
when p=(x1,x2).

The image of the plane under dlift, or

D := dlift(E2) = {dlift(x) | x a point},
is a paraboloid of revolution.

Say that a plane h is nonvertical if h does not contain two points with the same x and y coordinates. For any such plane, the idea of above and below is well-defined: a given point in 3-space is either above, on, or below a nonvertical plane.

Lemma. The projection of the intersection of a nonvertical plane h with D is either empty, a point, or a circle.

Proof. A nonvertical plane h can be written as

h={[x z] | x in E2, z = x·a + b},
or as
h={[x z] | x in E2, z = 2x·p - p·p + c}, w
here p := a/2 and c := b + p·p.

The vertical distance between a point [x z] on h and the paraboloid D is x·x - z, or

x·x - 2x·p + p·p - c = (x-p)·(x-p) - c = r2 - c,
where r2 is the squared distance from x to p.

The relation of h and D depends on c:
c < 0: the vertical distance from D to any point on h is positive, and h is below D everywhere.
c=0: the vertical distance is zero at x=p, and positive everywhere else: that is, h is tangent to D at [p, p2], and the intersection of h and D is the point [p p2].
c > 0: since the vertical distance is r2-c, it is zero when r=sqrt(c), that is, when the point x is at distance sqrt(c) from p. Otherwise:
vertical distance < 0: holds when r2 < c, so D is below h, and x is closer to p than sqrt(c).
vertical distance > 0: holds when r2 > c, so D is above h, and x is farther to p than sqrt(c).
QED
lower face: a face f of conv(dlift(S)) is a lower face if there is a witness plane h for f such that h is nonvertical and dlift(S) is above or on h.
proj(P): For a set P in E3, the set {(x1, x2) | there is z with (x1, x2, z) in P}

Theorem. Let f be a lower face of conv(dlift(S)). Then proj(f) is a face of the Delaunay triangulation DG(S), and every face of DG(S) is such a projection.

Proof. Since vert(f) is a subset of dlift(S), and f is a lower face, there is some witness plane h for f, which means that vert(f) is a subset of h as well. Hence there is some circle such that the projection of vert(f) is on the circle, and all other sites are outside the circle. Hence proj(f) is part of the Delaunay triangulation of S.

The converse is easily shown, using the circle containing the Delaunay cell vertices to define the needed nonvertical witness plane. QED

Voronoi diagrams and halfspace intersections

For a site s and point x in E2, define
hs(x) := 2 x·s - s2
hs := {[x z] | z>=hs(x)}
HS := {hs | s in S}
Theorem. Let f be a face of P(HS) with f a subset of hs for s in some set of sites S', then proj(f) is the Voronoi region VR(S').

Proof. For two sites s and s',

hs(x)-hs'(x) = 2 x·s - s·s - 2 x·s' - s'·s'
= (x-s)2 - (x-s')2.

That is, if the vertical distance of x to hs is greater than the vertical distance of x to hs', the x is farther to s than to s'.

Thus if x in on a face f of P(HS), where x is on the bounding planes of hs for s in some subset S', then x is is equidistant to all s in S', and farther to any site in S\S' than to any site in S'. QED

Duality, VD(S), and DG(S)

The duality between the Voronoi diagram and the Delaunay graph can now be expressed in a way similar to polarity.

For a nonvertical upper halfspace

h = {[x z] | z >= hn·z - hz},
and point p = [hn hz], let p#:= h and h# := p.

Then a point q is in h iff h# is in q#, and so a given halfspace h contains a set of points T iff h# is contained in each of the halfspaces T#; that is, h contains conv(T) iff h# is in P(T#).

Notice that the point [0 inf] is in any upper halfspace h; just at the condition that 0 be in polytope P is needed for polarity between polytopes, so the condition that the T include [0 inf] is needed for this duality mapping.

With this inclusion, the duality between DG(S) and VD(S) is then the "projection" of the duality between the lower hull of conv(dlift(S) union [0 inf]) and P(dlift(S)#).

A randomized incremental algorithm for convex hulls

We suppose for now that S is in general position: no d+1 sites are on a common hyperplane.

Given a set S in d dimensions, the randomized incremental algorithm maintains conv(Ri), for i=1..n, where Ri is a random subset of S.

As in the planar case, adding a site to conv(Ri) has a search phase and an update phase:
search: find a facet of conv(Ri-1) visible to pi
update: find all the new facets of conv(Ri) incident to pi.

The algorithm of Lecture 1, and that given in the text, maintain for each pk with k>j, a visible facet of Rj. This makes the search phase easy, but is offline: the sites pk must be considered throughout the algorithm.

The randomized incremental algorithm for trapezoidal diagrams maintains a data structure for speeding search; the data structure uses the history of TD(Ri): the trapezoids of TD(Rj) for j
Here we use the history of the convex hull construction in a particularly convenient way: we keep a triangulation TR(Ri) of Ri which encodes the history of the construction of conv(Ri), as follows: when pi is added, for each facet f of conv(Ri) visible to pi, a simplex t=conv(f union {pi}) is added to TR(Ri). The facet f will be called the base facet of t, and pi will be called t's peak vertex.

Representation

The triangulation TR(Ri) is represented so that for each simplex t in TR(Ri), t knows each of its vertices v , and the corresponding simplex tv opposite v, where tv has all the vertices of t except v.

This representation is enough to navigate through the triangulation.

There is an initial root simplex, which is the convex hull of Rd+1.

The facets of the convex hull of Ri can be represented as simplices of this triangulation: they will have a peak vertex which is null.

Search

To search TR(Ri-1) to find a facet, start at the root simplex, and do the following: 
repeat
if
	a simplex is the root,
	or has a peak vertex on the same side of it base facet as pi,
then
	search the neighbors of the simplex.
until the visited simplex has a null peak.

Say that a non-root simplex of TR(Ri-1) visited under these conditions has a visible base facet.

Suppose that a simplex has peak vertex pj. The search condition means that if pi had been added at time j, then the base facet would have been visible to pi. That is, the searched simplices correspond to facets of the hull at different times that would be visible to pi.

For this procedure to be correct in finding a currently visible facet, it's enough to show that such a facet is connected to the root in the adjacency graph on the simplices, where the path comprises simplices with visible base facets. Such a path can be found by walking the line segment between pi and some point inside the root simplex. Each simplex along the walk has a visible base facet.

Update

The updating procedure considers each simplex with a visible base facet and null peak vertex, and changes the peak vertex to be pi.

Some additional simplices are created: those corresponding to facets of conv(Ri) incident to pi. Suppose r is a ridge contained in a facet visible to pi and also contained in a facet not visible to pi. Then conv(r union {pi}) is a facet of conv(Ri); considered in the polar, ro is a line segment between two vertices, and pi* has a bounding plane cutting ro, yielding a new vertex.

Analysis of Update

Suppose the expected number of facets of conv(Ri) is fi. We are interested in bounding the total number of base facets generated up to time i.

Theorem. The expected number of base facets of TR(Ri) is no more than
sumj<=i dfi/i.

Proof. We count the expected number of convex hulls facets created when pi is added:

= the expected number of facets of conv(Ri) incident to pi
= the expected number of facets of conv(Ri) incident to a random element of Ri
= dfi/i, since each facet is incident to d elements of Ri

Summing up to i gives the theorem. QED

Analysis of search

Theorem. If the expected number of facets of conv(Ri) is fi, then the expected number of visible base facets for pi is no more than
sumj<=i d(d-1)fj/j(j-1).

Proof. Let M be the set of base facets for the sequence p1...pi-1, and let M' be the set of base facets for the sequence pi, p1...pi-1; that is M' is the set of base facets in the "alternative" history where pi is added to the random subset first.

The set M\M' is the set of base facets visible to pi (hence not in M').

The set M'\M is the set of base facets incident to pi in the alternative history.

We have

|M| + |M'\M| = |M'| + |M\M'|,
and we want to know E|M\M'|, which leads to the problem of bounding
E|M| - E|M'| + E|M'\M|.

From the analysis of updates, we have

E|M| - E|M'| = sumj < i dfi/i - sumj<=i dfi/i = -dfi/i.

To bound the expected number of base facets in M'\M, we consider the expected number created when pj is added. Let R'j denote Rj union {pi}. We want the expected number of facets of R'j incident to both pj and to pi.

Since pi and pj are random elements of R'j, this is

d(d-1)fi/i(i-1),
since each facet is incident to {d choose 2} pairs of sites, and there are {i choose 2} equally likely pairs of sites, and we are averaging over all random subsets R'j.

Summing up to i gives the theorem. QED

From this analysis, the following results are easy:

Theorem. A randomized incremental algorithm can be used to compute the convex hull of a set of n points in d dimensions in Theta(nfloor(d/2)) expected time, and the Delaunay triangulation in Theta(nceil(d/2)) expected time. If the point set is such that the expected size of the convex hull of a random subset is O(i) for i=1..n, then the expected time needed by the algorithm is O(n log n).