Is the enclosing circle algorithm correct?

[robinhouston]

I have some doubts about the correctness of the enclosing circle algorithm.

You’re using Welzl’s algorithm adapted for enclosing circles rather than points. It is known that this doesn’t work in general, as detailed in Chapter 5 of Kaspar Fischer’s thesis (2005). The fundamental problem is that the correctness proof (which works for an enclosing circle of a set of points) doesn’t carry over to enclosing differently-sized circles.

But you’re also using the move-to-front heuristic, and this makes matters a bit less certain. Fischer’s counterexample doesn’t appear to bite when move-to-front is in effect. (I have tried all permutations of his input circles, for example.)

So it’s possible that the algorithm is in fact correct, and that the move-to-front rule somehow prevents bad cases from happening. I don’t know. I do think the correctness or otherwise of this algorithm is an open question. (If there’s a proof of correctness out there I’d love to hear about it.)

The algorithm implemented in CGAL is provably correct, I believe.

[mbostock]

Initial condition from Fischer’s thesis:

Output of implementation:

The algorithm verifies that every circle (in L) is contained inside the enclosing circle, so it should be correct in the sense that it returns an enclosing circle. If you can demonstrate that’s not the case with a specific input, please let me know! I would be very surprised.

Whether or not that enclosing circle is the minimum enclosing circle, that is harder to prove and beyond my capacity.

I’m not really sure what else to do with this issue so I’m closing it, but please contribute if you can make it better. Thanks, Robin!

[robinhouston]

Sure, that’s probably best. I suspect it’s fine, actually – I just don’t completely understand why. I’ll let you know if I find a concrete problem!

[mbostock]

If you learn and can teach me that would be most excellent. Cheers!

[robinhouston]

I now pretty much believe that this algorithm is correct – which is interesting, since that is counter to the prevailing wisdom in the academy.

If I’m right, correctness doesn’t depend on the order elements are chosen in. It’s possible for a recursive call to encloseN to return undefined in examples like Fischer’s, but as long as we treat the undefined circle as not enclosing anything (which your implementation does) then it comes out right in the end.

If I can pin down the remaining details of the proof, would you be interested in co-authoring a short paper? I think this would be of some interest to the computational geometry people, and it’s a joint effort in that you developed the algorithm and implicitly conjectured its correctness.

[mbostock]

That would be amazing! I would love to.

[mbostock]

It’s possible for a recursive call to encloseN to return undefined in examples like Fischer’s, but as long as we treat the undefined circle as not enclosing anything (which your implementation does) then it comes out right in the end.

If I understand this correctly, the following change would cause the algorithm to throw an error (but it doesn’t, or at least I can’t get it to by throwing random inputs at it… although those inputs were non-overlapping circles, which might be a factor):

switch (B.length) {
  case 0: break;
  case 1: circle = enclose1(B[0]); break;
  case 2: circle = enclose2(B[0], B[1]); break;
  case 3: circle = enclose3(B[0], B[1], B[2]); break;
  default: throw new Error;
}

If encloseN were called with B.length > 3, then it would start behaving degenerately, recursing for each circle in L but in all cases returning undefined. So it seems like you’d want to short-circuit that and change the code to:

switch (B.length) {
  case 0: break;
  case 1: circle = enclose1(B[0]); break;
  case 2: circle = enclose2(B[0], B[1]); break;
  case 3: circle = enclose3(B[0], B[1], B[2]); break;
  default: return;
}

But if your theory is correct we should definitely look for an input that can reproduce this case so that we can understand what’s happening. I added a test script in 98f0f6a2353a7b84afcba0545cc8d83ca866cdd8, but a few thoughts:

• The test script runs d3.packSiblings before d3.packEnclose, but we might have better luck reproducing edge cases if we allow the input circles to be overlapping.

• d3.packEnclose is supposed to shuffle the circles before computing the enclosing circle, but it currently doesn’t (#83). This is convenient for reproducible test cases but we should keep this in mind for the future…

[robinhouston]

Ah yes. This is a really interesting question, and I have no idea what the answer is. I’ve mainly been analysing the algorithm without the move-to-front optimisation, because it’s easier to keep track of what’s where. (I think the algorithm gives the right answer whatever order you use, even if it’s like Welzl’s unoptimised version where you pick the circle randomly at each step.)

Like you I’ve noticed that these cases are hard to reproduce when move-to-front is in effect. So there’s another question here: does move-to-front actually prevent these cases from occurring at all?

Without move-to-front it isn’t hard to find a case where a recursive call returns undefined. Here’s the implementation I’ve been experimenting with. It assumes the circles have a name property for the trace output:

// Pretty-printers
function pp_circle(a) {
  return a.name;
  // return "(" + a.x + ", " + a.y + ", " + a.r + ")";
}

function pp_array(B) {
  return B.map(pp_circle).join(",");
}

function pp_list(L) {
  var s = [];
  for (var node = L.head; node; node = node.next) {
    s.push(node._);
  }
  return pp_array(s);
}

// Returns the smallest circle that contains circles L and intersects circles B.
function encloseN(L, B, level=0) {
  console.log("  ".repeat(level) + "encloseN", pp_list(L), "[" + pp_array(B) + "]");
  var circle,
      l0 = null,
      l1 = L.head,
      l2,
      p1;

  switch (B.length) {
    case 1: circle = enclose1(B[0]); break;
    case 2: circle = enclose2(B[0], B[1]); break;
    case 3: circle = enclose3(B[0], B[1], B[2]); break;
  }

  while (l1) {
    p1 = l1._, l2 = l1.next;
    if (!circle || !encloses(circle, p1)) {

      // Temporarily truncate L before l1.
      if (l0) l0.next = null;
      else L.head = null;

      B.push(p1);
      circle = encloseN(L, B, level + 1);
      B.pop();

      // Reconnect the truncated list L.
      if (l0) l0.next = l1;
      else L.head = l1;
    }

    l0 = l1;
    l1 = l2;
  }

  console.log("  ".repeat(level) + "returning", JSON.stringify(circle) || "undefined");
  return circle;
}

With this, if you define

const a = {name: "a", x: 365, y: 90, r: 5},
      b = {name: "b", x: 430, y: 90, r: 5},
      c = {name: "c", x: 400, y: 150, r: 10},
      d = {name: "d", x: 400, y: 100, r: 30},
      D = {name: "D", x: 400, y: 500, r: 10};

and run enclose([ d, D, c, b, a ]) the output is:

encloseN a,b,c,D,d []
  encloseN  [a]
  returning {"x":365,"y":90,"r":5}
  encloseN a [b]
    encloseN  [b,a]
    returning {"x":397.5,"y":90,"r":37.5}
  returning {"x":397.5,"y":90,"r":37.5}
  encloseN a,b [c]
    encloseN  [c,a]
    returning {"x":383.7596775638102,"y":122.15944725224608,"r":42.23110997362451}
    encloseN a [c,b]
      encloseN  [c,b,a]
      returning {"x":397.5,"y":114.42982663671509,"r":45.65791964058115}
    returning {"x":397.5,"y":114.42982663671509,"r":45.65791964058115}
  returning {"x":397.5,"y":114.42982663671509,"r":45.65791964058115}
  encloseN a,b,c [D]
    encloseN  [D,a]
    returning {"x":382.7126412472094,"y":297.4909403244528,"r":213.24559533559886}
    encloseN a [D,b]
      encloseN  [D,b,a]
      returning {"x":397.5,"y":296.23512452682576,"r":213.78021119970956}
    returning {"x":397.5,"y":296.23512452682576,"r":213.78021119970956}
  returning {"x":397.5,"y":296.23512452682576,"r":213.78021119970956}
  encloseN a,b,c,D [d]
    encloseN  [d,a]
    returning {"x":394.51904934551027,"y":98.43401409871723,"r":35.7002747232013}
    encloseN a [d,b]
      encloseN  [d,b,a]
      returning {"x":397.5,"y":92.80124908814332,"r":37.62049963525733}
    returning {"x":397.5,"y":92.80124908814332,"r":37.62049963525733}
    encloseN a,b [d,c]
      encloseN  [d,c,a]
      returning {"x":396.6489361702128,"y":114.8936170212766,"r":45.26595744680849}
      encloseN a [d,c,b]
        encloseN  [d,c,b,a]
        returning undefined
      returning undefined
    returning undefined
    encloseN a,b,c [d,D]
    returning {"x":400,"y":290,"r":220}
  returning {"x":400,"y":290,"r":220}
returning {"x":400,"y":290,"r":220}

[mbostock]

It returns the right answer even without the move-to-front heuristic? That is interesting. I was going to suggest that maybe the move-to-front heuristic does ensure correctness somehow.

[robinhouston]

Yes. I suspected at first that move-to-front was making it work – as I said earlier in this thread – but now I think it works for any ordering.

But now there is a separate question whether move-to-front has the further effect that no recursive call will ever return undefined. That seems harder to analyse.

[mbostock]

Yes, sorry, I forgot you mentioned that already!

[robinhouston]

Update: I’m actively working on this in my sadly limited free time, and I’m hoping to have a (very) rough draft ready in the next week.

[mbostock]

Great! Whenever you’re ready and want to create a private repo and add me as a collaborator I’d be happy to take a look.

[robinhouston]

Perfect! Yeah, that’s exactly what I am planning to do.

[VladanMajerech]

The first mentioned algorithm is incorrect, the proposed implementation with array of randomly permuted points repermuted whenever point is fixed is incorrect as well (it's propper implementation of the first one). Second proposal just with move front and no repermutation works well, but I am not sure it is proven to have the specified randomized complexity.

The algorithm works well with patch "the 4 last points are not selected unless they are the only remaining". I am going to write an article about the patch.

[robinhouston]

@VladanMajerech Yes. This algorithm has been entirely replaced, because we realised it was wrong. We’re now using the MSW algorithm: see @mbostock’s write-up at https://beta.observablehq.com/@mbostock/miniball

[VladanMajerech]

Actually in "A Subexponential Bound for Linear Programming" by Matousek, Sharir, Welzl. The need to fix d+2 points (basis(B,h) call to be sure the diameter increases) are mentioned. So I probably resign to writing a paper about it unless wikipedia forces me to have argument based on the published paper.

[robinhouston]

Reading back through this thread, it’s very misleading. Sorry! The discussion moved elsewhere (maybe to email?) after this thread was written, so reading this you’re only seeing the first part of the story.

I no longer think it’s plausible that Welzl’s algorithm works in general to find the minimum enclosing circle; I seem to remember @mbostock even found an explicit counterexample.

[robinhouston]

@VladanMajerech If you’re saying you can prove it’s correct with move-to-front and no shuffling, I would definitely be interested in reading that! I tried to prove that at one point (as mentioned in this thread), but I did not succeed.

[VladanMajerech]

@robinhouston OK, You have encouraged me to continue in writing the article. I have not finished writing the proof, but I am almost sure I have it. ... sketch... Last 3 points define the circle (of radius $r$). Enclosing circle radius for last 3 points could decrease, but not for last 4 ... as it should enclose defining points of previous circle. The enclosing circle radius for last 4 points only grows. When a point $p$ (when $k$ points remain) is going to be fixed, the circle defined by last 3 encloses the $k-1$ points. When $p$ is not among a set of the $k-1$ points, the enclosing circle's radius is at most $r$. This is why $p$ is defining for all enclosing circles for subsets of the "current" last $k$ points with radius bigger than $r$. Yes, move to front without reshufling is correct.

My implementation differs, I use array and swap selected index with $n$-th index when called for $n$ vertices. The change is the selection is not from all $n$ vertices, but from $n-4$ last vertices. When $4$ vertices remain, I select random nonfixed one. The randomized complexity upper bound have $d/(n-4)$ instead of $d/n$, but for $n>8$ it's at most twice the original bound (for $n<8$ the complexity is constant). So the randomized complexity is still $O(n)$ ... with constant at most $8$ times higher than in the original paper (when we have to fix at most $d=3$ points). (Actually I bet the probability is typically much smaller as often we have defining points even when we have not fixed them yet, only in the case we know, the radius is not sufficient, the probability chosen vertex is defining grows from $d/(n-f)$ to $d/(n-4)$)

Move to front and reshufling affecting last $n-4$ vertices would do the same job. I am not sure how the randomized complexity is affected without reshufling.

You can see in the already presented picture of counterexample of original Welzl. The enclosing circle diameter was decreasing ... what allowed the fixed point to leave the boundary of the enclosing circle.

[robinhouston]

@VladanMajerech In case it’s useful or interesting, there’s a very incomplete draft of a paper here that I started writing in 2017, and then stopped working on when it became clear that the MSW algorithm could be implemented more efficiently than Welzl’s algorithm, which meant it no longer seemed to have any practical significance.

[VladanMajerech]

BTW: MSW is equivalent to not selecting from last $4$ points. This implicitly maintains basis in last $3$ points whenever $4$ points are recomputed. The one extra point serves to recompute the basis "using Welzl's algorithm" (which is correct for 4 points).

[VladanMajerech]

I have updated the English wiki page. As my patch is equivalent to MSW, there is no reason to publish it as a paper 23 years after MSW corrected it.

[VladanMajerech]

I have run repeated experiments with 1000 and 5000 random points and move to front (of all) without reshufling wins ... it saves about 3% reccurent calls compared to reshufling of n-4 points when we are fixing point with n remaining. Interesting are statistics for the variant with x->0->1->2->3->x permutation without reshufling. The number of calls to trivial case is worse than move to front, but better to reshufling n-4 but the number of all recurent calls is worse than in both previous algorithms.

[VladanMajerech]

So finally as David Epstein rejected my counterexample, I have reread the Welzl's original paper and I must say the algorithm is correct (actually only my argument it is incorrect was wrong). There is trick that when 3 points are fixed it returns circle determined by them even when they do not enclose the whole set. When the enclosing circle's diameter was decremented, there is surely point "on the stack" before the last fix outside the circle, so there would be less fixed vertices and the vertex will be outside the subresult, so the wrong solution would be rejected, the point fixed ...

[robinhouston]

The algorithm is correct for a set of points, but not for a set of circles, I think. Which Wikipedia article is this?

[VladanMajerech]

So finally I have rerun the tests with corrected implementation and it is still failing. https://en.wikipedia.org/wiki/Smallest-circle_problem