3

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.

Stefan Falk
  • 23,898
  • 50
  • 191
  • 378
Hennich
  • 682
  • 3
  • 18
  • The question is, do you need more dimensions then `3`, is this or a CG Programm? How accurate do you want to be? – user1767754 Nov 23 '17 at 19:51
  • Yes, I need arbitrarily large dimensions. The problem comes from polynomial optimisation. My coordinates are all natural numbers and the answer is a list of indices, so an exact algorithm would be nice. – Hennich Nov 23 '17 at 20:23
  • It's going to be computationally complex, https://cw.fel.cvut.cz/wiki/_media/misc/projects/oppa_oi_english/courses/ae4m39vg/lectures/05-convexhull-3d.pdf if you look here, only quickhull behvaes kinda "fast" but for multi dimensions this can be tricky. – user1767754 Nov 23 '17 at 20:25
  • @user1767754 This talk, as far as I see, want to find the facets, which I don"t. My problem can be solved in polynomial time via the LPs above. – Hennich Nov 23 '17 at 20:28

1 Answers1

3

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.

sebigu
  • 164
  • 5