Add NumPy optimization guide#36
Conversation
Add intel-numpy mkl extension optimizations readme
|
Overall comment: this guide recommends setting the IOMP threading layter for MKL, but pretty much every other PyPI package outside of Intel-distributed NumPy will bundle LibGOMP and could potentially cause incompatibilities. Perhaps it could recommend setting |
|
Comment again that the guide specifically mentions AVX-512 as the highest level of SIMD instructions, but that will become outdated soon as hardware with avx10.2 gets released. |
| | `MKL_DYNAMIC` | `FALSE` | Disable automatic thread scaling | | ||
| | `KMP_AFFINITY` | `granularity=fine,compact,1,0` | Pin threads to physical cores (Intel OpenMP only) | | ||
|
|
||
| `KMP_AFFINITY` is an Intel OpenMP setting, so it applies only when oneMKL is on the Intel runtime (`MKL_THREADING_LAYER=INTEL`); under the GNU layer use `GOMP_CPU_AFFINITY` or `numactl` instead. `KMP_AFFINITY=granularity=fine,compact,1,0` is appropriate for single-socket systems or when running one process per socket. On multi-socket systems without `numactl` it may bind threads across sockets; verify the actual binding with `KMP_AFFINITY=verbose`. |
There was a problem hiding this comment.
What about OMP_PROC_BIND?
There was a problem hiding this comment.
oneMKL has a detailed guide for working with Intel OpenMP under different scenarios: https://www.intel.com/content/www/us/en/docs/onemkl/developer-guide-linux/2026-0/improving-performance-with-threading.html.
KMP_AFFINITY has existed in Intel OpenMP before the OpenMP spec introduced OMP_PROC_BIND and KMP_AFFINITY takes precedence over the OMP_PROC_BIND variables.
Reference: https://www.intel.com/content/www/us/en/docs/dpcpp-cpp-compiler/developer-guide-reference/2025-2/supported-environment-variables.html (check OMP_PROC_BIND)
| | Variable | Recommended value | Effect | | ||
| |---|---|---| | ||
| | `MKL_THREADING_LAYER` | `GNU` (mixed env) or `INTEL` (all-Intel) | Select MKL's OpenMP runtime; see note below | | ||
| | `MKL_NUM_THREADS` | physical core count | Cap MKL thread count | |
There was a problem hiding this comment.
Is this guaranteed to work as intended if you set MKL_NUM_THREADS to number of physical cores, then bind the threads to numbers from the system, but don't specify something like OMP_PLACES=threads? Wouldn't it potentially end up using hyperthreads if the system enumerates them in an interleaved order?
There was a problem hiding this comment.
This behavior depends on MKL_DYNAMIC as described in https://www.intel.com/content/www/us/en/docs/onemkl/developer-guide-linux/2026-0/mkl-dynamic.html.
The default value of MKL_DYNAMIC is TRUE,
regardless of OMP_DYNAMIC, whose default value may be FALSE.
When MKL_DYNAMIC is TRUE, Intel® oneAPI Math Kernel Library (oneMKL)
may use fewer OpenMP threads than the maximum number you specify.
For example, MKL_DYNAMIC set to TRUE enables optimal choice of
the number of threads in the following cases:
- If the requested number of threads exceeds the number of physical cores
(perhaps because of using the Intel® Hyper-Threading Technology),
Intel® oneAPI Math Kernel Library (oneMKL) scales down the number of
OpenMP threads to the number of physical cores.
This means if MKL_DYNAMIC=FALSE, user-side thread-binding is respected by oneMKL and hyperthreading can take place.
There was a problem hiding this comment.
Good catch,
@vmalia thanks for pitching in,
if you set MKL_DYNAMIC=FALSE (which this guide's table currently recommends), that automatic scale-down is off and oneMKL respects
| print(f"speedup : {stock_ms / mkl_ms:.1f}x") | ||
| ``` | ||
|
|
||
| Measured on AWS, Intel® Xeon® 6975P-C, 16 cores / 32 threads (HT on), 1 socket, Ubuntu 26.04 LTS. Numbers vary by hardware. |
There was a problem hiding this comment.
Have these numbers been through PDT?
There was a problem hiding this comment.
These are example snippets and freshly measured (FFT on the AWS Xeon 6975P-C) but have not been through PDT.
where as the benchmarking results section displayed numbers and speedups are been through the PCR review process and approved and published on Intel Distribution for Python page - https://www.intel.com/content/www/us/en/developer/tools/oneapi/distribution-for-python.html
|
|
||
| --- | ||
|
|
||
| ## Benchmark results |
There was a problem hiding this comment.
Have these numbers all been through PDT?
| conda activate idp_env | ||
| ``` | ||
|
|
||
| Pin python version to match your project if you need a specific interpreter. NumPy comes from conda-forge; the Intel channel supplies the `mkl_fft`/`mkl_random`/`mkl_umath` extensions and Intel's latest oneMKL builds. To add oneMKL to an *existing* environment that already has conda-forge NumPy installed, swap its BLAS to the MKL variant and add the extensions in place (this re-links the NumPy you already have, it does not reinstall NumPy). Intel-channel packages are built to be compatible with conda-forge but not with the Anaconda defaults channel: |
There was a problem hiding this comment.
Why does it need to pull them from the intel repository? They are also available in conda-forge:
https://anaconda.org/channels/conda-forge/packages/mkl_random/overview
| Use `--index-url`, not `--extra-index-url`: Intel's index is a partial mirror, and with `--extra-index-url` pip would see PyPI's higher-numbered OpenBLAS wheel and install that instead. Packages Intel does not mirror (for example `threadpoolctl`, used for [verification](#verifying-onemkl-is-active)) install normally from PyPI in a separate step. The Intel wheels target Linux and Windows; if `pip` reports no matching distribution, check that your platform and Python version are covered on the index. | ||
|
|
||
| Whichever path you take, choose the OpenMP threading layer and set it **before anything imports NumPy or MKL**. The variable is read once at MKL load time, so exporting it after the import has no effect. Which value to pick is explained under [Threads and NUMA](#threads-and-numa); the safe default for a typical pip or mixed environment is: | ||
| Whichever path you take, choose the OpenMP threading layer and set it **before anything imports NumPy, SciPy or MKL**. The variable is read once at MKL load time, so exporting it after the import has no effect. Intel OpenMP gives the best oneMKL performance on Intel hardware and is the recommended starting point; the environment-specific cases are covered under [Threads and NUMA](#threads-and-numa): |
There was a problem hiding this comment.
This part is not very descriptive. It could be more explicit about all-Intel vs. mixture of packages.
|
|
||
| ```bash | ||
| export MKL_THREADING_LAYER=GNU # share one OpenMP runtime (libgomp) with other packages | ||
| export MKL_THREADING_LAYER=INTEL # Intel OpenMP (libiomp5): best oneMKL performance |
There was a problem hiding this comment.
I think this is the default and setting it like this doesn't change anything. CC @vmalia .
| ### Linear Algebra: BLAS and LAPACK | ||
|
|
||
| This is the lever you get for free. Once oneMKL is the backend, `np.dot`, `np.matmul`, `np.linalg.*`, and everything built on them (covariances, distances, decompositions) run on oneMKL's BLAS and LAPACK with no code change and nothing to activate. These kernels dispatch at runtime to an optimized code path for the CPU's instruction set (e.g., AVX-512 on current Xeons). This is the largest single contributor to the geomean in [Benchmark results](#benchmark-results). | ||
| This is the lever you get for free. Once oneMKL is the backend, `np.dot`, `np.matmul`, `np.linalg.*`, and everything built on them (covariances, distances, decompositions) run on oneMKL's BLAS and LAPACK with no code change and nothing to activate. These kernels dispatch at runtime to an optimized code path for the CPU's instruction set (e.g., AVX-512 on current Xeons). This is the largest single contributor to the geomean speedup in [Benchmark results](#benchmark-results). |
There was a problem hiding this comment.
Has redundant comma after 'e.g.'.
| > **Reproducibility note:** `mkl_random` and `numpy.random` produce different sequences from the same seed. If your tests or simulations depend on specific random values, do not swap them. | ||
|
|
||
| `brng='MT19937'` selects the Mersenne Twister generator, matching `numpy.random`'s default algorithm. `method='BoxMuller'` is the faster of oneMKL's two normal-sampling methods; the alternative `'ICDF'` (inverse CDF) is slower but more accurate in the tails. | ||
| `brng='MT19937'` selects the Mersenne Twister generator, matching `numpy.random`'s default algorithm. For normal sampling, [`'BoxMuller'`](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform) is the faster method and the reason it is used here; the alternative [`'ICDF'`](https://en.wikipedia.org/wiki/Inverse_transform_sampling) (inverse-CDF sampling) is slower but maps each uniform to a normal one-to-one, which some workloads prefer for tail fidelity. For raw generation speed, keep `'BoxMuller'`. |
There was a problem hiding this comment.
NumPy's and SciPy's default is Ziggurat. I don't think that can be described as ICDF, and it involves a variable amount of random bit draws per output normal number.
| Alternatively, to partition a single generator instead of using the `MT2203` family, `skipahead(n)` advances a stream by `n` steps so each worker starts at a non-overlapping offset: | ||
|
|
||
| ```python | ||
| block = 1_000_000 |
There was a problem hiding this comment.
I think something like 2^63-1 (INT64_MAX) would be a more reasonable number. Can also be obtained programmatically from numpy np.iinfo(np.int64).
| `KMP_AFFINITY` is an Intel OpenMP setting, so it applies only when oneMKL is on the Intel runtime (`MKL_THREADING_LAYER=INTEL`); under the GNU layer use `GOMP_CPU_AFFINITY` or `numactl` instead. `KMP_AFFINITY=granularity=fine,compact,1,0` is appropriate for single-socket systems or when running one process per socket. On multi-socket systems without `numactl` it may bind threads across sockets; verify the actual binding with `KMP_AFFINITY=verbose`. | ||
| `KMP_AFFINITY` is an Intel OpenMP setting, so it applies only when oneMKL is on the Intel runtime (`MKL_THREADING_LAYER=INTEL`). For a runtime-agnostic alternative that also works under LLVM and GNU, use the standard OpenMP controls `OMP_PROC_BIND=close` and `OMP_PLACES=cores`; `OMP_PLACES=cores` is what keeps threads off hyperthread siblings, which core-count alone does not guarantee. In `KMP_AFFINITY=granularity=fine,compact,1,0`, `granularity=fine` pins each thread to a single logical CPU, `compact` places threads on adjacent cores, and the trailing `1,0` are the permute and offset ([Intel reference](https://www.intel.com/content/www/us/en/docs/dpcpp-cpp-compiler/developer-guide-reference/2025-0/thread-affinity-interface.html)). This is appropriate for single-socket systems or when running one process per socket; on multi-socket systems without `numactl` it may bind threads across sockets. To see the actual binding at startup, prepend `verbose`: `KMP_AFFINITY=verbose,granularity=fine,compact,1,0`. | ||
|
|
||
| `MKL_DYNAMIC` [defaults to `TRUE`](https://www.intel.com/content/www/us/en/docs/onemkl/developer-guide-linux/2026-0/mkl-dynamic.html), and while it is on, oneMKL scales the thread count down to the number of physical cores when a higher count is requested, which keeps work off hyperthread siblings automatically. Setting `MKL_DYNAMIC=FALSE` fixes the thread count and honors your binding as-is, so pair it with `OMP_PLACES=cores` (or `numactl`) to ensure threads land on distinct physical cores rather than hyperthreads. |
There was a problem hiding this comment.
Important to mention that it needs to be paired with MKL_NUM_THREADS, otherwise it will still end up using all threads.
| The one situation that hurts performance is loading two OpenMP runtimes in the same process, for example oneMKL on Intel OpenMP (`libiomp5`) alongside a package that bundles GNU's `libgomp`. Two runtimes over-subscribe the cores. `MKL_THREADING_LAYER` controls oneMKL's runtime only; it does not stop another package from loading its own. When your environment already standardizes on one OpenMP runtime, match oneMKL to it: | ||
|
|
||
| In conda, the runtime is selected through the `_openmp_mutex` metapackage. To pin LLVM's LibOMP (the upstreamed continuation of Intel OpenMP) on Linux: | ||
| - With LLVM's LibOMP present, use `MKL_THREADING_LAYER=INTEL`: oneMKL runs its Intel threading layer on LibOMP and the `KMP_*` controls apply. |
There was a problem hiding this comment.
Shouldn't it also set MKL_THREADING_LAYER=GNU when coupled with LLVM-OMP libraries? CC @vmalia
| - With LLVM's LibOMP present, use `MKL_THREADING_LAYER=INTEL`: oneMKL runs its Intel threading layer on LibOMP and the `KMP_*` controls apply. | ||
| - If other packages in the process pull in `libgomp` (some pip-installed scientific packages do), `MKL_THREADING_LAYER=GNU` selects oneMKL's GNU threading layer so the process shares a single runtime. | ||
|
|
||
| Platform defaults differ. Linux pip stacks often bundle `libgomp`; conda environments select the runtime through the `_openmp_mutex` metapackage (default LLVM LibOMP on recent conda-forge). On Windows conda, NumPy/SciPy/MKL resolve `llvm-openmp`, and oneMKL there reports the `intel` threading layer. To pin LLVM's LibOMP on Linux: |
There was a problem hiding this comment.
Default on linux is GOMP, unless MKL from conda-forge is installed, which switches it to LLVM-OMP. Not sure if the packages from the intel channel also set the OMP mutex though.
|
@jharlow-intel @vchamarthi If following the advise from this guide with intel packages, it ends up pulling GOMP instead of LLVM-OMP. Is that intended behavior? Are the intel packages meant to set the OMP mutex to LibOMP like the conda-forge ones do? See log below Details |
Adds a new tuning guide documenting how to run NumPy with Intel® oneMKL-backed performance (BLAS/LAPACK plus optional FFT/random/umath patching), and links it from the repository’s main README
Changes:
CC @xaleryb @jharlow-intel @napetrov for addition review