Repository: scikit-learn/scikit-learn
Pull request: FIX scoring != None for RidgeCV should used unscaled y for evaluation
closes https://github.com/scikit-learn/scikit-learn/issues/13998
closes https://github.com/scikit-learn/scikit-learn/pull/15648

While discussing with @jeromedockes, we recall to have observed something weird in the `RidgeCV` code. I check a bit closer and I open this PR to highlight what is the potential problem.

In `RidgeCV`, when having `sample_weight` we scale the data using the `sqrt(sample_weight)`:

https://github.com/scikit-learn/scikit-learn/blob/35164b3e8b605bc990eaf42a8a170082036dff59/sklearn/linear_model/_ridge.py#L2133-L2136

The idea is that the mean squared error can be expressed as:

https://github.com/scikit-learn/scikit-learn/blob/35164b3e8b605bc990eaf42a8a170082036dff59/sklearn/linear_model/_base.py#L212-L223

Those "centered" data are used to optimize the ridge loss. Later in the code, we want to compute a score that can be an arbitrary metric via a scorer.

https://github.com/scikit-learn/scikit-learn/blob/35164b3e8b605bc990eaf42a8a170082036dff59/sklearn/linear_model/_ridge.py#L2158-L2169

The problem here is that `predictions` is computed efficiently as provided in the GCV paper. But these predictions are in the "scaled" space and it seems incorrect to compute any metric in this space with an arbitrary metric. Instead, we should unscale these predictions and the scaled true targets to compute the metric in the original space.

This is what this PR is intended to. I did not add any non-regression test (I assume that using the MedAE should lead to some failures) because I wanted to be sure that what I'm saying is correct.

@jeromedockes @ogrisel @lorentzenchr Does the above description make sense to you?

Edit: It seems that it relates to https://github.com/scikit-learn/scikit-learn/issues/13998 and https://github.com/scikit-learn/scikit-learn/pull/15648

Probably, I should check the tests that were written in https://github.com/scikit-learn/scikit-learn/pull/15648

Problem statement:
`RidgeCV` `cv_values_` are for preprocessed data: centered and scaled by sample weights.
when `store_cv_values=True`, `RidgeCV` stores the leave-one-out squared errors,
when `scoring=None`, or the leave-one-out predictions, when `scoring` is
provided by the user, in its `cv_values_` attribute.

However, when `scoring` is not `None`, it stores the predictions for the
preprocessed data, i.e. rescaled by the square roots of the sample weights and
with the mean of `y` removed:

```python
import numpy as np
from sklearn.linear_model import RidgeCV
from sklearn.datasets import make_regression

x, y = make_regression(n_samples=6, n_features=2, random_state=0)
squared_error = RidgeCV(
    store_cv_values=True, alphas=[10.]).fit(x, y).cv_values_.ravel()
custom_scoring = RidgeCV(
    store_cv_values=True, scoring='neg_mean_squared_error',
    alphas=[10.]).fit(x, y)
# to get the actual predictions we need to add the y mean
custom = (y - (custom_scoring.cv_values_.ravel() + y.mean()))**2
assert np.allclose(squared_error, custom)

sw = np.arange(6) + 1
squared_error = RidgeCV(store_cv_values=True, alphas=[10.]).fit(
    x, y, sample_weight=sw).cv_values_.ravel()
custom_scoring = RidgeCV(
    store_cv_values=True, scoring='neg_mean_squared_error',
    alphas=[10.]).fit(x, y, sample_weight=sw)
# to get the actual predictions we need to rescale by inverse square root
# sample weights and add the y mean
custom = sw * (y
               - (custom_scoring.cv_values_.ravel() / np.sqrt(sw)
                  + np.average(y, weights=sw)))**2
assert np.allclose(squared_error, custom)
```

I think that for a user, it would be easier to get directly the predictions in
the original space, and not need to do this post-processing of `cv_values_`.

Should we rescale the cv values and add the intercept during `fit`?

Hints from issue discussion:
ping @agramfort 
yes `cv_values_` should be put in native space and fixed when using a custom scoring.

Commit message: address thomas review + n_target test

Patch to review:
diff --git a/doc/whats_new/v1.6.rst b/doc/whats_new/v1.6.rst
index 22a0d7acfd24e..bf141b582b365 100644
--- a/doc/whats_new/v1.6.rst
+++ b/doc/whats_new/v1.6.rst
@@ -257,6 +257,12 @@ Changelog
   for the calculation of test scores.
   :pr:`29419` by :user:`Shruti Nath <snath-xoc>`.
 
+- |Fix| :class:`linear_model.RidgeCV` now properly use predictions in the original
+  scale and store them in `cv_results_` when `scoring != None`.
+  :pr:`29842` by :user:`Guillaume Lemaitre <glemaitre>`,
+  :user:`Jérôme Dockes <jeromedockes>` and
+  :user:`Hanmin Qin <qinhanmin2014>`.
+
 - |API| Deprecates `copy_X` in :class:`linear_model.TheilSenRegressor` as the parameter
   has no effect. `copy_X` will be removed in 1.8.
   :pr:`29105` by :user:`Adam Li <adam2392>`.
diff --git a/sklearn/linear_model/_ridge.py b/sklearn/linear_model/_ridge.py
index 7b72efa810840..a634e84aa12fd 100644
--- a/sklearn/linear_model/_ridge.py
+++ b/sklearn/linear_model/_ridge.py
@@ -2129,6 +2129,7 @@ def fit(self, X, y, sample_weight=None, score_params=None):
 
         self.alphas = np.asarray(self.alphas)
 
+        unscaled_y = y
         X, y, X_offset, y_offset, X_scale = _preprocess_data(
             X,
             y,
@@ -2178,13 +2179,20 @@ def fit(self, X, y, sample_weight=None, score_params=None):
                     self.cv_results_[:, i] = squared_errors.ravel()
             else:
                 predictions = y - (c / G_inverse_diag)
+                # represent the predictions in the original scale
+                if predictions.ndim > 1:
+                    predictions /= sqrt_sw[:, None]
+                else:
+                    predictions /= sqrt_sw
+                predictions += y_offset
+
                 if self.store_cv_results:
                     self.cv_results_[:, i] = predictions.ravel()
 
                 score_params = score_params or {}
                 alpha_score = self._score(
                     predictions=predictions,
-                    y=y,
+                    y=unscaled_y,
                     n_y=n_y,
                     scorer=scorer,
                     score_params=score_params,
diff --git a/sklearn/linear_model/tests/test_ridge.py b/sklearn/linear_model/tests/test_ridge.py
index cb42e1a473647..b456b900da13e 100644
--- a/sklearn/linear_model/tests/test_ridge.py
+++ b/sklearn/linear_model/tests/test_ridge.py
@@ -859,7 +859,9 @@ def test_ridge_loo_cv_asym_scoring():
     loo_ridge.fit(X, y)
     gcv_ridge.fit(X, y)
 
-    assert gcv_ridge.alpha_ == pytest.approx(loo_ridge.alpha_)
+    assert gcv_ridge.alpha_ == pytest.approx(
+        loo_ridge.alpha_
+    ), f"gcv_ridge.alpha_: {gcv_ridge.alpha_}, loo_ridge.alpha_: {loo_ridge.alpha_}"
     assert_allclose(gcv_ridge.coef_, loo_ridge.coef_, rtol=1e-3)
     assert_allclose(gcv_ridge.intercept_, loo_ridge.intercept_, rtol=1e-3)
 
@@ -2252,6 +2254,52 @@ def test_ridge_cv_values_deprecated():
         ridge.cv_values_
 
 
+@pytest.mark.parametrize("with_sample_weight", [False, True])
+@pytest.mark.parametrize("fit_intercept", [False, True])
+@pytest.mark.parametrize("n_targets", [1, 2])
+def test_ridge_cv_results_predictions(with_sample_weight, fit_intercept, n_targets):
+    """Check that the predictions stored in `cv_results_` are on the original scale.
+
+    The GCV approach works on scaled data: centered by an offset and scaled by the
+    squared root of the sample weights. Thus, previous to compute scores, the
+    predictions need to be scaled back to the original scale. Those predictions are the
+    one stored in `cv_results_`.
+
+    In this test, we check that the internal predictions stored in `cv_results_` are
+    equivalent to a naive LOO-CV grid-search with a `Ridge` estimator.
+
+    Non-regression test for:
+    https://github.com/scikit-learn/scikit-learn/issues/13998
+    """
+    X, y = make_regression(
+        n_samples=100, n_features=10, n_targets=n_targets, random_state=0
+    )
+    sample_weight = np.ones(shape=(X.shape[0],))
+    if with_sample_weight:
+        sample_weight[::2] = 0.5
+
+    alphas = (0.1, 1.0, 10.0)
+
+    # scoring should be set to store predictions and not the squared error
+    ridge_cv = RidgeCV(
+        alphas=alphas,
+        scoring="neg_mean_squared_error",
+        fit_intercept=fit_intercept,
+        store_cv_results=True,
+    )
+    ridge_cv.fit(X, y, sample_weight=sample_weight)
+
+    # manual grid-search with a `Ridge` estimator
+    predictions = np.empty(shape=(*y.shape, len(alphas)))
+    cv = LeaveOneOut()
+    for alpha_idx, alpha in enumerate(alphas):
+        for idx, (train_idx, test_idx) in enumerate(cv.split(X, y)):
+            ridge = Ridge(alpha=alpha, fit_intercept=fit_intercept)
+            ridge.fit(X[train_idx], y[train_idx], sample_weight[train_idx])
+            predictions[idx, ..., alpha_idx] = ridge.predict(X[test_idx])
+    assert_allclose(ridge_cv.cv_results_, predictions)
+
+
 def test_ridge_cv_multioutput_sample_weight(global_random_seed):
     """Check that `RidgeCV` works properly with multioutput and sample_weight
     when `scoring != None`.


Code review comments:
diff_hunk: @@ -2178,13 +2179,20 @@ def fit(self, X, y, sample_weight=None, score_params=None):
                     self.cv_results_[:, i] = squared_errors.ravel()
             else:
                 predictions = y - (c / G_inverse_diag)
+                # represent the predictions in the original scale
path: sklearn/linear_model/_ridge.py
text: @user1:
```suggestion
                # Rescale predictions back to original scale
```

diff_hunk: @@ -859,7 +859,9 @@ def test_ridge_loo_cv_asym_scoring():
     loo_ridge.fit(X, y)
     gcv_ridge.fit(X, y)
 
-    assert gcv_ridge.alpha_ == pytest.approx(loo_ridge.alpha_)
+    assert gcv_ridge.alpha_ == pytest.approx(
+        loo_ridge.alpha_
+    ), f"gcv_ridge.alpha_: {gcv_ridge.alpha_}, loo_ridge.alpha_: {loo_ridge.alpha_}"
path: sklearn/linear_model/tests/test_ridge.py
text: @user1:
```suggestion
    ), f"{gcv_ridge.alpha_=}, {loo_ridge.alpha_=}"
```
if you want to make it shorter.

diff_hunk: @@ -2178,13 +2179,20 @@ def fit(self, X, y, sample_weight=None, score_params=None):
                     self.cv_results_[:, i] = squared_errors.ravel()
             else:
                 predictions = y - (c / G_inverse_diag)
+                # represent the predictions in the original scale
+                if predictions.ndim > 1:
path: sklearn/linear_model/_ridge.py
text: @user1:
Maybe add
```python 
if sample_weight is not None:
```

@author:
Let's add to have a small optimization.