Repository: sktime/sktime
Pull request #8274 (2025-05-21 14:10:53)

[BUG] fix recursive reduction tests
This fixes some bugs in recursive reduction tests, where instead of `assert` statements there were `return` statements or simply no statements, thus resulting in no check where there should have been one.

FYI @PranavBhatP, it would be appreciated if you could review this since you created these tests in #7770.

PranavBhatP: Hi @fkiraly , your changes are valid and I don't think they will need a separate review. This must have gone unnoticed while merging PR. Thanks for pointing this out and apologies for the mistake :)
fkiraly: Well, they seem to pass, so all is good!

Patch:
diff --git a/sktime/forecasting/reconcile.py b/sktime/forecasting/reconcile.py
index ee2695f8c..8355809f4 100644
--- a/sktime/forecasting/reconcile.py
+++ b/sktime/forecasting/reconcile.py
@@ -58,6 +58,9 @@ class ReconcilerForecaster(BaseForecaster):
         * If True, prediction data frames include total values at ``__total`` levels
         * If False, prediction data frames are returned without ``__total`` levels
 
+    alpha: float default=0
+        Optional regularization parameter to avoid singular covariance matrix
+
     See Also
     --------
     Aggregator
@@ -121,10 +124,11 @@ class ReconcilerForecaster(BaseForecaster):
     METHOD_LIST = ["mint_cov", "mint_shrink", "wls_var"] + TRFORM_LIST
     RETURN_TOTALS_LIST = [True, False]
 
-    def __init__(self, forecaster, method="mint_shrink", return_totals=True):
+    def __init__(self, forecaster, method="mint_shrink", return_totals=True, alpha=0):
         self.forecaster = forecaster
         self.method = method
         self.return_totals = return_totals
+        self.alpha = alpha
 
         super().__init__()
 
@@ -402,6 +406,7 @@ class ReconcilerForecaster(BaseForecaster):
             np.dot(
                 inv(
                     np.dot(np.transpose(self.s_matrix), np.dot(cov_mat, self.s_matrix))
+                    + self.alpha * np.eye(self.s_matrix.shape[1])
                 ),
                 np.dot(np.transpose(self.s_matrix), cov_mat),
             )
@@ -425,6 +430,13 @@ class ReconcilerForecaster(BaseForecaster):
     def get_test_params(cls):
         """Return testing parameter settings for the estimator.
 
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+            There are currently no reserved values for clusterers.
+
         Returns
         -------
         params : dict, default = {}
@@ -441,9 +453,11 @@ class ReconcilerForecaster(BaseForecaster):
             {
                 "forecaster": FORECASTER,
                 "method": x,
+                "alpha": alpha,
                 "return_totals": totals,
             }
             for x in cls.METHOD_LIST
             for totals in cls.RETURN_TOTALS_LIST
+            for alpha in [0, 1]
         ]
         return params_list
diff --git a/sktime/performance_metrics/forecasting/_functions.py b/sktime/performance_metrics/forecasting/_functions.py
index 624c18f9c..99c7078f1 100644
--- a/sktime/performance_metrics/forecasting/_functions.py
+++ b/sktime/performance_metrics/forecasting/_functions.py
@@ -1495,21 +1495,30 @@ def mean_absolute_percentage_error(
     horizon_weight=None,
     multioutput="uniform_average",
     symmetric=False,
+    relative_to="y_true",
     **kwargs,
 ):
-    """Mean absolute percentage error (MAPE) or symmetric version.
+    """Mean absolute percentage error (MAPE) or symmetric MAPE (sMAPE).
 
-    If ``symmetric`` is False then calculates MAPE and if ``symmetric`` is True
-    then calculates symmetric mean absolute percentage error (sMAPE). Both
-    MAPE and sMAPE output is non-negative floating point. The best value is 0.0.
+    - If ``symmetric`` is False, calculates MAPE (mean absolute percentage error).
+      By default, this is relative to the true values.
+      If ``relative_to="y_pred"``, then MAPE is relative to the predicted values.
+    - If ``symmetric`` is True, calculates sMAPE
+      (symmetric mean absolute percentage error).
+
+    All variants output a non-negative floating point value. The best value is 0.0.
+
+    - MAPE is measured in percentage error relative to the true values.
+    - sMAPE is measured in symmetric percentage error relative to the sum of
+      absolute true and predicted values.
 
     sMAPE is measured in percentage error relative to the test data. Because it
     takes the absolute value rather than square the percentage forecast
     error, it penalizes large errors less than MSPE, RMSPE, MdSPE or RMdSPE.
 
-    There is no limit on how large the error can be, particulalrly when ``y_true``
-    values are close to zero. In such cases the function returns a large value
-    instead of ``inf``.
+    There is no limit on how large the error can be, particularly when
+    denominator values are close to zero. In such cases the function returns
+    a large value instead of ``inf``.
 
     Parameters
     ----------
@@ -1524,8 +1533,8 @@ def mean_absolute_percentage_error(
     horizon_weight : array-like of shape (fh,), default=None
         Forecast horizon weights.
 
-    multioutput : {'raw_values', 'uniform_average'}  or array-like of shape \
-            (n_outputs,), default='uniform_average'
+    multioutput : {'raw_values', 'uniform_average'}  or array-like of shape (n_outputs,),\
+                  default='uniform_average'
         Defines how to aggregate metric for multivariate (multioutput) data.
         If array-like, values used as weights to average the errors.
         If 'raw_values', returns a full set of errors in case of multioutput input.
@@ -1534,14 +1543,21 @@ def mean_absolute_percentage_error(
     symmetric : bool, default=False
         Calculates symmetric version of metric if True.
 
+    relative_to : {'y_true', 'y_pred'}, default='y_true'
+        Determines the denominator of the percentage error.
+
+        * If 'y_true', the denominator is the true values,
+        * If 'y_pred', the denominator is the predicted values.
+
     Returns
     -------
     loss : float
         MAPE or sMAPE loss.
-        If multioutput is 'raw_values', then MAPE or sMAPE is returned for each
-        output separately.
-        If multioutput is 'uniform_average' or an ndarray of weights, then the
-        weighted average MAPE or sMAPE of all output errors is returned.
+
+        * If multioutput is 'raw_values', then MAPE or sMAPE is returned for each
+          output separately.
+        * If multioutput is 'uniform_average' or an ndarray of weights, then the
+          weighted average MAPE or sMAPE of all output errors is returned.
 
     See Also
     --------
@@ -1559,40 +1575,43 @@ def mean_absolute_percentage_error(
     >>> import numpy as np
     >>> from sktime.performance_metrics.forecasting import \
     mean_absolute_percentage_error
+    >>> # Univariate MAPE
     >>> y_true = np.array([3, -0.5, 2, 7, 2])
     >>> y_pred = np.array([2.5, 0.0, 2, 8, 1.25])
     >>> mean_absolute_percentage_error(y_true, y_pred, symmetric=False)  # doctest: +SKIP
     np.float64(0.33690476190476193)
+    >>> # Univariate sMAPE
     >>> mean_absolute_percentage_error(y_true, y_pred, symmetric=True)  # doctest: +SKIP
     np.float64(0.5553379953379953)
+    >>> # Multivariate MAPE
     >>> y_true = np.array([[0.5, 1], [-1, 1], [7, -6]])
     >>> y_pred = np.array([[0, 2], [-1, 2], [8, -5]])
     >>> mean_absolute_percentage_error(y_true, y_pred, symmetric=False)  # doctest: +SKIP
     np.float64(0.5515873015873016)
+    >>> # Multivariate sMAPE
     >>> mean_absolute_percentage_error(y_true, y_pred, symmetric=True)  # doctest: +SKIP
     np.float64(0.6080808080808081)
-    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput='raw_values', \
-        symmetric=False)  # doctest: +SKIP
+    >>> # Raw output per dimension
+    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput='raw_values', symmetric=False)
     array([0.38095238, 0.72222222])
-    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput='raw_values', \
-        symmetric=True)  # doctest: +SKIP
+    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput='raw_values', symmetric=True)
     array([0.71111111, 0.50505051])
-    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput=[0.3, 0.7], \
-    symmetric=False)  # doctest: +SKIP
+    >>> # Weighted multioutput
+    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput=[0.3, 0.7], symmetric=False)  # doctest: +SKIP
     np.float64(0.6198412698412699)
-    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput=[0.3, 0.7], \
-    symmetric=True)  # doctest: +SKIP
+    >>> mean_absolute_percentage_error(y_true, y_pred, multioutput=[0.3, 0.7], symmetric=True)  # doctest: +SKIP
     np.float64(0.5668686868686869)
+    >>> # Univariate MAPE relative to predicted values
+    >>> mean_absolute_percentage_error(y_true, y_pred, relative_to='y_pred')  # doctest: +SKIP
     """  # noqa: E501
     _, y_true, y_pred, multioutput = _check_reg_targets(y_true, y_pred, multioutput)
     if horizon_weight is not None:
         check_consistent_length(y_true, horizon_weight)
 
-    output_errors = np.average(
-        np.abs(_percentage_error(y_true, y_pred, symmetric=symmetric)),
-        weights=horizon_weight,
-        axis=0,
+    perc_err = _percentage_error(
+        y_true, y_pred, symmetric=symmetric, relative_to=relative_to
     )
+    output_errors = np.average(np.abs(perc_err), weights=horizon_weight, axis=0)
 
     if isinstance(multioutput, str):
         if multioutput == "raw_values":
@@ -1611,6 +1630,7 @@ def median_absolute_percentage_error(
     horizon_weight=None,
     multioutput="uniform_average",
     symmetric=False,
+    relative_to="y_true",
     **kwargs,
 ):
     """Median absolute percentage error (MdAPE) or symmetric version.
@@ -1654,6 +1674,12 @@ def median_absolute_percentage_error(
     symmetric : bool, default=False
         Calculates symmetric version of metric if True.
 
+    relative_to : {'y_true', 'y_pred'}, default='y_true'
+        Determines the denominator of the percentage error.
+
+        * If 'y_true', the denominator is the true values,
+        * If 'y_pred', the denominator is the predicted values.
+
     Returns
     -------
     loss : float
@@ -1707,12 +1733,21 @@ def median_absolute_percentage_error(
     _, y_true, y_pred, multioutput = _check_reg_targets(y_true, y_pred, multioutput)
     if horizon_weight is None:
         output_errors = np.median(
-            np.abs(_percentage_error(y_true, y_pred, symmetric=symmetric)), axis=0
+            np.abs(
+                _percentage_error(
+                    y_true, y_pred, symmetric=symmetric, relative_to=relative_to
+                )
+            ),
+            axis=0,
         )
     else:
         check_consistent_length(y_true, horizon_weight)
         output_errors = _weighted_percentile(
-            np.abs(_percentage_error(y_pred, y_true, symmetric=symmetric)),
+            np.abs(
+                _percentage_error(
+                    y_pred, y_true, symmetric=symmetric, relative_to=relative_to
+                )
+            ),
             sample_weight=horizon_weight,
         )
 
@@ -1734,6 +1769,7 @@ def mean_squared_percentage_error(
     multioutput="uniform_average",
     square_root=False,
     symmetric=False,
+    relative_to="y_true",
     **kwargs,
 ):
     """Mean squared percentage error (MSPE) or square root version.
@@ -1781,6 +1817,12 @@ def mean_squared_percentage_error(
     symmetric : bool, default=False
         Calculates symmetric version of metric if True.
 
+    relative_to : {'y_true', 'y_pred'}, default='y_true'
+        Determines the denominator of the percentage error.
+
+        * If 'y_true', the denominator is the true values,
+        * If 'y_pred', the denominator is the predicted values.
+
     Returns
     -------
     loss : float
@@ -1837,7 +1879,11 @@ def mean_squared_percentage_error(
         check_consistent_length(y_true, horizon_weight)
 
     output_errors = np.average(
-        np.square(_percentage_error(y_true, y_pred, symmetric=symmetric)),
+        np.square(
+            _percentage_error(
+                y_true, y_pred, symmetric=symmetric, relative_to=relative_to
+            )
+        ),
         weights=horizon_weight,
         axis=0,
     )
@@ -1863,6 +1909,7 @@ def median_squared_percentage_error(
     multioutput="uniform_average",
     square_root=False,
     symmetric=False,
+    relative_to="y_true",
     **kwargs,
 ):
     """Median squared percentage error (MdSPE)  or square root version.
@@ -1914,6 +1961,12 @@ def median_squared_percentage_error(
     symmetric : bool, default=False
         Calculates symmetric version of metric if True.
 
+    relative_to : {'y_true', 'y_pred'}, default='y_true'
+        Determines the denominator of the percentage error.
+
+        * If 'y_true', the denominator is the true values,
+        * If 'y_pred', the denominator is the predicted values.
+
     Returns
     -------
     loss : float
@@ -1967,7 +2020,9 @@ def median_squared_percentage_error(
     np.float64(0.7428571428571428)
     """  # noqa: E501
     _, y_true, y_pred, multioutput = _check_reg_targets(y_true, y_pred, multioutput)
-    perc_err = _percentage_error(y_true, y_pred, symmetric=symmetric)
+    perc_err = _percentage_error(
+        y_true, y_pred, symmetric=symmetric, relative_to=relative_to
+    )
     if horizon_weight is None:
         output_errors = np.median(np.square(perc_err), axis=0)
     else:
@@ -2738,7 +2793,7 @@ def _relative_error(y_true, y_pred, y_pred_benchmark):
     return (y_true - y_pred) / denominator
 
 
-def _percentage_error(y_true, y_pred, symmetric=False):
+def _percentage_error(y_true, y_pred, symmetric=False, relative_to="y_true"):
     """Percentage error.
 
     Parameters
@@ -2754,6 +2809,9 @@ def _percentage_error(y_true, y_pred, symmetric=False):
     symmetric : bool, default = False
         Whether to calculate symmetric percentage error.
 
+    relative_to : bool, default = "y_true"
+        Whether to calculate percentage error by forecast.
+
     Returns
     -------
     percentage_error : float
@@ -2764,13 +2822,10 @@ def _percentage_error(y_true, y_pred, symmetric=False):
     forecast accuracy", International Journal of Forecasting, Volume 22, Issue 4.
     """
     if symmetric:
-        # Alternatively could use np.abs(y_true + y_pred) in denom
-        # Results will be different if y_true and y_pred have different signs
-        percentage_error = (
-            2
-            * np.abs(y_true - y_pred)
-            / np.maximum(np.abs(y_true) + np.abs(y_pred), EPS)
-        )
+        denominator = np.maximum(np.abs(y_true) + np.abs(y_pred), EPS) / 2
+    elif relative_to == "y_pred":
+        denominator = np.maximum(np.abs(y_pred), EPS)
     else:
-        percentage_error = (y_true - y_pred) / np.maximum(np.abs(y_true), EPS)
+        denominator = np.maximum(np.abs(y_true), EPS)
+    percentage_error = np.abs(y_true - y_pred) / denominator
     return percentage_error
diff --git a/sktime/pipeline/pipeline.py b/sktime/pipeline/pipeline.py
index 6efc1d813..f8851b70a 100644
--- a/sktime/pipeline/pipeline.py
+++ b/sktime/pipeline/pipeline.py
@@ -165,7 +165,7 @@ class Pipeline(BaseEstimator):
             "such as TransformedTargetForecaster, ForecastingPipeline, "
             "ClassificationPipeline, etc., see for instance "
             "notebooks 01_forecasting.ipynb and "
-            "02_classification.ipynb at"
+            "02_classification.ipynb at "
             "https://github.com/sktime/sktime/blob/main/examples/.",
             stacklevel=1,
         )


Test patch:
diff --git a/sktime/forecasting/compose/tests/test_reduce.py b/sktime/forecasting/compose/tests/test_reduce.py
index 5cac6aa91..775494452 100644
--- a/sktime/forecasting/compose/tests/test_reduce.py
+++ b/sktime/forecasting/compose/tests/test_reduce.py
@@ -836,7 +836,7 @@ def test_recursive_reduction_with_X():
     manual_input = np.hstack([last_window, X_new.values])
     manual_pred = manual_lr.predict(manual_input)
 
-    np.allclose(y_pred, manual_pred)
+    assert np.allclose(y_pred, manual_pred)
 
 
 @pytest.mark.skipif(
@@ -867,7 +867,7 @@ def test_pooled_direct_reduction():
     last_win = np.array([[3, 4], [7, 8]])
     lr_pred = lr.predict(last_win)
 
-    return np.allclose(y_pred, lr_pred.reshape(-1, 1), rtol=1e-3)
+    assert np.allclose(y_pred, lr_pred.reshape(-1, 1), rtol=1e-3)
 
 
 @pytest.mark.skipif(
@@ -897,4 +897,4 @@ def test_pooled_recursive_reduction():
     last_win = np.array([[3, 4], [7, 8]])
     lr_pred = lr.predict(last_win)
 
-    return np.allclose(y_pred, lr_pred.reshape(-1, 1), rtol=1e-3)
+    assert np.allclose(y_pred, lr_pred.reshape(-1, 1), rtol=1e-3)
diff --git a/sktime/forecasting/tests/test_reconcile.py b/sktime/forecasting/tests/test_reconcile.py
index 5f3a4219a..361a2c8eb 100644
--- a/sktime/forecasting/tests/test_reconcile.py
+++ b/sktime/forecasting/tests/test_reconcile.py
@@ -5,15 +5,21 @@
 __author__ = ["ciaran-g"]
 
 import numpy as np
+import pandas as pd
 import pytest
 from pandas.testing import assert_frame_equal
 
 from sktime.forecasting.base import ForecastingHorizon
 from sktime.forecasting.exp_smoothing import ExponentialSmoothing
+from sktime.forecasting.naive import NaiveForecaster
 from sktime.forecasting.reconcile import ReconcilerForecaster
 from sktime.tests.test_switch import run_test_for_class
 from sktime.transformations.hierarchical.aggregate import Aggregator
-from sktime.utils._testing.hierarchical import _bottom_hier_datagen, _make_hierarchical
+from sktime.utils._testing.hierarchical import (
+    _bottom_hier_datagen,
+    _make_hierarchical,
+    _make_index,
+)
 from sktime.utils.dependencies import _check_soft_dependencies
 
 # get all the methods
@@ -173,3 +179,37 @@ def test_reconcilerforecaster_return_totals(method, return_totals):
     else:
         assert len(y_test) == len(y_pred)
         assert y_test.index.equals(y_pred.index)
+
+
+@pytest.mark.parametrize("alpha", [0, 1])
+def test_reconcilerforecaster_singular(alpha):
+    """
+    Test that ReconcilerForecaster handles highly correlated series,
+    where alpha=0 leads to a singular residual covariance matrix.
+    """
+    hierarchy_levels = (1, 2)
+    levels = [
+        [f"h{i}_{j}" for j in range(hierarchy_levels[i])]
+        for i in range(len(hierarchy_levels))
+    ]
+    level_names = [f"h{i}" for i in range(len(hierarchy_levels))]
+
+    time_index = _make_index(12, None)
+    index = pd.MultiIndex.from_product(
+        levels + [time_index], names=level_names + ["time"]
+    )
+    input_data = np.array([i for i in range(24)])
+    input_data = np.stack([input_data, input_data], 1)
+    n_columns = 2
+
+    df = pd.DataFrame(
+        data=input_data, index=index, columns=[f"c{i}" for i in range(n_columns)]
+    )
+    base_forecaster = NaiveForecaster(strategy="mean")
+
+    forecaster = ReconcilerForecaster(
+        forecaster=base_forecaster, method="mint_cov", alpha=alpha
+    )
+    if alpha == 0:
+        with pytest.raises(np.linalg.LinAlgError, match="Singular matrix"):
+            forecaster.fit(df, fh=[1]).predict()