Talip7/scikit-learn-issues-multilabel huggingface dataset:
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:
Related: apparently `LinearRegression` sometimes yield the minimal norm solution rather in a platform dependent manner as seen in this feedback from the MOOC:

- https://github.com/INRIA/scikit-learn-mooc/issues/641

So this behavior in the `Ridge` class might be either platform dependent or dependent on the version of scipy / LAPACK.

Labels: Bug, module:linear_model
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'm working on it :)

Labels: Bug, module:linear_model
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 think for `Ridge`, it's ok to offer the possibility to use solver that do not converge to the minimum solution when alpha is 0 as long as:

- the docstring for the `solver` parameter clearly documents which solvers are (numerically) safe to use without penalization,
- raise an informative `ConvergenceWarning` when calling Ridge with `alpha < np.finfo(X.dtype).eps` (or similar) and a value of `solver` that does not guarantee to remove the minimum norm solution.

For `LinearRegression` we should probably not expose solvers that are not guaranteed to converge to the minimum norm solution. This specific issue is being discussed in #22855.

Labels: Bug, module:linear_model
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:
Haha, I re-read @lorentzenchr description of the issue:

> Probable Cause
>
> For wide X, the least squares problem reads a bit different: 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 w0, this reads w = X'(XX' + 1 1')^-1 y, where 1 is a column vector of ones. w0 = 1'(XX' + 1 1')^-1 y.
This is incompatible with our mean centering approach.

So this is confirmed by my quick experiments.

Labels: Bug, module:linear_model
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
I'm no expert in sklearn's code base by any means, but as far as I read and have used the code w.r.t. `y_mean` it does the right thing because `WhiteKernel` has [kind of a special API](https://scikit-learn.org/stable/modules/gaussian_process.html#gaussian-process-kernel-api). While for all kernels `k(X) == k(X,X)`, `WhiteKernel` behaves as `w(X) == np.eye(...) * noise_level` while `w(X,X) == np.zeros(...)`. This behavior carries over to sum kernels.

[`fit()` calls](https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L308) `K = self.kernel_(self.X_train_)` while [`predict()` does](https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L392) `K_trans = self.kernel_(X, self.X_train_)`.


Consider this code:

```py
import numpy as np
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

N = 5
lamb = 3.45
X = np.random.rand(N, 3)
X_pred = np.random.rand(2*N, 3)

k = RBF(length_scale=1.23)
w = WhiteKernel(noise_level=lamb)

kw = k + w
I_lamb = np.eye(N) * lamb

assert (k(X) == k(X, X)).all()
assert (w(X) == I_lamb).all()
assert (w(X, X) == np.zeros((N, N))).all()
assert (kw(X, X) == k(X)).all()

# fit()
assert (kw(X) == k(X) + I_lamb).all()

# predict()
assert (kw(X_pred, X) == k(X_pred, X)).all()
```

Labels: module:gaussian_process
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
However [`predict()` also calculates `y_cov`](https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/gaussian_process/_gpr.py#L409). The code does `y_cov = self.kernel_(X) - V.T @ V`. If the above explanation is correct, then one might need to use `self.kernel_(X, X)` instead, which would be `K(X*, X*)` in R&W 2006 (eq. 2.24, 2.26, Alg. 2.1 line 6). The same probably applies to `y_std`.

Labels: module:gaussian_process
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
Nice contrib to solve the pb, hope that sklearn experts will tell us how they proceed and if one should be careful in using the code. Thanks

Labels: module:gaussian_process
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
I don't have the time currently, but it could be worthwhile to use some toy data and compare to `y_cov` from other implementations such as `gpytorch`.

Labels: module:gaussian_process
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
Hello, 
I have setup a Google Collab nb [HERE](https://colab.research.google.com/drive/1zWMxgLRSEO4FsT7dQhH5wCHT7GG6faA7?usp=sharing)
You will see 4 plots (beyond the dataset):
- one showing my personal GP small-and-not-otmized lib. with the implementation of the Algo2.1
- one using the current `Sklean` lib. 
- one using the current `GPy` library using the `predict` method
- and finally one using `GPy predict_noiseless` method

You will certainly conclude that :
- Sklearn and GPy gives the same results using the `predict`method but this is not the Algo2.1
- GPy and my GP lib gives very similar results when using the `GPy predict_noiseless` method which is conform to the Algo 2.1.

So, I think the libraries are correct in their implementations no doubt about that, but clearly `GPy` address to the user the responsibility to use the noisy or the noiseless prediction. My personal suggestion would be that Sklearn may be improved by providing the two options, and mentioning that the  `noiseless` prediction is conform to the Algo 2.1.

Hope that it is ok. I do not know gpytorch but one can add it to the nb if possible

Labels: module:gaussian_process
GaussianProcessRegressor (predict)

### Discussed in https://github.com/scikit-learn/scikit-learn/discussions/22925

<div type='discussions-op-text'>

<sup>Originally posted by **jecampagne** March 23, 2022</sup>
Hello,
I am questioning the code of `predict` of the  `GaussianProcessRegressor`. The code is based on Alg. 2.1 of C. E. Rasmussen & C. K. I. Williams (2003). I have the 2006 version and I do not know if there have been a modification between the two versions. Well, the algorithm is based on the following formula:

![image](https://user-images.githubusercontent.com/20539759/159709115-293c9cc3-aaf9-497f-9c0a-a7c9c10da73f.png)

![image](https://user-images.githubusercontent.com/20539759/159711732-5c3b305a-ecb0-4874-b79b-11680eb670a5.png)


Notice that the `K(X,X)` (ie, X=X_train) is the only one that contains the "noise" parameter while `K(X*, X)` and `K(X*,X*)` (ie, X*=X_test) do not get this additional part to the diagonal. And this is ok.

Looking now at the code it is ok for the default kernel (ie. RBF with fixed scale and length):
```python
self.kernel_ = C(1.0, constant_value_bounds="fixed") * RBF(
                1.0, length_scale_bounds="fixed"
            )
```
But, if the user pass a kernel composed with the `WhiteKernel`as
```python
kernel = C(1.0) * RBF() + WhiteKernel(0.5)
```
then it seems that `K(X*, X)` and `K(X*,X*)` will use the `WhiteKernel` part which is not what the Alg. 2.1 of Rasmussen & Williams does.




</div>

COMMENT:
The point that is rather strange is that if your model does not include the measurement noise as it is fixed (you know that your senor as an error of 0.01 for instance) then you will not include a WhiteNoise kernel and you have a `predict` result. If now you have a little idea on the noise level as it should be around 0.01 but you prefer the model to get a better estimate, than you introduce a WhiteNoise kernel, than let say it ends up to 0.01 with a very little error, so does the `predict`  will give you the same answer? It would be nice to get something clearly state for the user. But ok it may be stated that it is the responsibility of the users to not use library as a `black box`.

Labels: module:gaussian_process