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:
>> The great advantage of lstsq is it's stability w.r.t. badly conditioned X (e.g. singular).

> I think we need tests for this: in particular it would be great that we would call LinearRegression and Ridge on such badly conditioned learning problem and checking that warning with informative error message that suggest to use a more stable solver by its name are raised (or that the switch happens automatically with a warning message to inform the user).

I ran a quick experiment with a well specified, over-complete model with duplicated columns + machine precision noise to probe the behavior of different solvers with near degenerate problems:

```python
import numpy as np
from sklearn.linear_model import LinearRegression, Ridge
from time import perf_counter

rng = np.random.RandomState(42)
n_samples = 10_000_000
n_features = 2
y = rng.randn(n_samples)
X = np.repeat(y.reshape(-1, 1), n_features, axis=1)
X += rng.randn(*X.shape) * np.finfo(X.dtype).eps / 1e2


class PInvLinearRegression(LinearRegression):

    def fit(self, X, y):
        self.coef_ = np.linalg.pinv(X.T @ X) @ X.T @ y
        self.intercept_ = 0
        return self


models = [
    LinearRegression(fit_intercept=False),
    Ridge(fit_intercept=False, alpha=0, solver="lsqr"),
    Ridge(fit_intercept=False, alpha=0, solver="svd"),
    Ridge(fit_intercept=False, alpha=0, solver="cholesky"),
    PInvLinearRegression(),
]


print(f"n_samples: {n_samples}, n_features: {n_features}")
print()


for model in models:
    durations = []
    print(model)
    for _ in range(3):
        t0 = perf_counter()
        model.fit(X, y)
        durations.append(perf_counter() - t0)
    print(f"{model.coef_.max() - model.coef_.min() = }")
    print(f"duration: {np.mean(durations):.3f} ± {np.std(durations):.3f} s")
    print()
```

which outputs for different values of `n_samples` and `n_features`:

```
n_samples: 10000000, n_features: 2

LinearRegression(fit_intercept=False)
model.coef_.max() - model.coef_.min() = 0.9999999165424127
duration: 0.375 ± 0.007 s

Ridge(alpha=0, fit_intercept=False, solver='lsqr')
model.coef_.max() - model.coef_.min() = 0.0
duration: 0.119 ± 0.003 s

Ridge(alpha=0, fit_intercept=False, solver='svd')
model.coef_.max() - model.coef_.min() = 2.979069482632696
duration: 0.163 ± 0.002 s

Ridge(alpha=0, fit_intercept=False, solver='cholesky')
model.coef_.max() - model.coef_.min() = 2.979069482632696
duration: 0.205 ± 0.001 s

PInvLinearRegression()
model.coef_.max() - model.coef_.min() = 5.551115123125783e-17
duration: 0.047 ± 0.000 s
```

```
n_samples: 1000000, n_features: 100

LinearRegression(fit_intercept=False)
model.coef_.max() - model.coef_.min() = 6.938893903907228e-18
duration: 4.316 ± 0.134 s

Ridge(alpha=0, fit_intercept=False, solver='lsqr')
model.coef_.max() - model.coef_.min() = 0.0
duration: 0.257 ± 0.016 s

Ridge(alpha=0, fit_intercept=False, solver='svd')
model.coef_.max() - model.coef_.min() = 20.90600371073267
duration: 7.155 ± 0.320 s

Ridge(alpha=0, fit_intercept=False, solver='cholesky')
model.coef_.max() - model.coef_.min() = 20.90600371073267
duration: 7.013 ± 0.083 s

PInvLinearRegression()
model.coef_.max() - model.coef_.min() = 1.734723475976807e-17
duration: 0.501 ± 0.038 s
```

My conclusions:

- the current solver of `LinearRegression` [`scipy.linalg.lstsq`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html) is indeed one of the slowest for those tall and narrow problems,
- it is often more numerically stable than 'svd' and 'cholesky' but it's not always that stable: depending on the number of features, the random seed, the magnitude of the near-machine level perturbations and the number of samples, it sometimes fails to recover the minimum norm solution. I think this is surprising.
- 'svd' and 'cholesky' are not numericaly stable, but that was expected (at least for 'cholesky'),
- the [`scipy.sparse.linalg.lsqr`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsqr.html) solver of `Ridge` is both very fast and numerically stable: in my experiments it would always recover the minimum norm solution.
- my textbook naive pinv-based solution is even faster and also seems numerically stable.

Of course the pinv-based solution would be catastrophic from a memory usage point of view for large `n_features`. Not sure what are the space complexities of the other solvers.

So I think we should at least expose the `lsqr` in `LinearRegression` (by reusing the `_solver_lsqr` private function from the `Ridge` code base with alpha=0). We should further conduct a more systematic review of the performance profiles (both time and space) and numerical stability and ability to recover the minimum norm solution for different n_samples, n_features and condition numbers of `X.T @ X` or `X @ X.T`.

Then we might want to also add a pinv based solver dedicated to the case where `n_features >> n_samples` but I am not sure how to implement `fit_intercept=True` and sample weights efficiently.

Labels: Performance, module:linear_model