# Hierarchical Forecasting at Scale

Olivier Sprangers<EMAIL_ADDRESS>Wander Wadman<EMAIL_ADDRESS>Sebastian Schelter<EMAIL_ADDRESS>Maarten de Rijke<EMAIL_ADDRESS>AIRLab, University of Amsterdam, Science Park 900, 1098 XH Amsterdam, The
Netherlands  bol.com, Papendorpseweg 100, 3528 BJ Utrecht, The Netherlands
University of Amsterdam, Science Park 900, 1098 XH Amsterdam, The Netherlands

###### Abstract

Hierarchical forecasting techniques allow for the creation of forecasts that
are coherent with respect to a pre-specified hierarchy of the underlying time
series. This targets a key problem in e-commerce, where we often find millions
of products across many product hierarchies, and forecasts need to be made for
both individual products and product aggregations. However, existing
hierarchical forecasting techniques scale poorly when the number of time
series increases, which limits their applicability at a scale of millions of
products.

In this paper, we propose to learn a coherent forecast for millions of
products with a single bottom-level forecast model by using a loss function
that directly optimizes the hierarchical product structure. We implement our
loss function using sparse linear algebra, such that the number of operations
in our loss function scales quadratically rather than cubically with the
number of products and levels in the hierarchical structure. The benefit of
our sparse hierarchical loss function is that it provides practitioners a
method of producing bottom-level forecasts that are coherent to any chosen
cross-sectional or temporal hierarchy. In addition, removing the need for a
post-processing step as required in traditional hierarchical forecasting
techniques reduces the computational cost of the prediction phase in the
forecasting pipeline, as well as its deployment complexity.

In our tests on the public M5 dataset, our sparse hierarchical loss function
performs up to 10% better as measured by RMSE and MAE compared to the baseline
loss function. Next, we implement our sparse hierarchical loss function within
an existing gradient boosting-based forecasting model at bol, a large European
e-commerce platform. At bol, each day a forecast for the weekly demand of
every product for the next twelve weeks is required. In this setting our
sparse hierarchical loss resulted in an improved forecasting performance as
measured by RMSE of about 2% at the product level, as compared to the the
baseline model, and an improvement of about 10% at the product level as
measured by MAE. Finally, we found an increase in forecasting performance of
about 5–10% (both RMSE and MAE) when evaluating the forecasting performance
across the cross-sectional hierarchies that we defined. These results
demonstrate the usefulness of our sparse hierarchical loss applied to a
production forecasting system at a major e-commerce platform.

###### keywords:

Hierarchical forecasting, Large-scale forecasting, Efficiency in forecasting
methods

##  1 Introduction

In e-commerce, we are often faced with two forecasting challenges. First,
forecasts at the lowest granularity - often the individual product level - are
required but we also need forecasts at higher granularities, for example at
the category, department or regional level, as higher level forecasts are
often needed in logistics and financial planning. Second, forecasts at
different time granularities are required, for example daily or weekly
forecasts. It is common that separate forecast models are made for each
separate (temporal) granularity, and as such these forecasts may not be
coherent with each other. Hierarchical forecasting [1] and temporal
hierarchical forecasting techniques [2, 3, 4] aim to solve the problem of
creating forecasts that are coherent with respect to a pre-specified cross-
sectional and/or temporal hierarchy of the underlying time series.

#### Challenges with existing cross-sectional and temporal hierarchical
forecasting techniques

Reconciliation methods adjust the forecasts for each level in the hierarchy by
minimizing the errors at each forecast level. These methods are applied as a
post-processing step that requires a matrix inversion that scales cubically
with the number of products or product hierarchies [1, 2, 5]. In settings with
millions of products such as in e-commerce, this becomes computationally
expensive at prediction time. Neural network methods can optimize for the
hierarchy in an end-to-end manner, however, these are either multivariate
methods that scale poorly to millions of time series [6] or they can only
optimize for the temporal hierarchy [3].

#### Sparse loss function

In order to overcome these scaling issues, we design a sparse Hierarchical
Loss (HL) function that directly optimizes both cross-sectional and temporal
hierarchical structures. Our corresponding sparsity-aware implementation
ensures that the number of operations in our loss function scales
quadratically rather than cubically with the number of products and levels in
the hierarchical structure, enabling computationally efficient training. The
benefit of our sparse hierarchical loss function is that it provides
practitioners a method of producing bottom-level forecasts that are coherent
to any chosen cross-sectional and temporal hierarchy. In addition, removing
the need for a post-processing step as used in traditional hierarchical
forecasting techniques reduces the computational cost of the prediction phase
in the forecasting pipeline. Furthermore, this also reduces the deployment
complexity of the forecasting pipeline.

#### Evaluation

We evaluate our sparse HL function on a gradient-boosted forecasting system on
the public M5 dataset [7] and a proprietary dataset from our e-commerce
partner. For the M5 dataset, we demonstrate that our implementation provides
up to 10% better forecasting performance as measured by both RMSE and MAE
compared with (i) reconciliation methods and (ii) baseline bottom-level
forecasting methods that use a standard loss function. For the proprietary
dataset, we present the results of an offline test on the product-level
forecast system of bol, a European e-commerce company with a catalog of
millions of unique products. We find our sparse HL function improves the
forecasting performance by about 2% on RMSE and 10% on MAE as compared to the
baseline forecasting system. This demonstrates the usefulness of our sparse HL
function in a large-scale setting.

#### Contributions

In summary, the main contributions of this paper are:

  1. 1.

We design a sparse hierarchical loss function that enables direct end-to-end
training of cross-sectional and temporal hierarchical forecasts in large-scale
settings in Section 4.

  2. 2.

We empirically demonstrate that our sparse hierarchical loss function can
outperform existing hierarchical forecasting reconciliation methods by up to
10% in Section 5.1.

  3. 3.

We show how our sparse hierarchical loss function scales to large-scale
settings and demonstrate a reduction of both training and prediction time of
up to an order of magnitude compared to the best hierarchical forecasting
reconciliation methods (Section 5.1).

  4. 4.

We present the results of an offline test of our method for the primary
product demand forecasting model at bol, a European e-commerce company with a
catalogue of millions of unique products, demonstrating an improvement of 2%
on RMSE and 10% on MAE as compared to the baseline forecasting system, in
Section 5.2.

##  2 Related work

#### Forecasting for large-scale settings

Contemporary large-scale forecasting applications require forecasting many
time series concurrently [8]. In academia, there has been a surge in the use
of neural network-based forecasting methods, which are methods that commonly
learn a single forecast model that can produce forecasts for many time series.
We refer the interested reader to the recent survey of Benidis et al. [9] for
an overview of these methods. However, tree-based methods topped the M5
forecasting competition [7], which is believed to be due to the strong
implementations available of these algorithms [10], such as the LightGBM [11]
or XGBoost [12] packages. Our own experience within bol confirms this view:
the ease of use, execution speed and strong default performance are key
reasons a tree-based method is often the default choice when creating a new
forecasting model.

#### Hierarchical forecasting

Hierarchical forecasting [1, 13, 14, 15, 5] and temporal hierarchical
forecasting techniques [16, 2, 3, 4] aim to solve the problem of creating
forecasts that are coherent with respect to a pre-specified cross-sectional
and/or temporal hierarchy of the underlying time series. We divide
hierarchical forecasting methods into Reconciliation methods and Other
methods.

Reconciliation methods. These methods solve the hierarchical forecasting
problem as a post-processing step by reconciling the forecasts to a pre-
specified cross-sectional and/or temporal hierarchy [1, 13, 14, 15, 5, 17,
18]. Limitations of these approaches are (i) that they require a post-
processing step, (ii) computing the reconciliation may be computationally
expensive, as we show in Section 3.1, and (iii) approaches that are
computationally less expensive tend to perform worse, as we show in Section 5.
Recent work by Taieb [16], Ben Taieb and Koo [15] has improved forecasting
performance of previous reconciliation approaches but at the expense of even
higher computational costs, as we explain in Section 3.

Other methods. In [6, 3] neural network-based end-to-end hierarchical
probabilistic forecasting method are proposed to solve the hierarchical
forecasting problem. More recently and most closely related to our work, Han
et al. [19] introduced SHARQ, a method that reconciles probabilistic
hierarchical forecasts during training by employing a regularized loss
function that aims to improve hierarchical consistency of bottom-up forecasts
through regularization. However, the regularization does not strictly enforce
the cross-sectional hierarchy in this method.

##  3 Background

To understand our problem setting and the issues we identify with existing
hierarchical forecasting methods, we introduce the hierarchical forecasting
problem and common methods of solving the hierarchical forecasting problem.

Suppose we have a time series $\vec{y}_{t}$, where $t$ denotes the time stamp.
We are interested in finding estimates $\hat{y}_{t}$ of future values
${y}_{t}$ of the time series by employing a model $f$ based on past values
$y_{t-1},\dots,y_{t-T}$ of the time series and additional attributes $X$:

$\hat{y}_{t}=f(y_{t-1},\dots,y_{t-T},X_{t},X_{t-1},\dots,X_{t-T}).$ (1)

In our hierarchical forecasting setting, we aim to create forecasts for many
time series concurrently, whilst adhering to pre-specified hierarchical
relationships that exist between the time series. This can be formalized as
[20]:

$\tilde{\textbf{y}}_{t}=SP\hat{\textbf{y}}_{t}\;,$ (2)

where $\hat{\textbf{y}}_{t}\in\mathbb{R}^{m}$ denotes the vector of forecasts
for all $m$ time series in the hierarchy, $S\in\\{0,1\\}^{m\times n}$ is a
matrix that defines the hierarchical relationship between the $n$ bottom-level
time series and the $m^{*}=m-n$ aggregations, $P\in\mathbb{R}^{n\times m}$ is
a matrix that encapsulates the contribution of each forecast to the final
estimate, and $\tilde{\textbf{y}}_{t}\in\mathbb{R}^{m}$ is the vector of
forecasts adjusted for the hierarchy. We can use the matrix $P$ to define
various forecast contribution scenarios. Note that we can straightforwardly
extend Equation (2) to the setting of temporal hierarchies [2, 3] by
considering forecasts of different time granularities in our vector of base
forecasts $\hat{\textbf{y}}_{t}$ and using an appropriate choice of $S$ to
aggregate series of a different time granularity. We will show how cross-
sectional and temporal hierarchical forecasting can be combined in Section 4.

The optimal solution to the problem in Equation (2) can be found using
Reconciliation methods and Other methods.

#### Reconciliation methods

MinTShrink [5] and variants find the optimal $P$ matrix by solving the
following minimization problem for a particular choice of $W$:

$\displaystyle\min_{P}$
$\displaystyle(\hat{\textbf{y}}_{t}-SP\hat{\textbf{y}}_{t})^{T}W(\hat{\textbf{y}}_{t}-SP\hat{\textbf{y}}_{t})$
$\displaystyle\text{s.t.}\quad PS=I.$ (3)

Assuming $W$ is positive definite, this has the following solution (ref.
Theorem 1 of [5]):

$\displaystyle P$ $\displaystyle=(S^{T}W^{-1}S)^{-1}S^{T}W^{-1}$
$\displaystyle=(J-JWU(U^{T}WU)^{-1}U^{T})\;,$ (4)

in which $S$ is partitioned as $S^{T}=[C^{T}\;I_{n}]$, $J=[0_{n\times
m^{*}}\;I_{n}]$, $U^{T}=[I_{m^{*}}\;-C]$. In MinTShrink, $W$ is estimated
using the shrunk empirical covariance estimate of [21]. Simpler choices for
$W$, such as the identity matrix, reduce the solution to the Ordinary Least
Squares (OLS) solution of [1]. In ERM, Ben Taieb and Koo [15] note than
MinTShrink and variants rely on the assumption of unbiasedness of the base
forecasts. They relax this assumption by formulating the hierarchical
reconciliation problem as an Empirical Risk Minimization problem, introducing
the ERM method. In addition, they propose two regularized variants of ERM
aimed at reducing forecast variance.

#### Other methods

Hier-E2E [6] solves the problem of Equation (2) by learning a neural network
model that combines the forecasting and reconciliation step in a single model,
resulting in an end-to-end solution removing the need for a post-processing
step. Similarly, COPDeepVAR [3] is an end-to-end neural network method that
enforces temporal hierarchies, however this is a univariate method that is not
able to enforce structural hierarchies (i.e., cross-sectional hierarchies)
simultaneously, and therefore not suited to our task. SHARQ [19] also moves
the reconciliation step into the training phase and achieves reconciliation
using a regularized loss function, where the regularization enforces the
coherency. However, this method does not enforce absolute coherency to the
hierarchy.

###  3.1 Scaling issues of hierarchical forecasting methods

Our main motivation for this paper are the limitations of prior work for
problem settings with many time series.

#### Scaling issues with reconciliation methods

In reconciliation methods, we encounter the following issues when scaling to
many time series:

  * •

The reconciliation is performed as a post-processing step, and thus has to be
performed as an additional step after generating the base forecasts. Even
though $P$ in Eq. 2 needs to be computed only once using Eq. (4), the
reconciliation still needs to be performed after each base forecast is
produced. Also, $P$ ideally is sparse [15], but no reconciliation method
guarantees this, so computing Eq. (2) will generally be a dense matrix-vector
product that scales with the number of time series.

  * •

For MinTShrink [5], estimating $W$ according to the method of [21] is
computationally expensive, with a computational complexity of $O(Nn^{2})$,
where $N$ denotes the number of training samples used to compute the shrunk
covariance estimate. In addition, the shrunk covariance estimate of [21] is
not guaranteed to give consistent results in high-dimensional settings [22],
making it less applicable for problem settings with many time series. Finally,
the estimate for $W$ will generally be a dense matrix, so we cannot make use
of efficient sparse algorithms to solve Eq. (4). However, even for simpler,
sparse choices of $W$ (such as the identity matrix of OLS [1]), we still need
to invert a matrix of size $m^{*}\times m^{*}$ in order to solve Eq. (4),
which becomes computationally costly for problems with many aggregations,
which naturally arise in retail forecasting scenarios. For example, for the M5
retail forecasting competition [23], $m^{*}=12,350$, even though there are
only 3,049 unique products in this dataset.

  * •

For ERM and its regularized variants [15], we need to either invert multiple
dense matrices that scale quadratically with the number of time series, or we
need to compute a Kronecker product that scales quadratically with the number
of time series, followed by an expensive lasso search procedure. Improving the
computational complexity of the ERM methods is also mentioned in [15] as an
avenue for future work.

#### Scaling issues with other methods

Hier-E2E [6] is a multivariate method, which means both input and output of
the neural network scale with the number of time series. For neural networks,
this significantly adds to the training cost and parameter cost as a large
amount of parameters are required to handle all the separate time series. This
in turn requires GPUs with more memory to train these models, which increases
cost to operate them.

##  4 Sparse Hierarchical Loss

In this section we present our main technical contribution, the sparse
hierarchical loss. First, we show how cross-sectional and temporal
hierarchical forecasting can be combined. Then, we introduce our loss function
and demonstrate it via a toy example.

#### Combining cross-sectional and temporal hierarchical forecasting

We are interested in finding forecasts that can be aggregated according to a
pre-specified cross-sectional hierarchy $S^{c}\in\\{0,1\\}^{m^{c}\times
n^{c}}$ and temporal hierarchy $S^{t}\in\\{0,1\\}^{m^{t}\times n^{t}}$:

$\displaystyle\tilde{\textbf{y}}^{c}_{t}$
$\displaystyle=S^{c}P^{c}\hat{\textbf{y}}^{c}_{t}\;,$ (5)
$\displaystyle\tilde{\textbf{y}}^{t}$
$\displaystyle=S^{t}P^{t}\hat{\textbf{y}}^{t}\;.$ (6)

These equations can be interpreted as follows:

  * •

In Eq. (5), we aggregate $n^{c}$ bottom-level time series from the same
timestep $t$ across a set of $m^{c}=n^{c}+m^{*,c}$ cross-sectional
aggregations.

  * •

In Eq. (6), we aggregate each time series consisting of $n^{t}$ timesteps
across a set of $m^{t}=n^{t}+m^{*,t}$ temporal aggregations, hence we drop the
subscript $t$.

We will only create bottom-level forecasts, thus $P^{c}=[0_{n^{c}}\times
m^{*,c}\;I_{n^{c}}]$ and $P^{t}=[0_{n^{t}}\times m^{*,t}\;I_{n^{t}}]$,
yielding:

$\displaystyle\tilde{\textbf{y}}^{c}_{t}$
$\displaystyle=S^{c}\hat{\textbf{y}}^{n^{c}}_{t}\;,$ (7)
$\displaystyle\tilde{\textbf{y}}^{t}_{t}$
$\displaystyle=S^{t}\hat{\textbf{y}}^{n^{t}}_{t}\;.$ (8)

Considering only bottom-level forecasts has a number of benefits: (i) each
forecast is coherent to any hierarchy by design, and (ii) we reduce the number
of required forecasts from $m$ to $n$, which can be a significant reduction
(there is no need for a forecast for $m^{*}$ aggregations in the hierarchy).
We now construct a matrix of bottom-level forecasts
$\hat{\textbf{Y}}\in\mathbb{R}^{n^{c}\times n^{t}}$, in which the columns
represent the forecasts of the bottom-level time series at a timestep $t$.
This allows us to combine (7) and (8) as follows:

$\tilde{\textbf{Y}}=S^{c}\hat{\textbf{Y}}(S^{t})^{\intercal}\;,$ (9)

in which $\tilde{\textbf{Y}}\in\mathbb{R}^{m^{c}\times m^{t}}$ represents the
matrix of forecasts aggregated according to both cross-sectional and temporal
hierarchies. Equivalently, we can aggregate our bottom-level ground truth
values $\bar{\textbf{Y}}\in\mathbb{R}^{n^{c}\times n^{t}}$:

$\textbf{Y}=S^{c}\bar{\textbf{Y}}(S^{t})^{\intercal}\;.$ (10)

#### Sparse hierarchical loss

To find the best forecasts for the hierarchical forecasting problem (9), we
try to find a model $f(\cdot)$ by gradient-based optimization of the following
loss function:

$\displaystyle L$
$\displaystyle=\sum\left[\vbox{\hbox{$\displaystyle\frac{\frac{1}{2}\left(\textbf{Y}-\tilde{\textbf{Y}}\right)^{2}}{\underbrace{\left[l^{c}\sum_{j=1}^{n^{c}}S^{c}_{ij}\right]}_{d^{c}}\otimes\underbrace{\left[l^{t}\sum_{j=1}^{n^{t}}S^{t}_{ij}\right]^{\intercal}}_{(d^{t})^{\intercal}}}$}}\right]\;,$
(11)
$\displaystyle=\sum\left[\vbox{\hbox{$\displaystyle\frac{\frac{1}{2}\left(\textbf{Y}-\tilde{\textbf{Y}}\right)^{2}}{d^{c}\otimes(d^{t})^{\intercal}}$}}\right]\;,$
(12)

in which $l^{c}$ and $l^{t}$ denote the number of levels in hierarchies
$S^{c}$ and $S^{t}$, respectively, and $i,\,j$ are the (row, column) indices
of the elements of the $S$-matrices, and $\otimes$ denotes the outer-product
between the denominator components in the loss function. Note that (11) shares
similarities with the Weighted Root Mean Squared Error from the M5 competition
[7].

We can derive the gradient and hessian of (11) with respect to the bottom-
level forecasts $\hat{\textbf{Y}}$:

$\displaystyle\frac{\partial L}{\partial\tilde{\textbf{Y}}}$
$\displaystyle=\left[\frac{\tilde{\textbf{Y}}-\textbf{Y}}{d^{c}\otimes(d^{t})^{\intercal}}\right]\;,$
(13) $\displaystyle\frac{\partial L}{\partial\hat{\textbf{Y}}}$
$\displaystyle=(S^{c})^{\intercal}\left(\frac{\partial
L}{\partial\tilde{\textbf{Y}}}\right)S^{t}\;,$ (14)
$\displaystyle\frac{\partial^{2}L}{\partial\hat{\textbf{Y}}^{2}}$
$\displaystyle=(S^{c})^{\intercal}\left(\frac{1}{d^{c}\otimes(d^{t})^{\intercal}}\right)S^{t}\;.$
(15)

#### Analysis

The best possible forecast is achieved when the loss (11) is minimized, or
equivalently when the gradient (13) is zero:

$\displaystyle\frac{\partial L}{\partial\tilde{\textbf{Y}}}$
$\displaystyle=\left[\frac{\tilde{\textbf{Y}}-\textbf{Y}}{d^{c}\otimes(d^{t})^{\intercal}}\right]\;,$
$\displaystyle=\left[\frac{S^{c}\hat{\textbf{Y}}(S^{t})^{\intercal}-S^{c}\bar{\textbf{Y}}(S^{t})^{\intercal}}{d^{c}\otimes(d^{t})^{\intercal}}\right]\;,$
$\displaystyle=\left[\frac{S^{c}\left(\hat{\textbf{Y}}-\bar{\textbf{Y}}\right)(S^{t})^{\intercal}}{d^{c}\otimes(d^{t})^{\intercal}}\right]\;,$

which becomes zero when $\hat{\textbf{Y}}=\bar{\textbf{Y}}$. Thus, the best
forecast model is found when each bottom-level forecast equals the ground
truth. This is equivalent to the standard (i.e., non-hierarchical) squared
error loss often used in forecasting problems. We argue that our hierarchical
loss gradient can be seen as a smoothed gradient compared to the standard
squared error loss gradient (i.e., $\hat{\textbf{Y}}-\bar{\textbf{Y}}$). For
example, consider the canonical case where we have two bottom-level time
series ($n^{c}=2$) consisting of two timesteps ($n^{t}=2$). Furthermore,
suppose we have a single cross-sectional aggregation (the sum of the two time
series, thus $m^{*,c}=1$ and $m^{c}=m^{*,c}+n^{c}=3$), and a single temporal
aggregation (the sum of the two timesteps, thus $m^{*,t}=1$ and
$m^{t}=m^{*,t}+n^{t}=3$). Finally, there are two levels in our cross-sectional
hierarchy and in our temporal hierarchy, thus $l^{c}=2$ and $l^{t}=2$. The
standard squared error loss gradient for this problem is:

$\displaystyle\begin{bmatrix}\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{0}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{1}}\\\ \frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}\end{bmatrix}=\begin{bmatrix}e_{0}^{0}&e_{0}^{1}\\\
e_{1}^{0}&e_{1}^{1}\\\ \end{bmatrix}\;,$ (16)

in which $e^{i}_{j}$ denotes the bottom-level error
($\hat{y}^{i}_{j}-\bar{y}^{i}_{j}$). For our hierarchical loss, Eq. (9) reads:

$\displaystyle\tilde{\textbf{Y}}=\underbrace{\begin{bmatrix}1&1\\\ 1&0\\\
0&1\end{bmatrix}}_{S^{c}}\underbrace{\begin{bmatrix}\hat{\textbf{y}}^{0}_{0}&\hat{\textbf{y}}^{0}_{1}\\\
\hat{\textbf{y}}^{1}_{0}&\hat{\textbf{y}}^{1}_{1}\\\
\end{bmatrix}}_{\hat{\textbf{Y}}}\underbrace{\begin{bmatrix}1&1&0\\\ 1&0&1\\\
\end{bmatrix}}_{(S^{t})^{\intercal}}\;,$ (17)

and the gradient of the loss with respect to the bottom level time series (14)
reads (ref. A for the full derivation):

$\displaystyle\begin{bmatrix}\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{0}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{1}}\\\ \frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}\end{bmatrix}=(S^{c})^{\intercal}\left[\frac{S^{c}\left(\hat{\textbf{Y}}-\bar{\textbf{Y}}\right)(S^{t})^{\intercal}}{d^{c}\otimes(d^{t})^{\intercal}}\right]S^{t}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}\frac{9}{16}e_{0}^{0}+\frac{3}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{1}+\frac{1}{16}e_{1}^{1}&\frac{9}{16}e_{0}^{1}+\frac{3}{16}e_{0}^{1}+\frac{3}{16}e_{1}^{1}+\frac{1}{16}e_{1}^{0}\\\
\frac{9}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{0}+\frac{3}{16}e_{1}^{1}+\frac{1}{16}e_{0}^{1}&\frac{9}{16}e_{1}^{1}+\frac{3}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{1}+\frac{1}{16}e_{0}^{0}\\\
\end{bmatrix}\;.$

When we compare this result to the standard squared error loss gradient (16),
we find that we smooth the bottom-level gradient by adding to it portions of
the gradients of all cross-sectional and temporal aggregations the bottom-
level series belongs to. This derivation also shows the motivation of adding
the denominator parts $d^{c}$ and $d^{t}$ to the loss function (11): it is
neccessary to scale the aggregation gradients by the number of elements in the
aggregation, otherwise the magnitude of the gradient grows with the number of
time series and the number of levels in the hierarchy, which we found to be
undesirable when trying to facilitate stable learning.

#### Sparsity

$S^{c}$ and $S^{t}$ are highly sparse. For example, $S^{c}$ has only
$n^{c}l^{c}$ non-zero elements: the number of bottom-level time series
multiplied by the number of aggregations in the hierarchy. Hence, the overall
sparsity of $S^{c}$ is given by $1-\frac{n^{c}l^{c}}{m^{c}n^{c}}$. For the M5
dataset [23], $n^{c}=3,049$, $l^{c}=12$, $m^{c}=42,840$, corresponding to a
sparsity of 99.97%. Next, the matrix of bottom-level ground truth values
$\bar{\textbf{Y}}$ in (10) may be sparse too, for example in the case of
products that are not on sale for every timestep $n^{t}$ in the dataset. All
these sources of sparsity motivate the use of sparse linear algebra when
computing Equations (11)–(15).

#### Implementation

We implement the hierachical loss (11), the bottom-level gradient (13), (14)
and hessian (15) in Python using the sparse library from SciPy [24]. Note that
equations (13)-(14) can be rearranged:

$\frac{\partial
L}{\partial\hat{\textbf{Y}}}=\frac{(S^{c})^{\intercal}}{d^{c}}\left(\tilde{\textbf{Y}}-\textbf{Y}\right)\frac{S^{t}}{(d^{t})^{\intercal}}\;,$
(18)

such that the parts before and after the brackets can be precomputed as they
do not depend on the forecast values $\tilde{\textbf{Y}}$, avoiding a costly
division operation inside a training iteration. Also note that the hessian
(15) does not depend on the forecast values $\tilde{\textbf{Y}}$, so it can be
precomputed as well. Our implementation, including the code to reproduce the
experiments on public data from Section 5, is available on
GitHub.111https://github.com/elephaint/hfas/

##  5 Experiments

In this section we empirically verify the usefulness of our sparse
hierarchical loss. First, we evaluate forecasting accuracy using a set of
experiments on the public M5 dataset [23]. Then, we evaluate our sparse
hierarchical loss in an offline experiment on a proprietary dataset from our
e-commerce partner.

###  5.1 Public datasets

#### Task & Dataset

Our task is to forecast product demand. We use the M5 dataset [23] for our
offline, public dataset experiments. The M5 dataset contains product-level
sales from Walmart for 3,049 products across 10 stores in the USA.
Furthermore, the dataset contains 12 cross-sectional product aggregations
(e.g. department, region), which allow us to test hierarchical forecasting
performance. We preprocess the dataset resulting in a set of features as
described in B. We forecast 28 days into the future.

#### Baseline models

For our baseline forecasting model, we primarily use LightGBM [11], trained to
predict one-day ahead. We subsequently recursively generate predictions for 28
days. Tree-based models dominated the M5 forecasting competition due to their
strong performance and ease of use [7, 10]. Moreover, our e-commerce partner’s
primary product forecasting is a LightGBM-based model, so we expect results
from offline experiments on public datasets to transfer to our proprietary
setting when using the same base forecasting model. We compare the performance
of our LightGBM models against traditional statistical methods ARIMA [25], ETS
[26], Theta [27], SeasonalNaive [20], Naive [20] and Croston [28]. We note
that deep learning based approaches are becoming more prevalent in e-commerce
[29], especially with the rise of the Transformer-architecture in forecasting
models [30, 31]. We consider this for future work, and did not consider this
for our study as (i) the cloud cost to operate these models is 10x higher for
our e-commerce partner as compared to a tree-based model, and (ii) none of the
neural-network based methods are able to scale to the size of our e-commerce
partner, as explained in Section 3.1.

#### Experimental setup

To test our hierarchical sparse loss function against baseline forecasting
systems, we consider the following scenarios:

  1. 1.

Bottom-up. We train a single global model on only the bottom-level time
series.

  2. 2.

Separate aggregations. We train separate models for every aggregation in the
hierarchy, resulting in 12 models for the entire M5 dataset.

  3. 3.

Global. We train a single global model on all time series in the dataset,
including all the aggregations.

For the first scenario in our experiments (Bottom-up), we vary both the
objective (i.e. loss function that is optimized by LightGBM) and the
evaluation metric (i.e. the loss function that governs early-stopping during
hyperparameter optimization). For the objective, we consider the squared error
loss (SL), the Tweedie loss (TL) and our sparse hierarchical loss (HL). The
Tweedie loss is a loss function that assumes that the time series follow a
distribution somewhere in between a Poisson and a Gamma distribution, which is
useful in zero-inflated settings such as retail demand forecasting. It is a
loss function that was favored by contestants in the M5 forecasting
competition [10], and it is the loss also used in the primary forecasting
system of our e-commerce partner.

For the latter two scenarios, we will obtain non-coherent forecasts. Thus,
these methods require a reconciliation post-processing step to reconcile the
forecasts to the hierarchy. We employ the following cross-sectional
reconciliation methods:

  * •

Base. No reconciliation is performed.

  * •

OLS. Ordinary Least Squares (OLS) [1], where $W$ in Equation (3) is the
identity matrix.

  * •

WLS-struct. and WLS-var. Weighted Least Squares (WLS) [5], where $W$ in
Equation (3) is a diagonal matrix containing respectively the sum of the rows
of $S$ (WLS-struct) or the in-sample forecast errors (WLS-var).

  * •

MinT-shrink. Trace Minimization [5], where $W$ in Equation (3) is the shrunk
covariance matrix of in-sample forecast errors. We also experimented with
using the non-shrunk covariance matrix of the in-sample forecast errors (MinT-
cov), but this produced erroneous/high variance results, which we attribute to
precisely the motivation to shrink the covariance matrix in MinT-shrink: to
reduce the variance when the amount of time series considered becomes very
large.

  * •

ERM. The Empirical Risk Minimization (ERM) method [15]. Due to computational
issues explained in Section 3.1, we were not able to apply the regularized ERM
variants to our experiments, but only the unregularized variant.

We optimize the hyperparameters of each of the LightGBM models by Bayesian
hyperparameter optimization using Optuna [32]. The settings for the
hyperparameter optimization can be found in C. Each model is trained for 10
different random seeds, and our results are based on the mean and standard
deviation of those 10 rollouts.

#### Evaluation

We evaluate our results for every aggregation in the hierarchy using the Root
Mean Squared Error (RMSE) and Mean Absolute Error (MAE) [20]. In the results
section, we present the RMSE / MAE relative to the Bottom-up scenario using
the squared-loss objective with the squared-loss metric. For absolute values
and standard deviation of the results, see D.

#### Results - LightGBM as baseline model

For our first experiment, we only consider cross-sectional hierarchies (i.e.
$S^{t}=I_{n^{t}}$). We present our results on relative RMSE using LightGBM as
baseline model in Table 4 and conclude the following:

  * •

The best method is the Bottom-up-scenario combined with our sparse
hierarchical loss as objective, outperforming the baseline by 0–20% across
aggregations. This holds for both settings in which we use our sparse
hierarchical loss.

  * •

Even when we only use our sparse hierarchical loss as an evaluation metric
during training whilst optimizing the standard squared loss (the SL/HL
scenario), we already see a small improvement of $\pm$5% across aggregations.

  * •

Even though the Tweedie loss improves over the baseline loss, our sparse
hierarchical loss function still outperforms it by $\pm$5% across
aggregations.

  * •

From the reconciliation methods, MinT-shrink and WLS-var perform best in the
Separate aggregations-scenario, although the performance delta across
aggregations is still $\pm$5-30% as compared to the best (our) method.

For relative MAE, we present our results in Table 5. We find that overall, our
sparse hierachical loss still performs best by $\pm$5% compared to other loss
functions and scenarios. However, the results are more nuanced: we find that
MinT-shrink in the Separate aggregations-scenario performs strong as well. In
addition, we also find that the Tweedie loss (TL) performs relatively well.
This finding corroborates the usefullness of the TL in intermittent demand
settings, such as retail, where zero demand is often observed.

Next, we compare our findings against the forecasting results when employing
different baseline models in Table 1. For brevity, we only show the metrics
for all time series combined (incl. aggregations). In addition, we only show a
single reconciliation method for the other baseline models, as we found little
difference in results when employing different reconciliation methods. We then
find that on RMSE, our sparse hierachical loss outperforms the other baseline
models by at least 50%, and on MAE by at least 10%. This verifies that on this
dataset and with this type of problem, using a more complex model such as
LightGBM greatly improves forecasting performance, as was also shown in the M5
forecasting competition [7].

Table 1: Forecasting results for all time series (incl. aggregations) on the M5 dataset, using different baseline models. We show absolute and relative RMSE and MAE. Lower is better, and bold indicates the best performing method. |  | RMSE | MAE
---|---|---|---
Model | Reconciliation | Abs. | Rel. | Abs. | Rel.
LightGBM (SL/SL) | None | 22.39 | 2.20 | 1.00 | 1.00
LightGBM (HL/HL) | None | 19.54 | 2.10 | 0.87 | 0.95
LightGBM (HL/SL) | None | 19.59 | 2.10 | 0.88 | 0.95
ARIMA | MinT-shrink | 39.88 | 2.43 | 1.78 | 1.10
ETS | MinT-shrink | 36.48 | 2.35 | 1.63 | 1.07
Theta | MinT-shrink | 36.66 | 2.39 | 1.64 | 1.08
Croston | None | 39.40 | 2.76 | 1.76 | 1.25
Naive | None | 74.91 | 3.95 | 3.35 | 1.80
Seasonal Naive | None | 39.40 | 2.76 | 1.76 | 1.25
Table 2: Forecasting results for all time series (incl. aggregations) on the M5 dataset, ablating for the use of cross-sectional and temporal hierarchies. We show absolute and relative RMSE and MAE, with the standard deviation in brackets. Lower is better, and bold indicates the best performing method. Note that when not using cross-sectional nor temporal aggregations, the hierarchical loss is equal to the standard squared error loss. Hierarchies | RMSE | MAE
---|---|---
Cross-sectional | Temporal | Abs. | Rel. | Abs. | Rel.
No | No |  22.39 (0.16) |  2.20 (0.01) | 1.00 | 1.00
Yes | No |  19.54 (0.38) |  2.10 (0.01) | 0.87 | 0.95
Yes | Yes | 29.81 (1.52) | 2.47 (0.04) | 1.33 | 1.12
No | Yes | 26.65 (0.32) | 2.36 (0.01) | 1.19 | 1.07
Random | No | 25.62 (4.75) | 2.31 (0.15) | 1.14 | 1.05
Table 3: Forecasting results for all time series (incl. aggregations) on the M5 dataset, ablating for the use of cross-sectional and temporal hierarchies. We show relative RMSE for several forecasting day buckets of the forecast. Lower is better, and bold indicates the best performing method. Hierarchies | Forecast day
---|---
Cross-sectional | Temporal | 1–7 | 8–14 | 15–21 | 22–28
No | No | 1.00 | 1.00 | 1.00 | 1.00
Yes | No | 0.98 | 0.99 | 0.81 | 0.80
Yes | Yes | 1.75 | 1.50 | 0.96 | 1.66
No | Yes | 1.37 | 1.24 | 0.97 | 1.41
Random | No | 1.25 | 1.16 | 1.00 | 1.32

Table 4: Forecasting results for all stores on the M5 dataset, using LightGBM as baseline model. We report relative RMSE as compared to the baseline (shown in italic). Lower is better, and bold indicates best method for the aggregation, taking into account standard deviation of the best method across the 10 seeds. For absolute values and standard deviation of the results, see D. |  |  |  |  |  | Store | Product | State |  |
---|---|---|---|---|---|---|---|---|---|---
Scenario/Objective | Metric | Reconciliation | Product | Department | Category | Department | Category | Total | Store | State | Department | Category | Total | Total | All series
Bottom-up |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | None | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
SL | HL | None | 1.00 | 0.97 | 0.95 | 0.98 | 0.97 | 0.98 | 0.99 | 1.00 | 0.98 | 0.96 | 0.99 | 1.00 | 0.98
HL | HL | None | 1.00 | 0.88 | 0.80 | 0.93 | 0.89 | 0.93 | 0.97 | 0.99 | 0.89 | 0.84 | 0.88 | 0.87 | 0.87
HL | SL | None | 1.00 | 0.88 | 0.81 | 0.94 | 0.90 | 0.94 | 0.96 | 0.98 | 0.89 | 0.84 | 0.88 | 0.87 | 0.88
TL | HL | None | 1.00 | 0.95 | 0.96 | 0.99 | 1.00 | 1.00 | 0.99 | 1.00 | 0.97 | 0.98 | 0.98 | 0.96 | 0.97
TL | SL | None | 1.00 | 0.94 | 0.93 | 1.00 | 1.00 | 1.00 | 0.99 | 1.00 | 0.97 | 0.98 | 0.98 | 0.93 | 0.96
TL | TL | None | 1.17 | 2.72 | 2.81 | 1.76 | 1.83 | 1.73 | 1.52 | 1.33 | 2.15 | 2.18 | 2.08 | 2.71 | 2.38
Sep. agg. |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.00 | 1.44 | 1.29 | 1.19 | 1.14 | 1.14 | 1.01 | 0.99 | 1.23 | 1.34 | 1.27 | 1.60 | 1.35
SL | SL | OLS | 1.00 | 1.39 | 1.41 | 1.10 | 1.06 | 1.07 | 1.00 | 1.00 | 1.20 | 1.19 | 1.23 | 1.50 | 1.30
SL | SL | WLS-struct | 1.00 | 1.26 | 1.37 | 1.03 | 1.05 | 1.02 | 0.99 | 0.99 | 1.11 | 1.16 | 1.16 | 1.39 | 1.23
SL | SL | WLS-var | 1.00 | 1.12 | 1.23 | 0.99 | 1.02 | 0.99 | 0.99 | 0.99 | 1.03 | 1.09 | 1.07 | 1.22 | 1.12
SL | SL | MinT-shrink | 1.00 | 1.15 | 1.27 | 0.97 | 0.99 | 0.97 | 1.00 | 1.00 | 1.03 | 1.09 | 1.09 | 1.30 | 1.15
SL | SL | ERM | 1.22 | 1.25 | 1.29 | 1.07 | 1.03 | 1.07 | 1.17 | 1.22 | 1.17 | 1.14 | 1.22 | 1.49 | 1.26
Global |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.02 | 1.33 | 1.45 | 1.09 | 1.10 | 1.10 | 1.03 | 1.03 | 1.25 | 1.27 | 1.81 | 1.57 | 1.46
SL | SL | OLS | 1.01 | 1.32 | 1.39 | 1.07 | 1.09 | 1.16 | 1.02 | 1.02 | 1.20 | 1.25 | 1.38 | 1.49 | 1.34
SL | SL | WLS-struct | 1.01 | 1.38 | 1.54 | 1.08 | 1.13 | 1.11 | 1.03 | 1.02 | 1.19 | 1.28 | 1.27 | 1.55 | 1.36
SL | SL | WLS-var | 1.01 | 1.51 | 1.70 | 1.18 | 1.27 | 1.22 | 1.03 | 1.02 | 1.31 | 1.43 | 1.38 | 1.66 | 1.48
SL | SL | MinT-shrink | 1.03 | 1.26 | 1.41 | 1.05 | 1.10 | 1.15 | 1.06 | 1.05 | 1.11 | 1.17 | 1.24 | 1.54 | 1.30
SL | SL | ERM | 1.21 | 1.59 | 1.69 | 1.26 | 1.28 | 1.34 | 1.20 | 1.23 | 1.45 | 1.49 | 1.61 | 1.80 | 1.59

Table 5: Forecasting results for all stores on the M5 dataset, using LightGBM as baseline model. We report relative MAE as compared to the baseline (shown in italic). Lower is better, and bold indicates best method for the aggregation, taking into account standard deviation of the best method across the 10 seeds. For absolute values and standard deviation of the results, see D. |  |  |  |  |  | Store | Product | State |  |
---|---|---|---|---|---|---|---|---|---|---
Scenario/Objective | Metric | Reconciliation | Product | Department | Category | Department | Category | Total | Store | State | Department | Category | Total | Total | All series
Bottom-up |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | None | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00 | 1.00
SL | HL | None | 1.00 | 0.98 | 0.96 | 0.98 | 0.98 | 0.99 | 1.00 | 1.00 | 0.97 | 0.97 | 1.00 | 1.02 | 0.99
HL | HL | None | 0.99 | 0.81 | 0.78 | 0.90 | 0.89 | 0.94 | 0.98 | 0.99 | 0.85 | 0.83 | 0.89 | 0.88 | 0.95
HL | SL | None | 0.99 | 0.81 | 0.77 | 0.90 | 0.90 | 0.94 | 0.98 | 0.99 | 0.85 | 0.83 | 0.89 | 0.87 | 0.95
TL | HL | None | 0.98 | 0.83 | 0.82 | 0.93 | 0.94 | 0.96 | 0.97 | 0.98 | 0.88 | 0.89 | 0.92 | 0.88 | 0.95
TL | SL | None | 0.99 | 0.85 | 0.84 | 0.95 | 0.95 | 0.98 | 0.98 | 0.99 | 0.90 | 0.91 | 0.96 | 0.90 | 0.96
TL | TL | None | 1.02 | 2.06 | 2.39 | 1.47 | 1.61 | 1.65 | 1.14 | 1.07 | 1.69 | 1.88 | 1.97 | 2.86 | 1.29
Sep. agg. |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.00 | 1.11 | 1.06 | 1.02 | 1.10 | 1.08 | 0.96 | 0.97 | 1.03 | 1.15 | 1.20 | 1.55 | 1.02
SL | SL | OLS | 0.96 | 1.08 | 1.16 | 1.01 | 0.99 | 1.00 | 0.96 | 0.97 | 1.00 | 1.01 | 1.12 | 1.49 | 0.99
SL | SL | WLS-struct | 0.97 | 0.99 | 1.13 | 0.92 | 0.95 | 0.95 | 0.96 | 0.97 | 0.93 | 0.98 | 1.04 | 1.42 | 0.98
SL | SL | WLS-var | 0.98 | 0.91 | 1.00 | 0.91 | 0.92 | 0.92 | 0.96 | 0.97 | 0.90 | 0.93 | 0.95 | 1.16 | 0.96
SL | SL | MinT-shrink | 0.97 | 0.93 | 1.05 | 0.90 | 0.90 | 0.91 | 0.96 | 0.97 | 0.89 | 0.92 | 0.98 | 1.30 | 0.96
SL | SL | ERM | 1.18 | 1.17 | 1.21 | 1.09 | 1.06 | 1.07 | 1.19 | 1.22 | 1.11 | 1.12 | 1.19 | 1.57 | 1.18
Global |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.04 | 1.04 | 1.17 | 1.00 | 1.02 | 1.05 | 0.99 | 0.99 | 1.04 | 1.09 | 1.62 | 1.58 | 1.05
SL | SL | OLS | 0.98 | 1.09 | 1.17 | 1.05 | 1.07 | 1.13 | 0.99 | 1.00 | 1.09 | 1.13 | 1.27 | 1.50 | 1.03
SL | SL | WLS-struct | 0.99 | 1.09 | 1.26 | 0.98 | 1.01 | 1.04 | 1.00 | 0.99 | 1.00 | 1.07 | 1.15 | 1.57 | 1.02
SL | SL | WLS-var | 0.99 | 1.24 | 1.40 | 1.09 | 1.13 | 1.13 | 1.01 | 1.00 | 1.14 | 1.21 | 1.25 | 1.69 | 1.06
SL | SL | MinT-shrink | 0.99 | 1.13 | 1.24 | 1.06 | 1.08 | 1.15 | 1.02 | 1.01 | 1.05 | 1.07 | 1.15 | 1.60 | 1.04
SL | SL | ERM | 1.17 | 1.49 | 1.62 | 1.30 | 1.34 | 1.37 | 1.19 | 1.20 | 1.39 | 1.49 | 1.58 | 1.96 | 1.27

#### Analysis: impact of hierarchy

We investigate the impact of the choice of hierarchy.

Temporal hierarchies. As we noted before, we only used cross-sectional
aggregations in our first experiments. We now also include temporal
aggregations by aggregating our bottom-level time series across years, weeks
and months. We ablate for every setting and show the results in Table 2.
Interestingly, we find that using temporal hierarchies jointly with cross-
sectional hierarchies reduces forecasting performance by $\pm$35% (RMSE) and
$\pm$17% (MAE). This setting is even worse than only using temporal
hierarchies, which performs worse than using only cross-sectional hierarchies
by $\pm$26% (RMSE) and $\pm$12% (MAE). We further analyze these results by
studying the RMSE across the forecast days in Table 3. As noted before, we
forecast 28 days ahead, and each forecast is created by recursively applying
the one-step ahead LightGBM model. We find that as we forecast further into
the future, the setting with only using cross-sectional aggregations starts to
perform better by up to 20% as compared to the baseline where we do not use
any aggregations. Again, the setting where we employ temporal hierarchies too
shows relatively bad performance across all forecast day buckets.

Random hierarchies. In hierarchical forecasting problems, the aggregation
matrices $S^{c}$ and $S^{t}$ are commonly fixed a priori and considered
constant during training and prediction. As we are performing the
reconciliation in an end-to-end fashion during training, we can modify these
matrices at every iteration. This allows us to understand the robustness of
our solution to possible misspecification of the hierarchy, and more
generally, to what extent the choice of the hierarchy has an effect on
forecasting performance. We perform an experiment by randomly sampling an
$S^{c}$-matrix at every iteration of the LightGBM training process. At every
iteration, we sample uniformly at random (i) a number of levels for the cross-
sectional hierarchy and (ii) a number of maximum categories for the level and
construct a random $S^{c}$-matrix to be used in the gradient (14) and hessian
(15). We validate and test on the ‘true’ $S^{c}$-matrix. We present the
results in Table 2 and Table 3, under ‘Random’. We find that performance
deteriorates by about 30% as compared to the baseline, and we also find that
our results in this experiment show higher variance compared to the baseline
methods. Thus, misspecification of the hierarchy can severely deteriorate
forecasting performance.

#### Analysis: time complexity

We investigate the computational time complexity required to perform training
and prediction for each scenario and present the results in Table 6. The
training and prediction time complexity is indicated by how respectively the
training time and prediction time scales with respect to the default LightGBM
training and prediction time complexity. We first investigate the case where
we only consider cross-sectional hierarchies. This case is indicated by ‘HL’
in Table 6. First, we note that adding our hierarchical loss objective adds a
component to the time complexity that scales with $n^{c,3}$, as we need to
compute (13). However, our sparse implementation of the hierarchical loss
reduces this component from $n^{c,3}$ to $n^{c,2}l^{c}$, effectively reducing
the scaling from cubic to quadratic in the number of bottom-level time series,
as $l^{c}$ is generally small. In the reconciliation scenarios, we always need
to compute a matrix inversion to solve Eq. (4) that scales cubically with the
number of cross-sectional aggregations $m^{c,*}$ or with the total number of
time series $m^{c}$. The first is not problematic as generally $m^{c,*}\ll
n^{c}$ in large-scale settings, but methods with this time complexity
consequently trade in performance, as we observed in Table 4. To empirically
verify the differences in asymptotic time complexity, we recorded the training
and prediction time for each scenario. We show timings for training and
prediction for a single store of the M5 dataset (4M training samples) and for
the entire M5 dataset (52M training samples), to provide an indication of
scaling when the problem size increases by an order of magnitude. First, we
note that using our sparse implementation of the HL reduces training time by a
factor of $3\times$ when training for all stores. Second, our sparse HL has a
prediction time similar to the baseline (SL). Next, we find that the training
time of our sparse hierarchical loss is two orders of magnitude faster than
reconciliation methods in the Separate aggregations-scenario. This is mainly
due to the many individual models that need to be trained in this scenario and
thus shows a clear benefit of having just a single model. We observe an order
of magnitude difference in prediction time when comparing the sparse
hierarchical loss to the Separate aggregations-scenario when predictinf all
stores. Again, this shows a clear benefit of having just a single model for
this forecasting task. For the Global-scenario, we see that reconciliation
methods require a smaller training time when training for all stores (about
twice less), however that scenario also did not give strong forecasting
performance as we established in Table 4. Also, the prediction time using our
sparse HL is an order of magnitude lower. As ML costs in production systems
mainly consists of prediction costs, having a lower prediction time is
beneficial.222For example, Google designed its first TPU for inference:
https://techcrunch.com/2017/05/17/google-announces-second-generation-of-
tensor-processing-unit-chips. Finally, we also show the time complexity of
using both cross-sectional and temporal hierarchies jointly, as indicated by
‘$HL+$’ in Table 6. Adding temporal hierarchies adds another matrix
multiplication that scales with the number of timesteps to the complexity. In
our experiments, we find that adding temporal hierarchies results in a twice
higher training time when training for all stores and a 50% higher prediction
time when predicting for all stores. We view it as potential future work to
investigate how to perform this end-to-end learning of both cross-sectional
and temporal hierarchies even more efficiently.

To conclude, we showed that our sparse HL incurs some additional training
overhead but no additional prediction overhead as compared to the base case
SL, whereas it does not require the additional reconciliation step that
reconciliation methods require.

Table 6: Computational time complexity and observed timings in seconds for all
scenarios. The complexity is indicated by how respectively the training time
and prediction time scales with respect to the default LightGBM
training/prediction time $T$, where $n^{t}$ denotes the number of timesteps
per time series, $n^{c}$ denotes the number of bottom-level time series in the
hierarchy, $n^{c}_{l}$ the number of time series in each level in the
hierarchy and $l^{c}$ the number of levels in the cross-sectional hierarchy,
$m^{c}$ ($m^{t}$) the total number of cross-sectional (temporal) aggregations,
and $m^{*}=m-n$.

|  |  | Complexity | Training time (s) | Prediction time (s)
---|---|---|---|---|---
Scenario/Objective | Metric | Reconciliation | Training | Prediction | 1 store | All stores | 1 store | All stores
Bottom-up |  |  |  |  |  |  |  |
SL | SL | None | $O(T(n^{t}n^{c}))$ | $O(T(n^{t}n^{c})+n^{t}n^{c,3})$ | 8 | 173 | 1.1 | 11
SL | SL | None | $O(T(n^{t}n^{c}))$ | $O(T(n^{t}n^{c})+n^{t}n^{c,2}l^{c})$ | 10 | 174 | 0.1 | 10
HL (dense) | HL (dense) | None | $O(T(n^{t}n^{c}+n^{t}n^{c,3}))$ | $O(T(n^{t}n^{c})+n^{t}n^{c,3})$ | 14 | 1,185 | 1.1 | 10
HL (sparse) | HL (sparse) | None | $O(T(n^{t}n^{c}+n^{t}n^{c,2}l^{c}))$ | $O(T(n^{t}n^{c})+n^{t}n^{c,2}l^{c})$ | 12 | 318 | 0.1 | 11
HL+ (sparse) | HL+ (sparse) | None | $O(T(n^{t}n^{c}+n^{t}n^{c,2}l^{c}n^{c}n^{t,2}l^{t}))$ | $O(T(n^{t}n^{c})+n^{t}n^{c,2}l^{c}n^{c}n^{t,2}l^{t}))$ |  | 723 |  | 15
Sep. agg. |  |  |  |  |  |  |  |
SL | SL | Base | $O(l\cdot T(n^{c}_{l}n^{t}_{l}))$ | $O(l^{c}\cdot T(n^{c}_{l}n^{t}_{l}))$ | 11 | 36,018 | 4.4 | 103
SL | SL | OLS | $O(l^{c}\cdot T(n^{t}_{l}n_{l}^{c})+m^{c*,3})$ | 4.5 | 149
SL | SL | WLS-struct | $O(l^{c}\cdot T(n^{t}_{l}n_{l}^{c})+m^{c*,3})$ | 4.5 | 151
SL | SL | WLS-var | $O(l^{c}\cdot T(n^{t}_{l}n_{l}^{c})+m^{c*,3})$ | 4.5 | 151
SL | SL | MinT-shrink | $O(l^{c}\cdot T(n^{t}_{l}n_{l}^{c})+m^{c,3})$ | 5.8 | 305
SL | SL | ERM | $O(l^{c}\cdot T(n^{t}_{l}n_{l}^{c})+m^{c,3})$ | 6.0 | 239
Global |  |  |  |  |  |  |  |
SL | SL | Base | $O(T(m^{t}m^{c}))$ | $O(T(m^{t}m^{c}))$ | 4 | 173 | 2.4 | 71
SL | SL | OLS | $O(T(m^{t}m^{c})+m^{c*,3})$ | 2.5 | 118
SL | SL | WLS-struct | $O(T(m^{t}m^{c})+m^{c*,3})$ | 2.5 | 120
SL | SL | WLS-var | $O(T(m^{t}m^{c})+m^{c*,3})$ | 2.5 | 120
SL | SL | MinT-shrink | $O(T(m^{t}m^{c})+m^{c,3})$ | 4.2 | 274
SL | SL | ERM | $O(T(m^{t}m^{c})+m^{c,3})$ | 4.0 | 207

###  5.2 Proprietary datasets

At our e-commerce partner bol, a LightGBM-based forecasting model is used as
the primary product forecasting model. The model is used to forecast weekly
product demand for 12 weeks. Every day, 12 separate models are trained, each
tasked to forecast demand for a single week for every product. The model is
used to forecast the majority of the products on sale at any moment, which are
approximately 5 million unique items. We investigate the use of our sparse
hierarchical loss function as a drop-in replacement for the existing Tweedie
loss that is used within the company. The offline dataset consists of 36M
training samples from the period January 2017 to the end of June 2021. We test
on 55M samples from the period July 2021 to January 2022. The baseline model
for every weekly forecast model is a LightGBM model with a Tweedie loss (TL).
The baseline model uses 19 proprietary features.

#### Experimental setup

We replace the TL with our HL and investigate forecasting performance on the
test set. For the HL, we use the proprietary aggregations product_group and
seasonality_group, each containing respectively $\pm$70 and $\pm$6,000 unique
values. We have $n^{c}=\pm 5$M bottom-level time series and $m^{c,*}=\pm
6,070$ aggregated time series across $l^{c}=4$ levels: product (bottom-level),
product_group, seasonality_group and total.

#### Results

On average, we find our sparse HL outperforms the existing TL model by about
1–2% on RMSE and $\pm$10% on MAE. We further investigate the performance by
investigating how the RMSE and MAE vary across the 12 forecasting horizons and
weekly demand buckets, and present the results in Figure 1. We find that our
sparse HL performs best on both RMSE and MAE on the lower weekly demand
buckets (up to 100 products sold per week), where it outperforms the TL
averaged over all the forecasting horizons. The TL is clearly better for
higher weekly demand buckets, commonly outperforming the HL and SL by up to
5%. Next, we investigate forecasting performance across the cross-sectional
hierarchies that we defined. We show the results in Figure 2. We find that for
most forecasting horizons, the HL and SL outperform the TL, with an average
outperformance of the HL over the TL of $\pm$10% at the product level, $\pm$5%
at product group level and $\pm$4–7% at seasonality group level. Hence, we are
able to confirm some of the results we found in the M5 experiment, although
the baseline SL performed quite strong in this experiment as well. We believe
this is due to the M5 experiment having much more hierarchical levels (12 as
compared to the 4 we used for our proprietary dataset experiment), since the
HL is equal to the SL in the case of no hierarchies, and with fewer
hierarchies the HL thus becomes closer to the SL. Hence, we believe our HL is
most useful in settings with both many timeseries as well as many hierarchies.
To conclude, we find that this experiment demonstrates the usefulness of our
sparse HL applied to a production forecasting system at a major e-commerce
platform.

(a) RMSE (b) MAE

Figure 1: Forecasting results for the primary product forecasting model at our
e-commerce partner bol. We show RMSE (a, left) and MAE (b, right) by weekly
demand bucket relative to the Tweedie loss baseline for each forecasting
horizon (week). The Hierarchical loss outperforms the Tweedie loss on RMSE and
MAE on smaller weekly demand buckets. Figure 2: Forecasting results for the
primary product forecasting model at our e-commerce partner bol. We show RMSE
(left column of figures) and MAE (right column of figures) by aggregation
level relative to the Tweedie loss baseline for each forecasting horizon
(week). The Hierarchical loss commonly outperforms the Tweedie loss on every
aggregation level.

##  6 Conclusion

We introduced a sparse Hierarchical loss function to perform hierarchical
forecasting in large-scale settings. We demonstrated that we are able to
outperform existing hierarchical forecasting methods both in terms of
performance as measured by RMSE and MAE by up to 10% as well as in terms of
computational time required to perform the end-to-end hierarchical forecasting
in large-scale settings, reducing prediction time as compared to the best
hierarchical forecasting reconciliation method by an order of magnitude. We
empirically verified our sparse hierarchical loss in an offline test for bol,
where we confirmed the results from our offline test on the public M5 dataset.

In addition to our main contributions, one of our main learnings has been that
we could not find a benefit of having multiple models for separate
aggregations in the hierarchy, as the bottom-up scenario we employed
consistently outperformed other scenarios. Secondly, we did not find a benefit
of training a model whilst adhering to both cross-sectional and temporal
hierarchies jointly.

Limitations of our work are that we did not consider the probabilistic
forecasting setting, where reconciled forecasts are required across an entire
forecast distribution.

For future work, we aim to extend our work to the setting of probabilistic
forecasting by combining our sparse hierarchical loss with existing
probabilistic forecasting frameworks from, e.g., [33, 34, 35]. In addition, we
seek to further investigate solutions for efficiently combining cross-
sectional and temporal hierarchies.

## Acknowledgments

This research was (partially) funded by the Hybrid Intelligence Center, a
10-year program funded by the Dutch Ministry of Education, Culture and Science
through the Netherlands Organisation for Scientific
Research333https://www.hybrid-intelligence-centre.nl/, and Ahold
Delhaize.444https://www.aholddelhaize.com/

All content represents the opinion of the authors, which is not necessarily
shared or endorsed by their respective employers and/or sponsors.

## References

  * Hyndman et al. [2011] R. J. Hyndman, R. A. Ahmed, G. Athanasopoulos, H. L. Shang,  Optimal combination forecasts for hierarchical time series,  Computational Statistics & Data Analysis 55 (2011) 2579–2589. doi:10.1016/j.csda.2011.03.006.
  * Athanasopoulos et al. [2017] G. Athanasopoulos, R. J. Hyndman, N. Kourentzes, F. Petropoulos,  Forecasting with temporal hierarchies,  European Journal of Operational Research 262 (2017) 60–74. doi:10.1016/j.ejor.2017.02.046.
  * Rangapuram et al. [2023] S. S. Rangapuram, S. Kapoor, R. S. Nirwan, P. Mercado, T. Januschowski, Y. Wang, M. Bohlke-Schneider,  Coherent Probabilistic Forecasting of Temporal Hierarchies,  in: Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, PMLR, 2023, pp. 9362–9376.
  * Theodosiou and Kourentzes [2021] F. Theodosiou, N. Kourentzes, Forecasting with Deep Temporal Hierarchies, 2021. doi:10.2139/ssrn.3918315.
  * Wickramasuriya et al. [2019] S. L. Wickramasuriya, G. Athanasopoulos, R. J. Hyndman,  Optimal Forecast Reconciliation for Hierarchical and Grouped Time Series Through Trace Minimization,  Journal of the American Statistical Association 114 (2019) 804–819. doi:10.1080/01621459.2018.1448825.
  * Rangapuram et al. [2021] S. S. Rangapuram, L. D. Werner, K. Benidis, P. Mercado, J. Gasthaus, T. Januschowski,  End-to-End Learning of Coherent Probabilistic Forecasts for Hierarchical Time Series,  in: Proceedings of the 38th International Conference on Machine Learning, PMLR, 2021, pp. 8832–8843.
  * Makridakis et al. [2022] S. Makridakis, E. Spiliotis, V. Assimakopoulos,  M5 accuracy competition: Results, findings, and conclusions,  International Journal of Forecasting 38 (2022) 1346–1364. doi:10.1016/j.ijforecast.2021.11.013.
  * Böse et al. [2017] J.-H. Böse, V. Flunkert, J. Gasthaus, T. Januschowski, D. Lange, D. Salinas, S. Schelter, M. Seeger, Y. Wang,  Probabilistic demand forecasting at scale,  Proceedings of the VLDB Endowment 10 (2017) 1694–1705. doi:10.14778/3137765.3137775.
  * Benidis et al. [2023] K. Benidis, S. S. Rangapuram, V. Flunkert, Y. Wang, D. Maddix, C. Turkmen, J. Gasthaus, M. Bohlke-Schneider, D. Salinas, L. Stella, F.-X. Aubet, L. Callot, T. Januschowski,  Deep Learning for Time Series Forecasting: Tutorial and Literature Survey,  ACM Computing Surveys 55 (2023) 1–36. doi:10.1145/3533382. arXiv:2004.10240.
  * Januschowski et al. [2022] T. Januschowski, Y. Wang, K. Torkkola, T. Erkkilä, H. Hasson, J. Gasthaus,  Forecasting with trees,  International Journal of Forecasting 38 (2022) 1473–1481. doi:10.1016/j.ijforecast.2021.10.004.
  * Ke et al. [2017] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, T.-Y. Liu,  LightGBM: A Highly Efficient Gradient Boosting Decision Tree,  in: Advances in Neural Information Processing Systems 30, Curran Associates, Inc., 2017, pp. 3146–3154.
  * Chen and Guestrin [2016] T. Chen, C. Guestrin,  XGBoost: A Scalable Tree Boosting System,  in: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, ACM, San Francisco California USA, 2016, pp. 785–794. doi:10.1145/2939672.2939785.
  * Hyndman et al. [2016] R. J. Hyndman, A. J. Lee, E. Wang,  Fast computation of reconciled forecasts for hierarchical and grouped time series,  Computational Statistics & Data Analysis 97 (2016) 16–32. doi:10.1016/j.csda.2015.11.007.
  * Taieb et al. [2017] S. B. Taieb, J. W. Taylor, R. J. Hyndman,  Coherent Probabilistic Forecasts for Hierarchical Time Series,  in: International Conference on Machine Learning, 2017, pp. 3348–3357.
  * Ben Taieb and Koo [2019] S. Ben Taieb, B. Koo,  Regularized Regression for Hierarchical Forecasting Without Unbiasedness Conditions,  in: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, ACM, Anchorage AK USA, 2019, pp. 1337–1347. doi:10.1145/3292500.3330976.
  * Taieb [2017] S. B. Taieb,  Sparse and Smooth Adjustments for Coherent Forecasts in Temporal Aggregation of Time Series,  in: Proceedings of the Time Series Workshop at NIPS 2016, PMLR, 2017, pp. 16–26.
  * Panagiotelis et al. [2021] A. Panagiotelis, G. Athanasopoulos, P. Gamakumara, R. J. Hyndman,  Forecast reconciliation: A geometric view with new insights on bias correction,  International Journal of Forecasting 37 (2021) 343–359. doi:10.1016/j.ijforecast.2020.06.004.
  * Girolimetto and Di Fonzo [2023] D. Girolimetto, T. Di Fonzo, Point and probabilistic forecast reconciliation for general linearly constrained multiple time series, 2023. doi:10.48550/arXiv.2305.05330. arXiv:2305.05330.
  * Han et al. [2021] X. Han, S. Dasgupta, J. Ghosh,  Simultaneously Reconciled Quantile Forecasting of Hierarchically Related Time Series,  in: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, PMLR, 2021, pp. 190–198.
  * Hyndman and Athanasopoulos [2021] R. J. Hyndman, G. Athanasopoulos, Forecasting: Principles and Practice (3rd Ed), OTexts: Melbourne, Australia, 2021.
  * Schäfer and Strimmer [2005] J. Schäfer, K. Strimmer,  A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics,  Statistical Applications in Genetics and Molecular Biology 4 (2005). doi:10.2202/1544-6115.1175.
  * Touloumis [2015] A. Touloumis,  Nonparametric Stein-type Shrinkage Covariance Matrix Estimators in High-Dimensional Settings,  Computational Statistics & Data Analysis 83 (2015) 251–261. doi:10.1016/j.csda.2014.10.018. arXiv:1410.4726.
  * Makridakis et al. [2021] S. Makridakis, E. Spiliotis, V. Assimakopoulos,  The M5 competition: Background, organization, and implementation,  International Journal of Forecasting (2021). doi:10.1016/j.ijforecast.2021.07.007.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. Van Der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. Van Mulbregt, SciPy 1.0 Contributors, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. De Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko, Y. Vázquez-Baeza,  SciPy 1.0: Fundamental algorithms for scientific computing in Python,  Nature Methods 17 (2020) 261–272. doi:10.1038/s41592-019-0686-2.
  * Box and Pierce [1970] G. E. P. Box, D. A. Pierce,  Distribution of Residual Autocorrelations in Autoregressive-Integrated Moving Average Time Series Models,  Journal of the American Statistical Association 65 (1970) 1509–1526. doi:10.2307/2284333. arXiv:2284333.
  * Hyndman et al. [2008] R. Hyndman, A. B. Koehler, J. K. Ord, R. D. Snyder, Forecasting with Exponential Smoothing: The State Space Approach, Springer Science & Business Media, 2008.
  * Assimakopoulos and Nikolopoulos [2000] V. Assimakopoulos, K. Nikolopoulos,  The theta model: A decomposition approach to forecasting,  International Journal of Forecasting 16 (2000) 521–530. doi:10.1016/S0169-2070(00)00066-2.
  * Croston [1972] J. D. Croston,  Forecasting and Stock Control for Intermittent Demands,  Operational Research Quarterly (1970-1977) 23 (1972) 289–303. doi:10.2307/3007885. arXiv:3007885.
  * Kunz et al. [2023] M. Kunz, S. Birr, M. Raslan, L. Ma, Z. Li, A. Gouttes, M. Koren, T. Naghibi, J. Stephan, M. Bulycheva, M. Grzeschik, A. Kekić, M. Narodovitch, K. Rasul, J. Sieber, T. Januschowski, Deep Learning based Forecasting: A case study from the online fashion industry, 2023. doi:10.48550/arXiv.2305.14406. arXiv:2305.14406.
  * Lim et al. [2021] B. Lim, S. Ö. Arık, N. Loeff, T. Pfister,  Temporal Fusion Transformers for interpretable multi-horizon time series forecasting,  International Journal of Forecasting 37 (2021) 1748–1764. doi:10.1016/j.ijforecast.2021.03.012.
  * Li et al. [2019] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, X. Yan,  Enhancing the Locality and Breaking the Memory Bottleneck of Transformer on Time Series Forecasting,  in: Advances in Neural Information Processing Systems 32, Curran Associates, Inc., 2019, pp. 5244–5254.
  * Akiba et al. [2019] T. Akiba, S. Sano, T. Yanase, T. Ohta, M. Koyama,  Optuna: A Next-generation Hyperparameter Optimization Framework,  in: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD ’19, Association for Computing Machinery, New York, NY, USA, 2019, pp. 2623–2631. doi:10.1145/3292500.3330701.
  * Sprangers et al. [2021] O. Sprangers, S. Schelter, M. de Rijke,  Probabilistic Gradient Boosting Machines for Large-Scale Probabilistic Regression,  in: Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, KDD ’21, Association for Computing Machinery, New York, NY, USA, 2021, pp. 1510–1520. doi:10.1145/3447548.3467278.
  * Hasson et al. [2021] H. Hasson, B. Wang, T. Januschowski, J. Gasthaus,  Probabilistic Forecasting: A Level-Set Approach,  in: Advances in Neural Information Processing Systems, volume 34, Curran Associates, Inc., 2021, pp. 6404–6416.
  * Stankeviciute et al. [2021] K. Stankeviciute, A. M. Alaa, M. van der Schaar,  Conformal Time-series Forecasting,  in: Advances in Neural Information Processing Systems, volume 34, Curran Associates, Inc., 2021, pp. 6216–6228.

##  Appendix A Derivation of gradient

We have the following hierarchical forecasting problem:

$\displaystyle\tilde{\textbf{Y}}=\underbrace{\begin{bmatrix}1&1\\\ 1&0\\\
0&1\end{bmatrix}}_{S^{c}}\underbrace{\begin{bmatrix}\hat{\textbf{y}}^{0}_{0}&\hat{\textbf{y}}^{0}_{1}\\\
\hat{\textbf{y}}^{1}_{0}&\hat{\textbf{y}}^{1}_{1}\\\
\end{bmatrix}}_{\hat{\textbf{Y}}}\underbrace{\begin{bmatrix}1&1&0\\\ 1&0&1\\\
\end{bmatrix}}_{(S^{t})^{\intercal}}\;,$

which results in the following gradient and hessian:

$\displaystyle\begin{bmatrix}\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{0}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{0}_{1}}\\\ \frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}&\frac{\partial
L}{\partial\hat{\textbf{y}}^{1}_{1}}\end{bmatrix}$
$\displaystyle=(S^{c})^{\intercal}\left[\frac{S^{c}\left(\hat{\textbf{Y}}-\bar{\textbf{Y}}\right)(S^{t})^{\intercal}}{d^{c}\otimes(d^{t})^{\intercal}}\right]S^{t}$
$\displaystyle=(S^{c})^{\intercal}\left[\vbox{\hbox{$\displaystyle\frac{S^{c}\left(\hat{\textbf{Y}}-\bar{\textbf{Y}}\right)(S^{t})^{\intercal}}{\underbrace{\left[l^{c}\sum_{j=1}^{n^{c}}S^{c}_{ij}\right]}_{d^{c}}\otimes\underbrace{\left[l^{t}\sum_{j=1}^{n^{t}}S^{t}_{ij}\right]^{\intercal}}_{(d^{t})^{\intercal}}}$}}\right]S^{t}$
$\displaystyle=\begin{bmatrix}1&1&0\\\ 1&0&1\\\
\end{bmatrix}\left[\frac{\begin{bmatrix}1&1\\\ 1&0\\\
0&1\end{bmatrix}\left(\begin{bmatrix}\hat{\textbf{y}}^{0}_{0}&\hat{\textbf{y}}^{0}_{1}\\\
\hat{\textbf{y}}^{1}_{0}&\hat{\textbf{y}}^{1}_{1}\end{bmatrix}-\begin{bmatrix}\bar{\textbf{y}}^{0}_{0}&\bar{\textbf{y}}^{0}_{1}\\\
\bar{\textbf{y}}^{1}_{0}&\bar{\textbf{y}}^{1}_{1}\end{bmatrix}\right)\begin{bmatrix}1&1&0\\\
1&0&1\\\
\end{bmatrix}}{\left(2\begin{bmatrix}2&1&1\end{bmatrix}^{\intercal}\right)\otimes\left(2\begin{bmatrix}2&1&1\end{bmatrix}\right)}\right]\begin{bmatrix}1&1\\\
1&0\\\ 0&1\end{bmatrix}\;,$ $\displaystyle=\begin{bmatrix}1&1&0\\\ 1&0&1\\\
\end{bmatrix}\left[\frac{\begin{bmatrix}1&1\\\ 1&0\\\
0&1\end{bmatrix}\left(\begin{bmatrix}\overbrace{\hat{\textbf{y}}^{0}_{0}-\bar{\textbf{y}}^{0}_{0}}^{e^{0}_{0}}&\overbrace{\hat{\textbf{y}}^{0}_{1}-\bar{\textbf{y}}^{0}_{1}}^{e^{0}_{1}}\\\
\underbrace{\hat{\textbf{y}}^{1}_{0}-\bar{\textbf{y}}^{1}_{0}}_{e^{1}_{0}}&\underbrace{\hat{\textbf{y}}^{1}_{1}-\bar{\textbf{y}}^{1}_{1}}_{e^{1}_{1}}\end{bmatrix}\right)\begin{bmatrix}1&1&0\\\
1&0&1\\\ \end{bmatrix}}{\begin{bmatrix}16&8&8\\\ 8&4&4\\\ 8&4&4\\\
\end{bmatrix}}\right]\begin{bmatrix}1&1\\\ 1&0\\\ 0&1\end{bmatrix}\;,$
$\displaystyle=\begin{bmatrix}1&1&0\\\ 1&0&1\\\
\end{bmatrix}\begin{bmatrix}\frac{1}{16}(e^{0}_{0}+e^{1}_{0}+e^{0}_{1}+e^{1}_{1})&\frac{1}{8}(e^{0}_{0}+e^{1}_{0})&\frac{1}{8}(e^{0}_{1}+e^{1}_{1})\\\
\frac{1}{8}(e^{0}_{0}+e^{0}_{1})&\frac{1}{4}e^{0}_{0}&\frac{1}{4}e^{0}_{1}\\\
\frac{1}{8}(e^{1}_{0}+e^{1}_{1})&\frac{1}{4}e^{1}_{0}&\frac{1}{4}e^{1}_{1}\\\
\end{bmatrix}\begin{bmatrix}1&1\\\ 1&0\\\ 0&1\end{bmatrix}\;,$
$\displaystyle=\begin{bmatrix}\frac{1}{16}(e^{0}_{0}+e^{1}_{0}+e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{0}_{0}+e^{1}_{0})+\frac{1}{8}(e^{0}_{0}+e^{0}_{1})+\frac{1}{4}e^{0}_{0}&\frac{1}{16}(e^{0}_{0}+e^{1}_{0}+e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{0}_{0}+e^{0}_{1})+\frac{1}{4}e^{0}_{1}\\\
\frac{1}{16}(e^{0}_{0}+e^{1}_{0}+e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{0}_{0}+e^{1}_{0})+\frac{1}{8}(e^{1}_{0}+e^{1}_{1})+\frac{1}{4}e^{1}_{0}&\frac{1}{16}(e^{0}_{0}+e^{1}_{0}+e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{0}_{1}+e^{1}_{1})+\frac{1}{8}(e^{1}_{0}+e^{1}_{1})+\frac{1}{4}e^{1}_{1}\\\
\end{bmatrix}\;,$
$\displaystyle=\begin{bmatrix}\frac{9}{16}e_{0}^{0}+\frac{3}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{1}+\frac{1}{16}e_{1}^{1}&\frac{9}{16}e_{0}^{1}+\frac{3}{16}e_{0}^{1}+\frac{3}{16}e_{1}^{1}+\frac{1}{16}e_{1}^{0}\\\
\frac{9}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{0}+\frac{3}{16}e_{1}^{1}+\frac{1}{16}e_{0}^{1}&\frac{9}{16}e_{1}^{1}+\frac{3}{16}e_{1}^{0}+\frac{3}{16}e_{0}^{1}+\frac{1}{16}e_{0}^{0}\\\
\end{bmatrix}\;.$

##  Appendix B M5 Dataset

For each of the scenarios of our experiments in Section 5, we construct a set
of features for the LightGBM model as given in Table 7. To facilitate the most
‘fair’ comparison across methods, each model has the same features, and for
the time series aggregations in the hierarchy we construct the features taken
over the aggregation.

Table 7: Features used for the M5 dataset in our experiments. Feature | Description
---|---
Aggregation | Aggregation level in the hierarchy
Value | Identifier of time series of this aggregation
sales_lag1-7 | Lagged sales (target) (7 features)
sales_lag28 | Sales 28 days ago
sales_lag56 | Sales 56 days ago
sales_lag364 | Sales last year
sales_lag1_mavg7 | Moving average of sales last 7 days
sales_lag1_mavg28 | Moving average of sales last 28 days
sales_lag1_mavg56 | Moving average of sales last 56 days
dayofweek | Day of the week
dayofmonth | Day of the month
weekofyear | Week of year
monthofyear | Month of year
sell_price_avg | Sell price (average if aggregation)
sell_price_change | Day-to-day change in sell price
weeks_on_sale_avg | Weeks on sale
snap_CA | State indicator for California
snap_TX | State indicator for Texas
snap_WI | State indicator for Wyoming
event_type_1_enc | Encoded events
event_type_2_enc | Encoded events

##  Appendix C M5 model training & optimization

We optimize the hyperparameters of our LightGBM models using Optuna [32],
using the settings found in Table 8. The validation is performed on a rolling-
forward basis for 3 validation sets, where we use three years of data to
predict the next 28 days ahead. After the hyperparameter optimization
procedure, we use the average number of iterations at which the lowest
validation loss was achieved across the 3 validation sets as the number of
estimators to use in our final model. The final model uses the last three
years of data preceding the first day in the test set.

Table 8: Key hyperparameters used in our experiments. The parameters with a search range included are optimized in a hyperparameter search. Parameter | Description | Default value | Search range
---|---|---|---
n_estimators | Number of trees in each model | 2000 | Lowest validation loss
n_trials | Number of optimization trials to run | 100 |
learning_rate | Learning rate | 0.05 |
n_validation_sets | Number of validation sets | 3 |
n_days_test | Number of days in validation and test sets | 28 |
max_levels_random | Max. number of levels when using a random hierarchy | 2 |
max_categories_per _random_level | Max. categories per level in the random hierarchy | 1000 |
hier_freq | Frequency of performing the randomized hierarchical aggregation | 1 | uniform($1$, $10$)
lambda_l1 | L1-regularization | 0 | log_uniform($10^{-8}$, $10^{1}$)
lambda_l2 | L2-regularization | 0 | log_uniform($10^{-8}$, $10^{1}$)
num_leaves | Max. number of leaves per tree | 31 | uniform($2^{3}$, $2^{10}$)
feature_fraction | Fraction of features to use to build a tree | 1.0 | uniform($0.4$, $1.0$)
bagging_fraction | Fraction of training samples to use to build a tree | 1.0 | uniform($0.4$, $1.0$)
bagging_freq | Frequency at which to create a new bagging batch | 1.0 | uniform($1$, $7$)
min_child_samples | Minimum number of samples per leaf | 20 | log_uniform($5$, $5000$)

##  Appendix D Experiments

Table 9: Forecasting results for all stores on the M5 dataset. We report mean
RMSE scores across 10 random seeds with standard deviation in brackets. Lower
is better.

|  |  |  |  |  | Store | Product | State |  |
---|---|---|---|---|---|---|---|---|---|---
Scen./Obj. | Metric | Reconciliation | Product | Department | Category | Department | Category | Total | Store | State | Department | Category | Total | Total | All series
Bottom-up |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | None | 2.11 (0.00) | 580 (3.95) | 1096 (13.71) | 110 (0.60) | 207 (1.21) | 416 (2.07) | 8.8 (0.03) | 4.26 (0.01) | 263 (1.30) | 510 (3.81) | 997 (8.44) | 2086 (38.53) | 22.4 (0.16)
SL | HL | None | 2.11 (0.00) | 564 (5.58) | 1039 (14.64) | 107 (0.74) | 200 (1.69) | 408 (2.37) | 8.76 (0.06) | 4.25 (0.01) | 256 (1.30) | 490 (3.93) | 983 (5.93) | 2079 (33.46) | 21.8 (0.15)
HL | HL | None | 2.11 (0.00) | 511 (8.93) | 882 (20.54) | 102 (0.90) | 185 (1.48) | 388 (4.21) | 8.5 (0.02) | 4.21 (0.01) | 233 (2.42) | 427 (4.78) | 878 (17.17) | 1819 (78.42) | 19.5 (0.38)
HL | SL | None | 2.11 (0.00) | 511 (9.47) | 886 (24.67) | 103 (1.08) | 187 (2.12) | 390 (3.80) | 8.46 (0.04) | 4.19 (0.01) | 234 (3.15) | 430 (7.18) | 881 (15.27) | 1813 (64.6) | 19.6 (0.36)
TL | HL | None | 2.11 (0.00) | 549 (9.21) | 1055 (22.68) | 109 (0.58) | 207 (1.27) | 415 (2.58) | 8.73 (0.03) | 4.25 (0.01) | 256 (2.09) | 502 (5.04) | 978 (13.8) | 2007 (63.57) | 21.8 (0.34)
TL | SL | None | 2.11 (0.00) | 544 (7.66) | 1017 (25.23) | 110 (0.40) | 208 (1.20) | 416 (2.41) | 8.68 (0.04) | 4.25 (0.01) | 256 (1.55) | 497 (5.23) | 972 (7.67) | 1943 (41.21) | 21.5 (0.25)
TL | TL | None | 2.46 (0.00) | 1579 (9.71) | 3075 (20.5) | 193 (0.87) | 379 (1.84) | 717 (5.17) | 13.41 (0.05) | 5.68 (0.02) | 566 (3.11) | 1111 (6.42) | 2074 (18.68) | 5646 (60.40) | 53.3 (0.39)
Sep. agg. |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 2.11 (0.00) | 835 (23.47) | 1417 (69.78) | 130 (1.86) | 236 (5.25) | 474 (8.71) | 8.88 (0.02) | 4.23 (0.01) | 322 (11.48) | 682 (19.02) | 1269 (47.11) | 3339 (86.79) | 30.1 (0.21)
SL | SL | OLS | 2.11 (0.00) | 804 (14.22) | 1548 (22.24) | 120 (1.49) | 220 (1.95) | 443 (2.37) | 8.77 (0.02) | 4.25 (0.01) | 315 (4.15) | 605 (7.98) | 1222 (9.69) | 3125 (24.25) | 29.1 (0.22)
SL | SL | WLS-struct | 2.10 (0.00) | 731 (8.55) | 1505 (16.09) | 113 (0.58) | 217 (1.13) | 426 (2.17) | 8.71 (0.02) | 4.23 (0.01) | 291 (2.38) | 591 (4.82) | 1153 (7.53) | 2909 (23.52) | 27.6 (0.22)
SL | SL | WLS-var | 2.11 (0.00) | 647 (7.23) | 1345 (18.95) | 109 (0.51) | 211 (1.43) | 411 (2.06) | 8.72 (0.03) | 4.24 (0.01) | 272 (2.06) | 555 (5.58) | 1068 (8.77) | 2541 (31.06) | 25.1 (0.25)
SL | SL | MinT-shrink | 2.11 (0.00) | 667 (18.15) | 1397 (44.32) | 106 (1.15) | 206 (3.14) | 405 (5.45) | 8.76 (0.02) | 4.24 (0.01) | 272 (4.92) | 558 (12.59) | 1091 (22.52) | 2706 (80.28) | 25.8 (0.63)
SL | SL | ERM | 2.58 (0.01) | 728 (48.01) | 1411 (118.36) | 117 (3.66) | 213 (8.73) | 443 (20.5) | 10.30 (0.10) | 5.21 (0.04) | 307 (14.61) | 583 (33.15) | 1215 (75.3) | 3101 (258.29) | 28.1 (1.83)
Global |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 2.16 (0.00) | 771 (54.94) | 1587 (256.06) | 120 (2.11) | 227 (8.81) | 458 (18.45) | 9.09 (0.05) | 4.4 (0.01) | 329 (18.51) | 650 (31.09) | 1807 (1203.62) | 3273 (338.71) | 32.8 (6.60)
SL | SL | OLS | 2.13 (0.01) | 764 (58.91) | 1524 (181.38) | 117 (10.32) | 227 (27.24) | 484 (109.14) | 8.99 (0.05) | 4.34 (0.03) | 316 (42.99) | 637 (107.63) | 1371 (417.84) | 3114 (417.38) | 29.9 (2.52)
SL | SL | WLS-struct | 2.13 (0.00) | 800 (45.37) | 1686 (103.86) | 118 (2.52) | 235 (5.96) | 463 (11.90) | 9.09 (0.03) | 4.35 (0.01) | 314 (10.72) | 651 (24.97) | 1264 (55.94) | 3231 (279.96) | 30.4 (1.68)
SL | SL | WLS-var | 2.14 (0.01) | 876 (43.62) | 1867 (106.77) | 130 (2.31) | 263 (7.10) | 508 (18.51) | 9.10 (0.05) | 4.35 (0.01) | 345 (11.71) | 727 (30.30) | 1374 (80.66) | 3463 (283.05) | 33.2 (1.84)
SL | SL | MinT-shrink | 2.18 (0.02) | 730 (50.73) | 1545 (155.65) | 115 (2.42) | 228 (5.97) | 478 (19.61) | 9.29 (0.16) | 4.46 (0.06) | 293 (12.06) | 599 (33.27) | 1235 (89.41) | 3211 (432.02) | 29.1 (2.33)
SL | SL | ERM | 2.56 (0.14) | 923 (113.39) | 1857 (271.23) | 138 (16.52) | 265 (34.54) | 558 (105.87) | 10.52 (0.62) | 5.23 (0.34) | 382 (56.45) | 762 (118.54) | 1607 (371.34) | 3752 (669.36) | 35.7 (3.74)

Table 10: Forecasting results for all stores on the M5 dataset. We report mean
MAE scores across 10 random seeds with standard deviation in brackets. Lower
is better.

|  |  |  |  |  | Store | Product | State |  |
---|---|---|---|---|---|---|---|---|---|---
Scen./Obj. | Metric | Reconciliation | Product | Department | Category | Department | Category | Total | Store | State | Department | Category | Total | Total | All series
Bottom-up |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | None | 1.07 (0.00) | 420 (5.51) | 772 (12.94) | 70 (0.43) | 131 (0.84) | 312 (2.27) | 4.41 (0.01) | 2.2 (0.00) | 181 (1.45) | 343 (2.57) | 812 (7.46) | 1628 (24.81) | 2.2 (0.01)
SL | HL | None | 1.07 (0.00) | 410 (5.82) | 743 (12.96) | 68 (0.37) | 128 (0.75) | 309 (1.79) | 4.39 (0.01) | 2.19 (0.00) | 177 (1.29) | 332 (2.25) | 808 (5.34) | 1653 (40.22) | 2.2 (0.01)
HL | HL | None | 1.06 (0.00) | 342 (7.54) | 598 (18.24) | 63 (0.58) | 117 (1.36) | 293 (4.26) | 4.33 (0.01) | 2.18 (0.00) | 154 (2.01) | 283 (4.61) | 724 (14.31) | 1427 (81.18) | 2.1 (0.01)
HL | SL | None | 1.06 (0.00) | 340 (6.2) | 597 (18.44) | 63 (0.56) | 117 (1.19) | 294 (3.06) | 4.33 (0.01) | 2.18 (0.00) | 154 (1.91) | 284 (4.49) | 725 (13.35) | 1417 (57.53) | 2.1 (0.01)
TL | HL | None | 1.05 (0.00) | 350 (4.5) | 634 (10.36) | 65 (0.17) | 122 (0.49) | 299 (1.69) | 4.27 (0.01) | 2.15 (0.00) | 160 (0.69) | 305 (1.95) | 747 (9.07) | 1425 (50.6) | 2.1 (0.00)
TL | SL | None | 1.06 (0.00) | 359 (4.18) | 647 (10.08) | 66 (0.25) | 124 (0.52) | 305 (1.54) | 4.3 (0.01) | 2.17 (0.00) | 164 (1.27) | 311 (2.84) | 778 (4.19) | 1457 (28.24) | 2.1 (0.00)
TL | TL | None | 1.09 (0.00) | 866 (4.9) | 1843 (11.88) | 102 (0.44) | 210 (1.1) | 515 (4.95) | 5.03 (0.01) | 2.35 (0.00) | 306 (1.48) | 645 (3.75) | 1601 (17.55) | 4647 (59.66) | 2.8 (0.01)
Sep. agg. |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.07 (0.00) | 467 (13.86) | 819 (36.81) | 71 (0.62) | 144 (2.96) | 336 (9.36) | 4.23 (0.01) | 2.13 (0.00) | 187 (4.24) | 393 (11.11) | 971 (50.11) | 2523 (105.3) | 2.2 (0.00)
SL | SL | OLS | 1.03 (0.00) | 455 (8.08) | 895 (10.88) | 70 (0.73) | 129 (0.87) | 311 (1.76) | 4.24 (0.01) | 2.13 (0.00) | 181 (1.81) | 345 (2.84) | 910 (9.07) | 2430 (45.28) | 2.2 (0.01)
SL | SL | WLS-struct | 1.03 (0.00) | 417 (2.56) | 876 (8.99) | 64 (0.28) | 124 (0.66) | 298 (1.64) | 4.22 (0.00) | 2.12 (0.00) | 168 (1) | 334 (2.27) | 842 (6.11) | 2315 (22.33) | 2.1 (0.00)
SL | SL | WLS-var | 1.04 (0.00) | 383 (3.33) | 771 (10.62) | 63 (0.25) | 120 (0.66) | 286 (1.48) | 4.24 (0.00) | 2.13 (0.00) | 162 (1) | 318 (2.75) | 770 (6.87) | 1887 (32.1) | 2.1 (0.00)
SL | SL | MinT-shrink | 1.04 (0.00) | 393 (10.24) | 811 (27.92) | 62 (0.56) | 118 (1.64) | 283 (4.33) | 4.23 (0.01) | 2.13 (0.00) | 162 (2.44) | 316 (7.33) | 796 (18.35) | 2116 (81.24) | 2.1 (0.01)
SL | SL | ERM | 1.26 (0.00) | 490 (28.29) | 935 (91.26) | 76 (1.99) | 139 (5.8) | 335 (17.11) | 5.26 (0.03) | 2.67 (0.01) | 201 (7.9) | 383 (22.22) | 966 (63.05) | 2559 (296.35) | 2.6 (0.04)
Global |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
SL | SL | Base | 1.12 (0.01) | 436 (22.31) | 900 (156.94) | 70 (1.01) | 133 (3.84) | 326 (15.88) | 4.37 (0.01) | 2.18 (0.01) | 188 (8.04) | 373 (13.48) | 1314 (611.36) | 2569 (320.87) | 2.3 (0.04)
SL | SL | OLS | 1.05 (0.01) | 456 (32.4) | 906 (97.58) | 73 (11.48) | 139 (28.62) | 352 (93.72) | 4.38 (0.03) | 2.19 (0.04) | 197 (41.58) | 388 (103.17) | 1033 (336.46) | 2441 (389.89) | 2.3 (0.13)
SL | SL | WLS-struct | 1.05 (0.00) | 456 (17.42) | 970 (53.68) | 68 (0.83) | 133 (2.1) | 323 (7.31) | 4.4 (0.01) | 2.18 (0.01) | 181 (3.3) | 368 (10.47) | 937 (40.46) | 2552 (258.94) | 2.2 (0.02)
SL | SL | WLS-var | 1.06 (0.00) | 522 (18.96) | 1082 (59.2) | 76 (0.8) | 148 (3.13) | 353 (13.72) | 4.43 (0.02) | 2.19 (0.01) | 206 (4.75) | 415 (14.39) | 1017 (62.93) | 2744 (260.65) | 2.3 (0.02)
SL | SL | MinT-shrink | 1.06 (0.01) | 474 (24.59) | 960 (92.69) | 74 (1.56) | 141 (3.13) | 358 (11.28) | 4.49 (0.06) | 2.22 (0.02) | 190 (7) | 367 (22.7) | 935 (87.38) | 2605 (424.61) | 2.3 (0.03)
SL | SL | ERM | 1.25 (0.06) | 625 (60.87) | 1252 (148.15) | 90 (9.25) | 175 (24.13) | 426 (85.57) | 5.26 (0.29) | 2.65 (0.14) | 252 (34.4) | 512 (89.23) | 1285 (310.6) | 3190 (636.58) | 2.8 (0.19)