Hierarchically Regularized Deep Forecasting
Authors: Biswajit Paria
Rajat Sen
Amr Ahmed
Abhimanyu Das
Venue: ICLR 2022 Submitted (2022)
Keywords: hierarchical time series
deep learning

Abstract: Hierarchical forecasting is a key problem in many practical multivariate forecasting applications - the goal is to simultaneously predict a large number of correlated time series that are arranged in a pre-specified aggregation hierarchy. The main challenge is to exploit the hierarchical correlations to simultaneously obtain good prediction accuracy for time series at different levels of the hierarchy. In this paper, we propose a new approach for hierarchical forecasting which consists of two components. First, decomposing the time series along a global set of basis time series and modeling hierarchical constraints using the coefficients of the basis decomposition. And second, using a linear autoregressive model with coefficients that vary with time. Unlike past methods, our approach is scalable (inference for a specific time series only needs access to its own history) while also modeling the hierarchical structure via (approximate) coherence constraints among the time series forecasts. We experiment on several public datasets and demonstrate significantly improved overall performance on forecasts at different levels of the hierarchy, compared to existing state-of-the-art hierarchical models.

Under review as a conference paper at ICLR 2022
HIERARCHICALLY REGULARIZED DEEP FORECASTING
Anonymous authors Paper under double-blind review
ABSTRACT
Hierarchical forecasting is key problem in many practical multivariate forecasting applications the goal is to simultancously predict a large number of correlated time series that are arranged in a pre-specified aggregation hierarchy. The main challenge is to exploit the hierarchical correlations to simultaneously obtain good prediction accuracy for time series at different levels of the hierarchy. this paper, we propose new approach for hierarchical forecasting which consists of two components. First, decomposing the time series along a global set of basis time series and modeling hierarchical constraints using the coefficients of the basis decomposition. And second, using a linear autoregressive model with coefficients that vary with time. Unlike past methods, our approach is scalable (inference for a specific time series only needs access to its own history) while also modeling the hierarchical structure via (approximate) coherence constraints among the time series forecasts. We experiment on several public datasets and demonstrate significantly improved overall performance on forecasts different levels of the hierarchy, compared to existing state-of-the hicrarchical models.
INTRODUCTION
Multivariate time series forecasting is a key problem in many domains such as retail demand forecasting (Böse et al., 2017), financial predictions (Zhou et al., 2020), power grid optimization (Hyndman & Fan, 2009), road traffic modeling (Li et al., 2017), and online ads optimization (Cui al., 2011), 1). In many of these setting, the problem involves simultaneously forecasting a large number of possibly correlated time series for various downstream applications. In the retail domain, the time series may capture sales of items in a product inventory, and in power grids, the time series may correspond to energy consumption in household. Often, these time series are arranged in natural multi-level hierarchy- for example in retail forecasting, items are grouped into subcategories and categories, and arranged in product taxonomy. In the case of power consumption forecasting, individual households are grouped into neighborhoods, counties, and cities. The hierarchical structure among the time series can usually be represented as tree, with the leaf nodes corresponding to time series at the finest granularity, and the edges representing parent-child relationships. Figure I illustrates typical hierarchy in the retail forecasting domain for time series of product sales.
In such settings, it is often required to obtain good forecasts, not just for the leaf level time-series (fine grained forecasts), but also for the aggregated time-series corresponding to higher level nodes (coarse gained forecasts). Furthermore, for inter- pretability and business decision making purposes, it is often desirable to obtain predictions that are roughly coherent or con- sistent (Hyndman et al., 2011) with respect to the hierarchy tree. This means that the predictions for each parent time-series is equal to the sum of the predictions for its children time-series. More importantly, incorporating coherence constraints in a hi- erarchical forecasting model captures the natural inductive bias in most hierarchical datasets, where the ground truth parent and children time series indeed adhere to additive constraints. For example, total sales of a product category is equal to the sum of sales of all items in that category.
Root
Categary 1..Catepery
Categary 1... Categary
Nem
Rem-
Ro N_
Figure 1: An example hierarchy for retail demand forecasting. The blue triangle represents the sub- tree rooted the node Storel with leaves denoted by Item i.

---

Under review as a conference paper at ICLR 2022
Some standard approaches for hierarchical forecasting include bottom-up aggregation, or reconciliation based approaches. Bottom-Up Up aggregation involves training a model to obtain pre- dictions for the leaf nodes, and then aggregate up along the hierarchy tree to obtain predictions for higher-level nodes. Reconciliation methods (Ben Taieh & Koo, 2019; Taieb et al., 2017; Van Erven & Cugliari, 2015; Hyndman al., 2016; Wickramasuriya al., 2015; 2020; Panagiotelis et al., 2020) make use of a trained model to obtain predictions for all nodes of the tree, and then, in a separate obtain
post-processing phase, reconcile (or modify) them using various optimization formulations coherent predictions. Both of these approaches suffer from shorteomings in term of either aggregating noise as one moves to higher level predictions (bottom-up up aggregation), or not jointly optimizing the forecasting predictions along with the coherence constraints (for instance, reconciliation).
At the same time, there have been several recent advances on using Deep Neural Network models for multivariate forecasting, including Recurrent Neural Network (RNN), Convolutional Neural Network (CNN) architectures (Salinas et al., 2020; Oreshkin al., 2019; Rangapuram et al., 2018; Benidis et al., 2020), and models designed for multivariate time series based on dimensionality reduction techniques (Sen et al., 2019; Wang et al., 2019; Nguyen & Quanz, 2021; Salinas et al., 2019; Rasul et al., 2020; de Bézenac et al., 2020), that have been shown outperform classical time-series models such as autoregressive and exponential smoothing models (McKenzie, 1984; Hyndman et al., 2008; Hyndman & Athanasopoulos, 20, especially for large datasets. However, most of these approaches do not explicitly address the question of how to model the hierarchical relationships in the dataset. Deep forecasting models based on Graph Neural Networks (GNN) (Bai al., 2020; Cao et al., 2020; Yu et al., 2017; Li et al., 2017; Wu et al., 2020) do offer general framework for learning on graph-structured data. However it is well known (Bojchevski et al., 2020) that GNNs are hard to scale for learning on graphs with a very large number of nodes, which in real-world settings such as retail forecasting, could involve hundreds of thousands of time series. More importantly, a desirable practical feature for multi-variate forecasting models is to let the prediction of future values for particular time series only require as input historical data from that time series (along with covariates), without requiring access to historical data from all other time series in the hierarchy. This allows for scalable training and inference of such models using mini-batch gradient descent, without requiring each batch contain all the time series in the hierarchy. This is often not possible for GNN-based forecasting models, which require batch sizes of the order of the number of time series.
Based on the above motivations, goal is to design hierarchical forecasting model with the following requirements: 1) The model can be trained using a single-stage pipeline on all the time series data, without any separate post-processing. 2) The model captures the additive coherence constraints along the edges of the hierarchy, 3) The model is efficiently trainable on large datasets, without requiring, for instance, batch sizes that scale with the number of time series.
Problem Statement:
We propose a principled methodology to address all these above requirements for hierarchical fore- casting. Our model comprises of two components, both of which can support coherence constraints. The first component is purely function of the historical values of time series, without distin- guishing between the individual time series themselves in any other way. Coherence constraints on such a model correspond to imposing an additivity property on the prediction function which constrains the model to be a linear autoregressive model. However, crucially, our model uses time- varying autoregressive coefficients that can themselves be nonlinear functions of the timestamp and other global features (linear versions of varying AR have been historically used to deal with non-stationary signals (Sharman & Friedlander, 1984)). We will refer to this component as the time-varying autoregressive model.
The second component focuses on modeling the global temporal patterns in the dataset through identifying a small set of temporal global basis functions. The basis time-series, when combined in different ways, can express the individual dynamics of each time series. In our model, the basis time-series are encoded in a trained seq-2-seq model (Sutskever et al., 2014) model in a functional form. Each time series is then associated with a learned embedding vector that specifies the weights for dccomposition along these basis functions. Predicting a time series into the future using this model then just involves extrapolating the global basis functions and combining them using its weight vector, without explicitly using the past values of that time series. The coherence constraints therefore only impose constraints on the embedding vectors of each time series, which can be easily modeled by a hierarchical regularization function. We call this component basis decomposition model. As we will see, this part of the model is only approximately coherent unless the embedding constraints hold exactly. In particular, in this paper, we focus on improving model accuracy rather than preserving
2

---

Under review as a conference paper at ICLR 2022
exact coherency. In Section A.2, also provide theoretical justification for how such hierarchical regularization using basis decomposition results in improved prediction accuracy.
We experimentally evaluate our model on multiple publicly available hicrarchical forecasting datasets. We compare our approach to state-of-the-a models and reconciliation models, and show that our approach can obtain consistently more accurate predictions at all levels of the hierarchy tree.
monnhhirrarhicca) deep forecasting models, GNN-basce
2 RELATED WORK ON DEEP HIERARCHICAL, MODELS
In addition to the works referenced in the previous section, we now discuss a few papers that are more relevant to our approach. Specifically, we discuss some recent deep hierarchical forecasting methods that do not require post-processing reconciliation step. Hierarchical forecasting methods can be roughly divided into two categories: point forecasters and probabilistic forecasters. Mishchenko et al. (2019) propose point-forecasting approach which imposes coherency a base model via le regularization on the predictions. Gleason (2020) extend the idea further to impose the hierarchy embedding space rather than the predictions directly. SHARQ (Han et al., 2021) follows similar lg regularization based approach as Mishchenko et al. (2019), and also extends the idea to probabilistic forecasting. Their model is trained separately for each of the hierarchical levels starting from the leaf level, thus requiring a separate prediction model for each level.
Probabilistic forecasting methods include, Hier-E2E (Rangapuram al., 2021) which produces perfectly coherent forecasts by using projection operation on base predictions from a DeepVAR model (Salinas et al., 2019). It requires the whole hierarchy of time series to be fed as input to the model leading to a large number of parameters, and hence does not scale well to large hierarchies. Yanchenko al. (2021) take a fully Bayesian approach by modelling the hierarchy using conditional distributions.
3 PROBLEM SETTING
We are given a set of N coherent time series of length T. arranged in a pre-defined hierarchy consisting of nodes. At time step I, the time series data can be represented as a vector yr E R denoting the time series values of all N nodes. We compactly denote the set of time series for all
[yr..... yrr E RTTN, TPxNNANN Also define ( as the ith column of the matrix Y denoting all time steps of the th time series, and y, as the t th value of the i th time series. We compactly denote the H-step history of Y by Yu [y:-H, - E # and the H H-step history of y by VH V ". Vr- 1] R Similarly define the F-step future of Y as Y [y. ... --+F-1] 1T G FFN, We use the notation to denote predicted values, for example Y,, ÿr and 9r.
T steps as a matrix
Time series forecasts can often be improved by using features as input to the model along with historical time series. The features often evolve with time, for example, categorical features such as type of holiday, or continuous features such as time of the day. We denote the matrix of such features by X E RTXD where the f th row denotes the D-dimensional feature vector at the t time step. For simplicity, assume that the features are global, meaning that they are shared across all time series. We similarly define X and X5 y as above.
Hierarchically Coherent Time Series. We assume that the time series data are coherent, that is. they satisfy the sum constraints over the hierarchy. The time series at each node of the hierarchy is the equal to the sum of the time series of its children, or equivalently, equal to the sum of the leaf time series of the sub-tres rooted at that node. Figure 1 shows an example of a sub-tree rooted at node.
As a result of aggregation, the data can have widely varying scales with the values at higher level nodes being magnitudes higher than the leaf level nodes. It is well known that neural network training is more efficient if the data are similarly scaled. Hence, in this paper, we work with rescaled time series data. The time series at each node is downscaled by the number of leaves in the sub-tree rooted at the node, so that now they satisfy mean constraints rather than sum constraints described above. Denote by L(p), the set of leaf nodes of the sub-tree rooted at p. Hierarchically coherent data satisfy
3

---

Under review as a conference paper at ICLR 2022
the following data mean property,
y(m)
1
Σ y()) iEC(P)
(Data Mean Property).
(1)
[(p)
4 HIERARCHICALLY REGULARIZED DEEP FORECASTING - HIRED
We now introduce the two components in our model, namely the time-varying AR model and the basis decompasition model. As mentioned in the introduction combination of these two components satisfy the three requirements in our problem statement. In particular, we shall see that the coherence property plays central role in both the components. For simplicity, in this section, all our equations will be for forecasting one step into the future (F = 1), even though all the ideas trivially extend to multi-step forecasting. The defining equation of our model can be written as,
93 = f(y7 XH, Fx, 0.)
(1) YH al a(xx, X,, ZHe) + (0, H: X ZH!
Time varying AR (TVAR)
Basis decomposition (BD)
In the above equation, Zy is a latent state vector that contains some summary temporal informa- tion about the whole dataset, and 0; is the embeddinn/weeigh vector for time-series i in the basis decomposition model. Zu can be a relatively low-dimensional temporally evolving variable that represents some information about the global state of the dataset at a particular time. We use the Non-Negative Matrix Factorization (NMF) algorithm by Gillis & Vavasis (2013) to select a small set of representative time series that encode the global state. If the indices of the selected representative time-series is denoted by (i1,... ,in) (R denotes the rank of the factorization), then we define Z = -1) .... (+)) E Tx Note that we only feed the past values ZH as input to the model, since future values are not available during forecasting. Also, note that the final basis time-series can be a non-linear function of Zz. In our experiments, R is tuned but is always much much less than N. a and are functions not dependent on (1) and 0,. We will provide more details as we delve into individual components.
Time- Vore-narying... AR (TVAR): The first part of the expression in equation 2 denoted by Time Varying AR (TVAR) resembles a linear auto-regressive model with coefficients al X Ht Xx, Zx) E RH that are a function of the input features, and thus can change with time. The AR parameters of this model are shared across all time series and hence do not encode any time-series specific information, a drawback that is overcome by the Basis Decomposition part our model. This component is coherent by design because it is shared linear AR model. However, even though the AR weights are shared across all the time-series at a given time-point, they can crucially change with time, thus lending more flexibility to the model.
Implementation: In order to model the sequential nature of the data, we use an LSTM encoder to encode the past X and ZH- Then, we use fully connected (FC) decoder for predicting the auto-regressive weights. Similar to Wen et al. (2017)'s hoult--horizon approach, we use different head of the decoder for each future time step resulting in a F-headed decoder producing F-step predictions for TVAR weights. The decoder also takes as input the future covariates Xy: y if available. The produced weights are then multiplied (inner product) to the history to produce the final TVAR predictions. We illustrate this architecture in Figure 2 (right).
Basis Decomposition (BD) with Hierarchical Regularization: Now we come the second part of our model in equation 2. As discussed before, this part of the model has per time-series adaptivity, as different time-series can have different embeddings. It resembles an expansion of the time series on set of basis functions ((X Xx, Z+) R2K, with the basis weights/embedding for time series i denoted by 0₂ E I Both the basis functions and the time series specific weights are learned from the data, rather than fixing a specifie form such as Fourier or Wavelet basis.
The idea of using a basis has also been recently invoked in the time-series literature (Sen et al., 2019; Wang et al., 2019). The basis recovered in the implementation of Wang et al. (2019) is allowed to vary for each individual time-series and therefore is not a true basis. Sen et al. (2019) do explicitly recover approximate basis in the training set through low-rank matrix factorization regularized
4

---

Under review as a conference paper at ICLR 2022
by deep global predictive model alternatingly trained on the basis vectors, thus not amenable to end-to-end optimization. We shall see that our model can be trained in an end-to-end manner.
Embedding Regularization for Approximate Coherency: The TVAR part of our model is coherent by design due to its linearity. The BD model however requires the embeddings of the time-series to satisfy the mean property along the hierarchy. This directly translates coherency of the predictions due to linearity with respect to 9.
8
Σ (EC(M)
0: (Embedding Mean Property).
(3)
[(P)
We impose this constraint approximately via an &2 regularization on the embedding.
Ereg ())
10 -::12
p=1 (EC(p)
The purpose of this regularizer is two fold. Firstly, we observe that, when the leaf embeddings are kept fixed, the regularizer minimized when the embeddings satisfy the mean property (3), thus encouraging coherency in the predictions. Secondly, it also encodes the inductive bias present in the data corresponding to the hierarchical additive constraints. We provide some theoretical justification for this hierarchical regularization in Section 5.
Implementation: As before, we use an LSTM encoder to encode the past Xy. and Zy. Then, we use the encoding from the encoder along with the future features X, (sequential in nature) and pass them through an LSTM decoder to yield the F-step basis predictions which are then multiplied with the embedding (inner product) to produce the final BD predictions. We illustrate this architecture in the top right of Figure 2. Thus, a functional representation of the basis time-series is implicitly maintained within the trained weights of the basis generating seq-2-seq model. Note that the embeddings are also trained in our end-to-end model. We illustrate this architecture in Figure 2 (left).
We emphasize that the main ideas in our model are agnostic to the specific type of neural network architecture used. For our experiments, we specifically use an LSTM architecture (Hochreiter & Schmidhuber, 1997) for the encoder and decoder. Other types of architectures including transformers (Vaswani et al., 2017) and temporal convolution networks (Borovykh et al., 2017) can also be used.
Loss Function: During training, we minimize the mean absolute error (MAE) of the predictions along with the embedding regularization term introduced above (our method generalizes to other losses too, such as mean square error, or mean absolute percentage error). For regularization weight Ag, and y defined as Eq (2), and θ denoting the trainable parameters of a, b, our training loss function is,
E(0, 0) = ΣV y
E
(5)
Emhedking regularization
Prediction loss
Note that the time-series dependent part of the loss function can be easily mini-batched and the embeddings are not memory intensive.
5 THEORETICAL JUSTIFICATION FOR HIERARCHICAL MODELING
In this section, we theoretically analyze the benefits of modeling hierarchical constraints in a much simplified setting, and show how it can result provably improved accuracy, under some assumptions. Since analyzing our actual deep non-lineat model for an arbitrary hierarchical set of time series can be complex, we make some simplifying assumptions to the problem and model. We assume that all the time series in the dataset is linear combination of a small set of basis time series. That is, Tx denotes the set of basis vectors, = [0 Ov] E K denotes the set of weight vectors used in the linear combination for each time series, and E RV denotes the noise matrix sampled i.id as w ~ (O, o²) for the leaf nodes. A classical example of such basis set can be small subset of Fourier or Wavelet basis (Strang, 1993; van den Oord et al., 2016) that is relevant to the dataset. Note that we ignore the TVAR model for the sake of analysis and focus mainly the BD model which includes the hierarchical regularization.
Y = B0 + w. where
5

---

Under review as a conference paper at ICLR 2022
Decompostion
Varying
Preducted
(Pegg
Embecciings
Mu Millrraraoo cwaatdr
solies
Futre feaures
Futur
Fiu 8D Prddiie
by
Carsarr
Past
Fast ckibal stalus
Xx
states
ER
R
Figure 2: In this figure we show the architectures of our two model components separately. On the left we show the BD model, where the seq-2-seq model implicitly maintains the basis in a functional form. Note that the time-series specific weights [0:] are also trained. On the right, we show the TVAR model. The fully connected decoder has different prediction head for each future time-point.
Algorithm 1: Basis Recovery
Algorithm 2: Parameter Recovery
Input: Observed y, basis dict B, regularization parameter d Output: Estimated basis B 0g
Input: Observed time series y. estimated basis B, regularization parameter λ
Output: Estimated parameters 0 0o- argming, 7 Fllyo B0g l1 for L()) do
1 argmin 27 aER"
ll/0 Ball3 + d A||lllla
0, argming, +ly. Be,, 112 (100 1
Estimate support S {i |00| 0}
end
Estimate true basis B < Bg
In this section, we consider a 2-level hierarchy of time-series, consisting of single root node (indexed by 0) with L children (denoted by ((0)). We will also assume that instead of learning the K basis vectors B from scratch, the basis vectors are assumed come from a much larger dictionary B E R of D (> K) vectors that is fixed and known to the model. While the original problem learns the basis and the coefficients simultaneously, in this case the goal is to select the basis from among larger dictionary, and learn the coefficients 0.
We analyze this problem, and show that under the reasonable assumption of the parent embedding & being close to all the children embeddings O,, using the hierarchical constraints result in mean- square error at the leaf nodes that a multiplicative factor L smaller than the optimal mean-square error of any model that does not use the hierarchical constraints. Our proposed HIRED model, when applied in this setting would result in the following (hierarchically) regularized regression problem:
min 0
NT y B
+ A
Σ CEC
110 - 8,
(6)
For the sake of analysis, we instead consider a two-stage version, described in Algorithm and Algorithm 2: first recover the support of the basis using Basis Pursuit (Chen al., 2001). We then estimate the parameters of the root node, which is then plugged-in to solve for the parameters of the children node. We also define the baseline (unregularized) optimization problem for the leaf nodes that does not use any hierarchical information, as
0, = argmin
T llnn - B8, Vn E C(0).
(7)
O,,
The basis support recovery follows from standard analysis (Wainwright, 2009) detailed in Lemma in the Appendix. We focus on the performance of Algorithm 2 here. The following theorem bounds the error of the unregularized (0,, and the hicrarchically regullarie (0,,, see Algorithm 2) optimization
solutions. A proof of the theorem can be found in Appendix A.2.
Theorem 1. Suppose the rows of B are norm bounded as | 12 and |00,, - 00ll2 S B. Define R = BT B/T as the empirical covariance matrix. For Xg 7:22
Hn and 9, can be bounded as,
2KK T
0 2K
+ 6 2 K
0 - 13. <
E E - 0,1
<3
(8)
T
1
g2K
TL
Tr=3

---

Under review as a conference paper at ICLR 2022
In fixed design linear regression ||0,, - 6. = |B 0,)ll is the population squared error (see Appendix A.5 for a bound on the parameter estimation error). The gains due to the regularization can be understood by considering the case when 3 is upper bounded by sufficiently small quantity. Note that an upper bound on 33 essentially implies that the children time-series have structural similarities as further elaborated in Appendix A.5. We show that the above assumption yields a smaller upper bound on the error. In fact, if B = ((yK/T), then the numerator + Tr in Eq. (8) is (1)
o² K
resulting in BIJØ, - 1 = of - which decays faster than =²K K Furthermore, if /3 is even smaller as 3 = OVv K/ LT), then following similar calculations, E 0n 0 = 0 o² K smaller than the unregularized bound.
-
which again LT
EXPERIMENTS
We implemented our proposed model in Tensorflow (Abadi et al., 2016) and compared against multiple baselines on popular hierarchical time-series datasets.
Datasets. We experimented with three hierarchical forecasting datasets Two retail forecasting datasets, M5 (M5, 2020) and Favorita (Favorita, 2017); and the Tourism (Tourism. 2019) dataset consisting of tourist count data. The history length and forecast horizon (H, F) were set to (28, 7), (28, 7) and (24, 4), for Favorita, M5 and Tourism respectively. More information can be found in Appendix B.2. We divide each of the datasets into training, validation and test sets, with details on the splits provided in Appendix B.3.
Baselines. We compare our proposed approach HIRED with the following baseline models: (i) RNN we use seq-2-seq model shared across all the time series, (ii) DeepGLO (Sen et al., 2019), (iii) DCRNN (Li et al., 2017), a GNN based approach where we feed the hierarchy tree as the input graph, (iv) Deep Factors (DF) (Wang et al., 2019), (v) Ly (Gleason, 2020), which is an improvement over Mishchenko et al. (2019), (vi) SHARQ (Han et al., 2021), (vii) HierE2E (Rangapuram et al., 2021) the standard implementation produces probabilistic forecasts, however adapt it point forecasts by using their projection step on top of seq-2-seq model. We use code publicly released by the authors for DeepGLO' and DCRNN We implemented our own version of DeepFactors, HierE2E, and Emb for a fair comparison, since the official implementations either make rolling probabilistic forecasts, or use different set of covariates. Additionally, we also compare with the recent ERM (Ben Taieb & Koo, 2019) reconciliation method applied to the base forecasts from the RNN model, denoted as (vii) RNN+ It has been shown (Ben Taieb & Koo, 2019) to outperform many previous reconciliation techniques such as MinT (Wickramasuriya et al., 2019). For a fair comparison, we use the same Mean Absolute Error (MAE) loss for all the compared methods, and make sure that all the models have approximately the same number of parameters. Further details about the basclines and training parameters can be found in Appendix B.
Metrics. We compare the accuracy of the various approaches with respect to two metries: (i) weighted absolute percentage error (WAPE), and (ii) symmetric mean absolute percentage error (SMAPE). A description of these metries can be found in Appendix B.1. We report the metrics on the test data, for each level of the hierarchy (with level denoting the root) in Table 1. As a measure of the aggregate performance across all the levels, also report the mean of the metries in all the levels of the hicrarchy denoted by Mean.
6.1 RESULTS
Table I shows the averaged test metrics for M5, Favorita, and Tourism datasets. We present only the Mean metrics for the Tourism dataset due to lack of space. Complete results with confidence intervals for all the three datasets can be be found in Appendix B.6.
We find that for all three datasets, our proposed model either yields the smallest error or close the smallest error across most metrics and most levels. In particular, we find that our proposed method achieves the smallest errors in the mean column for all datasets in terms of WAPE and SMAPE, thus indicating good performance generally across all levels. We find that RNN+ERM in general, yields
'https://git httph//..ooorrajajaaaaaaea /deepglo
hhttps://g / thub. com wttptt.co/m/aaggaagg//o/o//// htw.////yyggggg//////////
7

---

Under review as a conference paper at ICLR 2022
Table 1: The tables show the WAPE/SMAPE test metrics for the M5, Tourism, and Favorita datasets, averaged over 10 runs. We present only the Mean metrics for the Tourism dataset due to lack of space. A complete set of results including the standard deviations can be found in Appendix B.6. In-coherent and coherent baselines are separated by horizontal line.
Level
Level
Level
Level 3
Mean
Tourism
Mean
0.048
0.055
0.072/
0.279/0.511
0.172
0.186/ 0.322
RNN
0059/0099
0.083 /0.083
0.085 /0.08
0.252/0.517
0.127 /
0.21//0.335
DF
0.055/0.056
0,061 /0,060
0.076 /0005
0,272 00.50
0.116/0 0.176
DF
0.204/0.334
DeepGLD
0.027/0.081
0,087 0.092
0.099 / 00.13
0.278 /0.338
0.135 / 0.2006
DeexGLO
0.199/ 0.346
DCRNN
0.078/0.079
0,096/00004 00009
0.165 / 090
0.282 /0.512
0.156 /0.219
DCRNN
0.281 00,39
Lg Kimb
0.055 /0.056
0,064 /0063
0.080/0.000
/..50
0.117/0.175
Emg 0.215/0.342
SHARQ
0.00/00006
0.071 70.062
0.009 /0.094
0.277 00.538
0.135/0195
SHARQ
0.229/0.378
RNN-ERM
0.052/
0,066
0.286 /0.520
0.122/ 0.117
RNN+ERM
0.211/04177
Hlee-E26
0.12/0..60
0..15/
0.152/0.181
0.396/0.615
0.13//02778
E2ce-E2E
0,208/0,340
Faworica
Level
Lovel
Level
Level
Mean
IIRED
0,061 /0.061
0.094/ 0.182
0.127 /0.267
0,210/ 0.302
0.123
RNN
0.067 0.068
0.114/0.197
0.134 00.290
0.203 /0.399
0.130/0.220
TNF
0.064/ / 0.064
0.110/0.194
0.135 / 00.91
0.213/0.543
0.130 0.230/0.223
DeepGLO
0.098 0.088
0.126/0.197
0.156 0.156/0.333 0,333
0,226
0.15L/0.255 70,256
DCRNN
0.080/0.080
0.120/0.212
0.134 0.134/0,328
/0.389
0.134 /0,252
LyEib
0.070 0.070
0.116/0.199
0.136 0,276
0.207 00.321
0.132 / 002216
SHARQ
0,058 0.085
0.142/0.199
0.156/0.335
0.230/0.404
0.154/0.256
RNN-ERM
0/856 0.058
0.103/0.185
0.129
0.220/0.344
0.127/0.219
Hien E2E
0.120
0,206 /0,334
0.247 /0.488
00.59/
0.245 /0,370
an improvement over the base RNN predictions for the higher levels closer to the root node (Levels 0 and 1), while, worsening at the lower levels. DCRNN, despite using the hierarchy as a graph also does not perform as well as our approach, especially in Tourism and M5 Datasets possibly due to the fact that a GNN is not the most effective way to model the tree hierarchies. We notice that HierE2E performs reasonably well for the smaller Tourism while performing badly for the larger M5 and Favorita datasets a possible explanation being that this is a VAR model that requires much more parameters to scale to large datasets. Therefore, for HierE2E perform experiments with 50x more parameters for M5 and Favorita and report the results in Table 4 in the appendix, showing that while the results improve, it still performs much worse than our model. Overall, find that our proposed method consistently works better or at par with the other baselines at all hierarchical levels.
Ablation study. Next, we perform an ablation study of our proposed model to understand the effects of its various components, the results of which are presented in Table 2. We compare our proposed model, to the same model without any regularization (set λø O in Eq (5)), and a model consisting of only TVAR. We find that both these components in our model are important, and result in improved accuracy in most metries.
Table 2: We report the test WAPE/SMAPE metrics for an ablation study on the M5 dataset, for each of the components in the HiReD model. We compare our model with two ablated variants: first, we remove the regularization (Ag 0), and second, we remove the BD component (TVAR only).
Level
Level
Level
3
Mena
IIRED
0.048/ 0.148
0.055 0.053
0,072/ 0.077
0.279 00.511
0.113/
>g
0,054 0.054
0.054
0.058 /0056
0,0744 0.078
0.279 00.533
00175
TVAR only
0,050/ 0.049
0.064/00065
0.084/00.80
0. 0.788/ 00.530
1.122./0.180
8

---

Under review as a conference paper at ICLR 2022
Levi3
HiRrD
Favoriaa
A
0.012
0013
0,000
RNN
0.045
0.044
0.042
HiReD
0,030
0.034
0.004
M5
Ag
0:035
0.040
0,099
0.042
0.047
IiReD
0,092
0.079
0,066
0,060
Tourism
0.097
0.082
0.083
Lovel
Figure Coherency metric for all our datasets, at all hierarchical levels. Leaf node metrics are identically zero, and hence not reported in the table. Leaf nodes for Favorita and M5 are denoted by L3. Tourism has hierarchical levels and hence L3 values reported in this case.
Figure 3: Left: Plots of the basis generated the validation set of the M5 dataset over 35 days. Right: We plot the true time series over the same time period, and compare it with the predicted time series, AR predictions and BD predictions.
Coherence. We also compare the coherence of our predictions that of the RNN model and an ablated model with λÆ O. Specifically, for each node p we measure the deviation of our forecast from c(p) = 1/L(p) EiE C(p) ((), the mean of the leaf node predictions of the corresponding sub- tree. Perfectly coherent predictions will have a zero deviation from this quantity. In Table 4, we report the WAPE metric between the predictions from our model ÿ and c, for each of the hierarchical levels. The leaf level predictions are trivially coherent. We find that our proposed model consistently produces more coherent predictions compared both the models, indicating that our hierarchical regularization indeed encourages coherency in predictions, in addition to improving accuracy.
Basis Visualization. We visualize the basis generated by the BD model for the M5 validation set in Figure 3 (left). We notice that the bases capture various global temporal patterns in the dataset. In particular, most of the bases have a period of 7, indicating that they represent weekly patterns. We also show the predictions made from the various components of our model at all hierarchical levels, in Figure 3 (right). We notice that the predictions from the BD part closely resemble the general patterns of the true time series values, where as the AR model adds further adjustments the predictions, including a constant bias, for most time series. For the leaf level (Level 3) predictions however, the final prediction is dominated by the AR model indicating that global temporal pattems may be less useful in this case.
7 CONCLUSION
In this paper, we proposed method for hierarchical time series forecasting, consisting of two components, the TVAR model, and the BD model. The TVAR model is coherent by design, whereas we regularize the BD model to impose approximate coherency. Our model is fully differentiable and is trainable via SGD, while also being scalable with respect to the number of nodes in the hierarchy. Furthermore, it also does not require any additional pre-processing steps.
We empirically evaluated our method on three benchmark datasets and showed that our model consistently improved over state of the art baselines for most levels of the hierarchy. We perform an ablation study to justify the important components of our model and also show empirically that our forecasts are more coherent than the RNN baseline, Lastly, we also visualize the learned basis and observe that they capture various global temporal patterns.
In this work, we aimed at maximizing the overall performance without emphasizing performance at individual hierarchical levels. For future work, plan to treat this is multi-objective problem with the aim of understanding the performance tradeoffs at various levels of the hierarchy. Finally, also plan extend our current model to probabilistic forecasting.

---

Under review as a conference paper at ICLR 2022
Reproducibility statement: We provide code for our model along with instructions in the supple- mentary file. The instructions are provided in the readme txt file inside the zipped folder. Our code includes links to public datasets and code for preprocessing them. We also fix the random seed in our code in order for to be reproducible. However, exact results may depend on the version of the python libraries, type of GPU and compute used.
REFERENCES
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving. Michael Isard, et al. Tensorflow: system for large-scale machine learning. In /2th USENIX symposium on operating systems design and implementation (OSDI 16), pp. 265-283, 2016.
Lei Bai, Lina Yao, Can Li, Xianzhi Wang, and Can Wang. Adaptive graph convolutional recurrent network for traffic forecasting. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 17804-17815. Cur- ran Associates, Inc., 2020. URL https: http://poocecee..cceeeeppoo neurips. cc cattaps.ccpppppeee wwtppppe///2///00/ http://eaaaaee992faaaaaaa0aa...0
file/celaad92b3 ww290332222c1111 www.:0occc10005541111210 -
Souhaib Ben Taieb and Bonsoo Koo. Regularized regression for hierarchical forecasting without unbiasedness conditions. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1337-1347, 2019.
Konstantinos Benidis, Syama Sundar Rangapuram, Valentin Flunkert, Bernie Wang, Danielle Maddix, Caner Turkmen, Jan Gasthaus, Michael Sohlke-Schneidee David Salinas, Lorenzo Stella, al. Neural forecasting: Introduction and literature overview. arXiv preprint arXiv: 10240, 2020.
Aleksandar Bojchevski, Johannes Klicpera, Bryan Perozzi, Amol Kapoor, Martin Blais, Benedek Rózemberczki, Michal Lukasik, and Stephan Günnemann. Scaling graph neural networks with approximate pagerank. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 2464-2473, 2020.
Anastasia Borovykh, Sander Bohte, and Comelis W Oosterlee. Conditional time series forecasting with convolutional neural networks. arXiv preprint arXiv:1703.04600 2017.
Joos-Hendrik Böse, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Dustin Lange, David Salinas, Sebastian Schelter, Matthias Seeger, and Yuyang Wang. Probabilistic demand forecasting at scale. Proceedings of the VLDB Endowment, 10(12):1694-1700 2017.
Defu Cao, Yujing Wang. Juanyong Duan, Ce Zhang, Xia Zhu, Congrui Huang. Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, and Qi Zhang. Spectral temporal graph neural network for multivariate time-series forecasting. In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 17766-17778. Curran As- rewriis.eeris.... cc/ paper/2020000 fiie// sociates, Inc., 2020, URL https: Paper. pdf. (tt.e5a.cc7cccceeeeeeeeeee33eee.
Scott Shaobing Chen, David I. Donoho, and Michael A Saunders. Atomic decomposition by basis pursuit. SIAM review, 43(1):129-159) 2001.
Ying Cui, Ruofei Zhang, Wei Li, and Jianchang Mao. Bid landscape forecasting in online ad exchange marketplace. In Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 265-273, 2011.
Emmanuel de Bézenac, Syama Sundar Rangapuram, Konstantinos Benidis, Michael Bohlke- Schneider, Richard Kurle, Lorenzo Stella, Hilaf Hasson, Patrick Gallinari, and Tim Januschowski. Normalizing kalman filters for multivariate time series analysis. In NewrlPS, 2020.
Favorita. Favorita forecasting dataset.
https: hetps::/www..r /ww. www.:wwwwaggaaa..o baaglel e waw.gglee.occ......o.oo/o////
favorita-grocery- www.y-sssll---sesssssseeessssses 2017.
Nicolas Gillis and Stephen A Vavasis. Fast and robust recursive algorithmsfor separable nonnegative matrix factorization. IEEE transactions on pattern analysis and machine intelligence, 36(4):
698-714, 2013.
10

---

Under review as a conference paper at ICLR 2022
Jeffrey L Gleason. Forecasting hierarchical time series with regularized embedding space. San Diego, 7, 2020.
Xing Han, Sambarta Dasgupta, and Joydeep Ghosh. Simultaneously reconciled quantile forecasting of hierarchically related time series. In International Conference on Artificial Intelligence and Stalistics, pp. 190-198. PMLR, 2021.
Sepp Hochreiter and Jürgen Schmidhuber, Long short-term memory. Neural computation, 9(8):
1997.
Daniel Hsu, Sham M Kakade, and Tong Zhang. Random design analysis of ridge regression. In Conference on learning theory, pp. 9-1. JMLR Workshop and Conference Proceedings, 2012.
Rob Hyndman, Anne B Koehler, J Keith Ord, and Ralph D Snyder. Forecasting with exponential smoothing: the state space approach. Springer Science & Business Media, 2008
Rob J Hyndman and George Athanasopoulos. Forecasting: principles and practice. OTexts, 2018.
Rob Hyndman and Shu Fan. Density forecasting for long-term peak electricity demand. IEEE Transactions on Power Systems, 25(2):1 142-1153, 2009.
Rob Hyndman, Roman A Ahmed, George Athanasopoulos, and Han Lin Shang. Optimal com- bination forecasts for hierarchical time series. Computational statistics & data analysis, 55(9): 2579-2589, 2011.
Rob J Hyndman, Alan Lee, and Earo Wang. Fast computation reconciled forecasts for hicrarchical and grouped time series. Computational statistics & data analysis, 97:16-32, 2016.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. Diffusion convolutional recurrent neural network: Data-driven traffic forecasting. arXiv preprint
2017.
M5.
M5 forecasting dataset.
https: / //www/ www..kaggle.co com/
m5- wwt.eastiingccccaaaa./.ooo 2020.
ED McKenzie. General exponential smoothing and the equivalent arma process. Journal of Forecast- ing, 333):333-344, 1984.
Konstantin Mishchenko, Mallory Montgomery, and Federico Vaggi. A self-supervised approach to hierarchical forecasting with applications to groupwise synthetic controls. arXiv preprint arXiv:1906.10538 2019.
Nam Nguyen and Brian Quanz. Temporal latent auto-encoder: method for probabilistic multivariate time series forecasting, In Proceedings of the AAAI Conference Artificial Intelligence, volume 35, pp. 9117-9125, 2021.
Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. N-beats: Neural basis expansion analysis for interpretable time series forecasting. arXiv preprint 1rXiv:1905.10433 10437, 2019.
Anastasios Panagiotelis, Puwasala Gamakumara, George Athanasopoulos, Rob Hyndman, al. Probabilistic forecast reconciliation: Properties, evaluation and score optimisation. Monash econometrics and business statistics working paper series, 26:20, 2020.
Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski. Deep state space models for time series forecasting. Advances neural information processing systems, 31:7785-7794, 2018.
Syama Sundar Rangapuram, Lucien D Werner, Konstantinos Benidis, Pedro Mercado, Jan Gasthaus, and Tim Januschowski. End-to-end learning of coherent probabilistic forecasts for hierarchical time series. In International Conference Machine Learning, pp. 8832-8843. PMLR, 2021.

---

Under review as a conference paper at ICLR 2022
Kashif Rasul, Abdul-Saboor Sheikh, Ingmar Schuster, Urs Bergmann, and Roland Vollgraf. Mul- tivariate probabilistic time series forecasting via conditioned normalizing flows. arXiv preprint arXiv:2002.066103 2020.
David Salinas, Michael Schnee-, Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank gaussian copula processes. arXiv preprint arXiv:1910.03002. 2019.
David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3): 1181-1191, 2020.
Rajat Sen, Hsiang-Fu Yu, and Inderjit Dhillon. Think globally, act locally; A deep neural network approach high-dimensional time series forecasting. arXiv preprint arXiv:/905.000
2019.
K Sharman and Benjamin Friedlander. Time-varying autoregressive modeling of class of nonsta- tionary signals. In ICASSP 84. IEEE International Conference on Acoustics, Speech, and Signal Processing, volume 9, pp. 227-230. IEEE, 1984.
Gilbert Strang. Wavelet transforms versus fourier transforms. Bulletin of the American Mathematical Society, 28(2):288-3)) -008 1993.
Ilya Sutskever, Oriol Vinyals, and Quoc V Le. Sequence to sequence learning with neural networks. Advances in Neural Information Processing Systems (NIPS 2014), 2014.
Souhaib Ben Taieb, James W Taylor, and Rob J Hyndman. Coherent probabilistic forecasts for hierarchical time series. In International Conference on Machine Learning, pp. 3348 -357. PMLR, 2017.
Tourism. Tourism forecasting dataset. mint/.201
Aäron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu. Wavenet: A generative model for raw audio. In 9th ISCA Speech Synthesis Workshop, pp. 125-125, 2016.
Tim Van Erven and Jairo Cugliari. thcretlco optimal reconciliation of contemporancous hierarchical time series forecasts. In Modeling and stochastic learning for forecasting in high dimensions, pp. 297-317. Springer, 2015.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention all you need. In NIPS, 2017.
Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using f Ej-constra quadratic programming (lasso). IEEE transactions on information theory, 55(5): 2183-2202, 2009.
Martin Wainwright, Higg-diienssonal statistics: non-asymptotic viewpoint, volume 48. Cam- bridge University Press, 2019.
Yuyang Wang, Alex Smola, Danielle Maddix, Jan Gasthaus, Dean Foster, and Tim Januschowski. Deep factors for forecasting. In International Conference on Machine Learning, pp. 6607-6617.
PMLR, 2019.
Ruofeng Wen, Kari Torkkola, Balakrishnan Narayanaswamy, and Dhruv Madeka. multi-horizon quantile recurrent forecaster. NeurIPS Time Series Workshop, 2017.
Shanika I. Wickramasuriya, George Athanasopoulos, Rob J Hyndman, et al. Forecasting hierarchical and grouped time series through trace minimization. Department of Econometrics and Business Statistics, Monash University, 105, 2015.
Shanika L Wickramasuriya, George Athanasopoulos, and Rob J Hyndman. Optimal forecast rec- onciliation for hierarchical and grouped time series through trace minimization. Journal of the American Statistical Association, 114(526)804-811) 2019.
12

---

Under review as a conference paper at ICLR 2022
Shanika L Wickramasuriya, Berwin A Turlach, and Rob J Hyndman. Optimal non-negative forecast reconciliation. Staristics and Computing. 3065)11177-1122, 2020.
Kevin W Wilson, Bhiksha Raj, and Paris Smaragdis. Regularized non-negative matrix factorization with temporal dependencies for speech denoising. In Interspeech, pp. 41 411-414, 2008
Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang. Xiaojun Chang, and Chengqi Zhang. Connecting the dots: Multivariate time series forecasting with graph neural networks. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp- 753-763, 2020.
Anna K Yanchenko, Di Daniel Deng, Jinglan Li, Andrew J Cron, and Mike West. Hierarchical dynamic modeling for individualized bayesian forecasting. arXiv preprint arXiv:2101.03008, 2021.
Bing Yu, Haoteng Yin, and Zhanxing Zhu. Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting. arXiv preprint arXiv:1709.04891 2017.
Dawei Zhou, Lecheng Zheng, Yada Zhu, Jianbo Li, and Jingrui He. Domain adaptive multi-modality neural attention network for financial forecasting. Proceedings The Web Conference 2020, pp- 2230-2240 2020.
13

---

Under review as a conference paper at ICLR 2022
A THEORY
A.1 SUPPORT RECOVERY
Lemma 1. Suppose B satisfies the lower eigenvalue condition (Assumption 1 in Appendix A.3) with parameter Cmin and the mutual incoherence condition (Assumption 2 in Appendix A.3) with parameter Also assume that the columns of the basis pool B are normalized so that maxjes= [B S VT, and the true parameter Bg of the root satisfies
1001|00 2 d
||
|
4a
(9)
VLC
where || = max; Σ, Aijl denotes matrix operator foo norm, and E B² B/T denotes the
empirical covariance matrix. Then for X 2 v 20 load with probability least 1 - 4 exp((=c7TXX
LT
(for some constant e1), the support S = {i Jaol Algorithm 1) is equal to the true support
0} recovered from the Lasso solution (see
Proof. We are given pool of basis vectors B from which the observed data is generated using a subset of K columns which we have denoted by B in the text. We denote the correct subset of columns by S and recover them from the observed data using basis pursuit also known as the support recovery problem in the literature. Given the observed data and the pool of basis vectors B, we recover the support from the following regression problem corresponding the root node time series.
(10) where is K -spaase with the non-zero indices at S, and the non-zero values equal 9g the true parameters of the root node. Here we have used the fact that the root node has a 1/ times smaller variance due aggregation. The true support S can be recovered from the observed data yo. by solving the sparse regression problem (Lasso) given in Algorithm 2. A number of standard Lasso assumptions are needed ensure that the support identifiable, and that the non-zero parameters are large enough to be estimated. Assuming that Bş (= B) and a satisfy all the assumptions of Theorem 2, the theorem ensures that the true support is recovered with high probability.
Vo Ba + wo: Wo ~ ~N(, I/I.),
A.2 PROOF OF THEOREM 1 ERROR BOUNDS FOR REGULARIZED ESTIMATORS
For this proof, assume that the support S is recovered and the true basis functions B are known with high probability (see Section A.I). We divide the proof into multiple steps.
Step I: By Corollary 1, the OLS estimate Og (see Algorithm 2) of parameters of the root node and the OLS estimate O, (see Eq. (7)) can be bounded as,
0² K TL
a 0 <
E - 6 1
El
Vn E C(0).
(11)
T
Step II: Next, using change of variables, we notice that the ridge regression loss for the child nodes (see Algorithm 2) is equivalent to the following.
Gn argmin
T Ilv B
Bri 12 All Vn L()),
(12)
where vi, = 0,, G.. The final estimate for the child parameters can be written sum of the the
estimate and the root node estimate, O, = 2h + Og. We also consider related problem that will help us in computing the errors bounds.
Hn argmin *h
+ llv
BO, B + Al Vn E L()).
(13)
Here we have replaced tn with the true value dg. Note that this regression problem cannot be solved in practice since do not have access to the true value of lg. We will only it to assist in the
14

---

Under review as a conference paper at ICLR 2022
analysis. Now will bound the difference between the estimates on and din- The closed form solution for ridge regression is well known in the literature.
on = ( - AD 'B" (yn Byg)
& (E) BB (yn Bob),
where Σ B" B/T, as defined earlier. The norm of the difference of the estimates can be bounded
as
(2) 0o)
Vn Un
11 - F = (0o X2 11)-11106 00) + -100 - 0o)
= (0o
X
= (0o - -))
V (00 - -6))
(A: + X)2
Here we have used eigen-decomposition of the symmetric sample covariance matrix Σ V VDDA.VV. We use the notation D [>ss] to denote diagonal matrix with values ) on the diagonal. The above can be further upper bounded using the fact that X; S X + X.
1 1 < (0. ".) (0, 0,) 10 0.1
(14)
Step III: Now we will bound the error on on and finally use in the next step with triangle inequality to prove our result. Note that yn B0p BO,, 0o) - w.. Therefore, we can see from
Eq. (13) that vin is estimate for On 0o. Using the fact that ||0, 0olla and corollary th can be bounded as,
BBOKK
E Elll 5 (0n. 00)1 <
(15)
Tr2 o²K
Finally using triangle inequality, we bound the error of our estimate 6n-
10,- - 0x = V + 0 - 0,113 (Using the decomposition from Step II).
S (11e - 3. | - (0,, 0o) + 10 - 0 (13
(Using triangle and Sauchy-Schwarne inequality)
< (l (O,, 0o) 12: + 2 0oll3) (Using Eq. (14)).
Taking the expectation of the both sides, and using Eq. (11) and (15), we get the desired result.
Ellō. 1 < 3 Tr232 + o2K + 6 o2K
0²6 K
TL
A.3 REVIEW OF SPARSE LINEAR REGRESSION
We consider the following sparse recovery problem. We are given data (X, y) E Raxd x R following the observation model XO + w, where ~ NN(0, o²) I). indices indexed by set (S-sparse). We estimate using the following Lagrangian Lasso program,
and 0" is supported in the
{ 1 - Fly
i E argmin
X0
x
+ Am 1011
(16)
GER"
2n
We consider the fixed design setting, where the matrix X is fixed and not sampled randomly. Following (Wainwright, 2009), we make the following assumptions required for recovery of the true support S of 0*
Assumption (Lower eigenvalue). The smallest eigenvalue of the sample covariance sub-matrix indexed by S is bounded below:
X
Amin
2 Cuin > 0
(17)
15

---

Under review as a conference paper at ICLR 2022
Assumption 2 (Mutual incoherence). There exists some (0, 1] such that
| Xx 'I1 -1-7,
(18)
where | = max; E, |iil denotes matrix operator foc norm.
Theorem 2 (Support Recovery, Wainwright (2009)). Suppose the design matrix satisfies assumptions I and 2. Also assume that the design matrix has its n-dimensional columns normalized so that maxjese ||x I vT. Then for An satisfying,
2
20 log
An 2
(19)
the Lasso solution satisfies the following properties with probability at least 1 - exeepppp1n_)
1. The Lasso has unique optimal solution with its support contained within the true support S(Ō) C S(0*) and satisfies the c. bound
X
Xs
40
100g 03 S An
(20)
n
VCmia
g(...)
2. If in addition, the minimum value of the regression vector O is lower bounded by g g(An). then it recovers the exact support.
A.4 REVIEW OF RIDGE REGRESSION
In this section review the relevant background from (Hsu al., 2012) fixed design ridge regression. As usual, we assume data (X,y) E Raad x Rn following the observation model
y X0 + where NN If Define the ridge estimator 0 as the minimizer of the & regularized mean squared error,
0 E argmin
Ellllk x012
(21)
EER
We denote the sample covariance matrix by E = X X /n. Then for any parameter 0, the expected
l2 prediction error given by, |00 l| = X (0- 0*) 0-) We also assume the standard ridge regression setting of bounded ||0- l2 < B. We have the following proposition from Hsu et al. (2012) on expected error bounds for ridge regression.
Proposition 1 (Hsu et al. (2012)). For any regularization parameter À > O, the expected prediction loss can be upper bounded as
011
Σ
A
0
A A, + X
2
(22)
(A//A+ 1)2
n
where >; denote the eigenvalues of the empirical covariance matrix
Using the fact that X; S tr (E), and (x + c) is increasing x for x 2 0, the above bound can be simplified as
tr Σ
o²d n
tr (E)
E -110--011 11
(tr (2) / 1 8 10-1
(x +
tr (2) B tr (2:) d2d/n
(tr +
Assuming that the covariate vectors X, norm bounded as |X,112 r, and using the fact that tr (E)= < gives us the following corollary.
Corollary 1. When choosing À o²d nB2 the prediction loss can be upper bounded as,
o2d
e 811 < 22 + odd
(23)
The usual ordinary least squares bound of z²d can be derived when considering the limit + oc, corresponding to
16

---

Under review as a conference paper at ICLR 2022
A.5 FURTHER DISCUSSION ON THEOREM 1
Upper bound |0,,- - 0 2 I A: The upper bound & essentially bounds the distance between sibling leaf embeddings belonging the same parent. This is directly related to an upper bound on the distance between the parent embedding 9n and the leaf embeddings 8,, as 0, is essentially the mean of the leaf nodes (mean property). In many practical scenarios, the children time series of parent may not have too different seasonal trends (for example power consumption of houses in the same neighborhood, or sales of items under the same category) resulting the parent time series following similar trends as well.
Bounding | 00 In most theoretical analyses of linear regression (Wainwright, 2019), the main quantity of interest is the prediction error ||0,, -0ol's rather than the parameter estimation error ||0, - fol 12.12. as the former is directly related to the performance metric of the model. However, a bound on the parameter estimation error can be easily established using the property that ||0, - Boll2 I ||0, Vollx:///minn VCmin, where Cmin is the lower bound on the smallest eigenvalue of the sample covariance matrix as defined in Assumption in Section A.3.
B FURTHER EXPERIMENTAL DETAILS
B.1 ACCURACY METRICS
In this section we define the evaluation metrics used in this paper. Denote the true values by y and the predicted values by ÿ. both n-dimensional vectors.
1. Symmetric mean absolute percent error SMAPE- = Σ 91
2. Weighted absolute percentage error WAPE Σ, Σ, 19
B.2 DATASET DETAILS
We use three publicly available benchmark datasets for our experiments.
1.
The M5 dataset' consists of time series data of product sales from 10 Walmart stores in three US states. The data consists of two different hierarchies: the product hierarchy and store location hierarchy. For simplicity, in our experiments we use only the product hierarchy consisting of 3k nodes and 1.8k time steps. The validation scores are computed using the predictions from time steps 1843 to 1877, and test scores on steps 1878 to 1913.
2. The Favorita dataset is a similar dataset, consisting of time series data from Corporación Favorita, a South-American grocery store chain. As above, we use the product hierarchy, consisting of 4.5k nodes and 1.7k time steps. The validation scores are computed using the predictions from time steps 1618 to 1652, and test scores on steps 1653 1687.
3. The Australian Tourism dataset consists of monthly domestic tourist count data in Australia across 7 states which are sub-divided into regions, sub-regions, and visit-type. The data consists of around 500 nodes and 230 time steps. The validation scores are computed using the predictions from time steps 122 to 156, and test scores on steps 157 to 192.
For the three datasets, all the time-series (corresponding both leaf and higher-level nodes) of the hierarchy that we used are present in the training data.
B.3 TRAINING DETAILS
Table 3 presents all the hyperparameters used in our proposed model. All models were trained via SGD using the Adam optimizer (Kingma & Ba, 2014), and the training data was standardized mean zero and unit variance. The datasets were split into train, validation, and test sets, the sizes of which
'nttps:/ :/ htw.aag.....coooooooooooooetttte
http///wwkkkgggg/e com/ wttp:////aavvrraaarrrrrreeeeeeee woss-forraattinnggg hhttpsss/rdbbyyyddaacccc/// hwt.:oo//ii.ooo////////////t//
17

---

Under review as a conference paper at ICLR 2022
Table 3: Final model hyperparameters for various datasets tuned using the Mean WAPE metric on the validation set.
Model hyperparameters
M5
Favorita
Tourism
LSTM hidden dim
42
24
14
Embedding dim K
8
6
NMF rank It
12
4
6
Multi-Horizon decoder hidden dim
24
16
12
Embedding regularization Äg
3.4e-6
4.644e-4
7.2498e-8
History length H and forecast horizon F
(28, 7)
(28, 7)
(24, 4)
No. of rolling val/test windows
5
3
Initial learning rate
0.004
0.002
0.07
Decay rate and decay interval
(0.5, 6)
(0.5, 6)
(0.5, 6)
Early stopping patience
10
10
10
Training epochs
40
40
40
Batch size
512
512
512
Total #params
80k
80k
8k
are given in the Table 3. We used learning rate decay and early stopping using the Mean WAPE score on the validation set, with a patience of 10 for all models. We tuned the model hyper-parameters using the same metric. The various model hyper-parameters
given in Table 3.
All our experiments were implemented in Tensorflow 2, and run on Titan Xp GPU with 12GB of memory. The computing server used, had 256GB of memory, and 32 CPU cores, however, code did not seem to use more than 10GB of memory and CPÚ cores.
MiniiBatchingg: During each training iteration, sample minibatch of time series for computing the gradients. For constructing a minibatch, first time window (of length H + F) is uniformly randomly selected from the training time steps, after which subset of nodes sampled from the hierarchy tree. The time series data corresponding to this subset of nodes and the sampled time window constitutes our minibatch.
Predicting for new time steps: Once trained, our model can be used to predict for new datapoints without any retraining. However, in practice it may be beneficial to retrain the model with newer data to improve performance, even though this is not constraint for our proposed approach.
Knowledge of the hierarchy: Our proposed approach requires the hierarchy to be known during training be able to regularize the embeddings. In scenarios where aggregations presented during test time, reasonable prediction can be produced using the following strategy: the embedding of the new aggregation is set the mean of the embeddings of nodes in that aggregation the rest of the model remains the same. Predictions for unseen aggregations are made in the same way with seen aggregations. This idea beyond the scope of the current paper, thus left for future exploration.
B.4 FURTHER DETAILS ABOUT GLOBAL STATE Z
In many practical scenarios, the evolving global state of the set of time series may not be captured by the global covariates X only. For instance, when there is an overall increase/decrease in sales across all time series, required to feed in past values of time series to the BD model. However, since we cannot feed in the whole set of time series (order of 1000s) without leading to scalability issues, we choose small set of representative time series using Non-negative Matrix Factorization, thus approximating the global state. NMF assumes that the columns of the time series matrix Y E IR T lies in the convex set spanned by small subset of columns. However, this assumption may not hold true for many datasets, and therefore, most NMF algorithms compute an approximate factorization. We use this aforementioned subset of columns approximation the global state. We would also like to
captured in the past values Y rather than X. As result, may be
18

---

Under review as a conference paper at ICLR 2022
emphasize that our prediction model only sees the past values Zn input since the future time series values are unknown.
While NMF is one of the choices for Z, it is definitely not the only choice. Another option, PCA (or cquivalently SVD), may lead to latent vectors which do not have any temporal dependencies thus requiring additional temporal regularization (Sen al., 2019; Wilson et al., 2008). One may also use more sophisticated methods such as Temporal Latent Auto-Encoders (Nguyen & Quanz, 2021). We leave this idea for future exploration.
B.5 BASELINES
Further details about the baselines we compare with are provided as below. For a fair comparison, we use the the same Mean Absolute Error loss functions for training, and ensure that all the baseline models have similar number of parameters.
1. RNN: We use an LSTM decoder and encoder to implement a seq-2-seq model shared across all the time series, trained using mean absolute error loss.
2. DeepGLO (Sen al., 2019): We the implementation released by the authors Github We modify the loss function, data handling and evaluation to adapt to our setting.
3. DCRNN (Li et al., 2017): DCRNN being a GNN based approach requires a correlation graph imput. We use the official implementation released by the authors' and provide the hierarchy tree as the input graph. The implementation uses MAE loss by default.
4. Deep Factors (DF) (Wang et al., 2019): The original implementation released by the authors makes rolling probabilistic forecasts. We implement our own version in Tensorflow using an LSTT encoder for the global model (producing point predictions) as in the original implementation, while leaving out the probabilistic local model. We manually tune the hyper-paramenn for each of the datasets on the validation set.
5. EnE (Gleason, 2020): We implement this model using an LSTM decoder and encoder with MAE as the main training loss. In addition, we also use node embeddings which are fed as input to the encoder and decoder, and regularized according to the hierarchy as described by Gleason (2020).
6. SHARQ (Han et al., 2021): We were not able to find an official release of SHARQ. We implemented it using an LSTM based seq-to-seq model, with layer-wise training as described in the paper. We used MAE as the data fit loss function and the default squared error reconciliation loss regularizer (See Han et al. (2021) for terminology).
7. HierE2E (Rangapuram al., 2021): The official implementation produces rolling proba- bilistic forecasts. We adapt it to point forecasts by using the proposed projection step on outputs from a seq-2-seq model.
8. RNN+ERM (Ben Taieb & Koo, 2019): This approach involves learning a sparse projection matrix from data, resulting in coherent predictions. We tune the sparsity parameter using the model accuracy on the validation
B.6 MORE RESULTS
Table 4 show the test metrics averaged over 10 independent runs on the three datasets along with the standard deviations.
o
"https://ge http::/gg.gtb....mmoooogggg/////
19

---

Under review as a conference paper at ICLR 2022
Table 4: WAPE/SMAPE test metrics for all the three datasets, averaged over 10 runs. The standard deviations are shown in the parenthesis. We bold the smallest mean in each column and anything that comes within two standard deviations.
M5
Level
Level
Level 2
Tevel 3
Mesn
HIRED
(0.0011)
0.048
0.088 (0.0006)
(0.0006)
(0.0007)
0-072
(0.0000)
000077
0.279
(0.0012)
(0.0005)
(0.0006)
RNN
(0.002)
0.059
(b.005)
0.059
(0.013)
0.085
(0.011)
0.083
(0.002)
0.085
(0.004)
0..09
(0.006)
0.282
(0.007)
0.517
(0.008)
0.127
(0.005)
0.189
DF
(0.001)
0.055
(0.001)
0.056
(0.001)
0.061
/ 0.080
(0.001)
(0.001)
0.076
(0.002)
0.085
(0.000)
0.272
(0.002)
0. 0.801
(0.001)
0.116
(0.001)
0.176
DeepGLO
(0.0003)
(0.0004)
0.081
(0.0003)
(0.0004)
0.092
(0.0003)
0.099
(0.0005)
(0.0001)
(0.0001)
(0.0008)
(0.0003)
DURKN
(0.008)
0.078
(0:007)
0.079
(0:005)
0.086
(0.001)
0. 092
(0.000)
0.166
(0-007)
0.188
(0.000)
0.282
(0.000)
0-512
(0.002)
0.156
(0.003)
0-219
42Emb
(0.0018)
0.035
(0.001)
0,056
(0.0014)
0.084
(0.001)
0.083
(0.0011)
0.080
(0.001)
0.092
(0:0603)
0:269
(0.005)
0.501
(a (0.0009)
0.117 / (0.001)
0.178
SHARQ
(0.002)
093
(0.002)
0.086
(0.004)
0.071 / 0 0.2
(0.003)
(0.002)
0.090
(0.001)
0.094
(0.000)
0.277
(0.000)
0.528
0.135 (0.001)
(0.001)
195
RNN-EKM
(0.001)
0.052
(0.001)
0.052
(0.001)
D-066 / 0-071
(0 (0.002)
(0.00 001)
D 0.084
(0.002)
0. 194
(0.002)
0. 286
(0-004)
0.522
(0.001)
0:122 0. 0.187
(0.001)
Hicr-RZE
(0.003)
0.152
/ 0.100
(0:002)
(0:002)
0.152
(0.002)
0.158
(0 003)
152
(0-002)
0.181
(0.001)
0.386
(0-002)
0.611
(0.002)
0:213
(0. 003)
275
Hicr-FFF Largo
(0.002)
(0.003)
(0.001)
001)
001)
(0.001)
(0,001)
(0.001)
(0.001)
0.130
Favorita
Level
Level
Level
Level
Mzon
HIRID
(0.002) (0. 00.0
(0. (0.001)
(0.002)
(010))
(0.003)
0000
(0-004)
(0.001)
(0.009)
0.087
9. pes
0.111 0.197
0.133
0.290
0.203
0.339
0.130
0.222
(0.004)
(0.008)
(0.003) (0.004)
(0 (0000))
(0:305)
(0. 00100)
(0.003)
(0.002)
(0.004)
DP
(0.003)
0.064
(0.004)
0.064
[0. (0.002)
0.110
(0.003)
194
(0 002)
0. 135
(0:007)
251
(0. 001)
0.213
(0:007)
0.333
(0.003)
0.130
(0.004)
0.223
DeepGi
(0.001)
0.086
(0.001)
0.008
(0.001) (0.001)
0.126
0.197
(0.001) (0.001)
0.156 0.338
(0.001)
222
(0.001)
0.404
(0.001)
0.161
(0.001)
0.256
IXCNNN
(0.004)
(0.005)
(0.001) (0.002)
000)
(0.000)
(0. 000)
(0-000)
(0.001)
(6.001)
F-gh
(0:003)
0.070
(0.008)
0.070
(0.002)
0.114
(0.004)
199
0 136
(010)
(0:006)
0.276
(0. 0010
0. 207
(0:007)
0.891
(0. (0.003)
0. 132
(0.004)
0.216
SHARQ
(0.002)
0.0ss
(0.002)
0.085
(0.001)
0.142
(0.001)
199
(0.001)
0. 158
(0.001)
335
(0.000)
0. 0.280
(0.000)
0.404
(0,000)
0.164
(0.000)
0.256
RNN-ERM
(0.002)
0. 0.058
(0.002)
0.058
(0.001)
0.108
(0.905)
0.185
(0. (0.001)
0.129
(0:085)
0.283
(0.001)
0,220
(0.008)
0.348
(0.001)
0.127
(0.008)
0.219
Hier B2B
(0.005)
0.320
(0.008)
0.125
(0.003)
0.208
(0.005)
0.334
(0 (0.002)
0 247
(0.006)
0.448
(0.007)
4.9
(0:014)
0.573
(0.004)
0.245
(0.007)
0.370
Hier E2E Large
(0.002)
0.052
(0.002)
0.077
(0.008)
0.105
(0.010)
0.268
(0.002)
0.190
(0.003)
0,360
(0.002)
0.314
(0.001)
0.440
(0.002)
0.189
(0. (0.008)
0.256
Tourism
Level
Level
Level
Level
Level 4
Meu
IIRED
0.059
0010
(0-001)
0. 0.081
(0.001)
0 0.125
(0.003)
0.102
(0:001)
0.172
[0 002)
222
(0:001)
0.229
0.378
00.0
(0:001)
0.847
(0. (0.007)
0.786
(0. (0.001)
0.188
(0-002)
0.322
RNN
(0. (0.001)
0.110
(0.001)
0.106
(0.001)
0.148
(0.002)
0.181
(0:001)
0.188
(0.001)
0.231
(0.000)
0.240
(0.005)
0.385
(0.001)
0-589
(0.012)
0.782
(0.001)
0:211
(0.002)
0.333
DF
(0.003)
0.097
(0-003)
0.086
(0.002)
0.141
(0.003)
0.170
(0.001)
0-187
(0.002)
0.240
(0-001)
0.201 / 3-380
(0.002)
(0-000)
0-355
(0:019)
0.783
(0.001)
0.204
(0-003)
0.334
DeepGLO
(0.0002)
(0.0002)
(0.0001)
(0.0001)
(0 (0.00033
(0 0003))
(0 (0.0001)
(0.0001)
(0.0001)
(0.0002)
0.500
(0.0001)
(0.0001)
DCKNN
(0.003)
0. 0.1887
(0-003)
0-171
(0:002)
0.231
(0.003)
0.248
(0:001)
0.258
[0
278
(0-001)
0-283
(0.001)
0 398
(0.000)
0.434
(0.000)
0. 868
(0.000)
0.281
(0-001)
0-392
F.gEinh
(0. (0.007)
0.114
(0:007)
0.115
(0:002)
0.153
(0.004)
0.180
(0:002)
0.192
[0 002)
0.244
(0:001)
0.245 / 0.388
((0.002)
(0 (0.002)
0.372
[0. (1.010)
0.7a9
(0 (0.0
0:215
(0. (0.003)
0.342
SHARQ
(0.005)
0.100
(0-005)
0.104
(0.002)
0.164
(0.001)
0.208
(0.003)
0.217
(0.002)
250
(0-003) (0.001))
0.265 / 3.8
(0-003)
0-388
(0.005)
0.931
(0.00)
0.228
(0.001)
0.378
RNN+ERM
(0.005)
0.078
(0.006)
0.079
(0.003)
0.165
(0.006)
0.206
(0.004)
0.225
(0.006)
0.291
(0.006)
0.307 / 0.498
[0.008)
(0.009)
0.488
(000100
110113
(0.006)
0.251
(0.006)
0.417
Hies-FR
(0.002)
0.110
(0.002)
0.113
(0.002)
0.143
(0.003)
0.161
(0.002)
0.187
(0.003)
0.232
(0.001)
0.240
[0. (0.004)
0.371
(0.001)
0.358
[0. (0.003)
0.882
(0.001)
0.208
(0.002)
34340
20


Reviews:
[{"review_id": "odZSbvRh9uc", "reviewer": "Reviewer_UnFM", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper’s subject, forecasting, is being studied regularly in academia and \nresearch groups, and hierarchical forecasting specifically, has many use \ncases in industry.\nThe paper is well written and organized overall and conveys a good “reader \njourney” experience. The figures and tables describing and emphasizing \nrelevant parts and have good explanations.\nThe ‘Introduction’ describing the problem, gives relevant use cases and \nserves as a great entrance gate for the paper. Figure1 gives an important \nvisual intuition. The ‘Problem Statement’ paragraph is a must have in every \npaper and helps to prepare the reader for the ‘Problem Setting’ section and \nthe suggested model section. For a weakness, I would like to point out that \nthe paper is dealing with ‘additive coherence constraints’ only, while one of \nthe use cases described in the first paragraph is financial predictions, that \nmight not have the additive constraint (for example in financial markets), but \na different form of aggregation. It’s fine that the paper doesn’t deal with it, \nbut it would be worth if the paper touched on this point.\nIn section 4, the ‘HiReD’ model is being very well explained and Figure2 \ncompletes the picture of the model perfectly. For the model itself, the \nstrength of it is the fact that it uses simple building blocks that can be \nimplemented easily with state-of-the-art gradient based libraries, also it can \ndeal with hierarchical structured data and scales linearly with more time \nseries. Moreover, it uses information from the entire hierarchy in training but \nnot in the forecasting phase, that should result in faster forecasts. Its \nweakness is the fact that you must know the hierarchy structure in advance \n(it’s OK for the use cases described in the paper, but it might not always be \nthe case). The predictions are point forecasts and not probabilistic, while in \nmany domains confidence level or intervals are very important for decision \nmaking (like in retail). Another issue that wasn’t discussed is the scalability of\nthe model to new time steps – do you have to train the model all over again \nwith all history plus the new time step to produce an improved prediction (for\nt+2)?\n\nIn the ‘Experiments’ section, the method is compared to a variety of \nmethods, that can be divided into: (i) neural network (NN) models for time \nseries - that are not necessarily optimized for hierarchical problems; (ii) graph\nNN models; (iii) NN for hierarchical models; (iv) reconciliation models. The \nexplanations, metrics and result tables are well written, but I want to point \nout that two methods that were introduced in the ‘Related Work on Deep \nHierarchical Models’ part, Han et al. (2021) and Yanchenko et al. (2021), were\nnot  used as baseline. As actionable feedback, I recommend including in\nthe baseline more NN for hierarchical models methods than “regular” NN \ntime series methods, as the main subject of the paper is hierarchical time \nseries forecasting.", "summary_of_the_paper": "The paper introduces a method for hierarchical time series forecasting.\nThe problem setting is: given historical hierarchical univariate time series \ndata and given historical and future features (like holidays, etc..), try to \npredict future values for all time series, while keeping the coherence \nconstraints of the hierarchy.\nThe forecasting is broken down into an autoregressive part with shared \nparameters for all time series, and a basis decomposition with different \nweights for each time series (regularized by the hierarchy structure).\nThe method is then being compared to different baseline methods on three \ndatasets.", "summary_of_the_review": "The paper is well written and gives great introduction for the subject. The \nsuggested method’s strength is in its easy to implement building blocks, the \nability to take as into account information from all hierarchy and to scale well \nfor larger hierarchies, though it could be better if the issue of dealing with \nnew time steps would be addressed. The datasets chosen for the \nexperiments are showing good results against relevant baseline, although \nthere are few more relevant methods that were not included in the baseline.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "2FskTIoi5MK", "reviewer": "Reviewer_iRbQ", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n- Dimension reduction is an important topic that has been overlooked in the literature on hierarchical forecasting.\n- Experiments with several publicly-available real-world datasets.\n\nWeaknesses:\n- Paper clarity should be improved. Reading the paper raised a lot of (unanswered) questions. See comments below.\n- About the proposed method:\n* The proposed method produces roughly/approximate coherent forecasts. It is hard to compare methods that produce coherent and incoherent forecasts.\n* The paper only considers point forecasting. Probabilistic forecasting is a more challenging and important problem.\n* It is not clear whether the basis decomposition should be applied on all series. A hierarchical time series is a multivariate time series with hierarchical aggregation constraints. As a result, some series are linear combinations of other series. Isn't it enough to just apply it to the bottom-level series? \n- I am wondering if the comparison between methods is fair. In fact, the proposed method minimizes (regularized) MAE while other baselines minimize MSE. Your method is more aligned with the evaluation metrics, WAPE/SMAPE, which are essentially scaled absolute errors.\n\n\nOther comments:\n- Related work: literature on hierarchical probabilistic forecasting should be discussed, as well as dimension reduction methods for multivariate time series.\n- \"Hierarchical coherency\" has nothing to do with coherence as discussed in Thomson et al., 2019.\n- Section 2 should be split in two different paragraphs. One on mean forecasting and the other on probabilistic/quantile forecasting. You cannot really compare these two types of methods.\n- Figure 4: it is not clear why Tourism has a non-zero L3 colum.\n- What does \"asymptotically better\" means? Please be more rigorous or avoid this term.\n- \"under the reasonable assumption of the parent embedding being close to all the children embeddings\" -> Why is it a reasonable assumption?\n- Why imposing the hierarchical constraints on the embedding is equivalent to constraining the forecasts? This should be clearly explained.\n- The motivation for NMF should be clearly stated in the paper.\n- More details on the baselines should be given, at least in the appendix.\n- What does minibatch means with hierarchical time series data? This should be clearly explained.\n- Given the simplified assumption, I am not sure the theoretical analysis really gives additional insight. Also, it only takes one paragraph in the whole paper. \n- While simplified assumption are often required for a theoretical analysis, these assumptions shoul be clearly stated. \n- \"A small beta implies that the children time-series have structural similarities which is common in hierarchical datasets.\" -> I am not convinced by this sentence. Please give more details on Beta.\n- Is Lemma 1 really useful? Replacing values in the Theorem should be enough.\n- Simulation experiments could be used to confirm the usefulness of the method when indeed the data is low-dimensional. This would complement the theoretical analysis.\n", "summary_of_the_paper": "The paper considers point forecasting of hierarchical time series, i.e. multivariate time series with hierarchical aggregation constraints. The authors propose a new approach based on decomposing the series along a global set of basis time series where (approximate) hierarchical constraints are applied on the coefficients of the basis decomposition. Forecasts are produced using a dynamic linear autoregressive model. Compared to existing state-of-the-art hierarchical models, the proposed approach improved overall performance on forecasts at different levels of the hierarchy on several public datasets.\n", "summary_of_the_review": "While dimension reduction in hierarchical forecasting is an important problem, the paper lacks clarity and justifications for all the design choices. Major references are missing (dimension reduction for multivariate time series, probabilistic hierarchical forecasting, etc). The experimental setup is also questionable (loss function, etc). As a result, I do not think the paper is ready for publication.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "1iL6LE-3FGR", "reviewer": "Reviewer_U9mT", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\n- The assumption $|| \\theta_n - \\theta_0 ||_2 \\leq \\beta$ seems to be restrictive. The embedding parameters could satisfy the hierarchical structure while the parent embedding parameter $\\theta_0$ would be different that the children embedding $\\theta_n$. Authors mention at the top of page 7 that as $\\beta \\rightarrow 0$, the hierarchically regularized estimator approaches an error $L$ times smaller when compared to the un-regularized estimator. That seems to be more related to the assumption than the developed algorithm. For very small $\\beta$, essentially, the sample size is increased from $T$ to $LT$ since $L$ additional time series are observed, but thanks to the assumption that $\\beta$ is very small, the number of parameters to be estimated did not increase. As a result, the consistency rate for estimating the embedding parameter will be of order $O(1/(LT))$ instead of $O(1/(T))$. Thus, the theoretical justification of the proposed methodology is unclear.\n\n- Results of Theorem 1 are related to prediction error consistency rather than parameter estimation consistency. It is appropriate to consider verifying the rate for $|| \\tilde{\\theta}_n - \\theta_n ||_2$ instead of  $|| \\tilde{\\theta}_n - \\theta_n ||_\\Sigma$.\n\n- It is not clear why the Non-Negative Matrix Factorization (NMF) algorithm is used to select a small set\nof representative time series that encode the global state, please elaborate further.\n\n- Typo in page 7: \"Tables 1 and 2 shows ...\"\n\n- Typo in page 9: \"... state of the art baselines for most levels of the hierarchically.\"", "summary_of_the_paper": "In this manuscript, the authors proposed a method for hierarchical time series forecasting, consisting of two\ncomponents, the TVAR model, and the BD model. Further, the performance of the proposed algorithm is empirically evaluated on three benchmark datasets and showed that the proposed model consistently improved over state of the art baselines.", "summary_of_the_review": "\n- The paper is well written and easy to read.\n- Theoretical results need clarification, while some assumptions need to be relaxed.\n- Code for reproducing the results is provided together with links to publicly available real data sets used in the manuscript.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a method for time series forecasting based on a hierarchical deep learning approach. Three reviewers submitted reviews, with two marginally accept and one marginally reject. The paper was therefore borderline, but the issues raised by the marginal reject reviewer on the justification for the design choice of a deep latent model and the experimental setup appear worth addressing in a revision resubmitted to another conference.
Decision: Reject