# How to generate higher dimensional convex hull?

I like the function ConvexHullMesh very much, since I often need to take a look at the convex hull of points in 2D and 3D, which as an example can be generated and tested with the following code

(*Dimension and number of points*)
d = 2;
np = 10^1;
(*Generate data*)
data = RandomReal[{-1, 1}, {np, d}];
(*Construct convex hull*)
ch = ConvexHullMesh@data;
(*Test new point*)
new = RandomReal[{-1, 1}, d];
Element[new, ch]


True

But just out of fun, I wanted to do this with higher dimensional points, e.g., with 5D points. But how do you construct in Mathematica the convex hull then out of the data and create a set for testing new points?

• Let's restrict ourselves to the 4D case, for simplicity. How do you expect the points making up the corners of your convex hull to be arranged? Apr 24, 2016 at 17:29
• Arranged? I dont see why they would have to have a specific order. Or what do you mean by arrangement? I think I dont understand. As far as I know, the convex hull is the set of all possible convex combinations of the data points, no specific order is imposed. Apr 24, 2016 at 17:31
• Consider for instance the 3D convex hull. The points are organized (alright, that may have been the more appropriate word to use) so that there is a polygon corresponding to each facet of the hull. In the 4D case, you have cells; how should the points be grouped in those individual cells? Or will you be satisfied with a routine that returns just the coordinates of the corners? Apr 24, 2016 at 17:35
• A routine with the corners would be fine for me, since then any point inside can be constructued as a convex combination. Apr 24, 2016 at 17:44
• The best algorithm I know to find the convex hull in R^n is the Quickhull. There is software, Qhull, that uses this algorithm for finding the convex hull and with the output format you can, by combinatorial methods, find the facets. Apr 24, 2016 at 20:24

As others mentioned, Qhull can do this. There are multiple ways to access Qhull from Mathematica.

One way is through the mPower package, now part of xCellerator. An answer discussing how to install the package is here: https://mathematica.stackexchange.com/a/18909/12

Another way is the qh-math package, described in this answer: https://mathematica.stackexchange.com/a/13445/12

Yet another way is through MATLAB and MATLink (MATLAB has Qhull built in).

<<MATLink
OpenMATLAB[]
convhulln = MFunction["convhulln"]/*Round

pts = RandomReal[1, {20,5}] (* 20 points in 5D *)

convhulln[pts]


This returns a list of facets encoded as point-index-lists.

• Cool, thank you. I dont have MATLAB on my laptop now, but I will take a look at the references you gave. Do you happen to know if the generated convex hull object is also usable for later optimization problems like Maximize[f,Element[{x1,x2,x3,x4,x5},conhull]]? This is the case for the Mathematica ConvexHullMesh which is kind of a nice property. Apr 26, 2016 at 15:47
• @MauricioLobos Well, it just gives you the facets, not anything Mathematica-specific. BoundaryMeshRegion and MeshRegion only support 1D, 2D and 3D regions... I'm not sure how to handle such a thing as a proper region. Maybe instead of that you could convert the constraint to a set of inequalities? Sorry, I don't know enough about what works well for this type of optimization problem ... Apr 26, 2016 at 15:53
• Ok, thanks, I will have to take a deeper look into it. Apr 26, 2016 at 15:54

New implementation (EDIT:2017-May-08) I just wanted to share a slightly different implementation (based on the old one, see below). The following function (almost the same as TestCandidate from the old implementation) tests if a candidate point cand can be described by a convex combination of the given corner points corners

ConvexDepenentQ[corners_, cand_] := Module[
{w, ws},
w = Array[ws, Length@corners];
1 == Length@
FindInstance[
w.corners == cand && Total[w] == 1 &&
And @@ Table[w[[i]] >= 0, {i, Length@w}], w]
];


The following function reduces a given data set data to its convex corners, which describe the convex hull. This is done by picking one data point out and testing if it is convex dependent of the rest.

ConvexReduce[data_] := Module[
{corners, ncorners, test},
corners = data;
Do[
ncorners = Delete[corners, Position[corners, data[[i]]]];
test = ConvexDepenentQ[ncorners, data[[i]]];
If[test, corners = ncorners];
, {i, Length@data}
];
corners
];


Naturally, you could get the convex corners just with ConvexReduce but we can speed this up a little by starting with the data points on the bounding rectangle (as in the old implementation), going through the data and only append points which are not convex dependent. At the end we can convex reduce the smaller set.

ConvexHullCorners[data_] := Module[
{corners, rd},
corners = {};
Do[
corners =
Join[corners,
Select[data,
Min[data[[;; , i]]] == #[[i]] ||
Max[data[[;; , i]]] == #[[i]] &]];
, {i, Length@data[[1]]}
];
corners = DeleteDuplicates@corners;
rd = Delete[data, First@Position[data, #] & /@ corners];
Do[
If[ConvexDepenentQ[corners, rd[[i]]], , AppendTo[corners, rd[[i]]]]
, {i, Length@rd}
];
ConvexReduce@DeleteDuplicates@corners
];


At least compared to the algorithms implemented by Loren Petrich, see http://lpetrich.org/Science/ and Jason's answer in this question, my new implementation seems to perform ok. But please, keep in mind that I obly look for the convex corners of the data (this is the easy part), while qhull also delivers the highdimensional facettes of the data (which is the veeeery hard and time expensive part, which I at least dont need). For increasing number of points, this seems to be ok.

d = 6;
np = 40;
ts = {};
Do[
data = RandomReal[{-10, 10}, {i, d}];
AppendTo[ts, {
AbsoluteTiming[ConvexHullCorners@data;][[1]]
, AbsoluteTiming[CHNChanShattering@data;][[1]]
, AbsoluteTiming[CHNIncremental@data;][[1]]
, AbsoluteTiming[CHNQuickHull@data;][[1]]
}];
, {i, 10, np, 10}
]
ListPlot[Transpose@ts, Joined -> True, PlotLegends -> {
"ConvexHullCorners"
, "CHNChan"
, "CHNInc"
, "CHNQuickH"
}]


But specially for increasing dimension, the bounding rectangle approach seems to be a good start.

d = 7;
np = 30;
ts = {};
Do[
data = RandomReal[{-10, 10}, {np, i}];
AppendTo[ts, {
AbsoluteTiming[ConvexHullCorners@data;][[1]]
, AbsoluteTiming[CHNChanShattering@data;][[1]]
, AbsoluteTiming[CHNIncremental@data;][[1]]
, AbsoluteTiming[CHNQuickHull@data;][[1]]
}];
, {i, 2, d}
]
ListPlot[Transpose@ts, Joined -> True, PlotLegends -> {
"ConvexHullCorners"
, "CHNChan"
, "CHNInc"
, "CHNQuickH"
}]


If you want to check, that the corners delivered by the algorithms are identical, use

Sort@ConvexHullCorners@data == Sort@data[[CHNQuickHull[data][[1]]]]


True

If you want to check, that all data points are convex dependent of the corners, use

And @@ (ConvexDepenentQ[ConvexHullCorners@data, #] & /@ data)


True

I am running these computations in Windows 10, Mathematica 11.0.1. Sadly, I was never able to connect to Matlab, so I dont know, if my implementation is slower (probably). Nevertheless, this might help somebody.

Old implementation

Thanks @physicien for the comment on Quickull and Qhull, see his comment on the question. Sadly, I did not understand the paper and so I wanted to give this a shot in Mathematica with some very rude implementation (sorry, I am extremly bad at programming). There is a lot of room for improvement (my solution is extremly slow in high dimensions and for large data), but may be this help also others like me, at least for a start.

For a candidate vector $v \in \mathbb{R}^n$ (cand) I will first test if the candidate is represtable by a convex combination of a given set of coners $\{c_1,c_2,...,c_m\}, c_i \in \mathbb{R}^n$ (corners)

$$\sum_{i=1}^m w_i c_i = v \ , \qquad \sum_{i=1}^m w_i = 1 \ , \qquad w_i \geq 0 \ \forall \ i$$

with the following 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
];


This functions returns True if the candidate is representable, since if so, an instance solving this can be found, False otherwise, i.e., the candidate is not represtable as a convex combination.

Now I will search the set of corners, which is able to represent any point of a given data set (data) by convex combinations as follows:

1. Pick initial corner points from data belonging to a $n$-dimensional estimate for a bounding box

2. Delete the picked points from 1. from the data

3. Check if any points picked in 1. are representable by a smaller set (this happens sometimes if all points are actually in a $(n-1)$-dimensional plane.

4. Run through the rest of the data points and add a candidate if it is NOT representable by the updated corners. If a new corner is added, check all subsets of corners in case a new corner allows for the representation of old corners and reduce momentary set of corners.

Code

ConvexHullCorners[data_] := Block[
{datatemp, dim, center, low, pos, upp, corners, cand, cornerstemp},
(***********************)
(*1*)
(*Generate initial bounding rectangle estimate points*)
dim = Length[data[[1]]];
center = Mean[data];
corners = {};
Do[
low = Max[center[[i]] - data[[;; , i]]];
pos = Position[data[[;; , i]], center[[i]] - low][[1, 1]];
AppendTo[corners, data[[pos]]];
upp = Max[data[[;; , i]] - center[[i]]];
pos = Position[data[[;; , i]], center[[i]] + upp][[1, 1]];
AppendTo[corners, data[[pos]]];
, {i, dim}
];
corners = DeleteDuplicates@corners;

(***********************)
(*2*)
(*Delete points from bounding box estimate from data and \
duplicates*)
pos = Flatten[
Table[Position[data, corners[[i]]], {i, Length@corners}], 1];
datatemp = DeleteDuplicates@Delete[data, pos];

(***********************)
(*3*)
(*Check for representable points,
and delete them from corners if so*)
cornerstemp = corners;
Do[
cand = corners[[j]];
pos = Position[cornerstemp, cand][[1]];
cornerstemp = Delete[cornerstemp, pos];
If[TestCandidate[cornerstemp, cand] == False,
PrependTo[cornerstemp, cand]];
, {j, Length[corners]}
];
corners = cornerstemp;

(***********************)
(*4*)
(*Find corners of convex hull*)
Do[
(*Test candidate from rest of data*)
cand = datatemp[[i]];
If[
(*If not representable, add to corners*)
TestCandidate[corners, cand] == False,
AppendTo[corners, cand];
(*Check if just added corner allows for representation of some \
former old corners, and delete them from corners if so*)
cornerstemp = corners;
Do[
cand = corners[[j]];
pos = Position[cornerstemp, cand][[1]];
cornerstemp = Delete[cornerstemp, pos];
If[TestCandidate[cornerstemp, cand] == False,
PrependTo[cornerstemp, cand]];
, {j, Length[corners] - 1}
];
corners = cornerstemp;
];
, {i, Length@datatemp}
];

(***********************)
(*OUT*)
(*Return corners of convex hull*)
corners
];


Some examples:

1) 2D example

Generate some random points in 2D (red points), compute the normal convex hull with Mathematica (blue region), compute the corners with the program given above (green points) and visualize results.

(*Dimension and number of points*)
d = 2;
np = 10^1*2;
(*Generate data*)
data = RandomReal[{-10, 10}, {np, d}];
ch = ConvexHullMesh@data;
datac = ConvexHullCorners@data;
pic = Show[
ch,
Graphics[{PointSize -> 0.05, Green, Point[datac],
Directive[Red, Opacity[0.4]], PointSize -> Large, Point[data]},
ImageSize -> Small]
]


2) 6D example:

Analogously, generate data and extract corners of convex hull with program. Check if all points of the data are really representable by the given corners of the program and check a false point.

(*Dimension and number of points*)
d = 6;
np = 10^1*3;
(*Generate data*)
int = {5, 10};
dis = 7;
data = RandomInteger[int, {np, d}];
(*Extract corner of convex hull*)
datac = ConvexHullCorners@data;
(*Test original data and point outside of convex hull*)
And @@ Table[TestCandidate[datac, data[[i]]], {i, Length@data}]
TestCandidate[datac, RandomReal[int + dis, d]]


True

False

Loren Petrich has done a wonderful job of implementing quite a few convex hull algorithms in Mathematica, and you can get all of his original code from his website. Since it is written under the MIT license I took the liberty of wrapping the code relevant to an N-dimensional convex hull into a package.

Included are the functions CHNGiftWrapping, CHNChanShattering, CHNQuickHull, and CHNIncremental, which use the gift wrapping algorithm, Chan's algorithm, quick hull, and the incremental algorithms, respectively.

All these functions work the same, they take a list of N-dimensional vertices and return a sorted list of vertices which make up the hull and a list of simplexes.

(* load the package *)
<<"https://gist.githubusercontent.com/jasondbiggs/c3d9410af3195da514a442be5b563ab8/raw/80ac5074077f0d4b1366540c8c710e35cb530ddd/NDConvexHull.m"


Here I apply the quick hull function to a list of 4-dimensional points:

SeedRandom[4];

list = RandomReal[5, {50, 4}];

CHNQuickHull[list]
(* {{2, 3, 7, 8, 11, 12, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
26, 27, 28, 29, 31, 32, 33, 34, 35, 36, 37, 38, 40, 41, 42, 45, 46,
47, 48, 49,
50}, {{2, 15, 21, 47}, {2, 15, 47, 50}, {2, 15, 50, 21}, {2, 21, 35,
47}, {2, 21, 50, 35}, {2, 22, 35, 50}, {2, 22, 47, 35}, {2, 22,
50, 47}, {3, 16, 17, 22}, {3, 16, 22, 49}, {3, 16, 49, 17}, {3, 17,
49, 22}, {7, 22, 29, 35}, {7, 22, 34, 49}, {7, 22, 35, 34}, {7,
22, 36, 29}, {7, 22, 49, 36}, {7, 29, 34, 35}, {7, 29, 36, 49}, {7,
29, 49, 34}, {8, 16, 17, 49}, {8, 16, 18, 24}, {8, 16, 24,
17}, {8, 16, 49, 18}, {8, 17, 20, 32}, {8, 17, 24, 20}, {8, 17, 32,
49}, {8, 18, 20, 24}, {8, 18, 23, 20}, {8, 18, 29, 37}, {8, 18,
36, 23}, {8, 18, 37, 36}, {8, 18, 49, 29}, {8, 20, 23, 33}, {8, 20,
33, 41}, {8, 20, 41, 32}, {8, 23, 36, 33}, {8, 29, 49, 37}, {8,
32, 33, 36}, {8, 32, 36, 49}, {8, 32, 41, 33}, {8, 36, 37,
49}, {11, 12, 15, 23}, {11, 12, 17, 40}, {11, 12, 18, 42}, {11, 12,
21, 45}, {11, 12, 23, 18}, {11, 12, 40, 21}, {11, 12, 42,
17}, {11, 12, 45, 15}, {11, 15, 18, 23}, {11, 15, 21, 26}, {11, 15,
26, 27}, {11, 15, 27, 18}, {11, 15, 45, 21}, {11, 16, 17,
24}, {11, 16, 19, 28}, {11, 16, 21, 40}, {11, 16, 24, 19}, {11, 16,
28, 21}, {11, 16, 40, 17}, {11, 17, 42, 24}, {11, 18, 19,
24}, {11, 18, 24, 42}, {11, 18, 27, 29}, {11, 18, 29, 19}, {11, 19,
29, 46}, {11, 19, 35, 28}, {11, 19, 46, 35}, {11, 21, 28,
35}, {11, 21, 35, 38}, {11, 21, 38, 26}, {11, 26, 38, 27}, {11, 27,
38, 46}, {11, 27, 46, 29}, {11, 35, 46, 38}, {12, 15, 23,
45}, {12, 17, 20, 42}, {12, 17, 40, 45}, {12, 17, 45, 20}, {12, 18,
42, 23}, {12, 20, 23, 42}, {12, 20, 45, 23}, {12, 21, 45,
40}, {15, 18, 23, 27}, {15, 21, 26, 50}, {15, 21, 47, 45}, {15, 23,
36, 27}, {15, 23, 45, 48}, {15, 23, 48, 36}, {15, 26, 27,
38}, {15, 26, 38, 50}, {15, 27, 36, 50}, {15, 27, 50, 38}, {15, 31,
36, 47}, {15, 31, 47, 50}, {15, 31, 50, 36}, {15, 36, 48,
47}, {15, 45, 47, 48}, {16, 17, 21, 40}, {16, 17, 22, 47}, {16, 17,
47, 21}, {16, 18, 19, 29}, {16, 18, 24, 19}, {16, 18, 29,
49}, {16, 19, 21, 35}, {16, 19, 28, 21}, {16, 19, 34, 29}, {16, 19,
35, 34}, {16, 21, 47, 35}, {16, 22, 34, 35}, {16, 22, 35,
47}, {16, 22, 49, 34}, {16, 29, 34, 49}, {17, 20, 32, 47}, {17, 20,
42, 24}, {17, 20, 47, 45}, {17, 21, 40, 47}, {17, 22, 32,
49}, {17, 22, 47, 32}, {17, 40, 45, 47}, {18, 20, 24, 42}, {18, 20,
42, 23}, {18, 23, 27, 37}, {18, 23, 37, 36}, {18, 27, 29,
37}, {19, 21, 35, 28}, {19, 29, 35, 34}, {19, 29, 46, 35}, {20, 23,
33, 41}, {20, 23, 41, 45}, {20, 32, 47, 41}, {20, 41, 47,
45}, {21, 26, 50, 38}, {21, 35, 38, 50}, {21, 40, 47, 45}, {22, 29,
35, 36}, {22, 31, 36, 50}, {22, 31, 47, 36}, {22, 31, 50,
47}, {22, 32, 36, 47}, {22, 32, 49, 36}, {22, 35, 50, 36}, {23, 27,
37, 36}, {23, 33, 41, 36}, {23, 36, 41, 48}, {23, 41, 45,
48}, {27, 29, 36, 38}, {27, 29, 37, 36}, {27, 29, 38, 46}, {27, 36,
50, 38}, {29, 35, 36, 38}, {29, 35, 38, 46}, {29, 36, 49,
37}, {32, 33, 36, 41}, {32, 36, 47, 48}, {32, 36, 48, 41}, {32, 41,
48, 47}, {35, 36, 38, 50}, {41, 45, 48, 47}}} *)


You can verify that they all produce the same results:

CHNAlgorithms = {CHNGiftWrapping, CHNChanShattering,
CHNQuickHull, CHNIncremental};

SameQ[Through[CHNAlgorithms][list]]
(* True *)

• Wow, cool! Thank you very much! Just one question, to check I understand the output. The first entry of the output is the list of positions of the corners in the original data, right? And the second entry contains the positions in the original data of the respective simplex, right? May 27, 2016 at 14:50
• Yes! The code he has for 3D and 2D gives answers that match those from ConvexHullMesh, but I couldn't check the higher dimensions. I'm looking for a contact email on his site to let him know I posted it here, but I don't see it May 27, 2016 at 14:52
• Yes, it works exactly like that, e.g., a 6D-problem list = RandomReal[5, {50, 6}]; chl = CHNQuickHull[list]; chp = list[[chl[[1]]]]; And @@ Table[TestCandidate[chp, list[[i]]], {i, Length@list}] returns True with my function from above, i.e., all points of the 6D-list are representable with the points addressed in the output of CHNQuickHull. Thank you! May 27, 2016 at 15:44
• Just as a remark, the algorithms work only for dimensions higher than 2. data = RandomReal[2, {10, 2}]; CHNQuickHull[data] returns empty lists, i.e., no results. But, of course, for 2D you can use Mathematica's ConvexHullMesh. May 29, 2016 at 7:20
• @MauricioLobos Yeah, I was surprised that the CHN` functions work in 3D as well, since the comments made it seem like it would only work for over 3 dimensions. If you go to the original source there are 2D and 3D versions of all the algorithms and even some spherical coordinate versions. May 29, 2016 at 9:58