Estimation of min_samples for DBSCAN I'm attempting to speed up some python code that is supposed to automatically pick the minimum samples argument in DBSCAN. Currently the execution time grows exponentially as the number of training samples increases:


0.128999948502 seconds for 100 training examples 
0.281999826431 seconds for 1000 training examples
11.5310001373 seconds for 10000 training examples



  min_samples is the number of samples (or total weight) in a
  neighborhood for a point to be considered as a core point. This
  includes the point itself.


Can anyone make some suggestions for how to speed up the code?

As a bonus, adding the n_jobs argument to DBSCAN doesn't seem to affect the execution time of DBSCAN.

import numpy as np
from sklearn.metrics.pairwise import euclidean_distances
from time import time

def autoselect_minsamps_orig(D_sym,eps,K,kmedian_thresh):
   min_samps = 10
   kmedian_dist_prev = 0
   for k in range(2,K):
      k_avg = []
      # For each row, compute median distance to k nearest neighbors
      for row in D_sym:
         flat_row = np.squeeze(np.asarray(row))
         row_nonzero = flat_row[np.flatnonzero(flat_row)]
         # Some nodes aren't connected to anything
         if len(row_nonzero) == 0: continue
         sorted_row = np.sort(row_nonzero)
         sorted_row_t = sorted_row[0:k]
         k_avg.append(np.median(sorted_row_t))
      kmedian_dist = np.median(k_avg)
      kmedian_delta = abs(kmedian_dist - kmedian_dist_prev)

      if kmedian_delta &lt;= kmedian_thresh:
         min_samps=k
         break
      kmedian_dist_prev = kmedian_dist
   return min_samps

num_training_examples = 100
num_features = 10
X = np.random.randint(5, size=(num_training_examples, num_features))

# generate a symmetric distance matrix
D = euclidean_distances(X,X)
eps = 0.25
kmedian_thresh = 0.005
K = num_training_examples

start = time()
result = autoselect_minsamps_orig(D,eps,K,kmedian_thresh)
end = time()
total_time = end - start
print('autoselect_minsamps took {} seconds for {} training examples'.format(total_time,num_training_examples))

num_training_examples = 1000
num_features = 10
X = np.random.randint(5, size=(num_training_examples, num_features))

# generate a symmetric distance matrix
D = euclidean_distances(X,X)
eps = 0.25
kmedian_thresh = 0.005
K = num_training_examples

start = time()
result = autoselect_minsamps_orig(D,eps,K,kmedian_thresh)
end = time()
total_time = end - start
print('autoselect_minsamps took {} seconds for {} training examples'.format(total_time,num_training_examples))

num_training_examples = 10000
num_features = 10
X = np.random.randint(5, size=(num_training_examples, num_features))

# generate a symmetric distance matrix
D = euclidean_distances(X,X)
eps = 0.25
kmedian_thresh = 0.005
K = num_training_examples

start = time()
min_samples = autoselect_minsamps_orig(D,eps,K,kmedian_thresh)
end = time()
total_time = end - start
print('autoselect_minsamps took {} seconds for {} training examples'.format(total_time,num_training_examples))

# min_samples is used here
db = DBSCAN(eps=eps,
            min_samples = min_samples,
            metric='precomputed',
            n_jobs=-1).fit(D)
Since the elements in each row grows, sorting the row only to take the k smallest is a bit wasteful. Using numpy.partition should always be faster:

sorted_row_t = np.partition(row_nonzero, k)[:k]


This is not needed with the examples you give, the row is already an array and flattened:

flat_row = np.squeeze(np.asarray(row))


Similarly, you can just use numpy.nonzero, no need to flatten the flattened row again:

row_nonzero = flat_row[np.nonzero(frow)]




However, this whole for loop can be done entirely in numpy:


First, mask all zeros by setting them to numpy.nan:

D = np.where(D == 0, np.nan, D)

Get the k smallest elements of each row. Since numpy.nan is not smaller than any number it will only appear at the beginning if (almost) the whole row was zeros:

k_smallest = np.partition(D, k, axis=1)[:, :k]

Compute the median of each row, and then the median of the median, ignoring numpy.nan values using numpy.nanmedian:

kmedian_dist = np.nanmedian(np.nanmedian(k_smallest, axis=1))



The first two steps can even be pulled outside of the for k in ... loop, which will save some more time. In that case we have to sort the whole array, though, because numpy.partition does not guarantee that the values are sorted, only that the k smallest values are in the slice up to k. Since K seems to be max(D_sym.shape) in all cases anyway, this should not make such a big difference. It will if the k at which the for loops stops is small enough, though, so you will have to test this.

There's no need to break and then immediately return afterwards, just return directly.

Python has an official style-guide, PEP8. It recommends using 4 spaces as indentation (not 3) and spaces after commas.



In summary, this should do the same as your code:

def autoselect_minsamps_graipher(D_sym, eps, K, kmedian_thresh):
    D = np.sort(np.where(D_sym == 0, np.nan, D_sym), axis=1)
    kmedian_dist_prev = 0
    for k in range(2, K):
        kmedian_dist = np.nanmedian(np.nanmedian(D[:,:k], axis=1))
        if abs(kmedian_dist - kmedian_dist_prev) &lt;= kmedian_thresh:
            return k
        kmedian_dist_prev = kmedian_dist
    return 10 # fall back value


This is faster than your implementation, but suffers from similar scaling behaviour:

n       autoselect_minsamps_orig    autoselect_minsamps_graipher
100     0.028 s                     0.0088 s
1000    0.32 s                      0.064 s
10000   9.5 s                       3.6 s




You could also make timing your code easier by using a decorator:

from time import perf_counter

class Timer:
    def __init__(self, name=""):
        self.name = ""
        self.start = None

    def __enter__(self):
        self.start = perf_counter()

    def __exit__(self, *args, **kwargs):
        delta = perf_counter() - self.start
        print(f"{self.name} finished in {delta:.2f} seconds")


Or, if you are stuck on Python 2:

from time import time

class Timer:
    def __init__(self, name=""):
        self.name = ""
        self.start = None

    def __enter__(self):
        self.start = time()

    def __exit__(self, *args, **kwargs):
        delta = time() - self.start
        print "{} finished in {:.2f} seconds".format(self.name, delta)


Which you can use like this:

D = euclidean_distances(X,X)
eps = 0.25
kmedian_thresh = 0.005
K = num_training_examples

with Timer(f"autoselect_minsamps_orig, K = {K}"):
    min_samples = autoselect_minsamps_orig(D, eps, K, kmedian_thresh)


Note that I used time.perf_counter, which is better suited for timing, but which is Python 3 only.