Speed of ConvexHullMesh

Question

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[3], {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. — , Sep 12 '14 at 03: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 03: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. — , Sep 12 '14 at 04:02
@IgorRivin, you posted a list of lists with some numbers; 40000 appears twice. Why is it so hard to post the exact code you run? I am not saying what you say is not correct but you make it much harder to reproduce than it must. — user21, Sep 12 '14 at 19:46
@IgorRivin, could you run this and tell me how that performs? "Needs["TetGenLink`"] experiment[n_] := AbsoluteTiming[TetGenConvexHull[Table[randpt[3], {n}]];][[1]] Table[{i, experiment[i]}, {i, 10000*{1, 2, 3, 4, 5}}]" — user21, Sep 12 '14 at 19:53
@user21 I am not at the same computer now, but can certainly try on the macbook I am on. — Igor Rivin, Sep 12 '14 at 19:53
@user21 Also sprach der MacBook Pro: {{10000, 1.407865}, {20000, 2.679615}, {30000, 3.890393}, {40000, 4.975892}, {50000, 6.142689}} — Igor Rivin, Sep 12 '14 at 19:58
@user21 my conjecture is all of the time is spent in the communication overhead (otherwise, I am unable to explain either why TetGen is so much slower than qhull, or why the times are actually growing slightly sublinearly) — Igor Rivin, Sep 12 '14 at 19:59
@user21 Do you mean "with what arguments"? Otherwise, I called it as a standalone program, though I believe the same effect can be achieved by using MATLink (since Matlab uses qhull internally, I believe). — Igor Rivin, Sep 15 '14 at 15:16
@IgorRivin, I'd like to know how you generated the input for qhull, how you exported it, how you called it on the command line, with which opetions; everything I need to know to reproduce your claim. — user21, Sep 15 '14 at 15:59
@user21 qconvex informs me that it took it 0.262secs to compute the hull. — Igor Rivin, Sep 15 '14 at 16:10
In version 13 I see {{10000, 0.233783}, {20000, 0.489283}, {30000, 0.701157}, {40000, 0.952774}, {50000, 1.21516}} — Greg Hurst, Apr 08 '22 at 15:07

Mr.Wizard · Answer 1 · 2014-09-12T14:38:55.303

10

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[3], {#}] &,
 "IncludeFits" -> True,
 TimeConstraint -> 30
]

enter image description here

Following your update I tried a plot with larger n values as indicated. The result:

enter image description here

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.

edited Sep 12 '14 at 14:38

answered Sep 12 '14 at 11:45

Mr.Wizard

271,378
34
587
1,371

1

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
So, should the question be closed? – Mark McClure Sep 12 '14 at 11:57
@Mark I'll wait to see if Igor can demonstrate a case with different apparent complexity. p.s. I voted for your answer of course. – Mr.Wizard Sep 12 '14 at 11:58
Ah, I see! I tend to use Module as a replacement for the Scheme let, but most of the time With is, in fact, the right thing if you program in a functional style (as I tend to). – Igor Rivin Sep 12 '14 at 11:59
@MarkMcClure I don't understand your question - I realize that it was before my updates, but is your a priori assumption that I am hallucinating? Just wondering... – Igor Rivin Sep 12 '14 at 16:24
@Mr.Wizard It is clearly a bug. In fact, if you look at Szabolcs' MATlink, the primary example for why you might want to link with MATLAB is Delaunay triangulation and its various off-shoots (like convex hull). It is a shame, since Mathematica's mesh primitives are quite nice to use, but this performance is not really acceptable. – Igor Rivin Sep 12 '14 at 16:26
@Igor Oliver's answer might change your mind about that. – Mr.Wizard Sep 13 '14 at 00:24

user21 · Answer 2 · 2014-09-16T11:17:38.687

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["NDSolve`FEM`"]
Needs["TetGenLink`"]
randpt[n_] := 
 With[{prept = RandomVariate[NormalDistribution[], 3]}, 
  prept/Norm[prept]]

Needs["GeneralUtilities`"]

BenchmarkPlot[{ToBoundaryMesh[#, "CheckIntersections" -> False] &, 
  TetGenConvexHull}, Table[randpt[3], {#}] &, 100, 
 "IncludeFits" -> True, TimeConstraint -> 60]

enter image description here

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.

score 7 · Answer 3 · answered Sep 12 '14 at 09:38

7

Quadratic?

randpt[n_] := Module[{prept = RandomVariate[NormalDistribution[], 3]}, 
  prept/Norm[prept]];
experiment[n_] := First[AbsoluteTiming[ConvexHullMesh[Table[randpt[3], {n}]];]];
ListPlot[experiment /@ Range[500]]

enter image description here

answered Sep 12 '14 at 09:38

Mark McClure

32,469
3
103
161

See the various edits to the original post... – Igor Rivin Sep 12 '14 at 15:20

Speed of ConvexHullMesh

3 Answers3