# Finding a known number of circle centers that maximize the number of points within a fixed distance

I have a set of 2-D data where I want to find the centers of a specified number of centers of circles ($N$) that maximize the total number of points within a specified distance ($R$).

e.g. I have 10,000 data points $(X_i, Y_i)$ and I want to find the centers of $N=5$ circles that capture as many points as possible within a radius of $R=10$. The 5 centers and radius of 10 are given beforehand, not derived from the data.

The presence of a data point within a circle is a binary either/or proposition. If $R=10$, there's no difference in value to a point 11 units away vs. 100 units away, as they are both > 10. Similarly for being within the circle, there's no extra value to being near the center vs. near the edge. A data point is either in one of the circles or out.

Is there a good algorithm that can be used to solve this problem? These seems related to clustering techniques, but rather than minimizing the average distance, the "distance" function is 0 if the point is within $R$ of any of the $N$ points, and 1 otherwise.

My preference would be to find a way to do this in R, but any approach would be appreciated.

• Is circle overlap allowed? Apr 11, 2013 at 17:48
• This is essentially a neighborhood (or focal) operation on a raster dataset. It would be good to check out the GIS site to see if it has been answered, and to examine R packages to conduct Raster analysis. Apr 11, 2013 at 18:18
• Circle overlap is allowed, but the data points covered by both circles would not be double-counted. Thanks for the pointer to neighborhood/focal operation on raster datasets. I will search for something along those lines. Apr 11, 2013 at 19:06
• @Andy W Although focal operations would naturally be involved in a solution, this question is beyond the expertise of the GIS community, IMHO, because it is really a (pretty hard) optimization problem. It's not a straightforward find-the-max-of-a-focal-mean-grid. I would recommend keeping it here for a while and then, if no satisfactory solution emerges, migrating to a programming-oriented site.
– whuber
Apr 12, 2013 at 20:49
• ....or migrating to math.overflow? They might have some insights on this too. Apr 13, 2013 at 5:40

This is a variation k-means problem. The radius of the centers doesn't matter, as long as their are assumed equal.

It will put the centers of the circles at locations of highest probability of the points.

Classic K-means Procedure:

1. set cluster count to 5
2. put each point in a random cluster
3. for each cluster, calculate the mean position
4. for each point, calculate the distance to each new mean position
5. associate membership with the nearest cluster
6. repeat until done (iterations, change in position, or other error metric)

Options:

• You can use some under-relaxation after 3, where you translate the mean position slowly toward the new position.
• this is a discrete system so it doesn't converge perfectly. Sometimes it does and you can end when points stop changing membership, but sometimes they just wiggle a bit.
• If you are making your own code (as most folks should) then you can use the POR k-means above as a starting point, and do some variation on EM informed by percent of points exclusively and completely encompassed by the circles.

Why K-means attacks the problem:

• It is the equivalent of fitting a Gaussian Mixture Model where the covariances of the components are equal. The centers of the mixture components are going to be located at the positions of highest expectation of points. The curves of constant probability are going to be circles. This is EM algorithm so it has the asymptotic convergence. The memberships are hard, not soft.
• I think that if the fundamental assumption of the equal variance components mixture model is reasonably "close", whatever that means, then this method is going to fit. If you just randomly distribute points, it is less likely to fit well.

There should be some analog of a "Zero Inflated Poisson" where there is a component that is non-gaussian that picks up the uniform distribution.

If you wanted to "tune" you model and were confident that there were enough sample points then you could initialize with the k-means, and then make an augmented k-means adjuster that removes points outside of the radii of the circles from competition. It would slightly perturb the circles you have, but it might have slightly improved performance given the data.

• Could you please be a little more explicit about how K-means solves this problem?
– whuber
Aug 7, 2013 at 20:40
• Thanks for the suggestion. It's still not clear to me that the K-means approach solves the problem? Consider the example of three clusters of normal (0,1) generated data, where the centers are offset by 5 units or so. The K-means centers would give the max density. Now cut out some of the points with "holes" such that data closer than 0.5 to the centers are removed. K-means will still show about the same centers, but if you've trying to get max coverage for N=3, R=0.5 that's clearly not the right answer (because the doughnut holes contain no data). Am I misunderstanding something? Aug 26, 2013 at 17:31
• Will look into your question more for a better answer when I have time. I like to allow negative weights. The can sometimes handle data doughnuts as well as radial rational polynomials. Aug 27, 2013 at 17:15

Someone probably has a better formal algorithm, but here's one brute force approach (a hack?). I'd use one of the hexagonal binning algorithms to compute a 2D histogram. Like hexbin in R.

I'd use a hexagon size that'd roughly circumscribe your circle of radius R and then sort on the top N bins. If you got N distinct far away bins, great. Now one way is to move about the circle locally on a 2*R scale (in x and y directions) from the center of the top density hexagons. Computing densities can roughly optimize position locally. This will account for the fact that the hexagons were not a moving window with respect to a fixed origin.

If all top bins are close by you'd have to have some smarter way of moving your circles in that vicinity.

Note that I can think of several corner cases where such a naive strategy will spectacularly fail. Yet, just a starting point.

Meanwhile, I hope someone has a better algorithm.

• Something like this might solve the problem, at least approximately, for one circle. (This can easily be done using focal counts with a GIS.) But it won't solve the multiple-circle problem.
– whuber
Apr 12, 2013 at 16:56
• @whuber: What about solving for one circle then dropping all points that lie within that circle and then repeating the original algorithm? Can you see situations where this would fail? Apr 12, 2013 at 18:23
• Yes, easily. (Yours is a "greedy algorithm.") Consider the case $R=10, N=2$ in one dimension with points at $0,1,2,20,21,28,29,30,31,32,39,40$. Your algorithm puts the first circle covering $28,29,30,31,32$ and the second covering $0,1,2$: eight points in toto. A better solution covers $20,21,28,29,30$ with one circle and $30,31,32,39,40$ with another: nine points.
– whuber
Apr 12, 2013 at 20:33
• @whuber: True. You are right. Though depending on the structure of the input points in some (many?) cases the greedy and non greedy solutions might be identical or close to? I don't know. Apr 13, 2013 at 5:36
• @whuber: The problem seems mostly at boundaries. What if (somewhat like I mentioned in my answer) one moves the window +R and -R and then puts all feasible solutions on a stack and selects among them. e.g. In your 1D example on hitting 28,29,30,31,32 it'd slide the window till 18-28 and 38-48 looking for all feasable solutions. Then within these one can look for maximum point yielding combinations. Not sure if that'd help? I'm trying to see if my naive algorithm can be salvaged any? :) Apr 13, 2013 at 6:03