LogisticRegression coefficients should be stored in column major order to improve sparse inference performance.

### Describe the workflow you want to enable

I want to enable the workflow of training a `LogisticRegression` model on sparse data, with substantial numbers of features and classes, with acceptable performance. Currently, inference for a `LogisticRegression` model slows rather substantially when using sparse features since the transpose of the coefficient matrix must be stored in row-major order to perform the sparse matrix multiply. Changing the ordering of the coefficient array improves the performance dramatically. This is illustrated in the below code snippet:

```python
import scipy
import numpy as np

from sklearn.linear_model import LogisticRegression

X = np.random.rand(10_000, 10_000)
y = np.random.randint(0, 500, size=10_000)
clf = LogisticRegression(max_iter=10).fit(X, y)

%timeit clf.predict_proba(scipy.sparse.rand(1, 10_000))
# 21.9 ms ± 973 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%prun clf.predict_proba(scipy.sparse.rand(1, 10_000))
#   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
#        2    0.019    0.010    0.019    0.010 {method 'ravel' of 'numpy.ndarray' objects}
#        1    0.003    0.003    0.022    0.022 _compressed.py:493(_mul_multivector)

clf.coef_ = np.asfortranarray(clf.coef_)

%timeit clf.predict_proba(scipy.sparse.rand(1, 10_000))
# 467 µs ± 11 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%prun clf.predict_proba(scipy.sparse.rand(1, 10_000))
#   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
#        1    0.000    0.000    0.000    0.000 {built-in method scipy.sparse._sparsetools.csr_matvecs}
#        1    0.000    0.000    0.000    0.000 {method 'choice' of 'numpy.random.mtrand.RandomState' objects}
```

Additional information can be found [here](https://stackoverflow.com/questions/74325646/why-doesnt-scikit-learns-logisticregression-classifier-use-column-major-for-co).

### Describe your proposed solution

My proposed solution is to ensure that the coefficient matrix in `LogisticRegression` is converted to column major order with `np.asfortranarray` once training has completed. For sparse input, this will dramatically improve performance, while for dense inference, it will not degrade performance too much, as demonstrated below:

```python
import scipy
import numpy as np

from sklearn.linear_model import LogisticRegression

X = np.random.rand(10_000, 10_000)
y = np.random.randint(0, 500, size=10_000)
clf = LogisticRegression(max_iter=10).fit(X, y)

X_test = scipy.sparse.rand(1, 10_000)
X_dense = X_test.A

%timeit clf.predict_proba(X_test)
# 23.4 ms ± 253 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit clf.predict_proba(X_dense)
# 2.81 ms ± 277 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

clf.coef_ = np.asfortranarray(clf.coef_)

%timeit clf.predict_proba(X_test)
# 201 µs ± 705 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit clf.predict_proba(X_dense)
# 3.05 ms ± 48.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```

### Describe alternatives you've considered, if relevant

Alternatives include simply ensuring the input to inference is always dense, which results in a 10x speedup. However, this is still a ~10x slowdown relative to the column major ordering solution.

### Additional context

_No response_

COMMENT:
I made a benchmark to evaluate in which regime, we actually have a problem:

<details>

```python
# %%
import collections
import itertools
import time
import warnings

import numpy as np
import pandas as pd

from sklearn.exceptions import ConvergenceWarning
from sklearn.linear_model import LogisticRegression
from sklearn.utils._testing import _convert_container

rng = np.random.default_rng(0)

n_repeats = 10
n_classes = [2, 10, 100]
n_samples_test = [1, 10, 100, 1_000, 10_000]
n_features = [10, 100, 1_000, 10_000]
make_sparse = True
as_fortran = [True, False]

results = collections.defaultdict(list)
for ns, nf, nc, is_f in itertools.product(n_samples_test, n_features, n_classes, as_fortran):
    local_results = []
    for _ in range(n_repeats):
        print(
            f"n_samples: {ns}, n_features: {nf}, repeat: {_ + 1}/{n_repeats}, "
            f"as_fortran: {is_f}, sparse: {make_sparse}, cl_asses: {nc}"
        )
        X_train = rng.normal(size=(1_000, nf))
        X_test = rng.normal(size=(ns, nf))
        y_train = rng.integers(0, nc, size=(1_000,))

        if make_sparse:
            # X_train = _convert_container(X_train, "sparse_csr")
            X_test = _convert_container(X_test, "sparse_csr")

        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=ConvergenceWarning)
            model = LogisticRegression(max_iter=1).fit(X_train, y_train)
        if is_f:
            model.coef_ = np.asfortranarray(model.coef_)
        start = time.time()
        model.predict_proba(X_test)
        end = time.time()

        local_results.append((end - start) * 1_000)  # in ms
    results["n_samples"].append(ns)
    results["n_features"].append(nf)
    results["n_classes"].append(nc)
    results["format"].append("Fortran" if is_f else "C")
    results["mean_time"].append(np.mean(local_results))
    results["std_time"].append(np.std(local_results))
results = pd.DataFrame(results)

# %%
results.to_csv("results.csv", index=False)
```
</details>

The results could be condensed into 2 figures:

![image](https://user-images.githubusercontent.com/7454015/201664775-04615f90-c98e-43bf-b987-71cbf83943db.png)

![image](https://user-images.githubusercontent.com/7454015/201664872-f695d0d3-5bd1-4130-b805-32417d98d377.png)

So from the figure above, the problem seems to happen only when `n_classes` is large (>100) and `n_samples` is really low (single data point). We also see that there are numerous cases where the current implementation is better or equivalent.

In practice, I am not sure that we will encounter a case where we have `n_classes` to be that large. Also, I don't think that we want to optimize the case for single-sample prediction. For this case, I am thinking that we have too much overhead on the data validation that we should not encourage it. Instead, serialization with ONNX should be preferred.

@ogrisel WDYT?

Labels: New Feature, module:linear_model