How to describe the convex hull of a set of points as an implicit region for optimization?

Question

Last edit: 2016-05-03

A.G.'s answer was the one helped me the most, see explanation at the end of this question.

Original question

In Mathematica you can create out of a set of points in 2D and 3D the convex hull and use it as a region for optimization problem, e.g.,

(*Dimension and number of points*)
d = 2;
np = 10;
(*Generate data and convex hull*)
int = {5, 10};
data = RandomReal[int, {np, d}];
ch = ConvexHullMesh@data;
(*Optimize on convex hull*)
var = Array[v, d];
max = Maximize[var.var, Element[var, ch]];
pic = Show[{ch, 
   Graphics[{Point[data], PointSize -> Large, Red, 
     Point[var /. max[[2]]]}]}, Axes -> True]

You can also performe optimization problems in higher dimensions, e.g.,

d = 7;
reg = Ball[ConstantArray[0, d], 3];
var = Array[v, d];
Maximize[var.var, Element[var, reg]]

(* {9, {v[1] -> -3, v[2] -> 0, v[3] -> 0, v[4] -> 0, v[5] -> 0, v[6] -> 0, v[7] -> 0}} *)

But now suppose we have a set of 5D vectors $\{c_1,c_2,...,c_m\}, c_i \in \mathbb{R}^5$ and we consider the set of all possible convex combinations of them (so we consider the convex hull $C$ of them).

\begin{equation} C = \left\{ v \in \mathbb{R}^5 \ : \ \sum_{i=1}^m w_i c_i = v \ , \quad \sum_{i=1}^m w_i = 1 \ , \quad w_i \geq 0 \ \forall \ i \right\} \end{equation}

I know how to obtain the corners of a given set of points (data), see this question

How to generate higher dimensional convex hull?

But even if we would not know how to, we could still just used the original data. Is there any way to define in Mathematica a region, I dont know, may be with ImplicitRegion, which then can be used for optimization problems as in the 2D and 7D examples above? By that I mean something like the following (not working) structure.

(*Dimension and number of points*)
d = 5;
np = 4;
(*Generate data and region*)
int = {5, 10};
data = RandomReal[int, {np, d}];
reg = ImplicitRegion[
   w1*data[[1, 1]] + w2*data[[2, 1]] + w3*data[[3, 1]] + w4*data[[4, 1]] == x1
    && w1*data[[1, 2]] + w2*data[[2, 2]] + w3*data[[3, 2]] + w4*data[[4, 2]] == x2
    && w1*data[[1, 3]] + w2*data[[2, 3]] + w3*data[[3, 3]] + w4*data[[4, 3]] == x3
    && w1*data[[1, 4]] + w2*data[[2, 4]] + w3*data[[3, 4]] + w4*data[[4, 4]] == x4
    && w1*data[[1, 5]] + w2*data[[2, 5]] + w3*data[[3, 5]] + w4*data[[4, 5]] == x5
    && w1 + w2 + w3 + w4 == 1
    && w1 >= 0 && w2 >= 0 && w3 >= 0 && w4 >= 0
   , {x1, x2, x3, x4, x5}
   ];
(*Optimize on convex hull*)
var = Array[v, d];
max = NMaximize[var.var, Element[var, reg]];

I just dont know how to formulate correctly the region. I can test candidate points cand for given data/corner points corners with the function

TestCandidate[corners_, cand_] := Block[
   {ws, w, eqs, fi},
   ws = Array[w, {Length[corners]}];
   eqs = Sum[ws[[i]]*corners[[i]], {i, Length@corners}] - cand;
   eqs = Flatten@Join[
      Table[eqs[[i]] == 0, {i, Length@eqs}]
      , {Total[ws] == 1}
      , Table[ws[[i]] >= 0, {i, Length@ws}]
      ];
   eqs = And @@ eqs;
   fi = FindInstance[eqs, ws];
   Length@fi != 0
   ];

e.g., I can test that the data can represent itself

And @@ Table[TestCandidate[data, data[[i]]], {i, Length@data}]

True

and also other candidate points cand. But I have not been able to formulate a region with this function. Any ideas on how to formulate the region for given data?

EDIT - 2016-05-03

Thank you very much for the answers. From my point of view, Chip's answer is excellent but A.G.'s answer is the best one, since it made me recognize that the optimization can be carried out a lot easier over the simplex he pointed out. You can treat then any problem for given data following his recommendation as follows

(*Dimension and number of points*)
d = 9;
np = 10^2;
(*Generate data*)
data = RandomReal[{1, 2}, {np, d}];
(*Weights and region*)
ws = Array[w, Length[data]];
reg = Simplex[Table[UnitVector[np, k], {k, np}]];
(*Function and optimization*)
v = Sum[ws[[i]]*data[[i]], {i, Length@data}];
f = Sin[v.v] + 1/v.v;
NMinValue[{f, Total[ws] == 1, Table[0 <= ws[[i]], {i, Length@ws}]}, 
   ws] // Quiet // AbsoluteTiming
NMinValue[f, Element[ws, reg]] // Quiet // AbsoluteTiming

{263.576, -0.942131}

{270.522, -0.942131}

The main problem is not the dimensionality of the vector variables v, as he pointed out, but the size of the data, and therefore, the number of ws. As you can see, you can either use the conditions for the ws directly or use the Simplex (see also Wikipedia), which is equivalent in this case. Both seem to be equally fast (around 270 seconds for the created data on 4 cores). If you want to reduce your data to the corners of the convex hull, see again the linked question

How to generate higher dimensional convex hull?

Answer 1

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\}$$

Answer 2

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}}

Answer 3

One way to treat this in some very special, low dimensional and friendly cases is to use ParametricRegion

(*Dimension and number of points*)
d = 5;
np = 4;
(*Generate data*)
data = RandomInteger[{-10, 10}, {np, d}];
(*Convex hull*)
ws = Array[w, Length[data]];
reg3 = ParametricRegion[
   {Sum[ws[[i]]*data[[i]], {i, Length@data}], Total[ws] == 1}
   , Evaluate@Table[{ws[[i]], 0, 1}, {i, Length@ws}]
   ];
(*Test data*)
Table[Element[data[[i]], reg3], {i, Length@data}]
(*Optimization in convex hull*)
v = Array[vc, d];
Minimize[v.v, Element[v, reg3]]

{True, True, True, True}

{27763033/486654, {vc[1] -> -(86579/243327), vc[2] -> -(57359/34761), vc[3] -> 632741/162218, vc[4] -> -(181833/81109), vc[5] -> -(2835953/486654)}}

Naturally, you will obtain other results, due to the random number generator. But it works for this friendly example. For higher dimensional problems with a lot of real valued corner points, the optimization is either reaaaaaally slow or just returns nothing.