n维点集的凸包顶点

Question

我有一个在n维空间中给出的点集。我想要找到其中构成凸包的顶点（角落）。

我希望使用Python来解决这个问题（但也可以调用其他程序）。

编辑：所有坐标都是自然数。输出应该是顶点的索引。

通常通过搜索引擎查找出来的答案只能解决2D情况，或者要求列出面，这在计算上更加困难。

我目前的尝试：

• scipy.spatial.ConvexHull(): Throws error for my current example. And I think, I have read, it does not work for dimension above 10. Also my supervisor advised against it.

• Normaliz (as part of polymake): works, but too slow. But maybe I did something wrong.

  import PyNormaliz
  def find_column_indices(A,B):
    return [i for i in range(A.shape[1]) if list(A[:,i]) in B.T.tolist()]
  
  def convex_hull(A):
    poly = PyNormaliz.Cone(polytope = A.T.tolist())
    hull_cone = poly.IntegerHull()
    hull_vertices = np.array([entry[:-1] for entry in hull_cone.VerticesOfPolyhedron()]).T
    hull_indices = find_column_indices(A, hull_vertices)
    return hull_indices
  

• Solve with linear programmes: works, but completely not optimised, so I think there must be a better solution.

  import subprocess
  from multiprocessing import Pool, cpu_count
  from scipy.optimize import linprog
  import numpy as np
  
  def is_in_convex_hull(arg):
    A,v = arg
    m = A.shape[1]
    A_ub = -np.eye(m,dtype = np.int)
    b_ub = np.zeros(m)
    res = linprog([1 for _ in range(m)],A_ub,b_ub,A,v)
    return res['success']
  
  def convex_hull2(A):
    pool = Pool(processes = cpu_count())
    res = pool.map(is_in_convex_hull,[(np.delete(A,i,axis=1),A[:,i]) for i in range(A.shape[1])])
    return [i for i in range(A.shape[1]) if not res[i]]
  

例子：

A = array([[  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1],
...:        [ 0,  0,  0,  0,  0,  0,  0,  0,  0,  2,  2,  2,  2,  4,  6,  6,  6,  8,  1,  2,  2,  2,  2,  3,  2,  1,  2,  1,  2,  1,  1,  1,  1,  2,  2,  2,  3,  1,  2,  2,  1,  2,  1,  1,  1,  2,  1,  2],
...:        [ 0,  0,  0,  0,  0,  2,  2,  4,  6,  0,  0,  2,  4,  0,  0,  2,  2,  2,  1,  2,  1,  1,  1,  1,  1,  2,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  0,  2,  0,  1,  0,  1],
...:        [ 0,  0,  2,  4,  4,  2,  2,  0,  0,  0,  6,  2,  0,  4,  0,  2,  4,  0,  1,  1,  2,  2,  2,  1,  1,  1,  2,  1,  1,  1,  2,  1,  1,  2,  1,  2,  1,  1,  1,  2,  1,  1,  1,  1,  1,  1,  1,  1],
...:        [ 0,  6,  0,  2,  4,  0,  6,  4,  2,  2,  0,  0,  8,  4,  8,  4,  0,  2,  4,  2,  2,  2,  2,  2,  2,  2,  2,  3,  3,  2,  2,  2,  2,  2,  4,  2,  2,  1,  1,  1,  2,  3,  2,  2,  2,  2,  1,  2],
...:        [ 0,  2, 14,  0,  4,  6,  0,  0,  4,  0,  2,  0,  4,  4,  4,  0,  0,  2,  2,  2,  2,  2,  2,  1,  2,  4,  1,  3,  2,  1,  1,  1,  1,  2,  1,  4,  1,  1,  0,  1,  1,  1,  2,  3,  1,  1,  1,  1],
...:        [ 0,  0,  0,  2,  6,  0,  4,  6,  0,  0,  6,  2,  2,  0,  0,  2,  2,  0,  1,  1,  2,  1,  2,  1,  1,  1,  1,  1,  1,  1,  2,  2,  1,  1,  1,  1,  1,  1,  1,  2,  1,  1,  1,  1,  1,  1,  0,  1],
...:        [ 0,  2,  2, 12,  2,  0,  0,  2,  0,  8,  2,  4,  0,  4,  0,  4,  0,  0,  2,  1,  2,  1,  1,  2,  1,  1,  4,  2,  1,  2,  3,  1,  3,  2,  2,  2,  1,  1,  2,  1,  1,  1,  1,  1,  3,  1,  1,  2],
...:        [ 0,  8,  2,  0,  0,  2,  2,  4,  4,  4,  2,  4,  0,  0,  2,  0,  0,  6,  2,  2,  1,  1,  1,  3,  2,  2,  1,  2,  2,  1,  2,  1,  2,  1,  3,  1,  2,  1,  1,  1,  1,  3,  1,  3,  2,  2,  0,  1]])

产出运行时间

In [44]: %timeit convex_hull(A)
1.79 s ± 16.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [45]: %timeit convex_hull2(A)
337 ms ± 3.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

对于稍微大一些的例子，比率更糟糕，因此也不能通过并行化来解释。

非常感谢任何帮助和改进。

- Hennich

---

问题是，您是否需要比3更多的维度，这是一个CG程序吗？您希望有多精确？ - user1767754

是的，我需要任意大的维度。这个问题来自于多项式优化。我的坐标都是自然数，答案是一个索引列表，因此最好使用精确算法。 - Hennich

这将是计算上复杂的，如果你在这里看 https://cw.fel.cvut.cz/wiki/_media/misc/projects/oppa_oi_english/courses/ae4m39vg/lectures/05-convexhull-3d.pdf ，只有quickhull表现得有点“快”，但对于多维来说可能会很棘手。 - user1767754

@user1767754，就我所见，此次演讲旨在寻找我未能找到的方面。通过上述线性规划问题，我可以在多项式时间内解决我的问题。 - Hennich

Answer 1

你可以通过以下方式修改你的is_in_convex_hull方法：

def convex_hull(A):
    vertices = A.T.tolist()
    vertices = [ i + [ 1 ] for i in vertices ]
    poly = PyNormaliz.Cone(vertices = vertices)
    hull_vertices = np.array([entry[:-1] for entry in poly.VerticesOfPolyhedron()]).T
    hull_indices = find_column_indices(A, hull_vertices)
    return hull_indices

该算法的规范化版本将运行得比原来快得多。

- sebigu