The Closest Pair Problem:

An Efficient Divide and Conquer Solution

Introduction

Closest Pair Problem

• Problem: Given n points in a plane, find the distance between the pair of points that is closest together
• Each point is an (x,y) pair


• Complexity goal: $T(n) = \Theta(n\lg n)$


• Approach: Divide and conquer


• Needed Recurrence: $T(n) = 2T(n/2) + \Theta(n)$

Outline

1. Review Master Method


2. Solutions for 1 dimension
• Brute Force: Θ(n²)
• DivCon: Θ(n + n lg n) = Θ(n lg n)


3. Solutions for 2 dimensions
• Brute Force: Θ(n²)
• DivCon - Naive: Θ(n²)
• DivCon - Slab - Naive: Θ(n²)
• DivCon - Slab - Many Sorted: Θ(n lg² n)


• DivCon - Slab - Sorted Once: Θ(n lg n)

Clickable Text: Yellow and Light Blue

• Yellow text is clickable to open an answer or explanation


• Light blue answers and explanation disappear when clicked

Review: Master Method

Master Method: Book version

• Solutions for recurrence $T(n) = aT(n/b) + n^d$
• $aT(n/b)$: cost of $a$ calls, each of size $n/b$
• $n^d$ for dividing data and then recombining subproblem solutions


• Three cases:

Case		$a :: b^d$		$T(n)$
1		$a \lt b^d$		$ n^d$
2		$a = b^d$		$n^d \lg n$
3		$a \gt b^d$		$n^{\log_b a}$

Master Method: Alternate Conditions

• Solutions for recurrence $T(n) = aT(n/b) + n^d$
• Alternate Conditions:
• Same cases, easier to understand $T(n)$
• Take $\log_b$ of both sides of $a :: b^d$ (remember: $\log_b b^d = d$)


• Three cases:

Case		$\log_b a :: d$		$T(n)$
1		$\log_b a \lt d$		$ n^d$
2		$\log_b a = d$		$n^d \lg n = $
				$n^{\log_b a}\lg n $
3		$\log_b a \gt d$		$n^{\log_b a}$

1D Closest Pair

1D Closest Pair Problem - Brute Force

• Problem: Given n integers, find the distance between the pair of ints that is the closest
• What do you think the brute force complexity is?
• What do you think the best case complexity is?
• Brute Force Code:

   type IntArray is Array(1 .. n) of Integer
   A: IntArray
   d: Float
   ...

   procedure closest_brute(A: IntArray, d: out Float)
      d := Float'last;

      -- Look at each pair, but only once
      for i in A'range
         for j in i + 1 .. A'last
            dis := distance (A(i), A(j)) -- A(j) - A(i)
            d := (if dis < d then dis else d)
            

• T(n) = $\displaystyle \sum_{i=1}^n i$ = n * (n+1) / 2 = Θ(n ²)
• PDF figure

1D Closest Pair Problem - Divide and Conquer

• Implementation:

   procedure closest(A: IntArray, d: out Float)
      sort(A)
      cp_helper(A, 1, A'last, d)

   procedure cp_helper(A: Array, l, r: Positive; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         mid := (l + r) / 2
         cp_helper (A, l,       mid, dL)
         cp_helper (A, mid + 1, r,   dR)
         d := (if dL < dR then dL else dR)
    

         --  Are we done? 
      

         across_dis := distance (A(mid), A(mid+1)) --  A(mid + 1) - A(mid) 

         d := (if across_dis < d then across_dis else d)
      

• T(n) = 2T(n/2) + Θ(1) = $ \Theta(n^{\log_b a})$ = $ \Theta(n^{\log_2 2})$ = $ \Theta(n^1) = \Theta(n)$


• Is this the actual performance? No. Must include the sort: $ \Theta(n) + \Theta(n\lg n) = \Theta(n\lg n)$


• Still pretty amazing!


• Hmmm, once the list is sorted, how can we solve the problem?

Now 2D

Closest Pair Problem

• Problem: Given n points in a plane, find the distance betweeen of points that is closest together
• Each point is an (x,y) pair


• Complexity goal: $T(n) = \Theta(n\lg n)$


• Needed Recurrence: $T(n) = 2T(n/2) + \Theta(n)$

Brute Force

   type Point is record x, y: Natural; end record
   type PointArray is Array(1 .. n) of Point

   A: Array(1 .. n) of Point
   d: Float
   ...

   procedure closest_brute(A: PointArray, d: out Float)
      d = Natural'last

      -- Look at each pair, but only once
      for i in 1 .. A'last loop
         for j in i + 1 .. A'last loop
            tmp_dis = distance(A(i), A(j))
            if tmp_dis < d then
               d := tmp_dis
    

• T(n) = $\displaystyle \sum_{i=1}^n i$ = n * (n+1) / 2 = Θ(n ²)


• Compared to 1D?? What's different?

Divide and Conquer: Naive

• PDF figure

   procedure closest(A: PointArray, d: out Float)
      sort_by_x(A)
      cp_helper(A, 1, A'last, d)

   procedure cp_helper(A: Array, l, r: Positive; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         -- Divide: find closest of left and right halves
         mid := (l + r) / 2
         cp_helper (A, l,       mid, dL)
         cp_helper (A, mid + 1, r,   dR)
         d := (if dL < dR then dL else dR)
    

         --  Are we done? 
      

         -- Conquer: Compare each point on left with each on right
         for i in l .. mid loop
            for j in mid + 1 .. r loop
               tmp_dis = distance(A(i), A(j))
               if tmp_dis < d then
                  d := tmp_dis
      

• T(n) = 2T(n/2) + n/2 * n/2 = Θ(n ²)

Improvement: Use a Slab: Naive Version

• Idea: Only look at pairs close to midline
• Implementation: Fill a slab with all points that are within distance d of midline

   procedure closest(A: PointArray, d: out Float)
      sort_by_x(A)
      cp_helper(A, 1, A'last, d)

   cp_helper(A: Array, l, r: Positive; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         -- Divide: find closest of left and right halves
         mid := (l + r) / 2

         cp_helper(A, l,     mid, dL)
         cp_helper(A, mid+1, r,   dR)

         d = min(dL, dR)

         -- Conquer: Only compare points close to midline

         -- Fill slab with all points within d left or right of mid
         slab : Array (l .. r) of Point
      

                fillslab(A, l, r, d, slab, slab_size)

         fillslab(A: Array; l, r: Positive, d: Float; slab: out PointArray, num_added: out Natural)
            num_added := 0
            midx := A((l + r) / 2).x
            for p of A(l..r)
               if abs(p.x - midx) < d then
                  num_added := @ + 1
                  slab(num_added) := p
      

         --  Search the slab
         for i in l .. slab_size loop
            for j in i + 1 .. slab_size loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis
    

• Q: What is complexity of searching the slab, if it contains n/10 elements from each side? A: $n^2/100$


• Complexity: T(n) = 2T(n/2) + $\Theta(n^2)$ = $\Theta(n^2)$
• And, in the worst case, the slab could contain all of the points!


• Q: How to avoid comparing all pairs in the slab?
• A: Compare each slab point with a fixed number of other slab points above
• But how to do that? See below

Checking the Slab More Efficiently

• Assume: the slab points are sorted by Y
• Later we look at Y-sorting them


• Then our slab-check only needs to look at the next 7 points, like this:

         --  Search the Y-sorted slab
         for i in l .. slab_size - 1 loop
            for j in i + 1 .. min(i + 7, slabsize) loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis
    

• Clearly this loop executes about $7\times \text{ slab_size}$ times
• This is a big win since it's $\Theta(n)$
• But why can you limit the inner loop to 7 repetitions?


• Clarifications:
• Only consider points above: pairs with points below have already been checked
• We don't care if points on opposite side of midline are checked
• The point is that no more than 7 points ever have to be checked!

Why 7?

• Consider point A(i)
• Clearly we don't have to look at points that are more than distance d above A(i)


• When checking point A(i), consider a rectangle of dimension $d\times 2d$ above point A(i)
• In other words, the bottom of the rectangle is at A(i).y
• The rectangle consists of two $d\times d$ squares, one on each side of the midline
• The bottom of one of the squares contains point A(i)


• Q: How many points can be within this $d\times 2d$ rectangle? A: No more than 8
• To see this, divide each square into 4 smaller squares
• The sides of the smaller squares are length $d/2$
• The diagonal of the smaller square is $d/\sqrt 2 \lt d/1.4 \lt d$
• $\sqrt 2 \approx 1.414$
• Thus, a smaller square can hold no more than one point
• Thus the $d\times 2d$ rectangle can hold no more than 8 points


• Clarification: The actual maximum number of points in the $d\times 2d$ rectangle is 6 (if sides and top are excluded), but this is harder to prove

Sort the Slab - Many Times

• Sorting each slab that is created makes searching each slab $\Theta(n)$
• But what is the performance?

   procedure closest(A: PointArray, d: out Float)
      sort_by_x(A)
      cp_helper(A, 1, A'last, d)

   procedure cp_helper(A: Array, l, r: Positive; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         -- Divide: find closest of left and right halves
         mid := (l+r)/2

         cp_helper(A, l,     mid, dL)
         cp_helper(A, mid+1, r,   dR)

         d = min(dL, dR)

         -- Conquer: Only compare points close to midline

         -- Fill slab with all points within d
         --   left or right of mid
         slab : Array (l .. r) of Point
      

                fillslab(A, l, r, d, slab, slab_size)

         procedure fillslab(A: Array; l, r: Positive, d: Float; slab: out PointArray, num_added: out Natural)
            num_added := 0
            midx := A((l + r) / 2).x
            for p of A(l..r)
               if abs(p.x - midx) <= d then
                  num_added := @ + 1
                  slab(num_added) := p
      

         --  Sort the slab (later we see how to avoid this sort!)
         sort(slab(1 .. slab_size));

         --  Search the slab
         for i in l .. slab_size - 7 loop
            for j in i + 1 .. min(i + 7, slabsize) loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis
    

• Complexity: T(n) = 2T(n/2) + $\Theta(n\lg n + 7n)$ = $2T(n/2) + \Theta(n\lg n)$ = $\Theta(n\lg^2 n)$
• Q: Where did this result come from?
• A: From a generalization of the Master Method.


• This is better than $T(n)=\Theta(n^2)$!
• Q: But can we do better?
• A: Yes! Sort by Y once, and distribute and pass the sorted array!

How to Avoid Sorting the Slab Each Time?

• The trick has two parts:
• Sort a copy of the entire array by Y just once, before the first recursive call
• Before each recursive call, distribute the Y-sorted array into left and right parts, BL and BR
• By distributing the array B, which is ordered, BL and BR end up being in order, without needing to be sorted!


• BL must be a copy of exactly the same points as in A(l   .. mid)
• BR must be a copy of exactly the same points as A(mid+1 .. r)
• What we want: BL := sort_by_y(copy(A(1 .. mid)))
• But this is not efficient
• Instead we distribute the points of B, which is in order, into either BL or BR, leaving them in order


• Remember to Assume the recursive call works!

   procedure closest(A: PointArray, d: out Float)
      B: Array(A'range) of Point := A
      sort_by_X(A)
      sort_by_Y(B)
      cp_helper(A, 1, A'last, B, d)

   procedure cp_helper(A: PointArray, l, r: Positive; B: PointArray; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         -- Divide: Distribute, then solve left and right halves
         mid := (l + r) / 2

         BL: Array(l .. mid) of Point
         BR: Array(mid + 1 .. r) of Point
   

                distribute(B, A(mid).x, BL, BR)

         procedure distribute(B: PointArray; midx: Natural; BL, BR: out PointArray)
            countL, countR: Natural := 0
            for p of B
               if p.x < midx then
                  countL := @ + 1
                  BL(countL) := p
               else
                  countR := @ + 1
                  BR(countR) := p
             assert(countL + countR = r - l + 1 AND countL = countR)
      

         cp_helper(A, l,       mid, BL, dL)
         cp_helper(A, mid + 1, r,   BR, dR)

         d = min(dL, dR)

         -- Conquer: Only compare points close to midline

         -- Fill slab with all points OF B within d
         --   left or right of mid of A
         slab : Array (l .. r) of Point
         fillslab(B, d, slab, slab_size)

         --  Search the slab
         for i in l .. slab_size - 1 loop
            for j in i + 1 .. min(i + 7, slabsize)  loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis

      

• Complexity:
• Distribute B into BL and BR: $O(n)$
• Two recursive calls: $2T(n/2)$
• Fill slab: $O(n)$
• Check slab: $O(7n)$


• $T(n) = 2T(n/2) + \Theta(9n) = n\lg n$

What Can go Wrong?

• What if two points have the same X location and that location is on a Mid?
• A: You must put the points into the correct BL or BR, depending on where they are in A


• How to solve this problem:
• Distribute very carefully, or
• B: Array(1 .. n) of Triples!
• The triple p contains x, y, and p's location in the sorted A

Array of Triples

   type Triple is ...
   type TripleArray is ...

   procedure closest(A: PointArray, d: out Float)
      sort_by_X(A)
      B: TripleArray(A'range) of Triple := makeTripleArray(A)
      sort_by_Y(B)
      cp_helper(A, 1, A'last, B, d)

   procedure cp_helper(A: PointArray, l, r: Positive; B: TripleArray; d: out Float)
      if r - l < 3 then --  Handle the base case
      else
         -- Divide: find closest of left and right halves
         mid := (l + r)/2

         BL: Array(l .. mid) of Point
         BR: Array(mid+1 .. r) of Point

         distribute(B, A(mid).x, BL, BR)

         cp_helper(A, BL, l,     mid, dL)
         cp_helper(A, BR mid+1, r,   dR)

         d = min(dL, dR)

         -- Conquer: Only compare points close to midline

         -- Fill slab with all points OF B within d left or right of mid of A
         slab : Array (l .. r) of Point
         fillslab(B, d, slab, slab_size)

         --  Search the slab
         for i in l .. slab_size - 7 loop
            for j in i + 1 .. min(i + 7, slabsize) loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis

      

Merge Solution

• Another way: Merge sort:
• each recursive call changes x-sorted to a y-sorted list
• After the recursive calls, the two y-sorted lists of size n/2 are merged to a y-sorted list of size n

   type Pair is ...
   type PairArray is ...

   procedure closest(A: PointArray, d: out Float)
      sort_by_X(A)
      cp_helper(A, 1, A'last, d)

   procedure cp_helper(A: PointArray, l, r: Positive; d: out Float)
         with pre => is_X_sorted(A(l .. r)), post => is_Y_sorted(A(l .. r))

      if r - l < 3 then --  Handle the base case
      else
         -- Divide: find closest of left and right halves
         mid := (l + r)/2

         cp_helper(A, l,     mid, dL)
         cp_helper(A, mid+1, r,   dR)
         -- AL and AR are each now in Y-sorted order!

         d = min(dL, dR)

         B: PointArray(l .. r)    -- Space for merge

         merge(A(l .. mid), A(mid + 1 .. r), B)
         -- Merge AL and AR into B, in Y-sorted order!

         A(l .. r) := B   
         -- when we're done, A must be sorted!

         -- Conquer: Only compare points close to midline

         -- Fill slab with all points OF A within d left or right of mid of A
         -- TODO: Hmmm, we'd better calculate and pass midx to fillslab 
         slab : Array (l .. r) of Point
         fillslab(A, d, slab, slab_size)

         --  Search the slab
         for i in l .. slab_size - 7 loop
            for j in i + 1 .. min(i + 7, slabsize) loop
               tmp_dis = distance(slab(i), slab(j)) 
               tmp_dis = distance(A(i), A(j)) 
               if tmp_dis < d then
                  d := tmp_dis

      

• Theoretical Performance: $T(n) = 2T(n/2) + \Theta(n)=\Theta(n \lg n)$
• But ... the assignment checker???
• Don't use this unless you can make it n lg n

And We're Done!