Nearest point in a convex hull

Specific Domains Optimization (Mathematical)
1

Background:

I’m using ScatteredInterpolation to approximate a higher-order value function as needed.

The interpolation works reasonably well for points that are in some sense “inside” prior points that I’ve used. My plan is to use LazySets to create a convex hull of existing interpolation points. Then when a new request comes in, I check to see if it’s interior to the convex hull. If it is not, I want to add it.

However, I think the best approach will not be to add the just-request point. Instead, I want to move some distance “away” from the hull, so that if the next call is epsilon in the same direction, it will be in the interior. I think the way to calculate the direction “away” is to find the nearest point in the hull, and then move some distance in the opposite direction.

Question:

Given a point in R^N, and a convex hull of R^N. Assume, the point is not in the hull. How do I find closes point that is in the hull. Some potential solutions:

  1. I think this is isomorphic to point in the minkowski_difference that has the smallest norm. But I don’t see how to find the point with the smallest norm.

  2. Sample 100 points from the convex hull and then use the smallest.

  3. I also wonder if there’s a way to iterate through all violated constraints, moving in the direction of the a_vector.

  4. Don’t worry about finding the closest point, just find the closest vertex, as if there’s more than one violation, the closet point is going to be a vertex. (This is false. See image at end of post. The Green Point violates 3 constraints, but the nearest point is the blue point which is not a vertex.)

3+4: check the number of violated constraints. If it’s 1, that constraint’s vector tells us the direction. If it’s 2 or more, just iterate over vertices.

Screenshot 2023-02-18 at 18.31.12
Screenshot 2023-02-18 at 18.31.12874×702 51.7 KB

1 Like
2

A similar question is addressed and answered on mathematics.stackexchange:

3

That is interesting. From the linked-to FAQ. It says, “For Delaunay triangulations” you can call the specific function. Are convex hulls equivalent to Delaunay triangulations?

This comment on the answer says that the function is not accessible in scipy, so presumably it’s not accessible in the Julia interface on top of scipy.

4

If you have a finite set of 2d points, P_k, k=1:n, you can define its convex hull, which a convex polygon. Delaunay triangulation associated to this set of points is a family of triangles having these points as vertices, and such that the open disk bounded by the circumcircle of a triangle does not contain any other point P.
In 3d the convex hull is a volume, and the Delaunay triangulation consists in tetrahedra with a similar property, that the open ball bounded by the circumsphere of a tetrahedron does not contain any point P from the initail set of points, P_k.
The convex hull of a set of points is the smallest convex set containing these points, while the Delaunay triangulation is a triangulation/tetrahedralization of the convex hull. Points, segments of lines, triangles, tetrahedra in dimension 3 are all simplices, and the Delaunay triangulation is a simplicial complex Simplicial complex - Wikipedia with the property outlined above.

LE: Triangulated convex hull:
delaunay-tr-cvxh

1 Like
5

Yes, in the sense that, for a given set of points \{ x_i : i = 1, \cdots, n\} of dimension d, one can first lift to the paraboloid surface graph \{ (x_i, || x_i ||_2^2) : i = 1, \cdots, n\} in dimension d+1 (one greater), find the convex hull of this set in the lifted space, then project the edges of the convex hull back to dimension d to obtain the Delaunay triangulation. Conversely, if you have a method of computing the Delaunay triangulation, you can find the convex hull of points on a paraboloid surface.

image

See also d-dimensional Delaunay on the Wikipedia page.

Note also that QHull.jl wraps scipy.spatial.ConvexHull.

6

Here is a different approach based on the efficient computation of the minimal volume ellipsoid approximation, an example of which I recently put together as a JuMP tutorial. This could be put to work if you’re willing to accept a more relaxed acceptance criterion of “some distance away” for new points, namely that the candidate point must lie outside the minimal volume ellipsoid of the current set. The advantage of this is that it is very simple to test whether a point satisfies being in the interior of an ellipsoid. Be interested to see if this can be modified to your problem in more specific ways, like only considering certain subsets of incumbent points to given tighter criteria.

7

Thank you everyone for the suggestions, and especially to Christian Schilling for recalling ConvexBodyProximityQueries.jl.

You are right that we can treat this as a convex optimization problem. Also there is a library that can solve this when I pass in the new point as a Singleton, and a guess at a direction (I suspect the difference between any vertex and the point would work.)

I have not compared the speeds, but considering ConvexBodyProximityQueries performs 0 allocations, I would expect it to be faster.

using Revise
using LazySets
using Polyhedra
using LinearAlgebra
import Convex
using SCS
using ConvexBodyProximityQueries
using StaticArrays

subN = float.([
    [0, 0, 0, 0],
    [0, 1, 0, 0],
    [1, 0, 0, 0],  
    [1, 1, 0, 0], 
    [.7, .2, 0, 0],
])

# will drop the final point
h4 = convex_hull(subN)
@show length(h4)

V = VPolytope(h4)

test = SA[.5, -.2, 0, 0]


# Using Convex
y = Convex.Variable(4)

problem = Convex.minimize(dot(y, y) - 2*dot(y, test))

for c in constraints(V)
    problem.constraints += dot(y, c.a) <= c.b
end

Convex.solve!(problem, SCS.Optimizer; silent_solver=true)

problem.status
problem.optval

@show Convex.evaluate(y)



# Using ConvexBodyProximityQueries
closest_points(Singleton(test), V, -test)

Outputs

length(h4) = 4
Convex.evaluate(y) = [0.49999688179102586, -0.00011282533880594664, 0.0, 0.0]
([0.5, -0.2, 0.0, 0.0], [0.5, 0.0, 0.0, 0.0])
2 Likes