Is there a way to make this Python kNN function more efficient?

Question

After having troubles with MATLAB I decided to try Python:

I wrote a function that calculates kNN when the samples are of my own class using my own distance function:

def closestK(sample, otherSamples, distFunc, k):
"Returns the closest k samples to sample based on distFunc"
    n = len(otherSamples)
    d = [distFunc(sample, otherSamples[i]) for i in range(0,n)]
    idx  = sorted(range(0,len(d)), key=lambda k: d[k])
    return idx[1:(k+1)]

def kNN(samples, distFunc, k):
    return [[closestK(samples[i], samples, distFunc, k)] for i in range(len(samples))]

and this is the distance function:

@staticmethod    
def distanceRepr(c1, c2):
    r1 = c1.repr
    r2 = c2.repr
    # because cdist needs 2D array
    if r1.ndim == 1:
        r1 = np.vstack([r1,r1])
    if r2.ndim == 1:
        r2 = np.vstack([r2,r2])

    return scipy.spatial.distance.cdist(r1, r2, 'euclidean').min()

But it still works amazingly slower compared to the "normal" kNN function, even when using "brute" algorithm. Am I doing something wrong?

UPDATE

I'm adding the constructor of the class. The attribute repr contains a set of vectors (from 1 to whatever) and the distance is calculated to be the minimal euclidean distance between the two sets of repr.

class myCluster:
    def __init__(self, index = -1, P = np.array([])):
        if index ==-1 :
            self.repr = np.array([])
            self.IDs = np.array([])
            self.n = 0
            self.center = np.array([])
        else:
            self.repr = np.array(P)
            self.IDs = np.array(index)
            self.n = 1
            self.center = np.array(P)

and the rest of relevant code (X is a matrix whose rows are samples and columns are variables):

level = [myCluster(i, X[i,:]) for i in range(0,n)]
kNN(level, myCluster.distanceRepr, 3)

UPDATE 2

I've made some measurements and the line that takes most of the time is

d = [distFunc(sample, otherSamples[i]) for i in range(0,n)]

So there is something with the distFunc. When I change it to return

np.linalg.norm(c1.repr-c2.repr)

i.e. "normal" vector calculation, with no sorting, the running time stays the same. So the problem lies in the calling of this function. Does it make sense that the use of classes changes the running time by a factor of 60?

Answer 1

You're just running into the slowness of Python (or rather the CPython interpreter I should say I guess). From wikipedia:

NumPy targets the CPython reference implementation of Python, which is a non-optimizing bytecode compiler/interpreter. Mathematical algorithms written for this version of Python often run much slower than compiled equivalents. NumPy seeks to address this problem by providing multidimensional arrays and functions and operators that operate efficiently on arrays. Thus any algorithm that can be expressed primarily as operations on arrays and matrices can run almost as quickly as the equivalent C code.

And from the Scipy FAQ:

Python’s lists are efficient general-purpose containers. They support (fairly) efficient insertion, deletion, appending, and concatenation, and Python’s list comprehensions make them easy to construct and manipulate. However, they have certain limitations: they don’t support “vectorized” operations like elementwise addition and multiplication, and the fact that they can contain objects of differing types mean that Python must store type information for every element, and must execute type dispatching code when operating on each element. This also means that very few list operations can be carried out by efficient C loops – each iteration would require type checks and other Python API bookkeeping.

_{Note this doesn't concern only Python; for more background see e.g. this and this question on SO.}

Due to the overhead from the dynamic type system and the interpreter, Python would be a lot less usefull for high performance number crunching, if it wouldn't be able to tap into all sorts of compiled C and Fortran libraries (e.g. Numpy). Also, there are JIT compilers like Numba and PyPy that try to get Python code to execute closer to the speeds of statically typed, compiled code.

Bottomline: You're doing to much in plain Python relative to the work that you're offloading to fast C code. I suppose you need to adopt more like an "array oriented" coding style rather than object oriented to achieve good performance with Numpy (MATLAB is a very similar story in this regard). On the other hand, if you would use a more efficient algorithm (see the answer by Ara) then the slowness of Python might not be such an issue.

Answer 2

Here are the points I can think of:

• you are computing the distance between one sample and every other every time you call closestK so you compute the distance between every sample twice (once distance(a,b) then distance(b,a)), this can be avoided by computing it once and for all
• you recompute r (that might involve a costfull vstack) 2 * (n - 1) times where n is len(sample), you can also compute it once and for all (and store it as an attribute of myCluster ?).
• you compute a sort on the full list while you only need the top-k (no need to sort after the k-th element)
• to compute the minium distance between the points of your set you create a matrix containing every distance then take it's min: you can certainly do better

My suggestion would be to implement a top-k class with an insert method that insert only when you are better than the current k-th element (and remove it) and to modify myCluster to include r. Then your code might look like

kNN = {i : TopK() for i in xrange(len(samples))}
for i, sample1 in enumerate(samples):
    for j, sample2 in enumerate(samples[:i]):
        dist = distanceRepr(sample1, sample2)
        kNN[i].insert(j, -dist)
        kNN[j].insert(i, -dist)
return kNN

Here is a possible implementation ok TopK:

import heapq

class TopK:
    def __init__(self, k):
        self.k = k
        self.content = []

    def insert (self, key, score):
        if len(self.content) < self.k:
            heapq.heappush(self.content, (score, key))
        else:
            heapq.heappushpop(self.content, (score, key))

    def get_keys(self):
        return [elem[1] for elem in self.content]

For distanceRepr you can use something like:

import scipy.spatial

def distanceRepr(set0 ,set1):
    if len(set0) < len(set1):
        min_set = set0
        max_set = set1
    else:
        min_set = set1
        max_set = set0
    if len(min_set) == 0:
        raise Exception("Empty set")

    min_dist = scipy.inf
    tree = scipy.spatial.cKDTree(max_set)

    for point in min_set:
        distance, _ = tree.query(point, 1, 0., 2, min_dist)
        if min_dist > distance:
            min_dist = min(min_dist, distance)

    return min_dist

It will be faster than your current method for medium and big entries (lets say with sample1 and 2 with a size > 5k), it will also a much smaller memory usage, allowing it to work with big samples (where cdit simply is out of memory).