DBSCAN "region query" too slow; implement a tree? My current DBSCAN in Python works...but its indexing is far too slow; its a linear scan:

def find(self, p):
    nearby = []
    for pt in self.data:
        # Haversine algorithm
        px, ptx, py, pty = map(math.radians, [p.x, pt.x, p.y, pt.y])
        dx = px - ptx
        dy = py - pty
        a = math.sin(dy/2)**2 + math.cos(py)*math.cos(pty)*math.sin(dx/2)**2
        c = 2*math.asin(math.sqrt(a))
        d = 6371.0*c
        if d &lt; self.eps:
            nearby.append(pt)
    return nearby


self.eps: The defined radius in km (hence the Haversine implementation)

self.data: Python list object of my Point() objects:

class Point:
    def __init__(self, lon, lat):
        self.x = lon
        self.y = lat
        self.cluster = None
        self.neighbors = 0


When I run the DBSCAN with say 27,000 points, the find() has to iterate through each point (technically minus 1 each time, but that's moot) for at worst O(N2) (assumes no points in the 27,000 are noise (i.e. they have valid minimum points within self.eps).

I have heard maybe a tree would be better for this? Not sure which and I would prefer to not keep adding modules (if possible!).
You don't provide much information, which is unfortunate, since this kind of optimization problem relies on knowing how to "skew" the code in order to get better performance.

For example, are your points constant? If you are always searching the same set of 27k points - for example, cities or road termini or railroad stations or whatever - then you can pre-compute lots of data about those points to make the code faster.

Are the points accumulated gradually over time? Could you perform an expensive computation on each point safely, knowing that the various points are collected at some interval that allows for that?

Are the points close together or far apart, relative to the likely neighborhood distances? That is, can you estimate the number of points that will be in the likely neighborhood queries?

Is your program memory-constrained? Can we trade memory allocation for performance? 

Are the eps values - the neighborhood distances - truly variable? Or is there a static list (1km, 10km, 25km, 50km) that the user chooses from?

All these questions, and more, could be used to tune the code for speed. When you're asking a speed question, there's no such thing as too much data!

Finally, what flavor of python are you using? PyPy? Cpython? Jython? Cython? 

With that out of the way, let me make a few "low hanging fruit" suggestions. I'll assume you're using "the python that came with my system" - meaning Cpython 2.7 or 3.4+.

The cpython compiler isn't nearly as smart as you might like. Let's look at the generated code from your function (I put it in test.py in a class named C):

$ python2
Python 2.7.13 (default, Apr 23 2017, 16:50:35)
[GCC 4.2.1 Compatible Apple LLVM 8.0.0 (clang-800.0.42.1)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
&gt;&gt;&gt; import test
&gt;&gt;&gt; import dis
&gt;&gt;&gt; dis.dis(test.C.find)
  3           0 BUILD_LIST               0
              3 STORE_FAST               2 (nearby)

  4           6 SETUP_LOOP             230 (to 239)
              9 LOAD_FAST                0 (self)
             12 LOAD_ATTR                0 (data)
             15 GET_ITER
        &gt;&gt;   16 FOR_ITER               219 (to 238)
             19 STORE_FAST               3 (pt)

  6          22 LOAD_GLOBAL              1 (map)
             25 LOAD_GLOBAL              2 (math)
             28 LOAD_ATTR                3 (radians)
             31 LOAD_FAST                1 (p)
             34 LOAD_ATTR                4 (x)
             37 LOAD_FAST                3 (pt)
             40 LOAD_ATTR                4 (x)
             43 LOAD_FAST                1 (p)
             46 LOAD_ATTR                5 (y)
             49 LOAD_FAST                3 (pt)
             52 LOAD_ATTR                5 (y)
             55 BUILD_LIST               4
             58 CALL_FUNCTION            2
             61 UNPACK_SEQUENCE          4
             64 STORE_FAST               4 (px)
             67 STORE_FAST               5 (ptx)
             70 STORE_FAST               6 (py)
             73 STORE_FAST               7 (pty)

  7          76 LOAD_FAST                4 (px)
             79 LOAD_FAST                5 (ptx)
             82 BINARY_SUBTRACT
             83 STORE_FAST               8 (dx)

  8          86 LOAD_FAST                6 (py)
             89 LOAD_FAST                7 (pty)
             92 BINARY_SUBTRACT
             93 STORE_FAST               9 (dy)

  9          96 LOAD_GLOBAL              2 (math)
             99 LOAD_ATTR                6 (sin)
            102 LOAD_FAST                9 (dy)
            105 LOAD_CONST               1 (2)
            108 BINARY_DIVIDE
            109 CALL_FUNCTION            1
            112 LOAD_CONST               1 (2)
            115 BINARY_POWER
            116 LOAD_GLOBAL              2 (math)
            119 LOAD_ATTR                7 (cos)
            122 LOAD_FAST                6 (py)
            125 CALL_FUNCTION            1
            128 LOAD_GLOBAL              2 (math)
            131 LOAD_ATTR                7 (cos)
            134 LOAD_FAST                7 (pty)
            137 CALL_FUNCTION            1
            140 BINARY_MULTIPLY
            141 LOAD_GLOBAL              2 (math)
            144 LOAD_ATTR                6 (sin)
            147 LOAD_FAST                8 (dx)
            150 LOAD_CONST               1 (2)
            153 BINARY_DIVIDE
            154 CALL_FUNCTION            1
            157 LOAD_CONST               1 (2)
            160 BINARY_POWER
            161 BINARY_MULTIPLY
            162 BINARY_ADD
            163 STORE_FAST              10 (a)

 10         166 LOAD_CONST               1 (2)
            169 LOAD_GLOBAL              2 (math)
            172 LOAD_ATTR                8 (asin)
            175 LOAD_GLOBAL              2 (math)
            178 LOAD_ATTR                9 (sqrt)
            181 LOAD_FAST               10 (a)
            184 CALL_FUNCTION            1
            187 CALL_FUNCTION            1
            190 BINARY_MULTIPLY
            191 STORE_FAST              11 (c)

 11         194 LOAD_CONST               2 (6371.0)
            197 LOAD_FAST               11 (c)
            200 BINARY_MULTIPLY
            201 STORE_FAST              12 (d)

 12         204 LOAD_FAST               12 (d)
            207 LOAD_FAST                0 (self)
            210 LOAD_ATTR               10 (eps)
            213 COMPARE_OP               0 (&lt;)
            216 POP_JUMP_IF_FALSE       16

 13         219 LOAD_FAST                2 (nearby)
            222 LOAD_ATTR               11 (append)
            225 LOAD_FAST                3 (pt)
            228 CALL_FUNCTION            1
            231 POP_TOP
            232 JUMP_ABSOLUTE           16
            235 JUMP_ABSOLUTE           16
        &gt;&gt;  238 POP_BLOCK

 14     &gt;&gt;  239 LOAD_FAST                2 (nearby)
            242 RETURN_VALUE


First, have a look at "line 10" - the block of code that starts with a "10" in the left-most column. That code corresponds to this python:

    c = 2*math.asin(math.sqrt(a))


The Python VM is a stack machine, so every opcode either takes something off the stack, puts something on the stack, or both. For example, the various LOAD_ ops put things onto the stack. The CALL_FUNCTION op calls a function (which is stored on the stack) with a number of arguments. And BINARY_MULTIPLY takes two arguments off the stack, multiplies them, and pushes the result back onto the stack.

There's a couple of optimization opportunities here that are good "rules of thumb." You can take this knowledge and apply it to any Cpython code you need to speed up (note: Cpython). 

First, there's that 2 * at the front of the expression. That's fairly harmless, except that on the very next line you multiply by 6371.0, and you never do anything else with the c variable - it's just a temp. So to squeeze out some performance, go ahead and move the 2* to the next line. It turns out that the compiler will go ahead and do constant math in the same expression, so you can say 6371.0 * 2 * c and it will collapse that to a single constant.

Second, note that there's this sequence of lookups all over the place:

            169 LOAD_GLOBAL              2 (math)
            172 LOAD_ATTR                8 (asin)


In general, Python believes that any operation might cause the result of a name lookup to change. Because you can return any type, and any type can overload pretty much any operation, so why shouldn't it be possible for 2 * math.asin(x) to change the function pointed to by math.sqrt? 

Because of this, and because name lookups like this really can be expensive (remember property methods and descriptors), you generally want to cache dotted-name lookups in a local variable if you're going to be looping with them:

asin = math.asin
sqrt = math.sqrt


Please note: This is not the same as doing from math import asin because that will still create a "global" name, and require a LOAD_GLOBAL operation, instead of a LOAD_FAST. Notice the FAST in the opcode name. ;-)

The same thing is true for almost every dotted name - if you are going to use it more than once, put it in a local variable.

One thing that might not be obvious is that you can cache "bound methods" - that is, method lookups like nearby.append. So you can also do:

nearby_append = nearby.append

...

nearby_append(pt)


And squeeze out another lookup.

But wait, there's more!

Because before you start doing that "peephole" optimization stuff, you need to look at improving your code in some other ways. Take a look at this:

def find(self, p):
    nearby = []
    for pt in self.data:
        # Haversine algorithm
        px, ptx, py, pty = map(math.radians, [p.x, pt.x, p.y, pt.y])


For the purposes of this function, the p object is a constant. So why do you keep doing the same conversion on it? Remember that map just applies the given function to each argument in turn. So you're calling math.radians on p.x and p.y 27000 times. Guess what? The answer isn't going to change!

Of course, I suspect that none of your points are going to change. That is, I suspect that once you create a point object, it doesn't get updated. So I think you would benefit from performing this conversion in the class initializer and storing the radian values in the object:

class Point:
    def __init__(self, lon, lat):
        self.x = lon
        self.x_radians = math.radians(lon)   # &lt;- HERE
        self.y = lat
        self.y_radians = math.radians(lat)   # &lt;- HERE
        self.cluster = None
        self.neighbors = 0


If you do that, you can eliminate that entire map call, and just retrieve the values. Yes, you pay an up-front cost of 27000 extra function calls, but it beats doing it in your find loop. (If you can't pay that cost, you could certainly do it as a separate loop at the front of your find function. You're going to do it anyway, and this lets you store the result.)

Also of note: math.cos(py) and math.cos(pty). You can certainly compute the cosine of py once at the top of the function. And if you're going to call this function a lot, it's probably worth storing the cosine(y) for all the points in the initializer. 

Heuristics, anyone?

Another suggestion is simply to guess what's a neighbor and what's not. You've got a nice calculation there, with lots of trigonometry. But keep in mind what it's doing: it's trying to compute distance on a sphere, which is basically distance along a circle. The key word is circle. 

At θ = 45°, cosθ is √2/2 as is sinθ. The sum of the two is √2, which is the maximum sum for a circle of radius 1. Thus, you know that for whatever distance N defines your neighborhood, the sum of cosine and sine will be less than or equal to √2 * N. 

If your points are widely separated - that is, if you're dealing with cities and not the entrances to buildings in the same town - you can implement a heuristic: change √2 into 2.5 to allow for spherical error, and compute the corners of a "square" centered at p with sides 2 * 2.5 * N. Then just compare pt.x and pt.y with the min/max x and y of this square. Only if the point is inside the square do you perform the "expensive" computations. This adds a set of 4 comparisons to your loop, but eliminates the trig calls for all the non-matching candidates. Hopefully, that's most of the points.

This is important! If your points are heavily clustered, there's no sense having a heuristic to guess if one point is a neighbor of another - because you know the answer is yes! In that case, it would make more sense to develop a reverse heuristic, that says "here's a square contained in the circle. If the point is inside this square you know it's a neighbor, otherwise do the computation."  The objective with heuristics it to produce a fast computation that handles 75% or more of cases. If your data is biased, make the heuristic align with the biases!

Finally, another simple optimization would be to take your self.eps value and divide it by the various factors you are presently using to multiply. This can be done outside the loop, once, and saves you some cycles.

Larger-scale changes

First, don't forget that if A is a neighbor of B, then B is a neighbor of A. Make sure you are only looping over the "higher" values in your N-squared neighbor search:

for i, p1 in enumerate(points):
    for p2 in points[i+1:]:
        # etc.


You just have to be sure and add A to B's neighborhood and B to A's neighborhood at the same time.

Depending on what you're doing, there are a couple of geospatial data structures that you might profitably use. Two obvious ones are a hash of some sort and a K-d tree. The hash approach will work if your points are clustered and if you can determine a key generation algorithm that will reflect this. For example, if you know that your points are GPS coordinates in major metropolitan areas (New York, Los Angeles, Dallas, Chicago) then a hash key that "rounded" the coordinates to the nearest city center might be adequate. (Or truncated the lat/lon to, say, the nearest 5 degrees.) In this case, you'd want to have a hash value with lots of collisions. The collisions would be the likely neighbors. So you would iterate over the hash keys (equal to, say, 10% of the points) and identify which buckets you would search in. On the other hand, if the points are more distributed, or the clustering isn't obvious, then go with the K-d tree.

The K-d tree is a binary tree (or binary-search-tree if you're a child of the 21st century) that can be used for locating any point. To find neighboring points using your distance function, you would simply do the distance computation at each node of the tree. The special case would be when a node's key value is within the heuristic neighborhood of your target point, you'll have to traverse both children of the node, to find any other nodes in the other half of the tree that might be in the neighborhood. For example, if a node is a lat node, with lat=A, and your target point has lat=A-2km, then for neighborhood=5km you'll have to search on either side of that node.

Also, note that you'll need to develop a strategy for nodes with equal values. For example, if two points have the same longitude, would you store a list-of-points in the node, or make a rule that says "sort by lowest latitude and put the other points in the right hand subtree"? It's up to you, but you need to handle it. (I suggest the latter.)