Convex hull in higher dimensions, finding the vertices of a polytope

Question

Suppose I have a point cloud given in 6-dimensional space, which I can make as dense as needed. These points turn out to lie on the surface of a lower-dimensional polytope (i.e. the point vectors (x1, x2, ... x6) appear to be coplanar).

I would like to find the vertices of this unknown polytope and my current attempt makes use of the qhull algorithm, via the scipy interface in Python. In the beginning I would only get error messages, apparently caused by the lower dimensional input and/or the many degenerate points. I have tried a couple of brute-force methods to eliminate the degenerate points, but not very successful and so in the end, I suppose that all these points must lie on the convex hull.

This question has been very helpful, as it suggests a dimension-reduction via Principal Component Analysis. If I project the points to a 4D hyperplane, the qhull algorithm runs without errors (for any higher dimension it does not run).

from scipy.spatial import ConvexHull
from sklearn.decomposition import PCA

model = PCA(n_components=4).fit(initial_points)
proj_points = model.transform(initial_points)
hull = ConvexHull(proj_points, qhull_options = "Qx")

The answer in above question mentions that the simplices need to be transformed back after one computes the convex hull of the projected points. But the qhull output consists of only the indices, and why would these not match the indices of the initial points?

Now my problem is that I do not know which precision to use to actually obtain the proper vertices. Regardless of how dense I make the point cloud, the obtained vertices differ with different precisions. For instance, for initial points in a (10000, 6) array I get (where E0.03 is the maximum for which this works):

hull1 = ConvexHull(proj_points, qhull_options = "Qx, E0.03")
print len(hull1.vertices)
print hull1.vertices

5
[ 437 2116 3978 7519 9381]

And plotting it in some (not terribly informative) projection of axes 0,1,2 (where the blue points represent a selection of the initial point cloud):

But for a higher precision (of course) I get a different set:

hull2 = ConvexHull(proj_points, qhull_options = "Qx, E0.003")
print len(hull2.vertices)
print hull2.vertices

29
[  74   75  436  437  756 1117 2116 2366 2618 2937 3297 3615 3616 3978 3979
 4340 4561 4657 4659 4924 5338 5797 6336 7519 7882 8200 9381 9427 9470]

Same projection (just slightly different angle):

I would suspect that the first picture has not enough vertices and that the second picture possibly has too many. Though of course I cannot extract rigorous information from these plots. But is there a good way of finding out which precision to use? Or perhaps a completely different approach to this problem (I tried a few already)?

Lavine

Fascinating question. I don't have a ready answer, but agree that the first example certainly looks (by eye) to have too few vertices. The later one, I guess, will always tend to have lots of vertices along the "edges" (excuse me if bad terminology, not my field of expertise) of your projected polytope just because the initial points are random - you're unlikely to get one on the "true" vertex of the polytope which you're saying exists. If you've time to experiment - have you tried the Q8 option which seems to ignore "nearly inside" points. – Feld

Thanks for commenting. Turns out that most of the different options in Qhull yield the same (varying) answers, as does Q8. The only one which gives a slightly different number (but still unstable with different precisions) is Q9. It's correct that the set will unlikely contain the expected "true" vertices, yet since they come so very close, I feel like there should be a way to obtain them. – Lavine

The more I think, the more I'm intrigued. It seems this is still the subject of papers in maths. This shows an approach (2D), where their alpha parameter seems to have a similar effect to your precision. The problem is, the hull is by definition the smallest polytope that can contain the points and yet we're assuming the "true" vertices might lie outside the sample set and that, in some sense, the polytope has a "simpler shape" than that produced by the high precision estimate. By eye, OK, algorithmically, difficult. – Feld

Yes, I see the resemblance in the linked paper, quite nice, even if the 'inapproachable-ness' only arises in higher dimensions. Moreover because every other approach of comparing the points has failed (i.e. dot products, etc). In some sense the convex hull approach to this problem is most likely not what the algorithm is optimized to do, if only because every point in this point cloud will actually be a part of the hull. Yet I am not a programmer, and so I thought I'd double check on here if I am missing something. Seeing that this may very well be open is a finding, too, I suppose. – Lavine

I hadn't quite grasped the significance of all the points being on the hull: maybe techniques to identify (hyper)planes in a point cloud could be used. The intersection of such planes might give you the simple hull you are looking for. Could check afterwards that all points were on or inside. I found the RANSAC algorithm cited for this. 1, 2 – Feld

Oh, great, thank you! I was in the process of finding the normal vectors to each plane by hand, but the RANSAC algorithm seems promising, too. I'll look into it. – Lavine

No problem - let me know if it works. If it does, we can edit the discussion into an answer in case it's useful for anyone in future. – Feld

Even though late, I must admit that I did not find a way to implement the RANSAC algorithm in higher dimensions. Perhaps there is a way to do it, by generalizing the ideas given, but I couldn't quite get a handle on it. – Lavine

No matter, I think the algorithm described by @timothyshields below does what you want using gradient descent instead. – Feld

Answer 1

In this answer, I will assume you have already used PCA to near-losslessly compress the data to 4-dimensional data, where the reduced data lies in a 4-dimensional polytope with conceptually few faces. I will describe an approach to solve for the faces of this polytope, which will in turn give you the vertices.

Let x_i in R⁴, i = 1, ..., m, be the PCA-reduced data points.

Let F = (a, b) be a face, where a is in R⁴ with a • a = 1 and b is in R.

We define the face loss function L as follows, where λ⁺, λ^- > 0 are parameters you choose. λ⁺ should be a very small positive number. λ^- should be a very large positive number.

L(F) = sum_i(λ⁺ • max(0, a • x_i + b) - λ^- • min(0, a • x_i + b))

We want to find minimal faces F with respect to the loss function L. The small set of all minimal faces will describe your polytope. You can solve for minimal faces by randomly initializing F and then performing gradient descent using the partial derivatives ∂L / ∂a_j, j = 1, 2, 3, 4, and ∂L / ∂b. At each step of gradient descent, constrain a • a to be 1 by normalizing.

∂L / ∂a_j = sum_i(λ⁺ • x_j • [a • x_i + b > 0] - λ^- • x_j • [a • x_i + b < 0]) for j = 1, 2, 3, 4

∂L / ∂b = sum_i(λ⁺ • [a • x_i + b > 0] - λ^- • [a • x_i + b < 0])

Note Iverson brackets: [P] = 1 if P is true and [P] = 0 if P is false.

Coomer

There are a couple of details that still elude me with respect to the actual application, but in theory I really like this idea. To begin with, I must admit that it is not quite clear to me how to obtain the faces. Can you perhaps explain (or point me to some source), where the condition a • a = 1 stems from? – Lavine

The vector a is the normal of the (directed) face F. The condition a • a = 1 just restricts the normal a to have a length of 1; that is, ||a|| = 1. – Coomer

Ah, right. That makes sense. What exactly is b determined by? I haven't fully gotten around to implement the idea, but I suppose it makes sense to accept this. – Lavine

b is the offset of the face F along its normal a. When b = 0 the face goes through the origin. When b = 5 the face is a distance of 5 away from the origin. – Coomer

@Lavine I added the partial gradients to the end of the answer. – Coomer