DOC: Implicit parallelism and how to disable it

[stuarteberg]

Questions:

Apparently OpenMP (or BLAS?) is used in pandas internals (or one of its dependencies), at least in some cases. For example, value_counts() uses all CPU cores on my machine (but does not benefit from it).

1. Is this expected behavior?
2. Is there any documentation on which functions are parallelized, and how to disable it if desired?

---

Details:

Here's a test program that repeatedly calls value_counts() on a large array:

# test.py
import numpy as np
import pandas as pd

a = np.random.randint(10, size=20_000_000)
s = pd.Series(a)

for _ in range(50):
    c = s.value_counts(sort=True)

From monitoring htop I discovered that value_counts() is using multiple CPU cores.

Apparently the parallelism can be controlled via OMP_NUM_THREADS, as shown here in the difference between real and user time when the number of threads is changed.

$ time OMP_NUM_THREADS=1 python test.py

real	0m8.714s
user	0m6.715s
sys	0m1.653s

$ time OMP_NUM_THREADS=10 python /test.py

real	0m7.921s
user	0m58.198s
sys	0m7.345s

This was a surprise to me. I know that some operations in np.linalg are implicitly parallelized via the BLAS implementation, but I would not have guessed that value_counts() used any of those functions. (Am I mistaken?)

This behavior is especially undesirable when one is already parallelizing one's code via multiprocessing (or, in my case, pyspark). In that case, the multiple threads in each process compete with each other, causing thrashing, etc.

Is setting the OMP_NUM_THREADS variable the recommended way to disable this behavior? Thanks in advance for any tips!

(Side note: As you can see, adding more threads apparently does not provide much benefit even in this simple case -- it just eats more CPU time without saving wall time.)

---

Appendix

This issue can be reproduced by simply installing the latest version of pandas with conda, but FWIW, here are the versions I'm using:

$ conda list | grep -E '(^python )|(^numpy)|(^pandas)|(blas)'
blas                      1.0                         mkl
numpy                     1.13.3           py36hdbf6ddf_4
openblas                  0.2.20                        7    conda-forge
pandas                    0.23.4           py36hf8a1672_0    conda-forge
python                    3.6.3                h1284df2_4

[TomAugspurger]

I would think mkl would be more likely to blame / thank here, rather than bias. Can you try in an environment with nomkl?

[jreback]

pandas doesn't explicity call any BLAS routines at all, or any offer an parallelism, except for dispatch to numexpr, which doesn't apply here. Likely some numpy routines are actually dispatching to multiple threads.

[stuarteberg]

Can you try in an environment with nomkl?

You're right! I tried creating a fresh environment and experimented with the openblas vs. mkl builds of blas. When mkl is active, I see strange parallelism as described above. But when mkl is not active (e.g. viaconda install nomkl), then I see single-threaded behavior.

pandas doesn't explicity call any BLAS routines at all [...] Likely some numpy routines are actually dispatching to multiple threads.

I am not familiar with pandas internals, but naively, I think this is the function of interest:

pandas/pandas/_libs/hashtable_func_helper.pxi.in

Lines 75 to 124 in f6da1f1

	@cython.wraparound(False)
	@cython.boundscheck(False)
	{{if dtype == 'object'}}
	cpdef value_count_{{dtype}}(ndarray[{{dtype}}] values, bint dropna):
	{{else}}
	cpdef value_count_{{dtype}}({{scalar}}[:] values, bint dropna):
	{{endif}}
	cdef:
	Py_ssize_t i=0
	kh_{{ttype}}_t *table
	{{if dtype != 'object'}}
	{{dtype}}_t[:] result_keys
	int64_t[:] result_counts
	{{endif}}
	Py_ssize_t k
	table = kh_init_{{ttype}}()
	{{if dtype == 'object'}}
	build_count_table_{{dtype}}(values, table, 1)
	{{else}}
	build_count_table_{{dtype}}(values, table, dropna)
	{{endif}}
	result_keys = np.empty(table.n_occupied, dtype=np.{{dtype}})
	result_counts = np.zeros(table.n_occupied, dtype=np.int64)
	{{if dtype == 'object'}}
	for k in range(table.n_buckets):
	if kh_exist_{{ttype}}(table, k):
	result_keys[i] = <{{dtype}}> table.keys[k]
	result_counts[i] = table.vals[k]
	i += 1
	{{else}}
	with nogil:
	for k in range(table.n_buckets):
	if kh_exist_{{ttype}}(table, k):
	result_keys[i] = table.keys[k]
	result_counts[i] = table.vals[k]
	i += 1
	{{endif}}
	kh_destroy_{{ttype}}(table)
	{{if dtype == 'object'}}
	return result_keys, result_counts
	{{else}}
	return np.asarray(result_keys), np.asarray(result_counts)
	{{endif}}

I don't see anything suspicious in there, so I'm at a loss as to where the threads are coming from. How is BLAS coming into play here? (Or is mkl used for anything besides BLAS?)

[stuarteberg]

Minor clarification: The pandas/numpy versions I mentioned in the OP are somewhat out-of-date, but this issue is easy to reproduce (on linux) using the most up-to-date packages from both conda-forge and also conda's defaults channel (since mkl is installed either way.)

Simply try conda create -n test pandas and then run the test program above to see the behavior. I guess this means anyone using multiprocessing, pyspark, or dask with a default conda install will see degraded performance if they call value_counts() (and maybe other functions)?.

Should this issue be retagged as a Bug or Performance issue?

[TomAugspurger]

Do you actually observe lower performance when using one of those libraries in an environment with mkl? I may be wrong, but I thought that mkl handled oversubscription.

[stuarteberg]

Do you actually observe lower performance when using one of those libraries in an environment with mkl?

Yes, it appears so.

In my case, I was running 16 processes on a 16-core machine. But since each process was using openmp, they had 16 threads each, so my total thread count was in the 100s. The code I was running happened to be dominated by value_counts(), so when I set OMP_NUM_THREADS=1, the thread count dropped down to normal, and my overall performance jumped up ~6x.

[stuarteberg]

pandas doesn't explicity call any BLAS routines at all [...] Likely some numpy routines are actually dispatching to multiple threads.

I did a little more digging, to see if functions besidesSeries.value_counts() also exhibit multi-core behavior. Here are my results. Maybe a core developer can take a look and discern a pattern.

(In the statements below, s is a large pd.Series with random data.)

---

On my test machine, the following functions use every CPU core available, (consuming 80-100% of cycles per core).

Fully multicore functions

s[:] = s[:]
_ = (s == s)
_ = (s != s)
_ = s.all()
_ = s.any()
_ = s.between(3, 7)
_ = s.copy()
_ = s.corr(s)
_ = s.cov(s)
_ = s.cummax()
_ = s.cummin()
_ = s.cumprod()
_ = s.cumsum()
_ = s.describe()
_ = s.drop_duplicates()
_ = s.dropna()
_ = s.duplicated()
_ = s.idxmax()
_ = s.idxmin()
_ = s.interpolate()
_ = s.kurt()
_ = s.kurtosis()
_ = s.mad()
_ = s.mask(cond)
_ = s.max()
_ = s.mean()
_ = s.median()
_ = s.min()
_ = s.mode()
_ = s.nlargest(10)
_ = s.nsmallest(10)
_ = s.pct_change()
_ = s.ptp()
_ = s.quantile()
_ = s.reindex()
_ = s.replace(1, 5) # Compare with 'inplace' call
_ = s.reset_index(drop=True)
_ = s.sem()
_ = s.skew()
_ = s.sort_index()
_ = s.value_counts()

---

The following functions use all 16 CPU cores, but not at 100%, indicating that they are probably calling one of the above functions, but the majority of their time is spent in a single thread.

Partially multi-core functions

_ = s.align(s2) 
_ = s.eq(s2)
_ = s.ne(s2)
_ = s.rank()
_ = s.sort_values()

---

The following functions use only 3 CPU cores (at 100%), so maybe they are using a different parallelization mechanism:

Moderately parallelized functions

# Exactly 3 cores
_ = s.div(s)
_ = s.divide(s)
_ = s.rdiv(s) # Exactly 3 cores

---

And the following functions use exactly 1 core:

Single-threaded functions

_ = s.abs()
_ = s.add(s)
_ = s.agg(lambda s: s+1)
_ = s.clip(3, 7)
_ = s.clip_lower(3)
_ = s.clip_upper(7)
_ = s.combine(s2, lambda a,b: a+b)
_ = s.combine_first(s2)
_ = s.compound()
_ = s.count()
_ = s.diff()
_ = s.dot(s) # Really???
_ = s.dot(s2) # Really???
_ = s.drop(3)
_ = s.equals(s) # Compare with s.eq(s)
_ = s.ge(s) # compare with s.eq(s)
_ = s.groupby(s % 2)
_ = s.groupby(s % 2).sum()
_ = s.gt(s) # compare with s.eq(s)
_ = s.isin([3,4,5])
_ = s.isna()
_ = s.isnull()
_ = s.le(s)
_ = s.lt(s)
_ = s.map(lambda x: x+2)
_ = s.mod(s)
_ = s.mul(s)
_ = s.multiply(s)
_ = s.nonzero()
_ = s.notna()
_ = s.notnull()
_ = s.nunique()
_ = s.prod()
_ = s.product()
_ = s.radd(s)
_ = s.replace(1, 5, inplace=True) # Compare with non-inplace call
_ = s.sample(1000)
_ = s.tolist()
_ = s.unique()

---

Some interesting observations:

• Apparently there is a big difference between (s == s) (uses all cores) and s.eq(s) (mostly single-threaded).
• There is a big difference in running s.replace(..., inplace=True) vs. s.replace(..., inplace=False), indicating that the multi-CPU usage in the latter case is due to the creation of a new Series for the result.
• Given that we think the multi-core behavior is related to mkl and/or BLAS functions, I expected np.dot to use multiple cores. But to my surprise, it uses only a single core.
• Interestingly, s.any() and s.all() use multiple cores, but the numpy equivalents np.all() and np.any() do not.

---

TEST SETUP:

The above statements were tested in Python, by inserting them into the for-loop shown below. Then I just looked at htop to see how many CPUs were being used.

import numpy as np
import pandas as pd

a = np.random.randint(10, size=20_000_000)
s = pd.Series(a)
s2 = s.iloc[::-1].reindex().copy()

for _ in range(50):
    {{STATEMENT}}

[jreback]

any kind of parallism is likely from numexpr
which we use in expressions and comparisons
eg df + 1 or df == 1

[stuarteberg]

any kind of parallism is likely from numexpr
which we use in expressions and comparisons
eg df + 1 or df == 1

OK, that explains why (s == s) uses all cores.

Even so, it seems strange that s.copy() uses all cores. Does Series.copy() use numexpr?

Also, does the numexpr theory explain the difference between the following two lines?

s.replace(1, 2, inplace=True)  # 1 core
s.replace(1, 2, inplace=False) # 16 cores

(Sorry if that's a dumb question; I don't know much about the pandas internals.)

[jreback]

.copy() is ultimately just np.copy at lowest level
ok multiple arrays (index and values)

all operations call .copy (except inplace)

[stuarteberg]

.copy() is ultimately just np.copy at lowest level

Nice tip, thanks. OK, I verified that, yes, np.copy(a) does indeed use all cores. So that explains many of the above results.

However, some are still mysterious, e.g. Series.any() is multicore. I thought the numexpr theory might explain it, but then I tried export NUMEXPR_NUM_THREADS=1 before running my test. It seems to make no difference -- s.any() still uses all cores, while np.any(a) does not.

In any case, I'm sure there's no way to avoid np.copy() in the pandas internals, even if you wanted to, so the best we can do is put a strong warning in the documentation to encourage users to use OMP_NUM_THREADS=1 before they run their code.

FWIW, there seems to be a related discussion in the numpy repo: numpy/numpy#11826

[stuarteberg]

the best we can do is put a strong warning in the documentation to encourage users to use OMP_NUM_THREADS=1 before they run their code.

I'll write something up. Which section(s) of the docs should I edit? I think the following sections are candidates:

• Install
• Enhancing Performance
• FAQ

[jreback]

a section in FAQ i think would be good