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).
Edit: All coordinates are natural numbers. As output I am looking for the indices of the vertices.
Googling usually yielded the problem in 2D, or asked for listing the faces, which is computationally much harder.
My own attempts so far
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]]
Example:
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]])
yields running time
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)
For a slightly larger example, the ratio was worse, so it also can't be explained by the parallelisation.
Any help or improvement is appreciated.
You can change your is_in_convex_hull method in the following way:
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
The Normaliz version of the algorithm will run much faster then.