OrthogonalMatchingPursuit reproducibility with OpenBLAS and different number of cores

[evax]

We've been experiencing reproducibility issues with OrthogonalMatchingPursuit in different computing environments.
The output has been very stable for us when running on hosts with two cores, but we noticed we were getting different results on hosts with eight cores.
The two result sets are stable, but different.

Searching the project's issues it looks like there was some OMP reproducibility issues in the CI environment around 2014, which led to the OMP tests being disabled.

Looking at the code, there isn't a lot of potential culprits explaining such a behavior. Parallel/delayed is used in OrthogonalMatchingPursuitCV, but can be controlled with n_jobs and this doesn't change the outcome.

Diving deeper brings us to scipy's use of OpenBLAS.
We noticed the sklearn test suite is running with OPENBLAS_NUM_THREADS=2 in CI. Setting this environment variable in the 8 cores environment solved our reproducibility issue.

So now to the questions:

• is this a known OpenBLAS issue?
• are the results with OPENBLAS_NUM_THREADS=2 more correct than without with more than two cores or are the results different but equally valid?
• what are the performance implications of running sklearn with OPENBLAS_NUM_THREADS=2?

[ogrisel]

is this a known OpenBLAS issue?

With multithreading enabled (that is OPENBLAS_NUM_THREADS > 1) the ordering of sub-operations (e.g. blocks of computations performed during a matrix matrix multiplication) is not deterministic (it depends on the thread scheduling decisions made by the operating system). From a mathematical point of view, this should not matter because arithmetic operations such as additions performed on the intermediate results are commutative.

However, computers work with finite precision floating point numbers instead of true real numbers. In consequence this can cause non deterministic rounding error changes in the results.

In case you have correlated data, those small changes can dramatically impact the sparsity pattern of the solution OMP converges to.

are the results with OPENBLAS_NUM_THREADS=2 more correct than without with more than two cores or are the results different but equally valid?

None is more correct than the others. The rounding errors that you get with sequential code have no particular reasons to be more correct than the rounding errors you get in non-deterministic orders of execution in general.

The sparsity pattern of the solution of the OMP optimization problem is fundamentally not stable when dealing with correlated features.

To quantify this (lack of) stability, it is a good idea to fit the same models many times with random resampling of the dataset, for instance using ShuffleSplit or RepeatedKFold cross validation and consider the resulting distribution of the coefficients.

You can have a look at this example in the documentation for instance:

https://scikit-learn.org/stable/auto_examples/inspection/plot_linear_model_coefficient_interpretation.html#checking-the-variability-of-the-coefficients

This above example considers Ridge and Lasso regression but you have similar (or even stronger) variability with OMP.

what are the performance implications of running sklearn with OPENBLAS_NUM_THREADS=2?

On multicore machines, OPENBLAS_NUM_THREADS=2 can be up to 2x faster than sequential code OPENBLAS_NUM_THREADS=1.

If you do not set OPENBLAS_NUM_THREADS explicitly, openblas will introspect the number of cores on your machine. If you have 8 cores, it can run up to 8x faster than sequential code. So setting OPENBLAS_NUM_THREADS=2 in this context it can slow down by a factor up to 8 / 2 = 4x compared to not setting OPENBLAS_NUM_THREADS explicitly.

You can introspect the number of threads used by openblas on your machine by default using the threadpoolctl library. For instance here is result on my machine:

% python -m threadpoolctl -i scipy
[
  {
    "filepath": "/Users/ogrisel/miniforge3/envs/dev/lib/libopenblasp-r0.3.12.dylib",
    "prefix": "libopenblas",
    "user_api": "blas",
    "internal_api": "openblas",
    "version": "0.3.12",
    "num_threads": 8,
    "threading_layer": "openmp"
  },
  {
    "filepath": "/Users/ogrisel/miniforge3/envs/dev/lib/libomp.dylib",
    "prefix": "libomp",
    "user_api": "openmp",
    "internal_api": "openmp",
    "version": null,
    "num_threads": 8
  }
]

[evax]

Thank you very much @ogrisel for this very thorough answer, it all makes a lot of sense!

Regarding performance, in your own experience, how much are the OpenBLAS and joblib thread pools detrimental one to the other when they both use the maximum number of cores?
Would it make sense to configure things in such a way that OPENBLAS_NUM_THREADS + n_jobs = number of cores?

[ogrisel]

It really depends on the workload.

Note that if you use the loky backend (default) or the older multiprocessing backend of joblib, there is already a builtin protection against over-subscription in nested OpenBLAS or OpenMP calls: joblib automatically sets the OPENBLAS_NUM_THREADS and OMP_NUM_THREADS variables to max(n_cores // n_jobs, 1) for the spawned joblib workers so that concurrent Python processes do not use too many openblas or openmp threads.

More details here: https://joblib.readthedocs.io/en/latest/parallel.html#avoiding-over-subscription-of-cpu-resources

If you use joblib and openblas (e.g. scipy) sequentially in your main program, then there is no problem. Each of them can use all the available cores, one after the other: there is no oversubcription problem in this case.

But feel free to tweak those parameters on intermediate sized problems (e.g. tasks that run in more than 10s but less than 1min) and measure the effective execution duration while monitoring the CPU usage of your machine to gain better insights.

[evax]

Thanks again @ogrisel, this is very clear and useful!