Repository: sktime/sktime
Pull request: [ENH] Add regularization hyperparameter in ReconcilerForecaster
<!--
Welcome to sktime, and thanks for contributing!
Please have a look at our contribution guide:
https://www.sktime.net/en/latest/get_involved/contributing.html
-->

#### Reference Issues/PRs
Fixes #7051 
<!--
Example: Fixes #1234. See also #3456.

Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests.
If no issue exists, you can open one here: https://github.com/sktime/sktime/issues
-->


#### What does this implement/fix? Explain your changes.
Add regularization parameter to avoid singular covariance matrix in `ReconcilerForecaster`
<!--
A clear and concise description of what you have implemented.
-->

#### What should a reviewer concentrate their feedback on?

<!-- This section is particularly useful if you have a pull request that is still in development. You can guide the reviews to focus on the parts that are ready for their comments. We suggest using bullets (indicated by * or -) and filled checkboxes [x] here -->

#### Did you add any tests for the change?
Yes
<!-- This section is useful if you have added a test in addition to the existing ones. This will ensure that further changes to these files won't introduce the same kind of bug. It is considered good practice to add tests with newly added code to enforce the fact that the code actually works. This will reduce the chance of introducing logical bugs.
-->

#### PR checklist
<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/sktime/sktime/blob/main/CONTRIBUTORS.md) with any new badges I've earned :-)
  How to: add yourself to the [all-contributors file](https://github.com/sktime/sktime/blob/main/.all-contributorsrc) in the `sktime` root directory (not the `CONTRIBUTORS.md`). Common badges: `code` - fixing a bug, or adding code logic. `doc` - writing or improving documentation or docstrings. `bug` - reporting or diagnosing a bug (get this plus `code` if you also fixed the bug in the PR).`maintenance` - CI, test framework, release.
  See here for [full badge reference](https://allcontributors.org/docs/en/emoji-key)
- [x] The PR title starts with either [ENH], [MNT], [DOC], or [BUG]. [BUG] - bugfix, [MNT] - CI, test framework, [ENH] - adding or improving code, [DOC] - writing or improving documentation or docstrings.



Problem statement:
[ENH] Add optional regularization hyperparameter to ReconcilerForecaster
**Is your feature request related to a problem? Please describe.**
When using optimal reconciliation methods, if more than 1 timeseries in the hierarchy has constant 0 values or is highly correlated with another, the residual covariance matrix can be singular and its inversion will break the code or produce unstable forecasts. If we add a possible regularization parameter, this issue can be easily solved.

**Describe the solution you'd like**
Adding a float hyperparameter `lambda` that would be added to the diagonal of the covariance matrix before inversion.

Hints from issue discussion:
FYI @ciaran-g 
I could work on this

Commit message: Add tests

Patch to review:
diff --git a/.all-contributorsrc b/.all-contributorsrc
index ad4440829e1..843d0aeebe7 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -3240,6 +3240,15 @@
         "code"
       ]
     },
+    {
+      "login": "HarshvirSandhu",
+      "name": "Harshvir Sandhu",
+      "avatar_url": "https://avatars.githubusercontent.com/HarshvirSandhu",
+      "profile": "https://github.com/HarshvirSandhu",
+      "contributions": [
+        "code"
+      ]
+    },
     {
       "login": "alyssadsouza",
       "name": "Alyssa D'Souza",
@@ -3261,7 +3270,7 @@
     },
     {
       "login": "TanviPooranmal",
-      "name": "Tanvi Pooranmal Meena"
+      "name": "Tanvi Pooranmal Meena",
       "avatar_url": "https://avatars.githubusercontent.com/TanviPooranmal",
       "profile": "https://github.com/TanviPooranmal",
       "contributions": [
diff --git a/sktime/forecasting/reconcile.py b/sktime/forecasting/reconcile.py
index 050ecc397b3..1d85c61c06a 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
+        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,12 @@ def _get_g_matrix_mint(self, shrink=False, diag_only=False):
             np.dot(
                 inv(
                     np.dot(np.transpose(self.s_matrix), np.dot(cov_mat, self.s_matrix))
+                    + self.alpha
+                    * np.eye(
+                        np.dot(
+                            np.transpose(self.s_matrix), np.dot(cov_mat, self.s_matrix)
+                        ).shape[0]
+                    )
                 ),
                 np.dot(np.transpose(self.s_matrix), cov_mat),
             )
diff --git a/sktime/forecasting/tests/test_reconcile.py b/sktime/forecasting/tests/test_reconcile.py
index 5f3a4219a45..e3224a1265a 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,33 @@ 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):
+    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()


Code review comments:
diff_hunk: @@ -402,6 +406,12 @@ def _get_g_matrix_mint(self, shrink=False, diag_only=False):
             np.dot(
                 inv(
                     np.dot(np.transpose(self.s_matrix), np.dot(cov_mat, self.s_matrix))
+                    + self.alpha
+                    * np.eye(
+                        np.dot(
path: sktime/forecasting/reconcile.py
text: @user1:
it is not optimal to do the computation in the `dot` just in order to get the size of hte matrix. The size of the matrix should be known from the objects already floating around - could you replace the "shape" with something that does not carry out unnecessary computations?

@author:
I think it can be replaced with `self.alpha * np.eye(self.s_matrix.shape[1])`

@user1:
yes, exactly