Maybe, this should also depend on n_jobs, e.g. n_threads = min(n_jobs, _openmp_effective_n_threads()).

We decided to restrict n_jobs to control joblib-level parallelism in other estimators (KMeans, HGBDT). So for consistency I would stick to that policy for now. To control the number of OpenMP (and possibly also OpenBLAS threads) users can rely on threadpoolctl which is a dependency of scikit-learn.

However, I suspect that the iterated calls to OpenMP cython code and OpenBLAS code can lead to performance slow downs when OpenBLAS using its own native threading layer instead of the same OpenMP runtime as scikit-learn. This can be confirmed by comparing the performance of this PR in 2 environments:

• one that installs scipy from PyPI that uses an OpenBLAS with its own threadpool (independent of the threadpool used by OpenMP)
• another env where scipy's openmp as explained here: Poor performance of sklearn.cluster.KMeans for numpy >= 1.19.0 #20642 (comment)

Profiling on a MNIST shaped problem shows that this does not seem a problem in practice. Maybe this is because the 2 OpenBLAS calls dominate the Cython calls that are comparatively negligible.

I only tried on macOS though. Maybe the OpenBLAS/OpenMP story is different on Linux.

What is the conclusion? No code change needed?

I only tried on macOS though. Maybe the OpenBLAS/OpenMP story is different on Linux.

@ogrisel: I can experiment with that, do you have your script at hand? If not, no problem.

No code change needed?

I don't think code change is required now to change the way Cython/prange/OpenMP threading is used in this PR but it might be important to keep this potential problem in mind if we observe slowe than expected performance on data with different shapes (e.g. smaller data size maybe where the interacting threadpool problem might still be present or in other locations in the scikit-learn code base whenever with have the following code pattern:

while not converged:
     step 1: Cython prange loop that uses an OpenMP threads
     step 2: numpy matrix multiply / scipy.linalg BLAS open that uses OpenBLAS threads

The long term fix for this kind of problems would ideally happen at the community level in upstream packages (scipy/scipy#15129) but if ever we need to workaround it for a specific in scikit-learn before that happens, we can always decide to use threadpoolctl (or just passing num_threads=1 to prange) to disable threading either for step 1 or step 2 when we detect that OpenBLAS is packaged in a way that it relies on its native threadpool instead of sharing the OpenMP threadpool.

OpenMathLib/OpenBLAS#3187 relates to the issue @ogrisel mentions regarding the following code pattern.

Actually, looking at this again, and independently of the OpenBLAS / OpenMP threadpool interaction hypothetical problem, there might be problematic case when n_jobs > 1 because we could get over-subscription between joblib parallelism and openmp parallelism. Maybe we should always hardcode n_threads=1 to stay on the same side for this PR.

We might want to add a TODO comment to say that we should refactor this all to avoid using joblib parallelism entirely when doing binary and multinomial multiclass classification and only use joblib only for for the one-vs-rest multiclass case. This is probably beyond the scope of this PR.

Done in dea9bf0 (sorry for the stupid commit message).

This is really special for multinomial case / n_classes >= 3. The hessian returned by HalfMultinomialLoss is only the diagonal part, diagonal in classes. Therefore, we call the method gradient_proba which exists only for HalfMultinomialLoss for exactly this reason: to define hessp such that it computes products of a vector with the full hessian for which gradient and predicted probabilties are enough.

This notation implies that all intercept values are contiguously packed at the end of the flat coef array. However in the code we have

            w = coef.reshape(self._loss.n_classes, -1)
            if self.fit_intercept:
                intercept = w[:, -1]
                w = w[:, :-1]

which means that if coef is contiguuous, neither w is not contiguous nor intercept are. Furthermore the intercept values are not packaged at the end of coef anymore

I tried to evaluate the performance impact of using non contiguous coef on a MNIST shaped problem and at least on my machine it does not have a significant impact:

In [1]: import numpy as np

In [2]: n_samples, n_features, n_classes = 60000, 784, 10

In [3]: X = np.random.randn(n_samples, n_features)

In [5]: flat_coef = np.random.normal(size=(n_classes * (n_features + 1)))

In [6]: coef = flat_coef.reshape(n_classes, -1)[:, :-1]

In [7]: intercept = flat_coef.reshape(n_classes, -1)[:, -1]

In [8]: coef.flags
Out[8]: 
  C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

In [9]: intercept.flags
Out[9]: 
  C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

In [10]: coef_contig = np.ascontiguousarray(coef)

In [11]: intercept_contig = np.ascontiguousarray(intercept)

In [12]: %timeit X @ coef.T + intercept
54.8 ms ± 2.1 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [13]: %timeit X @ coef_contig.T + intercept_contig
55.5 ms ± 2.8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Still, I would feel more comfortable we'd shaped the problem to use contiguous arrays everywhere.

You are right! D*m*d.

The scipy solvers only understand 1-d arrays x to solve for, i.e. minimize(func, x, ...). The x in logistic regression is initialized as

w0 = np.zeros(
            (classes.size, n_features + int(fit_intercept)), order="F", dtype=X.dtype
        )

Then this is ravelled for lbfgs and newton-cg:

w0 = w0.ravel()

w0 is now (C- and F-) contiguous. Then in _w_intercept_raw

w = coef.reshape(self._loss.n_classes, -1)

This is surprisingly C-contiguous.

Here is a full example on its own:

import numpy as np

n_classes, n_features = 3,5
w0 = np.zeros((n_classes, n_features), order="F")
w1 = w0.ravel()
w1.flags
#  C_CONTIGUOUS : True
#  F_CONTIGUOUS : True

w2 = w1.reshape(n_classes, -1)
w2.flags
#  C_CONTIGUOUS : True
#  F_CONTIGUOUS : False

w3 = w2[:, :-1]
w3.flags
# C_CONTIGUOUS : False
# F_CONTIGUOUS : False

Possible solution
Should we use np.ravel(..., order="F") and np.reshape(..., order='F')?

w0 = np.zeros(
(classes.size, n_features + int(fit_intercept)), order="F", dtype=X.dtype
)

Oh I missed the order="F" here! So my analysis was flawed, I think.

Maybe we should rename w0 to packed_coef when it is rectangular with the fortran layout as is the case above and then use the view flat_coef = np.ravel(packed_coef) for solvers that need to 1d view representation.

Then we can have a centralized documentation somewhere to explain the memory layout and contiguity assumptions for packed_coef and flat_coef and how to extract the contiguous views for the coef and intercept sub-components throughout this code.

39ffe9d is my attempt as solving this. It is a larger change:

• coef.shape is either (n_dof) or (n_classes, n_dof) or ravelled (n_classes * n_dof,)
  In the multiclass case, this enables to handle both, ravelled and unravelled coefeficients at the same time.
• In all steps involving ravel or reshape, the order is set explicitly, mostly "F".
• The benchmark script above gives 1.821 seconds in this PR vs 2.070 seconds on master.

We decided to restrict n_jobs to control joblib-level parallelism in other estimators (KMeans, HGBDT). So for consistency I would stick to that policy for now. To control the number of OpenMP (and possibly also OpenBLAS threads) users can rely on threadpoolctl which is a dependency of scikit-learn.

However, I suspect that the iterated calls to OpenMP cython code and OpenBLAS code can lead to performance slow downs when OpenBLAS using its own native threading layer instead of the same OpenMP runtime as scikit-learn. This can be confirmed by comparing the performance of this PR in 2 environments:

• one that installs scipy from PyPI that uses an OpenBLAS with its own threadpool (independent of the threadpool used by OpenMP)
• another env where scipy's openmp as explained here: Poor performance of sklearn.cluster.KMeans for numpy >= 1.19.0 #20642 (comment)

Profiling on a MNIST shaped problem shows that this does not seem a problem in practice. Maybe this is because the 2 OpenBLAS calls dominate the Cython calls that are comparatively negligible.

I only tried on macOS though. Maybe the OpenBLAS/OpenMP story is different on Linux.