# Selecting most points from a set of points with distance constraint

I am looking for an algorithm to select the largest subset of $$M$$ points from a set of $$N$$ points ($$M < N$$) such that no point is within a certain minimal distance d to any other point in $$M$$? I assume the points are in 3D space. I could exhaustively try every combination but I am wondering if there is a better method.

This can be stated as a formal optimization problem:

\begin{aligned} \max& \sum_i \color{darkred}x_i \\ & \mathit{\color{darkblue}{dist}}_{i,j} \ge \mathit{\color{darkblue}{mindist}}\cdot \color{darkred}x_i \cdot \color{darkred}x_j && \forall i\lt j\\ & \color{darkred}x_i \in \{0,1\} \end{aligned}

(red=decision variables, blue=constants) This is an MIQCP model: it has a (convex) quadratic constraint. It can be solved with solvers like Gurobi and Cplex. We can also easily linearize the model as:

\begin{aligned} \max& \sum_i \color{darkred}x_i \\ & \mathit{\color{darkblue}{dist}}_{i,j} \ge \mathit{\color{darkblue}{mindist}}\cdot (\color{darkred}x_i+\color{darkred}x_j-1) && \forall i\lt j\\ & \color{darkred}x_i \in \{0,1\} \end{aligned}

This is now a linear MIP model. A linear alternative is:

\begin{aligned} \max& \sum_i \color{darkred}x_i \\ & \color{darkred}x_i+\color{darkred}x_j \le 1 && \forall i\lt j|\mathit{\color{darkblue}{dist}}_{i,j} \lt \mathit{\color{darkblue}{mindist}} \\ & \color{darkred}x_i \in \{0,1\} \end{aligned}

The interpretation of the model is: $$\color{darkred}x_i$$ and $$\color{darkred}x_j$$ cannot be both in the solution if they are too close to each other. This model generates fewer constraints.

The linear MIP models can be solved with any MIP solver. Good solvers solve them very quickly, even for large $$N$$.

For a small 2d data set with minimum distance of 50, a solution can look like:

I could solve the linear model using a random data set with $$N=500$$ points to proven optimality in less than a second. It was solved during preprocessing, so zero branch&bound nodes were needed. $$N=1000$$ took 3 seconds.

Notes

• the solution is, in general, not unique. The $$N=50$$ dataset has 29 different, optimal solutions.
• we only need the upper-triangular part of the distance matrix
• this model works with 2d, 3d or other dimensional coordinates (it only needs the distances)
• depending on the solver, the linear formulation can be much faster than the quadratic formulation.