n维点集凸包的顶点

Question

我在维度 n 中有一组给定的点。在这些中，我想找到那些，它们是凸包的顶点（角）。我想用 Python 解决这个问题（但可能会调用其他程序）。

编辑：所有坐标均为自然数。作为输出，我正在寻找顶点的索引。

谷歌搜索通常会产生二维问题，或者要求列出面孔，这在计算上要困难得多。

到目前为止我自己的尝试

scipy.spatial.ConvexHull()：为我当前的示例抛出错误。而且我认为，我已经阅读过，它不适用于 10 以上的维度。我的主管也反对它。

Normaliz（作为 polymake 的一部分）：有效，但速度太慢。但也许我做错了什么。

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

用线性规划求解：可行，但完全没有优化，所以我认为一定有更好的解决方案。

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)

稍微大一点的例子，比例更差，所以也不能用并行来解释。

感谢任何帮助或改进。

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

算法的 Normaliz 版本运行会快得多。

n维点集凸包的顶点

Vertices of the convex hull of n-dimensional point set

python

convex-hull