【问题标题】:Vertices of the convex hull of n-dimensional point setn维点集凸包的顶点
【发布时间】:2018-05-07 19:43:14
【问题描述】:

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

编辑:所有坐标都是自然数。作为输出,我正在寻找顶点的索引。

谷歌搜索通常会在 2D 中产生问题,或者要求列出人脸,这在计算上要困难得多。

到目前为止我自己的尝试

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

对于一个稍微大一点的例子,比率更差,所以它也不能用并行化来解释。

感谢任何帮助或改进。

【问题讨论】:

  • 问题是,你需要比3更多的维度吗,这是还是CG程序?您希望准确到什么程度?
  • 是的,我需要任意大的尺寸。问题来自多项式优化。我的坐标都是自然数,答案是索引列表,所以一个精确的算法会很好。
  • 这将是计算复杂的,cw.fel.cvut.cz/wiki/_media/misc/projects/oppa_oi_english/… 如果你看这里,只有 quickhull 行为有点“快”,但对于多维这可能会很棘手。
  • @user1767754 这个谈话,据我所知,想找到我没有的方面。我的问题可以通过上面的 LP 在多项式时间内解决。

标签: python convex-hull


【解决方案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 版本的算法会运行得更快。

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 2015-03-09
    • 2018-06-29
    • 1970-01-01
    • 2011-03-06
    • 2011-01-08
    • 2017-05-01
    • 2013-05-09
    相关资源
    最近更新 更多