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:
I just played for a course that I will do tomorrow and indeed our `LinearRegression` is terrible.

Given the following dataset:

```python
from sklearn.datasets import make_regression

X, y = make_regression(
    n_samples=100_000,
    n_features=1_000,
    noise=10,
)
```

I did a quick comparison:

```python
%%time
from sklearn.linear_model import LinearRegression

linear_regression = LinearRegression()
linear_regression.fit(X, y)
```

```shell
CPU times: user 1min 50s, sys: 3.84 s, total: 1min 54s
Wall time: 18.5 s
```

and then I though about using the Normal equation:

```python
%%time
coef = np.linalg.inv(X_with_dummy.T @ X_with_dummy) @ X_with_dummy.T @ y
```

```shell
CPU times: user 33.1 s, sys: 850 ms, total: 33.9 s
Wall time: 6.36 s
```

and finally using a L-BFGS-B solver:

```python
def loss(coef, X, y):
    y_estimated = X @ coef
    error = y_estimated - y
    cost = (1 / 2 * len(y)) * np.sum(error ** 2)
    return cost


def grad(coef, X, y):
    y_estimated = X @ coef
    error = y_estimated - y
    gradient = (1 / len(y)) * X.T @ error
    return gradient

coef = np.random.randn(X_with_dummy.shape[1])
res = minimize(loss, coef, args=(X_with_dummy, y), method="L-BFGS-B", jac=grad)
coef = res.x
```

```shell
CPU times: user 6.84 s, sys: 3.73 s, total: 10.6 s
Wall time: 1.89 s

fun: 495382014981.9635
 hess_inv: <1001x1001 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 7.67078546e-07, -6.21855026e-07,  7.71316058e-07, ...,
        1.05001006e-06, -2.00795501e-06,  2.40024978e-08])
  message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 10
      nit: 7
     njev: 10
   status: 0
  success: True
        x: array([ 0.08592853, -0.03905589,  0.01284664, ...,  0.03827157,
        0.08447174, -0.03549184])
```

I recall a discussion in the book of @ageron regarding the complexity of the `lstsq` and the gradient descent: https://www.oreilly.com/library/view/hands-on-machine-learning/9781491962282/ch04.html

It might be nice to investigate as well the impact of the LAPACK driver because I am not sure that they are all equivalent in different `n_samples/n_features` regimes.

NB: I hope that I did not mess up the gradient computation, it was done while catching a train so better be prepared for a stupid error.

NB2: I missed to read the above discussion that already discuss the different solvers inclusion.

Labels: Performance, module:linear_model