# Stacking as Accelerated Gradient Descent

Naman Agarwal
Google DeepMind
<EMAIL_ADDRESS>
Pranjal Awasthi
Google Research
<EMAIL_ADDRESS>Satyen Kale
Google Research
<EMAIL_ADDRESS>Eric Zhao
UC Berkeley, Google Research
<EMAIL_ADDRESS>
###### Abstract

Stacking, a heuristic technique for training deep residual networks by
progressively increasing the number of layers and initializing new layers by
copying parameters from older layers, has proven quite successful in improving
the efficiency of training deep neural networks. In this paper, we propose a
theoretical explanation for the efficacy of stacking: viz., stacking
implements a form of Nesterov’s accelerated gradient descent. The theory also
covers simpler models such as the additive ensembles constructed in boosting
methods, and provides an explanation for a similar widely-used practical
heuristic for initializing the new classifier in each round of boosting. We
also prove that for certain deep linear residual networks, stacking does
provide accelerated training, via a new potential function analysis of the
Nesterov’s accelerated gradient method which allows errors in updates. We
conduct proof-of-concept experiments to validate our theory as well.

##  1 Introduction

Deep learning architectures are ubiquitous today and have been responsible
tremendous technological advances in machine learning. However, until 2006,
training deep architectures was extremely challenging. The deep learning
revolution of the past couple of decades was ushered in with the discovery
that a classical technique, viz. _greedy layer-wise pretraining_ , can be used
to train general deep architectures (see (Goodfellow et al., 2016, Section
15.1) for a historical account of this). Previously, only deep architectures
with special structure like convolutions or recurrences were known to be
feasible to train. Greedy layer-wise pretraining is a very intuitive technique
where a deep network is built and trained in a stagewise manner. Starting with
a small network that is easy to train, the technique prescribes adding new
layers over a number of stages, and in each stage training the newly added
layers (and, potentially, the older layers as well) for a certain number of
training steps. This process continues until the desired model depth is
reached.

More modern developments such as residual connections He et al. (2016) and
normalization layers Ioffe & Szegedy (2015) have made it possible to directly
train deep networks without using greedy layer-wise pretraining. However, in
recent years, the tremendous success of deep learning architectures based on
transformers (Vaswani et al., 2017) in domains such as language modeling, and
computer vision (Radford et al., 2019; Brown et al., 2020; Chen et al., 2022;
Tu et al., 2023) has led to the trend of scaling model capacity with ever
increasing model sizes for improved performance (Chowdhery et al., 2022;
Achiam et al., 2023; Team et al., 2023). This endeavour comes at a significant
cost as model training may often take months and require several million
dollars of compute resources (Chowdhery et al., 2022). As a result there has
been a surge of recent work aimed at faster training of large transformer
based models. These works include methods for sparse training such as mixture
of experts (MoE) models (Shazeer et al., 2017; Jiang et al., 2024), methods
for approximate sparse attention mechanisms (Beltagy et al., 2020) and better
optimization methods (Shazeer & Stern, 2018; Liu et al., 2023; Gupta et al.,
2018).

In the effort to reduce the massive costs of training these giant transformer
models, greedy layer-wise pretraining has re-emerged as a very effective
strategy in recent times. Specifically, a technique for initializing the new
layers known as _stacking_ Gong et al. (2019); Reddi et al. (2023) has been
shown to be very effective in speeding up training of deep transformer models.
Stacking prescribes a heuristic for _initializing_ the newly added layers.
Specifically, it prescribes that the newly added layers should be initialized
by copying parameters from the previously trained layers. In the original work
of Gong et al. (2019), the authors proposed to double the model depth at each
time by stacking an exact copy of the current model on top. The recent work of
Reddi et al. (2023) argues that doubling the model depth may be suboptimal and
a better approach is gradual stacking where a few new layers (say 3-4) are
added during each stage. These layers are initialized by copying the top most
layers from the existing model. See Figure 1 for an illustration of this
technique for training deep transformer models.

Figure 1: Stacking for stagewise training language models. In each stage, a
new transformer block is added, initialized with the parameters of the top
block from the previous stage, and then trained for a certain number of steps.

The classical greedy layer-wise pretraining strategy doesn’t have a specific
prescription for initializing the new layers. In general, they’re initialized
randomly in some standard fashion. Stacking initialization provides a clear
benefit over random initialization: Figure 2 shows one example of this effect,
for training the BERT Base Devlin et al. (2019) model with 4 stages of
stagewise training.

Figure 2: Stacking init vs random init for stagewise training of BERT Base
model. Four stages are used with 168,750 steps in each stage. Stage boundaries
are marked by vertical dashed lines. Stacking init provides a clear benefit
over random init.

Structurally, greedy layer-wise pretraining resembles another classical
technique, viz. _boosting_. In boosting, an additive ensemble of classifiers
is constructed via greedy stagewise training of classifiers in the same greedy
manner. Boosting algorithms such as AdaBoost Freund & Schapire (1997) and
Gradient Boosting Friedman (2001) have found tremendous practical application,
especially when using decision trees as base classifiers (e.g. XGBoost Chen &
Guestrin (2016)). A heuristic similar to stacking has also found practical
application in boosting algorithms. The heuristic is to initialize each new
classifier (e.g. a decision tree) by copying over the just-trained classifier
and then updating it using new training data. This process is illustrated in
Figure 3. Due to the similarity with stacking for training deep transformer
models, in the rest of the paper we use “stacking” to also refer to this
initialization strategy in the context of boosting.

Figure 3: Stacking for boosting. In each stage, a new classifier is added,
initialized with the parameters of the last trained classifier from the
previous stage, and then trained for a certain number of steps.

While stacking based methods lead to impressive speed up in training of
transformer models and additive models in boosting, we currently do not have a
good understanding of why this is the case. The recent work of Reddi et al.
(2023) provided a theoretical explanation based on the assumption that each
transformer block is a good few-shot learner. This assumption, along with a
few others, are then used to conclude copying parameters as stacking does
leads to fast learning. However, the assumptions made in the paper are fairly
strong and hard to verify.

In this work, we make progress on understanding the mechanism for the efficacy
of stacking by studying it from an optimization perspective. In particular,
our main contribution is that when viewed from the perspective of function
optimization, stacking speeds up stagewise training by enabling a form of the
_accelerated_ gradient descent method (AGD) developed by Nesterov (Nesterov,
1983). In other words, each stage of the stacking-initialized stagewise
training procedure will reduce the training loss at an accelerated rate.

In contrast, we also show that without using any form of initialization, or in
other words, initializing the new block/classifier to implement the zero
function, stagewise training simply recovers usual (non-accelerated) gradient
descent, whereas random initialization recovers stochastic gradient descent on
a smoothed version of the loss. Hence, stacking initialization accelerates
stagewise training over zero or random initialization. In more detail, our
contributions are as follows:

  * •

We propose a general theoretical framework towards learning a prediction
function $F$ via an ensemble, i.e., a sequence of functions
$(f_{1},f_{2},\ldots,f_{T})$ in a greedy stagewise manner. The generality of
our framework lets us unify classical approaches such as boosting (Freund &
Schapire, 1997; Friedman, 2001) that build the ensemble in an additive manner,
and modern approaches that build the ensemble via stagewise training of
_residual_ function compositions (e.g. ResNets He et al. (2016) and
Transformer models Vaswani et al. (2017)).

  * •

Our proposed framework lets us formally establish the connection between
various initialization strategies used for building the ensemble and the
convergence properties of the resulting overall learning procedure. In
particular, we show that the zero initialization strategy recovers the vanilla
functional gradient descent algorithm, for both the additive (i.e. boosting)
and residual compositional forms of learning, whereas random initialization
recovers _stochastic_ functional gradient descent (on a smoothed loss) for
both types of models. Furthermore, in the case of additive models, the use of
the popular stacking initialization exactly recovers Nesterov’s accelerated
functional gradient descent. The consequence is that for $T$ stages of
boosting with stacking initialization, loss reduces at a rate of $O(T^{-2})$
for smooth losses, or $\exp(-\Omega(T/\sqrt{\kappa}))$ for smooth and
strongly-convex losses with condition number $\sqrt{\kappa}$, as opposed to
rates of $O(T^{-1})$ and $\exp(-\Omega(T/\kappa))$ respectively for zero
initialization.

  * •

For the case of compositional updates we show that stacking initialization
results in updates that look remarkably similar to Nesterov’s accelerated
functional gradient descent. Proving an accelerated rate in the general non-
parametric functional setting seems intractable, so we analyze stacking in a
special parametric setting of deep linear networks with a convex loss
function. In this setting we prove that the stacking initialization
quantitatively leads to the same kind of convergence benefits over vanilla
gradient descent as is observed for Nesterov’s accelerated method. At the core
of our proof is a novel potential function based analysis of Nesterov’s method
with errors in the momentum term that may be of independent interest.

  * •

Finally, we also perform proof-of-concept experiments to validate our theory
on synthetic and real world data.

###  1.1 Related work

Boosting is a classical technique for constructing additive ensembles via
greedy stagewise training, and has a long and rich history of work. We refer
the interested reader to the excellent textbook of Schapire & Freund (2012)
for the literature on this topic.

The idea of training deep residual networks in a layer wise manner has been
explored in many prior works. In earlier studies (Hinton et al., 2006; Bengio
et al., 2006) the focus was on greedily adding trained layers to the model
while keeping the bottom layers frozen followed by a final fine tuning step
where the entire network is trained. In recent years progressive or gradual
stacking (Gong et al., 2019; Gu et al., 2020; Shen et al., 2022; Reddi et al.,
2023) has emerged as a powerful way to train deep networks especially
transformer based architectures.

The empirical insight in the work of (Gong et al., 2019) was that the
attention patterns in neighboring layers of trained transformer models show
remarkable similarity. Hence, by copying the parameters from the previous
layer one is providing a better initialization for the optimization procedure.
As mentioned previously, Reddi et al. (2023) developed the gradual stacking
approach based on the assumption that the trained transformer blocks are good
few-shot learners, and showed that gradual stacking leads to significant
wallclock improvements during training.

##  2 Stagewise training as functional gradient descent

#### Preliminaries.

We consider a fairly general supervised learning setting. Denote the input
space by $\mathcal{X}$ and the output space by $\mathcal{Y}$. Examples
$(x,y)\in\mathcal{X}\times\mathcal{Y}$ are drawn from a distribution $D$
(which may simply be the empirical data distribution in the case of empirical
risk minimization). We aim to model the input-output relationship via
_predictions_ in $\mathbb{R}^{d}$, for some dimension parameter $d$. Given an
example $(x,y)\in\mathcal{X}\times\mathcal{Y}$ the quality of a prediction
$\widehat{y}\in\mathbb{R}^{d}$ is measured via a loss function
$\ell:\mathbb{R}^{d}\times\mathcal{Y}\to\mathbb{R}$. Predictions are computed
using functions $f:\mathcal{X}\to\mathbb{R}^{d}$. We will assume that the
predictor functions $f$ are square integrable with respect to $D$, i.e.
$\operatorname*{\mathbb{E}}_{(x,y)\sim D}[\|f(x)\|^{2}]<\infty$. The space of
such functions forms a Hilbert space, denoted $\mathcal{L}_{2}$, with the
inner product defined as $\langle
f,g\rangle=\operatorname*{\mathbb{E}}_{(x,y)\sim D}[\langle
f(x),g(x)\rangle]$. Unless specified otherwise, all functions in the
subsequent discussion will be assumed to be in $\mathcal{L}_{2}$. The loss
function can then naturally be extended to predictor functions
$f\in\mathcal{L}_{2}$ by defining, with some abuse of notation,
$\ell(f):=\operatorname*{\mathbb{E}}_{(x,y)\sim D}[\ell(f(x),y)]$. The goal of
training is to obtain a function $f\in\mathcal{L}_{2}$ that minimizes
$\ell(f)$.

In the rest of this section, we perform the analysis in a purely functional
setting, which affords a convenient analysis. However, we note that, in
practice, functions are parameterized (say by neural networks) and hence
update rules for functions may not always be realizable via the specific
parameterization used. We ignore these realizability issues in the following
analysis in order to focus on the conceptual message that stacking
initialization enables accelerated updates.

We now define a general ensemble learning setup within the above setting. In
this setup, we aim to approximate the minimizer of $\ell$ on $\mathcal{L}_{2}$
via an _ensemble_ , which is a sequence of functions
$(f_{1},f_{2},\ldots,f_{T})$, where $T>0$ is a given parameter defining the
size of the ensemble. The functions in the ensemble are typically “simple” in
the sense that they are chosen from a class of functions that is easy to
optimize over. A predictor function can be obtained from an ensemble
$(f_{1},f_{2},\dots,f_{T})$ by aggregating its constituent functions into a
single function $F_{T}:\mathcal{X}\to\mathbb{R}^{d}$. The loss of an ensemble
can then be defined (again with some abuse of notation) in terms of its
aggregation as $\ell((f_{1},f_{2},\ldots,f_{T})):=\ell(F_{T})$. Two specific
aggregation operators we will consider are the following:

  1. 1.

Addition: Used in boosting, this results in a summation over ensemble outputs:
$F_{T}=f_{1}+f_{2}+\cdots f_{T}$.

  2. 2.

Residual composition: Used in deep residual neural networks, this results in
the composed function
$F_{T}=(I+f_{T})\circ(I+f_{T-1})\circ\cdots\circ(I+f_{1})$, where the domain
is $\mathcal{X}=\mathbb{R}^{d}$ and $I:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is the
identity mapping.

#### Greedy stagewise training.

_Stagewise training_ is a simple greedy procedure to train ensembles in a
progressive manner. Suppose we have already obtained a (partial) ensemble
$(f_{1},f_{2},\ldots,f_{t})$. Then, the next function in the ensemble,
$f_{t+1}$, is _ideally_ obtained by minimizing the loss of the new ensemble,
i.e. $f_{t+1}=\arg\min_{f}\ell((f_{1},f_{2},\ldots,f_{t},f))$.

However, in practice, this ideal is hard to implement, and instead two
heuristics are commonly used: (a) the new function to be trained is
_initialized_ in some carefully chosen manner, and (b) the optimization above
is done using _early stopping_ , i.e. a few steps of gradient descent, which
ensures that the new function stays close to initialization. To analyze these
practical considerations in a theoretically convenient framework, we model
them in a functional optimization setting as follows.

First, we assume that the function $f_{t+1}$ to be trained is initialized at
some carefully chosen value $f_{t+1}^{0}$. For notational convenience, we
denote the aggregation of the ensemble
$(f_{1},f_{2},\ldots,f_{t},f_{t+1}^{0})$ by $F_{t+1}^{0}$ and that of the
generic ensemble $(f_{1},f_{2},\ldots,f_{t},f)$ by $F$.

Next, we note that early stopping acts as a form of $\ell_{2}$ regularization
(see, e.g. (Goodfellow et al., 2016, Section 7.8)) which ensures that the
trained function $f_{t+1}$ remains close to its initialization, $f_{t+1}^{0}$.
Based on this intuition, for ease of theoretical analysis, we approximate
$\ell(F)$ by its linearization around the initialization:
$\ell(F)\approx\ell(F_{t+1}^{0})+\langle\nabla\ell(F_{t+1}^{0}),F-F_{t+1}^{0}\rangle$.
Here, $\nabla\ell(\cdot)$ is the Fréchet derivative, and
$\langle\cdot,\cdot\rangle$ denotes the inner product in $\mathcal{L}_{2}$.
Thus, $f_{t+1}$ is obtained by minimizing

$\displaystyle\ell(F_{t+1}^{0})+\langle\nabla\ell(F_{t+1}^{0}),F-F_{t+1}^{0}\rangle+\frac{\lambda}{2}\|F-F_{t+1}^{0}\|^{2}.$
(1)

While a characterization of $f_{t+1}$ may not be apparent from the above
discussion, in the cases of interest (additive and residual compositional
models) one can get useful insights via minimizing the above expression over
$F$ to directly obtain $F_{t+1}$:

$F_{t+1}=F_{t+1}^{0}-\frac{1}{\lambda}\nabla\ell(F_{t+1}^{0}).$ (2)

We can now consider specific initialization strategies (i.e. zero
initialization, random initialization, and stacking initialization) in the
context of additive and residual compositional models and see how these
initializations lead to various forms of functional gradient descent.

#### Stagewise training with zero initialization recovers functional gradient
descent.

First, consider stagewise training where functions are initialized to be zero
functions, i.e. $f_{t+1}^{0}=0$. It is easy to see that with this
initialization, for both additive and residual compositional models, we have
$F_{t+1}^{0}=F_{t}$. Thus, from (2), we have that the updated ensemble’s
predictor can be written as

$F_{t+1}=F_{t}-\frac{1}{\lambda}\nabla\ell(F_{t}).$

This exactly describes functional gradient descent with step size
$\frac{1}{\lambda}$. In the additive setting this is well-known: indeed,
boosting can be seen as a functional gradient descent Mason et al. (1999). The
result for the residual compositional setting appears to be new.

#### Stagewise training with random initialization recovers stochastic
functional gradient descent on smoothed loss.

We now consider stagewise training where functions are initialized randomly,
i.e. $f_{t+1}^{0}$ is a randomly drawn function, independent of all randomness
up to stage $t$. In the following, we will assume that
$\operatorname*{\mathbb{E}}[f_{t+1}^{0}]=0$, where the $0$ on the RHS denotes
the zero function. With this initialization, for both additive and residual
compositional models, we have $F_{t+1}^{0}=F_{t}+g_{t}$, where
$g_{t}=f_{t+1}^{0}$ for additive models, and $g_{t}=f_{t+1}^{0}\circ F_{t}$
for residual compositional models. In either case, note that
$\operatorname*{\mathbb{E}}[g_{t}]=0$. Now define the loss functional
$\ell_{t}(F):=\operatorname*{\mathbb{E}}[\ell(F+g_{t})]$. Since
$\operatorname*{\mathbb{E}}[g_{t}]=0$, we can interpret $\ell_{t}$ as a
randomized smoothing of $\ell$, similar to convolving with a Gaussian. Then,
from (2), we have that the updated ensemble’s predictor can be written as

$F_{t+1}=F_{t}-\frac{1}{\lambda}(\nabla\ell(F_{t}+g_{t})-\lambda g_{t}).$

Now, note that

$\operatorname*{\mathbb{E}}[\nabla\ell(F_{t}+g_{t})-\lambda
g_{t}]=\nabla\ell_{t}(F_{t}).$

Or in other words, the above update can be seen as a _stochastic_ functional
gradient descent step on the smoothed loss function $\ell_{t}$.

####  Stagewise training with stacking initialization recovers accelerated
functional gradient descent.

We now consider stagewise training where functions are initialized in a
stacking-like fashion with $f_{t+1}^{0}=f_{t}$, which we will refer to as the
_stacking initialization_. When the ensemble aggregation operator is addition,
we have $F_{t+1}^{0}=f_{t}+F_{t}=F_{t}+(F_{t}-F_{t-1})$ and hence (2) implies
that the updated ensemble’s predictor is

$F_{t+1}=F_{t}+(F_{t}-F_{t-1})-\frac{1}{\lambda}\nabla\ell(F_{t}+(F_{t}-F_{t-1})).$

The above formula essentially describes Nesterov’s accelerated gradient
descent, which has the following update rule:

$F_{t+1}=F_{t}+\beta(F_{t}-F_{t-1})-\frac{1}{\lambda}\nabla\ell(F_{t}+\beta(F_{t}-F_{t-1})).$
(3)

Here, $\beta\in[0,1)$ is a constant that can depend on $t$. In fact, we can
_exactly_ recover Nesterov’s accelerated gradient descent if we modify the
stacking initialization to $f_{t+1}^{0}=\beta f_{t}$. Thus, stacking enables
accelerated descent for training additive models.

When the ensemble aggregation operator is residual composition, stagewise
training with the stacking initialization $f_{t+1}^{0}=f_{t}$ results in

$F_{t+1}^{0}=(I+f_{t})\circ F_{t}=F_{t}+f_{t}\circ F_{t}.$

(2) therefore implies the updated ensemble’s predictor is

$F_{t+1}=F_{t}+f_{t}\circ F_{t}-\frac{1}{\lambda}\nabla\ell(F_{t}+f_{t}\circ
F_{t}).$ (4)

In contrast, Nesterov’s update rule (3), and the fact that for residual
compositional models,

$F_{t}-F_{t-1}=(I+f_{t})\circ F_{t-1}=f_{t}\circ F_{t-1},$

yields the following equation for $F_{t+1}$:

$F_{t+1}=F_{t}+\beta f_{t}\circ
F_{t-1}-\frac{1}{\lambda}\nabla\ell(F_{t}+\beta f_{t}\circ F_{t-1}).$ (5)

Comparing (4) and (5), barring the minor difference of the $\beta$ parameters,
which can be easily rectified exactly as in the case of the additive models by
setting $f_{t+1}^{0}=\beta f_{t}$, the major difference is that $f_{t}\circ
F_{t-1}$ is replaced by $f_{t}\circ F_{t}$. Although possibly intractable to
prove formally, we believe that the updates in (4) also provide an accelerated
convergence rate, since we expect $F_{t-1}$ to be close to $F_{t}$ as iterates
converge to the optimal function. In the following section, we show that in
certain deep linear networks, the above intuition is indeed correct and
stacking provides an accelerated convergence rate.

##  3 Accelerated convergence of deep linear networks by stacking

To demonstrate that stacking can provide a provably accelerated rate of
convergence, we now turn to studying the narrower setting of training deep
residual linear networks, which are fully connected feedforward neural
networks without non-linear activations and residual connections. Such
networks are a common subject of study in the theory of deep learning Saxe et
al. (2014); Kawaguchi (2016); Hardt & Ma (2017). As they have no non-linear
components, deep linear networks effectively compute a linear functions,
albeit via a parametrization as a product of the weight matrices.

#### Setup.

Consider again the general supervised learning setting from Section 2 and
suppose, as is often the case in modern neural networks, that examples consist
of inputs $x\in\mathbb{R}^{d}$ and outputs $y\in\mathcal{Y}$. The loss
function $\ell:\mathbb{R}^{d}\times\mathcal{Y}\to\mathbb{R}$ is assumed to be
convex in the first argument. Let the samples be drawn from a distribution $D$
over $\mathbb{R}^{d}\times\mathcal{Y}$. Then the expected loss of the linear
predictor $x\mapsto Wx$ for a matrix $W\in\mathbb{R}^{d\times d}$ is (with
some abuse of notation) $\ell(W):=\operatorname*{\mathbb{E}}_{(x,y)\sim
D}[\ell(Wx,y)]$. In the following, we will assume the expected loss $\ell(W)$
is $L$-smooth and $\mu$-strongly convex in $W$, by which we mean that the
following inequalities hold for any $W,V\in\mathbb{R}^{d\times d}$:

$\displaystyle\ell(W)+\langle\nabla\ell(W),V-W\rangle+\frac{\mu}{2}\left\lVert
W-V\right\rVert^{2}$
$\displaystyle\leq\ell(V)\leq\ell(W)+\langle\nabla\ell(W),V-W\rangle+\frac{L}{2}\left\lVert
W-V\right\rVert^{2}.$

Here, for matrices $W,V\in\mathbb{R}^{d}$, $\langle
W,V\rangle=\text{Tr}(W^{\top}V)$, and $\left\lVert W\right\rVert$ is the
Frobenius norm of $W$.

The _condition number_ $\kappa$ of the loss is defined as
$\kappa:=\frac{L}{\mu}$.

The deep residual neural networks we consider have $t$ layers with weight
matrices $w_{1},w_{2},\ldots,w_{t}$, and the function they compute is
$x\mapsto W_{t}x$, where

$W_{t}:=(I+w_{t})(I+w_{t-1})\dots(I+w_{1}).$

Here, $I\in\mathbb{R}^{d\times d}$ is the identity matrix providing the
residual connection. The expected loss of the neural network described above
on the data is $\ell(W_{t})$.

#### Derivation of stacking updates.

Suppose we train the deep residual linear network described above using
stacking initialization, but incorporating $\beta$-scaling: i.e., to train the
$(t+1)$-th layer, its weight matrix is initialized to $w_{t+1}^{0}=\beta
w_{t}$, for some constant $\beta\in[0,1]$, and then trained. Following the
exact same steps as in the derivation of stacking updates in the functional
setting of Section 2, we end up with the following formula for $W_{t+1}$:

$W_{t+1}=W_{t}+\beta w_{t}W_{t}-\frac{1}{\lambda}\nabla\ell(W_{t}+\beta
w_{t}W_{t}).$

When $W_{t-1}$ is non-singular, we have
$w_{t}=W_{t}W_{t-1}^{-1}-I=(W_{t}-W_{t-1})W_{t-1}^{-1}$, so the above equation
can be rewritten as

$\displaystyle W_{t+1}=$
$\displaystyle\;W_{t}+\beta(W_{t}-W_{t-1})W_{t-1}^{-1}W_{t}$
$\displaystyle-\frac{1}{\lambda}\nabla\ell(W_{t}+\beta(W_{t}-W_{t-1})W_{t-1}^{-1}W_{t}).$
(6)

As previously noted, (6) differs from Nesterov’s AGD method in form:
Nesterov’s AGD updates would be

$\displaystyle W_{t+1}$ $\displaystyle=W_{t}+\beta(W_{t}-W_{t-1})$
$\displaystyle\;-\frac{1}{\lambda}\nabla\ell(W_{t}+\beta(W_{t}-W_{t-1})).$ (7)

#### Accelerated convergence for stacking updates.

Despite the differences with Nesterov’s method, we can show that stacking
notably still yields a provably accelerated convergence rate. Let
$W^{*}:=\arg\min_{W}\ell(W)$. Suppose that $W^{*}$ is non-singular, i.e. its
smallest singular value $\sigma_{\min}(W^{*})>0$. Theorem 3.1 shows that as
long as the first two layers are initialized so that $W_{1}$ and $W_{2}$ are
close to optimal, stacking results in a suboptimality gap of
$\exp(-\widetilde{\Omega}(T/\sqrt{\kappa}))$ after $T$ stages of stacking.
This is of the same order as the rate obtained by Nesterov’s acceleration;
note that, in comparison, stagewise training with zero initialization results
would result in a suboptimality gap of $\exp(-\widetilde{\Omega}(T/\kappa))$.
All missing proofs are in the appendix.

######  Theorem 3.1.

Consider stagewise training with stacking initialization of a deep residual
linear network in the setup described above with
$\beta=\tfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$ and $\lambda=\tfrac{1}{L}$.
If the first two layers $w_{1},w_{2}$ are initialized so that
$W_{1},W_{2}\in\mathbb{R}^{d\times d}$ satisfy
$\max\\{\ell(W_{1}),\ell(W_{2})\\}-\ell(W^{*})\in
O(\tfrac{\mu^{5}}{dL{}^{4}}\sigma_{\min}(W^{*})^{2})$, then after $T$ stages
of stacking, we have

$\ell(W_{T})-\ell(W^{*})\leq\exp(-\widetilde{\Omega}(T/\sqrt{\kappa})).$

###### Proof.

We will show that $\ell(W_{T})-\ell(W^{*})\leq O(\sqrt{\kappa}\left\lVert
W_{0}-W^{*}\right\rVert\exp(-\Omega(T/\sqrt{\kappa})))$ whenever the
assumption

$\max\left\\{\ell(W_{0}),\ell(W_{1})\right\\}-\ell(W^{*})\leq\tfrac{4\sqrt{\kappa}-3}{\tfrac{4}{\mu}(4\sqrt{\kappa}-3)+1}(\tfrac{\alpha}{2\beta\sqrt{d}+\alpha}\sigma_{\min}(W^{*}))^{2}$

holds. One can rewrite this assumption more simply as
$\ell(W_{0})-\ell(W^{*})\leq\tfrac{\delta}{2(1+\rho)}$ and
$\ell(W_{1})-\ell(W^{*})\leq\tfrac{\delta}{2(1+\rho)}$, where we use $\rho$ as
defined in (12) and the shorthand
$\delta\coloneqq\tfrac{\mu}{2}(\tfrac{\alpha}{2\beta\sqrt{d}+\alpha}\sigma_{\min}(W^{*}))^{2}$.
We first observe that our deep linear networks $W_{1},\dots,W_{T}$, as defined
in (6), can be written as iterates of a perturbed version of Nesterov’s
acceleration in the form of (10) where the perturbation term $\Delta_{t}$ is
as defined in (22). This also means we set the gradient step size $\lambda$ to
match the usual step size used in Nesterov’s gradient method with
$\lambda=\tfrac{1}{L}$ and the momentum parameter $\beta$ to match its usual
setting of $\tfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$.

We will use induction to prove the following statement inductively: we have
$\left\lVert\Delta_{t}\right\rVert\leq\alpha\left\lVert
W_{t}-W_{t-1}\right\rVert$ for all $t\in[T-1]$ and $\Phi(t)\leq\delta$ for all
$t\in[T]$, where the potential functions $\Phi$ are as defined in (12). Once
we have shown this statement, the conditions for Lemma 3.5 are satisfied and
we immediately have that

$\displaystyle\ell(W_{t})-\ell(W^{*})\leq(\tfrac{2\sqrt{\kappa}-2}{2\sqrt{\kappa}-1})^{T}\left(L{}+\tfrac{8\sqrt{\kappa}\mu}{4(4\sqrt{\kappa}-3)}\right)\left\lVert
W_{0}-W^{*}\right\rVert^{2}.$

Using the shorthand
$C=\left(L{}+\tfrac{8\sqrt{\kappa}\mu}{4(4\sqrt{\kappa}-3)}\right)\left\lVert
W_{0}-W^{*}\right\rVert^{2}$, we can recover the original form of the
Theorem’s claim by rewriting the above in binomial form as the power series of
$\exp(-\Omega(\tfrac{T}{\sqrt{\kappa}}))$:

$\displaystyle\tfrac{1}{C}(\ell(W_{t})-\ell(W^{*}))$
$\displaystyle\leq(1-\tfrac{1}{2\sqrt{\kappa}-1})^{T}$
$\displaystyle=1-\tfrac{T}{2\sqrt{\kappa}-1}+\tfrac{T(T-1)}{2!}\tfrac{1}{(2\sqrt{\kappa}-1)^{2}}-\tfrac{T(T-2)(T-3)}{3!}\tfrac{1}{(2\sqrt{\kappa}-1)^{2}}+\dots$
$\displaystyle\in
1-\Omega(\tfrac{T}{\sqrt{\kappa}})+\tfrac{1}{2!}\Omega(\tfrac{T}{\sqrt{\kappa}})^{2}-\tfrac{1}{3!}\Omega(\tfrac{T}{\sqrt{\kappa}})^{3}+\dots$
$\displaystyle=\exp(-\Omega(\tfrac{T}{\sqrt{\kappa}})).$

Holding $\left\lVert W_{0}-W^{*}\right\rVert^{2}$, we note
$C\in\exp(-\Omega(\log(\tfrac{T}{\sqrt{\kappa}})))$ and thus
$\ell(W_{t})-\ell(W^{*})\in\exp(-\Omega(\tfrac{T}{\sqrt{\kappa}}+\log(\tfrac{T}{\sqrt{\kappa}})))$
as claimed.

#### Base case.

We first show the base case that $\Phi(1)\leq\delta,\Phi(0)\leq\delta$. We can
use the strong convexity of $\ell$ to upper bound each of these initial
potentials $\Phi(0),\Phi(1)$ by

$\displaystyle\Phi(0)$
$\displaystyle=\ell(W_{0})-\ell(W^{*})+(\tfrac{\mu}{2}+2\rho)\left\lVert
W_{0}-W^{*}\right\rVert^{2}\leq 2(1+\rho)(\ell(W_{0})-\ell(W^{*}))$ (8)
$\displaystyle\Phi(1)$
$\displaystyle=\ell(W_{1})-\ell(W^{*})+(\tfrac{\mu}{2}+2\rho)\left\lVert
W_{1}-W^{*}\right\rVert^{2}\leq 2(1+\rho)(\ell(W_{1})-\ell(W^{*})).$ (9)

Since we assume that $\ell(W_{0})\leq\ell(W^{*})+\tfrac{\delta}{2(1+\rho)}$
and $\ell(W_{1})\leq\ell(W^{*})+\tfrac{\delta}{2(1+\rho)}$, we thus have
$\Phi(0),\Phi(1)\leq\delta$.

#### Inductive step.

Fixing a $t\in[T-1]$, we now aim to prove that
$\left\lVert\Delta_{t}\right\rVert\leq\alpha\left\lVert
W_{t}-W_{t-1}\right\rVert$ and $\Phi(t+1)\leq\Phi(0)$. We can assume, as our
inductive hypothesis, that $\Phi(t)\leq\Phi(0)\leq\delta$, implying by strong
convexity that

$\delta\geq\Phi(t)\geq(1+\tfrac{\gamma}{2})^{t}(\ell(W_{t})-\ell(W^{*}))\geq(1+\tfrac{\gamma}{2})^{t}\tfrac{\mu}{2}\left\lVert
W_{t}-W^{*}\right\rVert^{2}.$

We will use the shorthand
$\eta_{t}\coloneqq\sqrt{\frac{2\delta}{\mu}(1+\tfrac{\gamma}{2})^{-t}}$ so
that we can write $\left\lVert W_{t}-W^{*}\right\rVert\leq\eta_{t}$. We can
similarly assume that $\Phi(t-1)\leq\Phi(0)$ and, by extension, $\left\lVert
W_{t-1}-W^{*}\right\rVert\leq\eta_{t-1}$. Finally, we can also assume that for
all $\tau\in[t-1]$ $\left\lVert\Delta_{\tau}\right\rVert\leq\alpha\left\lVert
W_{\tau}-W_{\tau-1}\right\rVert$ and therefore $\left\lVert
W_{\tau}W_{\tau-1}^{-1}-I\right\rVert\leq\frac{1}{\beta}\alpha$. This final
assumption means that, once we have shown the first half of our inductive
claim that $\left\lVert\Delta_{t}\right\rVert\leq\alpha\left\lVert
W_{t}-W_{t-1}\right\rVert$, then Lemma 3.5 directly implies the second part of
the inductive claim, i.e. that $\Phi(t+1)\leq\Phi(0)$. Thus, it only remains
to prove Claim 3.2.

######  Claim 3.2.

For any $t\in[T-1]$, if we assume that (1)
$\left\lVert\Delta_{\tau}\right\rVert\leq\alpha\left\lVert
W_{\tau}-W_{\tau-1}\right\rVert$ for all $\tau\in[t-1]$, (2)
$\Phi(t)\leq\Phi(0)$ and (3) $\Phi(t-1)\leq\Phi(0)$, then we can conclude that
$\left\lVert\Delta_{t}\right\rVert\leq\alpha\left\lVert
W_{t}-W_{t-1}\right\rVert$.

We begin by establishing that $W_{t-1}$ is indeed invertible and hence
$\Delta_{t}$ is well-defined.

######  Fact 3.3.

$W_{t-1}$ is invertible. Moreover, the smallest singular value of $W_{t-1}$ is
bounded by $\left\lVert
W_{t-1}\right\rVert^{-1}\leq\frac{\sqrt{d}}{\sigma_{\min{}}(W^{*})-\eta_{t-1}}$.

We can next bound by submultiplicativity the perturbation norm

$\left\lVert\Delta_{t}\right\rVert=\left\lVert\beta(W_{t}-W_{t-1})(W_{t-1}^{-1}W_{t}-I)\right\rVert\leq\beta\left\lVert
W_{t}-W_{t-1}\right\rVert\left\lVert W_{t-1}^{-1}W_{t}-I\right\rVert.$

We can further simplify right-most factor as follows:

$\displaystyle\left\lVert W_{t-1}^{-1}W_{t}-I\right\rVert$
$\displaystyle=\left\lVert W_{t-1}^{-1}(W_{t}-W_{t-1})\right\rVert$
$\displaystyle\leq\left\lVert W_{t}-W_{t-1}\right\rVert\left\lVert
W_{t-1}^{-1}\right\rVert$ (submultiplicative property)
$\displaystyle\leq(\left\lVert W_{t}-W^{*}\right\rVert+\left\lVert
W_{t-1}-W^{*}\right\rVert)\left\lVert W_{t-1}^{-1}\right\rVert$ (triangle
inequality)
$\displaystyle\leq(\eta_{t}+\eta_{t-1})\sqrt{d}\frac{1}{\sigma_{\min}(W^{*})-\eta_{t-1}}$
$\displaystyle\leq 2\eta_{0}\sqrt{d}\frac{1}{\sigma_{\min}(W^{*})-\eta_{0}}.$

Here, the second-to-last inequality applies the inductive hypothesis and Fact
3.3. The last inequality applies that $\eta_{1},\dots,\eta_{T}$ is a
monotonically decreasing sequence. We conclude by observing that, since we
chose
$\delta\leq\tfrac{\mu}{2}(\tfrac{\alpha}{2\beta\sqrt{d}+\alpha}\sigma_{\min}(W^{*}))^{2}$,
we guaranteed
$\eta_{0}\leq\tfrac{\alpha}{2\beta\sqrt{d}+\alpha}\sigma_{\min}(W^{*})$ and by
extension
$2\eta_{0}\sqrt{d}\frac{1}{\sigma_{\min}(W^{*})-\eta_{0}}\leq\tfrac{\alpha}{\beta}$.
This gives, as desired, that
$\left\lVert\Delta_{t}\right\rVert\leq\alpha\left\lVert
W_{t}-W_{t-1}\right\rVert$, completing our inductive step and proof. We now
turn to proving the previously stated Fact 3.3.

###### Proof of Fact 3.3.

Recall that our inductive hypothesis gives that
$\|W^{*}-W_{t-1}\|\leq\eta_{t-1}$. Since the Frobenius norm is just a norm of
singular values, we can upper bound the singular values of the matrix
$W^{*}-W_{t-1}$ by $\sigma_{\max{}}(W^{*}-W_{t-1})\leq\eta_{t-1}$, where we
use $\sigma_{\max{}}(W)$ to denote the largest singular values of a matrix
$W$. We can therefore invoke Weil’s inequality to argue that, for any
$i\in[d]$, the difference between the $i$th largest singular values of
$W_{t-1}$ and $W^{*}$ is upper bounded by

$\left|\sigma_{i}(W_{t-1})-\sigma_{i}(W^{*})\right|\leq\sigma_{\max}(W_{t-1}-W^{*})\leq\eta_{t-1}.$

The smallest singular value of $W_{t-1}$ is thus at least
$\sigma_{\min{}}(W_{t-1})\geq\sigma_{\min{}}(W^{*})-\eta_{t-1}$. Since our
choice of $\delta$ guarantees that $\eta_{t-1}\leq\sigma_{\min{}}(W^{*})$, we
have that $W_{t-1}>0$ and is invertible. This also implies that the largest
singular value of $W_{t-1}^{-1}$ is at most
$\sigma_{\max{}}(W_{t-1}^{-1})\leq\frac{1}{\sigma_{\min{}}(W^{*})-\eta_{t-1}}$.
We can therefore upper bound the Frobenius norm of $W_{t-1}^{-1}$ as claimed
by

$\left\lVert
W_{t-1}^{-1}\right\rVert\leq\sqrt{\sum_{i\in[d]}\sigma_{i}(W_{t-1}^{-1})^{2}}\leq\sqrt{d}\sigma_{\max{}}(W_{t-1}^{-1})\leq\sqrt{d}\frac{1}{\sigma_{\min{}}(W^{*})-\eta_{t-1}}.$

∎

∎

The $\widetilde{\Omega}(\cdot)$ notation above hides polylogarithmic
dependence on the problem parameters for clarity of presentation. Precise
expressions can be found in the proof.

Theorem 3.1 presents a _local_ convergence result: it assumes that the initial
two layers put the network in the vicinity of the optimal solution $W^{*}$. We
can also provide a _global_ convergence result by using a small warmup period
phase to stacking—that is, by training the first few (only a constant number)
stages of the deep linear network without stacking.

######  Corollary 3.4 (Corollary of Theorem 3.1).

Consider stagewise training of a deep residual linear network in the setup
described above with an initial warmup phase of zero initialization for
$\widetilde{O}(\kappa)$ stages followed by stacking initialization for the
remaining stages. Then after $T$ total stages, we have

$\ell(W_{T})-\ell(W^{*})\leq\exp(-\widetilde{\Omega}(T/\sqrt{\kappa})+\widetilde{O}(\sqrt{\kappa})).$

###### Proof.

We first recall the standard convergence rate of gradient descent. In the same
setting as Lemma 3.5 where we minimize a $L$-smooth and $\mu$-strongly convex
function on some Hilbert space $\mathcal{F}$, we can define gradient descent
iterates with the update rule $x_{t+1}=x_{t}-\tfrac{1}{L}\nabla\ell(x_{t})$.
For the resulting sequence of iterates $x_{0},\dots,x_{T}$, it is known that
$\ell(x_{T})-\ell(x^{*})\leq\exp(-T/\kappa)(\ell(x_{0})-\ell(x^{*}))$; see
e.g. Bansal & Gupta (2017).

We can further recall that, in the stagewise training of deep residual linear
networks, we can write networks as
$W_{t+1}=W_{t}+w_{t+1}^{0}W_{t}-\tfrac{1}{\lambda}\nabla\ell(W_{t}+w_{t+1}^{0}W_{t})$.
In the initial stages where new layers are initialized with zero weights, i.e.
$w_{t+1}^{0}=0$, we recover, as mentioned in Section 2, a gradient descent on
$\mathbb{R}^{d\times d}$:

$\displaystyle W_{t+1}=W_{t}-\tfrac{1}{\lambda}\nabla\ell(W_{t}).$

Putting these two pieces together, we have $\max_{t\leq
T_{0}}\ell(W_{t})-\ell(W^{*})\leq\exp(-\tfrac{T_{0}+1}{\kappa})(\ell(W_{0})-\ell(W^{*}))$.
By setting

$T_{0}\geq\kappa\log(\tfrac{dL^{4}}{C\mu^{5}\sigma_{\min}(W^{*})^{2}}(\ell(W_{0})-\ell(W^{*}))),$

the iterates $W_{T_{0}}$ and $W_{T_{0}-1}$ are close enough to the solution
that the local convergence result of Theorem 3.1 applies. Specifically, we
have that $\ell(W_{T_{0}-1})-\ell(W^{*})\leq
C\tfrac{\mu^{5}}{dL^{4}}\sigma_{\min}(W^{*})^{2}$ and
$\ell(W_{T_{0}-1})-\ell(W^{*})\leq
C\tfrac{\mu^{5}}{dL^{4}}\sigma_{\min}(W^{*})^{2}$. Theorem 3.1 thus implies
that performing $T^{\prime}$ rounds of stacking after $T_{0}$ rounds of warm-
start leads to the desired loss bound

$\displaystyle\ell(W_{T_{0}+T^{\prime}})-\ell(W^{*})$
$\displaystyle\leq\exp(-\widetilde{\Omega}(T^{\prime}/\sqrt{\kappa}))$
$\displaystyle=\exp(-\widetilde{\Omega}(T/\sqrt{\kappa})+\widetilde{O}(\sqrt{\kappa})).$

∎

#### Key ingredient of proof: robustness of Nesterov’s accelerated gradient
descent method.

The primary insight behind the result of Theorem 3.1 is that Nesterov’s
accelerated gradient method is relatively robust to perturbations in its
update rules. This robustness is formalized below in Lemma 3.5. The lemma is
described in a fairly general, standalone setting since it may be of
independent interest.

######  Lemma 3.5 (Robustness of Nesterov’s accelerated gradient method).

Let $\mathcal{F}$ be a Hilbert space and $\ell:\mathcal{F}\to\mathbb{R}$ an
$L{}$-smooth and $\mu$-strongly convex function to be minimized on
$\mathcal{F}$. Consider the iterates
$x_{0},y_{0},\dots,x_{T},y_{T}\in\mathcal{F}$ with $x_{0}=y_{0}$ chosen
arbitrarily, and the update rules

$\displaystyle y_{t+1}$ $\displaystyle\coloneqq
x_{t}-\frac{1}{L{}}\nabla\ell(x_{t})$ $\displaystyle x_{t+1}$
$\displaystyle\coloneqq y_{t+1}+\beta(y_{t+1}-y_{t})+\Delta_{t+1}$ (10)

where $\beta:=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$,
$\kappa=\tfrac{L}{\mu}$, and
$\Delta_{1},\Delta_{2},\dots,\Delta_{T-1}\in\mathcal{F}$ are error terms such
that $\left\lVert\Delta_{t}\right\rVert\in O(\kappa^{-2}\left\lVert
y_{t}-y_{t-1}\right\rVert)$ for all $t\in[T-1]$. Then the convergence rate of
the iterates to a suboptimality gap of $\varepsilon$ is of order
$\widetilde{O}(\sqrt{\kappa}\log(\kappa/\varepsilon))$. Specifically, for any
$T\geq 2$,

$\ell(y_{T})-\ell(x^{*})\leq(\tfrac{2\sqrt{\kappa}-2}{2\sqrt{\kappa}-1})^{T}\left(L{}+\tfrac{8\sqrt{\kappa}\mu}{4(4\sqrt{\kappa}-3)}\right)\left\lVert
x_{0}-x^{*}\right\rVert^{2}.$

###### Proof.

In this proof, we will assume that the norm of the perturbation term is
bounded by

$\displaystyle\|\Delta_{t}\|\leq\alpha\|y_{t}-y_{t-1}\|,\quad\text{ where
}\alpha^{-1}\coloneqq(\kappa-1)\sqrt{2\sqrt{\kappa}(\kappa-1)(\sqrt{\kappa}-3)}.$
(11)

Let us first define the coefficients $\tau:=\frac{1}{\sqrt{\kappa}+1}$,
$\gamma:=\frac{1}{\sqrt{\kappa}-1}$, and
$\rho=\frac{\mu\gamma}{4(4+\gamma)}=\frac{\mu}{4(4\sqrt{\kappa}-3)}$. $\tau$
can be understood as a momentum parameter, $\gamma$ as (proportional to) the
parameter of the exponential curve that is our convergence rate, and $\rho$ a
penalty for the presence of the perturbations $\Delta_{t}$. As is done in all
proofs of Nesterov’s method, we will define the iterates

$z_{t}:=\tfrac{1}{\tau}x_{t}-\tfrac{1-\tau}{\tau}y_{t}$

for all $t\in[T]$; we refer interested readers to Zhu & Orecchia (2017); Ahn &
Sra (2022) for interpretations of these iterates. We note that although our
definition of $z_{T}$ depends on $x_{T}$, and by extension $\Delta_{T}$, this
lemma’s claim about $y_{T}$ does not directly depend at all on $\Delta$. Thus,
we will freely choose to set $\Delta_{T}=0$, guaranteeing that
$\left\lVert\Delta_{T}\right\rVert\leq\alpha\left\lVert
y_{T}-y_{T-1}\right\rVert$.

To prove the statement of this lemma, it suffices to show the following claim.

######  Claim 3.6.

The following function $\Phi$ is a potential function, i.e. for all $t\in[T]$,
$\Phi(t)-\Phi(t-1)\leq 0$:

$\Phi(t)=(1+\tfrac{\gamma}{2})^{t}[(\ell(y_{t})-\ell(x^{*}))+\tfrac{\mu}{2}\|z_{t}-x^{*}\|^{2}+2\rho\|y_{t}-x^{*}\|^{2}].$
(12)

It suffices to show Claim 3.6 because the consequence that
$\Phi(T)\leq\Phi(0)$ directly implies the lemma’s statement via

$\displaystyle(1+\tfrac{\gamma}{2})^{T}\ell(y_{T})-\ell(x^{*})$
$\displaystyle\leq\ell(x_{0})-\ell(x^{*})+(\frac{\mu}{2}+2\rho)\left\lVert
x_{0}-x^{*}\right\rVert^{2}\leq(\tfrac{L{}+\mu}{2}+2\rho)\left\lVert
x_{0}-x^{*}\right\rVert^{2},$

with the last inequality following from the $L{}$-smoothness of $\ell$. We
therefore turn to showing Claim 3.6.

First, we note that the potential function (12) we are proving differs from
the usual potential function that is used to prove Nesterov’s acceleration,
which we will denote by $\Phi_{\mathrm{orig}}$:

$\displaystyle\Phi_{\mathrm{orig}}(t)=(1+\gamma)^{t}[(\ell(y_{t})-\ell(x^{*}))+\tfrac{\mu}{2}\left\lVert
z_{t}-x^{*}\right\rVert^{2}].$

Fact 3.7 says roughly that, for any $t\in[T]$, we can hypothetically recover a
large part of the usual proof of Nesterov’s acceleration by bounding the
difference $\Phi_{\mathrm{orig}}(t)-\Phi_{\mathrm{orig}}(t-1)\leq 0$, if only
we could remove the perturbation at timestep $t$, i.e. set $\Delta_{t}=0$ (but
keep in place the perturbations $\Delta_{1},\dots,\Delta_{t-1}$ from previous
timesteps).

######  Fact 3.7.

Fix any $t\in[T]$ and let
$\widetilde{z}_{t}\coloneqq\tfrac{1}{\tau}(y_{t}+\beta(y_{t}-y_{t-1}))-\tfrac{1-\tau}{\tau}y_{t}$.
Then,

$(1+\gamma)(\ell(y_{t})-\ell(x^{*}))-(\ell(y_{t-1})-\ell(x^{*}))+\tfrac{\mu}{2}\left((1+\gamma)\|\widetilde{z}_{t}-x^{*}\|^{2}-\|z_{t-1}-x^{*}\|^{2}\right)\leq
0$ (13)

Next, in Fact 3.8, we show that, since we can argue that
$\Phi_{\mathrm{orig}}(t)\leq\Phi_{\mathrm{orig}}(t-1)$ if we ignore the
perturbation $\Delta_{t}$, we can argue that $\Phi(t)\leq\Phi(t-1)$ even
taking the perturbation $\Delta_{t}$ into account. That is, Fact 3.8 shows
that the left-hand side of (13) upper bounds $\Phi(t)-\Phi(t-1)$, so long as
our stated assumption on the perturbation norm
$\left\lVert\Delta_{t}\right\rVert$ holds.

######  Fact 3.8.

For any $t\in[T]$, the left-hand side of (13) upper bounds
$\Phi(t)-\Phi(t-1)$; equivalently,

$\displaystyle(1+\tfrac{\gamma}{2})(\tfrac{\mu}{2}\left\lVert
z_{t}-x^{*}\right\rVert^{2}+2\rho\left\lVert
y_{t}-x^{*}\right\rVert^{2})-2\rho\left\lVert
y_{t-1}-x^{*}\right\rVert^{2}\leq\tfrac{\gamma}{2}(\ell(y_{t})-\ell(x^{*}))+(1+\gamma)\tfrac{\mu}{2}\left\lVert\widetilde{z}_{t}-x^{*}\right\rVert^{2}.$
(14)

Plugging Fact 3.8’s 14 into Fact 3.7’s 13, we recover our main claim that
$\Phi(t)-\Phi(t-1)\leq 0$. We conclude by turning to prove Facts 3.7 and 3.8.

###### Proof of Fact 3.7.

The proof of this fact closely follows Bansal & Gupta (2017)’s proof of
Nesterov’s accelerated gradient method in smooth strongly convex settings.

We begin by upper bounding the first summand in the left-hand side of (13).
Since $\ell$ is smooth and $y_{t}$ is a gradient step on $\ell$ from
$x_{t-1}$, we have that
$\ell(y_{t})\leq\ell(x_{t-1})-\frac{1}{2L{}}\left\lVert\nabla
x_{t-1}\right\rVert^{2}$ and thus

$\displaystyle(1+\gamma)(\ell(y_{t})-\ell(x^{*}))-(\ell(y_{t-1})-\ell(x^{*}))\leq\ell(x_{t-1})-\ell(y_{t-1})+\gamma(\ell(x_{t-1})-\ell(x^{*}))-\frac{1+\gamma}{2L{}}\left\lVert\nabla\ell(x_{t-1})\right\rVert^{2}.$
(15)

Using the $\mu$-strong convexity of $\ell$, we can further bound part of the
right-hand side by

$\displaystyle\ell(x_{t-1})-\ell(y_{t-1})+\gamma(\ell(x_{t-1})-\ell(x^{*}))\leq\;$
$\displaystyle\left<\nabla\ell(x_{t-1}),x_{t-1}-y_{t-1}\right>$
$\displaystyle\;+\gamma\left(\left<\nabla\ell(x_{t-1}),x_{t-1}-x^{*}\right>-\frac{\mu}{2}\left\lVert
x_{t-1}-x^{*}\right\rVert^{2}\right).$ (16)

Plugging (16) and the identity
$z_{t-1}=\frac{1}{\tau}x_{t-1}-\frac{1-\tau}{\tau}y_{t-1}$ into (15), direct
algebra yields an upper bound on the first summand of (13):

$\displaystyle(1+\gamma)(\ell(y_{t})-\ell(x^{*}))-(\ell(y_{t-1})-\ell(x^{*}))$
$\displaystyle\leq\frac{1}{1+\gamma}\left<\nabla\ell(x_{t-1}),\gamma(z_{t-1}-x^{*})+\gamma^{2}(x_{t-1}-x^{*})\right>-\frac{\mu\gamma}{2}\left\lVert
x_{t-1}-x^{*}\right\rVert^{2}-\frac{1+\gamma}{2L{}}\left\lVert\nabla\ell(x_{t-1})\right\rVert^{2}.$
(17)

Next, we turn to upper bounding the second summand in (13). Plugging the
iterate definitions $y_{t}=x_{t-1}-\frac{1}{L{}}\nabla\ell(x_{t-1})$ and
${z}_{t-1}=\frac{1}{\tau}{x}_{t-1}-\frac{1-\tau}{\tau}{y}_{t-1}$ into our
definition of
$\widetilde{z}_{t}=\frac{1}{\tau}\left((1-\tau)y_{t}-(1-2\tau)y_{t-1}\right)$,
we can recover the identity

$\widetilde{z}_{t}=\frac{1}{1+\gamma}z_{t-1}+\frac{\gamma}{1+\gamma}x_{t-1}-\frac{\gamma}{\mu(1+\gamma)}\nabla\ell(x_{t-1}).$

Plugging this identity for $\widetilde{z}_{t}$ into the expression
$\tfrac{\mu}{2}\left((1+\gamma)\|\widetilde{z}_{t}-x^{*}\|^{2}-\|z_{t-1}-x^{*}\|^{2}\right)$,
direct algebra yields the identity

$\displaystyle\tfrac{\mu}{2}\left((1+\gamma)\|\widetilde{z}_{t}-x^{*}\|^{2}-\|z_{t-1}-x^{*}\|^{2}\right)=\;$
$\displaystyle\frac{1+\gamma}{2L{}}\left\lVert\nabla\ell(x_{t-1})\right\rVert^{2}-\frac{1}{1+\gamma}\left<\nabla\ell(x_{t-1}),\gamma(z_{t-1}-x^{*})+\gamma^{2}(x_{t-1}-x^{*})\right>$
$\displaystyle\;+\frac{\mu\gamma}{2}\left\lVert
x_{t-1}-x^{*}\right\rVert^{2}-\frac{\mu\gamma}{2(1+\gamma)}\left\lVert
z_{t-1}-x_{t-1}\right\rVert^{2}.$ (18)

Summing (17) and (18) yields the following upper bound on the left-hand side
of (13):

$\displaystyle(1+\gamma)(\ell(y_{t})-\ell(x^{*}))-(\ell(y_{t-1})-\ell(x^{*}))+\tfrac{\mu}{2}\left((1+\gamma)\|\widetilde{z}_{t}-x^{*}\|^{2}-\|z_{t-1}-x^{*}\|^{2}\right)$
$\displaystyle\leq-\frac{\mu\gamma}{2(1+\gamma)}\left\lVert
z_{t-1}-x_{t-1}\right\rVert^{2}\leq 0.$

∎

###### Proof of Fact 3.8.

To show this fact, we can observe from direct algebra that it suffices to
prove the following inequalities:

$\displaystyle\tfrac{\gamma}{2}(\ell(y_{t})-\ell(x^{*}))\geq(2+\tfrac{\gamma}{2})2\rho\left\lVert
y_{t}-x^{*}\right\rVert^{2},$ (19)
$\displaystyle(1+\gamma)\left\lVert\widetilde{z}_{t}-x^{*}\right\rVert^{2}\geq(1+\tfrac{\gamma}{2})\left\lVert
z_{t}-x^{*}\right\rVert^{2}-\tfrac{4\rho}{\mu}(\left\lVert
y_{t-1}-x^{*}\right\rVert^{2}+\left\lVert y_{t}-x^{*}\right\rVert^{2}).$ (20)

The inequality (19) follows directly from the $\mu$-convexity of $\ell$, and
the identity $\rho=\frac{\mu\gamma}{4(4+\gamma)}$, as

$\tfrac{\gamma}{2}(\ell(y_{t})-\ell(x^{*}))\geq\tfrac{\mu\gamma}{4}\|y_{t}-x^{*}\|^{2}=(2+\tfrac{\gamma}{2})\cdot
2\rho\|y_{t}-x^{*}\|^{2}.$

To show (20), we first recall the general fact that $c(1+\delta)\|v\|^{2}\geq
c\|v+w\|^{2}-c(1+\frac{1}{\delta})\|w\|^{2}$ for any $\delta,c>0$. Choosing
$v=\widetilde{z}_{t}-x^{*}$, $w=z_{t}-\widetilde{z}_{t}$,
$c=1+\tfrac{\gamma}{2}$, and $\delta=\frac{\gamma}{\gamma+2}>0$, this yields

$\displaystyle
1+\gamma\|\widetilde{z}_{t}-x^{*}\|^{2}\geq(1+\tfrac{\gamma}{2})\|z_{t}-x^{*}\|^{2}-(2+\tfrac{2}{\gamma})(1+\tfrac{\gamma}{2})\|\tfrac{1}{\tau}\Delta_{t}\|^{2}.$
(21)

Using the perturbation norm bound, we can therefore bound

$(2+\tfrac{2}{\gamma})(1+\tfrac{\gamma}{2})\|\tfrac{1}{\tau}\Delta_{t}\|^{2}\leq\tfrac{2\rho}{\mu}\|y_{t}-y_{t-1}\|^{2}\leq\tfrac{4\rho}{\mu}(\|y_{t}-x^{*}\|^{2}+\|y_{t-1}-x^{*}\|^{2}),$

where the last inequality uses the triangle inequality; plugged into (21),
this yields (20) as desired. ∎

∎

Lemma 3.5 means that even though stagewise training with stacking differs from
Nesterov’s method, their similar form allows us to express the former as a
perturbation of the latter. In particular, if we write our stacking update (6)
for deep residual linear networks as a perturbation of Nesterov’s method, as
in (10), the perturbation term is exactly

$\displaystyle\Delta_{t}=\beta(W_{t}-W_{t-1})W_{t-1}^{-1}W_{t}-\beta(W_{t}-W_{t-1})$
(22)

for all $t\in[T-1]$ and $\Delta_{0}=0$. That is, if we examine the sequence of
iterates $W_{1},\dots,W_{T}$ produced by stagewise training with stacking
initialization, the term $\Delta_{t}$ is a measure of the disagreement between
the realized iterate $W_{t+1}$ and what Nesterov’s method says the iterate
$W_{t+1}$ should be conditioned on the iterates $W_{1},\dots,W_{t}$ from
previous timesteps. Note that, even when these perturbation terms
$\Delta_{1},\dots,\Delta_{T-1}$ are small in norm, it is possible for
Nesterov’s method to describe an iterate sequence that diverges significantly
in norm from the iterates realized by stagewise training.

Rewriting (22) as
$\Delta_{t}=\beta(W_{t}-W_{t-1})W_{t-1}^{-1}(W_{t}-W_{t-1})$, we immediately
see that to satisfy the requirement of Lemma 3.5 that
$\left\lVert\Delta_{t}\right\rVert\in O(\kappa^{-2}\left\lVert
W_{t}-W_{t-1}\right\rVert)$, we simply need $\left\lVert\beta
W_{t-1}^{-1}(W_{t}-W_{t-1})\right\rVert\in O(\kappa^{-2})$. This is satisfied
when $W_{t}\approx W_{t-1}$ and $W_{t-1}$ is reasonably non-singular. A
sufficient condition for this is that the iterates are sufficiently close to
the ground-truth solution $W^{*}$ and $W^{*}$ is non-singular, which explains
the conditions of Theorem 3.1.

##  4 Experiments

In this section we empirically validate our theoretical results.

###  4.1 Deep Linear Networks and Squared Losses

As our main theoretical results in Section 3 apply to the case of deep linear
networks, we consider the same function class in our experiments on synthetic
data with the square loss. Formally, the output space
$\mathcal{Y}=\mathbb{R}^{d}$, and for a predictor $x\mapsto W_{t}x$, we
consider the loss

$\displaystyle\ell(W)=\mathbb{E}_{(x,y)\sim
D}\left[\tfrac{1}{2}\|Wx-y\|^{2}\right].$ (23)

We consider a data distribution $D$ where the samples $(x,y)$ are drawn as
follows. Let $W^{*}$ and be the “ground truth” positive definite matrix and
let $\Sigma$ be the data covariance matrix. We first sample $x\sim
N(0,\Sigma)$ and then conditioned on $x$, the output $y$ is generated as
$y=W^{*}x+\xi$, where $\xi$ is a mean zero random variable.

We can then write the expected square loss explicitly as

$\displaystyle\ell(W)$ $\displaystyle=\mathbb{E}_{(x,y)\sim
D}\left[\tfrac{1}{2}\|Wx-y\|^{2}\right]$
$\displaystyle=\tfrac{1}{2}\text{Tr}((W-W^{*})\Sigma(W-W^{*})^{\top})+\tfrac{1}{2}\mathbb{E}\left[\|\xi\|^{2}\right].$
(24)

Note that for the case of the squared loss described above, the condition
number of the expected loss depends on the covariance matrix $\Sigma$, i.e.,
$\kappa=\frac{\sigma_{\text{max}}(\Sigma)}{\sigma_{\text{min}}(\Sigma)}$.

#### Stacking updates.

For the specific case of the squared loss we get the following closed form
expression of the stacking updates (see Eq. (6))

$\displaystyle W_{t+1}=W_{t}+\beta w_{t}W_{t}+\nabla\ell(W_{t}+\beta
w_{t}W_{t})$
$\displaystyle=(W_{t}+\beta(W_{t}-W_{t-1})W_{t-1}^{-1}W_{t})(I-\tfrac{1}{L}{}\Sigma)+\tfrac{1}{L}{}W^{*}\Sigma.$
(25)

Here $L$ is the smoothness of the loss which in turn depends on the largest
singular value of $\Sigma$, i.e. $L=\sigma_{\text{max}}(\Sigma)$.

#### Nesterov’s updates.

Similarly, we get the following closed form expression for obtaining
Nesterov’s updates (see Eq. (7)) for the case of squared loss

$\displaystyle W_{t+1}=W_{t}+\beta w_{t}W_{t}+\nabla\ell(W_{t}+\beta
w_{t}W_{t})$
$\displaystyle=(W_{t}+\beta(W_{t}-W_{t-1}))(I-\tfrac{1}{L}{}\Sigma)+\tfrac{1}{L}{}W^{*}\Sigma.$
(26)

For both stacking and Nesterov’s updates we set
$\beta=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$.

###  4.2 Synthetic data experiments

We compare the performance of the three updates namely vanilla gradient
descent, stacking updates (Eq. 25) and exact Nesterov updates (Eq. 26). Here
at each stacking stage only the last layer is updated which matches our
theoretical setup faithfully. Later we will also consider the effect of
training all the layers in each stacking stage which is closer to the manner
in which stacking is applied in practice.

We consider points in $d=20$ dimensions. We generate the ground truth $W^{*}$
to be of the form $I+\sigma S$, where $S$ is a random positive semi-definite
matrix of spectral norm $1$ and $\sigma$ is a parameter defining the closeness
of $W^{*}$ to identity. For a given value $\kappa>1$, we generate a random
covariance matrix $\Sigma$. Finally, we sample the noise $\xi$ in the output
from a mean zero Gaussian with a standard deviation of $0.1$.

Figure 4: Mean squared error (MSE) vs. number of stacking stages. We observe
that as the data becomes more ill conditioned both the stacking updates and
Nesterov’s updates demonstrate faster convergence than vanilla gradient
descent.

Figure 4 shows the performance of the three types of updates as the problem
becomes more ill conditioned, i.e. as a function of $\kappa$. As expected, at
small values of the condition number there is no advantage of the stacking
updates over vanilla gradient descent. However, for ill conditioned data the
stacking updates converge much faster than gradient descent. We also observe
that the convergence of the stacking updates mirrors very closely the
convergence behavior of the exact Nesterov’s updates.

To further understand the relationship between the stacking updates and
Nesterov’s updates, in Figure 5 we show the performance of the two as the
distance of $W^{*}$ from identity increases. As can be seen from the figures,
when $W^{*}$ is farther from identity the stacking updates behave
qualitatively different from Nesterov’s updates during the initial phase where
the loss for stacking updates explodes before converging in later stages. This
suggests that in practice there may be a better way to initialize a stacking
stage by making the initialization closer to the ideal Nesterov’s updates.
While in the case of deep linear networks and the squared loss we have a
closed form expression for such an initialization, in general this is a hard
problem.

Next we consider the case where in each stacking stage we train all the layers
of the deep linear network. We use the same data generation procedure as
described above. We perform $10$ stages of stacking where in each stage we
perform $2$ steps of gradient descent with a learning rate of $1/L$ where $L$
is the smoothness of the loss function. We train on $1024$ examples with batch
size of $32$ and test on $1024$ examples.

We consider two types of stacking based initialization schemes. The first one
namely Stacking Init. initializes the next layer’s weight matrix $w^{0}_{t+1}$
as $\beta w_{t}$. The second scheme namely Nesterov Init. initializes
$w^{0}_{t}$ such that we recover the precises Nesterov’s updates at
initialization, i.e., Eq. 26. From the analysis in Section 2 the
initialization that achieves this amounts to setting $w^{0}_{t+1}$ as $\beta
w_{t}(I+w_{t})^{-1}$.

Figure 5: Mean squared error (MSE) vs. number of stacking stages. The figure
shows the comparison of stacking updates and Nesterov’s updates as $W^{*}$
becomes farther from Identity, i.e. $\sigma$ increases. We observe that for
higher values of $\sigma$ the stacking updates display a diverging behavior in
the initial stages.

Figure 6 shows the performance of the two stacking initialization schemes as
compared to the random baseline where we initialize the next layer’s weight
matrix to be a random one. We again observe that both the stacking schemes
outperform the baseline particularly when the data is ill conditioned.

Figure 6: Mean squared error (MSE) vs. number of stacking stages when training
all the layers.

###  4.3 Stacking for BERT Base with $\beta$ parameters

The theory developed in Section 2 requires the initialization at the
$(t+1)$-th stage to be $f_{t+1}^{0}=\beta f_{t}$ for some $\beta\in[0,1)$. The
introduction of $\beta$ is crucial to get the accelerated convergence rate in
Nesterov’s method, but the standard stacking initialization doesn’t use a
$\beta$ parameter. We performed sanity check experiments on BERT Base to
ensure that the introduction of the $\beta$ parameter doesn’t affect the
efficacy of stacking. We introduced a trainable parameter, $\beta$, that
multiplies the output of the newly added transformer block in stacking, which
is initialized to the values $0.9$ and $0.99$, which are standard settings for
momentum parameters. Figure 7 shows that introduction of the $\beta$ parameter
doesn’t hurt the efficacy of stacking. The plot also shows that the final log
perplexity improves a bit when using trainable $\beta$.

Figure 7: Stacking initialization with trainable $\beta$ parameter multiplying
the output of the newly added transformer block. Experimental runs with
$\beta$ initialized to $0.99$ and $0.9$ are provided.

##  5 Conclusions and Future Work

This paper develops the perspective that the effectiveness of stacking
initialization, compared to other forms of initialization such as zero or
random, is because it enables a form of accelerated gradient descent in
function space. There are several directions for future work. While this work
provides a formal proof of accelerated convergence for a particular parametric
setting (deep residual linear networks), such a proof in the general
functional setting for deep residual networks is still open, and will probably
require some additional assumptions. From a practical standpoint, a very
intriguing and potentially impactful question is whether it is possible to
come up with an efficiently implementable initialization scheme that leads to
Nesterov’s AGD updates exactly for deep residual networks.

## References

  * Achiam et al. (2023) Achiam, J., Adler, S., Agarwal, S., Ahmad, L., Akkaya, I., Aleman, F. L., Almeida, D., Altenschmidt, J., Altman, S., Anadkat, S., et al.  Gpt-4 technical report.  _arXiv preprint arXiv:2303.08774_ , 2023.
  * Ahn & Sra (2022) Ahn, K. and Sra, S.  Understanding nesterov’s acceleration via proximal point method.  In Bringmann, K. and Chan, T. (eds.), _5th Symposium on Simplicity in Algorithms, SOSA@SODA 2022, Virtual Conference, January 10-11, 2022_ , pp. 117–130. SIAM, 2022.
  * Bansal & Gupta (2017) Bansal, N. and Gupta, A.  Potential-function proofs for first-order methods, 2017.
  * Beltagy et al. (2020) Beltagy, I., Peters, M. E., and Cohan, A.  Longformer: The long-document transformer.  _arXiv preprint arXiv:2004.05150_ , 2020.
  * Bengio et al. (2006) Bengio, Y., Lamblin, P., Popovici, D., and Larochelle, H.  Greedy layer-wise training of deep networks.  _Advances in neural information processing systems_ , 19, 2006.
  * Brown et al. (2020) Brown, T., Mann, B., Ryder, N., Subbiah, M., Kaplan, J. D., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., et al.  Language models are few-shot learners.  _Advances in neural information processing systems_ , 33:1877–1901, 2020.
  * Chen & Guestrin (2016) Chen, T. and Guestrin, C.  Xgboost: A scalable tree boosting system.  _CoRR_ , abs/1603.02754, 2016.  URL http://arxiv.org/abs/1603.02754.
  * Chen et al. (2022) Chen, X., Wang, X., Changpinyo, S., Piergiovanni, A., Padlewski, P., Salz, D., Goodman, S., Grycner, A., Mustafa, B., Beyer, L., et al.  Pali: A jointly-scaled multilingual language-image model.  _arXiv preprint arXiv:2209.06794_ , 2022.
  * Chowdhery et al. (2022) Chowdhery, A., Narang, S., Devlin, J., Bosma, M., Mishra, G., Roberts, A., Barham, P., Chung, H. W., Sutton, C., Gehrmann, S., et al.  Palm: Scaling language modeling with pathways.  _arXiv preprint arXiv:2204.02311_ , 2022.
  * Devlin et al. (2019) Devlin, J., Chang, M., Lee, K., and Toutanova, K.  BERT: pre-training of deep bidirectional transformers for language understanding.  In Burstein, J., Doran, C., and Solorio, T. (eds.), _NAACL-HLT_ , pp. 4171–4186. Association for Computational Linguistics, 2019.  doi: 10.18653/V1/N19-1423.  URL https://doi.org/10.18653/v1/n19-1423.
  * Freund & Schapire (1997) Freund, Y. and Schapire, R. E.  A decision-theoretic generalization of on-line learning and an application to boosting.  _Journal of computer and system sciences_ , 55(1):119–139, 1997.
  * Friedman (2001) Friedman, J. H.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ , pp. 1189–1232, 2001.
  * Gong et al. (2019) Gong, L., He, D., Li, Z., Qin, T., Wang, L., and Liu, T.  Efficient training of bert by progressively stacking.  In _International conference on machine learning_ , pp. 2337–2346. PMLR, 2019.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., and Courville, A.  _Deep Learning_.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * Gu et al. (2020) Gu, X., Liu, L., Yu, H., Li, J., Chen, C., and Han, J.  On the transformer growth for progressive bert training.  _arXiv preprint arXiv:2010.12562_ , 2020.
  * Gupta et al. (2018) Gupta, V., Koren, T., and Singer, Y.  Shampoo: Preconditioned stochastic tensor optimization.  In Dy, J. G. and Krause, A. (eds.), _ICML_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 1837–1845. PMLR, 2018\.  URL http://proceedings.mlr.press/v80/gupta18a.html.
  * Hardt & Ma (2017) Hardt, M. and Ma, T.  Identity matters in deep learning.  In _ICLR_. OpenReview.net, 2017.  URL https://openreview.net/forum?id=ryxB0Rtxx.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _CVPR_ , pp. 770–778. IEEE Computer Society, 2016.  doi: 10.1109/CVPR.2016.90.  URL https://doi.org/10.1109/CVPR.2016.90.
  * Hinton et al. (2006) Hinton, G. E., Osindero, S., and Teh, Y.-W.  A fast learning algorithm for deep belief nets.  _Neural computation_ , 18(7):1527–1554, 2006\.
  * Ioffe & Szegedy (2015) Ioffe, S. and Szegedy, C.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In Bach, F. R. and Blei, D. M. (eds.), _ICML_ , volume 37 of _JMLR Workshop and Conference Proceedings_ , pp. 448–456. JMLR.org, 2015\.  URL http://proceedings.mlr.press/v37/ioffe15.html.
  * Jiang et al. (2024) Jiang, A. Q., Sablayrolles, A., Roux, A., Mensch, A., Savary, B., Bamford, C., Chaplot, D. S., Casas, D. d. l., Hanna, E. B., Bressand, F., et al.  Mixtral of experts.  _arXiv preprint arXiv:2401.04088_ , 2024.
  * Kawaguchi (2016) Kawaguchi, K.  Deep learning without poor local minima.  In Lee, D. D., Sugiyama, M., von Luxburg, U., Guyon, I., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 29_ , pp. 586–594, 2016.
  * Liu et al. (2023) Liu, H., Li, Z., Hall, D., Liang, P., and Ma, T.  Sophia: A scalable stochastic second-order optimizer for language model pre-training.  _arXiv preprint arXiv:2305.14342_ , 2023.
  * Mason et al. (1999) Mason, L., Baxter, J., Bartlett, P. L., and Frean, M. R.  Boosting algorithms as gradient descent.  In Solla, S. A., Leen, T. K., and Müller, K. (eds.), _NeurIPS_ , pp. 512–518. The MIT Press, 1999.  URL http://papers.nips.cc/paper/1766-boosting-algorithms-as-gradient-descent.
  * Nesterov (1983) Nesterov, Y.  A method for unconstrained convex minimization problem with the rate of convergence o (1/k^ 2).  In _Doklady an ussr_ , volume 269, pp. 543–547, 1983.
  * Radford et al. (2019) Radford, A., Wu, J., Child, R., Luan, D., Amodei, D., Sutskever, I., et al.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9, 2019.
  * Reddi et al. (2023) Reddi, S. J., Miryoosefi, S., Karp, S., Krishnan, S., Kale, S., Kim, S., and Kumar, S.  Efficient training of language models using few-shot learning.  In _International Conference on Machine Learning_ , pp. 14553–14568. PMLR, 2023.
  * Saxe et al. (2014) Saxe, A. M., McClelland, J. L., and Ganguli, S.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  In Bengio, Y. and LeCun, Y. (eds.), _2nd International Conference on Learning Representations, ICLR 2014, Conference Track Proceedings_ , 2014.
  * Schapire & Freund (2012) Schapire, R. E. and Freund, Y.  _Boosting: Foundations and Algorithms_.  The MIT Press, 2012.  ISBN 0262017180.
  * Shazeer & Stern (2018) Shazeer, N. and Stern, M.  Adafactor: Adaptive learning rates with sublinear memory cost.  In _International Conference on Machine Learning_ , pp. 4596–4604. PMLR, 2018.
  * Shazeer et al. (2017) Shazeer, N., Mirhoseini, A., Maziarz, K., Davis, A., Le, Q., Hinton, G., and Dean, J.  Outrageously large neural networks: The sparsely-gated mixture-of-experts layer.  _arXiv preprint arXiv:1701.06538_ , 2017.
  * Shen et al. (2022) Shen, S., Walsh, P., Keutzer, K., Dodge, J., Peters, M., and Beltagy, I.  Staged training for transformer language models.  In _International Conference on Machine Learning_ , pp. 19893–19908. PMLR, 2022.
  * Team et al. (2023) Team, G., Anil, R., Borgeaud, S., Wu, Y., Alayrac, J.-B., Yu, J., Soricut, R., Schalkwyk, J., Dai, A. M., Hauth, A., et al.  Gemini: a family of highly capable multimodal models.  _arXiv preprint arXiv:2312.11805_ , 2023.
  * Tu et al. (2023) Tu, T., Azizi, S., Driess, D., Schaekermann, M., Amin, M., Chang, P.-C., Carroll, A., Lau, C., Tanno, R., Ktena, I., et al.  Towards generalist biomedical AI.  _arXiv preprint arXiv:2307.14334_ , 2023.
  * Vaswani et al. (2017) Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Zhu & Orecchia (2017) Zhu, Z. A. and Orecchia, L.  Linear coupling: An ultimate unification of gradient and mirror descent.  In Papadimitriou, C. H. (ed.), _8th Innovations in Theoretical Computer Science Conference, ITCS 2017, January 9-11, 2017, Berkeley, CA, USA_ , volume 67 of _LIPIcs_ , pp. 3:1–3:22. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2017.


# Short-Term Load Forecasting Using Time Pooling Deep Recurrent Neural Network

Elahe Khoshbakhti Vaygan  Faculty of ECE
Qom University of Technology
Qom, Iran
<EMAIL_ADDRESS>Roozbeh Rajabi  Faculty of ECE
Qom University of Technology
Qom, Iran
<EMAIL_ADDRESS>Abouzar Estebsari  School of the Built Environment and
Architecture
London South Bank University
London, United Kingdom
<EMAIL_ADDRESS>
###### Abstract

Integration of renewable energy sources and emerging loads like electric
vehicles to smart grids brings more uncertainty to the distribution system
management. Demand Side Management (DSM) is one of the approaches to reduce
the uncertainty. Some applications like Nonintrusive Load Monitoring (NILM)
can support DSM, however they require accurate forecasting on high resolution
data. This is challenging when it comes to single loads like one residential
household due to its high volatility. In this paper, we review some of the
existing Deep Learning-based methods and present our solution using Time
Pooling Deep Recurrent Neural Network. The proposed method augments data using
time pooling strategy and can overcome overfitting problems and model
uncertainties of data more efficiently. Simulation and implementation results
show that our method outperforms the existing algorithms in terms of RMSE and
MAE metrics.

###### Index Terms:

Load forecasting, Time pooling, Deep learning, Smart grids.

##  I Introduction

One of the main objectives of distribution systems is to ensure continuity of
service to final customers and make a balance between supply and demand.
Therefore, it is crucial to estimate both generation and load in advance to
manage system in operation more efficiently [1]. Load forecasting has been a
key element in Distribution Management Systems (DMS) in operation and planning
analyses as storing electric energy in large amount is very costly. An
accurate estimation of load help match supply to avoid imbalances due to power
deficit or surplus.

The measurements of electrical load over a long period of time not only show
different values of loads over a periodical cycle (e.g. 24 hours), but also
realize some trends of load variation in midterm or long-term horizons [2],
[3]. Therefore, the load forecasting methods aim at predicting the load values
in different time horizons using historic measurements, load trends analysis,
experimental rules and mathematical models [4]. Some early load forecasting
methods included Experimental Smoothing Models, Moving Average, and Auto
Regression. The forecasting methods could be classified in three groups based
on the techniques used: Gray Prediction Models, Statistical Analysis Models
and Non-linear Intelligent Models [5], [6]. The most common practical load
forecasting methods in traditional passive distribution networks are Auto
Regressive Integration Moving Average (ARIMA) and Support Vector Regression
(SVR) [5] , [7]. Integration of more electronic intelligent devices and
communication technologies to distribution systems from one side and changing
supply paradigms from bulk generation based models to distributed generation
create o-called Smart Grids. In Smart Grids, there is a great opportunity for
many advanced operation and analysis methods to be applied and utilized. These
methods include Artificial Intelligence (AI) application like Machine Learning
(ML) techniques and Deep Learning (DL) methods [8], [9], [10], [11].

In This paper, an AI method based on deep learning and pooling strategy
techniques is proposed to increase the accuracy of load forecasting of the
most volatile loads including single residential households. Similar to
supervised multi-layer neural networks, the historic data is divided into two
groups containing test and training sets, pooling test and pooling train,
respectively. Before presenting our model, the complexity of load forecast in
single residential customers and the parameters involved is briefed, and a
literature review of ARIMA, SVR, Recurrent Neural Network (RNN) and Deep RNN
is presented. The main contribution of this paper is to use time pooling
strategy to enhance the performance of deep recurrent neural networks,
overcoming the overfitting problem, and better modeling of the data.

The rest of the paper is organized as follows. Section II reviews the
challenges regarding load forecasting and its applications. In this section
previous methods used for short-term load forecasting also summarized. The
proposed method using pooling strategy and deep learning is presented in III.
Section IV describes the used dataset and experiments to evaluated the
proposed TPRNN. Finally section V concludes the paper.

##  II Literature Review

The residential appliances cover a wide range of devices with different
characteristics, duty cycles, switching requirements, and power consumption.
They are also different in terms of load types, from resistive types to
inductive loads. This would make the aggregated profile very volatile with
sharp increases or decreases of power in short periods of time. Moreover, the
consumption behavior of residential customers is not so predictable and there
are always a high level of uncertainty. The uncertainty is due to different
cultural and social consumption patterns, weather conditions which could be
intermittent, and residual loads considered as noise. Since residential loads
have a substantial share in total demand of distribution systems, the
prediction of such stochastic loads is essential to manage the system. The
importance of more accurate load forecasting is stressed when the system faces
intermittent behavior of renewable energy sources (RES) widely integrated into
the network.

Due to high volatility and uncertainty of single residential loads, we believe
advanced ML-based methods should be applied to get a high accuracy in the
results. DL techniques have been proposed in literature and in several
projects as efficient methods, however the accuracy and performance of the
applying DLs directly are arguable. In these methods, increase of deep layers
is a fundamental approach to improve the accuracy of the results compare to
Artificial Neural Networks (ANN), but in single residential load forecasting,
the increase of DL layers exponentially increases the intrinsic parameters of
the network in the training set. This would eventually grow the system noise
which can affect the forecasting results. In summary, it is important to
carefully model and adapt the DL network for for particular types of loads and
specific forecast time horizons in order to achieve desired accuracy and
performance [12], [13]. Some parameters which may be used as extra input to
the model include the economical and environmental conditions of the loads
such as temperature, humidity, lighting, population growth, wind speed,
economical situation, time and policies (Fig. 1).

Figure 1: Important factors having impact on load consumption.

In the following we briefly review 4 general methods to forecast loads in
short-term: Auto Regressive Integration Moving Average (ARIMA), Support Vector
Regression (SVR), Recurrent Neural Network (RNN) and Deep Recurrent Neural
Network (DRNN).

### Auto Regressive Integration Moving Average

Time series models are the most common and popular methods to perform short-
term forecasting. The Auto Regressive Integration Moving Average model is
derived from Regressive Moving Average method which is time-series based
method. In this model, time-series $\mathrm{y}(\mathrm{t})$ is defined based
on the previous values $y(t-1),y(t-2),\ldots,y(t-p)$ and stochastic noise. The
order of the process depends on the oldest previous value to which
$\mathrm{y}(\mathrm{t})$ is returned. This model is formulated as follow:

$\displaystyle\mathrm{y}(\mathrm{t})=\delta+\phi_{1}\mathrm{y}(\mathrm{t}+1)+\cdots+$
(1)
$\displaystyle\phi_{\mathrm{p}}\mathrm{y}(\mathrm{t}-\mathrm{p})+\varepsilon(\mathrm{t})$

where $y(t-1),y(t-2),\ldots,y(t-p)$ are the previous values of time series,
$\varepsilon(t)$ is the stochastic noise, and $\delta$ is defined as follow
[5]:

$\delta\equiv\left(1-\sum_{i=1}^{p}\varphi_{i}\right)\mu$ (2)

where $\mu$ is the average value of the process.

### Support Vector Regression

SVR is based on Vapnik-chervonenkis theory and aim at expanding data sets to
those not seen yet [7], [14]. SVM can be extended to be used as SVR. For this,
a zone named $\varepsilon$ -insensitivities defined around the function as
$\varepsilon$ -tube; then the SVR is formulated as follow [7]:

$\begin{gathered}\mathrm{f}(\mathrm{x})=\left(\begin{array}[]{c}\mathrm{x}\\\
1\end{array}\right)\left(\begin{array}[]{l}\mathrm{W}\\\
\mathrm{b}\end{array}\right)^{\mathrm{T}}=w^{T}+\mathrm{b}\\\
\mathrm{x},\mathrm{w}\in\mathrm{R}^{\mathrm{m}+1}\end{gathered}$ (3)

### Recurrent Neural Network

The Recurrent Neural Network is based on Long Short-Term Memory (LSTM) methods
which is able to exploit the interdependence in time-series over long periods
of time. This makes the predictions of electrical loads more accurate.
Considering
$\mathrm{x}=\left\\{\mathrm{x}_{1},\mathrm{x}_{2},\ldots,\mathrm{x}_{\mathrm{T}}\right\\}$
as input time series, and
$\mathrm{h}=\left\\{\mathrm{h}_{1},\mathrm{~{}h}_{2},\ldots,\mathrm{h}_{\mathrm{T}}\right\\}$
and
$\mathrm{y}=\left\\{\mathrm{y}_{1},\mathrm{y}_{2},\ldots,\mathrm{y}_{\mathrm{T}}\right\\}$
as output time series, which are created by updating the states of memory
cell, LSTM decides what dat to be removed from the cell gate [15]:

$\displaystyle\mathrm{f}(\mathrm{x})=$ (4)
$\displaystyle\sigma\left(w_{fx}x_{t}+w_{fh}h_{t-1}+w_{fc}c_{t-1}+b_{f}\right)$

In the next step, LSTM stores the new data as new cell state. this would be
the input gate layer named $i_{t}$ and formulated as the following:

$\mathrm{i}_{\mathrm{t}}=\sigma\left(\mathrm{w}_{\mathrm{ix}}\mathrm{x}_{\mathrm{t}}+\mathrm{w}_{\mathrm{ih}}\mathrm{h}_{\mathrm{t}-1}+\mathrm{w}_{\mathrm{ic}}\mathrm{c}_{\mathrm{t}-1}+\mathrm{b}_{\mathrm{i}}\right)$
(5)

A new vector is created, $\mathrm{U}_{\mathrm{t}}$, in order to store the new
data that are going to be added as the new state to the cell. This vector is
formed as the following equation:

$\mathrm{U}_{\mathrm{t}}=\mathrm{g}\left(\mathrm{w}_{\mathrm{cx}}\mathrm{x}_{\mathrm{t}}+\mathrm{w}_{\mathrm{ch}}\mathrm{h}_{\mathrm{t}-1}+\mathrm{b}_{\mathrm{c}}\right)$
(6)

Then, the old cell state $\mathrm{C}_{\mathrm{t}-1}$ is replaced by the new
cell state $\mathrm{c}_{\mathrm{t}}$, using $\mathrm{f}_{\mathrm{t}}$ and
$\mathrm{U}_{\mathrm{t}}$:

$c_{\mathrm{t}}=U_{t}{\mathrm{i}}_{\mathrm{t}}+\mathrm{c}_{\mathrm{t}-1}\mathrm{f}_{\mathrm{t}}$
(7)

The output of the output gate is represented by the following equation:

$o_{t}=\sigma\left(w_{0x}x_{t}+w_{oh}h_{t-1}+w_{0c}c_{t-1}+b_{0}\right)$ (8)

The new data, which is $\mathrm{h}_{\mathrm{t}}$, would be as follow:

$\mathrm{h}_{\mathrm{t}}=\mathrm{o}_{\mathrm{t}}\mathrm{l}\left(\mathrm{c}_{\mathrm{t}}\right)$
(9)

And the final output would equal $\mathrm{y}_{\mathrm{t}}$:

$y_{t}=k\left(w_{yh}h_{t}+b_{y}\right)$ (10)

Where $\mathbf{w}_{\mathrm{ix}}$, $\mathbf{w}_{\mathrm{fx}}$,
$\mathbf{w}_{\mathbf{ox}}$, and $\mathbf{w}_{\mathrm{cx}}$ are input matrices;
$\mathrm{w}_{\mathrm{ih}}$, $\mathrm{w}_{\mathrm{fh}}$,
$\mathrm{w}_{\mathrm{oh}}$, and $\mathrm{w}_{\mathrm{ch}}$ are recurrent
weight matrices; and $\mathbf{w}_{\mathrm{yh}}$ is the hidden output weight
matrix. The matrices $\mathbf{w}_{\mathbf{ic}}$, $\mathbf{w}_{\mathrm{fc}}$,
and $\mathrm{w}_{\mathrm{oc}}$ are connection weight matrices, and the vectors
$\mathrm{b}_{\mathrm{i}}$, $\mathrm{b}_{\mathrm{f}}$,
$\mathrm{b}_{\mathrm{o}}$, $\mathrm{b}_{\mathrm{c}}$, and
$\mathrm{b}_{\mathrm{y}}$ are based on bias vector [15].

### Deep Recurrent Neural Network

Deep Recurrent Neural Network (DRNN) has more layers in its network compare to
RNN and has more benefits comparatively. Figure 2 shows the structure of DRNN
and the functional process of a Deep RNN [12]. For a time-step, $\mathrm{t}$,
the parameters of the network neuron on the layer $\mathrm{l}$ updates their
values according to the following equations [12]:

$\mathrm{a}_{l}^{(\mathrm{t})}=\mathrm{b}_{l}+\mathrm{w}_{l}\mathrm{~{}h}_{l}^{(\mathrm{t}-1)}+\mathrm{U}_{l}\mathrm{x}^{(\mathrm{t})}$
(11)

$\mathrm{h}_{\mathrm{l}}^{(\mathrm{t})}=\text{activation
function}\left(\mathrm{a}_{\mathrm{l}}^{(\mathrm{t})}\right){\text{ for
}\mathrm{l}=1,2,\ldots,\mathrm{N}}$ (12)

$\mathrm{a}_{\mathrm{l}}^{(\mathrm{t})}=\mathrm{b}_{\mathrm{l}}+\mathrm{w}_{\mathrm{l}}\mathrm{h}_{\mathrm{l}}^{(\mathrm{t}-1)}+\mathrm{U}_{l}\mathrm{~{}h}_{l-1}^{(t)}{\text{
for }\mathrm{l}=2,3,...,N}$ (13)

$\mathrm{y}^{(\mathrm{t})}=\mathrm{b}_{\mathrm{N}}+\mathrm{w}_{\mathrm{N}}\mathrm{h}_{\mathrm{N}}^{(\mathrm{t}-1)}+\mathrm{U}_{\mathrm{N}}\mathrm{h}_{\mathrm{N}}^{(\mathrm{t})}$
(14)

$\mathrm{L}=\operatorname{loss}\mathrm{f}\left(\mathrm{y}^{(\mathrm{t})},\mathrm{y}_{\mathrm{target}}^{(\mathrm{t})}\right)$
(15)

where, $\mathbf{x}^{(\mathrm{t})}$ is the input of the time step $t$,
$\mathrm{y}^{(\mathrm{t})}$ is the predicted value, and
$\mathrm{y}_{\text{target }}^{(\mathrm{t})}$ is the output value.
$\mathrm{h}_{\mathrm{l}}^{(\mathrm{t})}$ is the activation function at layer
$l$ of the network in time step $t$. $a_{l}^{(t)}$ is the real input value on
layer $l$ in time step $t$. The characteristics of the DRNN enables it to
learn the uncertainties which are not very precedent.

Figure 2: Structure of DRNN.

##  III Proposed Method

In this section the proposed methodology is presented. Pooling strategy [12]
is used in the proposed method to tackle two main challenges in short-time
load forecasting for single house residential applications. Overfitting is one
of the major limits in deep learning-based load forecasting. Because of many
layers in deep networks, training of these networks needs a large amount of
data. Pooling stage can increase the amount of data greatly and avoid
overfitting. The second challenge is uncertainty of load consumption data. Due
to probabilistic nature of data, it is usually very difficult to model the
data. Although some of these uncertainties are due to extrinsic factors such
as weather conditions, diversity in training data can improve learning of the
uncertainties. So, in the proposed method, pooling strategy is used to augment
data and improve learning capabilities of the proposed method. A block diagram
of the proposed method is illustrated in Fig. 3. The first step is loading
data and data completion as discussed in IV-A.

After preprocessing stage, the data for each week is partitioned in a time-
based manner. Total number of 10,080 measurements are available for each week
with the resolution of one minute. This amount of data is divided to half day
intervals (N=24*60/2=720) that will generate M number of groups (M=14). Each
group will divided into training and test parts generating pools of data to
train and test the model. In the next stage, deep recurrent neural network
(DRNN) is trained to learn the load consumption profiles and finally the
performance of the proposed method is evaluated using test dataset and metrics
described in IV-C.

Figure 3: The proposed time pooling DRNN method.

##  IV Experiments and Results

###  IV-A Dataset and Preprocessing

In this study, the individual household electric power consumption dataset is
used to validate the proposed algorithm [16] 111Available online:
https://archive.ics.uci.edu/ml/datasets/individual+
household+electric+power+consumption. This dataset consists of 2,075,259
instance points that are measured over a period of almost 4 years. The global
active power (GAP) in kW, date and time columns of this dataset are used in
this paper. Preprocessing step based on averaging over available data is used
to complete missing data. Fig. 4 shows GAP measurements of two consecutive
weeks from the dataset exhibiting roughly similar pattern during a week. This
similarity is exploited in the proposed method.

Figure 4: Sample measurements from the dataset for two consecutive weeks
starting from 2006-12-18 and ending to 2006-12-31.

###  IV-B Experiments

All experiments are done using Python programming language on a machine with 8
GB of RAM and Intel core i5-5200U up to 2.7 GHz CPU. The parameters are set as
follows: N=720, M=14, training to test ratio: 67/23. Fig. 5 and 6 show the
training and test results of two methods DRNN and TPRNN (proposed method)
respectively.

Figure 5: Results of DRNN method for 2 weeks period.

Figure 6: Results of proposed method (TPRNN) for 4 weeks period.

###  IV-C Results and Discussion

To evaluate the proposed algorithm three widely used metrics are employed:
Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE). These metrics
are defined as follows:

$\mbox{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\hat{y_{i}})^{2}}$ (16)

$\mbox{MAE}=\frac{1}{N}\sum_{i=1}^{N}|y_{i}-\hat{y_{i}}|$ (17)

where $N$ is the total number of sample points. $y_{i}$ and $\hat{y_{i}}$ are
the true and estimated values respectively.

Results are compared with different algorithms including SVR, ARIMA, RNN and
DRNN in Table. I. The proposed method (TPRNN) exhibits better performance in
comparison with other state-of-the-art methods.

TABLE I: Comparison Results Based on RMSE and MAE. Method | RMSE | MAE
---|---|---
SVR | 0.96 | 0.77
ARIMA | 0.81 | 0.75
RNN | 0.75 | 0.55
DRNN | 0.39 | 0.20
TPRNN | 0.37 | 0.19

##  V Conclusion

Short-term load forecasting is an input or prerequisite for many applications
in power systems, specially in smart grids. Depending on the application,
different types of data are provided to run load forecasting algorithms. The
accuracy of the forecasting is important for different applications too. In
this paper, we presented a new short-term load forecasting method based on
time pooling DRNN to predict the power demand of single residential customers
for applications which require very short-term forecast. An application which
could use this method is nonintrusive load monitoring (NILM). A high
resolution historic data is used for training the model.

Load profiles demonstrate similarities in time. Based on these similarities,
in this paper time-pooling strategy in conjunction with deep recurrent neural
networks is used to solve the short-term load forecasting problem. Using
pooling strategy has two fold advantages. Firstly it satisfies the need for
large amount of data in deep learning methods and helps to model uncertainties
in the dataset precisely. Results demonstrates that the proposed method can
effectively forecast load consumption profile. Further research could be done
including more precise time pooling using clustering algorithms to put similar
power consumption patterns in a group.

## References

  * [1] A. Garulli, S. Paoletti, and A. Vicino, “Models and techniques for electric load forecasting in the presence of demand response,” _IEEE Transactions on Control Systems Technology_ , vol. 23, no. 3, pp. 1087–1097, 2014\.
  * [2] A. Estebsari and R. Rajabi, “Single residential load forecasting using deep learning and image encoding techniques,” _Electronics_ , vol. 9, no. 1, p. 68, 2020.
  * [3] R. Rajabi and A. Estebsari, “Deep learning based forecasting of individual residential loads using recurrence plots,” in _2019 IEEE Milan PowerTech_. IEEE, Conference Proceedings, pp. 1–5.
  * [4] H. N. Rafsanjani, “Factors influencing the energy consumption of residential buildings: a review,” in _Construction Research Congress 2016_ , Conference Proceedings, pp. 1133–1142.
  * [5] F. Mahia, A. R. Dey, M. A. Masud, and M. S. Mahmud, “Forecasting electricity consumption using ARIMA model,” in _2019 International Conference on Sustainable Technologies for Industry 4.0 (STI)_. IEEE, Conference Proceedings, pp. 1–6.
  * [6] A. Srivastava, A. S. Pandey, and D. Singh, “Notice of violation of ieee publication principles: Short-term load forecasting methods: A review,” in _2016 International conference on emerging trends in electrical electronics & sustainable energy systems (ICETEESES)_. IEEE, Conference Proceedings, pp. 130–138.
  * [7] M. Awad and R. Khanna, _Support vector regression_. Springer, 2015, pp. 67–80.
  * [8] T. Hossen, S. J. Plathottam, R. K. Angamuthu, P. Ranganathan, and H. Salehfar, “Short-term load forecasting using deep neural networks (DNN),” in _2017 North American Power Symposium (NAPS)_. IEEE, Conference Proceedings, pp. 1–6.
  * [9] A. Almalaq and G. Edwards, “A review of deep learning methods applied on load forecasting,” in _2017 16th IEEE international conference on machine learning and applications (ICMLA)_. IEEE, Conference Proceedings, pp. 511–516.
  * [10] E. Busseti, I. Osband, and S. Wong, “Deep learning for time series modeling,” _Technical report, Stanford University_ , pp. 1–5, 2012.
  * [11] H. Chen, C. A. Canizares, and A. Singh, “ANN-based short-term load forecasting in electricity markets,” in _2001 IEEE power engineering society winter meeting. Conference proceedings (Cat. No. 01CH37194)_ , vol. 2. IEEE, Conference Proceedings, pp. 411–415.
  * [12] H. Shi, M. Xu, and R. Li, “Deep learning for household load forecasting—a novel pooling deep rnn,” _IEEE Transactions on Smart Grid_ , vol. 9, no. 5, pp. 5271–5280, 2017.
  * [13] X. Sun, P. B. Luh, K. W. Cheung, W. Guan, L. D. Michel, S. S. Venkata, and M. T. Miller, “An efficient approach to short-term load forecasting at the distribution level,” _IEEE Transactions on Power Systems_ , vol. 31, no. 4, pp. 2526–2537, 2016.
  * [14] X. Sun, P. B. Luh, K. W. Cheung, W. Guan, L. D. Michel, S. Venkata, and M. T. Miller, “An efficient approach to short-term load forecasting at the distribution level,” _IEEE Transactions on Power Systems_ , vol. 31, no. 4, pp. 2526–2537, 2015.
  * [15] J. Zheng, C. Xu, Z. Zhang, and X. Li, “Electric load forecasting in smart grids using long-short-term-memory based recurrent neural network,” in _2017 51st Annual Conference on Information Sciences and Systems (CISS)_. IEEE, Conference Proceedings, pp. 1–6.
  * [16] G. Hebrail and A. Berard, “Individual household electric power consumption data set,” _É. d. France, Ed., ed: UCI Machine Learning Repository_ , 2012.