From 8cd41c070e71851e34f3c991961318a6da62a348 Mon Sep 17 00:00:00 2001 From: Christopher Titchen <109701765+christophertitchen@users.noreply.github.com> Date: Tue, 3 Sep 2024 16:39:09 +0100 Subject: [PATCH 1/3] ENH: :sparkles: Add sparse non-negative reconciliation heuristic. --- hierarchicalforecast/_modidx.py | 2 + hierarchicalforecast/methods.py | 65 ++++++++++++++++++++++++++++++++ nbs/methods.ipynb | 67 ++++++++++++++++++++++++++++++++- 3 files changed, 133 insertions(+), 1 deletion(-) diff --git a/hierarchicalforecast/_modidx.py b/hierarchicalforecast/_modidx.py index 237d1ac..f8f8f64 100644 --- a/hierarchicalforecast/_modidx.py +++ b/hierarchicalforecast/_modidx.py @@ -102,6 +102,8 @@ 'hierarchicalforecast/methods.py'), 'hierarchicalforecast.methods.MinTraceSparse._get_PW_matrices': ( 'methods.html#mintracesparse._get_pw_matrices', 'hierarchicalforecast/methods.py'), + 'hierarchicalforecast.methods.MinTraceSparse.fit': ( 'methods.html#mintracesparse.fit', + 'hierarchicalforecast/methods.py'), 'hierarchicalforecast.methods.OptimalCombination': ( 'methods.html#optimalcombination', 'hierarchicalforecast/methods.py'), 'hierarchicalforecast.methods.OptimalCombination.__init__': ( 'methods.html#optimalcombination.__init__', diff --git a/hierarchicalforecast/methods.py b/hierarchicalforecast/methods.py index 35eb049..4a34156 100644 --- a/hierarchicalforecast/methods.py +++ b/hierarchicalforecast/methods.py @@ -1098,6 +1098,71 @@ def get_P_action(y): return P, W + def fit(self, + S: sparse.csr_matrix, + y_hat: np.ndarray, + y_insample: Optional[np.ndarray] = None, + y_hat_insample: Optional[np.ndarray] = None, + sigmah: Optional[np.ndarray] = None, + intervals_method: Optional[str] = None, + num_samples: Optional[int] = None, + seed: Optional[int] = None, + tags: Dict[str, np.ndarray] = None, + idx_bottom: Optional[np.ndarray] = None): + # Clip the base forecasts if required to align them with their use in practice. + if self.nonnegative: + self.y_hat = np.clip(y_hat, 0) + else: + self.y_hat = y_hat + # Get the reconciliation matrices. + self.P, self.W = self._get_PW_matrices( + S=S, + y_hat=self.y_hat, + y_insample=y_insample, + y_hat_insample=y_hat_insample, + idx_bottom=idx_bottom, + ) + + if self.nonnegative: + # Get the number of leaf nodes. + _, n_bottom = S.shape + # Although it is now sufficient to ensure that all of the entries in P are + # positive, as it is implemented as a linear operator for the iterative + # method to solve the sparse linear system, we need to reconcile to find + # if any of the coherent bottom level point forecasts are negative. + y_tilde = self._reconcile( + S=S, P=self.P, y_hat=self.y_hat, level=None, sampler=None + )["mean"][-n_bottom:] + # Find if any of the forecasts are negative. + if np.any(y_tilde < 0): + # Clip the negative forecasts. + y_tilde = np.clip(y_tilde, 0) + # Force non-negative coherence by overwriting the base forecasts with + # the aggregated, clipped bottom level forecasts. + self.y_hat = S @ y_tilde + # Overwrite the attributes for the P and W matrices with those for + # bottom-up reconciliation to force projection onto the non-negative + # coherent subspace. + self.P, self.W = BottomUpSparse()._get_PW_matrices(S=S, idx_bottom=None) + + # Get the sampler for probabilistic reconciliation. + self.sampler = self._get_sampler( + S=S, + P=self.P, + W=self.W, + y_hat=self.y_hat, + y_insample=y_insample, + y_hat_insample=y_hat_insample, + sigmah=sigmah, + intervals_method=intervals_method, + num_samples=num_samples, + seed=seed, + tags=tags, + ) + # Set the instance as fitted. + self.fitted = True + return self + # %% ../nbs/methods.ipynb 55 class OptimalCombination(MinTrace): """Optimal Combination Reconciliation Class. diff --git a/nbs/methods.ipynb b/nbs/methods.ipynb index c74c0ec..e23fd4d 100644 --- a/nbs/methods.ipynb +++ b/nbs/methods.ipynb @@ -1708,7 +1708,72 @@ " )\n", " W = sparse.spdiags(W_diag, 0, W_diag.size, W_diag.size)\n", "\n", - " return P, W" + " return P, W\n", + "\n", + " def fit(self,\n", + " S: sparse.csr_matrix,\n", + " y_hat: np.ndarray,\n", + " y_insample: Optional[np.ndarray] = None,\n", + " y_hat_insample: Optional[np.ndarray] = None,\n", + " sigmah: Optional[np.ndarray] = None,\n", + " intervals_method: Optional[str] = None,\n", + " num_samples: Optional[int] = None,\n", + " seed: Optional[int] = None, \n", + " tags: Dict[str, np.ndarray] = None,\n", + " idx_bottom: Optional[np.ndarray] = None):\n", + " # Clip the base forecasts if required to align them with their use in practice.\n", + " if self.nonnegative:\n", + " self.y_hat = np.clip(y_hat, 0)\n", + " else:\n", + " self.y_hat = y_hat\n", + " # Get the reconciliation matrices.\n", + " self.P, self.W = self._get_PW_matrices(\n", + " S=S, \n", + " y_hat=self.y_hat, \n", + " y_insample=y_insample, \n", + " y_hat_insample=y_hat_insample, \n", + " idx_bottom=idx_bottom,\n", + " )\n", + "\n", + " if self.nonnegative:\n", + " # Get the number of leaf nodes.\n", + " _, n_bottom = S.shape\n", + " # Although it is now sufficient to ensure that all of the entries in P are \n", + " # positive, as it is implemented as a linear operator for the iterative \n", + " # method to solve the sparse linear system, we need to reconcile to find \n", + " # if any of the coherent bottom level point forecasts are negative.\n", + " y_tilde = self._reconcile(\n", + " S=S, P=self.P, y_hat=self.y_hat, level=None, sampler=None\n", + " )[\"mean\"][-n_bottom:]\n", + " # Find if any of the forecasts are negative.\n", + " if np.any(y_tilde < 0):\n", + " # Clip the negative forecasts.\n", + " y_tilde = np.clip(y_tilde, 0)\n", + " # Force non-negative coherence by overwriting the base forecasts with \n", + " # the aggregated, clipped bottom level forecasts.\n", + " self.y_hat = S @ y_tilde\n", + " # Overwrite the attributes for the P and W matrices with those for \n", + " # bottom-up reconciliation to force projection onto the non-negative \n", + " # coherent subspace.\n", + " self.P, self.W = BottomUpSparse()._get_PW_matrices(S=S, idx_bottom=None) \n", + "\n", + " # Get the sampler for probabilistic reconciliation.\n", + " self.sampler = self._get_sampler(\n", + " S=S,\n", + " P=self.P,\n", + " W=self.W,\n", + " y_hat=self.y_hat,\n", + " y_insample=y_insample,\n", + " y_hat_insample=y_hat_insample,\n", + " sigmah=sigmah,\n", + " intervals_method=intervals_method,\n", + " num_samples=num_samples,\n", + " seed=seed,\n", + " tags=tags,\n", + " )\n", + " # Set the instance as fitted.\n", + " self.fitted = True\n", + " return self" ] }, { From 59a2b6bcf73c8a8cd687e9aa4372d60ae613afc9 Mon Sep 17 00:00:00 2001 From: Christopher Titchen <109701765+christophertitchen@users.noreply.github.com> Date: Tue, 3 Sep 2024 17:22:43 +0100 Subject: [PATCH 2/3] BUG: :bug: Remove legacy absolute tolerance for BiCGSTAB convergence. --- hierarchicalforecast/methods.py | 2 +- nbs/methods.ipynb | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/hierarchicalforecast/methods.py b/hierarchicalforecast/methods.py index 4a34156..f760d57 100644 --- a/hierarchicalforecast/methods.py +++ b/hierarchicalforecast/methods.py @@ -1087,7 +1087,7 @@ def get_P_action(y): (b.size, b.size), matvec=lambda v: R @ (S @ v) ) - x_tilde, exit_code = sparse.linalg.bicgstab(A, b, atol="legacy") + x_tilde, exit_code = sparse.linalg.bicgstab(A, b) return x_tilde diff --git a/nbs/methods.ipynb b/nbs/methods.ipynb index e23fd4d..f7c07b0 100644 --- a/nbs/methods.ipynb +++ b/nbs/methods.ipynb @@ -1699,7 +1699,7 @@ " (b.size, b.size), matvec=lambda v: R @ (S @ v)\n", " )\n", "\n", - " x_tilde, exit_code = sparse.linalg.bicgstab(A, b, atol=\"legacy\")\n", + " x_tilde, exit_code = sparse.linalg.bicgstab(A, b)\n", "\n", " return x_tilde\n", "\n", From 73355cecfe31007101dde9a207e2963680936807 Mon Sep 17 00:00:00 2001 From: Christopher Titchen <109701765+christophertitchen@users.noreply.github.com> Date: Tue, 3 Sep 2024 17:38:45 +0100 Subject: [PATCH 3/3] TST: :white_check_mark: Add tests for sparse MinT. --- hierarchicalforecast/methods.py | 15 +++++---------- nbs/methods.ipynb | 31 ++++++++++++++++++++++++------- 2 files changed, 29 insertions(+), 17 deletions(-) diff --git a/hierarchicalforecast/methods.py b/hierarchicalforecast/methods.py index f760d57..a3d07fc 100644 --- a/hierarchicalforecast/methods.py +++ b/hierarchicalforecast/methods.py @@ -1024,11 +1024,6 @@ def _get_PW_matrices( "Only the methods with diagonal W are supported as sparse operations" ) - if self.nonnegative: - raise NotImplementedError( - "Non-negative MinT is currently not implemented as sparse" - ) - S = sparse.csr_matrix(S) if self.method in res_methods and y_insample is None and y_hat_insample is None: @@ -1111,7 +1106,7 @@ def fit(self, idx_bottom: Optional[np.ndarray] = None): # Clip the base forecasts if required to align them with their use in practice. if self.nonnegative: - self.y_hat = np.clip(y_hat, 0) + self.y_hat = np.clip(y_hat, 0, None) else: self.y_hat = y_hat # Get the reconciliation matrices. @@ -1136,7 +1131,7 @@ def fit(self, # Find if any of the forecasts are negative. if np.any(y_tilde < 0): # Clip the negative forecasts. - y_tilde = np.clip(y_tilde, 0) + y_tilde = np.clip(y_tilde, 0, None) # Force non-negative coherence by overwriting the base forecasts with # the aggregated, clipped bottom level forecasts. self.y_hat = S @ y_tilde @@ -1163,7 +1158,7 @@ def fit(self, self.fitted = True return self -# %% ../nbs/methods.ipynb 55 +# %% ../nbs/methods.ipynb 56 class OptimalCombination(MinTrace): """Optimal Combination Reconciliation Class. @@ -1197,7 +1192,7 @@ def __init__(self, super().__init__(method=method, nonnegative=nonnegative, num_threads=num_threads) self.insample = False -# %% ../nbs/methods.ipynb 64 +# %% ../nbs/methods.ipynb 65 @njit def lasso(X: np.ndarray, y: np.ndarray, lambda_reg: float, max_iters: int = 1_000, @@ -1229,7 +1224,7 @@ def lasso(X: np.ndarray, y: np.ndarray, #print(it) return beta -# %% ../nbs/methods.ipynb 65 +# %% ../nbs/methods.ipynb 66 class ERM(HReconciler): """Optimal Combination Reconciliation Class. diff --git a/nbs/methods.ipynb b/nbs/methods.ipynb index f7c07b0..a51b5bd 100644 --- a/nbs/methods.ipynb +++ b/nbs/methods.ipynb @@ -1636,11 +1636,6 @@ " \"Only the methods with diagonal W are supported as sparse operations\"\n", " )\n", "\n", - " if self.nonnegative:\n", - " raise NotImplementedError(\n", - " \"Non-negative MinT is currently not implemented as sparse\"\n", - " )\n", - "\n", " S = sparse.csr_matrix(S)\n", "\n", " if self.method in res_methods and y_insample is None and y_hat_insample is None:\n", @@ -1723,7 +1718,7 @@ " idx_bottom: Optional[np.ndarray] = None):\n", " # Clip the base forecasts if required to align them with their use in practice.\n", " if self.nonnegative:\n", - " self.y_hat = np.clip(y_hat, 0)\n", + " self.y_hat = np.clip(y_hat, 0, None)\n", " else:\n", " self.y_hat = y_hat\n", " # Get the reconciliation matrices.\n", @@ -1748,7 +1743,7 @@ " # Find if any of the forecasts are negative.\n", " if np.any(y_tilde < 0):\n", " # Clip the negative forecasts.\n", - " y_tilde = np.clip(y_tilde, 0)\n", + " y_tilde = np.clip(y_tilde, 0, None)\n", " # Force non-negative coherence by overwriting the base forecasts with \n", " # the aggregated, clipped bottom level forecasts.\n", " self.y_hat = S @ y_tilde\n", @@ -1917,6 +1912,28 @@ " )" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "for method in [\"ols\", \"wls_struct\"]:\n", + " for nonnegative in [False, True]:\n", + " cls_min_trace = MinTraceSparse(method=method, nonnegative=nonnegative)\n", + " test_close(\n", + " cls_min_trace(\n", + " S=S,\n", + " y_hat=S @ y_hat_bottom,\n", + " y_insample=S @ y_bottom,\n", + " y_hat_insample=S @ y_hat_bottom_insample,\n", + " idx_bottom=idx_bottom if nonnegative else None,\n", + " )[\"mean\"],\n", + " S @ y_hat_bottom,\n", + " )" + ] + }, { "cell_type": "markdown", "metadata": {},