Smallest Enclosing Circle, in depth: Welzl's randomized algorithm, why it runs in expected linear time, and a robust implementation

The two shapes of the answer

Three facts that make it solvable

Three facts make the problem tractable. First, the smallest enclosing circle is unique. If two different circles of the same minimum radius r both contained every point, the points would lie in their intersection, a lens that fits inside a circle centred midway between the two centres with radius strictly less than r, contradicting minimality.

Second, the circle is determined by at most three points on its boundary, called the support set. Either two points lie at opposite ends of a diameter, or three points lie on the circle and the triangle they form is not obtuse, so it contains the centre. If every boundary point lay within some open half of the circle, you could shift the centre towards them and shrink the radius, so the boundary points must surround the centre. Every other point is strictly inside or on the boundary and does not matter.

Third, and this is what the algorithm exploits: if a point lies outside the smallest circle of some subset, it must lie on the boundary of the smallest circle of that subset plus the point. Adding one point either leaves the circle unchanged, when the point is already inside, or forces the new circle to pass through the new point. That turns one search into a smaller search with a known boundary point, and the recursion bottoms out once three boundary points are known, because three points fix a circle.

These facts give a slow algorithm: try every pair and triple as a support set and keep the smallest candidate that contains all points, O(n⁴) time. It is a useful test oracle and nothing more.

Welzl's algorithm

Welzl's algorithm applies the third fact recursively. It processes points in random order and carries a set R of points already known to lie on the boundary. The recursive version is elegant, but its recursion depth is linear in n, which overflows the call stack in most languages for large inputs. The iterative form, which follows the same logic after a single shuffle, is what production code uses: the outer loop builds the circle of the first i points, and whenever a point falls outside, the inner loops rebuild the circle of the earlier points with that point pinned to the boundary.

# Recursive form (Welzl 1991). P: points not yet processed, R: points known
# to lie on the boundary, |R| <= 3.
welzl(P, R):
    if P is empty or |R| == 3:
        return trivial_circle(R)       # 0, 1, 2 or 3 boundary points
    p = remove a random point from P
    D = welzl(P, R)
    if p is inside D:
        return D
    return welzl(P, R + {p})           # p must be on the boundary

# Iterative form after one random shuffle: three nested loops.
shuffle(points)
C = circle of radius 0 at points[0]
for i in 1..n-1:
    if points[i] outside C:
        C = circle of radius 0 at points[i]          # points[i] on boundary
        for j in 0..i-1:
            if points[j] outside C:
                C = circle with diameter points[i], points[j]
                for k in 0..j-1:
                    if points[k] outside C:
                        C = circumcircle(points[i], points[j], points[k])
return C

Read it as three nested instances of one idea: the outer loop pins nothing, the middle loop pins i, the inner loop pins i and j, and a third outside point pins all three, fixing the circumcircle for that prefix.

Why it runs in expected linear time

The inner loops look quadratic or cubic, yet the expected running time is linear. The proof is backwards analysis: instead of asking how likely the next point is to change the circle, fix the set of the first i points and ask how likely it is that the last of them, in random order, was a support point. The circle of the first i points is determined by at most three of them. Since the order of those i points is uniformly random, the probability that the i-th one is among the at most three support points is at most 3/i.

The outer loop does O(1) work per point, plus O(i) work for the middle loop with probability at most 3/i. The expected cost is the sum of 3/i times O(i), which is O(1) per point and O(n) overall. The same argument applies one level down: inside the middle loop with point i pinned, the circle of the first j points plus i has at most two free support points, so point j triggers the inner loop with probability at most 2/j, and the inner loop's O(j) cost again averages to O(1) per step. Linearity of expectation adds the levels, giving expected O(n) total. The constant is small: in practice the algorithm touches each point a handful of times.

The bound is about the random order, not the input. No input is bad, only an unlucky shuffle, and the probability of running far above the expectation falls off quickly. If an adversary can choose the order, for example by feeding points already sorted along a spiral, skipping the shuffle can make the loops quadratic. Always shuffle, and for reproducible results seed the generator. The same randomized incremental pattern, analysed backwards, appears in treaps and in randomized convex hull and triangulation algorithms.

A robust implementation

The implementation below is the iterative form with three robustness measures that the textbook pseudocode leaves out: duplicates are removed, the circumcircle is computed relative to one of its points, and the containment test has a tolerance relative to the radius.

import math
import random

def _two(a, b):
    cx, cy = (a[0] + b[0]) / 2, (a[1] + b[1]) / 2
    return (cx, cy, math.hypot(a[0] - cx, a[1] - cy))

def _three(a, b, c):
    # Work relative to a to reduce cancellation for large coordinates.
    bx, by = b[0] - a[0], b[1] - a[1]
    cx, cy = c[0] - a[0], c[1] - a[1]
    d = 2 * (bx * cy - by * cx)
    scale = bx * bx + by * by + cx * cx + cy * cy
    if abs(d) <= 1e-12 * scale:
        # Nearly collinear: no stable circumcircle, use the widest pair.
        return max(_two(a, b), _two(a, c), _two(b, c), key=lambda t: t[2])
    b2, c2 = bx * bx + by * by, cx * cx + cy * cy
    ux = (cy * b2 - by * c2) / d
    uy = (bx * c2 - cx * b2) / d
    return (a[0] + ux, a[1] + uy, math.hypot(ux, uy))

def _inside(circle, p, rel=1e-9):
    cx, cy, r = circle
    return math.hypot(p[0] - cx, p[1] - cy) <= r * (1 + rel) + 1e-12

def smallest_enclosing_circle(points, seed=None):
    pts = list({(float(x), float(y)) for x, y in points})  # drop duplicates
    if not pts:
        raise ValueError("no points")
    random.Random(seed).shuffle(pts)
    c = (pts[0][0], pts[0][1], 0.0)
    for i in range(1, len(pts)):
        if _inside(c, pts[i]):
            continue
        c = (pts[i][0], pts[i][1], 0.0)
        for j in range(i):
            if _inside(c, pts[j]):
                continue
            c = _two(pts[i], pts[j])
            for k in range(j):
                if not _inside(c, pts[k]):
                    c = _three(pts[i], pts[j], pts[k])
    return c   # (centre_x, centre_y, radius)

The containment tolerance matters more than it looks. Points that lie exactly on the circle are common, because every support point does, and the distance to them is computed with rounding error. Without a tolerance a support point can test as outside its own circle and trigger needless rebuilds. A relative tolerance of about 10^-9 is far above double-precision rounding for well-scaled inputs and far below any meaningful geometric difference. For a guaranteed enclosure, inflate the final radius by the same tolerance.

Worked example: five points by hand

Take five points: (0,0), (4,0), (0,3), (1,1) and (2,2). The first three form a right triangle, so the answer should be the circle on the hypotenuse from (4,0) to (0,3): centre (2, 1.5), radius 2.5. Suppose the shuffle produced the order (1,1), (2,2), (4,0), (0,3), (0,0). The trace shows each time a point falls outside and what replaces the circle.

Points landed exactly on the boundary twice, (2,2) at step j=1 and (0,0) at the end, which is the case the tolerance exists for. And the circumcircle through (0,3), (2,2) and (1,1) was correct for its prefix but discarded one step later, which is normal: each loop level answers only for the points it has seen.

Numerical robustness

Most failures in practice are numerical. Large coordinates, such as projected map coordinates in millions of metres, lose precision when squared; translate the points so the bounding-box centre is at the origin, and translate the answer back. Nearly collinear triples give a huge radius from a tiny determinant; exact arithmetic never requests that circumcircle, but rounding can, hence the widest-pair fallback. Duplicates produce a zero determinant.

When the answer must be exact, use integer or rational coordinates with a division-free in-circle test, or adaptive-precision predicates like those used for robust convex hulls; the orientation and in-circle tests fail the same ways under floating point.

On latitude and longitude, a circle in degrees is an ellipse on the ground; project to a local metric system first.

Alternatives and generalisations

The problem is an instance of what Sharir and Welzl called LP-type problems: like a low-dimensional linear program, its solution is fixed by a small basis, here at most three points, and the same randomized algorithm solves both. In d dimensions the minimum enclosing ball has a basis of at most d+1 points and Welzl's method still runs in expected linear time in n, but the constant grows very fast with d, so in high dimensions approximate methods based on core sets are used instead. For very large inputs, only hull points can be support points, so computing the hull first, for example with a Graham scan, often shrinks the input sharply.

Using it in practice

Use the circle where its definition matches the question. For facility placement it is the 1-centre problem: it minimises the worst distance, so one far outlier moves the answer; decide whether outliers are errors or customers. For collision detection it is a broad-phase filter: pairs whose circles do not overlap cannot collide, and for many objects a spatial index such as a k-d tree finds the candidate pairs. For clustering, report the radius beside a percentile-based one, because it reacts to single points.

Failure modes

• Using the centroid or bounding box. Both give a valid enclosing circle that is not minimal; measure the difference on your data before accepting it.
• No shuffle. Pre-sorted input, such as points ordered by time along a track, can make the loops quadratic.
• Exact comparisons. Support points test as outside their own circle and the loops thrash.
• Recursive version on large input. Recursion depth n overflows the stack.

What to do next

1. Implement the iterative algorithm with a shuffle, duplicate removal and a relative containment tolerance.
2. Write a brute-force oracle over pairs and triples and compare the two on thousands of random small inputs, including collinear and duplicate points.
3. Add tests for the hand-traced example above and for points on a circle, on a line and at one location.
4. Centre large coordinates near the origin before computing, and project geographic data to a metric system.
5. Decide how outliers should affect the result, and remove or keep them deliberately before computing.
6. If you need a deterministic or exact answer, move to exact predicates or a geometry library rather than tightening the tolerance.