How to describe the convex hull of a set of points as an implicit region for optimization?
For 2D, just find the Polygon
representing the convex hull and use RegionMember
:
(* fake data *)
rand = Round[RandomReal[{0, 1}, {10, 2}], 1/100];
prims = MeshPrimitives[ConvexHullMesh[rand], 2][[1, 1]];
Refine[RegionMember[Polygon[Round[prims, 1/100]], {x, y}], {x, y} ∈ Reals]
1/25 (1/20 - x) + 18/25 (-1/50 + y) >= 0 && -23/25 (-77/100 + x) + 3/20 (-3/50 + y) >= 0 && 4/25 (-23/25 + x) - 87/100 (-49/50 + y) >= 0 && 4/5 (-1/20 + x) >= 0
In general, a convex hull in nD can be thought of as the intersection of half spaces that are generated by the (n-1)st dimensional faces of the hull.
In 3D we can compute the convex hull and define the half spaces, using the inequality in the documentation for HalfSpace
.
(* fake data *)
rand = Round[RandomReal[{0, 1}, {10, 3}], 1/100];
(* faces of the convex hull's boundary *)
prims = Round[
MeshPrimitives[RegionBoundary[ConvexHullMesh[rand]], 2][[All, 1]],
1/100
];
(* intersect half spaces *)
res = And @@ (Cross[#2 - #1, #3 - #1].({x, y, z} - #1) <= 0 & @@@ prims)
Let's plot this implicit description to make sure it's correct:
{
RegionPlot3D[res, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, PlotPoints -> 60,
Axes -> False, Boxed -> False, Mesh -> None],
ConvexHullMesh[rand]
}
In higher dimensions, if we can find the boundary of the convex hull, we can use the exact same idea.
To compute the convex hull in higher dimensions, I will use brute force and exploit this fact:
A face is on the boundary of the hull if and only if all points in the given set lie in its half space.
This method scales poorly, but for small enough data it works fine.
First write a function to determine if a given face is on the boundary of the hull:
hullBoundaryQ[data_, pts_] :=
Block[{p1 = First[pts], n},
n = Cross @@ Transpose[Transpose[Rest[pts]] - p1];
Abs[Subtract @@ MinMax[Sign[n.(# - p1)& /@ data]]] < 2
]
Now create data and find the convex hull. I'll work in 4D here:
(* fake data *)
rand = Round[RandomReal[{0, 1}, {10, 4}], 1/100];
(* all possible faces to test *)
candidatefaces = Pick[#, UnsameQ @@@ #]&[
DeleteDuplicates[Sort /@ Tuples[rand, Length[First[rand]]]]
];
(* find all faces on the hull *)
hullfaces = Select[candidatefaces, hullBoundaryQ[rand, #] &];
Now that we have the faces, find the implicit description of each and we're done:
halfSpacenD[pts_, test_, vars_] :=
Block[{p, n, sgn},
p = First[pts];
n = Cross @@ Transpose[Transpose[Rest[pts]] - p];
sgn = Sign[n.(p - test)];
sgn n.(vars - p) <= 0
]
Here's the equation for the 4D convex hull of our data:
vars = {x, y, z, w};
testpoints = First[Complement[rand, #]] & /@ hullfaces;
res = And @@ MapThread[halfSpacenD[##, vars] &, {hullfaces, testpoints}]
Finally, it seems we can optimize over this. For instance:
NMaximize[{vars.vars, res && 0 < x < 1 && 0 < y < 1 && 0 < z < 1 && 0 < w < 1}, vars]
{2.7231, {x -> 0.59, y -> 0.99, z -> 0.85, w -> 0.82}}
Finding the convex hull of points in $\Re^d$ and expressing it as a set of (in)equalities is hard. However, I would suggest you transform the problem by writing feasible points as convex combinations of the given points, i.e.
$$x=\sum_{i=1}^{d} w(i)\, x(i)$$
and then optimize over the simplex
$$\left\{ 0\leq w(i) \leq 1, \sum_{i=1}^{d} w(i) = 1\right\}$$
Implementation should be rather straightforward. Works in any number of dimensions.
BTW, if the sole objective function you want to maximize is the distance to some given point (the origin in you example) then the solution is just... one of the points that generate the convex hull. In that case all that is needed is
$$\max\{ \|x(i)\|_2\ ,\ 1\leq i\leq d\}$$