text:
\section{Introduction} \label{sec:intro}
Mean-flow-based linear analyses have been used since the 1970s to understand and model the dynamics of coherent structures in turbulent shear flow. The original idea behind this is that turbulence generates a mean flow that can be seen as an equivalent laminar flow on which disturbances evolve \citep{crighton_jfm_1976}. The mean flow includes some, but not all, of the effects of the non-linear flow dynamics. Approaches of this kind have evolved considerably in recent years as global stability \citep{lesshafft_jfm_2006,akervik_ejm_2008,sipp_amr_2010,mantic_prl_2014} and input-output, or resolvent analyses \citep{mckeon_jfm_2010,hwang_jfm_2010,towne_jfm_2018,schmidt_jfm_2018,cavalieri_amr_2019}, have become feasible thanks to progress in computational methods. Such analysis has been applied to a broad variety of fluid-mechanics problems (compressible, incompressible, wall-bounded and free shear flows). For the example of a turbulent jet considered later in this paper, recent reviews  have been provided by \citet{jordan_arm_2013} and \citet{cavalieri_amr_2019}.

Many of these studies use the mean-flow-based linear operator to analyse the flow dynamics in an input-output, or resolvent, framework \citep{mckeon_jfm_2010,towne_jfm_2018}. In this case, the forcing terms are considered to drive the response through the resolvent operator. It is shown that the resolvent operator fully represents the flow dynamics if forcing is white  \citep{towne_ctr_2016,lesshafft_prf_2019}. For coloured forcing, which is generally the case for turbulent flows, flow dynamics depend also on the spectral content of the forcing \citep{zare_jfm_2017}. \cite{beneddine_jfm_2016} discussed some situations in which flow dynamics may be partially governed by the linear operator.

In other studies, eddy viscosity models have been included in the linear operator. This can enhance the extent to which the resolvent modes match the observed turbulence structure \citep{cossu_jfm_2009,hwang_jfm_2010,morra_jfm_2019,pickering_aiaa_2019}. However, physical interpretation of the remaining, unmodelled forcing terms then becomes unclear, as the use of an eddy viscosity amounts to a partial modelling of the effects of nonlinear forcing from Reynolds stress fluctuations.

We consider mean-flow-based linear analyses from the point of view of their uniqueness. It is obvious that modifying the linear operator, either through eddy-viscosity-based modelling, or by directly changing the linearisation point, will change the characteristics of the linearised system. A less obvious ambiguity is investigated in the present paper: we aim to show that, by choosing two different, nonlinearly related, variable sets that define a given flow, one may obtain two linearised systems with different characteristics, even when the transformation between the variables is bijective. Part of this ambiguity comes from linearisation around the mean, which is a not a fixed point of the system for turbulent flows. Another part of the ambiguity is due to the non-equivalence of the means obtained for these two variable sets. A common example is the nonlinear transformation between primitive and conservative variables. We explore the effect of this choice on resolvent analysis.

We first provide (\S\ref{sec:nstransform}) the mathematical framework that relates two linear operators obtained via nonlinear transformation of the dependent variables. We illustrate (\S\ref{sec:testcases}) the analysis by considering a number of Large-Eddy Simulation (LES) datasets for turbulent jets ranging from isothermal subsonic to heated supersonic using conservative and primitive variable sets. Based on our analysis and observations, we argue (\S\ref{sec:conc}) that the properties of the resolvent operator cannot be regarded as universal, but instead depend on the choice of variables used to define the mean, and that the overall model becomes independent of this choice only if the nonlinear forcing terms are appropriately modelled.
\section{Nonlinear transformation of Navier-Stokes equations} \label{sec:nstransform}
This study originated from an issue faced by the authors while trying to use the forcing data from an LES database constructed using a primitive-like variable set \citep{towne_phd_2016,bres_aiaa_2017,bres_jfm_2018}, with a resolvent analysis tool written in conservative variables \citep{bugeat20193d}. The problem can be illustrated by considering the effect of a nonlinear transformation of variables on the Navier-Stokes equations (though the analysis in fact applies to any nonlinear dynamical system).
\subsection{Derivation}\label{subsec:derivation}
Navier-Stokes equations (NSE) are written,
\begin{align} \label{eq:qeq}
\partial_t \mathbf{q} = \mathcal{N}(\mathbf{q}),
\end{align}
where $\mathbf{q}\in \mathbb{R}^m$ is the vector of state variables with $m=5$ for three dimensional flows and $\mathcal{N}:\mathbb{R}^m\to\mathbb{R}^m$ is the nonlinear Navier-Stokes operator. Performing a Taylor-series expansion of $\mathbf{q}$ about a linearisation point, $\overline{\mathbf{q}}$ in state space, which is typically the mean of $\mathbf{q}$, (\ref{eq:qeq}) can be written,
\begin{align} \label{eq:qeqlin}
\partial_t \mathbf{q}^\prime -\mathcal{A}(\mathbf{q}^\prime) = \mathbf{f},
\end{align}
where $(\cdot)^\prime$ denotes fluctuations, $\mathcal{A}\triangleq D\mathcal{N}(\overline{\mathbf{q}})$ is the Fr\'{e}chet derivative of $\mathcal{N}$ at $\overline{\mathbf{q}}$ (which becomes the Jacobian $\mathbf{A}\triangleq\partial_{q}\mathcal{N}|_{\overline{\mathbf{q}}}$ for the discretised NSE) and $\mathbf{f}$ denotes all the remaining nonlinear terms, i.e., $\mathbf{f}=\mathcal{N}(\mathbf{q})-\mathcal{A}(\mathbf{q}^\prime)$. Here, an analogy with the theory of linear time-invariant (LTI) systems can be made by considering $\mathbf{A}$, as the system matrix, and $\mathbf{q}^\prime$ as the response of this system to a forcing, $\mathbf{f}$. Resolvent analysis adopts this analogy to determine a forcing function $\mathbf{f}$ that maximizes the linear gain associated with the resolvent operator, $\mathbf{R}\triangleq(i \omega \mathbf{I} - \mathbf{A})^{-1}$.

To investigate the impact of the choice of variables, we define a nonlinear bijective transformation, $\mathcal{H}:\mathbb{R}^m\to\mathbb{R}^m$, as,
\begin{align} \label{eq:trans}
\mathbf{q}_T \triangleq \mathcal{H}(\mathbf{q})
\end{align}
that maps the state variable from the original to a new set of variables, e.g., from primitive to conservative variables. The governing equation (\ref{eq:qeq}) can be represented in terms of the transformed variable $\mathbf{q}_T$ by defining another operator, $\mathcal{N}_T$ that satisfies,
\begin{align} \label{eq:veq}
\partial_t \mathbf{q}_T = \mathcal{N}_T(\mathbf{q}_T),
\end{align}
which can also be linearized around its mean, $\overline{\mathbf{q}}_T$, yielding
\begin{align} \label{eq:veqlin}
\partial_t \mathbf{q}_T^\prime -\mathcal{A}_T(\mathbf{q}_T^\prime) = \mathbf{f}_T.
\end{align}
In what follows, we express the $D\mathcal{N}_T(\mathbf{q}_T)$ in terms of $D\mathcal{N}(\mathbf{q})$. Differentiating (\ref{eq:trans}) with respect to time and inserting (\ref{eq:qeq}) and (\ref{eq:veq}) in the result gives,
\begin{align} \label{eq:vvsq}
\mathcal{N}_T(\mathbf{q}_T)=D\mathcal{H}(\mathbf{q})\mathcal{N}(\mathbf{q}),
\end{align}
where the multiplication with a Fr\'{e}chet derivative defines an inner product. To calculate the Jacobian, we need to take the derivative of the LHS of (\ref{eq:vvsq}) with respect to $\mathbf{q}_T$, which is equivalent to taking the derivative of the RHS with respect to $\mathcal{H}(\mathbf{q})$. Given that $DF(\mathcal{H}(\mathbf{q}))=D(F\circ \mathcal{H})(\mathbf{q})\left(D \mathcal{H}(\mathbf{q})\right )^{-1}$ for a smooth $F$, the Jacobian can be written as,
\begin{align} \label{eq:vvsqfinal}
D\mathcal{N}_T(\mathbf{q}_T)= D\mathcal{H}(\mathbf{q})D\mathcal{N}(\mathbf{q})\left(D \mathcal{H}(\mathbf{q})\right )^{-1} +
 D^2\mathcal{H}(\mathbf{q}) \left(D \mathcal{H}(\mathbf{q})\right )^{-1}\mathcal{N}(\mathbf{q}).
\end{align}
The LHS of (\ref{eq:vvsqfinal}) should be calculated at $\overline{\mathbf{q}}_T$ in order to obtain $\mathcal{A}_T$. The state vectors $\overline{\mathbf{q}}$ and $\overline{\mathbf{q}}_T$ are not equivalent in the sense that they do not satisfy (\ref{eq:trans}), i.e., $\overline{\mathbf{q}}_T\neq \mathcal{H}(\overline{\mathbf{q}})$. The equivalent expansion point for the RHS of (\ref{eq:vvsqfinal}), which is not equal to $\overline{\mathbf{q}}$, is then defined as,
\begin{align} \label{eq:qtilde}
\tilde{\mathbf{q}}\triangleq\mathcal{H}^{-1}(\overline{\mathbf{q}}_T).
\end{align}
Calculation of $\overline{\mathbf{q}}_T$ requires the dynamic state data. In case only the mean flow statistics are available, a first-order prediction of $\tilde{\mathbf{q}}$ can be achieved as described in the appendix.

Having defined $\tilde{\mathbf{q}}$, (\ref{eq:vvsqfinal}) can then be re-written as,
\begin{align} \label{eq:avsau}
{\mathcal{A}}_{T}= D\mathcal{H}(\tilde{\mathbf{q}})\tilde{\mathcal{A}}\left(D \mathcal{H}(\tilde{\mathbf{q}})\right )^{-1} +
 D^2\mathcal{H}(\tilde{\mathbf{q}}) \left(D \mathcal{H}(\tilde{\mathbf{q}})\right )^{-1}\mathcal{N}(\tilde{\mathbf{q}}),
\end{align}
where $\tilde{\mathcal{A}}\triangleq D\mathcal{N}(\tilde{\mathbf{q}})$. Note that the second term on the RHS of (\ref{eq:vvsqfinal}) is zero if $\mathcal{H}$ is linear ($D^2\mathcal{H}(\mathbf{q})=0$ for any $\mathbf{q}$) or if linearisation is performed about a `fixed' point of the system, which, by definition, satisfies $\mathcal{N}(\overline{\mathbf{q}})=0$. In that case, the first term in (\ref{eq:avsau}) indicates a similarity transformation between ${\mathcal{A}}_T$ and $\tilde{\mathcal{A}}$, which now reduces to $\mathcal{A}$ (due to having $\overline{\mathbf{q}}_T= \mathcal{H}(\overline{\mathbf{q}})$ satisfied, by definition in case of linear $\mathcal{H}$, and through (\ref{eq:vvsq}) in case of expansion around the fixed point), and the two linear operators share the same eigenvalues. In mean-flow-based linear analysis, however, the second term on the RHS (\ref{eq:avsau}) is non-zero, which implies a difference in the eigenvalues of $\tilde{\mathcal{A}}$ and ${\mathcal{A}}_T$ even if the corresponding linearisation points are equivalent.

We now derive an expression relating $\mathcal{A}$ and $\mathcal{A}_T$. It is more convenient to compare these two operators since one would use $\mathcal{A}$, not $\tilde{\mathcal{A}}$ when performing mean-flow linear analysis. Although (\ref{eq:avsau}) provides a means to compute $\mathcal{A}_T$, which can then be compared to $\mathcal{A}$, an exact general expression relating the two operators cannot be derived for nonlinear $\mathcal{H}$. However an approximate relation is given by replacing (\ref{eq:atilde}) into (\ref{eq:avsau}) as,

\begin{align} \label{eq:avsauapr}
{\mathcal{A}}_{T}&\approx D\mathcal{H}(\tilde{\mathbf{q}}){\mathcal{A}}\left(D\mathcal{H}(\tilde{\mathbf{q}})\right )^{-1} \nonumber \\
&+ D\mathcal{H}(\tilde{\mathbf{q}}) \big(\mathcal{A}D\epsilon(\overline{\mathbf{q}}) + D^2\mathcal{N}(\overline{\mathbf{q}})\epsilon(\overline{\mathbf{q}})\big)\left(D\mathcal{H}(\tilde{\mathbf{q}})\right )^{-1}
+ D^2\mathcal{H}(\tilde{\mathbf{q}}) \left(D\mathcal{H}(\tilde{\mathbf{q}})\right )^{-1}\mathcal{N}(\tilde{\mathbf{q}}).
\end{align}
The second term on the RHS of \eqref{eq:avsauapr} accounts for account the change in the linearisation point, i.e., non-equivalence of the mean flows, $\overline{\mathbf{q}}$ and $\overline{\mathbf{q}}_T$ due to taking the mean after nonlinear transformation, while the third term appears once again since the mean flow is not a fixed point of the Navier-Stokes system for turbulent flows. It is obvious that these extra terms do not cancel each other for arbitrary $\mathcal{H}$. As a result, expanded about the mean-flow, the system characteristics are modified when subjected to a nonlinear transformation. An immediate implication is that the stability characteristics of the linear operator are not unique and the linearised systems (\ref{eq:qeqlin}) and (\ref{eq:veqlin}) are equivalent only if the right-hand sides are maintained and the equal sign is respected.
\subsection{Implications for resolvent analysis} \label{subsec:res}
Resolvent analysis involves taking the Fourier transform of (\ref{eq:qeqlin}), or similarly (\ref{eq:veqlin}), and re-organizing the result as,
\begin{align} \label{eq:res}
\hat{\mathbf{q}}=\mathbf{R}\hat{\mathbf{f}},
\end{align}
where $\hat{\mathbf{q}}$ and $\hat{\mathbf{f}}$ are Fourier transforms of $\mathbf{q}$ and $\mathbf{f}$, respectively. Ignoring the nonlinear relation between $\hat{\mathbf{f}}$ and $\hat{\mathbf{q}}$, (\ref{eq:res}) can be seen as a forcing-response relation. The optimal forcing that would maximize the energy in the response can be be found by maximising the Rayleigh quotient, $\sigma \triangleq {\langle \hat{\mathbf{q}}, \hat{\mathbf{q}} \rangle}/{\langle \hat{\mathbf{f}}, \hat{\mathbf{f}}\rangle}$, where the inner product is defined as
\begin{align} \label{eq:inner}
\langle \mathbf{a},\mathbf{b}\rangle=\int_V\mathbf{a}^*\mathbf{M}\mathbf{b}dV,
\end{align}
with $\mathbf{M}$ defining a suitable energy norm. The Rayleigh quotient for the transformed system, $\sigma_T$, can be similarly defined using $\hat{\mathbf{q}}_T$ and $\hat{\mathbf{f}}_T$. To compare $\sigma$ and $\sigma_T$, the energy norms in the original and the transformed systems should be equivalent, which amounts to
\begin{equation} \label{eq:normequiv}
\int_V {\mathbf{q}^\prime}^*\mathbf{M}\mathbf{q}^\prime dV=\int_V {\mathbf{q}^\prime}^*_T\mathbf{M}_T\mathbf{q}^\prime_TdV.
\end{equation}
To leading order in fluctuation amplitude, $\mathbf{q}^\prime_T$ and $\mathbf{q}^\prime$ are related by the expression,
\begin{equation} \label{eq:fluc}
\mathbf{q}^\prime_T=\partial_q\mathcal{H}|_{\overline{q}}\mathbf{q}^\prime.
\end{equation}
Replacing $\mathbf{q}^\prime_T$ terms in (\ref{eq:normequiv}) using (\ref{eq:fluc}) yields a leading-order norm equivalence,
\begin{equation} \label{eq:normcorr}
\mathbf{M}=\partial_q\mathcal{H}|_{\overline{q}}^*\mathbf{M}_T\partial_q\mathcal{H}|_{\overline{q}}.
\end{equation}
Note that the norm correction is required even when the variable transformation is linear, in which case the resolvent operators are connected through a similarity transform. Equation (\ref{eq:normcorr}) can alternatively be obtained by enforcing that the Rayleigh quotients of the original and linearly transformed systems be equal.
\section{Application to test cases} \label{sec:testcases}
The above analysis is tested using four different LES databases of ideally expanded round jets: three isothermal jets at $M_j=0.4$, 0.9 and 1.5 and a hot jet at $M_j=1.5$ are investigated to observe the effect of switching from conservative to primitive variable sets on resolvent analysis. The operating conditions are given in terms of jet exit Mach number, $M_j$, nozzle pressure ratio, $NPR\triangleq P_{0,j}/P_\infty$, nozzle temperature ratio, $NTR\triangleq T_{0,j}/T_\infty$, and Reynolds number, $Re\triangleq\rho_jU_jD/\mu_j$, where $D$ is the nozzle diameter, $\mu_j$ is the dynamic viscosity at the jet exit, and the subscript 0 is used to denote stagnation quantities. Further details about the LES databases can be found in \cite{bres_aiaa_2017,bres_jfm_2018}.

\begin{table}
\vspace{-10px}
\caption{The list of axisymmetric jets investigated.}
\centering
\label{tab:2}
\begin{tabular}{p{2cm} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}}}
& $M_j$ & $T_j/T_\infty$ & $NPR$ & $NTR$ & $Re$ \rule{0pt}{4ex}\\
 \emph{Jet-1} & 0.4 & 1.0 & 1.117 & 1.032 & $0.45\times 10^6$ \rule{0pt}{3ex}\\
 \emph{Jet-2} & 0.9 & 1.0 & 1.691 & 1.162 & $1.01\times 10^6$\\
 \emph{Jet-3} & 1.5 & 1.0 & 3.671 & 1.450 & $1.69\times 10^6$\\
 \emph{Jet-4} & 1.5 & 1.79 & 3.671 & 2.596 & $0.95\times 10^6$\\
\end{tabular}
\end{table}

For each flow case, we performed resolvent analysis and spectral proper orthogonal decomposition (SPOD) of the axisymmetric mode using two different variable sets: $[\nu\; u_x\; u_r\; p]^T$ and $[\rho\; \rho u_x\; \rho u_r\; \rho E]^T$, which we refer to as primitive and conservative, respectively. Here, $\nu=1/\rho$, $\rho$, $u_x$, $u_r$, $p$ and $\rho E$ denote the specific volume, density, axial and radial velocities, pressure and total energy, respectively. In order to perform a singular value decomposition (SVD) of the resolvent, the Jacobian matrix of the compressible NSE, i.e., the linear operator $\mathbf{A}$ defined in \S\ref{subsec:derivation}, is computed using conservative variables by linearising the discretised equation as described by \cite{mettot2014computation}. At this point, it is possible to switch to primitive variables by first finding the linearisation point that is equivalent to the mean in primitive variables using equation \eqref{eq:qtilde}, and then applying equation \eqref{eq:avsau} to get $\mathbf{A}_T$. The SVD of the resolvent is then calculated based on Krylov methods as detailed in \cite{bugeat20193d}.
The open libraries PETSc \citep{petsc-efficient} and SLEPc \citep{Hernandez2005SSF} are used to solve the linear systems by direct LU decomposition and eigenvalue calculation by the Krylov-Schur algorithm \citep{hernandez2007krylov}.
Computations are carried out using an orthogonal mesh of $748\times 229$ points for the subsonic cases, and $636\times 229$ points for the supersonic cases, which were determined after a convergence study, in the numerical domain $x/D \in [0;30]$ and $r/D \in [0;12]$, where $x/D=0$ is the location of the nozzle exit and $r/D=0$ is the centre line.
Sponge zones are located at $x/D>20$ and $r/D>5$.
The Chu norm \citep{chu1965energy} is used for both response and forcing over $x/D \in [0;20]$ and $r/D \in [0;5]$.

To calculate the SPOD modes, we followed the method introduced by \cite{towne_jfm_2018}. Once again, the Chu norm is used. We calculated the Fast Fourier Transform (FFT) of $\mathbf{q}$ with 128 FFT points, and overlap ratio of 0.75, yielding 310 FFT blocks for the isothermal cases, and 154 blocks for the hot jet case.

\begin{figure}
  \centerline{\resizebox{1\textwidth}{!}{\includegraphics{fig1.eps}}
  \caption{Comparison of the leading SPOD and optimal resolvent modes at $St=0.6$ for different jets (\emph{Jet-1} to \emph{Jet-4}: top to bottom). Pressure fields are shown. For the analyses using conservative variables (right), pressure is obtained through linearised transformation. }
\label{fig:optmodes}
\vspace{-15px}
\end{figure}

The optimum response modes at $St=0.6$ obtained using primitive and conservative variables, respectively, are plotted in comparison to the optimum SPOD modes in primitive variables in Figure \ref{fig:optmodes} for all jet cases. The figure shows the pressure field, which is directly computed for the primitive-variable cases, while it is reconstructed using (\ref{eq:fluc}) for the conservative-variable cases. The reconstructed SPOD modes are identical to the directly computed ones (see Table \ref{tab:5}), so are not shown here. The mode shapes are very similar between the directly computed and the reconstructed response modes. However, the global weight distribution among the variables slightly differs in \emph{Jet-4} (higher for pressure in the reconstructed mode).

The transformation from primitive to conservative variables, which are denoted by the subscripts $p$ and $c$, respectively, is defined as
$\mathbf{q}_c=\mathcal{H}(\mathbf{q}_p),$ with
\begin{align} \label{eq:transfun}
\mathcal{H}(\mathbf{q})=\left[\frac{1}{\mathbf{q}_1},\;\frac{\mathbf{q}_2}{\mathbf{q}_1},\;\frac{\mathbf{q}_3}{\mathbf{q}_1},\; \frac{\mathbf{q}_4}{(\gamma-1)}+\frac{1}{2}\frac{1}{\mathbf{q}_1}\left(\mathbf{q}_2^2 + \mathbf{q}_3^2 \right)\right]^T,
\end{align}
where $\gamma=1.4$ is the heat capacity ratio. This function becomes linear for density and momentum terms when the flow is incompressible. The jet cases are listed in the order of increasing variable-density effects. Therefore we expect larger deviation in the singular values of the original and transformed resolvent operators as we move from \emph{Jet-1} to \emph{Jet-4}. To quantify the effect of density fluctuations on the mean flow, we define the percentage measure $\Delta\overline{\mathbf{q}}\triangleq\left|\mathcal{H}^{-1}(\overline{\mathbf{q}}_c) - \overline{\mathbf{q}}_p\right|/\overline{\mathbf{q}}_p\times100$. We integrate this quantity over the jet domain where turbulent kinetic energy (TKE) is greater than 1\% of its maximum value and normalize the result with the integration domain. To quantify the change in the singular values, we define, $\Delta\sigma_1 \triangleq (\sigma_{1,c}-\sigma_{1,p})/\sigma_{1,p}\times100$, and $\Delta G \triangleq (G_c-G_p)/G_p\times100$, where $G=\sigma_1/\sigma_2$ is the gain separation between the optimal and the first suboptimal singular values. The modification of the mean and the resulting changes in the singular values are tabulated in Table \ref{tab:3}. It can be seen from the table that the difference between the true mean, $\overline{\mathbf{q}}_p$ and the transformed mean, $\mathcal{H}^{-1}(\overline{\mathbf{q}}_c)$ does not exceed 0.6\% except in the radial velocity field where 7.6\% overall difference is seen in \emph{Jet-4}. The radial profiles of the mean flow calculated directly in primitive variables, $\overline{\mathbf{q}}_p$, and reconstructed from conservative variables, $\tilde{\mathbf{q}}_p$ are shown for \emph{Jet-4} in Figure \ref{fig:meanx2}, respectively. Once again, the difference between $\overline{\mathbf{q}}_p$ and $\tilde{\mathbf{q}}_p$ is mostly visible in radial velocity and pressure around the shear layer, with the largest percentage change in the former. This difference indicates a relatively strong correlation between density and radial velocity, and similarly between density and pressure, in this region. Despite the small modification of the mean, the optimal singular value and gain separation are modified by up to 40\% and 35\%, respectively, in \emph{Jet-4}.

\begin{table}
\vspace{-10px}
\caption{The change in the mean fields and the singular values calculated at $St=0.6$ due to transformation from primitive to conservative variables.}
\centering
\label{tab:3}
\begin{tabular}{p{2cm} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} }
& $\Delta\overline{\nu}$ [\%]& $\Delta\overline{u_x}$ [\%]& $\Delta\overline{u_r}$ [\%]& $\Delta\overline{p}$ [\%]& $\Delta\sigma_1$ [\%]& $\Delta G$ [\%] \rule{0pt}{4ex}\\
 \emph{Jet-1} & 0.00 & 0.01 & 0.31 & 0.04 & 0.05 & $-0.05$ \rule{0pt}{3ex}\\
 \emph{Jet-2} & 0.02 & 0.02 & 1.21 & 0.18 & $-2.81$ & $-3.25$ \\
 \emph{Jet-3} & 0.05 & 0.06 & 1.43 & 0.57 & $-6.02$ & $-7.45$ \\
 \emph{Jet-4} & 0.24 & 0.16 & 7.65 & 0.61 & $-40.58$ & $-35.87$ \\
\end{tabular}
\end{table}

\begin{figure}
  \centerline{\resizebox{1\textwidth}{!}{\includegraphics{fig2.eps}}
  \caption{Radial mean flow profiles calculated in primitive variables (solid) and reconstructed from conservative variables (dashed) at $x/D=2$ in \emph{Jet-4}. }
\label{fig:meanx2}
\vspace{-5px}
\end{figure}

For the two supersonic cases, \emph{Jet-3} and \emph{Jet-4} in which noticeable changes in the singular values are observed, the analysis is repeated at different frequencies. The optimal gain spectra obtained using primitive variables and their modification under variable transformation are plotted in Figure \ref{fig:spec}. It is seen that maximum change in singular values is obtained at $St=0.4$ and $St=0.6$ for \emph{Jet-3} and \emph{Jet-4}, respectively.

\begin{figure}
  \centerline{\resizebox{1\textwidth}{!}{\includegraphics{fig3.eps}}
  \caption{Optimal gain spectra computed using primitive variables (left), and their modification under variable transformation (right). }
\label{fig:spec}
\vspace{-10px}
\end{figure}

We furthermore investigate the alignment of the optimum response mode of the resolvent operator with the leading SPOD mode of the response, $\mathbf{q}$. The alignment is quantified as the absolute value of the scalar product, $|\langle\mathbf{\psi}_S,\mathbf{\psi}_R\rangle|$, where the subscripts $R$ and $S$ stand for the resolvent optimal response and SPOD modes, respectively. The results at $St=0.6$ are tabulated in Table \ref{tab:4} for all jets. It is seen that the alignment is slightly better when primitive variables are used, except for \emph{Jet-4}, where the alignment is improved by $\sim$24\% when conservative variables are used, although the gain separation is lowered by 35\%. This is counter-intuitive since the alignment of SPOD and resolvent modes is often associated with high gain separation in the resolvent operator \citep{beneddine_jfm_2016}, considering that the non-linear terms are close to white noise and thus have similar projection coefficients onto the forcing mode basis. However, the presence of coherent, low-rank forcing, as recently observed for turbulent channel flow \citep{morra_arxiv_2020}, may change this picture. The overall weight of a response mode in $\hat{\mathbf{q}}$ is determined by the corresponding gain together with the projection of the actual forcing, $\hat{\mathbf{f}}$ on to the corresponding forcing mode of the resolvent operator. Therefore, the observed increase in the alignment despite of the reduction in gain separation indicates enhanced projection of forcing onto the optimal forcing mode, which is sufficient to outweigh the reduced gain separation.

\begin{table}
\vspace{-10px}
\caption{Alignment coefficients between the leading SPOD and response modes for different jet cases at $St=0.6$. Results obtained using primitive and conservative variables, respectively, are compared.}
\centering
\label{tab:4}
\begin{tabular}{p{4.5cm} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} }
& \emph{Jet-1} & \emph{Jet-2} & \emph{Jet-3} & \emph{Jet-4} \rule{0pt}{4ex}\\
 Primitive variables & 0.593 & 0.571 & 0.578 & 0.420 \rule{0pt}{3ex}\\
 Conservative variables & 0.587 & 0.542 & 0.555 & 0.521  \\
\end{tabular}
\rule{\textwidth}{.25pt}
\caption{Alignment coefficients between the directly computed and the reconstructed modes for different jet cases at $St=0.6$.}
\centering
\label{tab:5}
\begin{tabular}{p{4.5cm} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} c@{\hspace{.5cm}} }
& \emph{Jet-1} & \emph{Jet-2} & \emph{Jet-3} & \emph{Jet-4} \rule{0pt}{4ex}\\
 Leading SPOD mode & 1.000 & 1.000 & 0.999 & 0.981 \rule{0pt}{3ex}\\
 Optimal resolvent mode & 0.996 & 0.987 & 0.983 & 0.462  \\
\end{tabular}
\end{table}

Finally, we investigate the change in SPOD and resolvent modes that occurs under the variable transformation. We transform the modes obtained using conservative variables to their primitive-variable counterparts using (\ref{eq:fluc}). The inner-product definition given in (\ref{eq:inner}) is now used to measure the alignment between the reconstructed and directly calculated modes. The results at $St=0.6$ are tabulated in Table \ref{tab:5}. The alignment of the leading SPOD modes calculated in conservative and primitive variables is above 98\% for all the jets. The similarity of directly computed and reconstructed optimal resolvent modes are again above 98\% for all the cases except for the heated, compressible \emph{Jet-4}, where we see a drop to 46\%. This change, consistent with the significant change observed in the singular vales, shows how the choice of variables may lead to quite different results and subsequent conclusions about the underlying flow physics.
\section{Discussion and concluding remarks} \label{sec:conc}
We have demonstrated a non-uniqueness issue associated with mean-flow-based linear analysis. We show that the characteristics of the linear operator that is obtained by linearising a non-linear dynamical system about its mean depend on the state variables considered. In the framework of resolvent analysis, this means that two studies of a given flow, using non-linearly related dependent variables (conservative and primitive, for instance), may produce different results.

This raises questions regarding the interpretation of such analyses, now widely used to study coherent structures in turbulent shear flow. The application of such analysis for jets is often justified on the basis of their weak non-parallelism, which may result in a large gain separation between leading and sub-optimal resolvent modes, and an associated spatial separation between forcing and response modes.
By means of resolvent analyses of turbulent jets with different operating conditions, we demonstrate how primitive- and conservative-variable-based analyses may differ significantly (up to 40\% change in gain; 35\% change in gain separation). This implies that mean-flow-based linear analyses of amplifier flows, where  flow dynamics are strongly dependent on forcing, cannot be regarded as universal, but are instead dependent on the specific form (and model) considered for the forcing terms

A similar issue has arisen in the context of aeroacoustics. Acoustic analogies of different forms have been proposed over the years (\citet{lighthill_prs_1952,lilley1974noise,howe_jfm_1975,doak_ap_1995,goldstein_2003}; etc.) based on different linear operators considered to describe sound propagation. For each wave operator there exists a corresponding `source' term; and it is this that confounds attempts to uniquely define what is meant by a `source' of sound in turbulent flows (cf. \cite{jordan2008sja}). We see that the same situation holds for mean-flow-based stability analysis. Just as there is no unique acoustic wave propagator, there is no unique resolvent operator. However, as these approaches are all exact rearrangements of the governing equations, with residual non-linear terms  treated as external `force' or `source', all analyses will lead to the same result if the residual terms are retained. This highlights the importance of not neglecting, or over-simplifying, the external forcing term in resolvent analyis.

While this ambiguity may be unsettling, it opens the door to optimisation of resolvent-based approaches for various applications. How to optimise the linear framework could vary depending on the problem considered. For instance, for a supersonic jet where a rank-1 model is sufficient to describe peak jet-noise \citep{sinha_jfm_2014,cavalieri_amr_2019}, finding a nonlinear transformation that maximises gain separation, whilst keeping the nonlinear forcing maximally aligned with the leading input mode, would constitute an interesting optimisation problem. For a subsonic jet, on the other hand, multiple input-output modes are necessary to correctly describe sound generation \citep{cavalieri_jfm_2014,towne_aiaa_2015,cavalieri_amr_2019}. Finding a transformation that maximises the gains and projections onto the forcing modes for the first $n$ in that case, where $n\geq2$ is to be determined as well, may help improve the modelling strategies.\\

This work has received funding from the Clean Sky 2 Joint Undertaking (JU) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 785303. Results reflect only the authors' view and the JU is not responsible for any use that may be made of the information it contains.

\vspace{-0.42cm}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Deep neural networks have been highly successful in
recent years \cite{graves2013generating,he2016identity,silver2016mastering,vaswani2017attention,zaremba2014recurrent}.
To achieve state-of-the-art performance, they often have to
leverage the computing power of multiple
machines during training \cite{dean-12,xing2015petuum,zinkevich2010parallelized,chen2016revisiting}.
Popular approaches include distributed synchronous SGD and its momentum variant
SGDM,
in
which the computational load for evaluating a mini-batch gradient is distributed
among the workers.
Each worker performs local computation,
and these local informations are then merged by the server for final update on the model parameters.
However, its scalability
is limited by the
possibly
overwhelming
cost due to
communication of the gradient and model parameter \cite{li2014communication}.
Let $d$ be the
gradient/parameter
dimensionality, and $M$ be the number of workers.
$64Md$ bits need to be transferred between the workers and server
in each iteration.

To mitigate this communication bottleneck,
the two common approaches are gradient sparsification
and gradient quantization.
Gradient sparsification only sends the most significant,
information-preserving
gradient entries.
A heuristic algorithm is first
introduced in
\cite{seide20141},
in which only the large
entries
are transmitted.
On training a neural machine translation model
with 4 GPUs,
this greatly reduces the communication overhead and achieves 22\% speedup
\cite{aji2017sparse}.
Deep gradient compression \cite{lin2017deep} is another heuristic method that combines gradient sparsification with other techniques such as momentum correction,
local gradient clipping, and momentum factor masking, achieving significant
reduction on communication cost.
Recently, a stochastic sparsification method was proposed in \cite{wangni2018gradient} that
balances sparsity and variance
by solving a constrained linear programming.
MEM-SGD \cite{stich2018sparsified}
combines top-$k$ sparsification with error correction.
By keeping track of the accumulated
errors, these can be added back to the gradient estimator before each transmission.
MEM-SGD converges at the same rate as SGD on convex problems, whilst
reducing the communication overhead by a factor equal to the problem
dimensionality.

On the other hand,
gradient quantization mitigates the communication bottleneck
by lowering the gradient's floating-point precision with a smaller bit width.
1-bit SGD achieves state-of-the-art results
on acoustic modeling
while dramatically reducing the communication cost
\cite{seide20141,strom2015scalable}.
TernGrad \cite{wen2017terngrad}
quantizes the gradients to ternary levels $\{-1, 0, 1\}$.
QSGD \cite{alistarh2017qsgd} employs stochastic randomized rounding
to ensure unbiasedness of the estimator.
Error-compensated quantized SGD (ECQ-SGD) was proposed in \cite{wu2018error}, wherein a similar stochastic quantization function used in
QSGD is employed, and an error bound is obtained for quadratic loss functions.
Different from the error-feedback mechanism proposed in MEM-SGD, ECQ-SGD requires
two more hyper-parameters and its quantization errors are decayed exponentially. Thus, error feedback is limited to a small number of iterations.
Also, ECQ-SGD uses
all-to-all broadcast (which may involve large network traffic and idle time), while we consider parameter-server architecture.
Recently, Bernstein {\em et al.\/}
proposed
signSGD with majority vote \cite{pmlr-v80-bernstein18a},
which only
transmits the
1-bit
gradient sign
between workers and server.
A variant
using momentum, called
signum with majority vote,
is also introduced though without convergence analysis \cite{bernstein2019signsgd} .
Using the majority vote,
signSGD
achieves a notion of Byzantine fault tolerance \cite{bernstein2019signsgd}.
Moreover, it
converges at the same rate as
distributed SGD,
though it has to rely on the unrealistic
assumptions of having a
large mini-batch and
unimodal symmetric
gradient noise.
Indeed,
signSGD can diverge in some simple cases  when these assumptions are violated
\cite{karimireddy2019error}.
With only a single worker,
this
divergence issue
can be fixed by
using the error correction technique in MEM-SGD, leading to
SGD with error-feedback (EF-SGD)
\cite{karimireddy2019error}.

While only a single worker is considered in EF-SGD,
we study in this paper the more interesting distributed setting.
An extension of MEM-SGD and EF-SGD with parallel computing was proposed in
\cite{cordonnier2018convex} for all-to-all broadcast.
Another related architecture is allreduce.
Compression at the server can be implemented between the reduce and broadcast
steps in tree allreduce, or between the reduce-scatter and allgather steps in ring allreduce.
However, allreduce requires repeated gradient aggregations, and the compressed gradients need to be first
decompressed before they are summed. Hence, heavy overheads may be incurred.

In this paper,
we study
the distributed setting with a parameter server architecture.
To ensure
efficient
communication,
we consider two-way
gradient compression, in which gradients in both directions (server to/from workers) are compressed.  Note that existing works (except signSGD/signum with majority vote \cite{pmlr-v80-bernstein18a,bernstein2019signsgd}) do not compress the aggregated gradients before sending
back to
workers.
Moreover,
as gradients in a deep network typically have similar magnitudes in each
layer,
each layer-wise gradient can be sufficiently represented using a sign vector and
its average $\ell_1$-norm.
This layer-wise (or blockwise in general) compressor
achieves nearly $32$x reduction in communication cost.
The resulant procedure is called communication-efficient distributed SGD with
error-feedback
(dist-EF-SGD).
Analogous to SGDM, we also propose a stochastic variant
dist-EF-SGDM
with Nesterov's momentum
\cite{nesterov1983method}.
The convergence properties of
dist-EF-SGD(M)
are studied theoretically.

Our contributions are: (i)
We provide a bound
on dist-EF-SGD
with general stepsize schedule for a class of compressors
(including the commonly used sign-operator and top-$k$ sparsification).
In particular, without relying on the unrealistic assumptions in
\cite{pmlr-v80-bernstein18a,bernstein2019signsgd}, we show that dist-EF-SGD
with constant/decreasing/increasing stepsize
converges at an
$\mathcal{O}(1/\sqrt{MT})$
rate, which matches that of
distributed synchronous SGD;
(ii) We study gradient compression with Nesterov's momentum in a parameter server.
For dist-EF-SGDM
with constant stepsize,
we obtain an
$\mathcal{O}(1/\sqrt{MT})$ rate.
To the best of our knowledge, these are the first convergence results on two-way gradient compression with Nesterov's momentum;
(iii) We propose a general blockwise compressor and show its theoretical properties.
Experimental results show that
the proposed algorithms
are efficient
without losing prediction accuracy.
After our paper has appeared,
we note
a similar idea was independently proposed in \cite{tang2019doublesqueeze}. Different from ours, they do not consider changing stepsize, blockwise compressor and Nesterov's momentum.

{\bf Notations}. For a vector $x$, $\|x\|_1$ and $\|x\|_2$ are its $\ell_1$- and $\ell_2$-norms, respectively.
$\text{sign}(x)$ outputs a vector in which each element is the sign of the corresponding entry of $x$.
For two vectors $x, y$, $\langle x, y \rangle$ denotes the dot product.
For a function $f$, its gradient is $\nabla f$.

\section{Related Work:
SGD with Error-Feedback}

In machine learning, one is often interested in minimizing the expected
risk
$F(x) = {\mathbb E}_{\xi}[f(x, \xi)]$.
which directly measures the generalization error \cite{bottou2004large}.
Here,
$x \in {\mathbb R}^d$ is the model parameter,
$\xi$ is drawn from some unknown distribution,
and
$f(x, \xi)$ is the
possibly nonconvex
risk due to $x$.
When the expectation is taken over a training set of size $n$, the expected risk
reduces to empirical risk.

Recently, Karimireddy {\em et al.\/}
\cite{karimireddy2019error} introduced
SGD with error-feedback (EF-SGD),
which combines gradient compression with error correction
(Algorithm~\ref{alg:ef-sgd}).
A single machine
is considered,
which keeps the gradient difference
that is not used for parameter update in the current iteration.
In the next iteration $t$, the accumulated residual
$e_t$
is added to the current gradient. The corrected gradient
$p_t$
is then fed into
an $\delta$-approximate compressor.

\begin{definition} \label{definition:ef_compressor}
\cite{karimireddy2019error}
An operator $\mathcal{C}: {\mathbb R}^d \rightarrow {\mathbb R}^d$ is an $\delta$-approximate compressor for $\delta \in (0, 1]$ if
$\|\mathcal{C}(x) - x\|^2_2 \leq (1 - \delta)\|x\|_2^2$.
\end{definition}
Examples
of $\delta$-approximate compressors
include the scaled sign operator
$\mathcal{C}(v) = \|v\|_1/d\cdot\text{sign}(v)$ \cite{karimireddy2019error}
and top-$k$ operator (which only preserves the $k$ coordinates with the largest absolute values) \cite{stich2018sparsified}.
One can also have randomized compressors that only satisfy
Definition~\ref{definition:ef_compressor}
in expectation.
Obviously, it is desirable to have
a large
$\delta$
while achieving low communication cost.

\begin{algorithm}[ht]
\caption{SGD with Error-Feedback (EF-SGD) \cite{karimireddy2019error}}
\label{alg:ef-sgd}
\begin{algorithmic}[1]
    \STATE {\bfseries Input:} stepsize $\eta$; compressor $\mathcal{C}(\cdot)$.
    \STATE {\bfseries Initialize:} $x_0 \in \mathbb{R}^d$; $e_{0} = 0 \in \mathbb{R}^d$
    \FOR{$t = 0, \ldots, T-1$}

        \STATE $p_{t} = \eta g_{t} +  e_{t}$ \algorithmiccomment{stochastic gradient $g_{t} = \nabla f(x_t, \xi_t)$}
        \STATE $\Delta_{t} = \mathcal{C}(p_{t})$
		  \COMMENT{compressed value output}
        \STATE $x_{t+1} = x_t - \Delta_{t}$
        \STATE $e_{t+1} = p_{t} - \Delta_{t}$
    \ENDFOR
\end{algorithmic}
\end{algorithm}

EF-SGD achieves the same $\mathcal{O}(1/\sqrt{T}))$ rate as SGD.
To obtain this convergence guarantee,
an important observation
is that the error-corrected iterate ${\tilde{x}}_t = x_t - e_t$ satisfies the recurrence:
${\tilde{x}}_{t+1} = {\tilde{x}}_t -  \eta g_t$,
which is similar to that of SGD.
This allows utilizing the convergence proof of SGD to bound the gradient difference $\|\nabla F({\tilde{x}}_t) - \nabla F(x_t)\|_2$.

\section{Distributed Blockwise Momentum SGD with Error-Feedback}

\subsection{Distributed SGD with Error-Feedback}
\label{sec:dist-ef-sgd}

The proposed procedure,
which extends EF-SGD to the distributed setting.
is shown in Algorithm~\ref{alg:dist-ef-sgd}.
The computational workload
is distributed
over $M$ workers. A local accumulated error vector $e_{t, i}$ and a local corrected gradient vector $p_{t,i}$
are stored in the memory of worker $i$.
At iteration $t$, worker $i$ pushes the compressed signal $\Delta_{t,i} = \mathcal{C}(p_{t, i})$ to the parameter server.
On the server side, all workers' $\Delta_{t,i}$'s are aggregated and used to update
its global error-corrected vector ${\tilde{p}}_t$. Before sending back the final update
direction ${\tilde{p}}_{t}$ to each worker,
compression
is performed
to ensure a comparable amount of communication costs between the push and pull operations.
Due to gradient compression on the server, we also employ a global accumulated error vector $\tilde{e}_t$.
Unlike EF-SGD in Algorithm~\ref{alg:ef-sgd}, we do not multiply gradient $g_{t, i}$
by the stepsize $\eta_t$ before compression. The two cases make no difference when
$\eta_t$ is constant. However, when the stepsize is changing over time, this would
affect convergence.
We also rescale the local accumulated error $e_{t, i}$ by $\eta_{t-1}/\eta_t$.
This modification,
together with the use of error correction on both workers and server,
allows us to obtain Lemma~\ref{lemma:recurrence_ec_iterate}.
Because of these differences,
note that dist-EF-SGD does not reduce to EF-SGD when $M=1$. When $\mathcal{C}(\cdot)$ is
the identity mapping,
dist-EF-SGD reduces to full-precision distributed SGD.

\begin{algorithm}[ht]
\caption{Distributed SGD with Error-Feedback (dist-EF-SGD)}
\label{alg:dist-ef-sgd}
\begin{algorithmic}[1]
    \STATE {\bfseries Input:}  stepsize sequence $\{\eta_t\}$ with $\eta_{-1} = 0$; number of workers $M$; compressor $\mathcal{C}(\cdot)$.
    \STATE {\bfseries Initialize:} $x_0 \in \mathbb{R}^d$; $e_{0,i} = 0 \in
	 \mathbb{R}^d$ on each worker $i$; $\tilde{e}_{0} = 0 \in \mathbb{R}^d$ on
	 server
    \FOR{$t = 0, \ldots, T-1$}
        \STATE \textbf{on each worker $i$}

        \STATE \qquad $p_{t,i} = g_{t,i} + \frac{\eta_{t-1}}{\eta_t}e_{t,i}$ \algorithmiccomment{stochastic gradient $g_{t,i} = \nabla f(x_t, \xi_{t, i})$}
        \STATE \qquad \textbf{push} $\Delta_{t,i} = \mathcal{C}(p_{t,i})$ \textbf{to server}
        \STATE \qquad $x_{t+1} = x_t - \eta_t\tilde{\Delta}_{t}$ \algorithmiccomment{$\tilde{\Delta}_{t}$ \textbf{is pulled from server}}
        \STATE \qquad $e_{t+1,i} = p_{t,i} - \Delta_{t,i}$

        \STATE \textbf{on server}
        \STATE \qquad \textbf{pull $\Delta_{t,i}$ from each worker $i$} and ${\tilde{p}}_{t} = \frac{1}{M} \sum_{i=1}^M \Delta_{t,i} + \frac{\eta_{t-1}}{\eta_t}{\tilde{e}}_{t}$
        \STATE \qquad \textbf{push $\tilde{\Delta}_t = \mathcal{C}({\tilde{p}}_{t})$ to each worker}
        \STATE \qquad $\tilde{e}_{t+1} = {\tilde{p}}_t - {\tilde{\Delta}}_t$
    \ENDFOR
\end{algorithmic}
\end{algorithm}

In the following,
we investigate the convergence of dist-EF-SGD.
We make the following assumptions, which are common in the stochastic approximation
literature.

\begin{assumption} \label{assumption:ef_smooth}
$F$ is lower-bounded (i.e., $F_* = \inf_{x \in {\mathbb R}^d} F(x) > -\infty$) and
$L$-smooth (i.e.,
$F(x) \leq F(y) + \langle \nabla F(y), x -y \rangle + \frac{L}{2}\|x - y\|_2^2$
for $x, y \in {\mathbb R}^d$).
\end{assumption}

\begin{assumption} \label{assumption:ef_bounded_var}
The stochastic gradient $g_{t, i}$ has bounded variance: ${\mathbb E}_t\left[\|g_{t, i} - \nabla F(x_t)\|_2^2\right] \leq \sigma^2$.
\end{assumption}

\begin{assumption} \label{assumption:ef_bounded_grad}
The full gradient $\nabla F$ is uniformly bounded: $\|\nabla F(x_t)\|_2^2 \leq \omega^2$.
\end{assumption}
This implies the
second moment
is bounded, i.e., ${\mathbb E}_t\left[\|g_{t, i}\|_2^2\right] \leq G^2 \equiv \sigma^2 + \omega^2$.

\begin{lemma} \label{lemma:recurrence_ec_iterate}
Consider the error-corrected iterate
${\tilde{x}}_t = x_t - \eta_{t-1}\left({\tilde{e}}_t + \frac{1}{M} \sum_{i=1}^M e_{t,i}\right)$,
where $x_t$, ${\tilde{e}}_t$, and $e_{t,i}$'s are generated from Algorithm~\ref{alg:dist-ef-sgd}. It satisfies the recurrence:
${\tilde{x}}_{t+1} = {\tilde{x}}_t -  \eta_t \frac{1}{M} \sum_{i=1}^Mg_{t,i}$.
\end{lemma}
The above Lemma shows that ${\tilde{x}}_t$ is very similar to the distributed SGD iterate except that the stochastic gradients are evaluated at $x_t$ instead of ${\tilde{x}}_t$.
This connection allows us to utilize the analysis of full-precision distributed SGD. In particular, we have the following Lemma.
\begin{lemma} \label{lemma:wm_error_bound}
${\mathbb E}\left[\left\|{\tilde{e}}_{t} + \frac{1}{M} \sum_{i=1}^M e_{t,i}\right\|_2^2\right] \leq  \frac{8(1-\delta)G^2}{\delta^2}\left[1 + \frac{16}{\delta^2}\right]$
for any $t \geq 0$.
\end{lemma}
This implies that $\nabla F({\tilde{x}}_t) \approx \nabla F(x_t)$ by Assumption~\ref{assumption:ef_smooth}.
Given the above results, we can prove convergence of the proposed method by
utilizing tools used on the full-precision distributed SGD.

\begin{theorem} \label{theorem:general_distefsgd_convergence}
Suppose that Assumptions~\ref{assumption:ef_smooth}-\ref{assumption:ef_bounded_grad} hold.
Assume that $0 < \eta_t < 3/(2L)$ for all $t$.
For the $\{x_t\}$
sequence generated from Algorithm~\ref{alg:dist-ef-sgd}, we have
\begin{eqnarray*}
{\mathbb E}\left[\left\|\nabla F(x_o)\right\|^2_2\right]
&\leq&  \frac{4}{\sum_{k=0}^{T-1}\eta_k\left(3 - 2L\eta_k \right)}[F(x_0) - F_*]  + \frac{2L\sigma^2}{M}\sum_{t=0}^{T-1}\frac{\eta_t^2}{\sum_{k=0}^{T-1}\eta_k\left(3 - 2L\eta_k \right)} \\
&& + \frac{32L^2(1 - \delta)G^2}{\delta^2}\left[1 + \frac{16}{\delta^2}\right]\sum_{t=0}^{T-1}\frac{\eta_t\eta_{t-1}^2}{\sum_{k=0}^{T-1}\eta_k\left(3 - 2L\eta_k \right)},
\end{eqnarray*}
where $o \in \{0, \dots, T-1\}$ is an index such that $P(o = k) = \frac{\eta_k\left(3 - 2L\eta_k \right)}{\sum_{t=0}^{T-1}\eta_t\left(3 - 2L\eta_t \right)}$, $\forall k = 0, \dots, T-1$.
\end{theorem}
The first term on the RHS
shows decay of the initial value.
The second term
is related to the variance,
and the proposed algorithm enjoys variance reduction with more workers.
The last term
is due to gradient compression.
A large
$\delta$ (less
compression)
makes this term smaller and thus faster convergence.
Similar to the results in \cite{karimireddy2019error}, our bound
also holds
for unbiased compressors (e.g., QSGD \cite{alistarh2017qsgd}) of the form $\mathcal{C}(\cdot) = cU(\cdot)$, where ${\mathbb E}[U(x)] = x$ and
${\mathbb E}[\|U(x)\|_2^2]\leq\frac{1}{c}\|x\|_2^2$ for some $0 < c < 1$. Then, $cU(\cdot)$ is a $c$-approximate compressor in expectation.

The following Corollary
shows that dist-EF-SGD has a convergence rate of
$\mathcal{O}(1/\sqrt{MT})$, leading to a $\mathcal{O}(1/(M\epsilon^4))$ iteration complexity for satisfying ${\mathbb E}[\|\nabla F(x_o)\|_2^2] \leq \epsilon^2$.

\begin{corollary} \label{corollary:convergence-dist-ef-sgd}
Let
stepsize
$\eta = \min(\frac{1}{2L}, \frac{\gamma}{\sqrt{T}/\sqrt{M} + (1 - \delta)^{1/3}\left(1/\delta^2 + 16/\delta^4\right)^{1/3}T^{1/3}})$
for some $\gamma > 0$.
Then,
\begin{eqnarray*}
{\mathbb E}[\|\nabla F(x_o)\|_2^2]
 & \leq &  \frac{4L}{T}[F(x_0) - F_*] + \left[\frac{2}{\gamma}[F(x_0) - F_*] +  L\gamma\sigma^2\right]\frac{1}{\sqrt{MT}} \\
&& + \frac{2(1 - \delta)^{1/3}\left[\frac{1}{\gamma}[F(x_0) - F_*] + 8L^2\gamma^2G^2\right]}{\delta^{2/3}T^{2/3}}\left[1 + \frac{16}{\delta^2}\right]^{1/3}.
\end{eqnarray*}
In comparison,
under the same assumptions,
distributed synchronous SGD
achieves
\begin{eqnarray*} \label{eq:dist_sgd_rate}
{\mathbb E}[\|\nabla F(x_o)\|^2_2]
 \leq  \frac{8L}{3T}[F(x_0) - F_*] + \left[\frac{2}{\gamma}[F(x_0) - F_*]  + L\gamma\sigma^2\right]\frac{2}{3\sqrt{MT}}.
\end{eqnarray*}
\end{corollary}
Thus, the convergence rate of
dist-EF-SGD
matches that of distributed synchronous SGD (with full-precision gradients)
after $T\geq O(1/\delta^2)$ iterations,
even though gradient compression is used.
Moreover, more
workers
(larger $M$) leads to faster convergence.
Note that the bound above does not reduce to that of EF-SGD when $M=1$, as we have two-way compression.
When $M=1$, our bound also differs from Remark~4 in \cite{karimireddy2019error} in that
our last term is
$O((1-\delta)^{1/3}/(\delta^{4/3}T^{2/3}))$, while theirs is
$O((1-\delta)/(\delta^2T))$ (which is for single machine with one-way compression). Ours is worse
by a factor of
$O(T^{1/3}\delta^{2/3}/(1-\delta)^{2/3})$, which is the price to pay for two-way
compression and a linear
speedup of using $M$ workers.
Moreover, unlike signSGD with majority vote \cite{pmlr-v80-bernstein18a}, we
achieve a convergence rate of $\mathcal{O}(1/\sqrt{MT})$ without assuming
a large mini-batch size ($=T$)
and
unimodal symmetric gradient noise.

Theorem~\ref{theorem:general_distefsgd_convergence}
only requires $0 < \eta_t < 3/(2L)$ for all $t$. This
thus allows the use
of any decreasing, increasing, or hybrid stepsize schedule.
In particular, we have the following Corollary.

\begin{corollary} \label{corollary:convergence-decrease-dist-ef-sgd}
Let $\eta_t = \frac{\gamma}{((t+1)T)^{1/4}/(\sqrt{M}) + (1 -
\delta)^{1/3}\left(1/\delta^2 + 16/\delta^4\right)^{1/3}T^{1/3}}$ (decreasing stepsize) with $T \geq 16L^4\gamma^4M^2$ or $\eta_t = \frac{\gamma\sqrt{t+1}}{T/\sqrt{M} + (1 - \delta)^{1/3}\left(1/\delta^2 + 16/\delta^4\right)^{1/3}T^{5/6}}$ (increasing stepsize) with $T \geq 4L^2\gamma^2M$. Then, dist-EF-SGD converges to a stationary point at a rate of $\mathcal{O}(1/\sqrt{MT})$.
\end{corollary}
To the best of our knowledge, this is the first such result for distributed
compressed SGD
with decreasing/increasing stepsize on nonconvex problems. These two stepsize
schedules can also be used together. For example, one can
use an increasing stepsize at the beginning of training as warm-up, and then a
decreasing stepsize afterwards.

\subsection{Blockwise Compressor}
\label{sec:blc_cp}

\begin{algorithm}[t]
\caption{Distributed Blockwise SGD with Error-Feedback (dist-EF-blockSGD)}
\label{alg:dist-bef-sgd}
\begin{algorithmic}[1]
    \STATE \textbf{Input:} stepsize sequence $\{\eta_t\}$ with $\eta_{-1} = 0$; number of workers $M$; block partition $\{\mathcal{G}_1, \dots, \mathcal{G}_B\}$.
    \STATE \textbf{Initialize:} $x_0 \in \mathbb{R}^d$; $e_{0,i} = 0 \in \mathbb{R}^d$ on each worker $i$; $\tilde{e}_{0} = 0 \in \mathbb{R}^d$ on server
    \FOR{$t = 0, \ldots, T-1$}
        \STATE \textbf{on each worker $i$}

         \STATE \qquad $p_{t,i} = g_{t,i} +  \frac{\eta_{t-1}}{\eta_t}e_{t,i}$ \algorithmiccomment{stochastic gradient $g_{t,i} = \nabla f(x_t, \xi_{t, i})$}
        \STATE \qquad \textbf{push} $\Delta_{t,i} = \left[\frac{\|p_{t, i, \mathcal{G}_1}\|_1}{d_1}\text{sign}(p_{t, i, \mathcal{G}_1}), \dots, \frac{\|p_{t, i, \mathcal{G}_B}\|_1}{d_B}\text{sign}(p_{t, i, \mathcal{G}_B})\right]$ \textbf{to server}

        \STATE \qquad $x_{t+1} = x_t - \eta_t\tilde{\Delta}_{t}$ \algorithmiccomment{$\tilde{\Delta}_{t}$ \textbf{is pulled from server}}
        \STATE \qquad $e_{t+1,i} = p_{t,i} - \Delta_{t,i}$

        \STATE \textbf{on server}
        \STATE \qquad \textbf{pull} $\Delta_{t,i}$ \textbf{from each worker} $i$ and ${\tilde{p}}_{t} = \frac{1}{M} \sum_{i=1}^M \Delta_{t,i} +  \frac{\eta_{t-1}}{\eta_t}{\tilde{e}}_{t}$
        \STATE \qquad \textbf{push} $\tilde{\Delta}_t = \left[\frac{\|{\tilde{p}}_{t,  \mathcal{G}_1}\|_1}{d_1}\text{sign}({\tilde{p}}_{t, \mathcal{G}_1}), \dots, \frac{\|{\tilde{p}}_{t,  \mathcal{G}_B}\|_1}{d_B}\text{sign}({\tilde{p}}_{t, \mathcal{G}_B})\right]$ \textbf{to each worker}
        \STATE \qquad $\tilde{e}_{t+1} = {\tilde{p}}_t - {\tilde{\Delta}}_t$
    \ENDFOR
\end{algorithmic}
\end{algorithm}

A commonly used compressor
is
\cite{karimireddy2019error}:
\begin{equation} \label{eq:orig}
\mathcal{C}(v) = \|v\|_1/d \cdot\text{sign}(v).
\end{equation}
Compared to using only the sign operator as in signSGD, the factor $\|v\|_1/d$ can preserve the gradient's magnitude.
However, as shown in
\cite{karimireddy2019error}, its
$\delta$ in Definition~\ref{definition:ef_compressor}
is  $\|v\|_1^2/(d\|v\|_2^2)$,
and can be
particularly small when $v$ is sparse.
When
$\delta$ is closer to $1$,
the bound
in Corollary~\ref{corollary:convergence-dist-ef-sgd}
becomes smaller  and thus
convergence is
faster.
In this section,  we achieve this
by
proposing a blockwise extension of
(\ref{eq:orig}).

Specifically, we
partition the
compressor input
$v$
into $B$ blocks, where
each block $b$ has $d_b$ elements indexed by
$\mathcal{G}_b$.
Block $b$
is then compressed with
scaling factor
$\|v_{\mathcal{G}_b}\|_1/d_b$ (where $v_{\mathcal{G}_b}$ is the subvector of $v$ with elements in
block $b$),
leading to:
$\mathcal{C}_B(v) = [\|v_{\mathcal{G}_1}\|_1/d_1 \cdot \text{sign}(v_{\mathcal{G}_1}), \dots,
\|v_{\mathcal{G}_B}\|_1/d_B\cdot\text{sign}(v_{\mathcal{G}_B})]$.
A similar
compression scheme,
with each layer being a block,
is considered in the experiments
of \cite{karimireddy2019error}. However, they provide no
theoretical justifications.
The following Proposition first shows that $\mathcal{C}_B(\cdot)$ is also an approximate compressor.

 \begin{proposition} \label{prop:block_density}
  Let $[B] = \{1, 2, \dots, B\}$.
 $\mathcal{C}_B$

 is a
 $\phi(v)$-approximate compressor, where

 $\phi(v)= \min_{b \in [B]} \frac{\|v_{\mathcal{G}_b}\|_1^2}{d_b\|v_{\mathcal{G}_b}\|_2^2} \geq \min_{b \in [B]} \frac{1}{d_b}$.

 \end{proposition}
The resultant algorithm
will be
called dist-EF-blockSGD
(Algorithm~\ref{alg:dist-bef-sgd}) in the sequel.
As can be seen, this is a special case of Algorithm~\ref{alg:dist-ef-sgd}.
By replacing $\delta$ with $\phi(v)$ in
  Proposition~\ref{prop:block_density},
the convergence results of dist-EF-SGD in Section~\ref{sec:dist-ef-sgd} can be directly applied.

There are many ways to partition the gradient into blocks. In practice, one can simply consider
each parameter
tensor/matrix/vector in the deep network as a block. The intuition is that (i)
gradients in the same parameter tensor/matrix/vector typically have similar
magnitudes,
and (ii) the corresponding scaling factors can thus be tighter than the scaling factor obtained on the whole parameter, leading to a larger $\delta$.
As an illustration of (i),
Figure~\ref{fig:cv_b16_cifar100} shows
the
coefficient of variation
(which
is defined as the ratio of the standard deviation to the mean)
of $\{|g_{t, i}|\}_{i \in \mathcal{G}_b}$
averaged over all blocks and iterations in an epoch,
obtained from ResNet-20 on the CIFAR-100 dataset (with a mini-batch size of 16 per worker).\footnote{The detailed experimental setup is in  Section~\ref{sec:distef-cifar100-expt}.}
A value smaller than $1$ indicates that the absolute gradient values in each block concentrate around the mean.
As for point (ii) above,
consider the case where all the blocks
are of the same size ($d_b = \tilde{d}, \forall b$), elements in the
same block have the same magnitude
($\forall i \in \mathcal{G}_b,
|v_i| = c_b$ for some $c_b$),
and
the magnitude is increasing across blocks ($c_b/c_{b+1}=\alpha$ for some $\alpha < 1$).
For the standard compressor in
(\ref{eq:orig}),
$\delta=\frac{\|v\|_1^2}{d\|v\|_2^2} = \frac{(1+\alpha)(1-\alpha^B)}{B(1-\alpha)(1+\alpha^B)} \approx
\frac{(1+\alpha)}{B(1-\alpha)}$ for a sufficiently large $B$;
whereas for the proposed blockwise compressor,
$\phi(v) = 1
\gg \frac{(1+\alpha)}{B(1-\alpha)}$.
 Figure~\ref{fig:delta_b16_cifar100}
shows the empirical estimates
of
$\|v\|_1^2/(d\|v\|_2^2)$ and $\phi(v)$
in the ResNet-20 experiment.
As can be seen, $\phi(v)
\gg \|v\|_1^2/(d\|v\|_2^2)$.

\begin{figure}[ht]
\vskip -.2in
\begin{center}
\subfigure[Coefficient of variation of $\{|g_{t, i}|\}_{i \in \mathcal{G}_b}$.]{\includegraphics[width=0.4\columnwidth]{cv} \label{fig:cv_b16_cifar100}}
\subfigure[$\delta$ for blockwise and non-block versions.]{\includegraphics[width=0.4\columnwidth]{delta_block} \label{fig:delta_b16_cifar100}}
\vskip -.1in
\caption{Illustrations using the ResNet-20 in
Section~\ref{sec:distef-cifar100-expt}.
Left: Averaged coefficient of variation of $\{|g_{t, i}|\}_{i \in \mathcal{G}_b}$.
Right: Empirical estimates of $\delta$ for the blockwise ($\phi(v)$ in Proposition~\ref{prop:block_density}) and non-block versions ($\|v\|_1^2/(d\|v\|_2^2)$).
Each point is the minimum among all iterations in an epoch.
The lower bounds,
$\min_{b \in [B]} 1/d_b$ and $1/d$,
are also shown.
Note that the ordinate is in log scale.
}
\end{center}
\end{figure}

The
per-iteration communication costs of the various distributed algorithms
are shown in
Table~\ref{tbl:bit_comp}.
Compared to signSGD with majority vote \cite{pmlr-v80-bernstein18a},
dist-EF-blockSGD requires an extra
$64MB$ bits
for transmitting the blockwise scaling factors
(each factor $\|v_{\mathcal{G}_b}\|_1/d_b$ is stored
in float32 format and transmitted twice in each iteration).
By treating each vector/matrix/tensor parameter as a block,
$B$ is typically in the order of hundreds. For most
problems of interest, $64MB/(2Md) < 10^{-3}$.
The reduction in communication cost
compared to full-precision distributed SGD
is thus nearly 32x.

 \begin{table}[h]
 \vskip -.1in
 \begin{center}
 \caption{Communication costs of the various distributed gradient compression
 algorithms and SGD.}
 \vskip .1in
 \begin{tabular}{c|c}
 \hline
 algorithm & \#bits per iteration \\
  \hline
 full-precision SGD & $64Md$ \\
  signSGD with majority vote & $2Md$ \\
  dist-EF-blockSGD & $2Md + 64MB$\\
  \hline
 \end{tabular}
  \label{tbl:bit_comp}
 \end{center}
 \vskip -.2in
 \end{table}

\subsection{Nesterov's Momentum}
\label{sec:dist-ef-sgdm}

Momentum has been widely used in deep networks
\cite{sutskever2013importance}.
Standard distributed SGD with
Nesterov's momentum
\cite{nesterov1983method}
and full-precision gradients
uses the update: $m_{t,i} = \mu m_{t-1,i} + g_{t, i}, \forall i \in [M]$ and $x_{t+1} =x_t -  \eta_t \frac{1}{M} \sum_{i=1}^M(\mu m_{t,i} + g_{t, i})$,
where
$m_{t, i}$
is a local momentum vector maintained by each worker $i$ at time $t$ (with
$m_{0,i} = 0$), and $\mu \in [0, 1)$ is the momentum parameter.
In this section, we extend
the proposed dist-EF-SGD
with momentum.
Instead of sending the compressed $g_{t, i} +
\frac{\eta_{t-1}}{\eta_t}e_{t, i}$
to the server,
the compressed
$\mu m_{t, i} + g_{t, i} + \frac{\eta_{t-1}}{\eta_t}e_{t, i}$ is sent.
The server merges all the workers's results
and sends it back to each worker.
The resultant procedure with blockwise compressor is called dist-EF-blockSGDM (Algorithm~\ref{alg:dist-bef-sgdm}),
and has the same communication cost as dist-EF-blockSGD.
The
corresponding non-block variant is analogous.

\vskip -.1in
\begin{algorithm}[ht]
\caption{Distributed Blockwise Momentum SGD with Error-Feedback (dist-EF-blockSGDM)}
\label{alg:dist-bef-sgdm}
\begin{algorithmic}[1]
    \STATE \textbf{Input:} stepsize sequence $\{\eta_t\}$ with $\eta_{-1} = 0$; momentum parameter $0 \leq \mu < 1$; number of workers $M$; block partition $\{\mathcal{G}_1, \dots, \mathcal{G}_B\}$.
    \STATE \textbf{Initialize:} $x_0 \in \mathbb{R}^d$; $m_{-1, i} = e_{0,i} = 0 \in \mathbb{R}^d$ on each worker $i$; $\tilde{e}_{0} = 0 \in \mathbb{R}^d$ on server
    \FOR{$t = 0, \ldots, T-1$}
        \STATE \textbf{on each worker $i$}

         \STATE \qquad $m_{t,i} = \mu m_{t - 1,i} +  g_{t,i}$ \algorithmiccomment{stochastic gradient $g_{t,i} = \nabla f(x_t, \xi_{t, i})$}
        \STATE \qquad $p_{t,i} = \mu m_{t,i} + g_{t,i} +  \frac{\eta_{t-1}}{\eta_t}e_{t,i}$
        \STATE \qquad \textbf{push} $\Delta_{t,i} = \left[\frac{\|p_{t, i, \mathcal{G}_1}\|_1}{d_1}\text{sign}(p_{t, i, \mathcal{G}_1}), \dots, \frac{\|p_{t, i, \mathcal{G}_B}\|_1}{d_B}\text{sign}(p_{t, i, \mathcal{G}_B})\right]$ \textbf{to server}

        \STATE \qquad $x_{t+1} = x_t - \eta_t\tilde{\Delta}_{t}$ \algorithmiccomment{$\tilde{\Delta}_{t}$ \textbf{is pulled from server}}
        \STATE \qquad $e_{t+1,i} = p_{t,i} - \Delta_{t,i}$

        \STATE \textbf{on server}
        \STATE \qquad \textbf{pull} $\Delta_{t,i}$ \textbf{from each worker} $i$ and ${\tilde{p}}_{t} = \frac{1}{M} \sum_{i=1}^M \Delta_{t,i} +  \frac{\eta_{t-1}}{\eta_t}{\tilde{e}}_{t}$
        \STATE \qquad \textbf{push} $\tilde{\Delta}_t = \left[\frac{\|{\tilde{p}}_{t,  \mathcal{G}_1}\|_1}{d_1}\text{sign}({\tilde{p}}_{t, \mathcal{G}_1}), \dots, \frac{\|{\tilde{p}}_{t,  \mathcal{G}_B}\|_1}{d_B}\text{sign}({\tilde{p}}_{t, \mathcal{G}_B})\right]$ \textbf{to each worker}
        \STATE \qquad $\tilde{e}_{t+1} = {\tilde{p}}_t - {\tilde{\Delta}}_t$
    \ENDFOR
\end{algorithmic}
\end{algorithm}
\vskip -.1in

Similar to Lemma~\ref{lemma:recurrence_ec_iterate}, the following Lemma shows that
the error-corrected iterate ${\tilde{x}}_t$ is very similar to Nesterov's accelerated
gradient iterate, except that the momentum is computed based on $\{x_t\}$.
\begin{lemma} \label{lemma:recurrence_mec_iterate}
The error-corrected iterate ${\tilde{x}}_t = x_t - \eta_{t-1}({\tilde{e}}_t + \frac{1}{M} \sum_{i=1}^M e_{t,i})$, where
$x_t$, ${\tilde{e}}_t$, and $e_{t,i}$'s are generated from
Algorithm~\ref{alg:dist-bef-sgdm},
satisfies the recurrence:
${\tilde{x}}_{t+1} = {\tilde{x}}_t -  \eta_t \frac{1}{M} \sum_{i=1}^M(\mu m_{t,i} + g_{t, i})$.
\end{lemma}
As in Section~\ref{sec:dist-ef-sgd},
it can be shown that
$\|{\tilde{e}}_t + \frac{1}{M} \sum_{i=1}^M e_{t,i}\|_2$ is bounded and
$\nabla F({\tilde{x}}_t) \approx \nabla F(x_t)$.
The following Theorem shows
the convergence rate of the proposed
dist-EF-blockSGDM.

 \begin{theorem} \label{theorem:general_distefsgdm_convergence}
Suppose that Assumptions~\ref{assumption:ef_smooth}-\ref{assumption:ef_bounded_grad} hold.
Let $\eta_t = \eta$ for some $\eta > 0$.
For any $\eta \leq \frac{(1 - \mu)^2}{2L}$, and the $\{x_t\}$ sequence generated from Algorithm~\ref{alg:dist-bef-sgdm}, we have
\begin{eqnarray}
{\mathbb E}\left[\left\|\nabla F(x_o)\right\|^2_2\right]
& \leq &  \frac{4(1 - \mu)}{\eta T}[F(x_0) - F_*]
+ \frac{2L\eta\sigma^2}{(1 - \mu)M}\left[1 + \frac{2L\eta\mu^4}{(1 - \mu)^3}\right] \label{bound:general_distefsgdm_convergence} \\
&& + \frac{32L^2\eta^2(1-\delta)G^2}{\delta^2(1-\mu)^2}\left[1 + \frac{16}{\delta^2}\right].  \nonumber
\end{eqnarray}
\end{theorem}
Compared to Theorem~\ref{theorem:general_distefsgd_convergence}, using a larger momentum parameter
$\mu$ makes the first term (which depends on the initial condition) smaller but
a worse variance
term (second term) and error term due to gradient compression
(last term).
Similar to Theorem~\ref{theorem:general_distefsgd_convergence}, a larger $\eta$ makes the third term larger.
The following Corollary shows that the proposed dist-EF-blockSGDM achieves a
convergence rate of $\mathcal{O}(((1 - \mu)[F(x_0) - F_*] + \sigma^2/(1-\mu))/\sqrt{MT})$.

\begin{corollary} \label{corollary:convergence-dist-ef-sgdm}
Let $\eta = \frac{\gamma}{\sqrt{T}/\sqrt{M} + (1 - \delta)^{1/3}\left(1/\delta^2 + 16/\delta^4\right)^{1/3}T^{1/3}}$ for some $\gamma > 0$.
For any $T \geq \frac{4\gamma^2L^2M}{(1 - \mu)^4}$,
${\mathbb E}\left[\left\|\nabla F(x_o)\right\|^2_2\right]
\leq   \left[\frac{2(1 - \mu)}{\gamma}[F(x_0) - F_*]  + \frac{L\gamma\sigma^2}{1-\mu}\right]\frac{2}{\sqrt{MT}}
+ \frac{4L^2\gamma^2\mu^4\sigma^2}{(1 - \mu)^4T}
+ \frac{4(1 - \delta)^{1/3}\left[\frac{(1 - \mu)}{\gamma}[F(x_0) - F_*] + 8L^2\gamma^2G^2/(1-\mu)^2\right]}{\delta^{2/3}T^{2/3}}\left[1 + \frac{16}{\delta^2}\right]^{1/3}$.
\end{corollary}

 \section{Experiments}

\subsection{Multi-GPU Experiment on CIFAR-100}
\label{sec:distef-cifar100-expt}

In this experiment, we demonstrate that the proposed
 dist-EF-blockSGDM

 and dist-EF-blockSGD ($\mu = 0$ in Algorithm~\ref{alg:dist-bef-sgdm}),
though using fewer bits for gradient transmission,
still has good convergence.
For faster experimentation,
we use a
single node
with multiple GPUs (an AWS P3.16 instance with 8 Nvidia V100 GPUs,
each GPU being a worker)
instead of a distributed setting.

Experiment is performed on
the CIFAR-100 dataset, with
50K training images and 10K test images.
We use a 20-layer ResNet \cite{he2016identity}.
Each parameter tensor/matrix/vector is treated as a block
in dist-EF-blockSGD(M).
They are compared with
 (i) distributed synchronous SGD (with full-precision gradient);
 (ii) distributed synchronous SGD
 (full-precision gradient)
 with momentum (SGDM);
 (iii) signSGD with majority vote \cite{pmlr-v80-bernstein18a}; and
 (iv) signum with majority vote \cite{bernstein2019signsgd}.
All the algorithms are implemented in MXNet.
We vary the mini-batch size per worker in $\{8, 16, 32\}$.
Results are averaged over 5 repetitions.
More details of the experiments are shown
in
Appendix~\ref{app:cifar}.

Figure~\ref{fig:ef_cifar100} shows
convergence of the testing accuracy w.r.t.
the number of epochs.
As can be seen,
dist-EF-blockSGD converges as fast as SGD and has slightly better accuracy, while
signSGD performs poorly. In particular, dist-EF-blockSGD is robust to the mini-batch size, while the performance of signSGD
degrades with smaller mini-batch size (which agrees with the results in
\cite{pmlr-v80-bernstein18a}).
Momentum makes SGD and dist-EF-blockSGD faster with mini-batch size of $16$ or $32$ per worker,
particularly before epoch $100$.
At epoch 100, the learning rate is reduced, and the difference is less obvious.
This is because a larger
mini-batch
means smaller
variance
$\sigma^2$,
so the initial optimality gap
$F(x_0) - F_*$ in (\ref{bound:general_distefsgdm_convergence})
is more dominant.
Use of momentum
($\mu > 0$) is then beneficial.
On the other hand, momentum significantly
improves signSGD.
However, signum is still much worse than dist-EF-blockSGDM.

 \begin{figure}[t]
\begin{center}
\subfigure{\includegraphics[width=0.32\columnwidth]{dist-ef-sgd-cifar100-b8.pdf}}
\subfigure{\includegraphics[width=0.32\columnwidth]{dist-ef-sgd-cifar100-b16.pdf}}
\subfigure{\includegraphics[width=0.32\columnwidth]{dist-ef-sgd-cifar100-b32.pdf}} \\
\setcounter{subfigure}{0}
\vskip -.18in
\subfigure[Mini-batch size: $8$ per worker.]{\includegraphics[width=0.32\columnwidth]{dist-ef-sgdm-cifar100-b8.pdf}}
\subfigure[Mini-batch size: $16$ per worker.]{\includegraphics[width=0.32\columnwidth]{dist-ef-sgdm-cifar100-b16.pdf}}
\subfigure[Mini-batch size: $32$ per worker.]{\includegraphics[width=0.32\columnwidth]{dist-ef-sgdm-cifar100-b32.pdf}}
\vskip -.1in
\caption{Testing accuracy on CIFAR-100. Top: No momentum; Bottom: With
momentum. The solid curve is the mean accuracy
over five repetitions. The shaded region
spans one standard deviation.}
\label{fig:ef_cifar100}
\end{center}
\vskip -.3in
\end{figure}

\subsection{Distributed Training on ImageNet}
\label{sec:distef-imagenet-expt}

In this section, we perform distributed optimization on ImageNet \cite{russakovsky2015imagenet}
using
a 50-layer ResNet.
Each worker is an AWS P3.2 instance with 1 GPU,
and the parameter server is housed in one node.
We use the publicly available
code\footnote{\url{https://github.com/PermiJW/signSGD-with-Majority-Vote}} in
\cite{bernstein2019signsgd}, and the default communication library
Gloo communication
library in PyTorch.
As in
\cite{bernstein2019signsgd},
we use its allreduce implementation
for SGDM, which is faster.

As momentum accelerates the training for large mini-batch size in Section~\ref{sec:distef-cifar100-expt},
we only compare the momentum variants here.
 The proposed dist-EF-blockSGDM is
 compared with
(i) distributed synchronous SGD with momentum (SGDM); and
 (ii) signum with majority vote \cite{bernstein2019signsgd}.
The number of workers $M$ is varied in $\{7, 15\}$.
With an odd number of workers, a majority vote will not
produce zero, and so signum does not lose accuracy by using 1-bit compression.
More details of the setup are in
Appendix~\ref{app:imagenet}.

\begin{figure}[t]
\begin{center}
\subfigure[Test accuracy w.r.t. epoch.]{\includegraphics[width=0.32\columnwidth]{imagenet-epoch-7x.pdf}}
\subfigure[Test accuracy w.r.t. time. \label{fig:t1}]{\includegraphics[width=0.32\columnwidth]{imagenet-time-7x.pdf}}
\subfigure[Workload breakdown. \label{fig:load1}]{\includegraphics[width=0.32\columnwidth]{imagenet-bar-7x.pdf}}
\vskip -.15in
\subfigure[Test accuracy w.r.t. epoch.]{\includegraphics[width=0.32\columnwidth]{imagenet-epoch-15x.pdf}}
\subfigure[Test accuracy w.r.t. time. \label{fig:t2}]{\includegraphics[width=0.32\columnwidth]{imagenet-time-15x.pdf}}
\subfigure[Workload breakdown. \label{fig:load2}]{\includegraphics[width=0.32\columnwidth]{imagenet-bar-15x.pdf}}
\vskip -.1in
\caption{Distributed training results on the ImageNet dataset.
Top:
7 workers;
Bottom:
15 workers.}
\label{fig:dist_imagenet}
\end{center}
\vskip -.3in
\end{figure}

Figure~\ref{fig:dist_imagenet} shows
the testing accuracy w.r.t. the number of epochs and wall clock time.
As in Section~\ref{sec:distef-cifar100-expt},
SGDM and dist-EF-blockSGDM
have comparable accuracies,
while signum is inferior.
When 7 workers are used, dist-EF-blockSGDM
has higher accuracy than SGDM
(76.77\% vs
76.27\%).
dist-EF-blockSGDM
reaches SGDM's highest accuracy
in around 13 hours,
while SGDM takes 24 hours (Figure~\ref{fig:t1}),
leading to a $46\%$ speedup.
With 15 machines, the improvement is smaller (Figure~\ref{fig:t2}). This is because the burden on the parameter server is heavier. We expect
comparable speedup with the 7-worker setting can be obtained by using more
parameter servers. In both cases, signum converges fast but the test
accuracies are about $4\%$ worse.

Figures~\ref{fig:load1}
and \ref{fig:load2}
show a breakdown
of wall clock time
into computation and communication time.\footnote{Following
\cite{bernstein2019signsgd}, communication time includes the extra computation time for
error feedback and compression.}
All methods have comparable computation costs,
but
signum and dist-EF-blockSGDM have lower communication costs than SGDM.
The communication costs
for signum and dist-EF-blockSGDM are comparable
for 7 workers,
but for 15 workers
signum is lower.
We speculate that
it is because
the sign vectors and scaling factors are sent separately to the server in our
implementation, which causes more latency on the server with more workers.
This may be alleviated if
the two operations
are fused.

 \section{Conclusion}

In this paper, we
proposed  a
distributed blockwise SGD algorithm with error feedback
and momentum.
 By partitioning the gradients into blocks, we can transmit each block of gradient using 1-bit quantization with its average $\ell_1$-norm.
 The proposed methods are communication-efficient and have the same convergence rates as
 full-precision distributed SGD/SGDM for nonconvex objectives.
Experimental results show that the proposed methods have
 fast convergence and achieve the same test accuracy as SGD/SGDM, while
 signSGD and signum only achieve much worse accuracies.