# Weighted closest-pair-of-points problem

I want to solve the following optimisation problem (an approximation or heuristic would be helpful as well).

I have two sets of points in the plane:
$P=\left\{ p_{1},p_{2},\dots,p_{N}\right\}$ and $Q=\left\{ q_{1},q_{2},\dots,q_{M}\right\}$.

Each point has a value/weight, i.e there is a function $v:P\cup Q\rightarrow\mathbb{R}$ which assigns a weight/value to each point.

I want to find a pair of points, one in $P$ and one in $Q$, which maximizes the function $f:P\times Q\rightarrow \mathbb{R}$,
$f\left(p,q\right)=v\left(p\right)+v\left(q\right)-\alpha\cdot d\left(p,q\right)$

where $d$ is the Manhattan distance between the points. So on one hand I want points with a high value, on the other hand I want them to be close together. $\alpha$ is just a constant which I will determine according to the specific application and it assigns a "relative importance" of the distance relative to the weights.

The choice for Manhattan distance is because I want to use this with actual geographical points which are connected by actual road networks.

Does anyone have a suggestion how to solve this without an exhaustive search over all pairs? even an approximation or heuristic suggestion would be helpful.

Thank you!

• Have you tried adapting the standard $O(n \lg n)$ divide-and-conquer algorithm for finding closest-pair-of-points? Normally you choose a vertical line that splits the set of points in two, solve the left half and the right half recursively, and then the merge step looks at all points that fall within some distance on either side of the splitting line. You might be able to use the values on the points to adjust the distance (to adjust how wide a strip around the splitting line you look at), in the merge step. Have you tried to see if you can make something like this work? – D.W. Apr 6 '15 at 18:18
• I have attempted it on paper. It seems like one possibility. The value of $\alpha$ will influence the size of the strip. It seems that for a realisitic value of $\alpha$, the strip is too wide, so that pretty much all pairs are considered. I guess the width could be lowered to get a heuristic solution, so that's one possibility – fiftyeight Apr 8 '15 at 4:59
• You could also look at adapting a $O(n)$ time randomized algorithm for finding the closest pair of points, e.g., these slides or Rabin's algorithm. I don't know if that will lead anywhere useful or not. It might depend upon the degree of variation in $v(\cdot)$ values, compared to the size of $\alpha$ and $\min d(p,q)$. – D.W. Apr 8 '15 at 6:25
• I'll check those links, thank you. If it helps, I anticipate that $v(\cdot)$ would span from $1$ to $10$ with the probability of smaller values being higher. So the probability distribution looks somwhat like $1\over{v}$. I also anticipate $\alpha \approx {1\over5}$ where the distance is measured in kilometers. And the points themselves will lie in a square of up to $50 \times 50$ in size (in kilometers) – fiftyeight Apr 8 '15 at 14:41
• Hi @fiftyeight, did my answer below meet your needs? Was it what you were looking for? – D.W. Apr 16 '15 at 9:44

## Suggested solution: use a k-d tree

I recommend you use a k-d tree to store the points of $Q$, with $k=3$ (3 dimensions). The three dimensions of each point are the $x$-coordinate, the $y$-coordinate, and the $v$-value (the value/weight of the point).

Store all of the points of $Q$ in a k-d tree. Then given a point $p$, you can use the structure of the $k$-d tree to find the point $q \in Q$ that maximizes $f(p,q)$ in a relatively efficient way -- more efficient than naively testing each possible point $q \in Q$. Then, you can iterate through the points of $P$, and for each point $p \in P$, use this procedure to find the point $q \in Q$ that maximizes $f(p,q)$. Keep track of the best pair you've seen, and this will solve your problem.

I guess I'd better describe how to use the k-d tree. Recall that each node in the k-d tree splits by a binary condition of the form $x \le c$ or $y \le c$ or $v \le c$, where $c$ is some constant (all points that satisfy the condition are stored in the left subtree of that node and all points that don't satisfy the condition are stored in the right subtree). Each point $q \in Q$ is stored somewhere in the tree. So, any node $s$ in the k-d tree corresponds to a subset $S_s \subseteq Q$ of points from $Q$, namely, the points that appear somewhere in the subtree rooted at $s$. Due to the splitting criteria, each set $S_s$ has the form $Q \cap R$ where $R$ is some axis-aligned 3-dimensional rectangle. For a particular node $s$, define

\begin{align*} s.x_\text{max} &= \max\{q.x : q \in S_s\}\\ s.x_\text{min} &= \min\{q.x : q \in S_s\}\\ s.y_\text{max} &= \max\{q.y : q \in S_s\}\\ s.y_\text{min} &= \min\{q.y : q \in S_s\}\\ s.v_\text{max} &= \max\{v(q) : q \in S_s\}\\ s.v_\text{min} &= \min\{v(q) : q \in S_s\} \end{align*}

As you form the k-d tree, precalculate the above 6 values for every node and store them in that node of the k-d tree, so they're available for ready lookup.

Define $R_s = [s.x_\text{min},s.x_\text{max}] \times [s.y_\text{min},s.y_\text{max}] \times [s.v_\text{min},s.v_\text{max}]$. We can see that $S_s = Q \cap R_s$. It will be helpful to define

$$d(p,R_s) = \min \{d(p,r) : r \in R_s\}.$$

Note that you can efficiently calculate $d(p,R_s)$, given the point $p$ and the 6 values associated with $s$. (Basically, this involves a case analysis through 27 cases, depending upon where $p$ is relative to $R_s$. The running time to compute $d(p,R_S)$ is $O(1)$.)

Let's fix a point $p$. I'm going to describe an algorithm to find the point $q \in Q$ that maximizes $f(p,q)$, by recursively searching the k-d tree for $Q$ in a suitable way. The idea is that we'll recursively visit the nodes of the k-d tree, but pruning the traversal when we reach a node that cannot possibly hold any point $q \in Q$ that will be better than the best seen so far.

To help us prune the traversal, notice that we can bound the maximum possible value of $f(p,q)$ over all $q \in S_s$, as follows:

$$f(p,q) \le v(p) + s.v_\text{max} - \alpha \cdot d(p,R_s) \text{ for all } q \in S_s.$$

Thus, whenever our recursive traversal reaches a node $s$ of the k-d tree, we'll use the above bound to figure out whether there is any possibility that exploring the subtree rooted at $s$ could possibly yield some improvement over the best result seen so far. If not, we'll prune the traversal and won't visit any of the children or descendants of $s$.

Thus, the algorithm looks something like the following:

def findbestpoint(p):
bestsofar := -infinity
visit(p, root of k-d tree for Q)
return bestsofar

def visit(p, s):
if v(p) + s.v_max - alpha * d(p,R_s) <= bestsofar:
return
update bestsofar based upon the point q stored in node s (if any)
if s is a leaf:
return
visit(p, s.rightchild)
visit(p, s.leftchild)


Given a k-d tree for $Q$ and a point $p \in P$, this algorithm computes $\max \{f(p,q) : q \in Q\}$. It will visit some of the nodes of the k-d tree but not all of them. So, iterate through the points $p \in P$ and call findbestpoint() on each one.

The above is correct but not necessarily optimal in performance. As an optimization, I suggest that at nodes that split on the $x$-coordinate or $y$-coordinate, you should choose the order of the recursive traversal to prefer the subtree that contains $p$; and for nodes that split on $v$-value, always prefer the subtree with the larger $v$-value. In other words, I suggest you replace the last two lines of visit(p, s) with

    if the condition associated with s is a split on the v-value:
visit(p, s.rightchild)
visit(p, s.leftchild)
else if p satisfies the condition associated with s (i.e., p is contained within the subtree rooted at s.leftchild):
visit(p, s.leftchild)
visit(p, s.rightchild)
else:
visit(p, s.rightchild)
visit(p, s.leftchild)


This will choose a traversal ordering that is more likely to help you prune many unpromising nodes.

As another optimization, rather than resetting bestsofar to $-\infty$ on each call to findbestpoint, I suggest you initialize it to $-\infty$ once and then don't reset it again. This will enable you to prune more effectively, and since you only want to find one pair $p,q$ that maximizes $f(p,q)$, this optimization is valid.

As yet another optimization, you can choose the order in which you iterate through the points of $P$ based on some heuristic. For instance, you could try sorting $P$ by value and iterating through $P$ in order of decreasing value; or you could use value minus $\alpha$ times the distance to the centroid of $Q$, or something like that. Basically, you want to first try values of $P$ that are more likely to yield you a large $f(\cdot,\cdot)$ value, as that will make the pruning more effective.

I have no theoretical guarantees on the asymptotic running time of this algorithm. However, when the set $Q$ is large, I expect this will give significant improvements. Once you understand how k-d trees, it should be relatively easy to code this up and give it a try.

Note that you can always choose whether to store $Q$ in the k-d tree and iterate through $P$ (as outlined above) or whether to store $P$ in the k-d tree and iterate through $Q$. I suggest you store the larger set in the k-d tree and iterate through the smaller set.

## A better algorithm using k-d trees

The idea above can be generalized further, to obtain an algorithm that might be even more efficient, at the cost of making the algorithm more complicated and more tedious to implement.

The algorithm above explores pairs $p,s$ where $p$ is a point and $s$ is a node in the tree (summarizing some subset of points from $Q$). In the refinement, we'll explore pairs $s,t$ where $s$ is a node in the k-d tree for $P$ and $t$ is a node in the k-d tree for $Q$ (so $s$ summarizes some subset of $P$ and $t$ summarizes some subset of $Q$). This will expose further opportunities for pruning.

Assume we have built a k-d tree for $P$ and a k-d tree for $Q$, augmented with the information listed above. For simplicity of exposition, assume we're using one of the variants of the k-d tree data structure where the points of $P,Q$ are stored in the leaves of the tree; the internal nodes don't store points. (You can generalize the ideas here to the case where each internal node also stores a point from $P$ or $Q$; the extension is messy but not conceptually difficult.)

Also, if $s$ is a node in the k-d tree for $P$ and $t$ is a node in the k-d tree for $Q$, define

$$d(R_s,R_t) = \min \{d(r_1,r_2) : r_1 \in R_s, r_2 \in R_t\}$$

to be the distance between the two rectangles $R_s,R_t$. This is the distance between the "closest point of approach" between the two rectangles. For example, if the two rectangles $R_s,R_t$ overlap, $d(R_s,R_t)=0$. This distance is a property solely of the rectangles $R_s,R_t$ and doesn't have anything to do with the points $p,q$ that are contained within them. Given the descriptions of $R_s,R_t$ (i.e., the 6 values stored in the nodes $s,t$), you can compute $d(R_s,R_t)$ in $O(1)$ time via an ugly but straightforward case analysis.

Now notice that

$$f(p,q) \le s.v_\text{max} + t.v_\text{max} - \alpha d(R_r,R_s)$$

for all $p \in R_s, q \in R_t$. This gives us an upper bound on the value of $f(p,q)$ over all class of possible pairs $p,q$ of points. Conceptually, the node $s$ summarizes a bunch of points $p$ (the set $S_s \subseteq P$), and the node $t$ summarizes a bunch of points $q$ (the set $S_t \subseteq Q$); the above bound lets us upper-bound the best possible $f(\cdot,cdot)$ value, among all pairs of points $p,q$ in these sets.

Based on this insight, we can get a recursive algorithm to find the best pair of points:

def findbestpair():
bestsofar := -infinity
visitpair(root of k-d tree for P, root of k-d tree for Q)
return bestsofar

def visitpair(s, t):
if s.v_max + t.v_max - alpha * d(R_s,R_t) <= bestsofar:
return
if s is a leaf and t is a leaf:
let p := the point contained within s
let q := the point contained within t
if f(p,q) > bestsofar:
bestsofar := f(p,q)
else if s is a leaf:
visit(s, t.rightchild)
visit(s, t.leftchild)
else if t is a leaf:
visit(s.rightchild, t)
visit(s.leftchild, t)
else:
visit(s.rightchild, t.rightchild)
visit(s.leftchild, t.leftchild)
visit(s.leftchild, t.rightchild)
visit(s.leftchild, t.leftchild)


We can see that this explores pairs $s,t$ (corresponding to a set $S_s \times S_t$ of points), pruning when we can prove that none of the pairs $(p,q) \in S_s \times S_t$ will have a larger value of $f(p,q)$ than the best we've seen so far.

You can apply similar optimizations to those mentioned above (e.g., when $s$ is a leaf that splits on $x$-coordinate or $y$-coordinate, first try the child of t that $s$ is contained in; when $s$ is a leaf that splits on $v$-value, first try the child with the larger $v$-value; and symmetrically when $t$ is a leaf).

I would conjecture that if the size of the both sets is large enough, this might perform better than the solution mentioned above. However, I have no proof, and the only way to know for sure is to implement it and give it a try on your data set.

## A connection to computational geometry

I also noticed that your problem is related to the following classical problem:

• Given a set of discs in the 2-D plane, determine whether any pair of discs overlap.

and the following version of the problem:

• Given a set of discs in the 2-D plane, determine the pairs of discs that are closest (where the distance between two discs is given by the closest point of approach, i.e., the minimum distance between any point in the first disc and any point in the second disc).

These are classic problems in computational geometry which I think have sub-quadratic time algorithms (by sweepline methods or something, if I recall correctly). Unfortunately your problem seems harder those questions about discs, so I don't know whether those methods will transfer over.

Here is the connection. Define $w:P \cup Q \to \mathbb{R}$ by $w(p) = v(p)/\alpha$ and $w(q) = v(q)/\alpha$. Define $g:P \times Q \to \mathbb{R}$ by

$$g(p,q) = d(p,q) - w(p) - w(q).$$

Your problem is equivalent to asking for the pair of points $p \in P,q \in Q$ that minimizes $g(p,q)$, as $f(p,q) = - \alpha \cdot g(p,q)$.

Next, treat each point $p$ as being associated with a disc centered at $p$ and having radius $w(p)$; and similarly for each point $q$. In this way if all of the $w(\cdot)$ values are positive, we obtain a collection of discs in the two-dimensional plane. Notice that if two discs $p,q$ don't overlap, then $g(p,q)$ denotes exactly the distance between those two discs. This condition is equivalent to requiring that $g(p,q)$ be positive for all $p,q$.

In summary, if we were promised that $v(p)$ was positive for all $p,q$ and that the maximum value of $f(p,q)$ was negative, then your problem would be equivalent to finding the closest pair of discs in the 2-D plane.

Of course, in your problem, we're not promised any of those things. Therefore, your problem is harder. However, if you were really eager on trying to find the best possible algorithm, you could explore the literature to find algorithms for finding the closest pair of discs and see if any of them can be generalized to your setting.