Repository: scikit-learn/scikit-learn
Pull request: ENH add force_finite to control the behaviour of constant target in explained variance and R2
#### Reference Issues/PRs

Fixes #16857

#### What does this implement/fix? Explain your changes.

Clarified both in the documentation and in the implementation, the fact that the default behaviour of `r2_score` is different from the actual R² formula. Added an explicit `fix_when_y_true_is_constant` arg (name to debate) with conservative default (preserving current sklearn's behaviour) allowing users to remove the fix and use the actual R² formula if they wish.

EDIT: I just found out that the same goes for Explained Variance so I'll propose to improve it too.


Problem statement:
`r2_score`: improved support of edge case where `y_true` is constant
As for many users, when I first encountered sklearn a few years ago I was surprised to see that there is a hardcoded behaviour in the `r2_score` to handle the case where training data is constant.

```python
import sys
from sklearn.metrics import r2_score

r2_constant_train_nonexact_preds = r2_score(y_true=[-1, -1, -1, -1], 
                                            y_pred=[-1, -1, -1, -1 + sys.float_info.epsilon])
print("Constant train, non-exact predictions. Score should be -Inf. Found: %s" % r2_constant_train_nonexact_preds)

r2_constant_train_exact_preds = r2_score(y_true=[-1, -1, -1, -1], 
                                         y_pred=[-1, -1, -1, -1])
print("Constant train, exact predictions. Score should be NaN. Found: %s" % r2_constant_train_exact_preds)
```

yields:

```
Constant train, non-exact predictions. Score should be -Inf. Found: 0.0
Constant train, exact predictions. Score should be NaN. Found: 1.0
```

Of course this has been done to make sure that the scoring function always return finite scores. Still, I believe that this should be improved for readability and usability. 

I suggest to add a new `ensure_finite=True` argument to the `r2_score` function. The documentation of this argument would explicitly state that when `True` (default), the score is turned to `1.0` or `0.0` when the data is constant so that model comparisons do not fail (typically in cross-validation). This new argument would allow users to set it to `False` to get the "true" R² definition, for other kind of analysis.

What do you think ? I can provide a PR if you think that this brings value.

Hints from issue discussion:
Update on this ticket before pushing my PR: there also seem to be an ambiguity in the docstring

> A constant model that always predicts the expected value of y, disregarding the input features, would get a R^2 score of 0.0.

For me this is incorrect, but maybe this is because the sentence is not clear (does "expected value" mean average ? and even so, if true data is constant this does not hold). I'll propose a reformulation in the PR as well.

Commit message: Update sklearn/metrics/_regression.py

Co-authored-by: Julien Jerphanion <git@jjerphan.xyz>

Patch to review:
diff --git a/doc/modules/model_evaluation.rst b/doc/modules/model_evaluation.rst
index 3dd7313bb9b59..48e495b719898 100644
--- a/doc/modules/model_evaluation.rst
+++ b/doc/modules/model_evaluation.rst
@@ -2004,6 +2004,18 @@ then the explained variance is estimated as follow:
 
 The best possible score is 1.0, lower values are worse.
 
+Note: when the prediction residuals have zero mean (perfectly unbiased
+model), the Explained Variance score and the :ref:`r2_score` are identical.
+
+In the particular case where the true y is constant, the Explained Variance
+score is not finite: it is either ``NaN`` (perfect predictions) or ``-Inf``
+(imperfect predictions). Such non-finite scores may prevent correct model
+optimization such as grid-search cross-validation to be performed correctly.
+For this reason the default behaviour of :func:`explained_variance_score` is
+to replace them with 1.0 (perfect predictions) or 0.0 (imperfect predictions).
+You can set the ``force_finite`` parameter to ``False`` to prevent this
+fix from happening and fallback on the original Explained Variance score.
+
 Here is a small example of usage of the :func:`explained_variance_score`
 function::
 
@@ -2018,6 +2030,18 @@ function::
     array([0.967..., 1.        ])
     >>> explained_variance_score(y_true, y_pred, multioutput=[0.3, 0.7])
     0.990...
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2]
+    >>> explained_variance_score(y_true, y_pred)
+    1.0
+    >>> explained_variance_score(y_true, y_pred, force_finite=False)
+    nan
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2 + 1e-8]
+    >>> explained_variance_score(y_true, y_pred)
+    0.0
+    >>> explained_variance_score(y_true, y_pred, force_finite=False)
+    -inf
 
 .. _max_error:
 
@@ -2240,8 +2264,11 @@ predicted by the model, through the proportion of explained variance.
 As such variance is dataset dependent, R² may not be meaningfully comparable
 across different datasets. Best possible score is 1.0 and it can be negative
 (because the model can be arbitrarily worse). A constant model that always
-predicts the expected value of y, disregarding the input features, would get a
-R² score of 0.0.
+predicts the expected (average) value of y, disregarding the input features,
+would get a :math:`R^2` score of 0.0.
+
+Note: when the prediction residuals have zero mean (perfectly unbiased
+model), the :math:`R^2` score and the :ref:`explained_variance_score` are identical.
 
 If :math:`\hat{y}_i` is the predicted value of the :math:`i`-th sample
 and :math:`y_i` is the corresponding true value for total :math:`n` samples,
@@ -2256,6 +2283,14 @@ where :math:`\bar{y} = \frac{1}{n} \sum_{i=1}^{n} y_i` and :math:`\sum_{i=1}^{n}
 Note that :func:`r2_score` calculates unadjusted R² without correcting for
 bias in sample variance of y.
 
+In the particular case where the true y is constant, the :math:`R^2` score is
+not finite: it is either ``NaN`` (perfect predictions) or ``-Inf`` (imperfect
+predictions). Such non-finite scores may prevent correct model optimization
+such as grid-search cross-validation to be performed correctly. For this reason
+the default behaviour of :func:`r2_score` is to replace them with 1.0 (perfect
+predictions) or 0.0 (imperfect predictions). If ``force_finite``
+is set to ``False``, this score fallbacks on the original :math:`R^2` definition.
+
 Here is a small example of usage of the :func:`r2_score` function::
 
   >>> from sklearn.metrics import r2_score
@@ -2275,7 +2310,18 @@ Here is a small example of usage of the :func:`r2_score` function::
   array([0.965..., 0.908...])
   >>> r2_score(y_true, y_pred, multioutput=[0.3, 0.7])
   0.925...
-
+  >>> y_true = [-2, -2, -2]
+  >>> y_pred = [-2, -2, -2]
+  >>> r2_score(y_true, y_pred)
+  1.0
+  >>> r2_score(y_true, y_pred, force_finite=False)
+  nan
+  >>> y_true = [-2, -2, -2]
+  >>> y_pred = [-2, -2, -2 + 1e-8]
+  >>> r2_score(y_true, y_pred)
+  0.0
+  >>> r2_score(y_true, y_pred, force_finite=False)
+  -inf
 
 .. topic:: Example:
 
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
index 1a9f773ce08df..f0b89313e45ce 100644
--- a/doc/whats_new/v1.1.rst
+++ b/doc/whats_new/v1.1.rst
@@ -63,6 +63,12 @@ Changelog
 :mod:`sklearn.metrics`
 ......................
 
+- |Feature| :func:`r2_score` and :func:`explained_variance_score` have a new
+  `force_finite` parameter. Setting this parameter to `False` will return the
+  actual non-finite score in case of perfect predictions or constant `y_true`,
+  instead of the finite approximation (`1.0` and `0.0` respectively) currently
+  returned by default. :pr:`17266` by :user:`Sylvain Marié <smarie>`.
+
 - |API| :class:`metrics.DistanceMetric` has been moved from
   :mod:`sklearn.neighbors` to :mod:`sklearn.metric`.
   Using `neighbors.DistanceMetric` for imports is still valid for
diff --git a/sklearn/metrics/_regression.py b/sklearn/metrics/_regression.py
index dd4c4c52f55fd..7183abd8e75e9 100644
--- a/sklearn/metrics/_regression.py
+++ b/sklearn/metrics/_regression.py
@@ -22,6 +22,7 @@
 #          Christian Lorentzen <lorentzen.ch@googlemail.com>
 #          Ashutosh Hathidara <ashutoshhathidara98@gmail.com>
 #          Uttam kumar <bajiraouttamsinha@gmail.com>
+#          Sylvain Marie <sylvain.marie@se.com>
 # License: BSD 3 clause
 
 import warnings
@@ -594,12 +595,30 @@ def median_absolute_error(
 
 
 def explained_variance_score(
-    y_true, y_pred, *, sample_weight=None, multioutput="uniform_average"
+    y_true,
+    y_pred,
+    *,
+    sample_weight=None,
+    multioutput="uniform_average",
+    force_finite=True,
 ):
     """Explained variance regression score function.
 
     Best possible score is 1.0, lower values are worse.
 
+    In the particular case when the true y is constant, the explained variance
+    score is not finite: it is either ``NaN`` (perfect predictions) or
+    ``-Inf`` (imperfect predictions). To prevent such non-finite numbers to
+    pollute higher-level experiments such as a grid search cross-validation,
+    by default these cases are replaced with 1.0 (perfect predictions) or 0.0
+    (imperfect predictions) respectively. If ``force_finite``
+    is set to ``False``, this score falls back on the original :math:`R^2`
+    definition.
+
+    Note: when the prediction residuals have zero mean (perfectly unbiased
+    model), the Explained Variance score is identical to the
+    :func:`R^2 score <r2_score>`.
+
     Read more in the :ref:`User Guide <explained_variance_score>`.
 
     Parameters
@@ -628,6 +647,14 @@ def explained_variance_score(
             Scores of all outputs are averaged, weighted by the variances
             of each individual output.
 
+    force_finite : boolean, optional
+        Flag indicating if ``NaN`` and ``-Inf`` scores resulting from constant
+        data should be replaced with real numbers (``1.0`` if prediction is
+        perfect, ``0.0`` otherwise). Default is ``True``, a convenient setting
+        for model optimization procedures (e.g. grid search cross-validation).
+
+        .. versionadded:: 1.1.0
+
     Returns
     -------
     score : float or ndarray of floats
@@ -648,6 +675,18 @@ def explained_variance_score(
     >>> y_pred = [[0, 2], [-1, 2], [8, -5]]
     >>> explained_variance_score(y_true, y_pred, multioutput='uniform_average')
     0.983...
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2]
+    >>> explained_variance_score(y_true, y_pred)
+    1.0
+    >>> explained_variance_score(y_true, y_pred, force_finite=False)
+    nan
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2 + 1e-8]
+    >>> explained_variance_score(y_true, y_pred)
+    0.0
+    >>> explained_variance_score(y_true, y_pred, force_finite=False)
+    -inf
     """
     y_type, y_true, y_pred, multioutput = _check_reg_targets(
         y_true, y_pred, multioutput
@@ -662,35 +701,41 @@ def explained_variance_score(
     y_true_avg = np.average(y_true, weights=sample_weight, axis=0)
     denominator = np.average((y_true - y_true_avg) ** 2, weights=sample_weight, axis=0)
 
-    nonzero_numerator = numerator != 0
-    nonzero_denominator = denominator != 0
-    valid_score = nonzero_numerator & nonzero_denominator
-    output_scores = np.ones(y_true.shape[1])
-
-    output_scores[valid_score] = 1 - (numerator[valid_score] / denominator[valid_score])
-    output_scores[nonzero_numerator & ~nonzero_denominator] = 0.0
-    if isinstance(multioutput, str):
-        if multioutput == "raw_values":
-            # return scores individually
-            return output_scores
-        elif multioutput == "uniform_average":
-            # passing to np.average() None as weights results is uniform mean
-            avg_weights = None
-        elif multioutput == "variance_weighted":
-            avg_weights = denominator
-    else:
-        avg_weights = multioutput
-
-    return np.average(output_scores, weights=avg_weights)
+    return _assemble_r2_explained_variance(
+        numerator=numerator,
+        denominator=denominator,
+        nb_outputs=y_true.shape[1],
+        multioutput=multioutput,
+        force_finite=force_finite,
+    )
 
 
-def r2_score(y_true, y_pred, *, sample_weight=None, multioutput="uniform_average"):
+def r2_score(
+    y_true,
+    y_pred,
+    *,
+    sample_weight=None,
+    multioutput="uniform_average",
+    force_finite=True,
+):
     """:math:`R^2` (coefficient of determination) regression score function.
 
     Best possible score is 1.0 and it can be negative (because the
-    model can be arbitrarily worse). A constant model that always
-    predicts the expected value of y, disregarding the input features,
-    would get a :math:`R^2` score of 0.0.
+    model can be arbitrarily worse). In the general case when the true y is
+    non-constant, a constant model that always predicts the average y
+    disregarding the input features would get a :math:`R^2` score of 0.0.
+
+    In the particular case when the true y is constant, the :math:`R^2` score
+    is not finite: it is either ``NaN`` (perfect predictions) or ``-Inf``
+    (imperfect predictions). To prevent such non-finite numbers to pollute
+    higher-level experiments such as a grid search cross-validation, by default
+    these cases are replaced with 1.0 (perfect predictions) or 0.0 (imperfect
+    predictions) respectively. You can set ``force_finite`` to ``False`` to
+    prevent this fix from happening.
+
+    Note: when the prediction residuals have zero mean (perfectly unbiased
+    model), the :math:`R^2` score is identical to the
+    :func:`Explained Variance score <explained_variance_score>`.
 
     Read more in the :ref:`User Guide <r2_score>`.
 
@@ -725,6 +770,14 @@ def r2_score(y_true, y_pred, *, sample_weight=None, multioutput="uniform_average
         .. versionchanged:: 0.19
             Default value of multioutput is 'uniform_average'.
 
+    force_finite : boolean, optional
+        Flag indicating if ``NaN`` and ``-Inf`` scores resulting from constant
+        data should be replaced with real numbers (``1.0`` if prediction is
+        perfect, ``0.0`` otherwise). Default is ``True``, a convenient settings
+        for model optimization procedures (e.g. grid search cross-validation).
+
+        .. versionadded:: 1.1.0
+
     Returns
     -------
     z : float or ndarray of floats
@@ -770,6 +823,18 @@ def r2_score(y_true, y_pred, *, sample_weight=None, multioutput="uniform_average
     >>> y_pred = [3, 2, 1]
     >>> r2_score(y_true, y_pred)
     -3.0
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2]
+    >>> r2_score(y_true, y_pred)
+    1.0
+    >>> r2_score(y_true, y_pred, force_finite=False)
+    nan
+    >>> y_true = [-2, -2, -2]
+    >>> y_pred = [-2, -2, -2 + 1e-8]
+    >>> r2_score(y_true, y_pred)
+    0.0
+    >>> r2_score(y_true, y_pred, force_finite=False)
+    -inf
     """
     y_type, y_true, y_pred, multioutput = _check_reg_targets(
         y_true, y_pred, multioutput
@@ -791,25 +856,50 @@ def r2_score(y_true, y_pred, *, sample_weight=None, multioutput="uniform_average
     denominator = (
         weight * (y_true - np.average(y_true, axis=0, weights=sample_weight)) ** 2
     ).sum(axis=0, dtype=np.float64)
-    nonzero_denominator = denominator != 0
-    nonzero_numerator = numerator != 0
-    valid_score = nonzero_denominator & nonzero_numerator
-    output_scores = np.ones([y_true.shape[1]])
-    output_scores[valid_score] = 1 - (numerator[valid_score] / denominator[valid_score])
-    # arbitrary set to zero to avoid -inf scores, having a constant
-    # y_true is not interesting for scoring a regression anyway
-    output_scores[nonzero_numerator & ~nonzero_denominator] = 0.0
+
+    return _assemble_r2_explained_variance(
+        numerator=numerator,
+        denominator=denominator,
+        nb_outputs=y_true.shape[1],
+        multioutput=multioutput,
+        force_finite=force_finite,
+    )
+
+
+def _assemble_r2_explained_variance(
+    numerator, denominator, nb_outputs, multioutput, force_finite
+):
+    """Common part used by explained variance score and :math:`R^2` score"""
+
+    if not force_finite:
+        # Standard formula, that may lead to NaN or -Inf
+        output_scores = 1 - (numerator / denominator)
+    else:
+        nonzero_denominator = denominator != 0
+        nonzero_numerator = numerator != 0
+        # Default = Zero Numerator = perfect predictions. Set to 1.0
+        # (note: even if denominator is zero, thus avoiding NaN scores)
+        output_scores = np.ones([nb_outputs])
+        # Non-zero Numerator and Non-zero Denominator: use the formula
+        valid_score = nonzero_denominator & nonzero_numerator
+        output_scores[valid_score] = 1 - (
+            numerator[valid_score] / denominator[valid_score]
+        )
+        # Non-zero Numerator and Zero Denominator:
+        # arbitrary set to 0.0 to avoid -inf scores
+        output_scores[nonzero_numerator & ~nonzero_denominator] = 0.0
+
     if isinstance(multioutput, str):
         if multioutput == "raw_values":
             # return scores individually
             return output_scores
         elif multioutput == "uniform_average":
-            # passing None as weights results is uniform mean
+            # Passing None as weights to np.average results is uniform mean
             avg_weights = None
         elif multioutput == "variance_weighted":
             avg_weights = denominator
-            # avoid fail on constant y or one-element arrays
-            if not np.any(nonzero_denominator):
+            if force_finite and not np.any(nonzero_denominator):
+                # Avoid failing on constant y or one-element arrays.
                 if not np.any(nonzero_numerator):
                     return 1.0
                 else:
diff --git a/sklearn/metrics/tests/test_regression.py b/sklearn/metrics/tests/test_regression.py
index b66ce18ec8da4..77b7dd92905eb 100644
--- a/sklearn/metrics/tests/test_regression.py
+++ b/sklearn/metrics/tests/test_regression.py
@@ -50,7 +50,11 @@ def test_regression_metrics(n_samples=50):
     assert mape > 1e6
     assert_almost_equal(max_error(y_true, y_pred), 1.0)
     assert_almost_equal(r2_score(y_true, y_pred), 0.995, 2)
+    assert_almost_equal(r2_score(y_true, y_pred, force_finite=False), 0.995, 2)
     assert_almost_equal(explained_variance_score(y_true, y_pred), 1.0)
+    assert_almost_equal(
+        explained_variance_score(y_true, y_pred, force_finite=False), 1.0
+    )
     assert_almost_equal(
         mean_tweedie_deviance(y_true, y_pred, power=0),
         mean_squared_error(y_true, y_pred),
@@ -135,8 +139,30 @@ def test_multioutput_regression():
     error = r2_score(y_true, y_pred, multioutput="uniform_average")
     assert_almost_equal(error, -0.875)
 
+    # Setting force_finite=False results in the nan for 4th output propagating
+    error = r2_score(
+        y_true, y_pred, multioutput="variance_weighted", force_finite=False
+    )
+    assert_almost_equal(error, np.nan)
+    error = r2_score(y_true, y_pred, multioutput="uniform_average", force_finite=False)
+    assert_almost_equal(error, np.nan)
+
+    # Dropping the 4th output to check `force_finite=False` for nominal
+    y_true = y_true[:, :-1]
+    y_pred = y_pred[:, :-1]
+    error = r2_score(y_true, y_pred, multioutput="variance_weighted")
+    error2 = r2_score(
+        y_true, y_pred, multioutput="variance_weighted", force_finite=False
+    )
+    assert_almost_equal(error, error2)
+    error = r2_score(y_true, y_pred, multioutput="uniform_average")
+    error2 = r2_score(y_true, y_pred, multioutput="uniform_average", force_finite=False)
+    assert_almost_equal(error, error2)
+
 
 def test_regression_metrics_at_limits():
+    # Single-sample case
+    # Note: for r2 and d2_tweedie see also test_regression_single_sample
     assert_almost_equal(mean_squared_error([0.0], [0.0]), 0.0)
     assert_almost_equal(mean_squared_error([0.0], [0.0], squared=False), 0.0)
     assert_almost_equal(mean_squared_log_error([0.0], [0.0]), 0.0)
@@ -146,7 +172,17 @@ def test_regression_metrics_at_limits():
     assert_almost_equal(median_absolute_error([0.0], [0.0]), 0.0)
     assert_almost_equal(max_error([0.0], [0.0]), 0.0)
     assert_almost_equal(explained_variance_score([0.0], [0.0]), 1.0)
+
+    # Perfect cases
     assert_almost_equal(r2_score([0.0, 1], [0.0, 1]), 1.0)
+
+    # Non-finite cases
+    # R² and explained variance have a fix by default for non-finite cases
+    for s in (r2_score, explained_variance_score):
+        assert_almost_equal(s([0, 0], [1, -1]), 0.0)
+        assert_almost_equal(s([0, 0], [1, -1], force_finite=False), -np.inf)
+        assert_almost_equal(s([1, 1], [1, 1]), 1.0)
+        assert_almost_equal(s([1, 1], [1, 1], force_finite=False), np.nan)
     msg = (
         "Mean Squared Logarithmic Error cannot be used when targets "
         "contain negative values."
@@ -270,6 +306,9 @@ def test_regression_multioutput_array():
     mape = mean_absolute_percentage_error(y_true, y_pred, multioutput="raw_values")
     r = r2_score(y_true, y_pred, multioutput="raw_values")
     evs = explained_variance_score(y_true, y_pred, multioutput="raw_values")
+    evs2 = explained_variance_score(
+        y_true, y_pred, multioutput="raw_values", force_finite=False
+    )
 
     assert_array_almost_equal(mse, [0.125, 0.5625], decimal=2)
     assert_array_almost_equal(mae, [0.25, 0.625], decimal=2)
@@ -277,6 +316,7 @@ def test_regression_multioutput_array():
     assert_array_almost_equal(mape, [0.0778, 0.2262], decimal=2)
     assert_array_almost_equal(r, [0.95, 0.93], decimal=2)
     assert_array_almost_equal(evs, [0.95, 0.93], decimal=2)
+    assert_array_almost_equal(evs2, [0.95, 0.93], decimal=2)
 
     # mean_absolute_error and mean_squared_error are equal because
     # it is a binary problem.
@@ -300,17 +340,38 @@ def test_regression_multioutput_array():
         [[0, -1], [0, 1]], [[2, 2], [1, 1]], multioutput="raw_values"
     )
     assert_array_almost_equal(evs, [0, -1.25], decimal=2)
+    evs2 = explained_variance_score(
+        [[0, -1], [0, 1]],
+        [[2, 2], [1, 1]],
+        multioutput="raw_values",
+        force_finite=False,
+    )
+    assert_array_almost_equal(evs2, [-np.inf, -1.25], decimal=2)
 
     # Checking for the condition in which both numerator and denominator is
     # zero.
-    y_true = [[1, 3], [-1, 2]]
-    y_pred = [[1, 4], [-1, 1]]
+    y_true = [[1, 3], [1, 2]]
+    y_pred = [[1, 4], [1, 1]]
     r2 = r2_score(y_true, y_pred, multioutput="raw_values")
     assert_array_almost_equal(r2, [1.0, -3.0], decimal=2)
     assert np.mean(r2) == r2_score(y_true, y_pred, multioutput="uniform_average")
+    r22 = r2_score(y_true, y_pred, multioutput="raw_values", force_finite=False)
+    assert_array_almost_equal(r22, [np.nan, -3.0], decimal=2)
+    assert_almost_equal(
+        np.mean(r22),
+        r2_score(y_true, y_pred, multioutput="uniform_average", force_finite=False),
+    )
+
     evs = explained_variance_score(y_true, y_pred, multioutput="raw_values")
     assert_array_almost_equal(evs, [1.0, -3.0], decimal=2)
     assert np.mean(evs) == explained_variance_score(y_true, y_pred)
+    evs2 = explained_variance_score(
+        y_true, y_pred, multioutput="raw_values", force_finite=False
+    )
+    assert_array_almost_equal(evs2, [np.nan, -3.0], decimal=2)
+    assert_almost_equal(
+        np.mean(evs2), explained_variance_score(y_true, y_pred, force_finite=False)
+    )
 
     # Handling msle separately as it does not accept negative inputs.
     y_true = np.array([[0.5, 1], [1, 2], [7, 6]])
@@ -332,6 +393,9 @@ def test_regression_custom_weights():
     mapew = mean_absolute_percentage_error(y_true, y_pred, multioutput=[0.4, 0.6])
     rw = r2_score(y_true, y_pred, multioutput=[0.4, 0.6])
     evsw = explained_variance_score(y_true, y_pred, multioutput=[0.4, 0.6])
+    evsw2 = explained_variance_score(
+        y_true, y_pred, multioutput=[0.4, 0.6], force_finite=False
+    )
 
     assert_almost_equal(msew, 0.39, decimal=2)
     assert_almost_equal(rmsew, 0.59, decimal=2)
@@ -339,6 +403,7 @@ def test_regression_custom_weights():
     assert_almost_equal(mapew, 0.1668, decimal=2)
     assert_almost_equal(rw, 0.94, decimal=2)
     assert_almost_equal(evsw, 0.94, decimal=2)
+    assert_almost_equal(evsw2, 0.94, decimal=2)
 
     # Handling msle separately as it does not accept negative inputs.
     y_true = np.array([[0.5, 1], [1, 2], [7, 6]])


Code review comments:
diff_hunk: @@ -725,6 +770,14 @@ def r2_score(y_true, y_pred, *, sample_weight=None, multioutput="uniform_average
         .. versionchanged:: 0.19
             Default value of multioutput is 'uniform_average'.
 
+    force_finite : boolean, optional
path: sklearn/metrics/_regression.py
text: @author:
```suggestion
    force_finite : bool, default=True
```

diff_hunk: @@ -628,6 +647,14 @@ def explained_variance_score(
             Scores of all outputs are averaged, weighted by the variances
             of each individual output.
 
+    force_finite : boolean, optional
path: sklearn/metrics/_regression.py
text: @author:
```suggestion
    force_finite : bool, default=True
```