Finding the minimum point in the convex hull of a finite set of points
Based on the work of Philip Wolf [1] and the recursive algorithm of Kazuyuki Sekitani and Yoshitsugu Yamamoto [2].
The algorithm in [2] has 3 epsilon to avoid comparison problems in three parts of the algorithm. The code below has few changes and only one epsilon. The aim of the change is to avoid infinite loops.
Code
1 from numpy import array, matrix, sin, sqrt, dot, cos, ix_, zeros, concatenate, abs, log10, exp, ones
2 from numpy.linalg import norm
3
4 from mpmath import mpf, mp
5 mp.dps=80
6
7 def find_min_point(P):
8 # print "Calling find_min with P: ", P
9
10 if len(P) == 1:
11 return P[0]
12
13 eps = mpf(10)**-40
14
15 P = [array([mpf(i) for i in p]) for p in P]
16
17 # Step 0. Choose a point from C(P)
18 x = P[array([dot(p,p) for p in P]).argmin()]
19
20 while True:
21
22 # Step 1. \alpha_k := min{x_{k-1}^T p | p \in P}
23 p_alpha = P[array([dot(x,p) for p in P]).argmin()]
24
25 if dot(x,x-p_alpha) < eps:
26 return array([float(i) for i in x])
27
28 Pk = [p for p in P if abs(dot(x,p-p_alpha)) < eps]
29
30 # Step 2. P_k := { p | p \in P and x_{k-1}^T p = \alpha_k}
31 P_Pk = [p for p in P if not array([(p == q).all() for q in Pk]).any()]
32
33 if len(Pk) == len(P):
34 return array([float(i) for i in x])
35
36 y = find_min_point(Pk)
37
38
39 p_beta = P_Pk[array([dot(y,p) for p in P_Pk]).argmin()]
40
41 if dot(y,y-p_beta) < eps:
42 return array([float(i) for i in y])
43
44
45 # Step 4.
46 P_aux = [p for p in P_Pk if (dot(y-x,y-p)>eps) and (dot(x,y-p)!=0)]
47 p_lambda = P_aux[array([dot(y,y-p)/dot(x,y-p) for p in P_aux]).argmin()]
48 lam = dot(x,p_lambda-y) / dot(y-x,y-p_lambda)
49
50 x += lam * (y-x)
51
52
53
54 if __name__ == '__main__':
55 print find_min_point( [array([ -4.83907292e+00, 2.22438863e+04, -2.67496763e+04]), array([ 9.71147604, -351.46404195, -292.18064276]), array([ 4.60452808e+00, 1.07020174e+05, -1.25310230e+05]), array([ 2.16080134e+00, 5.12019937e+04, -5.96167833e+04]), array([ 2.65472146e+00, 6.70546443e+04, -7.71619656e+04]), array([ 1.55775358e+00, -1.34347516e+05, 1.53209265e+05]), array([ 13.22464295, 1869.01251292, -2137.61850989])])
56
57
58 print find_min_point( [array([ -4.83907292e+00, 2.22438863e+04, -2.67496763e+04]), array([ 9.71147604, -351.46404195, -292.18064276]), array([ 4.60452808e+00, 1.07020174e+05, -1.25310230e+05]), array([ 2.16080134e+00, 5.12019937e+04, -5.96167833e+04]), array([ 2.65472146e+00, 6.70546443e+04, -7.71619656e+04]), array([ 1.55775358e+00, -1.34347516e+05, 1.53209265e+05]), array([ 13.22464295, 1869.01251292, -2137.61850989]), array([ 12273.18670123, -1233.32015854, 61690.10864825])])