Evaluate the weights of a convex combination

795 Views Asked by At

I'm using the scipy.spatial.ConvexHull API for evaluating the convex hull of a set of points and it works well. Given the code below:

p1 = points[11]
hull = ConvexHull(points)
vertices = [points[i] for i in hull.vertices] 

how can I evaluate the coefficients(weights) of the convex combination of vertices that equals to p1?

Thank you very much, Moshe

1

There are 1 best solutions below

2
Stef On BEST ANSWER

Note: if the number of vertices is higher than d+1, where d is the dimension, then the combination will not be unique. In the remainder of the answer I assume that d=2 for simplicity.

Input:

  • vertices v0 = (x0, y0), v1 = (x1, y1), ..., vn = (xn, yn);
  • a point p1 = (x,y);

Output:

  • a combination a0, a1, ..., an;

such that:

  • x = a0 * x0 + a1 * x1 + ... + an * xn;
  • y = a0 * y0 + a1 * y1 + ... + an * yn;
  • 1 = a0 + a1 + ... + an;
  • forall i, ai >= 0.

This is a system of three linear equations in the unknown (a0, a1, ..., an), along with n+1 linear inequations. We can solve this system using scipy's linear programming module scipy.optimize.linprog

Example solution:

import random                      # generate points
import scipy.spatial as scispa     # ConvexHull
import scipy.optimize as sciopt    # linprog
from scipy.sparse import identity  # identity matrix

points = [(random.randrange(0,100), random.randrange(0,100)) for _ in range(20)]
p1 = points[11]
hull = scispa.ConvexHull(points)
vertices = [points[i] for i in hull.vertices]

c = [1 for _ in vertices]
A = [[x for x,y in vertices], [y for x,y in vertices], [1 for _ in vertices]]
b = [p1[0], p1[1], 1]
s = sciopt.linprog(c, A_eq = A, b_eq = b)

>>> s.x
array([0.13393774, 0.06470577, 0.07367599, 0.09523271, 0.18924727,
       0.26909487, 0.17410566])
>>> p1
(36, 57)
>>> [sum(a*xi for a,(xi,_) in zip(s.x, vertices)), sum(a*yi for a,(_,yi) in zip(s.x, vertices))]
[36.00000000719907, 57.00000000671608]

Important note: I began this answer by warning that if there are more than 3 vertices in the plane, or more than d+1 vertices in dimension d, then there will be more than 1 solution to our system of equations. The code above only returns one solution, meaning that an arbitrary choice was made. You have control over this arbitrary choice: scipy.optimize.linprog is an optimisation tool; here it is minimizing the dot-product of the vector c we gave as parameter, and the solution vector s.x. Here I set all coefficients of c to 1, meaning linprog will find a solution that minimizes the sum of coefficients in the solution; try supplying a different vector c and you'll get a different solution.