scikit-learn
diff --git a/‎doc/whats_new/v0.24.rst
Lines changed: 8 additions & 0 deletions b/‎doc/whats_new/v0.24.rst
Lines changed: 8 additions & 0 deletions
diff --git a/‎sklearn/metrics/_regression.py
Lines changed: 8 additions & 3 deletions b/‎sklearn/metrics/_regression.py
Lines changed: 8 additions & 3 deletions
diff --git a/‎sklearn/metrics/tests/test_regression.py
Lines changed: 2 additions & 2 deletions b/‎sklearn/metrics/tests/test_regression.py
Lines changed: 2 additions & 2 deletions
@@ -61,6 +61,14 @@ Changelog
   change since `None` was defaulting to these values already.
   :pr:`16493` by :user:`Darshan N <DarshanGowda0>`.
 
+:mod:`sklearn.metrics`
+......................
+
+- |Fix| Fixed a bug in :func:`metrics.mean_squared_error` where the
+  average of multiple RMSE values was incorrectly calculated as the root of the
+  average of multiple MSE values.
+  :pr:`17309` by :user:`Swier Heeres <swierh>`
+
 :mod:`sklearn.model_selection`
 ..............................
 
 
@@ -244,6 +244,8 @@ def mean_squared_error(y_true, y_pred, *,
     >>> y_pred = [[0, 2],[-1, 2],[8, -5]]
     >>> mean_squared_error(y_true, y_pred)
     0.708...
+    >>> mean_squared_error(y_true, y_pred, squared=False)
+    0.822...
     >>> mean_squared_error(y_true, y_pred, multioutput='raw_values')
     array([0.41666667, 1.        ])
     >>> mean_squared_error(y_true, y_pred, multioutput=[0.3, 0.7])
@@ -255,15 +257,18 @@ def mean_squared_error(y_true, y_pred, *,
     check_consistent_length(y_true, y_pred, sample_weight)
     output_errors = np.average((y_true - y_pred) ** 2, axis=0,
                                weights=sample_weight)
+
+    if not squared:
+        output_errors = np.sqrt(output_errors)
+
     if isinstance(multioutput, str):
         if multioutput == 'raw_values':
-            return output_errors if squared else np.sqrt(output_errors)
+            return output_errors
         elif multioutput == 'uniform_average':
             # pass None as weights to np.average: uniform mean
             multioutput = None
 
-    mse = np.average(output_errors, weights=multioutput)
-    return mse if squared else np.sqrt(mse)
+    return np.average(output_errors, weights=multioutput)
 
 
 @_deprecate_positional_args
 
@@ -76,7 +76,7 @@ def test_multioutput_regression():
     assert_almost_equal(error, (1. / 3 + 2. / 3 + 2. / 3) / 4.)
 
     error = mean_squared_error(y_true, y_pred, squared=False)
-    assert_almost_equal(error, 0.645, decimal=2)
+    assert_almost_equal(error, 0.454, decimal=2)
 
     error = mean_squared_log_error(y_true, y_pred)
     assert_almost_equal(error, 0.200, decimal=2)
@@ -258,7 +258,7 @@ def test_regression_custom_weights():
     evsw = explained_variance_score(y_true, y_pred, multioutput=[0.4, 0.6])
 
     assert_almost_equal(msew, 0.39, decimal=2)
-    assert_almost_equal(rmsew, 0.62, decimal=2)
+    assert_almost_equal(rmsew, 0.59, decimal=2)
     assert_almost_equal(maew, 0.475, decimal=3)
     assert_almost_equal(rw, 0.94, decimal=2)
     assert_almost_equal(evsw, 0.94, decimal=2)