Introduction

Given nn points in the plane, there are (n2)=O(n2)\binom{n}{2} = O(n^2) pairs of points — and naively computing or storing all their distances costs quadratic time and space. For large nn that is simply too expensive.

Well-Separated Pair Decomposition (WSPD), introduced by Paul Callahan and S. Rao Kosaraju in 1995, offers a remarkable shortcut. It finds a collection of O(n)O(n) pairs of point sets {Ai,Bi}\{A_i, B_i\} such that:

  1. Every pair of individual points (p,q)(p, q) is represented by at least one pair (Ai,Bi)(A_i, B_i) with p∈Aip \in A_i and q∈Biq \in B_i.
  2. Each pair is well-separated: the distance between the two clusters is at least ss times the diameter of either cluster, for a user-chosen separation factor s>0s > 0.

Because the clusters in a well-separated pair are far apart relative to their own sizes, every point in AiA_i is roughly the same distance from every point in BiB_i — so one representative distance stands in for all ∣Ai∣⋅∣Bi∣|A_i| \cdot |B_i| individual distances. This is the core geometric insight: spatial proximity structure lets you collapse a quadratic table into a linear one.

The construction runs in O(nlog⁡n)O(n \log n) time using a split tree (a compressed quadtree), and the number of pairs is O(sd⋅n)O(s^d \cdot n) in dd dimensions — linear for any fixed ss and dd.

Build the WSPD

Click inside the canvas to add points (up to 20), then press Build WSPD to decompose them. Each colored arc connects the centers of a well-separated pair; the shaded halos show the cluster diameters. Drag the separation factor ss slider and rebuild — a larger ss demands tighter separation, which forces more pairs.

<!-- Well-Separated Pair Decomposition (WSPD) interactive demo -->
<div class="controls">
  <label class="lbl">Separation factor s = <strong id="sval">2</strong>
    <input type="range" id="slider" min="1" max="8" step="0.5" value="2">
  </label>
  <button id="buildBtn" type="button">Build WSPD</button>
  <button id="clearBtn" type="button" class="ghost">Clear</button>
</div>
<canvas id="canvas" width="480" height="300" title="Click to add points"></canvas>
<div id="status" class="status">Click on the canvas to add points, then press Build WSPD.</div>

/* Layout and visual style for the WSPD demo */
* { box-sizing: border-box; margin: 0; padding: 0; }
body { font-family: system-ui, sans-serif; color: #222; background: #fff; }
.controls { display: flex; align-items: center; gap: .6rem; flex-wrap: wrap; margin-bottom: .5rem; }
.lbl { font-size: .88rem; display: flex; align-items: center; gap: .35rem; }
#slider { width: 110px; cursor: pointer; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 8px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
#canvas { border: 1px solid #cdd9e3; border-radius: 10px; cursor: crosshair;
          display: block; max-width: 100%; }
.status { font-size: .9rem; margin-top: .45rem; min-height: 1.3em; color: #444; font-weight: 500; }
.status.ok { color: #0a7d33; }
.status.info { color: #1d3557; }

const fmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
// WSPD: decompose n points into O(n) well-separated pairs {A,B} such that
// every point-pair (p,q) is covered by some pair (A,B) with p in A, q in B.

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const sliderEl = document.getElementById('slider');
const svalEl = document.getElementById('sval');

// Point store — user adds up to MAX_PTS points by clicking
let pts = [];
const MAX_PTS = 20;

// Color palette: one color per well-separated pair
const COLORS = [
  '#e63946','#457b9d','#2a9d8f','#e9c46a','#f4a261',
  '#6a4c93','#1982c4','#8ac926','#ff595e','#6a994e',
  '#4361ee','#f77f00','#06d6a0','#ef476f','#118ab2',
  '#ffb703','#3a86ff','#fb5607','#8338ec','#219ebc'
];

// draw(): render background, cluster halos, arcs, and point dots
function draw(pairs) {
  const W = canvas.width, H = canvas.height;
  ctx.clearRect(0, 0, W, H);

  // Draw shaded cluster halos for each well-separated pair
  if (pairs) {
    pairs.forEach((pair, i) => {
      const col = COLORS[i % COLORS.length];
      const alpha = '44';
      [pair[0], pair[1]].forEach(cluster => {
        if (cluster.length === 0) return;
        const cx = cluster.reduce((s, p) => s + p.x, 0) / cluster.length;
        const cy = cluster.reduce((s, p) => s + p.y, 0) / cluster.length;
        let rad = 0;
        cluster.forEach(p => { const d = Math.hypot(p.x - cx, p.y - cy); if (d > rad) rad = d; });
        rad = Math.max(rad + 10, 14);
        ctx.beginPath();
        ctx.arc(cx, cy, rad, 0, Math.PI * 2);
        ctx.fillStyle = col + alpha;
        ctx.fill();
        ctx.strokeStyle = col;
        ctx.lineWidth = 1.5;
        ctx.stroke();
      });
      // Connect centroids with a dashed arc to show the pair relationship
      const c0 = centroid(pair[0]), c1 = centroid(pair[1]);
      ctx.beginPath();
      ctx.moveTo(c0.x, c0.y);
      ctx.lineTo(c1.x, c1.y);
      ctx.strokeStyle = col;
      ctx.lineWidth = 1.8;
      ctx.setLineDash([5, 4]);
      ctx.stroke();
      ctx.setLineDash([]);
    });
  }

  // Draw each point as a solid circle on top
  pts.forEach(p => {
    ctx.beginPath();
    ctx.arc(p.x, p.y, 5, 0, Math.PI * 2);
    ctx.fillStyle = '#1d3557';
    ctx.fill();
  });
}

function centroid(cluster) {
  return {
    x: cluster.reduce((s, p) => s + p.x, 0) / cluster.length,
    y: cluster.reduce((s, p) => s + p.y, 0) / cluster.length
  };
}

// Compute axis-aligned bounding box and diameter for a cluster
function bbox(cluster) {
  let minX = Infinity, maxX = -Infinity, minY = Infinity, maxY = -Infinity;
  cluster.forEach(p => {
    if (p.x < minX) minX = p.x; if (p.x > maxX) maxX = p.x;
    if (p.y < minY) minY = p.y; if (p.y > maxY) maxY = p.y;
  });
  return { minX, maxX, minY, maxY,
           diam: Math.hypot(maxX - minX, maxY - minY) };
}

// Test whether clusters A and B are s-well-separated
function areSeparated(A, B, s) {
  const bA = bbox(A), bB = bbox(B);
  const diam = Math.max(bA.diam, bB.diam);
  // Minimum inter-point distance between the two clusters
  let minDist = Infinity;
  A.forEach(a => B.forEach(b => {
    const d = Math.hypot(a.x - b.x, a.y - b.y);
    if (d < minDist) minDist = d;
  }));
  return minDist >= s * diam;
}

// Core WSPD recursion: pair A and B if well-separated,
// otherwise split the larger cluster and recurse.
function wspd(A, B, s, pairs) {
  if (A.length === 0 || B.length === 0) return;
  if (areSeparated(A, B, s)) {
    pairs.push([A, B]);
    return;
  }
  // Split the cluster with the larger diameter
  const bA = bbox(A), bB = bbox(B);
  if (bA.diam >= bB.diam) {
    // Split cluster A along the x-axis midpoint (fall back to y if degenerate)
    const mid = (bA.minX + bA.maxX) / 2;
    const left = A.filter(p => p.x <= mid);
    const right = A.filter(p => p.x > mid);
    if (left.length && right.length) {
      wspd(left, B, s, pairs);
      wspd(right, B, s, pairs);
    } else {
      const midY = (bA.minY + bA.maxY) / 2;
      wspd(A.filter(p => p.y <= midY), B, s, pairs);
      wspd(A.filter(p => p.y > midY), B, s, pairs);
    }
  } else {
    // Split cluster B along the x-axis midpoint (fall back to y if degenerate)
    const mid = (bB.minX + bB.maxX) / 2;
    const left = B.filter(p => p.x <= mid);
    const right = B.filter(p => p.x > mid);
    if (left.length && right.length) {
      wspd(A, left, s, pairs);
      wspd(A, right, s, pairs);
    } else {
      const midY = (bB.minY + bB.maxY) / 2;
      wspd(A, B.filter(p => p.y <= midY), s, pairs);
      wspd(A, B.filter(p => p.y > midY), s, pairs);
    }
  }
}

// Build WSPD for all O(n^2) point pairs by starting each recursion as singletons
function buildWSPD(points, s) {
  const pairs = [];
  for (let i = 0; i < points.length; i++) {
    for (let j = i + 1; j < points.length; j++) {
      wspd([points[i]], [points[j]], s, pairs);
    }
  }
  return pairs;
}

function setStatus(msg, cls) {
  const el = document.getElementById('status');
  el.textContent = msg;
  el.className = 'status ' + (cls || '');
}

// Button handler: compute and render the WSPD
document.getElementById('buildBtn').addEventListener('click', () => {
  if (pts.length < 2) { setStatus('Add at least 2 points first.', ''); return; }
  const s = parseFloat(sliderEl.value);
  const pairs = buildWSPD(pts, s);
  draw(pairs);
  setStatus(fmt('Found %length% well-separated pairs for %length2% points.', { length: pairs.length, length2: pts.length }), 'ok');
});

document.getElementById('clearBtn').addEventListener('click', () => {
  pts = [];
  draw(null);
  setStatus('Click on the canvas to add points, then press Build WSPD.', '');
});

sliderEl.addEventListener('input', () => { svalEl.textContent = sliderEl.value; });

// Canvas click handler: add a new point
canvas.addEventListener('click', (e) => {
  if (pts.length >= MAX_PTS) { setStatus('Maximum 20 points reached.', 'info'); return; }
  const rect = canvas.getBoundingClientRect();
  const scaleX = canvas.width / rect.width;
  const scaleY = canvas.height / rect.height;
  pts.push({ x: (e.clientX - rect.left) * scaleX, y: (e.clientY - rect.top) * scaleY });
  draw(null);
  setStatus(fmt('Added point %length% — press Build WSPD.', { length: pts.length }), 'info');
});

draw(null);
setStatus('Click on the canvas to add points, then press Build WSPD.', '');
Result
Console

Watch how two clusters that are far apart collapse into a single pair, while nearby clusters need several pairs to cover all combinations. The pair count shown is always O(n)O(n) — linear in the number of points regardless of their arrangement.

The Real Complexity

Building the WSPD takes O(nlog⁡n)O(n \log n) time using a fair split tree (also called a compressed quadtree). The algorithm:

  1. Build the split tree in O(nlog⁡n)O(n \log n).
  2. For every internal node, recursively pair its two sub-trees when they are well-separated; otherwise recurse on the larger one. This produces O(sd⋅n)O(s^d \cdot n) pairs in dd dimensions.

Why O(n)O(n) pairs? In two dimensions with fixed ss, the number of pairs is at most O(s2⋅n)O(s^2 \cdot n). For typical values like s=2s = 2 that is still O(n)O(n) — a constant times the input size, not quadratic.

What you get for free — once you have the WSPD with separation factor ss:

  • Approximate nearest neighbor: for each point pp, its nearest neighbor in the entire set is represented by a well-separated pair containing pp. Scan those pairs: O(nlog⁡n)O(n \log n) total.
  • tt-spanners: connect one representative edge per pair to get a graph on nn points with O(n)O(n) edges whose shortest-path distances approximate Euclidean distances within factor t=(s+2)/(s−2)t = (s+2)/(s-2). See the closest pair article for a related geometric divide-and-conquer.
  • kk-nearest neighbors and n-body potential computations: both reduce to scanning WSPD pairs.

The key insight is that well-separation is a proxy for distance equivalence: all (p,q)(p,q) with p∈Ai,q∈Bip \in A_i, q \in B_i have distances in the range [di,di⋅(1+2/s)][d_i, d_i \cdot (1 + 2/s)], where did_i is the inter-cluster distance. Larger ss tightens this ratio.

Where It Matters

Any algorithm that touches all O(n2)O(n^2) pairs of points is a candidate for WSPD acceleration:

  • Nearest-neighbor search: GPS devices, recommender systems, and image retrieval all need to find the closest match in a large point set. WSPD reduces this to O(nlog⁡n)O(n \log n).
  • Geometric spanners: network design wants a sparse graph that still approximates all distances. A WSPD spanner has O(n)O(n) edges and a user-controlled stretch factor — invaluable for road-network approximation and sensor placement.
  • Collision detection: in game engines and robotics, objects are replaced by bounding spheres; WSPD clusters distant groups so only nearby pairs are checked each frame.
  • nn-body simulations: gravitational or electrostatic forces between nn particles naively cost O(n2)O(n^2); WSPD (and its cousin the Barnes–Hut tree) group distant particles and approximate their combined influence in O(nlog⁡n)O(n \log n).
  • Clustering and dimensionality reduction: minimum spanning trees of point sets (used in single-linkage clustering) can be built from a WSPD in O(nlog⁡n)O(n \log n). Compare with closest pair algorithms that share the same divide-and-conquer spirit.

In short, WSPD is the geometric analogue of fast multipole methods: it replaces a quadratic interaction table with a linear summary, trading exactness for a controllable approximation ratio.

Conclusion

The Well-Separated Pair Decomposition is a beautiful example of how geometric structure pays off algorithmically. Points in space are not arbitrary: nearby points form clusters, and clusters that are far apart behave almost identically with respect to inter-cluster distances. WSPD formalizes this intuition and turns it into a O(nlog⁡n)O(n \log n) construction that compresses a quadratic table of distances into a linear one.

The separation factor ss is your lever: crank it up for tighter approximations at the cost of more pairs; lower it for fewer pairs with a looser ratio. Either way the pair count stays O(n)O(n), and downstream algorithms — nearest neighbor, spanners, nn-body forces — inherit that linear efficiency.

Next time you wonder why your GPS can find a route across a continent in milliseconds, or why a physics engine can simulate thousands of particles in real time, remember: somewhere under the hood, a decomposition very much like WSPD is quietly keeping the pair count linear.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/well-separated-pair-decomposition/Content licensed under CC BY-NC 4.0.