From 9bb4ebcb6b5fc00471dbc0fafc1092129f093927 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Fri, 11 Aug 2023 16:59:02 +0200 Subject: [PATCH 1/6] [Hierarchical] Allow computing projection based on a non-Euclidean norm --- .../model/deepvar_hierarchical/_estimator.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py index 8532f30eba..50feec0a3e 100755 --- a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py +++ b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py @@ -73,7 +73,10 @@ def constraint_mat(S: np.ndarray) -> np.ndarray: return A -def null_space_projection_mat(A: np.ndarray) -> np.ndarray: +def null_space_projection_mat( + A: np.ndarray, + D: Optional[np.ndarray], +) -> np.ndarray: """ Computes the projection matrix for projecting onto the null space of A. @@ -82,6 +85,12 @@ def null_space_projection_mat(A: np.ndarray) -> np.ndarray: A The constraint matrix A in the equation: Ay = 0 (y being the values/forecasts of all time series in the hierarchy). + D + Positive definite matrix (typically a diagonal matrix). Optional. + If provided then the distance between the reconciled and unreconciled + forecasts is calculated based on the norm induced by D. Useful for + weighing the distances differently for each level of the hierarchy. + By default Euclidean distance is used. Returns ------- @@ -89,8 +98,11 @@ def null_space_projection_mat(A: np.ndarray) -> np.ndarray: Projection matrix, shape (total_num_time_series, total_num_time_series) """ num_ts = A.shape[1] - return np.eye(num_ts) - A.T @ np.linalg.pinv(A @ A.T) @ A - + if D is None: + return np.eye(num_ts) - A.T @ np.linalg.pinv(A @ A.T) @ A + else: + D_inv = np.linalg.inv(D) + return np.eye(num_ts) - D_inv @ A.T @ np.linalg.pinv(A @ D_inv @ A.T) @ A class DeepVARHierarchicalEstimator(DeepVAREstimator): """ From f3b17a00aebbfa5dc2bcbf51dfa3977e54bd3c36 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Fri, 11 Aug 2023 17:21:07 +0200 Subject: [PATCH 2/6] Fix tests --- .../mx/model/deepvar_hierarchical/_estimator.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py index 50feec0a3e..d6bf49ef89 100755 --- a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py +++ b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py @@ -75,7 +75,7 @@ def constraint_mat(S: np.ndarray) -> np.ndarray: def null_space_projection_mat( A: np.ndarray, - D: Optional[np.ndarray], + D: Optional[np.ndarray] = None, ) -> np.ndarray: """ Computes the projection matrix for projecting onto the null space of A. @@ -88,9 +88,9 @@ def null_space_projection_mat( D Positive definite matrix (typically a diagonal matrix). Optional. If provided then the distance between the reconciled and unreconciled - forecasts is calculated based on the norm induced by D. Useful for - weighing the distances differently for each level of the hierarchy. - By default Euclidean distance is used. + forecasts is calculated based on the norm induced by D. Useful for + weighing the distances differently for each level of the hierarchy. + By default Euclidean distance is used. Returns ------- @@ -102,7 +102,10 @@ def null_space_projection_mat( return np.eye(num_ts) - A.T @ np.linalg.pinv(A @ A.T) @ A else: D_inv = np.linalg.inv(D) - return np.eye(num_ts) - D_inv @ A.T @ np.linalg.pinv(A @ D_inv @ A.T) @ A + return ( + np.eye(num_ts) - D_inv @ A.T @ np.linalg.pinv(A @ D_inv @ A.T) @ A + ) + class DeepVARHierarchicalEstimator(DeepVAREstimator): """ From 6062a5774eee0621e5b69ebaccae5fdfb09337f3 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Fri, 11 Aug 2023 17:39:20 +0200 Subject: [PATCH 3/6] Address comments --- src/gluonts/mx/model/deepvar_hierarchical/_estimator.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py index d6bf49ef89..14b7fbcdad 100755 --- a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py +++ b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py @@ -124,6 +124,12 @@ class DeepVARHierarchicalEstimator(DeepVAREstimator): Length of the prediction horizon S Summation or aggregation matrix. + D + Positive definite matrix (typically a diagonal matrix). Optional. + If provided then the distance between the reconciled and unreconciled + forecasts is calculated based on the norm induced by D. Useful for + weighing the distances differently for each level of the hierarchy. + By default Euclidean distance is used. num_samples_for_loss Number of samples to draw from the predicted distribution to compute the training loss. @@ -210,6 +216,7 @@ def __init__( freq: str, prediction_length: int, S: np.ndarray, + D: Optional[np.ndarray] = None, num_samples_for_loss: int = 200, likelihood_weight: float = 0.0, CRPS_weight: float = 1.0, @@ -288,7 +295,7 @@ def __init__( ), "Cannot project only during training (and not during prediction)" A = constraint_mat(S.astype(self.dtype)) - M = null_space_projection_mat(A) + M = null_space_projection_mat(A=A, D=D) ctx = self.trainer.ctx self.M = mx.nd.array(M, ctx=ctx) self.A = mx.nd.array(A, ctx=ctx) From 5c03a335c0a9925911b466d508d531b02b4e8199 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Mon, 14 Aug 2023 17:07:11 +0200 Subject: [PATCH 4/6] Add tests --- .../model/deepvar_hierarchical/_estimator.py | 3 +++ .../test_reconcile_samples.py | 27 ++++++++++++++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py index 14b7fbcdad..ea4e93892e 100755 --- a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py +++ b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py @@ -101,6 +101,9 @@ def null_space_projection_mat( if D is None: return np.eye(num_ts) - A.T @ np.linalg.pinv(A @ A.T) @ A else: + assert np.all( + np.linalg.eigvals(D) > 0 + ), "`D` must be positive definite." D_inv = np.linalg.inv(D) return ( np.eye(num_ts) - D_inv @ A.T @ np.linalg.pinv(A @ D_inv @ A.T) @ A diff --git a/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py b/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py index 8e67a74931..9c874c14ee 100644 --- a/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py +++ b/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py @@ -38,7 +38,6 @@ num_bottom_ts = S.shape[1] A = constraint_mat(S) -reconciliation_mat = null_space_projection_mat(A) @pytest.mark.parametrize( @@ -53,6 +52,28 @@ 100.0 + 2.0 * np.random.standard_normal(size=(10, 32, S.shape[0])), ], ) +@pytest.mark.parametrize( + "D", + [ + None, + # Root gets the maximum weight and the two aggregated levels get + # more weight than the leaf level. + np.diag( + [4, 2, 2, 1, 1, 1, 1] + ), + # Random diagonal matrix + np.diag( + np.random.rand(S.shape[0]) + ), + # Random positive definite matrix + np.diag( + np.random.rand(S.shape[0]) + ) + np.dot( + np.array([[4, 2, 2, 1, 1, 1, 1]]).T, + np.array([[4, 2, 2, 1, 1, 1, 1]]) + ) + ], +) @pytest.mark.parametrize( "seq_axis", [ @@ -63,9 +84,9 @@ [1, 0], ], ) -def test_reconciliation_error(samples, seq_axis): +def test_reconciliation_error(samples, D, seq_axis): coherent_samples = reconcile_samples( - reconciliation_mat=mx.nd.array(reconciliation_mat), + reconciliation_mat=mx.nd.array(null_space_projection_mat(A=A, D=D)), samples=mx.nd.array(samples), seq_axis=seq_axis, ) From f072bb323c8fc81c7f610928a0148e43a3218bb2 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Mon, 14 Aug 2023 17:10:20 +0200 Subject: [PATCH 5/6] Fix black --- .../test_reconcile_samples.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py b/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py index 9c874c14ee..f80ec0f0a8 100644 --- a/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py +++ b/test/mx/model/deepvar_hierarchical/test_reconcile_samples.py @@ -58,20 +58,15 @@ None, # Root gets the maximum weight and the two aggregated levels get # more weight than the leaf level. - np.diag( - [4, 2, 2, 1, 1, 1, 1] - ), + np.diag([4, 2, 2, 1, 1, 1, 1]), # Random diagonal matrix - np.diag( - np.random.rand(S.shape[0]) - ), + np.diag(np.random.rand(S.shape[0])), # Random positive definite matrix - np.diag( - np.random.rand(S.shape[0]) - ) + np.dot( + np.diag(np.random.rand(S.shape[0])) + + np.dot( np.array([[4, 2, 2, 1, 1, 1, 1]]).T, - np.array([[4, 2, 2, 1, 1, 1, 1]]) - ) + np.array([[4, 2, 2, 1, 1, 1, 1]]), + ), ], ) @pytest.mark.parametrize( From eb75ddec4fbbb99dea1322c68f2932927da5d436 Mon Sep 17 00:00:00 2001 From: Syama Sundar Rangapuram Date: Tue, 15 Aug 2023 14:31:29 +0200 Subject: [PATCH 6/6] Add symmetry check --- src/gluonts/mx/model/deepvar_hierarchical/_estimator.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py index ea4e93892e..7dd70587bb 100755 --- a/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py +++ b/src/gluonts/mx/model/deepvar_hierarchical/_estimator.py @@ -86,7 +86,8 @@ def null_space_projection_mat( The constraint matrix A in the equation: Ay = 0 (y being the values/forecasts of all time series in the hierarchy). D - Positive definite matrix (typically a diagonal matrix). Optional. + Symmetric positive definite matrix (typically a diagonal matrix). + Optional. If provided then the distance between the reconciled and unreconciled forecasts is calculated based on the norm induced by D. Useful for weighing the distances differently for each level of the hierarchy. @@ -101,9 +102,11 @@ def null_space_projection_mat( if D is None: return np.eye(num_ts) - A.T @ np.linalg.pinv(A @ A.T) @ A else: + assert np.all(D == D.T), "`D` must be symmetric." assert np.all( np.linalg.eigvals(D) > 0 ), "`D` must be positive definite." + D_inv = np.linalg.inv(D) return ( np.eye(num_ts) - D_inv @ A.T @ np.linalg.pinv(A @ D_inv @ A.T) @ A