# Speed of ConvexHullMesh

Pick $N$ random points on the surface of the sphere, then use ConvexHullMesh[] to compute their convex hull. Empirically, this takes time quadratic time in $N,$ whereas the state of the art as of 1996 is $N \log N.$ Is Mathematica doing something extra, or is its algorithm just 18 years out of date?

In the interest of "politeness".

randpt[n_]:= Module[{prept = RandomVariate[NormalDistribution[], 3]}, prept/Norm[prept]]
experiment[n_]:= ConvexHullMesh[Table[randpt, {n}]//AbsoluteTiming


EDIT Here is the experimental data (the first element in each pair is the number of points, the second the AbsoluteTiming [this on a very fast Windows PC])

{{10000, 6.133357},
{20000, 17.791034},
{30000, 50.460951},
{40000, 106.783207},
{50000, 182.384594}}


This looks awfully quadratic to me. I point out that experiments for small $n$ are meaningless, since there is linear overhead, which has a none-too-small constant, so swamps the quadratic growth. $50000$ points is not that many (and certainly would be a modest number for any FEM mesh), so this is clearly a major performance bug.

ANOTHER EDIT for comparison purposes, I also installed qhull on the very same machine, and ran qconvex on 50000 cospherical random points. The running time was, umm, well, it started spewing output at me the moment I hit return, so I am seeing a roughly three order of magnitude speedup.

• It would be polite to include your code in the question so that it's easier for other people to attempt to reproduce your observations. I tested on my end (v10, OS X) for $N=2, 4, 8, \ldots, 2^{14}, 2^{15}$, and it seems to behave subquadratically for up to $2^{13}\approx 16{,}000$ points, after which it becomes quadratic. Maybe it has an upper limit on working memory for the $N\log N$ algorithm, beyond which it switches to a quadratic time algorithm. – Rahul Sep 12 '14 at 3:44
• @RahulNarain I am not sure how this would be "polite", since the code is completely trivial, but ok. – Igor Rivin Sep 12 '14 at 3:57
• It is trivial to open a door, but it's polite for me to hold it open for someone. Thanks for adding the code. – Rahul Sep 12 '14 at 4:02
• Can you show some code that backs your claim up? – user21 Sep 12 '14 at 8:22
• @user21 Yes, watch this space... – Igor Rivin Sep 12 '14 at 11:46

It appears to be just as it should be here (v10.0.0 under Windows):

randpt[n_] := With[{prept = RandomVariate[NormalDistribution[], 3]}, prept/Norm[prept]]

Needs["GeneralUtilities"]

BenchmarkPlot[{ConvexHullMesh}, Table[randpt, {#}] &,
"IncludeFits" -> True,
TimeConstraint -> 30
] Following your update I tried a plot with larger n values as indicated. The result: Indeed it seems there is an issue with scale. I am not familiar with the algorithm so I don't know if this should be considered a designed trade-off or a bug.

• BenchmarkPlot, who knew! – Igor Rivin Sep 12 '14 at 11:46
• @Igor I first learned about it here: (54865) – Mr.Wizard Sep 12 '14 at 11:48
• @Igor By the way I used With here because for the time being Module has problems: (58375). I don't know if it would have affected the results but I thought it best to eliminate that factor. – Mr.Wizard Sep 12 '14 at 11:49
• Thanks! Is the semantics of With exactly the same as Module (if yes, it is weird that it does not inherit Module's bugs...)? – Igor Rivin Sep 12 '14 at 11:53
• @Igor Please see: What are the use cases for different scoping constructs? -- With, Block, and Module all have the same syntax but work in very different ways. In this case the value of prept does not need to be changed therefore With is a better choice. – Mr.Wizard Sep 12 '14 at 11:55

The best you may hope for is the speed of TetGen as that is the underlying engine of the convex hull computation. The increase of timing you see comes from the fact that a MeshRegion needs to check the consistency. You can switch that consistency check of for ToBoundaryMesh (and there are a few other checks you can switch off, please see the documentation for ElementMesh here)

Needs["NDSolveFEM"]
randpt[n_] :=
With[{prept = RandomVariate[NormalDistribution[], 3]},
prept/Norm[prept]]

Needs["GeneralUtilities"]

BenchmarkPlot[{ToBoundaryMesh[#, "CheckIntersections" -> False] &,
TetGenConvexHull}, Table[randpt, {#}] &, 100,
"IncludeFits" -> True, TimeConstraint -> 60] TetGen is used for the convex hull computation. Unfortunately, TetGen will compute a Delaunay tetrahedralization and the hull will be extracted from that. That is why you see a speed discrepancy compared with 'qconvex'. This may change in a future version.

randpt[n_] := Module[{prept = RandomVariate[NormalDistribution[], 3]},
` 