The fit performance of LinearRegression is sub-optimal

It seems that the performance of Linear Regression is sub-optimal when the number of samples is very large. 

[sklearn_benchmarks](https://github.com/mbatoul/sklearn_benchmarks) measures a [speedup of 48](https://mbatoul.github.io/sklearn_benchmarks/results/github_ci/master/20220315T132336/scikit_learn_intelex_vs_scikit_learn.html#linearregression) compared to an optimized implementation from scikit-learn-intelex on a `1000000x100` dataset. For a given set of parameters and a given dataset, we compute the speed-up `time scikit-learn` / `time sklearnex`. A speed-up of 48 means that sklearnex is 48 times faster than scikit-learn on the given dataset.

[Profiling](https://mbatoul.github.io/sklearn_benchmarks/results/github_ci/improve_cli/20220314T105014/profiling/sklearn_fit_0c037d7431f8fb8e89db7b3f3f29b083_ddd328af83cd4ef15164ddaa0fa01ba4.json.gz) allows a more detailed analysis of the execution of the algorithm. We observe that most of the execution time is spent in the [`lstsq` solver of scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html).

<img width="1538" alt="image" src="https://user-images.githubusercontent.com/36322423/158272733-92f107f4-92cf-4e26-887d-2a0818113d2d.png">

The profiling reports of sklearn_benchmarks can be viewed with [Perfetto UI](https://ui.perfetto.dev/).

<details>
  <summary>See benchmark environment information</summary>
  
<img width="946" alt="image" src="https://user-images.githubusercontent.com/36322423/158276716-2ab08b7d-079d-4cb8-82fc-054e6b9e8f40.png">
</details>

It seems that the solver could be better chosen when the number of samples is very large. Perhaps Ridge's solver with a zero penalty could be chosen in this case. On the same dimensions, [it shows better performance](https://mbatoul.github.io/sklearn_benchmarks/results/github_ci/improve_cli/20220314T105014/scikit_learn_intelex_vs_scikit_learn.html#Ridge).

Speedups can be reproduced with the following code:

```bash
conda create -n lr_perf -c conda-forge scikit-learn scikit-learn-intelex numpy jupyter
conda activate lr_perf
```

```python
from sklearn.linear_model import LinearRegression as LinearRegressionSklearn
from sklearnex.linear_model import LinearRegression as LinearRegressionSklearnex
from sklearn.datasets import make_regression
import time
import numpy as np

X, y = make_regression(n_samples=1_000_000, n_features=100, n_informative=10)

def measure(estimator, X, y, n_executions=10):
    times = []
    while len(times) < n_executions:
        t0 = time.perf_counter()
        estimator.fit(X, y)
        t1 = time.perf_counter()
        times.append(t1 - t0)

    return np.mean(times)

mean_time_sklearn = measure(
    estimator=LinearRegressionSklearn(),
    X=X,
    y=y
)

mean_time_sklearnex = measure(
    estimator=LinearRegressionSklearnex(),
    X=X,
    y=y
)

speedup = mean_time_sklearn / mean_time_sklearnex
speedup
```

COMMENT:
Thanks for the feedback. I guess we need to do some more systematic benchmarking (including with various combinations for dense vs sparse, with or without sample weight, fit intercept or not, on well specified / mispecified, well-conditioned / ill conditioned).

But in my experience (with `make_regression` with various `n_samples` / `n_features`, with or without biased data, ...), I found that `lsqr` and `lsmr` with `tol=1e-12` are almost always faster than `cholesky` (sometimes both have similar speed and I never saw a case with `cholesky` being more than 2x faster). The `mean_squared_error` are the same the 12 significant digits.

Since `lsqr` and `lsmr` are also numerically stable (and recover the minimum norm solution in case of singular `X.T @ X`) I think either of them would be fool-proof default, at least for `LinearRegression`. `lsmr` is sometimes a tiny bit faster than `lsqr` but not always.

I also tried `lsqr` and `lsmr` with `tol=1e-5` and they are even faster (typically around 2-3x) but then this is another hyper-parameter to adjust and we probably don't want to expose that for `LinearRegression`.

It seems that for biased problems, and with `fit_intercept=True`, the `y_mean` centering in `_preprocess_data` can make the convergence a tiny bit faster for `lsqr`. So we might want to keep it, but I think we would then converge to the minium `np.linalg.norm(coef_)` solution instead of the minimum `np.linalg.norm(np.concatenate(coef_, [intercept_]))` solution. I think either are fine as long as we document and properly test those behavior to avoid surprises in the future.

One think that I have not evaluated yet is the memory usage induced by the different solvers.

Labels: Performance, module:linear_model