Evaluate multi-threading `LogisticRegression`

When https://github.com/scikit-learn/scikit-learn/pull/32073 is merged, `n_jobs` will have no effect in `LogisticRegression` so should be deprecated.

It's also a good time to consider enabling multi-threading to compute the logistic regression path here https://github.com/scikit-learn/scikit-learn/blob/21e0df780772e9567b09249a05a42dfde6de465d/sklearn/linear_model/_logistic.py#L1373-L1388

So far it was disabled because parallelism was happening at a higher level using joblib. Now we should benchmark the impact of enabling multi-threading to see if it's positive for all the solvers and a wide range of problems.
Like for other estimators using OpenMP-based multi-threading, the number of thread should not be controlled by `n_jobs` but instead use all available core (`n_threads = _openmp_effective_n_threads()`).

COMMENT:
Here are some results of a benchmark comparing `n_threads=1` (denoted as n_thread False in the legends) to `n_trheads=_openmp_effective_n_threads()` (denoted as n_thread True in the legends) for fitting a logistic regression using the code from #32073. I ran the benchmark locally on my machine, limiting the number of available CPUs using either `threadpool_limits` in the script or a `OMP_NUM_THREADS` environment variable. 
I only tested the `lbfgs`, `newton-cg` and `newton-cholesky` solvers as the other ones don't take `n_threads` as an argument.

<details>
<summary>Selection of  plots</summary>

<img width="785" height="455" alt="lbfgs 100 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/7031dca7-5186-4778-baad-7c334d34ba79" />

<img width="776" height="455" alt="newton-cg 100 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/4da22f43-6711-40e0-8a6b-79657128c51f" />

<img width="777" height="455" alt="newton-cholesky 100 samples 100 features 8 classes" src="https://github.com/user-attachments/assets/8ff03e36-604f-4dff-a33c-019dea25cad0" />

<img width="776" height="455" alt="lbfgs 1e4 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/958d8aca-8fb2-40a9-90e8-30b2d7bc4082" />

<img width="768" height="455" alt="newton-cg 1e4 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/9ec3fa86-fd00-4036-be16-4d400b2539f4" />

<img width="776" height="455" alt="newton-cholesky 1e4 samples 100 features 8 classes" src="https://github.com/user-attachments/assets/7e0ed611-4f86-4e02-af77-a4146d0e7865" />

<img width="772" height="455" alt="lbfgs 1e6 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/b028c8b9-265d-4fba-88c6-5b2bb76ae83e" />

<img width="763" height="455" alt="newton-cg 1e6 samples 100 features 32 classes" src="https://github.com/user-attachments/assets/93d61af9-6880-43f0-b1b3-fd323f3b9de3" />

<img width="776" height="455" alt="newton-cholesky 1e6 samples 100 features 8 classes" src="https://github.com/user-attachments/assets/7499001e-288f-4791-8be0-1fa4a4eb9e6a" />

</details>

<details>
<summary>All lbfgs plots</summary>

<img width="576" height="455" alt="Image" src="https://github.com/user-attachments/assets/cfbd2c7c-8a19-40ab-afea-fc61a896115a" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/d0fa9048-0fa9-476b-abfa-bb6881d4a79d" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/033cb9cb-3a0e-4e47-8c35-cd0cbb979e8b" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/bfb3e042-f7ca-4550-b7ec-db606138ac2c" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/72692b6d-5e89-4e79-80c2-605d35f3b7af" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/f0395429-555f-4c16-9645-4acbafbd10e6" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/081376cd-c44e-4b58-bbe6-b370b1bf524a" />

<img width="794" height="455" alt="Image" src="https://github.com/user-attachments/assets/78e0d1c4-39de-4658-810c-69dfda83e66d" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/bf86fc6a-97d0-417a-8a44-971a74db2fd2" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/c892489c-7dc0-47cb-805c-e3070352af75" />

<img width="777" height="455" alt="Image" src="https://github.com/user-attachments/assets/fc2e02fb-b75d-4fc8-936b-56435e15379d" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/54091f2d-5879-4fc2-8d10-8d1ec92606e7" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/81d49c3b-4a2d-459b-a921-d26325e5fdf6" />

<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/ab87175d-a883-4903-a8e4-a551bc471b89" />

<img width="562" height="455" alt="Image" src="https://github.com/user-attachments/assets/ada8f9d9-8e85-44e8-966d-453190d419c2" />

<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/7c98cc97-e8ec-4107-96d8-d93170ba8984" />

<img width="763" height="455" alt="Image" src="https://github.com/user-attachments/assets/ac187b1f-3db4-4c70-bd33-51fc52329a9a" />

<img width="772" height="455" alt="Image" src="https://github.com/user-attachments/assets/4b39e820-80f1-4955-b7d2-3d0f67f7a8b9" />

</details>

<details>
<summary>All newton-cg plots</summary>

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/ae297d27-1284-4791-af25-47c8eb9b4c18" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/aad6949e-1f10-4965-8fae-8851064dea30" />

<img width="576" height="455" alt="Image" src="https://github.com/user-attachments/assets/1b6a844f-8bff-40b4-9fb2-019fd07ca846" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/d032d719-2ed4-4475-92c8-24a5107a7af8" />

<img width="576" height="455" alt="Image" src="https://github.com/user-attachments/assets/60933936-3f7c-4719-a988-a162d09361c0" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/de98cb73-dad2-4622-80c0-34cc448e7c55" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/869ae28c-3d11-4df2-aa61-d163e1d04f7e" />

<img width="794" height="455" alt="Image" src="https://github.com/user-attachments/assets/65f39145-a945-460e-b200-fd32232cadb7" />

<img width="777" height="455" alt="Image" src="https://github.com/user-attachments/assets/e32b5be8-4eed-434f-ab6a-bc669363dfe2" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/c2324824-d105-40ac-88e1-31ea4f181b17" />

<img width="777" height="455" alt="Image" src="https://github.com/user-attachments/assets/80456508-67c9-47dd-87ff-1945501cde98" />


<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/02825b44-2fb5-448c-a408-31b3e7ea76cc" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/a85599d3-31d1-4a27-9c97-4636bc5a6ac4" />

<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/3a07d9e7-9a7b-4815-b33a-486cf35862fb" />

<img width="763" height="455" alt="Image" src="https://github.com/user-attachments/assets/5bd16cef-3be3-45fe-b7c5-d2c6248b2e07" />

<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/bd09c301-8e2c-4d30-a199-9061cca4399f" />

<img width="763" height="455" alt="Image" src="https://github.com/user-attachments/assets/0fb651c0-d782-4288-b120-9df8eb480f0e" />

<img width="763" height="455" alt="Image" src="https://github.com/user-attachments/assets/84e5d42b-5736-4b9f-b3b7-c7814bca0b16" />
</details>

<details>
<summary>All newton-cholesky plots</summary>

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/9665a8b6-7644-4fd9-ad21-3af82dbc911f" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/ba4314d7-df63-4386-adc3-4bf83a7b668c" />

<img width="576" height="455" alt="Image" src="https://github.com/user-attachments/assets/d5955cf2-e32c-4492-8874-77ac437c17b8" />

<img width="777" height="455" alt="Image" src="https://github.com/user-attachments/assets/bbbeda72-2115-4d98-99ca-12c58b02c57e" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/c7d628ae-2707-4ae3-b622-669d941f7bd9" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/b4661b84-6857-4cf0-9737-5736d04ed4b6" />

<img width="785" height="455" alt="Image" src="https://github.com/user-attachments/assets/47b0eee3-287d-4d82-9923-bdfa55756b50" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/2236acb3-86c1-47af-bb24-338c441c0754" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/2939cc1e-6635-4da6-bcf2-6f1eec278963" />

<img width="777" height="455" alt="Image" src="https://github.com/user-attachments/assets/72972bb7-b825-4e08-8c9f-a03ba78a459f" />

<img width="768" height="455" alt="Image" src="https://github.com/user-attachments/assets/8bf491a6-cad2-43de-a187-7ff5f92070da" />

<img width="776" height="455" alt="Image" src="https://github.com/user-attachments/assets/8659a914-2057-41d3-a0e5-5b658d116977" />

</details>

The conclusion I draw from the plots are :

- for a small sample size (100) the overhead of using more CPUs increases the run time proportionally to the number of CPUs, for all solvers
- for a medium sample size (1e4) there is a moderate speedup for `lbfgs` when using more CPUs, with the fastest results with 4 CPUs and then we see some overhead slowdown for higher CPU counts. For `newton-cg` and `newton-cholesky` there seem to be no significant difference.
- for a large sample size (1e6) there is a significant speedup for both `lbfgs` and `newton-cg` with increasing CPU counts, and a saturation around 8 CPUs. For `newton-cholesky`, there are noticeable speedups only for the smallest values for number of featuers and classes, otherwise there is no clear difference.

Additionally, it seems that when there are speedups, they diminish with increasing number of features (10 or 100) and increase with increasing number of classes (2, 8 or 32). 


<details>
<summary>Here is the script I used for the benchmark</summary>

```python
import os
import time
import itertools
import argparse
import warnings
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.exceptions import ConvergenceWarning

from threadpoolctl import threadpool_limits
from rich.progress import (
    Progress,
    SpinnerColumn,
    BarColumn,
    MofNCompleteColumn,
    TextColumn,
    TimeElapsedColumn,
    TimeRemainingColumn,
)

warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", category=FutureWarning)


def main(cpu_limit_type, isolated, n_threads):
    results = []

    n_cpu_list = range(1, 11)
    n_samples_list = (100, 10000, int(1e6))
    n_features_list = (10, 100)
    n_classes_list = (2, 8, 32)
    solver_list = (
        "lbfgs",
        "newton-cg",
        "newton-cholesky",
    )

    param_lists = [n_cpu_list] if cpu_limit_type == "threadpool_limits" else []
    param_lists += [
        n_samples_list,
        n_features_list,
        n_classes_list,
        solver_list,
    ]

    with Progress(
        SpinnerColumn(),
        TextColumn("[progress.description]{task.description}"),
        BarColumn(),
        MofNCompleteColumn(),
        TimeElapsedColumn(),
        TimeRemainingColumn(),
        transient=False,
    ) as progress:
        task = progress.add_task(
            "Benchmarking", total=np.prod([len(l) for l in param_lists])
        )

        for params in itertools.product(*param_lists):
            if cpu_limit_type == "threadpool_limits":
                n_cpu, n_samples, n_features, n_classes, solver = params
            else:
                n_samples, n_features, n_classes, solver = params
                n_cpu = os.environ["OMP_NUM_THREADS"]

            progress.update(
                task,
                description=f"[magenta]CPU:[/] {n_cpu}"
                f" [yellow]Samples:[/] {n_samples} [red]Feat:[/] {n_features}"
                f" [green]Classes:[/] {n_classes} [cyan]Solver:[/] {solver}",
            )

            if n_classes == 32 and solver == "newton-cholesky":
                # takes too long so gets skipped
                progress.update(task, advance=1)
                continue

            n_informative = int(0.6 * n_features)
            n_redundant = int(0.2 * n_features)

            X, y = make_classification(
                n_samples=n_samples,
                n_features=n_features,
                n_informative=n_informative,
                n_redundant=n_redundant,
                n_classes=n_classes,
                random_state=0,
            )
            clf = LogisticRegression(solver=solver, random_state=0)

            start = time.perf_counter()
            if cpu_limit_type == "threadpool_limits":
                with threadpool_limits(limits=n_cpu):
                    clf.fit(X, y)
            else:
                clf.fit(X, y)
            duration = time.perf_counter() - start
            reached_max_iter = int(clf.n_iter_[0]) < clf.max_iter

            results.append(
                {
                    "n_cpu": n_cpu,
                    "n_samples": n_samples,
                    "n_features": n_features,
                    "n_classes": n_classes,
                    "n_informative": n_informative,
                    "n_redundant": n_redundant,
                    "solver": solver,
                    "duration": duration,
                    "cpu_limit_type": cpu_limit_type,
                    "isolated": isolated,
                    "n_threads": n_threads,
                    "reached_max_iter": reached_max_iter,
                }
            )
            progress.update(task, advance=1)

    df = pd.DataFrame.from_records(results)
    df.to_csv(
        "logistic_regression_bench_results.csv",
        mode="a",
        index=False,
        header=not Path("logistic_regression_bench_results.csv").exists(),
    )


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Benchmark the Logistic Regression.")
    parser.add_argument(
        "--isolated",
        action="store_true",
        help="Indicate that this script is the only thing running on the machine.",
    )
    parser.add_argument(
        "--cpu_limit_type",
        choices=("threadpool_limits", "OMP_NUM_THREADS"),
        help="How the CPU limitation is handled.",
    )
    parser.add_argument(
        "--n_threads",
        action="store_true",
        help="Indicate that the version of the logistic regression used is using n_threads > 1.",
    )
    args = parser.parse_args()
    main(args.cpu_limit_type, args.isolated, args.n_threads)

```
</details>



**To summarize** : the parallelization seems to improve the performance at higher sample sizes, especially for `lbfgs` and to a lesser degree also for `newton-cg`. There are slowdowns at low sample sizes due to overhead. The speedups seem greater for lower number of features and higher number of classes.

Labels: API, help wanted, Needs Benchmarks