Scanpy: regress_out n_jobs question

[zhangguy]

...

Hi,

I notice that even if I set n_jobs=1 in sc.pp.regress_out, all cpus are utilized. This also happens if I set n_jobs to other numbers. Basically there isn't a noticeable difference in cpu usage no matter what number of n_jobs I set. I'm using CentOS6.8 on a machine with 28 physical cores and hyper threading on (appears as 56 cores in the os).

Is this an intended behavior, or just my installation? I understand that some numpy functions are naturally multi-threaded because of the setup of BLAS libraries.

Thanks.

[JBreunig]

I have the same issue on Ubuntu 20.04.

[ivirshup]

AFAICT, I think the parallelization you're seeing will be due to the underlying calls in statsmodels. If you turn down the number of threads blas can use, do you see the same utilization?

[zhangguy]

AFAICT, I think the parallelization you're seeing will be due to the underlying calls in statsmodels. If you turn down the number of threads blas can use, do you see the same utilization?

FYI more n_jobs seems to be slower for regress_out if I don't disable the BLAS multi threading:

sc.pp.regress_out(adata, ['percent_mito'], n_jobs=1)

regressing out ['percent_mito']
    sparse input is densified and may lead to high memory use
    finished (0:04:05)



md5-ea96e5b62c806a31e5cc55de9d263362



sc.pp.regress_out(adata, ['percent_mito'], n_jobs=24)



md5-ea96e5b62c806a31e5cc55de9d263362



regressing out ['percent_mito']
    sparse input is densified and may lead to high memory use
    finished (0:07:41)



md5-5c3a7270f838abc0df7dda06ef3f36c5



export MKL_NUM_THREADS=1
export NUMEXPR_NUM_THREADS=1
export OMP_NUM_THREADS=1



md5-1882d585f4bc7e45f82a4880b38469ac



sc.pp.regress_out(adata, ['percent_mito'], n_jobs=24)



md5-ea96e5b62c806a31e5cc55de9d263362



regressing out ['percent_mito']
    sparse input is densified and may lead to high memory use
    finished (0:00:23)

[ivirshup]

More interesting is that regress_out becomes lightning fast when n_jobs = 24 and with BLAS multi threading disabled:

Thats not too surprising to me. This must be significantly over scheduling your machine.

---

This got me doing a little more digging into this, and it look's like there's actually a solution now! We can use threadpoolctl to dynamically manage the number of threads BLAS uses via the threadpool_limits context manager.

I'm definitely interested in using this inside scanpy to manage the number of threads used here. Not quite sure yet what the right behaviour/ api is. Some options:

• Should all calls use 1 blas thread by default, so parallelization will only happen through the number of jobs?
• Do we only limit the number of threads if n_jobs is specified?
• Do we try and be fancy, with something like n_threads = n_cpus // n_jobs?

Minor update

If we use joblib (with the loky backend) instead of multiprocessing, the fancy solution will be used by default.

I think this is what the code would look like inside of regress_out:

    from joblib import Parallel, delayed
    res = Parallel(n_jobs=n_jobs)(delayed(_regress_out_chunk)(task) for task in tasks)

[JBreunig]

Is there a general fix for this besides zhangguy's great suggestion? (It works but as was stated, but I'm not sure how this alters other functions in the scanpy). I continue to get the putative over scheduling and sometimes a crash (on big datasets) when using all cores versus the super fast completion when using 1/2 cores) on a 32 core/64 thread threadripper. (RAM doesn't seem to be a problem as I'm barely touching 10% of the 256GB.)