2

I am trying to find the euclidian distance between a d-dimensional polytope defined as A x <= b and a point, called point in the code below, that I am certain is outside of the polytope. This is a convex quadratic minimization problem. I would like to use cvxpy for this purpose since I have this dependency for the project anyways. Below, I have a 2 dimensional polytope (the third dimension is actually constant) and a point, and the output of cvxpy confuses me greatly.

Questions: why is the minimum found by cvxpy smaller than the manually computed 2-norm? We can clearly see that the point found by cvxpy is not the one that minimizes the Euclidian distance to the polytope, why is this? If this is the wrong approach for finding the minimum distance, how should I instead go about it?

NB pypoman is a package for finding the vertices of the polytope, which is convenient for plotting.

import numpy as np
import cvxpy as cp
import pypoman
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull, convex_hull_plot_2d

A = np.array([
    [ 0.,  0., -1.],
    [ 0.,  0.,  1.],
    [ 0., -1.,  0.],
    [ 0.,  1.,  0.],
    [ 0.,  0., -1.],
    [ 0.,  0.,  1.],
    [ 0.,  3., -1.],
    [ 0., -3.,  1.],
    [-1.,  0.,  0.],
    [ 1.,  0.,  0.],
    [ 1.,  0., -1.],
    [-1.,  0.,  1.],
    [-1.,  0.,  0.],
    [ 1.,  0.,  0.],
    [ 1.,  3., -1.],
    [-1., -3.,  1.]
])
b = np.array([-10.,  10.,   0.,   2.,   0., 100.,   0., 100.,   0., 100.,   0., 100.,   0., 100.,   0., 100.])

vertices = pypoman.duality.compute_polytope_vertices(A, b)
vertices = np.array(vertices)

x_cp = cp.Variable(A.shape[1])
constraints = [
    A @ x_cp <= b,
]
point = cp.Constant([8.0, 2.0, 10.0])
objective = cp.Minimize(cp.norm(x_cp - point, 2))
problem = cp.Problem(objective=objective, constraints=constraints)

print(f'cvxpy minimum: {problem.solve()}')
print(f'point coordinates: {x_cp.value}, CLEARLY NOT THE CLOSEST POINT!')
print(f'distance: {sum((point.value - x_cp.value)**2)}, DIFFERENT FROM CVXPY MINIMUM???')

hull = ConvexHull(vertices[:, :2])
_ = convex_hull_plot_2d(hull)
plt.scatter(*x_cp.value[:2], c='r', label='cvxpy closest')
plt.scatter(*point.value[:2], c='g', label='point')
plt.legend(loc=3, framealpha=1.0)
plt.show()

enter image description here

Patrickens
  • 321
  • 3
  • 14

0 Answers0