text:
\section{Introduction}
Discovering the underlying topics or discourses in large text corpora is
a challenging task in natural language processing (NLP).  A statistical
approach often starts by determining the frequency of occurrence of
terms across the corpus, and using the term frequency as a criterion for
word significance---a thesis put forward in a seminal paper by Luhn
\cite{Luhn:1958}.  From the list of terms ranked by frequency, terms
that are either too rare or too common are usually dropped, for they are
of little use.  For a domain-specific corpus, the top ranked terms in
the trimmed list often nicely summarize the main topics of the corpus,
as will be illustrated below.

For more detailed corpus analysis, such as discovering the subtopics
covered by the documents in the corpus, the term frequency list by
itself is, however, of limited use.  The main problem is that within a
given frequency range, function words, which primarily have an
organizing function and carry little or no meaning, appear together with
content words, which represent central features of texts and carry the
meaning of the context.  In other words, the rank of a term in the
frequency list is by itself not indicative of meaning \cite{Luhn:1958}.

This problem can be tackled by replacing the corpus-wide term frequency
with a more refined weighting scheme based on document-specific term
frequency \cite{Aizawa:2000}.  In such a scheme, a document is taken as
the context in which a word appears.  Since key words are typically
repeated in a document, they tend to cluster and to be less evenly
distributed across a text corpus than function words of the same
frequency.  The fraction of documents containing a given term can then
be used to distinguish them.  Much more elaborate statistical methods
have been developed to further explore the distribution of terms in
collections of documents, such as topic modeling \cite{Blei:2003} and
spacing statistics \cite{Ortuno:2002}.

An even more refined weighting scheme is obtained by reducing the
context of a word from the document in which it appears to a window of
just a few words.  Such a scheme is suggested by Harris' distributional
hypothesis \cite{harris54} which states ``that it is possible to define
a linguistic structure solely in terms of the `distributions' (=
patterns of co-occurrences) of its elements'', or as Firth famously put
it \cite{firth1957} ``a word is characterized by the company it keeps''.

Word co-occurrence is at the heart of several machine learning
algorithms, including the recently introduced word2vec by Mikolov and
collaborators \cite{MCCD13,MikolovSCCD13}.  Word2vec is a neural network
with a single hidden layer that uses word co-occurrence for learning a
relatively low-dimensional vector representation of each word in a
corpus, a so-called distributed representation \cite{hinton:1986}.  The
dimension is typically chosen of order 100 or 1000.  This is easily
orders of magnitude smaller than the size of a vocabulary, which would
be the dimension when a one-hot representation of words is chosen
instead.  Given the words appearing in a context, the neural network
learns by predicting (the representation of) the word in the middle, or
\textit{vice versa}.  During training, words that appear in similar
contexts are grouped together in the same direction by this unsupervised
learning algorithm.  The distributed representation thus ultimately
captures semantic similarities between words.  This has been
impressively demonstrated by a series of experiments in the original
word2vec papers, where semantic similarity was measured by the dot
product between \textit{normalized} vectors.

In this brief report, we consider the problem of identifying significant
terms that give information about content in text corpora made up of
short texts, such as abstracts of scientific papers, or news summaries.
It is proposed to use the L$_2$ norm, or length of a word vector, in
combination with the term frequency, as measure of word significance.

In a discussion forum dedicated to
word2vec,\footnote{http://groups.google.com/forum/\#!forum/word2vec-toolkit}
it has been argued by some that the length of a vector merely reflects
the frequency with which a word appears in the corpus, while others
argued that it in addition reflects the similarity of the contexts in
which a word appears.  According to this thesis, a word that is
consistently used in a similar context will be represented by a longer
vector than a word of the same frequency that is used in different
contexts.  Below, we provide experimental support for this thesis.  It
is this property that justifies measuring significance by word vector
length, for words represented by long vectors refer to a distinctive
context.

It is further proposed that the scatter plot of word vector length
versus term frequency of all the words in the vocabulary provides a
useful two-dimensional visualization of a text corpus.

The paper is organized as follows.  The next section introduces
the language corpus used and gives a global characterization based on
term frequency.  Section \ref{sec:exp} describes the experiments carried
out using word2vec, and presents the main results.  Section
\ref{sec:disc} concludes the paper with a short discussion.

\section{Dataset}
\label{sec:corpus}
For our experiments we use a dataset from the
arXiv,\footnote{http://arxiv.org/} a repository of scientific preprints.
The dataset consists of about 29k papers from one single subject class
in the arXiv, \textit{viz.}\ the hep-th section on theoretical
high-energy physics posted in the period from January 1992 to April
2003.  Although full papers are available, we consider only title and
abstract of the papers, which have about 100 word tokens on
average.\footnote{The dataset is available from the KDD Cup 2003
  homepage
  http://www.sigkdd.org/kdd-cup-2003-network-mining-and-usage-log-analysis}

LaTeX (or TeX) commands are removed from input through use of the detex
program.\footnote{http://www.cs.purdue.edu/homes/trinkle/detex/} The
input text is further converted to lowercase, and punctuation marks and
special symbols are separated from words, as was done in the
preprocessing step of a word2vec experiment by Mikolov on the IMDB
dataset of movie
reviews.\footnote{https://groups.google.com/forum/\#!msg/word2vec-toolkit/Q49FIrNOQRo/J6KG8mUj45sJ}

\subsection{Term Frequency List}
\begin{table}[t]
\begin{center}
\begin{tabular}{lrr}
term & $v$ & $tf$ \\ \hline
theory & 1.90 & 27702 \\
field & 2.00 & 17510 \\
gauge & 2.13 & 15536 \\
string & 2.33 & 13523 \\
model & 2.19 & 12389 \\
quantum & 2.14 & 12307 \\
theories & 2.22 & 10528 \\
space & 2.18 & 8035 \\
also & 1.67 & 7907 \\
models & 2.39 & 7313 \\
two & 2.03 & 7286 \\
fields & 2.11 & 7261 \\
solutions & 2.49 & 7129 \\
show & 1.87 & 7125 \\
action & 2.39 & 6602 \\
one & 1.82 & 6440 \\
black & 3.24 & 6011 \\
dimensions & 2.34 & 5953 \\
symmetry & 2.35 & 5792 \\
group & 2.46 & 5696 \\
equations & 2.54 & 5509 \\
algebra & 2.70 & 5461
\end{tabular}
\caption{Top ranked words in the term frequency list of the hep-th corpus
  with their vector length $v$ (included for later convenience) and term
  frequency $tf$.  Punctuation marks and stop words are removed from the
  list.}
\label{table:top}
\end{center}
\end{table}
After removing stop words and punctuation marks, the list of 50 most
frequently used words in the corpus reduces to the one given in
Table~\ref{table:top}.  Deriving from a domain-specific dataset, this
list indeed gives a succinct and fairly precise characterization of the
hep-th corpus, which is primarily about ``gauge theory'', ``quantum
field theories'', and ``string theory''.  It correctly reveals the
importance of ``models'' in this research area, as well as the
importance of the concepts ``space'', ``solutions'', ``action'',
``dimensions'', ``symmetry group'', ``equations'', and ``algebra''.  The
term ``black'' refers to ``black holes'', which play a distinctive role
in this corpus.  The term ``also'' is not filtered out by the
NLTK\footnote{www.nltk.org/} stop word list we use.  Finally, ``show''
(used exclusively as verb in this corpus, not as noun) appears mostly in
the context ``we show that'' and reveals that a large portion of the
corpus consists of research papers.  Note that ``show'' is the only verb
besides the stop word ``be'' that made it into the top 50 list.

\section{Experiments}
\label{sec:exp}
We next turn to word2vec.\footnote{The code is available for download at
  https:// code.google.com/p/word2vec} For training the neural network,
we use the same parameter settings as advertised for the IMDB dataset
referred to above.\footnote{Specifically, the parameters used are:
  \texttt{word2vec -train \$inputfile -output \$outputfile -cbow 0 -size
    100 -window 10 -negative 5 -hs 0 -sample 1e-4 -threads 40 -binary 0
    -iter 20 -min-count 1}} With these settings, the vector dimension is
100, the (maximum) context window size is 10, and the algorithm makes 20
passes through the dataset for learning.  The total number of tokens
processed by the algorithm is 3.2M.  As is typical for a highly specific
domain, the vocabulary is relatively small, containing about 44k terms,
of which about half is used only once.

\subsection{Similarity Distribution}
During training, similar words are grouped together in the same
direction by the learning algorithm, so that after training the vectors
encode word semantics.  One of the most popular measures of semantic
similarity in NLP is the cosine similarity given by the dot product
between two \textit{normalized} vectors.  Denoting the cosine of the
angle between the two vectors, the cosine similarity can take values in
the interval $[-1,1]$.

To analyze the hep-th corpus, we built a histogram of the cosine
similarity between arbitrarily chosen pairs of word vectors.  The words
are randomly selected from the vocabulary irrespective their frequency.
We have, however, discarded terms that appear only once.  The result,
given in Fig.~\ref{fig:sim}, is a bell-shaped distribution.
\begin{figure}
\begin{center}
  \includegraphics[width=.4\textwidth]{cos_sim_out_word2vecs.pdf}
\end{center}
\caption{Cosine similarity between arbitrarily chosen pairs of word
  vectors with $tf>1$.}
\label{fig:sim}
\end{figure}
To our surprise, the distribution is not centered around zero, but
around a positive value, $0.23$.  This means that word vectors in the
hep-th corpus have on average a certain similarity.  Closely related to
this is that the average word vector is non-zero, having a small length,
$v = 1.37$ ($v=1.51$ when words that only appear once are excluded).
This vector marks the center of the word cloud spanned in the word
vector space by all the words in the vocabulary.

To see if this behavior is shared by general purpose corpora, we
considered Wikipedia by way of example.\footnote{For a cleaned version
  of the Wikipedia corpus from October 2013, see
  https://blog.lateral.io/2015/06/the-unknown-perils-of-mining-wikipedia/}
For that corpus, covering diverse topics, we found, as expected, the
histogram to be centered around zero (and slightly right skewed).  The
non-zero value found for the hep-th corpus is therefore probably a sign
of the homogeneity of this dataset.

The reason for excluding terms that only appear once, which after all
make up half of the vocabulary, is that they have their own bell-shaped
distribution, slightly more peaked than the one shown in
Fig.~\ref{fig:sim} and centered at a higher cosine similarity of about
$0.5$.

Note the outlier at zero cosine similarity in the distribution in
Fig.~\ref{fig:sim}.  The reason for this outlier eludes us.

\subsection{Vector Length as Significance Factor}
To demonstrate that, besides depending on consistent use, the length of
a word vector also depends on term frequency, we consider the months of
the year, see Table~\ref{table:months}.  Apart from the word token
``may'', which in addition denotes a verbal auxiliary and is therefore
used in many different contexts, these terms consistently appear in the
abstracts of the hep-th corpus to indicate the time of a school or
conference where the paper was presented.  The data clearly show that
for fixed context, the vector length increases with term frequency, see
Fig.~\ref{fig:months}.
\begin{table}
\begin{center}
\begin{tabular}{lrr}
month     &   $v$  & $tf$ \\ \hline
january   &   4.16 & 16 \\
february  &   4.19 & 15 \\
march     &   4.90 & 37 \\
april     &   4.07 & 13 \\
may       &   2.23 & 2229 \\
june      &   5.95 & 73 \\
july      &   5.54 & 54 \\
august    &   5.10  & 31 \\
september &   5.54  & 51 \\
october   &   3.83 & 11 \\
november  &   4.35 & 15 \\
december  &   4.39 & 17
\end{tabular}
\caption{The months of the year with their word vector length $v$ and term
  frequency $tf$.}
\label{table:months}
\end{center}
\end{table}

\begin{figure}
\begin{center}
  \includegraphics[width=.4\textwidth]{months.pdf}
\end{center}
\caption{Word vector length as a function of frequency of appearance of
  the months of the year excluding ``may''.  The line through the data
  points serves as a guide to the eye.}
\label{fig:months}
\end{figure}
The three terms with the largest vector length are besides ``june'',
``school'' ($v=5.97, \, tf=114$) and ``conference'' ($v=5.89, \,
tf=93$).  If we take vector length as a measure of word significance,
this finding surprisingly supports another thesis by Luhn
\cite{Luhn:1958}, which states that:
\begin{quote}
The more often certain words are found in each other's company within a
sentence, the more significance may be attributed to each of these
words.
\end{quote}
Here, the phrase ``certain words'' refers to words that
\begin{quote}
a writer normally repeats [$\cdots$] as he advances or varies his
arguments and as he elaborates on an aspect of a subject.
\end{quote}

Table~\ref{table:months} also nicely demonstrates that term frequency
alone does not determine the length of a word vector.  The term ``may''
has a much higher frequency than the other terms in the table, yet it is
represented by the shortest vector.  This is because it is used in the
corpus mostly as a verbal auxiliary in opposing contexts.  When a word
appears in different contexts, its vector gets moved in different
directions during updates.  The final vector then represents some sort
of weighted average over the various contexts.  Averaging over vectors
that point in different directions typically results in a vector that
gets shorter with increasing number of different contexts in which the
word appears.  For words to be used in many different contexts, they
must carry little meaning.  Prime examples of such insignificant words
are high-frequency stop words, which are indeed represented by short
vectors despite their high term frequencies, see Table~\ref{table:stop}.
\begin{table}
\begin{center}
\begin{tabular}{lrr}
term     &   $v$  & $tf$ \\ \hline
the & 1.49 & 257866 \\
of & 1.51 & 148549  \\
. & 1.43 & 131022   \\
, & 1.41 & 84595    \\
in & 1.59 & 80056   \\
a & 1.61 & 72959    \\
and & 1.51 & 71170  \\
to & 1.61 & 53265   \\
we & 1.88 & 49756   \\
is & 1.69 & 49446   \\
for & 1.62 & 34970  \\
) & 2.03 & 28878
\end{tabular}
\caption{Top 12 terms in the term frequency list of the hep-th
  corpus with their word vector length $v$ and term frequency $tf$.  In
  addition to punctuation marks, this list exclusively features stop
  words.}
\label{table:stop}
\end{center}
\end{table}

\subsection{Vector Length vs. Term Frequency}
To study to what extent term frequency and word vector length can serve
as indicators of a word's significance, we represent all words in the
vocabulary in a two-dimensional scatter plot using these variables as
coordinates.  Figure~\ref{fig:norm_vs_freq} gives the result for the
hep-th corpus.
\begin{figure}[t]
\begin{center}
  \includegraphics[width=.4\textwidth]{length_vs_freq_word2vec.pdf}
\end{center}
\caption{Word vector length $v$ versus term frequency $tf$ of all words
  in the hep-th vocabulary.  Note the logarithmic scale used on the
  frequency axis.  The dark symbols denote bin means with the $k$th bin
  containing the frequencies in the interval $[2^{k-1}, 2^k-1]$ with
  $k=1,2,3,\ldots$.  These means are included as a guide to the eye.
  The horizontal line indicates the length $v = 1.37$ of the mean
  vector.}
\label{fig:norm_vs_freq}
\end{figure}
For given term frequency, the vector length is seen to take values only
in a narrow interval.  That interval initially shifts upwards with
increasing frequency.  Around a frequency of about 30, that trend
reverses and the interval shifts downwards.

Both forces determining the length of a word vector are seen at work
here.  Small-frequency words tend to be used consistently, so that the
more frequently such words appear, the longer their vectors.  This
tendency is reflected by the upwards trend in
Fig.~\ref{fig:norm_vs_freq} at low frequencies.  High-frequency words,
on the other hand, tend to be used in many different contexts, the more
so, the more frequently they occur.  The averaging over an increasing
number of different contexts shortens the vectors representing such
words.  This tendency is clearly reflected by the downwards trend in
Fig.~\ref{fig:norm_vs_freq} at high frequencies, culminating in
punctuation marks and stop words with short vectors at the very end.

Words represented by the longest vectors in a given frequency bin often
carry the content of distinctive contexts.  Typically, these contexts
are topic-wise not at the core of the corpus but more on the outskirts.
For example, the words with the longest vector in the high-frequency
ranges $[2^{k-1}, 2^k-1]$ with $k=9,10,11$, ``inflation'' ($v=4.64, \,
tf=571$), ``sitter'' ($v=3.81, \, tf=1490$) as in ``de Sitter'', and
``holes'' ($v=3.41, \, tf=2465$) as in ``black holes'', all refer to
general relativity.  General relativity, having its own subject class
(gr-qc) in the arXiv, is not one of the main subjects of the hep-th
corpus.  It takes a distinctive position in this corpus as it mostly
appears in studies that aim at reconciling general relativity with the
laws of quantum mechanics.
\subsection{POS Tagging}
To further assess the ability of word vector length to measure word
significance, we assign part-of-speech (POS) tags to each word in the
corpus.  For this task, we use the Stanford POS
tagger\footnote{Specifically, we use the
  english-caseless-left3words-distsim tagger.  For details on this model
  and for downloading, see
  http://nlp.stanford.edu/software/tagger.shtml.} \cite{Toutanova:2003}.
The final tag assigned to a word in the vocabulary is decided by
majority vote.

By the way that word2vec learns word representations, we expect nouns
(excluding proper nouns) and adjectives to be similarly distributed in
the $v$-$tf$ plane.  This is indeed what we observe, see
Fig.~\ref{fig:nn_vs_jj}.  Note that these word types pervade almost the
entire region covered in the full plot in Fig.~\ref{fig:norm_vs_freq} by
the complete vocabulary.
\begin{figure}
\begin{center}
  \includegraphics[width=.4\textwidth]{nn_vs_jj.pdf}
\end{center}
\caption{Word vector length $v$ versus term frequency $tf$ of all words
  in the hep-th vocabulary labeled nouns (red dots) or adjectives (blue
  dots).}
\label{fig:nn_vs_jj}
\end{figure}

We also find verbs and adverbs to be similarly distributed in the
$v$-$tf$ plane.  Again this was to be expected given that word2vec
learns word representations from word co-occurrences.  Somewhat
surprisingly, we observe in Fig.~\ref{fig:vb_vs_fw} that the
distribution of verbs also overlaps with that of function
words.\footnote{As function words we classify: prepositions (IN),
  pronouns (PRP, PRP\$, WP, WP\$), determiners (DT, PDT, WDT),
  conjunctions (CC), modal verbs (MD), and particles (RP). In brackets,
  we included here the tags used by the Stanford POS tagger.} These word
types no longer pervade the entire region covered in the full plot, but
are confined to the bottom band, corresponding to short vectors.  The
fact that function words are represented by short vectors underscores
the ability of vector length to measure word significance.
\begin{figure}[t]
\begin{center}
  \includegraphics[width=.4\textwidth]{vb_vs_fw.pdf}
\end{center}
\caption{Word vector length $v$ versus term frequency $tf$ of all words
  in the hep-th vocabulary labeled verbs (red dots) or function words
  (blue dots).}
\label{fig:vb_vs_fw}
\end{figure}

This proficiency is even more brought out by comparing the distribution
of function words and proper nouns, which typically are indicative of
distinctive contexts in the hep-th corpus.  The plot in
Fig.~\ref{fig:nnp_vs_fw} shows a clear separation of proper nouns and
function words for sufficiently large term frequencies.

\subsection{Visualization}
These results suggest an interesting technique for visualizing text
corpora.  By labeling the data points in the $v$-$tf$ scatter plot with
their terms, one obtains a two-dimensional visualization of all the
words in the vocabulary.  One advantage of a $v$-$tf$ plot is that words
are ranked by significance, thus allowing for effective exploration of a
corpus.  To deal with the large number of data points, an interactive
visualization tool can be build that allows the user to navigate a mouse
pointer over the plot, and that shows only the labels of the data points
in the vicinity of the pointer.

There exist other techniques for visualizing high-dimensional data, such
as the popular t-distributed stochastic neighbor embedding (t-SNE)
\cite{vdmaaten:2008}.  That machine learning algorithm, being an example
of multidimensional scaling, projects high-dimensional data points onto
a plane such that the distances, or similarities between them are
preserved as well as possible.  Words of similar meaning thus tend to be
projected together by the t-SNE algorithm.  Since the cosine similarity
is independent of vector lengths, word significance is ignored when
using this measure.  The t-SNE algorithm therefore arranges the data
points entirely differently from our proposal.  Moreover, in contrast to
the axes in the $v$-$tf$ scatter plot, those in the t-SNE plot have no
direct meaning.

\begin{figure}[t]
\begin{center}
  \includegraphics[width=.4\textwidth]{nnp_vs_fw.pdf}
\end{center}
\caption{Word vector length $v$ versus term frequency $tf$ of all words
  in the hep-th vocabulary labeled proper nouns (red dots) or function
  words (blue dots).}
\label{fig:nnp_vs_fw}
\end{figure}

\section{Discussion}
\label{sec:disc}
Most applications of distributed representations of words obtained
through word2vec so far centered around semantics.  A host of
experiments have demonstrated the extent to which the direction of word
vectors captures semantics.  In this brief report, it was pointed out
that not only the direction, but also the length of word vectors
carries important information.  Specifically, it was shown that word
vector length furnishes, in combination with term frequency, a useful
measure of word significance.  Also an alternative to the t-SNE
algorithm for projecting word vectors onto a plane was introduced, where
words are ordered by significance rather than similarity.

We have restricted ourselves to unigrams in this exploratory study.  For
more extended experiments and applications, including important bi- and
tri-grams into the vocabulary will certainly improve results.  We have
also restricted ourselves to running word2vec using parameters that were
recommended by the developers, and have not attempted to optimize them.

Finally, the question arises whether word vectors produced by other
highly scalable machine learning algorithms built on top of word
co-occurrences, such as the log bilinear model \cite{Mnih:2013} and GloVe
\cite{pennington:2014}, also encode word significance in their length.

\bibliographystyle{acl}


text:
\section{Introduction}

Stochastic gradient descent (SGD) is a widely used optimization technique in machine learning applications.
In distributed SGD, all $n$ workers collaboratively learn a global model $\mathbf{w}$ by minimizing the empirical loss with their local data:
\begin{align}
\label{eq:initialObjective}
\min_{\mathbf{w}\in \mathbb{R}^d} f(\mathbf{w}) :=
\frac{1}{n}\sum_{j=1}^n F_{j}(\mathbf{w}),
\end{align}
where $F_j(\cdot)$ is the local loss function of worker $j$. Traditionally, workers send their models or gradients to the central server after one local iteration, which is inefficient due to frequent model aggregations. This may cost high communication latency in practice.
Recently, a form of distributed SGD was proposed to
reduce communication cost by allowing multiple local iterations during one communication round \cite{McMahan2017a}.
This is referred to as local SGD \cite{stich2018local,Lin2020Don't}.
However, since data on workers can be non-IID, to avoid model divergence, the number of local iterations $P$ cannot be too large, which limits its advantage of reducing communication cost.

In practical scenarios, networks often have a hierarchical structure in nature, such as edge computing systems \cite{8270639} and software defined networks (SDN) \cite{6461195}. In these networks, workers often directly communicate with their local server rather than global server. Motivated by this hierarchical structure, a few works \cite{castiglia2020multi,liu2020client} proposed hierarchical SGD (H-SGD). In H-SGD, workers are partitioned into $N$ groups, where each group has a local server. To reduce communication cost, workers send their models to local servers to do several local aggregations before communicating with global server. They perform $I$ local iterations (referred to as local period) between local aggregations, and $G$ (referred to as global period) local iterations ($G>I$) between global aggregations.
In this paper, we call the connection between local servers and workers as ``downward'' network and the connection between local servers and the global server as ``upward'' network.

Recently, there have been a few works analyzing the convergence behavior of H-SGD.  A structure where the upward network is a peer-to-peer network while downward network is a server-worker network was considered by \citet{castiglia2020multi}. However, it only considers IID data. The work by \citet{liu2020client} considers non-IID data but it uses full-batch (non-stochastic) gradient descent and the convergence bound is an exponential function of $G$. In both works, a comprehensive comparison between H-SGD and local SGD is missing and none of them analyzes the effect of local aggregation on overcoming data heterogeneity, which is the key merit of H-SGD. Therefore, a more general and tighter analysis for H-SGD with non-convex objective function, non-IID data and stochastic gradient descent is needed. The effect of local aggregation on overcoming data heterogeneity lacks theoretical understanding.

\begin{table*}[t]
\centering
\caption{A summary of convergence bounds in the literature.}
\label{tab: comparison}
\small
\renewcommand{\arraystretch}{2}
\begin{tabular}{
>{\centering\arraybackslash}p{0.25\linewidth}|
>{\centering\arraybackslash}p{0.36\linewidth} | >{\centering\arraybackslash}p{0.05\linewidth} |
>{\centering\arraybackslash}p{0.03\linewidth} |
>{\centering\arraybackslash}p{0.08\linewidth} |
>{\centering\arraybackslash}p{0.08\linewidth}}
\hline
Paper&Convergence Bound& \!\!\!Non-IID\!\!\! & \!\!SGD\!\! & Type & %
\!\!\!Assumption\!\!\!\\
\hline
\citet{yu2019linear}& $O\left(\frac{1+\sigma^2}{\sqrt{nT}}+\frac{n}{T}(P\sigma^2+P^2\tilde{\epsilon}^2) \right)$ & \CheckmarkBold &\CheckmarkBold & \!\!\! Local SGD \!\!\!& $N=1$  \\
 \hline
\citet{liu2020client} & $O\left(\frac{1+B^G\tilde{\epsilon}^2}{\sqrt{nT}}\right)$ & \CheckmarkBold & \XSolidBrush & H-SGD & $\sigma^2=0$ \\
 \hline
\!\!\!\citet{castiglia2020multi}\!\!\!  &   $O\left(\frac{1+\sigma^2}{\sqrt{nT}}+ \frac{n}{T}\frac{G^2}{I}\sigma^2\right)$ & \XSolidBrush & \CheckmarkBold &  H-SGD & $\tilde{\epsilon}^2=0$ \\
 \hline
 Ours  &
$
O\left(\frac{1+\sigma^2}{\sqrt{nT}}+ \frac{(N-1)(G\sigma^2+G^2\tilde{\epsilon}^2) +(n-N)(I\sigma^2 + I^2\tilde{\epsilon}^2)}{T}\right)
    $
 & \CheckmarkBold & \CheckmarkBold &  H-SGD & None  \\
\hline
\end{tabular}
 \begin{tablenotes}
        \footnotesize
        \item \textsuperscript{1} $P$: aggregation period in local SGD; $G$: global aggregation period in H-SGD; $I$: local aggregation period in H-SGD; $N$: number of groups in H-SGD; $\tilde{\epsilon}^2$: global divergence (Assumption~\ref{assumption:divergence-sf-sgd}); $\sigma^2$: stochastic gradient noise
        \item \textsuperscript{2} $B$ is a constant and $B>2$. $n$ is the number of nodes and $N$ is the number of groups. $T$ is the total number of local iterations.
        \item \textsuperscript{3} Our bound can reduce to the IID case by setting $\tilde{\epsilon}^2=0$.
 \end{tablenotes}
\end{table*}

In this paper, to provide a better theoretical understanding to H-SGD,
we devise a novel characterization of the data heterogeneity of H-SGD by ``upward'' divergence and ``downward'' divergence, respectively. Furthermore, we show that the global divergence can be partitioned into upward and downward divergences.
H-SGD can be seen as performing distributed SGD both within a group and across groups. Within a group, workers perform multiple local iterations (on each worker) before local aggregation (within each group). Across groups, each group performs multiple local aggregations
before global aggregation.
With this characterization, we conduct a novel convergence analysis for H-SGD with non-IID data, non-convex objective function and stochastic gradients.
Furthermore,
we show that although data can be highly non-IID, local aggregation can help global convergence even when grouping is random.
Then by a detailed comparison with local SGD, we show that our convergence upper bound for H-SGD lies in between the convergence upper bounds of two single-level local SGD settings with aggregation periods of $I$ and $G$, respectively. This is referred to as the ``sandwich'' behavior, which reveals the fundamental impact of local aggregation.

Our convergence analysis shows that better convergence of H-SGD can be achieved with local aggregation when the number of groups, together with the global and local periods, are chosen appropriately to control the combined effect of upward  and downward divergences.
In general, we show that since local aggregation is more frequent, grouping strategies with a smaller upward divergence can strengthen the benefit of local aggregation. To reduce the communication cost while maintaining similar or better convergence, it can be beneficial to increase the global period $G$ and decrease the local period $I$.
We also extend our results to multi-level cases where there are multiple levels of local servers, where our characterization with upward and downward divergences can be applied to each level, from which we derive the convergence bound for general multi-level cases.

A comparison of our result with existing results is shown in Table~\ref{tab: comparison}.
When setting $N=1$ and $P=I=G$, our result recovers the well-known result for single-level local SGD by \citet{yu2019linear}. We can also see that choosing $I<G=P$ for H-SGD gives a smaller convergence upper bound, which shows the benefit of the hierarchical structure.
The result by \citet{liu2020client} only considers full gradient descent. In this case, the stochastic noise $\sigma^2 = 0$. Even when setting $\sigma^2=0$, our result is tighter. The result by \citet{castiglia2020multi} only considers IID case where the global divergence $\tilde{\epsilon}^2 =0$. Setting $\tilde{\epsilon}^2=0$, our result is still tighter than theirs since $I<G$. It can be seen that our result is the most general and tightest.

\textbf{Main Contributions.}
\begin{itemize}
    \item We introduce the new notion of ``upward'' and ``downward'' divergences to characterize data heterogeneity of H-SGD. We show that it can be extended %
    to multi-level cases.
    \item We derive a general convergence bound for two-level H-SGD with non-IID data, non-convex objective functions and stochastic gradient descent.
    \item We provide a novel convergence analysis for random grouping and show how local aggregation helps global convergence by a ``sandwich'' behavior.
    \item We extend our analysis to multi-level H-SGD and the result shows similar properties as the two-level case. To our knowledge, this is the first analysis which can be extended to multi-level cases.
\end{itemize}
We also conduct experiments on CIFAR-10, FEMNIST, and CelebA datasets. The results of experiments validate our theoretical results.

\section{Related Works}
Traditional distributed SGD is introduced by \citet{zinkevich2010parallelized}, which can be seen as a special case of local SGD with only one local iteration during one communication round.
There have been a large volume of works analyzing the convergence of local SGD, for convex objective functions \citep{Li2020On,WangJSAC2019}, non-convex objective functions \citep{NEURIPS2019_c17028c9,YuAAAI2019}, and their variants \citep{karimireddy2020scaffold,li2019federated,reddi2020adaptive,CooperativeSGD,yu2019linear}.
Local SGD can be regarded as a special case of H-SGD with only one level.
Our theoretical results for H-SGD recover results for both local SGD and traditional distributed SGD. There have been a few works analyzing the convergence of H-SGD, including \citet{castiglia2020multi,zhou2019distributed} for IID data, and \citet{liu2020client} for non-IID data with full-batch gradient descent, as described earlier. There are also works on system design for H-SGD without theoretical guarantees \citep{abad2020hierarchical,9127160}. %
In addition,
there are works on decentralized SGD \citep{CooperativeSGD,bellet2021d}, where workers exchange their models based on a doubly stochastic mixing matrix. However, the analysis for decentralized SGD cannot be applied to H-SGD, since the second largest eigenvalue of the mixing matrix in the case of H-SGD is one and the decentralized SGD analysis requires this second largest eigenvalue to be strictly less than one.

There are also some works focusing on practical aspects such as model compression and sparsification \citep{han2020adaptive,Peng2018,jiang2019model,Konecny2016} and partial worker participation \citep{Bonawitz2019,chen2020fedcluster}.
These algorithms and techniques are orthogonal to our work and may be applied together with H-SGD.

\section{H-SGD Setup}
In two-level H-SGD, all workers are partitioned into $N$ groups $\mathcal{V}_1,\mathcal{V}_2,\ldots ,\mathcal{V}_N$. The number of workers in each group is denoted by $n_i := |\mathcal{V}_i|$ ($i=1,2,\ldots ,N$). Then we have $n = \sum_{i=1}^N n_i$.
With this grouping, the objective function (\ref{eq:initialObjective}) is equivalent to
\begin{align}\label{eq:objective}
\min_{\mathbf{w}\in \mathbb{R}^d} f(\mathbf{w}) :=
\sum_{i=1}^N \frac{n_i}{n} f_{i}(\mathbf{w}),
\end{align}
where
$f_i(\cdot)$ is the averaged loss function of workers in group $i$ that is defined as follows:
\begin{align}
f_i(\mathbf{w}):=\frac{1}{n_i}\sum_{j \in \mathcal{V}_i} F_{j}(\mathbf{w}).
\end{align}
During each local iteration $t$, each worker $j$ updates its own model using SGD:
\begin{align}
\mathbf{w}_{j}^{t+1} = \mathbf{w}_{j}^t - \gamma \mathbf{g}(\mathbf{w}_{j}^t, \zeta^t_{j}),
\end{align}
where $\gamma$ is the learning rate, $\mathbf{g}(\mathbf{w}_{j}^t, \zeta^t_{j})$ is the stochastic gradient of $F_{j}(\mathbf{w})$, and $\zeta^t_{j}$ represents random data samples from the local dataset $\mathcal{D}_{j}$ at worker $j$.
We assume that $\mathbb{E}_{\zeta^t_{j}\sim \mathcal{D}_{j}}[\mathbf{g}(\mathbf{w}_{j}^t, \zeta^t_{j})] = \nabla F_{j}(\mathbf{w}_{j}^t)$.

During one communication round, local models are first averaged within group $i$ ($i = 1, 2,\ldots , N$) after every $I_i$ local iterations.
In particular, at local iteration $t \in \{I_i,2I_i,3I_i,\ldots \}$, we compute $\bar{\mathbf{w}}_i^{t}:=\frac{1}{n_i} \sum_{j \in \mathcal{V}_i} \mathbf{w}_{j}^{t}$. This can be done at a \textit{local server} (e.g., a computational component in close proximity) for group $i$.
After several rounds of ``intra-group'' aggregations, models from the $N$ groups are averaged globally. Let a global aggregation be performed for every  $G$ local iterations.
Then, at local iteration $t$, we have $\bar{\mathbf{w}}^{t}:=\frac{1}{N} \sum_{i=1}^{N} \bar{\mathbf{w}}_i^{t}$ for $t=G,2G,3G,\ldots $.
Note that we assume that all workers perform synchronous updates and let $G$ be a common multiple of $\{I_1, \ldots , I_N\}$. Therefore, distributed SGD is conducted both at the local level within a group  and the global level across groups.  We summarize the algorithm in Algorithm~\ref{alg:HF-SGD}, where $a\mid b$ (or $a\nmid b$) denotes that $a$ divides (or does not divide) $b$, i.e., $b$ is (or is not) an integer multiple of $a$.

In order to understand the fundamental convergence behavior of H-SGD, we will mainly focus on the two-level model introduced above. We will extend our analysis to %
more than two levels in a later section.

\begin{algorithm}[t]
\small
 \caption{Hierarchical SGD (H-SGD)}
 \label{alg:HF-SGD}
\KwIn{$\gamma$, $\bar{\mathbf{w}}^0$, $G$, $\{\mathcal{V}_i : i\in \{1,2,\ldots ,N\}\}$, $\{I_i : i\in \{1,2,\ldots ,N\}\}$}
\KwOut{Global aggregated model $\mathbf{\bar w}^{T}$}
\For{$t=0$ to $T-1$}{
\For{Each group $i\in \{1,2,\ldots ,N\}$, in parallel}{
\For{Each worker $j \in \mathcal{V}_i$, in parallel}{
Compute $\mathbf{g}(\mathbf{w}_{j}^t,\zeta^t_{j})$;\\
$\mathbf{w}_{j}^{t+1} \leftarrow \mathbf{w}_{j}^t - \gamma \mathbf{g}(\mathbf{w}_{j}^t,\zeta^t_{j})$;\\
}
\If{$I_i \mid t+1$}{
Local aggregate: $\mathbf{\bar w}_i^{t+1} \leftarrow \frac{1}{n_i} \sum_{j \in \mathcal{V}_i} \mathbf{w}_{j}^{t+1}$; \\
\If{$G \nmid t+1$}{
    Distribute: $\mathbf{w}_{j}^{t+1} \leftarrow \mathbf{\bar w}_i^{t+1}, \forall j\in\mathcal{V}_i$;
}
}
}
\If{$G \mid t+1$}{
Global aggregate: $\mathbf{\bar w}^{t+1} \leftarrow \frac{1}{N} \sum^N_{i=1} \mathbf{\bar w}_i^{t+1} $;\\
Distribute: $\mathbf{w}_{j}^{t+1} \leftarrow \mathbf{\bar w}^{t+1}, \forall j\in\mathcal{V}$;
}
}
\end{algorithm}

\section{Convergence Analysis for Two-level H-SGD}
\label{sec: Convergence Analysis}
We begin with a description of the assumptions made in our convergence analysis. Then, we present our results for general two-level H-SGD with fixed groupings. The number of workers and the number of local period during one communication round can vary among groups.

\subsection{Assumptions}
We make the following minimal set of assumptions that are common in the literature.

\begin{assumption}
\label{assumption:hf-sgd}
For H-SGD, we assume the following.

a) \textbf{Lipschitz gradient}
\begin{align}
\|\nabla F_{j}(\mathbf{w}) - \nabla F_{j}(\mathbf{w}')\| \le L\|\mathbf{w}-\mathbf{w}' \|, \forall j, \mathbf{w}, \mathbf{w}'
\end{align}
where $L$ is some positive constant.

b) \textbf{Bounded variance}
\begin{align}
\mathbb{E}_{\zeta^t_j\sim \mathcal{D}_{j} }\|\mathbf{g}(\mathbf{w};\zeta^t_j)-\nabla F_{j}(\mathbf{w})\|^2 \le \sigma^2, \forall j,\mathbf{w}
\end{align}
where $\mathcal{D}_{j}$ is the dataset at worker $j\in\mathcal{V}$.

c) \textbf{Bounded upward divergence (H-SGD)}
\begin{align}
\label{eq: div global}
	\frac{n_i}{n}\sum_{i=1}^N\|\nabla f_i(\mathbf{w}) - \nabla f(\mathbf{w})\|^2 \le \epsilon^2, \forall \mathbf{w}
\end{align}

d) \textbf{Bounded downward divergence (H-SGD)}
\begin{align}
\label{eq: bounded second-level divergence}
		\frac{1}{n_i}\sum_{j\in \mathcal{V}_i}\|\nabla F_{j}(\mathbf{w}) - \nabla f_i(\mathbf{w})\|^2 \le \epsilon_i^2,\forall i, \mathbf{w}
\end{align}
\end{assumption}

Note that the Lipschitz gradient assumption also applies to the group objective $f_i(\mathbf{w})$ and overall objective $f(\mathbf{w})$. For example, $\|\nabla f_i(\mathbf{w}) - \nabla f_i(\mathbf{w}')\|
= \|\frac{1}{n_i}\sum_{j\in \mathcal{V}_i}\nabla F_{j}(\mathbf{w}) -\frac{1}{n_i}\sum_{j\in \mathcal{V}_i}\nabla F_{j}(\mathbf{w}') \|
\le  \frac{1}{n_i} \sum_{j\in \mathcal{V}_i}\|\nabla F_{j}(\mathbf{w}) - \nabla F_{j}(\mathbf{w}')\|
\le  L\|\mathbf{w}-\mathbf{w}'\|$. Global divergence is often used to describe the data heterogeneity of single-level local SGD \cite{yu2019linear}, which is as the following.
\begin{assumption}{\textbf{Bounded global divergence}}
\label{assumption:divergence-sf-sgd}
\begin{align}
 \frac{1}{n}\sum_{j=1}^n\|\nabla F_{j}(\mathbf{w}) - \nabla f (\mathbf{w}) \|^2 \le \tilde{\epsilon}^2 ,\forall \mathbf{w}.
\label{eq:div_SF}
\end{align}

\end{assumption}
Note that when $N=1$, H-SGD reduces to local SGD and downward divergence becomes the global divergence while upward divergence becomes zero. When $N=n$, H-SGD also reduces to local SGD but upward divergence becomes the global divergence while downward divergence becomes zero. Here, we will use global divergence to show our results for H-SGD can encompass the results for local SGD.

Note for the global divergence, it is easy to show that
\begin{align}\label{eq:divergence partition}
\frac{1}{n}\sum_{j=1}^n &\|\nabla F_{j}(\mathbf{w}) \!-\! \nabla f (\mathbf{w}) \|^2 =\sum_{i=1}^N \frac{n_i}{n}\|\nabla f_i(\mathbf{w}) \!-\! \nabla f(\mathbf{w})\|^2 \notag \\
&+ \sum_{i=1}^N\frac{n_i}{n} \frac{1}{n_i}\sum_{j\in \mathcal{V}_i}\|\nabla F_{j}(\mathbf{w}) - \nabla f_i(\mathbf{w})\|^2.
\end{align}
In (\ref{eq:divergence partition}), we see that the upward and downward divergences are in fact a partition of the global divergence, which implies that upward and downward divergences do not increase at the same time. It is the hierarchical structure that makes this partition possible. Later we will show that this partition can exactly explain the benefits of local aggregation. We will discuss more about the relationship between upward/downward  and global divergences in a later section.

\subsection{Convergence Analysis}
\textbf{Technical Challenge.} The main challenge of the analysis is that workers only perform local iterations before global aggregation in local SGD, whereas in H-SGD, workers not only perform local iterations but also do aggregations with other workers in the same group. If we directly apply the analysis for local SGD in upward part, the effect of local aggregations would be neglected and the resulted bound would be loose. To address this, our key idea is to construct the local parameter drift $\|\mathbf{w}_j - \bar{\mathbf{w}}_i \|, j \in \mathcal{V}_i$ in the analysis of global parameter drift $\|\bar{\mathbf{w}}_i - \mathbf{\bar{w}}\|$, so that the analysis for the downward part can be incorporated in the analysis for the upward part.
\begin{theorem}\label{thm1}
Consider the problem in (\ref{eq:objective}). For any fixed worker grouping that satisfies Assumption~\ref{assumption:hf-sgd}, if  the learning rate in Algorithm $1$ satisfies $\gamma < \frac{1}{2\sqrt{6}GL}$,
then for any $T\ge 1$,
we have
\begin{subequations}
\begin{align}
&\frac{1}{T} \sum_{t=0}^{T-1}\mathbb{E}\|\nabla f(\mathbf{\bar w}^t) \|^2 \le \frac{2(f^0-f^*)}{\gamma T} +  \gamma L\frac{1}{n}\sigma^2\label{eq1a} \\
&\quad +2C\gamma^2G\frac{N-1}{n}\sigma^2
+3C\gamma^2G^2\epsilon^2
 \label{eq1b1} \\
&\quad +2C\gamma^2\sigma^2 \sum_{i=1}^N \frac{(n_i-1)I_i}{n} + 3C\gamma^2 \sum_{i=1}^N\frac{n_i}{n} I_i^2\epsilon_i^2 ,
\label{eq1d}
\end{align}
\end{subequations}
where $C=40/3$.
\end{theorem}

\textbf{Remark 1.}
Note that the bound (\ref{eq1a})--(\ref{eq1d}) can be partitioned into three parts. The terms in (\ref{eq1a}) are the original SGD part. If we set $N=n_i=1$ ($\epsilon$ and $\epsilon_i$ become zero), then only (\ref{eq1a}) will remain, which is the same as convergence bound for non-convex SGD in \cite{bottou2018optimization}. The terms in (\ref{eq1b1}) are the upward part, which consists of noise and upward divergence associated with $G$. The terms in (\ref{eq1d}) are the downward part, which has a similar form as (\ref{eq1b1}) associated with $I_i$. Divergence plays a more important role than noise since the coefficients in front of the SGD noise $\sigma^2$ include $G$ and $I_i$, while the coefficients in front of the divergences $\epsilon$ and $\epsilon_i$ include $G^2$ and $I_i^2$. Since $G>I_i,\forall i$, the upward part has a stronger influence on convergence.

\textbf{Remark 2.}
Note that all the divergences of H-SGD can be written as $O(\gamma^2G^2\epsilon^2 + \gamma^2 \sum_{i=1}^N \frac{n_i}{n} I_i^2 \epsilon_i^2)$  while the corresponding part in local SGD is $O(\gamma^2 G^2 \tilde{\epsilon}^2)$ \cite{yu2019linear}. This exactly shows how local aggregation overcomes data heterogeneity: global divergence is partitioned into two parts and local aggregation weakens the effects of the downward part. This also brings us some insights on how to group workers. %
The grouping strategy should have the smallest upward divergence since this can make full use of benefits of local aggregation. We will validate this property in the experiments.

\textbf{Remark 3.}
Let $\gamma = \sqrt{\frac{n}{T}} $ with $T\ge \frac{1}{24G^2L^2\sqrt{n}}$, when $T$ is sufficiently large, $\frac{1}{T} \sum_{t=0}^{T-1}\mathbb{E}\|\nabla f(\mathbf{\bar w}^t) \|^2
= O\left(\sqrt{\frac{1}{nT}}\right)+ O\left(\frac{1}{T}\right)$,
which achieves a linear speedup in $n$.

In the following corollary, we show that our results can encompass local SGD with local period $P$.
\begin{corollary}\label{corollary:single-level}
\textbf{(Degenerate to local SGD)}
Let $N=1$ and $\gamma \le \frac{1}{2\sqrt{6}PL}$,
from Theorem~\ref{thm1}, we obtain
\begin{align}
&\frac{1}{T} \sum_{t=0}^{T-1}\mathbb{E}\|\nabla f(\mathbf{\bar w}^t) \|^2 \le \frac{2(f^0-f^*)}{\gamma T} + \frac{\gamma L\sigma^2}{n} \notag\\
&\quad + 2C\gamma^2L^2\sigma^2 \left(1-\frac{1}{n}\right)P
+  3C\gamma^2L^2 P^2 \tilde{\epsilon}^2 \notag \\
&= O\left(\frac{1}{\gamma T}\right) + O\left(\frac{\gamma \sigma^2}{n}\right) \notag\\
&\quad + O\left(\gamma^2 P \sigma^2\left(1-\frac{1}{n}\right)\right)+ O\left(\gamma^2 P^2 \tilde{\epsilon}^2\right).
\label{eq: SF_SGD 1}
\end{align}
\end{corollary}
While (\ref{eq: SF_SGD 1}) is similar to the bound  by \citet{yu2019linear}, %
 we note that the third term $O\left(\gamma^2 P \sigma^2\left(1-\frac{1}{n}\right)\right) $ in the last equality in (\ref{eq: SF_SGD 1}) has an additional term of $\left(1-\frac{1}{n}\right)$ compared to the bound by \citet{yu2019linear}. This term can potentially make our bound tighter. Another important observation is that the techniques used to obtain this term is the key to make our H-SGD bound in Theorem~\ref{thm1} encompass original SGD cases.

\subsection{H-SGD with Random Grouping}
\label{sec: HF-SGD random grouping}
We now state our convergence results for H-SGD with random grouping.
We will show that the convergence upper bound of H-SGD with random grouping takes a value that is between the convergence upper bounds of two single-level local SGD settings with local and global periods of $I$ and $G$, respectively. This will provide insights on when and why local aggregation helps.

For worker grouping, we consider all possible grouping strategies with the constraint that $n_i = n/N$, $\forall i$.
Then, we uniformly select one grouping strategy at random.
Let the random variable $\mathsf{S}$ denote the uniformly random grouping strategy. This means that each realization of $\mathsf{S}$ corresponds to one grouping realization. First, we introduce two key lemmas for our convergence analysis.

\begin{lemma} \label{lemma1}
Using the uniformly random grouping strategy $\mathsf{S}$, for any $\mathbf{w}$, the average upward divergence is %
\begin{align}
\label{eq: lemma 1}
\mathbb{E}_{\mathsf{S}} \left[\frac{1}{N}\sum_{i=1}^N \left\| \nabla f(\mathbf{w}) \!-\! \nabla f_i(\mathbf{w}) \right\|^2\right]
\le  \left(\frac{N\!-\!1}{n\!-\!1}\right) \tilde{\epsilon}^2,
\end{align}
where $\tilde{\epsilon}$ is given in (\ref{eq:div_SF}).
\end{lemma}

\begin{lemma} \label{lemma2}
Using the uniformly random grouping strategy $\mathsf{S}$, for any $\mathbf{w}$, the average downward divergence is %
\begin{align}
\label{eq: lemma 2}
\mathbb{E}_{\mathsf{S}}\left[\frac{1}{n}\sum_{i=1}^N \sum_{\mathsf{k} \in \mathcal{V}_i} \left\|\nabla f_i(\mathbf{w}) \!-\! \nabla F_{\mathsf{k}}(\mathbf{w})\right\|^2 \right] \le \left(1\!-\!\frac{N\!-\!1}{n\!-\!1}\right)\tilde{\epsilon}^2.
\end{align}
\end{lemma}

Similar to (\ref{eq:divergence partition}), Lemma~\ref{lemma1} and Lemma~\ref{lemma2} show that the sum of upper bounds for upward and downward divergences is equal to the global divergence. Furthermore, when the number of groups increases, upward divergence becomes larger while downward divergence becomes smaller. When grouping is random, $N$ is the key parameter that determines how global divergence is partitioned. For simplicity, we let each group have the same local period $I_i =I, \forall i$. Then, we can obtain the following theorem.

\begin{theorem}
\label{thm2}
Using the uniformly random grouping strategy $\mathsf{S}$, let $\gamma \le \frac{1}{2\sqrt{6}GL}$, then we have
\begin{align}
&\mathbb{E}_{\mathsf{S}}\left[\frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\|\nabla f(\mathbf{\bar w}^t) \|^2\right]
\le \frac{2(f^0-f^*)}{\gamma T} + \frac{\gamma L\sigma^2}{n}  \notag\\
&\quad+2C\gamma^2L^2\left[\left(\frac{N-1}{n}\right)G+\left(1-\frac{N}{n}\right)I\right]\sigma^2 \notag\\
&\quad+3C \gamma^2L^2\!\left[\!\left(\frac{N\!-\!1}{n\!-\!1}\right)\!G^2\! +\! \left(\!1-\frac{N\!-\!1}{n\!-\!1}\right)\!I^2\right] \tilde{\epsilon}^2,
\label{eq:thm2}
\end{align}
where $C = 40/3$.
\end{theorem}
\textbf{Remark 4.} Note that both multiplicative factors of the noise term with $\sigma^2$ and the divergence term with $\tilde{\epsilon}^2$ are composed of two parts, where the upward part is ``modulated'' by $G$ while the downward part is ``modulated'' by $I$. As $N$ becomes larger, the upward part has a stronger influence on the convergence bound since both $\frac{N-1}{n}$ and $\frac{N-1}{n-1}$ become larger. Note that the convergence upper bound of H-SGD can be sandwiched by the convergence upper bounds of two local SGD. To see this,
we consider the following three scenarios: 1) local SGD with    aggregation period $P=G$, 2) local SGD with aggregation period $P=I$, and 3) H-SGD with local aggregation period $I$ and global aggregation period $G$.
We let them all start with the same $\mathbf{\bar w}^0 $ and choose learning rate $\gamma$ such that $\gamma \le \frac{1}{2\sqrt{6}GL}$. We can see that all these three cases have the same first two terms in (\ref{eq:thm2}). However, for the third and fourth terms in (\ref{eq:thm2}), we have
\begin{align}
\left(1-\frac{1}{n}\right)&I
\le \left(\frac{N\!-\! 1}{n}\right)G + \left(1-\frac{N}{n}\right)I
\le \left(1-\frac{1}{n}\right)G,    \label{eq: remark 21} \\
I^2  \le &\left(\frac{N-1}{n-1}\right)G^2 + \left(1-\frac{N-1}{n-1}\right)I^2 \le G^2.  \label{eq: remark 22}
\end{align}
Equations~(\ref{eq: remark 21}) and (\ref{eq: remark 22}) show that the convergence upper bound of H-SGD has a value between the convergence upper bounds of local SGD with aggregation periods of $P=I$ and $P=G$, respectively. The grouping approach can explicitly characterize how much the convergence bound moves towards the best case, i.e., local SGD with $P=I$. However, %
this case incurs the highest global communication cost, so grouping can adjust the trade-off between convergence and communication cost.

\textbf{Remark 5.} Even with a large $G$, if we make $I$ sufficiently small, the convergence bound of H-SGD can be improved. To see this, consider the last two terms of (\ref{eq:thm2}).
Suppose $G = mI, m=1,2, \cdots$. For a non-trivial worker grouping, i.e., $1<N<n$, if we increase $G$ to $G' = lG, 1<l<\sqrt{\frac{1}{m^2}\frac{n-N}{N}+1}$ and decrease $I$ to $I' = q I, q\le \sqrt{1-m^2(l^2-1)\frac{N}{n-N}}$, then one can  show that the bound (\ref{eq:thm2}) using $G'$ and $I'$ can be lower than %
that using $G$ and $I$.
A similar behavior can also be seen for fixed grouping, as shown empirically in Figure~\ref{fig:largeG} in the experiments.

\section{H-SGD with More Than Two Levels}
\label{sec:multi-level}
In the two-level H-SGD setting discussed in the previous sections, there is only one level of local servers between the global server and workers.
In this section, we  extend our analysis to H-SGD with more than one level of local servers.
Specifically, as shown in Figure~\ref{fig:3-level-case}, we consider a total of $M \geq 2$ levels.
For $\ell=1,2,\ldots ,M-1$, each server at level $\ell-1$ is connected to $N_\ell$ servers at level $\ell$, where we assume that the global server is at a ``dummy'' level $\ell=0$ for convenience.
For $\ell = M$, each server at level $M-1$ directly connects to $N_M$ workers.
As a result, we have a total of $n = \Pi_{\ell=1}^M N_\ell$ workers.
With this notation, a sequence of indices $(k_1, k_2, \ldots , k_\ell)$ denotes a ``path'' from the global server at level $0$ to a local server/worker at level $\ell$, where this path traverses the $k_1$-th server at level $1$, the $k_2$-th server at level $2$ connected to the $k_1$-th server at level $1$, and so on. Hence, the index sequence $(k_1, k_2, \ldots , k_\ell)$ uniquely determines servers (and workers if $\ell=M$) down to level $\ell$.

\begin{figure}[t]
\centering
\includegraphics[width=1\columnwidth]{multi-level.pdf}
\caption{{Three-level example with $N_1=N_2=2, N_3=3$.
}}
\label{fig:3-level-case}
\end{figure}

We assume in the multi-level case that each server at level $\ell -1$ connects to the same number ($N_\ell$) of servers/workers in the next level $\ell$, for ease of presentation. Our results can be extended to the general case with different group sizes.

Let $F_{k_1\ldots  k_M}(\mathbf{w})$  denote the objective function of worker $k_1  \ldots  k_M$ and $f_{k_1 \ldots  k_\ell}(\mathbf{w})$ denote the averaged objective function of workers connected to server $k_1 \ldots  k_\ell$,
\begin{align}
\!\!\!\!f_{k_1 \ldots  k_\ell}(\mathbf{w})\! :=\! \frac{1}{\prod_{j=\ell+1}^M N_j}\sum_{k_{\ell+1}=1}^{N_{\ell+1}}\!\!\!\cdots\!\!\!\sum_{k_M=1}^{N_M} \!\!\!F_{k_1 \ldots  k_M}(\mathbf{w}).  \!
\end{align}
Then the global objective function can be rewritten as
\begin{align}
\label{eq: f}
f(\mathbf{w}) &= \frac{1}{\prod_{j=1}^M N_j}\sum_{k_1=1}^{N_1}\cdots\sum_{k_M=1}^{N_M} F_{k_1 \ldots  k_M}(\mathbf{w})    \notag\\
&= \frac{1}{\prod_{j=1}^\ell N_j}\sum_{k_1=1}^{N_1}\cdots\sum_{k_\ell=1}^{N_{\ell}} f_{k_1 \ldots  k_\ell}(\mathbf{w}).
\end{align}
Let $P_{\ell}$ ($\ell=1,2,\ldots ,M$) denote the period (expressed as the number of local iterations)
that the parameters at servers at level $\ell$ are aggregated by their parent server at level $\ell - 1$.
We require that $P_1\!>\!P_2\!>\!\ldots \!>\!P_M$ and $P_{\ell}$ is an integer multiple of $P_{\ell-1}$. The algorithm is in Appendix~\ref{appendix:MultiLevel}.

In the following, we provide  results for the random grouping case. Results for the fixed grouping case can be found in Appendix~\ref{appendix:MultiLevel}.
We apply the divergence partition idea to each level, so we extend Lemma~\ref{lemma1} and Lemma~\ref{lemma2} to obtain the following lemma.
\begin{lemma}
\label{lemma3}
Using the uniformly random grouping strategy $\mathsf{S}$, the $\ell$-th level averaged upward and downward divergences are given by
\begin{align}
&\mathbb{E}_\mathsf{S} \bigg[ \frac{1}{n_\ell}\!\sum_{k_1=1}^{N_1}\!\!\!\!\cdots\!\!\!\sum_{k_\ell=1}^{N_\ell}\!\!\|\nabla f_{k_1 \ldots  k_\ell}(\mathbf{w})\!-\!\nabla f(\mathbf{w}) \|^2\bigg] \le\! \bigg(\frac{n_\ell-1}{n-1} \bigg)\tilde{\epsilon},\\
&\mathbb{E}_\mathsf{S} \bigg[\frac{n_\ell}{n}\!\sum_{k_{\ell+1}=1}^{N_{\ell+1}}\!\!\!\!\cdots\!\!\!\sum_{k_M=1}^{N_M}\!\!\!\|\nabla f_{k_1 \ldots  k_\ell}(\mathbf{w})\!-\!\nabla F_{k_1 \ldots  k_M}(\mathbf{w}) \|^2 \bigg]\notag\\
&\le\! \bigg(1-\frac{n_\ell-1}{n-1} \bigg)\tilde{\epsilon}^2 \label{eq:multi_level_local_divergence}
\end{align}
respectively,
$\forall\mathbf{w}$, $\forall k_1,\ldots ,k_\ell$ in (\ref{eq:multi_level_local_divergence}), where $n_\ell = \Pi_{j=1}^\ell N_j$ and $\tilde{\epsilon}$ is the global divergence.
\end{lemma}
When the number of servers in the $\ell$th-level $n_\ell$ increases, the $\ell$th-level's upward divergence becomes larger while its downward divergence becomes smaller. The sum of the $\ell$th-level upward and downward divergences are upper bounded by the global divergence. Based on this lemma, we derive the convergence bound for multi-level case as the following.
\begin{theorem}
\label{thm3}
Consider uniform random grouping strategy $\mathsf{S}$, if $\gamma\le \frac{1}{2\sqrt{6}P_1L}$, then for multi-level case,
\begin{align}
\label{eq: coro 3}
&\frac{1}{T}\sum_{t=0}^{T-1} \mathbb{E}_{\mathsf{S}} \bigg\|\nabla f(\bar{\mathbf{w}}^t) \bigg\|^2
\le  \frac{2}{\gamma T}\bigg(f^0-f^* \bigg) + \frac{\gamma L \sigma^2}{n} \notag\\
&\quad + C\gamma^2L^2\frac{1}{M-1}\! \cdot \!\sum_{\ell=1}^{M-1}\left[2A_1(\ell)\sigma^2 +3A_2(\ell)\tilde{\epsilon}^2\right]
\end{align}
where, $A_1(\ell) := P_1\Big(\frac{1}{\Pi_{j=\ell}^M N_j}-\frac{1}{n}\Big)
+P_{\ell}\Big(1-\frac{1}{\Pi_{j=\ell}^M N_j}\Big)$, $A_2(\ell) := P_1^2\Big(\frac{n_\ell - 1}{n-1}\Big)+P_{\ell}^2\Big(1-\frac{n_\ell - 1}{n-1}\Big)$.
\end{theorem}
\textbf{Remark 6.} Note that when $M=2$, this bound reduces to the bound given in Theorem~\ref{thm2}. The effect of the $\ell$-th level can be represented as $2A_1(\ell)\sigma^2 +3A_2(\ell)\tilde{\epsilon}^2$, where each $A_i(\ell), i=1,2$ is composed of an upward part ``modulated'' by $P_1$ and a downward part ``modulated" by $P_\ell$. %
Similar to
(\ref{eq: remark 21}) and (\ref{eq: remark 22}), we can obtain the sandwich behavior as:
\begin{align}
\Big(1 - \frac{1}{n}\Big) P_M &\leq \frac{1}{M - 1}\sum_{\ell=1}^{M-1} A_1(\ell) \leq \Big(1 - \frac{1}{n}\Big) P_1, \!\\
P_M^2 &\leq \frac{1}{M-1}\sum_{\ell=1}^{M-1} A_2(\ell) \leq P_1^2.
\end{align}
It can be seen that the convergence upper bound of H-SGD with more than two levels also takes a value that is between the convergence upper bounds of local SGD with local periods $P_1$ and $P_M$, respectively. Compared to the two-level case, this provides greater freedom to choose  parameters in $A_1(\ell)$ and $A_2(\ell)$, for $\ell = 1,\ldots ,M-1$.

\section{Experiments}
\label{sec:experiments}

\begin{table*}[t]
\centering
\caption{Total time (s) needed to achieve 50\% test accuracy for VGG-11 with CIFAR-10.}
\vspace{-0.2cm}
\label{tab: wallclocktime}
\small
\begin{tabular}{cccccc}
\hline
Case
& $P=5$  & $P=10$ & $P=50$ &$G=50, I=5$& $G=50, I=10$ \\
\hline
Time (s) & $673.5$ & $690.2$ & $944.3$ & $160.3$ & $381.1$ \\
\hline
\end{tabular}
\end{table*}

\begin{figure*}[h]
\centering
\begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{new-femnist.png}

  \caption{FEMNIST with CNN.}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{new-celeba.png}
  \caption{CelebA with VGG-11.}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{new-cifar-comm.png}
  \caption{CIFAR-10 with VGG-11.}
  \end{subfigure}
\vspace{-0.2cm}
\caption{Test accuracy v.s. communication time ($N=2$).}
\label{fig:commtime}

\end{figure*}

\begin{figure*}[tb]
  \centering
  \begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{new-cifar-noniid.png}
  \caption{Non-IID ($G=50$ for two-level case).}
  \label{fig:nIID}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{new-cifar-largerG.png}
  \caption{Non-IID ($G=200$ for two-level case).}
  \label{fig:largeG}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
  \centering
  \includegraphics[width=4.8cm]{cifar-giid-test-newv2.png}
  \caption{Group-IID/Group-non-IID.}
  \label{fig:gIID}
  \end{subfigure}
  \vspace{-0.1cm}
  \caption{Results with CIFAR-10. Test accuracy v.s. local iterations.
  By default, $N=2$.}
  \label{fig: cifar 10}
  \vspace{-0.3cm}
\end{figure*}

In this section, we validate our theoretical results with experiments on training the VGG-11 model over CIFAR-10~\citep{krizhevsky2009learning} and CelebA~\citep{liu2015faceattributes}, and training a convolutional neural network (CNN) over FEMNIST~\citep{cohen2017emnist}, all with non-IID data partitioning across workers. The communication time is emulated by measuring the round-trip time of transmitting the model between a device (in a home) and local (near)~/ global (far) Amazon EC2 instances.
For the computation time, we measured the averaged computation time needed for VGG-11 during one iteration, where we ran SGD with VGG-11 on a single GPU for 100 times and then computed the averaged computation time, which is approximately 4 ms per iteration on each worker.
The experiments presented in this section are for \textit{two-level} H-SGD. Additional details of the setup and experiments with \textit{more levels} and \textit{partial worker participation} can be found in Appendix~\ref{appendix:AdditionalExperiments}. Our code is available at \url{https://github.com/C3atUofU/Hierarchical-SGD.git}.

In all plots, curves with $P$ (denoting the aggregation period) are for local SGD and curves with $G$ (global period), $I$ (local period), and $N$ (number of groups) are for H-SGD. Figure~\ref{fig:commtime} shows
that for all the datasets and models evaluated in the experiments, H-SGD can achieve a better accuracy than local SGD given the same communication time.
In other words, for any given accuracy, H-SGD can achieve this accuracy using shorter or similar communication time. In Table~\ref{tab: wallclocktime}, we provide the total time (including both communication and computation) needed for VGG-11 with CIFAR-10 to achieve $50\%$ test accuracy (since $50\%$ is nearly the best accuracy for $P=50$). We can see H-SGD cases ($G=50, I=5$ and $G=50, I=10$) perform much better than single-level local SGD cases. This is because communication time plays a much more important role and H-SGD can achieve 50\% accuracy within less amount of time in total.

To show the ``sandwich'' behavior, we examine the test accuracy as a function of the total number of local iterations for the non-IID scenario in Figure~\ref{fig:nIID}.  We can observe that the performance of H-SGD with global period $G$ and local period $I$ is between that of local SGD with $P=G$ and with $P=I$. We can also observe that as $N$ becomes larger, the performance becomes worse since the upward divergence becomes larger, which is consistent with our analysis in Remark~4. Comparing Figures~\ref{fig:nIID} and \ref{fig:largeG}, we see that decreasing $I$ while increasing $G$ can even improve the performance, which is consistent with Theorem~\ref{thm2}.
For example, $G=200,I=2$ gives a better performance than $G=50,I=5$, while $G=50,I=5$ is similar to $P=10$.
This shows that by allowing more local aggregations, H-SGD can reduce the number of global aggregations by $95\%$ ($G=200$ vs. $P=10$) while maintaining a similar performance to local SGD. In Figure~\ref{fig:gIID}, we show the effects of grouping corresponding to our analysis for Theorem~\ref{thm1}. Group-IID is a grouping strategy which makes upward divergence nearly zero while group-non-IID is a grouping strategy with a large upward divergence. We can see that group-IID performs as well as group-non-IID after reducing $I$ by half.

\section{Conclusion}
We have studied H-SGD with multi-level model aggregations. In particular, we have provided a thorough theoretical analysis on the convergence of H-SGD over non-IID data, under both fixed and random worker grouping. We have successfully answered the important question on how local aggregation affects convergence of H-SGD.
Based on our novel analysis of the local and global divergences, we established
explicit comparisons of the convergence rate
for H-SGD and local SGD.
Our theoretical analysis provides valuable insights into the design of practical H-SGD systems, including the choice of global and local aggregation periods. Different grouping strategies are considered to best utilize the divergences to reduce communication costs while accelerating learning convergence. In addition, we have extended the analysis approach to H-SGD with arbitrary number of levels.
Future work could theoretically analyze the effect of partial worker participation in H-SGD.

\clearpage

\section*{Acknowledgment}
    This research was partly sponsored by the U.S. Army Research Laboratory and the U.K. Ministry of Defence under Agreement Number W911NF-16-3-0001. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Army Research Laboratory, the U.S. Government, the U.K. Ministry of Defence or the U.K. Government. The U.S. and U.K. Governments are authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation hereon.

\input{main.bbl}

\clearpage

\onecolumn

\newcommand{\mathbb{E}}{\mathbb{E}}
\newcommand{\mathbf{w}}{\mathbf{w}}
\newcommand{\mathbf{g}}{\mathbf{g}}

\begin{center}
\LARGE \textbf{Appendix}
\end{center}

redpajama_set_name: RedPajamaArXiv