BUG unpenalized Ridge does not give minimum norm solution

#### Describe the bug
As noted in #22910, `Ridge(alpha=0, fit_intercept=True)` does not give the minimal norm solution for wide data, i.e. `n_features > n_samples`.

Note that we nowhere guarantee that we provide the **minimum norm solution**.

Edit: Same seems to hold for `LinearRegression`, see #26164.

#### Probable Cause
For wide X, the least squares problem reads a bit different: $\mathrm{min} ||w||_2$ subject to $Xw = y$ with solution $w = X'(XX')^{-1} y$, see e.g. http://ee263.stanford.edu/lectures/min-norm.pdf.
With explicit intercept $w_0$, this reads $w = X'(XX' + 1 1')^{-1} y$, where 1 is a column vector of ones. $w_0 = 1'(XX' + 1 1')^{-1} y$.
**This is incompatible with our mean centering approach.**

#### Example

<details>

```python
import numpy as np
from numpy.testing import assert_allclose
from scipy import linalg
from sklearn.datasets import make_low_rank_matrix
from sklearn.linear_model import Ridge

n_samples, n_features = 4, 12  # wide data
k = min(n_samples, n_features)
rng = np.random.RandomState(42)
X = make_low_rank_matrix(
    n_samples=n_samples, n_features=n_features, effective_rank=k
)
X[:, -1] = 1  # last columns acts as intercept
U, s, Vt = linalg.svd(X)
assert np.all(s) > 1e-3  # to be sure X is not singular
U1, U2 = U[:, :k], U[:, k:]
Vt1, _ = Vt[:k, :], Vt[k:, :]
y = rng.uniform(low=-10, high=10, size=n_samples)
# w = X'(XX')^-1 y = V s^-1 U' y
coef_ols = Vt1.T @ np.diag(1 / s) @ U1.T @ y

model = Ridge(alpha=0, fit_intercept=True)
X = X[:, :-1]  # remove intercept
intercept = coef_ols[-1]
coef = coef_ols[:-1]
model.fit(X, y)

# Check that we have found a solution => residuals = 0
assert_allclose(model.predict(X), y)
# Check that `coef`, `intercept` also provide a valid solution
assert_allclose(X @ coef + intercept, y)

# Ridge does not give the minimum norm solution. (This should be equal.)
np.linalg.norm(np.r_[model.intercept_, model.coef_]) > np.linalg.norm(
    np.r_[intercept, coef]
)
```
This last statement should be be `False`. It proves that Ridge does not give the miminum norm solution. 

</details>

COMMENT:
> I suspect that to recover the minimum norm solution for iterative solvers such as "lbfgs" or "lsqr", we need to do a dump init of the all the parameters (including the intercept) at zero.

Actually, for `lsqr`, this is not really an init problem but more of the way we pose the problem. We preprocess the data based on `fit_intercept` to mean center `y` and then we use that assumption to simplify the handling of the sample weight:

https://github.com/scikit-learn/scikit-learn/blob/30bf6f39a7126a351db8971d24aa865fa5605569/sklearn/linear_model/_base.py#L276-L277

https://github.com/scikit-learn/scikit-learn/blob/30bf6f39a7126a351db8971d24aa865fa5605569/sklearn/linear_model/_ridge.py#L148-L175

If instead of mean centering `y`, we pass it raw to `lsqr`, we actually recover the minimum norm solution. Here is the proof of concept estimator I used to check this:

```python
class LsqrLinearRegression(LinearRegression):

    def __init__(self, fit_intercept=True, tol=1e-5):
        self.fit_intercept = fit_intercept
        self.tol = tol

    def fit(self, X, y):
        if self.fit_intercept:
            X = np.hstack([X, np.ones((X.shape[0], 1))])
        result = sparse_linalg.lsqr(
            X, y, damp=0, atol=self.tol, btol=self.tol, iter_lim=None
        )
        fitted_params = result[0]
        self.n_iter_ = result[2]
        if self.fit_intercept:
            self.coef_ = fitted_params[:-1]
            self.intercept_ = fitted_params[-1]
        else:
            self.coef_ = fitted_params
            self.intercept_ = 0
        return self
```

So I think we should try to avoid doing the `_preprocess_data` trick and instead rely on some linear operators when needed to do the `sample_weight` handling on the fly, both for sparse and dense data.

Still we would need to also empirically study the speed of convergence, both in terms of wall clock time and `n_iter_` and the impacts in terms of test R^2 do decide about which approach is the best. Maybe the inductive bias induced by the current `_preprocess_data` is not that bad.

cc @agramfort @lorentzenchr.

Labels: Bug, module:linear_model