title: Vertices of the convex hull of n-dimensional point set
body: <p>I have a given set of points in dimension n. Of these I want to find those, which are the vertices (corners) of the convex hull.
  I want to solve this with Python (but may call other programmes).</p>
  
  <p>Edit: All coordinates are natural numbers. As output I am looking for the indices of the vertices.</p>
  
  <p>Googling usually yielded the problem in 2D, or asked for listing the faces, which is computationally much harder.</p>
  
  <p>My own attempts so far</p>
  
  <ul>
  <li><code>scipy.spatial.ConvexHull()</code>: 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.</li>
  <li><p><code>Normaliz</code> (as part of polymake): works, but too slow. But maybe I did something wrong.</p>
  
  <pre><code>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
  </code></pre></li>
  <li><p>Solve with linear programmes: works, but completely not optimised, so I think there must be a better solution.</p>
  
  <pre><code>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]]
  </code></pre></li>
  </ul>
  
  <p>Example:</p>
  
  <pre><code>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]])
  </code></pre>
  
  <p>yields running time</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>For a slightly larger example, the ratio was worse, so it also can't be explained by the parallelisation.</p>
  
  <p>Any help or improvement is appreciated.</p>