Approximate Convex Hulls: sketching the convex hull using curvature
Robert Graham Adam M. Oberman

Approximate Convex Hulls]Approximate Convex Hulls: sketching the convex hull using curvature


Department of Mathematics and Statistics, McGill University.  This work was supported by an NSERC Engage grant 2015


Convex hulls are fundamental objects in computational geometry.  In moderate dimensions or for large numbers of vertices, computing the convex hull can be  impractical due to the computational complexity of convex hull algorithms.  In this article we approximate the convex hull in using a scalable algorithm which finds high curvature vertices with high probability.   The algorithm is particularly effective for approximating convex hulls which have a relatively small number of extreme points. 





[
    Adam M. Oberman
    December 30, 2023
=====================




















§ INTRODUCTION


Computing the convex hull of points is a fundamental problem in the field of computational geometry <cit.>.
In moderate dimensions or for large numbers of vertices, computing the exact convex hull can be computationally impractical.  Even the vertex redundancy problem, which computes the extreme points without the full geometric structure of the convex hull is impractical 

Even relatively modern algorithms <cit.> <cit.> break down or are too expensive for high dimensional problems.  For practical purposes, in moderate dimensions, the computational obstruction is not with the algorithm, but with the convex hull itself.  However, in moderate or high dimensions, the combinatorial structure of high dimensional data sets can be very different from the geometrical intuition obtained from studying low dimensional simplices <cit.> <cit.>.

Nevertheless, there are many practical problems which are well represented by sampling a large number of high dimensional points, where we expect that either: (i) the number of extreme points is small or, more generally, (ii) the data is contained in a set which is the convex hull of a small number of points.   The second case applies, for example, to data sampled uniformly from a simplex in high dimensions.  This problem occurs hyperspectral data analysis <cit.> <cit.> <cit.>. 


In figure <ref>, we give an illustration of this problem, and of our algorithm.  Here 200 points were sampled uniformly from the 2 dimensional simplex. We used the hyperplane compression algorithm to approximate the original simplex.

Our work is motivated by an application in data reduction.  The algorithm was used to reduce the number of points in a helicopter flight test of Bell Helicopter Model 505 Jet Ranger X <cit.>.
Approximately twenty million load vectors were recorded during the helicopter certification flight test. The goal
was to extract a small number of extremal loads to be applied on a flight representative fatigue test of the helicopter tailboom assembly. The
exact convex hull contained approximately two thousand points. The algorithm
identified approximately 200 high curvature points. These points were clustered
into a smaller number of points, and finally, using strain response of the tailboom to the load vectors, six load vectors were selected for the fatigue test. This robust data reduction method proved to be effective in identifying extremal loads and was instrumental for the timely certification of the recently released Bell Helicopter 505 Jet Ranger X.

The key idea in our approach is geometrically intuitive. Suppose we have a finite collection of points X⊂ℝ^n for which we want to approximate the convex hull, denoted CH(X). Given any unit vector d∈ S^n-1, any extreme point

    x_d ∈{x^⊺ d  |  x∈ X}


is a vertex of CH(X), moreover x^⊺ d ≤ x_d^⊺ d defines a  supporting hyperplane. Thus if we perform such computations for a large number of unit vectors, we obtain both a collection of vertices whose hull is contained in CH(X) and a collection of hyperplanes whose intersection contains CH(X). These determine inner and outer approximations to CH(X). 

However, if we return to the example of points sampled from a simplex, near a vertex, there can be a high number of extreme points.   Our goal is to reduce the number of extreme points without introducing too large of an error in the convex hull.   The influential Pixel Purity Index (PPI) algorithm <cit.> keeps only vertices which are the extremal for more unit vectors, d.   In this article, we justify and refine the algorithm using high dimensional curvature concepts. 
The connection with curvature has been hinted at in <cit.>, where it was observed that points which maximize many direction vectors are “presumed to be closer to the “corners" of the data".  As far as we know, we are the first to make this explicit.  Using this observation we can prove consistency <ref> and convergences <ref> of the algorithm. As an aside we give a theoretical answer to a question posed in <cit.>.


We extend the algorithm by also giving a method for reducing the number of hyperplanes <ref>. We also briefly touch on the endmember detection problem: given a collection of points uniformly sampled from some polytope, how do we find the corners of the polytope (as opposed to the corners of the convex hull of the points).

Our paper is structured as follows. In Section 2 we explain how our method computes approximate curvature. Section 3 describes how to compute the error in our approximations. Section 4 describes our algorithm and its extensions. Section 5 uses  the results of Section 2 and 3 to show our algorithm is consistent and converges. Finally in section 6 we give some examples.


Note there is a large body of work on approximately convex bodies, see for example <cit.> and <cit.>.  There the idea is to, for example, find the best ellipsoid inside or containing a convex body, where a scaling factor is allowed.  The goal of this work it to select extreme vertices in a computationally practical fashion, which is somewhat different from those works. 	







§ CONSISTENCY OF THE CURVATURE APPROXIMATION


Here we define the basics and show how one can compute the curvature of polytopes. This is the essence of our algorithm defined in a later section.


 A supporting hyperplane at v  for the convex set P is a hyperplane with normal n such that 
 
    n^⊺ y ≤ n^⊺ v,    for all y∈ P.


 Let S^n-1 be the unit sphere in 

Let N_v⊂ S^n-1 be the set of all unit normals of supporting hyperplanes at v for the convex set P . We view the normal vectors as being points on the sphere. The curvature of P at v is the spherical volume of N_v. The relative curvature of P at v is 

    K(v) = (N_v) /(S^n-1).





In two dimensions the curvature is a measure of the exterior angle at a point. In general, curvature is viewed as a n dimensional notion of spherical angle.	 An introduction to curvature of polytopes can be found in <cit.>. We note this definition applies just as well to any convex body <cit.>. For smooth bodies, the curvature of a subset V of a convex body is ( ⋃_v∈ V N_v).  


We can now show how to approximately measure the curvature of polytopes.


Let D ⊂ S^n-1 be finite and let V ⊂ also be finite.  
Let 

    D_v := { d∈ D |  v^⊺ d  ≥ w^⊺ d   for all w ∈ V}
 
be the set of extremal directions in D for v, and let 

    CH_D(V) = { v ∈ V |D_v > 0 }

be the D-convex hull of V. 
Define the relative D-curvature of S at a D-extremal point v to be

    K_D(v) = D_v/D




Let  D be a finite set of vectors uniformly sampled from S^n-1.  
Then for any ϵ >0

    (  K_D(v) - K(v)  >ϵ ) ≤K(v)(1-K(v))/Dϵ ^2





First note that since every direction vector is maximized by some vertex of a polytope we know the sets { N_v | v∈ V } cover the sphere. On the other hand for w≠ v, N_v ∩ N_w has measure zero. 
This follows since for a given normal vector in N_v ∩ N_w, the corresponding supporting hyperplane contains the line segment between v and w. Consequently the dimension of this set must have measure zero. 

Hence, any d_i has probability K(v) of satisfying d_i ∈ N_v.  Moreover, for each d_i belongs to N_v for some v (or with probability 0 it is maximized by more than one). Therefore D_v follows a multinomial distribution; every d_i has probability K(v) of belonging to N_v. Now, given that K_D (v) = D_v/D, Chebyshev's inequality states that
 
    ( | K_D(v) - E(K_D(v)) | >ϵ ) ≤var( K_D(v))/ϵ ^2

and since D_v follows a multinomial distribution E(K_D(v))=K(v) and var(K_D(v)) = K(v) (1- K(v))/D. The result follows immediately. 


Figure <ref> shows the convergence of K_D (v) for a given vertex as |D| increases.







§ SPARSE APPROXIMATION OF POLYTOPES


Here we show that removing low curvature extreme points from a polytope does not significantly change its shape (more precisely stated below).   The main result is related to the Aleksandrov's maximum principle <cit.>




Let d(S,S') denote the Hausdorff distance between two sets

    d(S,S') ≡max{sup_x∈ S d(x,S'), sup_y∈ S' d(y,S)   }


We need a simple lemma about the Hausdorff distance between two polytopes.
 Let S = CH({v_1,...v_m }) and let S' be convex and compact then 

    sup_x∈ S d(x,S') = sup_1≤ i ≤ m d(v_i,S')



Let W be a point in S such that d := d(W,S') is maximal. Let A be a point in S' such that dist(W,A)= d (which exists by compactness).  Consider the hyperplane ℋ at A with normal WA. We claim this is a supporting hyperplane at A with respect to the polytope S'. Suppose otherwise, then there exists a point B∈ S' that belongs to the same side of the hyperplane as W. But then the line AB belongs to S' and this line will intersect the ball of radius d centered at W and so there is a point in S' that is closer to W then A: contradiction. The following figure makes this clear (the circle pictured is of radius d).


    < g r a p h i c s >


On the other hand there must be some w_i on or above (i.e. not on the same side as A) the hyperplane at W with normal WA since w_i are the extreme points of S. If w_i is on the hyperplane then d(w_i,S')≥ d and we are done. Otherwise, if it is above, then since d(w_i,S')≤ d there exists a point A_0 in S' such that d(W,A_0) ≤ d and so A_0 is on the wrong side of ℋ giving the above contradiction.

 
	Recall for an angle θ∈ [0,π ] and v∈ S^n-1 a spherical cap with angle theta about v is the set 
	
    {w∈ S^n-1| w^⊺ v ≥cos (θ ) }

	
	Let S_cap(θ) denotes the volume of the n dimensional spherical cap with angle θ (about any v since this will not change the volume). 
	
	.  Let θ∈ [0,π ], then in S^n-1
	
    1/2( sin( θ/2 ) 
    )^n-1≤ S_cap (θ)

	

 A straightforward rewriting of  <cit.>


Let B_r denote the ball of radius r in . 

 Let V = {v_1,...v_m } and W = {w_1,...,w_k,} be subsets of  B_r. Let S = CH(V ∪ W) and S' = CH(V), suppose neither are degenerate.  Let ω be the sum of the curvatures of all the points w ∈ W. Then 

    S_cap(arcsin(d(S,S')/2r))≤ω

and 

    d(S,S') ≤√(2)π r(2ω)^1/n-1




The result above gives an upper bound on the distance of order ( ω (S^n-1))^1/n-1 where ω is the total relative curvature removed from the set. In Figure <ref> we plot this function for n=2,3,4,5 and r = 1. Note since 2 is a trivial worst case this shows how the estimate is only useful for small ω








Let d = d(S,S'). First we show how to obtain (<ref>) from (<ref>).   From lemma <ref>

    1/2( sin( arcsin(d/2r)/2) )^n-1≤ S_cap(arcsin( d/2r) ) ≤ω
 
hence: 

    sin( arcsin(d/2r)/2) ≤ 2 (2ω)^1/n-1

Finally we use the inequalities
 
    sin( arcsin(d/2r)/2) ≥4/√(2)π(arcsin(d/2r)/2)
 
 and 

    4/√(2)πarcsin( d/2r) ≥4/√(2)π(d/2r)
 
to complete the proof of this step. 




Next we prove (<ref>).  First we show it suffices to consider the case where k=1.  Suppose we proved the theorem for k=1. Let w_i be the point in S furthest away from conv({v_1 ,...v_m}) which exists by Lemma <ref>. But by assumption, S_cap(arcsin(d(w_i,S')/2r) ) ≤ω'  where ω' is the curvature of w_i with respect to S”:=conv(w_i,v_1,...v_n). Thus it would be enough to show ω' ≤ω. But note every supporting hyperplane for some w_j with respect to S will either be a supporting hyperplane for some v_ℓ or for w_i with respect to S”, on the other hand any supporting hyperplane for some v_j with respect to S will remain a supporting hyperplane of v_j with respect to S”. From this we conclude that ω' ≤ω as required. 

It remains to show the case k=1. Rename w_1 to W for clarity, so S=conv(W,v_1,...v_n), note d = d(W,S').  Let A be the unique point in S' such that dist(A,W)=d (which exists by compactness of convex hulls). Let  ℋ be the hyperplane at A with normal WA. Note that exactly as in the proof of  <ref> ℋ is a supporting hyperplane for S'

Consider the cone C_0 consisting of vertex W and base S∩ℋ. In other words the portion of S on the same side of the hyperplane as W. Let C_1 be the cone consisting of vertex W and base the ball of radius  √( (2r)^2 - d^2) centered at A in  ℋ. We claim C_0⊂ C_1, it suffices to show S∩ℋ is contained in the ball around A. Note by assumption all points of S are in a ball of radius r, in particular all points of  S∩ℋ are within 2r of W. Moreover the slant of C_1 is length 2r. So if any point in S∩ℋ were outside the ball around A it would be a distance >2r from W: contradiction.

Finally, since C_0⊂ C_1 it follows that the curvature of W with respect to C_0 is less then or equal to the curvature of W with respect to C_1 but this is precisely:

    S_cap(arcsin ( d/2r)) ≤ω










§ THE ALGORITHM





Algorithm <ref> is our basic approximate convex hull algorithm. In short, we generate many direction vectors and use these to compute curvature as described above. We then keep only the high curvature points.





Suppose we run algorithm <ref> on sets V and D as above. Let V' be all the vectors that were kept. We call CH(V') the inner hull. Now for each d∈ D let v_d∈ V be a vector that maximizes the dot product. Consider the collection of linear constraints 

    d^⊺ x ≤ d^⊺ v_d

This determines a convex body we call the outer hull. 

Note the inner hull is contained in the actual convex hull of V, which is contained in the outer hull. Our algorithm also gives the constraints for the outer hull. In step 2 one simply needs to keep track of the value d ^⊺ v_d (i.e a_v,d) for each d∈ D when computing the maximums. There is a possibility the constraints will not define a finite polytope if |D| is too small but this is exceedingly unlikely for large 
|D|. An even larger concern is that the outer hull contains a large number of constraints, we will attempt to remedy this later. 

We claim that both the inner and outer hull approach the actual convex hull if one increases the number of direction vectors. Thus we get an approximate convex hull in the vertex format and another approximate convex hull in the constraint format. The following notion of error will help make this precise. 


Under the same assumptions as the previous definition let A be the inner hull, B the actual hull and C the outer hull. We define the inner error as 

    sup_x∈ B d(x,A)
 
and the outer error as 

    sup_x∈ C d(x,B)

In both cases these are simply the Hausdorff distance defined above.



In <cit.> we find an elegant method for producing direction vectors (called skewers). Not only does this speed up the production of direction vectors but more importantly it speeds up the computation of the dot products in the algorithm.  We are confident choosing vectors uniformly from the sphere works (Since it gives a precise measure of the curvature as we saw above). However alternative methods seem to work just as well experimentally (intuitively they are `uniform enough'). In <cit.> they have tested various methods similar to <cit.> and have singled out what they found to be the best approach. 




 §.§ Sparse Approximation

For many applications the above algorithm finds too many extreme vectors in the inner hull, and for almost any application the outer hull has far too many hyperplanes (there will be one hyperplane for each direction vectors). In this section we discuss how to deal with this problem, in particular how to reduce the following two ratios:


Let V ⊂ℝ^n be a finite collection of points. Suppose we have an algorithm A that outputs vertices of a convex hull A_CH that approximates CH(V). Then the vertex compression ratio of algorithm A is

    Number of vertices in  A_CH/Number of vertices in   CH(V)

On the other hand suppose A outputs hyperplanes of a convex hull A_CH that approximates CH(V). Then the hyperplane compression ratio of algorithm A is 

    Number of hyperplanes in  A_CH/Number of hyerplanes in  CH(V)




To reduces the vertex compression ratio of algorithm <ref>  we run algorithm <ref>. We simply find vectors that are clustered together and keep only one of them. 






Let V,V',V”,β be as in algorithm <ref> then 

    d(CH(V'),CH(V”)) < β



This is an immediate application of Lemma  <ref>

Next we show how to reduce the hyperplane compression ration of our original algorithm. Let V,D be as above, we start by running our original algorithm followed by the vertex compression algorithm. Suppose when running our original algorithm that for each v∈ V' we keep track of 

    D_v := { d∈ D |  v^⊺ d  ≥ w^⊺ d   for all w ∈ V}
 
moreover suppose when running the vertex compression algorithm that for each v∈ V” we keep track of 

    E_v := { w ∈ V' | w  was removed during the step involving  v }∪{v}

(i.e. E_v is all elements in V' that were clustered around v). For v∈ V” let 

    F_v := ∪_w∈ E_v D_v


See algorithm <ref> below.




If desired, algorithm <ref> works without running the vertex compression algorithm (i.e. if β = 0). If one doesn't run the vertex compression algorithm and finds the true convex hull of the projected F_v then the result is equivalent to the outer hull found originally (this merely removes redundant hyperplanes)




We end this section by noting that hyperplane compression algorithm is a potential solution to the problem of finding endmembers (the vertices of the outputed hull are potential endmembers). This holds even without a pure pixel in the data set. In particular if the data is uniformly generated from some polytope and we wish to recover the polytope this method can be used.




§ APPLICATION OF THE CURVATURE ESTIMATE TO THE ALGORITHM



We can use the consistency result, Theorem <ref>, to show convergence of algorithm <ref>.

Let V ⊂ℝ^n and let D be uniformly sampled from S^n-1.  Assume that the C(v) ≥ω > 0 for all extremal points in V.  
Then

    CH_DV = CH(V),    with probability P ≥ 1-p

provided 

    D≥log (ω p)/log(1-ω).



 
We have seen from the proof above that the probability of not finding a point with relative curvature ω is (1-ω)^D.Since there are at most 1/ω such points we have by subadditivity of probabilities that the probability of missing one of them is ≤1/ω (1-ω)^D (for an exact answer use inclusion-exclusion) Suppose p is less then or equal to this: p≤1/ω (1-ω)^D. This is trivially equivalent to D≥log (ω p)/log(1-ω) as required



	The function log (ω p)/log(1-ω) in the estimate above is O(1/ω) as case be seen in   Figure <ref>.





Another consistency result is as follows.



Let V ⊂ℝ^n. Let D be sampled uniformly from S^n-1. 
Let  O_k and I_k be the outer and inner error respectively. Then {O_k }_k∈ℕ and {I_k }_k∈ℕ are non increasing sequences and with probability 1 they both converge to zero.


It is clear that {O_k }_k∈ℕ and {I_k }_k∈ℕ are non increasing sequences since each new direction vector adds one more constraint (possibly lowering O_k and may or may not find a new extreme vector (possibly lowering I_k).

It is also clear from above that I_k approaches 0. Indeed each extreme vector has non zero curvature and hence a non zero probability of being found so for large enough k we expect I_k = 0 with increasing probability.

Finally we can use theorem <ref> to show with probability 1 that O_k converges to zero. As before we know that with probability 1 all extreme points will be found. The extreme points of the outer hulls are made up of V and some other vectors. For each extreme point v∈ V, we claim the curvature of v with respect to the actual hull CH(V) approaches the curvature of v with respect to the outer hull as k→∞ with probability 1. To show this let E be the set of all sequences (d_1,d_2,...) ∈ (S^n-1)^∞ such that the claim fails. We wish to show the measure of E is zero. Let N_v ⊂ S^n-1 be the collection of direction vectors for which v is extremal and let { d_1 ',d_2 ' ... } be the subsequence of d_i that belong to N_v. Now if each point in N_v is a limit point of ∪_n CH({ d_1 ' ,...d_n '} ) the claim clearly follows. Therefore for each sequence in E there exists a rational point q∈ N_v and a rational r>0 such that B_r(a)∩∪_n ({ d_1 ',...d_n '} ) = ∅. Let E_a,r be all sequences for which this holds, now since N_v ∩ B_r(a) has positive measure then by definition of product measure the measure of E_a,r is zero. Since E belongs to the union of all E_a,r the claim follows from countable subadditivity of measure.

 Now since the sum of the curvatures of each extreme point of any polytope add up to (S^n-1) the claim implies that all extreme points of the outer hulls that are not in V must have a vanishing proportion of the total curvature as k→∞. Theorem  <ref> completes the proof.

It is easy to compute some controls on the convergence of the inner error. Equivalenty this gives a worst case calculation for choosing the appropriate number of  direction vectors to achieve a desired error. This is a theoretical answer to a problem raised in <cit.>, how to choose the number of direction vectors for PPI algorithms. It's worth mentioning the practical solutions from <cit.> for choosing the number of direction vectors. Essentially they suggest computing the maximum for only a small block of direction vectors and then repeating this process until no new extreme vectors are found. This way there is no need for human input about the choice of how many direction vectors to use. 

Number of direction vectors needed for a given inner error (simple messy bound)
We can use the above to find the number of direction vectors needed to achieve a 
particular inner error ϵ >0. Let V be a set of points in ball of radius r and let X be the true extreme points of (V).   Suppose we run our algorithm with some set of direction vectors uniformly chosen and we keep all extreme points found.  The error depends on the total amount of curvature of all the extreme points we've missed. Using <ref> we can ensure with probability ≥ 1-p that we have all points of relative curvature ≥ω for any 0≤ω≤ 1. Then the total missing curvature would be ≤ |X|ω(S^n-1). By <ref> this will give an error 

    ≤√(2)π r (2|X|ω(S^n-1) )^1/n-1
. 

We can set this error ≤ϵ and solve for ω to get 

    ω≤(ϵ/√(2)π r)^n-1/ 2|X|  ((S^n-1) )

Denote the right hand side by C. From the result above to achieve this ω we would require 

    log (Cp) /log (1-C)

direction vectors. To sum up, if the number of directions is more than log (Cp) / log(1-C) then with probability 1-p the error is bounded by ϵ. 

The growth of this is comparable to (S^n-1)/(ϵ/√(2)π r)^n-1. This unfortunately gives results that are very large. For example keeping p=.05,ϵ = .1, r=1,  |X|=10000,  and varying the dimension n = 3,4,5,6,7 we get roughly 10^10,10^11,10^13,10^15,10^17 respectively.




Improvement on previous remark
The previous remark doesn't relate to the inner error we expect in practice. For one thing, of course we don't claim the constants found above are in any way optimal but more importantly if we remove many low curvature points we don't expect this to be equivalent to the worst case where all the curvature is concentrated in one point. This is especially true since we are taking our directions uniformly from the sphere. In many cases if we have found all points with relative curvature ≥ω then removing all other points results in an error closer to ≤√(2)π r (2 ω(S^n-1) ) ^1/n-1 (the curvatures don't add up). This means we can set |X|=1 in the calculation above. So if we consider as before  p=.05,ϵ = .1, r=1, and vary the dimension n = 3,4,5,6,7 we get roughly 10^5,10^7,10^9,10^11,10^13 respectively. 






§ COMPUTATIONAL EXAMPLES



 §.§ Synthetic data

In this section we generate synthetic data to test the algorithm. 

For points generated uniformly from a simplex most curvature is in the corners.  So already a random uniform direction vector will likely find a corner. Our procedure is effective in this case. 


Consider Figure <ref>. This is 1000 direction vectors and a million points. The algorithm found 61 extreme vectors (unfilled dots) and kept 34(filled dots). The actual convex hull has around 300 extreme points. Clearly clustering (algorithm <ref> would be very effective in this case.

On the other hand for a million points generated uniformly from the sphere we do very poorly (Figure <ref>). Here all the direction vectors found different extreme points. Only 102 were kept. (This can be slightly rectified by modifying step 4 at the expense of worse performance on examples with sharp corners).





Figure <ref> is an example of a mix of low curvature points with one very high curvature point. In practice much data is such a mix, so we expect high performance in some region and low performance in others.

Figure <ref> studies the accuracy of the algorithm as a function of the number of direction vectors. We generated 10 thousand 3 dimensional points and applied the algorithm with an increasing number of direction vectors, keeping every extreme point we found. The figure shows how the outer error (yellow) and the inner error (blue) decrease as the number of found extreme vectors increases. The points in the first two images were generated randomly from the cube and from the sphere. The points in the third  randomly generated points from a simplex, and the we applied a fixed linear transformation. (Note for the sphere the outer error is only computed to the first digit due to the large number of computations involved).  Figure <ref>(d) shows a similar computation of the inner error using 1 million points in 3 dimensions. This time the error is compared to the number of extreme vectors found.























 §.§ Algorithm Applied to Helicopter Flight Test data set.

In this section, we demonstrate the results of the vertex compression algorithm on a real dataset consisting of 23 million data points in dimension 5. We used 70000 direction vectors  In Figure <ref>(a) we compare the 1491 found points (red filled) with the 2495 points on the true convex hull. There is an error of around 300, for context the width of the shape is around 66000. Moreover the average distance from the true extreme points to the mean of the extreme points is 32000. The figure has been projected into 2 dimensions.    After  running the vertex compression algorithm we were left with 29 points <ref>(b). The error here is 6800, for context recall the information above.  
















alpha






