Hdbscan: Limitation on the number of cores actively used

[albertofavaro]

Hi!

I have run HDBSCAN on a machine with 16 CPUs and 32GB of RAM. This produces good
results, but I noticed that only 4 cores out of 16 are used for the heavy parts
of the calculation. It would be fantastic if this limitation could be removed.

Gradually increasing core_dist_n_jobs from 1 to 4 reduces the computational
time as one would expect. Setting core_dist_n_jobs to be larger than 4,
however, does not generate further improvements. Running top on the terminal
verifies that the additional cores are not employed.

I suspect that the limitation on the number of cores being used arises because
of lines 397-416 in _hdbscan_boruvka.pyx:

# A shortcut: if we have a lot of points then we can split the points
# into four piles and query them in parallel. On multicore systems
# (most systems) this amounts to a 2x-3x wall clock improvement.
if self.tree.data.shape[0] > 16384 and self.n_jobs > 1:
    datasets = [
        np.asarray(self.tree.data[0:self.num_points//4]),
        np.asarray(self.tree.data[self.num_points//4:
                                                  self.num_points//2]),
        np.asarray(self.tree.data[self.num_points//2:
                                                  3*(self.num_points//4)]),
        np.asarray(self.tree.data[3*(self.num_points//4):
                                                  self.num_points])
        ]

    knn_data = Parallel(n_jobs=self.n_jobs)(
        delayed(_core_dist_query,
                      check_pickle=False)
        (self.core_dist_tree, points,
         self.min_samples + 1)
        for points in datasets)

To test the hypothesis that the above snippet causes the problem, I have
written a Python script that resembles the code from HDBSCAN while performing
a trivial operation. Here it is:

import os
import numpy as np
import time

from sklearn.externals.joblib import Parallel, delayed

os.environ["JOBLIB_TEMP_FOLDER"] = "/tmp"

NUM_POINTS = 10**2
UP_TO = 10**8
MAX_N_JOBS = 16


def add_total(array, up_to):
    total = 0
    for i in range(up_to):
        total += 1
    return array + total


def partition_data(data):
    num_points = data.shape[0]
    datasets = [
        np.asarray(data[0:num_points//4]),
        np.asarray(data[num_points//4: num_points//2]),
        np.asarray(data[num_points//2: 3*(num_points//4)]),
        np.asarray(data[3*(num_points//4): num_points])
        ]
    return datasets


def main():
    data = np.ones(NUM_POINTS)
    datasets= partition_data(data)

    timing = []
    for n_jobs in range(1, MAX_N_JOBS+1):

        print("Trying to use {} CPUs.".format(n_jobs))
        start = time.time()
        Parallel(n_jobs=n_jobs)(
            delayed(add_total, check_pickle=False)(points, UP_TO)
            for points in datasets
            )
        end = time.time()
        timing.append(end-start)

    print(timing)

if __name__ == '__main__':
    main()

One finds that joblib manages to parallelise this script on up to 4 CPUs. As
before, adding more cores yields no apparent benefits.

If useful, I could have a go at making a pull request. The idea would be to
separate the data in n_chunks, where this number is passed as a keyword
argument, and need not be equal to 4 (as it effectively is now). As a bonus
feature, one could have that setting n_chunks=None splits the data in as
many piles as there are core_dist_n_jobs. Would this be a good strategy (@lmcinnes )?

[lmcinnes]

Great observation. I think the n_chunks approach should work well -- there is an issue with diminishing returns, but as long as you can override to avoid that (which you suggested approach does) it will be fine. I would love to see a PR for this if you have the time. Thank you!

[hammadmazhar1]

I don't supppose any progress has been made on this regard? Hitting some long runtimes with 300K points and 26ish dimensions.