Convex hull in n-dimensions or linear programming to find the vertices of a polytope point cloud.

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. 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.

At first, it a dimension-reduction via Principal Component Analysis was suggested to me. If I project the points to a 4D hyperplane, the qhull algorithm runs without errors (for any higher dimension it does not run).

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")

5
[ 437 2116 3978 7519 9381]

(Please note, I have been using different programs, so this is not matlab code, but hopefully still makes sense, especially given that they all make use of quickhull).

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")

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):

Alternative approach.

Someone suggested the following, which I believe is a linear programming problem, and seems to give a nice way of solving this. Unfortunately I know nothing about LP, and thus I do not know how to fill in the details and transfer the idea into matlab code. Could someone perhaps provide some pointers, so that even a LP-novice can compute this? For instance, I do not know if the a and b are unknowns and/or how to randomly initialize the faces F.

• 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 xi in R4, i = 1, ..., m, be the PCA-reduced data points.
• Let F = (a, b) be a face, where a is in R4 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) = sumi(λ+ • max(0, a • xi + b) - λ- • min(0, a • xi + 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 / ∂aj, j = 1, 2, 3, 4, and ∂L / ∂b. At each step of gradient descent, constrain a • a to be 1 by normalizing.
• ∂L / ∂aj = sumi(λ+ • xj • [a • xi + b > 0] - λ- • xj • [a • xi + b < 0]) for j = 1, 2, 3, 4
• ∂L / ∂b = sumi(λ+ • [a • xi + b > 0] - λ- • [a • xi + b < 0])
• Note Iverson brackets: [P] = 1 if P is true and [P] = 0 if P is false.

[jim as]

A possible solution maybe is provided with the algorithm, cMinMax, described in https://arxiv.org/abs/2011.14035