Robustness of accelerated first-order algorithms for strongly convex optimization problems

TL;DR

This paper investigates how accelerated first-order optimization algorithms perform under stochastic gradient noise in strongly convex problems, providing bounds on their robustness and analyzing the influence of problem and network structure.

Contribution

It offers tight bounds on mean-squared error under noise, analyzes the spectral influence on robustness, and applies findings to distributed averaging networks.

Findings

Upper bounds on mean-squared deviation are tight up to constants.

Spectral properties of the Hessian influence robustness.

Distributed averaging robustness depends on network size and topology.

Abstract

We study the robustness of accelerated first-order algorithms to stochastic uncertainties in gradient evaluation. Specifically, for unconstrained, smooth, strongly convex optimization problems, we examine the mean-squared error in the optimization variable when the iterates are perturbed by additive white noise. This type of uncertainty may arise in situations where an approximation of the gradient is sought through measurements of a real system or in a distributed computation over a network. Even though the underlying dynamics of first-order algorithms for this class of problems are nonlinear, we establish upper bounds on the mean-squared deviation from the optimal solution that are tight up to constant factors. Our analysis quantifies fundamental trade-offs between noise amplification and convergence rates obtained via any acceleration scheme similar to Nesterov's or heavy-ball methods. To gain additional analytical insight, for strongly convex quadratic problems, we explicitly evaluate the steady-state variance of the optimization variable in terms of the eigenvalues of the Hessian of the objective function. We demonstrate that the entire spectrum of the Hessian, rather than just the extreme eigenvalues, influence robustness of noisy algorithms. We specialize this result to the problem of distributed averaging over undirected networks and examine the role of network size and topology on the robustness of noisy accelerated algorithms.

Tables

Table 1. TABLE I: Conventional values of parameters and the corresponding rates for f ∈ ℱ m L 𝑓 superscript subscript ℱ 𝑚 𝐿 f\in\mathcal{F}_{m}^{L} , ‖ x t − x ⋆ ‖ ≤ c ​ ρ t ​ ‖ x 0 − x ⋆ ‖ , norm superscript 𝑥 𝑡 superscript 𝑥 ⋆ 𝑐 superscript 𝜌 𝑡 norm superscript 𝑥 0 superscript 𝑥 ⋆ \|x^{t}-x^{\star}\|\leq c\,\rho^{t}\,\|x^{0}-x^{\star}\|, where κ : = L / m \kappa\mathrel{\mathop{:}}=L/m and c > 0 𝑐 0 c>0 is a constant   [ 9 , Theorems 2.1.15, 2.2.1] . The heavy-ball method does not offer acceleration guarantees for all f ∈ ℱ m L 𝑓 superscript subscript ℱ 𝑚 𝐿 f\in\mathcal{F}_{m}^{L} .

Method	Parameters	Linear rate
Gradient	$\alpha =\frac{1}{L}$	$\rho =\sqrt{\mathrm{\hspace{0.17em}1}-\frac{2}{\kappa +1}}$
Nesterov	$\alpha =\frac{1}{L}$, $\beta =\frac{\sqrt{\kappa }-\mathrm{\hspace{0.17em}1}}{\sqrt{\kappa }+\mathrm{\hspace{0.17em}1}}$	$\rho =\sqrt{1-\frac{1}{\sqrt{\kappa }}}$

Table 2. TABLE II: Optimal parameters and the corresponding convergence rates for a strongly convex quadratic objective function f ∈ ℱ m L 𝑓 superscript subscript ℱ 𝑚 𝐿 f\in\mathcal{F}_{m}^{L} with λ max ​ ( ∇ 2 f ) = L subscript 𝜆 superscript ∇ 2 𝑓 𝐿 \lambda_{\max}(\nabla^{2}f)=L and λ min ​ ( ∇ 2 f ) = m subscript 𝜆 superscript ∇ 2 𝑓 𝑚 \lambda_{\min}(\nabla^{2}f)=m , and κ : = L / m \kappa\mathrel{\mathop{:}}=L/m   [ 49 , Proposition 1] .

Method	Optimal parameters	Rate of linear convergence
Gradient	$\alpha =\frac{2}{L+m}$	$\rho =\frac{\kappa -1}{\kappa +1}$
Nesterov	$\alpha =\frac{4}{3L+m}$, $\beta =\frac{\sqrt{3\kappa +1}-2}{\sqrt{3\kappa +1}+2}$	$\rho =\frac{\sqrt{3\kappa +1}-2}{\sqrt{3\kappa +1}}$
Heavy-ball	$\alpha =\frac{4}{{\left(\sqrt{L}+\sqrt{m}\right)}^2},\beta =\frac{{\left(\sqrt{\kappa }-1\right)}^2}{{\left(\sqrt{\kappa }+1\right)}^2}$	$\rho =\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}$

Equations

$$\underset{x}{\mathrm{minimize}}f(x)$$

$$x^{t+1}\;=\;x^{t}\;-\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\nabla f(x^{t})\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma w^{t}}$$

$$x^{t+2}\;=\;x^{t+1}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}(x^{t+1}\,-\,x^{t})\;-\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\nabla f(x^{t+1})\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma w^{t}}$$

$$\begin{array}[]{rcl}x^{t+2}&\!\!\!=\!\!\!&x^{t+1}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}(x^{t+1}\,-\,x^{t})\,-\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\nabla f\!\left(x^{t+1}\,+\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}(x^{t+1}\,-\,x^{t})\right)\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma w^{t}}\end{array}$$

$$f\in {\mathcal{F}}_m^L\iff mI⪯{\nabla }^{2}f(x)⪯LI,\forall x\in {\mathbb{R}}^n.$$

$$\Vert x^t-x^{\star }\Vert \le c{\rho }^t\Vert x^0-x^{\star }\Vert$$

$$J:=\underset{t\to \infty }{\mathrm{limsup}}\frac{1}{t}\sum _{k=0}^{t}E\left(\Vert x^k-x^{\star }{\Vert }^2\right).$$

$$f(x)=\frac{1}{2}{x}^{T}Qx-q^{T}x$$

$$L={\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_n=m>0.$$

$$\begin{array}[]{rcl}\psi^{t+1}&\!\!\!=\!\!\!&A\,\psi^{t}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma}B\,w^{t}\\[0.0pt] z^{t}&\!\!\!=\!\!\!&C\,\psi^{t}\end{array}$$

$$A=I-\alpha Q,B=C=I$$

$$A\;=\;\left[\begin{array}[]{cc}0&I\\ -\beta I&(1+\beta)I-\alpha Q\end{array}\right],~{}~{}A\;=\;\left[\begin{array}[]{cc}0&I\\ -\beta(I-\alpha Q)&(1+\beta)(I-\alpha Q)\end{array}\right]$$

$$\displaystyle B^{T}\;=\;\left[\begin{array}[]{cc}0&I\end{array}\right],~{}~{}C\;=\;\left[\begin{array}[]{cc}I&0\end{array}\right].$$

$$P^{t+1}\;=\;A\,P^{t}A^{T}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}BB^{T}$$

$$P:=\underset{t\to \infty }{\lim }E\left({\psi }^t({\psi }^t{)}^T\right)$$

$$P\;=\;A\,PA^{T}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}BB^{T}.$$

$$J=\underset{t\to \infty }{\lim }\frac{1}{t}\sum _{k=0}^{t}E\left(\Vert z^k{\Vert }^2\right)=\underset{t\to \infty }{\lim }E\left(\Vert z^t{\Vert }^2\right)$$

$$J=\underset{t\to \infty }{\lim }\frac{1}{t}\sum _{k=0}^t\mathrm{trace}\left(Z^k\right)=\mathrm{trace}\left(Z\right)$$

$${\hat{x}}^t:=V^T{x}^t,{\hat{x}}^{\star }:=V^T{x}^{\star },{\hat{w}}^t:=V^T{w}^t$$

$${\hat{\psi }}_i^t={\hat{x}}_i^t-{\hat{x}}_i^{\star },{\hat{A}}_i=1-\alpha {\lambda }_i,{\hat{B}}_i={\hat{C}}_i=1.$$

$$\displaystyle\hat{A}_{i}\;=\;\left[\begin{array}[]{cc}0&1\\ -\beta&1+\beta-\alpha\lambda_{i}\end{array}\right],~{}~{}$$

$$\displaystyle\hat{A}_{i}\;=\;\left[\begin{array}[]{cc}0&1\\ -\beta(1-\alpha\lambda_{i})&(1+\beta)(1-\alpha\lambda_{i})\end{array}\right],~{}~{}$$

$$\hat{P}_{i}\;=\;\hat{A}_{i}\,\hat{P}_{i}\,\hat{A}_{i}^{T}\;+\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}\hat{B}_{i}\hat{B}_{i}^{T},\quad i\,=\,1,\ldots,n$$

$$J=\sum _{i=1}^n\hat{J}({\lambda }_i):=\sum _{i=1}^n\mathrm{trace}({\hat{C}}_i{\hat{P}}_i{\hat{C}}_i^T)$$

$$\begin{array}[]{rl}\mbox{Gradient:}&\hat{J}_{\mathrm{gd}}(\lambda)\;=\;\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}}{\alpha\lambda\left(2\,-\,\alpha\lambda\right)}\\[5.69046pt] \mbox{Polyak:}&\hat{J}_{\mathrm{hb}}(\lambda)\;=\;\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}(1\,+\,\beta)}{\alpha\lambda\left(1\,-\,\beta\right)\left(2(1\,+\,\beta)\,-\,\alpha\lambda\right)}\\[7.11317pt] \mbox{Nesterov:}&\hat{J}_{\mathrm{na}}(\lambda)\;=\;\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}(1\,+\,\beta(1\,-\,\alpha\lambda))}{\alpha\lambda\left(1\,-\,\beta(1\,-\,\alpha\lambda)\right)\left(2(1\,+\,\beta)\,-\,(2\beta\,+\,1)\alpha\lambda\right)}.\end{array}$$

$$J^{\prime }=\underset{t\to \infty }{\lim }E\left((x^t-x^{\star }{)}^{T}Q(x^t-x^{\star })\right).$$

$$\frac{J_{\mathrm{hb}}}{J_{\mathrm{gd}}}=\frac{(\sqrt{\kappa }+1{)}^4}{8\sqrt{\kappa }(\kappa +1)}.$$

$$\frac{{\hat{J}}_{\mathrm{hb}}(\lambda )}{{\hat{J}}_{\mathrm{gd}}(\lambda )}=\frac{1}{1-{\beta }^2}=\frac{(\sqrt{\kappa }+1{)}^4}{8\sqrt{\kappa }(\kappa +1)},\forall \lambda \in [m,L].$$

$$\begin{array}[]{rclcl}\max\limits_{\lambda\,\in\,[m,L]}\;\hat{J}_{\mathrm{gd}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{gd}}(m)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{gd}}(L)\;=\;\dfrac{(\kappa\,+\,1)^{2}}{4\kappa}\\[7.11317pt] \min\limits_{\lambda\,\in\,[m,L]}\;\hat{J}_{\mathrm{gd}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{gd}}({1}/{\alpha})&\!\!\!=\!\!\!&1\end{array}$$

$$\begin{array}[]{rclcl}&&\hat{J}_{\mathrm{na}}(L)&\!\!\!=\!\!\!&\dfrac{9\,\bar{\kappa}^{2}\!\left(\bar{\kappa}\,+\,2\sqrt{\bar{\kappa}}\,-2\,\right)}{32\left(\bar{\kappa}\,-\,1\right)\!\left(\bar{\kappa}\,-\,\sqrt{\bar{\kappa}}\,+\,1\right)\!\left(2\sqrt{\bar{\kappa}}\,-\,1\right)}\\[9.95863pt] \max\limits_{\lambda\,\in\,[m,L]}\;\hat{J}_{\mathrm{na}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{na}}(m)&\!\!\!=\!\!\!&\dfrac{\bar{\kappa}^{2}\!\left(\bar{\kappa}\,-\,2\sqrt{\bar{\kappa}}\,+\,2\right)}{32\!\left(\sqrt{\bar{\kappa}}\,-\,1\right)^{3}}\\[4.26773pt] \min\limits_{\lambda\,\in\,[m,L]}\;\hat{J}_{\mathrm{na}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{na}}({1}/{\alpha})&\!\!\!=\!\!\!&1\end{array}$$

$$\begin{array}[]{rcl}\max\limits_{\lambda\,\in\,[m,L]}\hat{J}_{\mathrm{gd}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{gd}}(m)\;=\;\hat{J}_{\mathrm{gd}}(L)\;=\;\Theta(\kappa)\\[4.26773pt] \max\limits_{\lambda\,\in\,[m,L]}\hat{J}_{\mathrm{hb}}(\lambda)&\!\!\!=\!\!\!&\hat{J}_{\mathrm{hb}}(m)\;=\;\hat{J}_{\mathrm{hb}}(L)\;=\;\Theta(\kappa\sqrt{\kappa}).\end{array}$$

$$\underset{\lambda \in [m,L]}{\max }{\hat{J}}_{\mathrm{na}}(\lambda )={\hat{J}}_{\mathrm{na}}(m)=\Theta (\kappa \sqrt{\kappa }),{\hat{J}}_{\mathrm{na}}(L)=\Theta (\sqrt{\kappa }).$$

$$J=\sum _{i=1}^n\frac{{\sigma }^2(1+{\hat{\lambda }}_i{\hat{\lambda }}_i^{\prime })}{(1-{\hat{\lambda }}_i{\hat{\lambda }}_i^{\prime })(1-{\hat{\lambda }}_i)(1-{\hat{\lambda }}_i^{\prime })(1+{\hat{\lambda }}_i)(1+{\hat{\lambda }}_i^{\prime })}$$

$$\frac{{\hat{J}}_{\mathrm{na}}(m)+(n-1){\hat{J}}_{\mathrm{na}}(L)}{{\hat{J}}_{\mathrm{gd}}(m)+(n-1){\hat{J}}_{\mathrm{gd}}(L)}\le \frac{J_{\mathrm{na}}}{J_{\mathrm{gd}}}\le \frac{{\hat{J}}_{\mathrm{na}}(L)+(n-1){\hat{J}}_{\mathrm{na}}(m)}{{\hat{J}}_{\mathrm{gd}}(L)+(n-1){\hat{J}}_{\mathrm{gd}}(m)}.$$

$$\frac{(3\kappa +1{)}^{\frac{3}{2}}}{32}\le {\hat{J}}_{\mathrm{na}}(m)\le \frac{(3\kappa +1{)}^{\frac{3}{2}}}{8},\frac{9\sqrt{3\kappa +1}}{64}\le {\hat{J}}_{\mathrm{na}}(L)\le \frac{9\sqrt{3\kappa +1}}{8}.$$

Full text

I Introduction

First-order algorithms are well-suited for solving a broad range of optimization problems that arise in statistics, signal and image processing, control, and machine learning [1, 2, 3, 4, 5]. Among these algorithms, accelerated methods enjoy the optimal rate of convergence and they are popular because of their low per-iteration complexity. There is a large body of literature dedicated to the convergence analysis of these methods under different stepsize selection rules [6, 7, 8, 5, 2, 9]. In many applications, however, the exact value of the gradient is not fully available, e.g., when the objective function is obtained via costly simulations (e.g., tuning of hyper-parameters in supervised/unsupervised learning [10, 11, 12] and model-free optimal control [13, 14, 15]), when evaluation of the objective function relies on noisy measurements (e.g., real-time and embedded applications), or when the noise is introduced via communication between different agents (e.g., distributed computation over networks). Another related application arises in the context of (batch) stochastic gradient, where at each iteration the gradient of the objective function is computed from a small batch of data points. Such a batch gradient is known to be a noisy unbiased estimator for the gradient of the training loss. Moreover, additive noise may be introduced deliberately in the context of non-convex optimization to help the iterates escape saddle points and improve generalization [16, 17].

In all above situations, first-order algorithms only have access to noisy estimates of the gradient. This observation has motivated the robustness analysis of first-order algorithms under different types of noisy/inexact gradient oracles [18, 19, 20, 21, 22, 23]. For example, in a deterministic noise scenario, an upper bound on the error in iterates for accelerated proximal gradient methods was established in [24]. This study showed that both proximal gradient and its accelerated variant can maintain their convergence rates provided that the noise is bounded and that it vanishes fast enough. Moreover, it has been shown that in the presence of random noise, with the proper diminishing stepsize, acceleration can be achieved for general convex problems. However, in this case optimal rates are sub-linear [25].

In the context of stochastic approximation, while early results suggest to use a stepsize that is inversely proportional to the iteration number [19], a more robust behavior can be obtained by combining larger stepsizes with averaging [20, 26, 27, 28]. Utility of these averaging schemes and their modifications for solving quadratic optimization and manifold problems has been examined thoroughly in recent years [29, 30, 31]. Moreover, several studies have suggested that accelerated first-order algorithms are more susceptible to errors in the gradient compared to their non-accelerated counterparts [32, 33, 24, 21, 22, 34].

One of the basic sources of error that arises in computing the gradient can be modeled by additive white stochastic noise. This source of error is typical for problems in which the gradient is being sought through measurements of a real system [35] and it has a rich history in analysis of stochastic dynamical systems and control theory [36]. Moreover, in many applications including distributed computing over networks [37, 38], coordination in vehicular formations [39], and control of power systems [40], additive white noise is a convenient abstraction for the robustness analysis of distributed control strategies [38] and of first-order optimization algorithms [41, 42]. Motivated by this observation, in this paper we consider the scenario in which a white stochastic noise with zero mean and identity covariance is added to the iterates of standard first-order algorithms: gradient descent, Polyak’s heavy-ball method, and Nesterov’s accelerated algorithm. By confining our attention to smooth strongly convex problems, we provide a tight quantitative characterization for the mean-squared error of the optimization variable. Since this quantity provides a measure of how noise gets amplified by the dynamics resulting from optimization algorithms, we also refer to it as noise (or variance) amplification. We demonstrate that our quantitative characterization allows us to identify fundamental trade-offs between the noise amplification and the rate of convergence obtained via acceleration.

This work builds on our recent conference papers [43, 44]. In a concurrent work [45], a similar approach was taken to analyze the robustness of gradient descent and Nesterov’s accelerated method. Therein, it was shown that for a given convergence rate, one can select the algorithmic parameters such that the steady-state mean-squared error in the objective value of a Nesterov-like method becomes smaller than that of gradient descent. This is not surprising because gradient descent can be viewed as a special case of Nesterov’s method with a zero momentum parameter. Using this argument, similar assertions have been made about the variance amplification of the iterates. This observation has been used to design an optimal multi-stage algorithm that does not require any information about the variance of the noise [46]. On the contrary, we demonstrate that there are fundamental differences between these two robustness measures, i.e., objective values and iterates, as the former does not capture the negative impact of acceleration in the presence of noise.

Focusing on the error in the iterates, we show that any choice of parameters for Nesterov’s or heavy-ball methods that yields an accelerated convergence rate increases variance amplification relative to gradient descent. More precisely, for the problem with the condition number $\kappa$, an algorithm with accelerated convergence rate of at least $1-c/\sqrt{\kappa}$, where $c$ is a positive constant, increases the variance amplification in the iterates by a factor of $\sqrt{\kappa}$. The robustness problem was also studied in [47] where the authors show a similar behavior of Nesterov’s method and gradient descent in an asymptotic regime in which the stepsize goes to zero. In contrast, we focus on the non-asymptotic stepsize regime and establish fundamental differences between gradient descent and its accelerated variants in terms of noise amplification.

More recently, the problem of finding upper bounds on the variance amplification was cast as a semidefinite program [48]. This formulation provided numerical results that are consistent with our theoretical upper bounds in terms of the condition number. In [48], structured objective functions (e.g., diagonal Hessians) that arise in distributed optimization were also studied and the problem of designing robust algorithms were formulated as a bilinear matrix inequality (which, in general, is not convex).

Contributions

The effect of imperfections on the performance and robustness of first-order algorithms has been studied in [30, 22] but the influence of acceleration on stochastic gradient perturbations has not been precisely characterized. We employ control-theoretic tools suitable for analyzing stochastic dynamical systems to quantify such influence and identify fundamental trade-offs between acceleration and noise amplification. The main contributions of this paper are:

We start our analysis by examining strongly convex quadratic optimization problems for which we can explicitly characterize variance amplification of first-order algorithms and obtain analytical insight. In contrast to convergence rates, which solely depend on the extreme eigenvalues of the Hessian matrix, we demonstrate that the variance amplification is influenced by the entire spectrum. 2. 2.

We establish the relation between the noise amplification of accelerated algorithms and gradient descent for parameters that provide the optimal convergence rate for strongly convex quadratic problems. We also explain how the distribution of the eigenvalues of the Hessian influences these relations and provide examples to show that acceleration can significantly increase the noise amplification. 3. 3.

We address the problem of tuning the algorithmic parameters and demonstrate the existence of a fundamental trade-off between the rate of convergence and noise amplification: for problems with condition number $\kappa$ and bounded dimension $n$, we show that any choice of parameters in accelerated methods that yields the linear convergence rate of at least $1-{c}/\sqrt{\kappa}$, where $c$ is a positive constant, increases noise amplification in the iterates relative to gradient descent by a factor of at least $\sqrt{\kappa}$. 4. 4.

We extend our analysis from quadratic objective functions to general strongly convex problems. We borrow an approach based on linear matrix inequalities from control theory to establish upper bounds on the noise amplification of both gradient descent and Nesterov’s accelerated algorithm. Furthermore, for any given condition number, we demonstrate that these bounds are tight up to constant factors. 5. 5.

We apply our results to distributed averaging over large-scale undirected networks. We examine the role of network size and topology on noise amplification and further illustrate the subtle influence of the entire spectrum of the Hessian matrix on the robustness of noisy optimization algorithms. In particular, we identify a class of large-scale problems for which accelerated Nesterov’s method achieves the same order-wise noise amplification (in terms of condition number) as gradient descent.

Paper structure

The rest of our presentation is organized as follows. In Section II, we formulate the problem and provide background material. In Section III, we explicitly evaluate the variance amplification (in terms of the algorithmic parameters and problem data) for strongly convex quadratic problems, derive lower and upper bounds, and provide a comparison between the accelerated methods and gradient descent. In Section IV, we extend our analysis to general strongly convex problems. In Section V, we establish fundamental trade-offs between the rate of convergence and noise amplification. In Section VI, we apply our results to the problem of distributed averaging over noisy undirected networks. We highlight the subtle influence of the distribution of the eigenvalues of the Laplacian matrix on variance amplification and discuss the roles of network size and topology. We provide concluding remarks in Section VII and technical details in appendices.

II Preliminaries and background

In this paper, we quantify the effect of stochastic uncertainties in gradient evaluation on the performance of first-order algorithms for unconstrained optimization problems

[TABLE]

where $f$: $\mathbb{R}^{n}\rightarrow\mathbb{R}$ is strongly convex with Lipschitz continuous gradient $\nabla f$. Specifically, we examine how gradient descent,

[TABLE]

amplify the additive white stochastic noise $w^{t}$ with zero mean and identity covariance matrix, $\operatorname{\mathbb{E}}\left[w^{t}\right]=0,$ $\operatorname{\mathbb{E}}\left[w^{t}(w^{\tau})^{T}\right]=I\,\delta(t-\tau)$. Here, $t$ is the iteration index, $x^{t}$ is the optimization variable, $\alpha$ is the stepsize, $\beta$ is an extrapolation parameter used for acceleration, $\sigma$ is the noise magnitude, $\delta$ is the Kronecker delta, and $\operatorname{\mathbb{E}}$ is the expected value. When the only source of uncertainty is a noisy gradient, we set $\sigma=\alpha$ in (2).

The set of functions $f$ that are $m$-strongly convex and $L$-smooth is denoted by $\mathcal{F}_{m}^{L}$; $f\in\mathcal{F}_{m}^{L}$ means that $f(x)-\frac{m}{2}\|x\|^{2}$ is convex and that the gradient $\nabla f$ is $L$-Lipschitz continuous. In particular, for a twice continuously differentiable function $f$ with the Hessian matrix $\nabla^{2}f$, we have

[TABLE]

In the absence of noise (i.e., for $\sigma=0$), for $f\in\mathcal{F}_{m}^{L}$, the parameters $\alpha$ and $\beta$ can be selected such that gradient descent and Nesterov’s accelerated method converge to the global minimum $x^{\star}$ of (1) with a linear rate $\rho<1$, i.e.,

[TABLE]

for all $t$ and some $c>0$. Table I provides the conventional values of these parameters and the corresponding guaranteed convergence rates [9]. Nesterov’s method with the parameters provided in Table I enjoys the convergence rate $\rho_{\mathrm{na}}=\sqrt{1-{1}/{\sqrt{\kappa}}}\leq 1-1/{(2\sqrt{\kappa})}$, where $\kappa\mathrel{\mathop{:}}=L/m$ is the condition number associated with $\mathcal{F}_{m}^{L}$. This rate is orderwise optimal in the sense that no first-order algorithm can optimize all $f\in\mathcal{F}_{m}^{L}$ with the rate $\rho_{\mathrm{lb}}=(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$ [9, Theorem 2.1.13]. Note that $1-\rho_{\mathrm{lb}}=O(1/\sqrt{\kappa})$ and $1-\rho_{\mathrm{na}}=\Omega(1/\sqrt{\kappa})$. In contrast to Nesterov’s method, the heavy-ball method does not offer any acceleration guarantees for all $f\in\mathcal{F}_{m}^{L}$. However, for strongly convex quadratic $f$, parameters can be selected to guarantee linear convergence of the heavy-ball method with a rate that outperforms the one achieved by Nesterov’s method [49]; see Table II.

To provide a quantitative characterization for the robustness of algorithms (2) to the noise $w^{t}$, we examine the performance measure,

[TABLE]

For quadratic objective functions, algorithms (2) are linear dynamical systems. In this case, $J$ quantifies the steady-state variance amplification and it can be computed from the solution of the algebraic Lyapunov equation; see Section III. For general strongly convex problems, there is no explicit characterization for $J$ but techniques from control theory can be utilized to compute an upper bound; see Section IV.

Notation

We write $g=\Omega(h)$ (or, equivalently, $h=O(g)$) to denote the existence of positive constants $c_{i}$ such that, for any $x>c_{2}$, the functions $g$ and $h$: $\mathbb{R}\rightarrow\mathbb{R}$ satisfy $g(x)\geq c_{1}h(x)$. We write $g=\Theta(h)$, or more informally $g\approx h$, if both $g=\Omega(h)$ and $g=O(h)$.

III Strongly convex quadratic problems

Consider a strongly convex quadratic objective function,

[TABLE]

where $Q$ is a symmetric positive definite matrix and $q$ is a vector. Let $f\in\mathcal{F}_{m}^{L}$ and let the eigenvalues $\lambda_{i}$ of $Q$ satisfy

[TABLE]

In the absence of noise, the constant values of parameters $\alpha$ and $\beta$ provided in Table II yield linear convergence (with optimal decay rates) to the globally optimal point $x^{\star}=Q^{-1}q$ for all three algorithms [49]. In the presence of additive white noise $w^{t}$, we derive analytical expressions for the variance amplification $J$ of algorithms (2) and demonstrate that $J$ depends not only on the algorithmic parameters $\alpha$ and $\beta$ but also on all eigenvalues of the Hessian matrix $Q$. This should be compared and contrasted to the optimal rate of linear convergence which only depends on $\kappa\mathrel{\mathop{:}}=L/m$, i.e., the ratio of the largest and smallest eigenvalues of $Q$.

For constant $\alpha$ and $\beta$, algorithms (2) can be described by a linear time-invariant (LTI) first-order recursion

[TABLE]

where $\psi^{t}$ is the state, $z^{t}\mathrel{\mathop{:}}=x^{t}-x^{\star}$ is the performance output, and $w^{t}$ is a white stochastic input. In particular, choosing $\psi^{t}\mathrel{\mathop{:}}=x^{t}-x^{\star}$ for gradient descent and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi^{t}}\mathrel{\mathop{:}}=[\,(x^{t}-x^{\star})^{T}\;(x^{t+1}-x^{\star})^{T}\,]^{T} for accelerated algorithms yields state-space model (5) with

[TABLE]

for gradient descent and

[TABLE]

for the heavy-ball and Nesterov’s methods, respectively, with

[TABLE]

Since $w^{t}$ is zero mean, we have $\operatorname{\mathbb{E}}\left(\psi^{t+1}\right)=A\operatorname{\mathbb{E}}\left(\psi^{t}\right)$. Thus, $\operatorname{\mathbb{E}}\left(\psi^{t}\right)=A^{t}\operatorname{\mathbb{E}}\left(\psi^{0}\right)$ and, for any stabilizing parameters $\alpha$ and $\beta$, $\lim_{t\,\to\,\infty}\operatorname{\mathbb{E}}\left(\psi^{t}\right)=0$, with the same linear rate as in the absence of noise. Furthermore, it is well-known that the covariance matrix $P^{t}\mathrel{\mathop{:}}=\operatorname{\mathbb{E}}\left(\psi^{t}(\psi^{t})^{T}\right)$ of the state vector satisfies the linear recursion

[TABLE]

where $Z=CPC^{T}$ is the steady-state limit of the output covariance matrix $Z^{t}\mathrel{\mathop{:}}=\operatorname{\mathbb{E}}\left(z^{t}(z^{t})^{T}\right)=CP^{t}C^{T}.$

We next provide analytical solution $P$ to (6c) that depends on the parameters $\alpha$ and $\beta$ as well as on the spectrum of the Hessian matrix $Q$. This allows us to explicitly characterize the variance amplification $J$ and quantify the impact of additive white noise on the performance of first-order optimization algorithms.

III-A Influence of the eigenvalues of the Hessian matrix

We use the modal decomposition of the symmetric matrix $Q=V\Lambda V^{T}$ to bring $A$, $B$, and $C$ in (5) into a block diagonal form, $\hat{A}=\mathrm{diag}\,(\hat{A}_{i}),$ $\hat{B}=\mathrm{diag}\,(\hat{B}_{i}),$ $\hat{C}=\mathrm{diag}\,(\hat{C}_{i}),$ with $i=1,\ldots,n.$ Here, $\Lambda=\mathrm{diag}\,(\lambda_{i})$ is the diagonal matrix of the eigenvalues and $V$ is the orthogonal matrix of the eigenvectors of $Q$. More specifically, the unitary coordinate transformation

[TABLE]

brings the state-space model of gradient descent into a diagonal form with

[TABLE]

This block diagonal structure allows us to explicitly solve Lyapunov equation (6c) for $P$ and derive an analytical expression for $J$ in terms of the eigenvalues $\lambda_{i}$ of the Hessian matrix $Q$ and the algorithmic parameters $\alpha$ and $\beta$. Namely, under coordinate transformation (7) and a suitable permutation of variables, equation (6c) can be brought into an equivalent set of equations,

[TABLE]

where $\hat{P}_{i}$ is a scalar for the gradient descent method and a $2\times 2$ matrix for the accelerated algorithms. In Theorem 1, we use the solution to these decoupled Lyapunov equations to express the variance amplification as

[TABLE]

where $\hat{J}(\lambda_{i})$ determines the contribution of the eigenvalue $\lambda_{i}$ of the matrix $Q$ to the variance amplification. In what follows, we use subscripts $\mathrm{gd}$, $\mathrm{hb}$, and $\mathrm{na}$ (e.g., $J_{\mathrm{gd}}$, $J_{\mathrm{hb}}$, and $J_{\mathrm{na}}$) to denote quantities that correspond to gradient descent (2a), heavy-ball method (2b), and Nesterov’s accelerated method (2c).

Theorem 1

For strongly convex quadratic problems, the variance amplification of noisy first-order algorithms (2) with any constant stabilizing parameters $\alpha$ and $\beta$ is determined by $J=\sum_{i\,=\,1}^{n}\hat{J}(\lambda_{i})$, where $\lambda_{i}$ is the $i$th eigenvalue of $Q=Q^{T}\succ 0$ and the modal contribution to the variance amplification $\hat{J}(\lambda)$ is given by

[TABLE]

Proof:

See Appendix -A. ∎

For strongly convex quadratic problems, Theorem 1 provides exact expressions for variance amplification of the first-order algorithms. These expressions not only quantify the dependence of $J$ on the algorithmic parameters $\alpha$ and $\beta$ and the impact of the largest and smallest eigenvalues, but also capture the effect of all other eigenvalues of the Hessian matrix $Q$. We also observe that the variance amplification $J$ is proportional to $\sigma^{2}$. Apart from Section V, where we examine the role of parameters $\alpha$ and $\beta$ on acceleration/robustness tradeoff and allow the dependence of $\sigma$ on $\alpha$, without loss of generality we choose $\sigma=1$ in the rest of the paper.

Remark 1

The performance measure $J$ in (6d) quantifies the steady-state variance of the iterates of first-order algorithms. Robustness of noisy algorithms can be also evaluated using alternative performance measures, e.g., the mean value of the error in the objective function [45],

[TABLE]

This measure of variance amplification can be characterized using our approach by defining $C=Q^{1/2}$ for gradient descent and $C=[\,Q^{1/2}\;\,0\,]$ for accelerated algorithms in state-space model (5). Furthermore, repeating the above procedure for the modified performance output $z^{t}$ yields $J^{\prime}=\sum_{i\,=\,1}^{n}\lambda_{i}\hat{J}(\lambda_{i}),$ where the respective expressions for $\hat{J}(\lambda_{i})$ are given in Theorem 1.

III-B Comparison for the parameters that optimize the convergence rate

We next examine the robustness of first-order algorithms applied to strongly convex quadratic problems for the parameters that optimize the linear convergence rate; see Table II. For these parameters, the eigenvalues of the matrix $A$ are inside the open unit disk, implying exponential stability of system (5). We first use the expressions presented in Theorem 1 to compare the variance amplification of the heavy-ball method to gradient descent.

Theorem 2

Let the strongly convex quadratic objective function $f$ in (4) satisfy $\lambda_{\max}(Q)=L$, $\lambda_{\min}(Q)=m>0$, and let $\kappa\mathrel{\mathop{:}}=L/m$ be the condition number. For the optimal parameters provided in Table II, the ratio between the variance amplification of the heavy-ball method and gradient descent with equal values of $\sigma$ is given by

[TABLE]

Proof:

For the parameters provided in Table II we have $\alpha_{\mathrm{hb}}=(1+\beta)\alpha_{\mathrm{gd}}$, where $\beta=(\sqrt{\kappa}-1)^{2}/(\sqrt{\kappa}+1)^{2}$ is the momentum parameter for the heavy-ball method. It is now straightforward to show that the modal contributions $\hat{J}_{\mathrm{hb}}$ and $\hat{J}_{\mathrm{gd}}$ to the variance amplification of the iterates given in Theorem 1 satisfy

[TABLE]

Thus, the ratio $\hat{J}_{\mathrm{hb}}(\lambda)/\hat{J}_{\mathrm{gd}}(\lambda)$ does not depend on $\lambda$ and is only a function of the condition number $\kappa$. Substitution of (12) into $J=\sum_{i}\hat{J}(\lambda_{i})$ yields relation (11). ∎

Theorem 2 establishes the linear relation between the variance amplification of the heavy-ball algorithm $J_{\mathrm{hb}}$ and the gradient descent $J_{\mathrm{gd}}$. We observe that the ratio $J_{\mathrm{hb}}/J_{\mathrm{gd}}$ only depends on the condition number $\kappa$ and that acceleration increases variance amplification: for $\kappa\gg 1$, $J_{\mathrm{hb}}$ is larger than $J_{\mathrm{gd}}$ by a factor of $\sqrt{\kappa}$. We next study the ratio between the variance amplification of Nesterov’s accelerated method and gradient descent. In contrast to the heavy-ball method, this ratio depends on the entire spectrum of the Hessian matrix $Q$. The following proposition, which examines the modal contributions $\hat{J}_{\mathrm{na}}(\lambda)$ and $\hat{J}_{\mathrm{gd}}(\lambda)$ of Nesterov’s accelerated method and gradient descent, is the key technical result that allows us to establish the largest and smallest values that the ratio $J_{\mathrm{na}}/J_{\mathrm{gd}}$ can take for a given pair of extreme eigenvalues $m$ and $L$ of $Q$ in Theorem 3.

Proposition 1

Let the strongly convex quadratic objective function $f$ in (4) satisfy $\lambda_{\max}(Q)=L$, $\lambda_{\min}(Q)=m>0$, and let $\kappa\mathrel{\mathop{:}}=L/m$ be the condition number. For the optimal parameters provided in Table II, the ratio $\hat{J}_{\mathrm{na}}(\lambda)/\hat{J}_{\mathrm{gd}}(\lambda)$ of modal contributions to variance amplification of Nesterov’s method and gradient descent is a decreasing function of $\lambda\in[m,L]$. Furthermore, for $\sigma=1$, the function $\hat{J}_{\mathrm{gd}}(\lambda)$ satisfies

[TABLE]

where $\bar{\kappa}\mathrel{\mathop{:}}=3\kappa+1$.

Proof:

See Appendix -A. ∎

For all three algorithms, Proposition 1 and Theorem 2 demonstrate that the modal contribution to the variance amplification of the iterates at the extreme eigenvalues of the Hessian matrix $m$ and $L$ only depends on the condition number $\kappa\mathrel{\mathop{:}}=L/m$. For gradient descent and the heavy-ball method, $\hat{J}$ achieves its largest value at $m$ and $L$, i.e.,

[TABLE]

Remark 2

Theorem 1 provides explicit formulas for variance amplification of noisy algorithms (2) in terms of the eigenvalues $\lambda_{i}$ of the Hessian matrix $Q$. Similarly, we can represent the variance amplification in terms of the eigenvalues $\hat{\lambda}_{i}$ of the dynamic matrices $\hat{A}_{i}$ in (8). For gradient descent, $\hat{\lambda}_{i}=1-\alpha\lambda_{i}$ and it is straightforward to verify that $J_{\mathrm{gd}}$ is determined by the sum of reciprocals of distances of these eigenvalues to the stability boundary, J_{\mathrm{gd}}=\sum_{i\,=\,1}^{n}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}/(1-\hat{\lambda}_{i}^{2}). Similarly, for accelerated methods we have,

[TABLE]

where $\hat{\lambda}_{i}$ and $\hat{\lambda}^{\prime}_{i}$ are the eigenvalues of $\hat{A}_{i}$. For Nesterov’s method with the parameters provided in Table II, the matrix $\hat{A}_{n}$, which corresponds to $\lambda_{n}=m$, admits a Jordan canonical form with repeated eigenvalues $\hat{\lambda}_{n}=\hat{\lambda}^{\prime}_{n}=1-2/\sqrt{3\kappa+1}$. In this case, $\hat{J}_{\mathrm{na}}(m)=\sigma^{2}(1+\hat{\lambda}_{n}^{2})/(1-\hat{\lambda}_{n}^{2})^{3}$, which should be compared and contrasted to the above expression for gradient descent. Furthermore, for both $\lambda_{1}=L$ and $\lambda_{n}=m$, the matrices $\hat{A}_{1}$ and $\hat{A}_{n}$ for the heavy-ball method with the parameters provided in Table II have eigenvalues with algebraic multiplicity two and incomplete sets of eigenvectors.

We next establish the range of values that the ratio $J_{\mathrm{na}}/J_{\mathrm{gd}}$ can take.

Theorem 3

For the strongly convex quadratic objective function $f$ in (4) with $x\in\mathbb{R}^{n}$, $\lambda_{\max}(Q)=L$, and $\lambda_{\min}(Q)=m>0$, the ratio between the variance amplification of Nesterov’s accelerated method and gradient descent, for the optimal parameters provided in Table II and equal values of $\sigma$ satisfies

[TABLE]

Proof:

See Appendix -A. ∎

Theorem 3 provides tight upper and lower bounds on the ratio between $J_{\mathrm{na}}$ and $J_{\mathrm{gd}}$ for strongly convex quadratic problems. As shown in Appendix -A, the lower bound is achieved for a quadratic function in which the Hessian matrix $Q$ has one eigenvalue at $m$ and $n-1$ eigenvalues at $L$, and the upper bound is achieved when $Q$ has one eigenvalue at $L$ and the remaining ones at $m$. Theorem 3 in conjunction with Proposition 1 demonstrate that for a fixed problem dimension $n$, $J_{\mathrm{na}}$ is larger than $J_{\mathrm{gd}}$ by a factor of $\sqrt{\kappa}$ for $\kappa\gg 1$.

This trade-off is further highlighted in Theorem 4 which provides tight bounds on the variance amplification of iterates in terms of the problem dimension $n$ and the condition number $\kappa$ for all three algorithms. To simplify the presentation, we first use the explicit expressions for $\hat{J}_{\mathrm{na}}(m)$ and $\hat{J}_{\mathrm{na}}(L)$ in Proposition 1 to obtain the following upper and lower bounds on $\hat{J}_{\mathrm{na}}(m)$ and $\hat{J}_{\mathrm{na}}(L)$ (see Appendix -A)

[TABLE]

Theorem 4

For the strongly convex quadratic objective function $f$ in (4) with $x\in\mathbb{R}^{n}$, $\lambda_{\max}(Q)=L$, $\lambda_{\min}(Q)=m>0$, and $\kappa\mathrel{\mathop{:}}=L/m$, the variance amplification of the first-order optimization algorithms, with the parameters provided in Table II and $\sigma=1$, is bounded by

[TABLE]

Proof:

As shown in Proposition 1, the functions $\hat{J}(\lambda)$ for gradient descent and Nesterov’s algorithm attain their largest and smallest values over the interval $[m,L]$ at $\lambda=m$ and $\lambda={1}/{\alpha}$, respectively. Thus, fixing the smallest and largest eigenvalues, the variance amplification $J$ is maximized when the other $n-2$ eigenvalues are all equal to $m$ and is minimized when they are all equal to $1/\alpha$. This combined with the explicit expressions for $\hat{J}_{\mathrm{gd}}(m)$, $\hat{J}_{\mathrm{gd}}(L)$, and $\hat{J}_{\mathrm{gd}}(1/\alpha)$ in (13a) leads to the tight upper and lower bounds for gradient descent. For the heavy-ball method, the bounds follow from Theorem 2 and for Nesterov’s algorithm, the bounds follow from (16). ∎

For problems with a fixed dimension $n$ and a condition number $\kappa\gg n$, there is an $\Omega(\sqrt{\kappa})$ difference in both upper and lower bounds provided in Theorem 4 for the accelerated algorithms relative to gradient descent. Even though Theorem 4 considers only the values of $\alpha$ and $\beta$ that optimize the convergence rate, in Section V we demonstrate that this gap is fundamental in that it holds for any parameters that yield an accelerated convergence rate. It is worth noting that both the lower and upper bounds are influenced by the problem dimension $n$ and the condition number $\kappa$. For large-scale problems, there may be a subtle relation between $n$ and $\kappa$ and the established bounds may exhibit different scaling trends. In Section VI, we identify a class of quadratic optimization problems for which $J_{\mathrm{na}}$ scales in the same way as $J_{\mathrm{gd}}$ for $\kappa\gg 1$ and $n\gg 1$.

Before we elaborate further on these issues, we provide two illustrative examples that highlight the importance of the choice of the performance metric in the robustness analysis of noisy algorithms. It is worth noting that an $O(\kappa)$ upper bound for gradient descent and an $O(\kappa^{2})$ upper bound for Nesterov’s accelerated algorithm was established in [24]. Relative to this upper bound for Nesterov’s method, the upper bound provided in Theorem 4 is tighter by a factor of $\sqrt{\kappa}$. Theorem 4 also provides lower bounds, reveals the influence of the problem dimension $n$, and identifies constants that multiply the leading terms in the condition number $\kappa$. Moreover, in Section IV we demonstrate that similar upper bounds can be obtained for general strongly convex objective functions with Lipschitz continuous gradients.

III-C Examples

We next provide illustrative examples to (i) demonstrate the agreement of our theoretical predictions with the results of stochastic simulations; and (ii) contrast two natural performance measures, namely the variance of the iterates $J$ in (6d) and the mean objective error $J^{\prime}$ in (10), for assessing robustness of noisy optimization algorithms.

Example 1

Let us consider the quadratic objective function in (4) with

[TABLE]

For all three algorithms, the performance measures $J$ and $J^{\prime}$ are given by

[TABLE]

As shown in (14), $\hat{J}(m)$ and $\hat{J}(L)$ only depend on the condition number $\kappa$ and the variance amplification of the iterates satisfies

[TABLE]

• •

for $m\ll 1$ and $L=O(1)$

[TABLE]

• •

for $L\gg 1$ and $m=O(1)$

[TABLE]

Relation (18a) reveals the detrimental impact of acceleration on the variance of the optimization variable. On the other hand, (18b) and (18c) show that, relative to gradient descent, the heavy-ball method increases the mean error in the objective function while Nesterov’s method reduces it. Thus, if the mean value of the error in the objective function is to be used to assess performance of noisy algorithms, one can conclude that Nesterov’s method significantly outperforms gradient descent both in terms of convergence rate and robustness to noise. However, this performance metric fails to capture large variance of the mode associated with the smallest eigenvalue of the matrix $Q$ in Nesterov’s algorithm. Theorem 2 and Proposition 1 show that the modal contributions to the variance amplification of the iterates for gradient descent and the heavy-ball method are balanced at $m$ and $L$, i.e., $\hat{J}_{\mathrm{gd}}(m)=\hat{J}_{\mathrm{gd}}(L)=\Theta(\kappa)$ and $\hat{J}_{\mathrm{hb}}(m)=\hat{J}_{\mathrm{hb}}(L)=\Theta(\kappa\sqrt{\kappa})$. On the other hand, for Nesterov’s method there is a $\Theta(\kappa)$ gap between $\hat{J}_{\mathrm{na}}(m)=\Theta(\kappa\sqrt{\kappa})$ and $\hat{J}_{\mathrm{na}}(L)=\Theta(\sqrt{\kappa})$. While the performance measure $J^{\prime}$ reveals a superior performance of Nesterov’s algorithm at large condition numbers, it fails to capture the negative impact of acceleration on the variance of the optimization variable; see Fig. 1 for an illustration.

Figure 2 shows the performance outputs $z^{t}=x^{t}$ and $z^{t}=Q^{1/2}x^{t}$ resulting from $10^{5}$ iterations of noisy first-order algorithms with the optimal parameters provided in Table II for the strongly convex objective function $f(x)=0.5\,x_{1}^{2}+0.25\times 10^{-4}\,x_{2}^{2}$ ($\kappa=2\times 10^{4}$). Although Nesterov’s method exhibits good performance with respect to the error in the objective function (performance measure $J^{\prime}$), the plots in the first row illustrate detrimental impact of noise on both accelerated algorithms with respect to the variance of the iterates (performance measure $J$). In particular, we observe that: (i) for gradient descent and the heavy-ball method, the iterates $x^{t}$ are scattered uniformly along the eigen-directions of the Hessian matrix $Q$ and acceleration increases variance equally along all directions; and (ii) relative to gradient descent, Nesterov’s method exhibits larger variance in the iterates $x^{t}$ along the direction that corresponds to the smallest eigenvalue $\lambda_{\min}(Q)$.

Example 2

Figure 3 compares the results of twenty stochastic simulations for a strongly convex quadratic objective function (4) with $q=0$ and a Toeplitz matrix $Q\in\mathbb{R}^{50\times 50}$ with the first row $[\,2\,-\!1~{}\,0\;\cdots\;0\;\;0\,]^{T}$. This figure shows the time-dependence of the variance of the performance outputs $z^{t}=x^{t}$ and $z^{t}=Q^{1/2}x^{t}$ for the algorithms subject to additive white noise with zero initial conditions. The plots further demonstrate that the mean error in the objective function does not capture detrimental impact of noise on the variance of the iterates for Nesterov’s algorithm. The bottom row also compares variance obtained by averaging outcomes of twenty stochastic simulations with the corresponding theoretical values resulting from the Lyapunov equations.

IV General strongly convex problems

In this section, we extend our results to the class $\mathcal{F}_{m}^{L}$ of $m$-strongly convex objective functions with $L$-Lipschitz continuous gradients. While a precise characterization of noise amplification for general problems is challenging because of the nonlinear dynamics, we employ tools from robust control theory to obtain meaningful upper bounds. Our results utilize the theory of integral quadratic constraints [50], a convex control-theoretic framework that was recently used to analyze optimization algorithms [49] and study convergence and robustness of the first-order methods [51, 52, 53, 54]. We establish analytical upper bounds on the mean-squared error of the iterates (3) for gradient descent (2a) and Nesterov’s accelerated (2c) methods. Since there are no known accelerated convergence guarantees for the heavy-ball method when applied to general strongly convex functions, we do not consider it in this section.

We first exploit structural properties of the gradient and employ quadratic Lyapunov functions to formulate a semidefinite programing problem (SDP) that provides upper bounds on $J$ in (3). While quadratic Lyapunov functions yield tight upper bounds for gradient descent, they fail to provide any upper bound for Nesterov’s method for large condition numbers ($\kappa>100$). To overcome this challenge, we present a modified semidefinite program that uses more general Lyapunov functions which are obtained by augmenting standard quadratic terms with the objective function. This type of generalized Lyapunov functions has been introduced in [55, 52] and used to study convergence of optimization algorithms for non-strongly convex problems. We employ a modified SDP to derive meaningful upper bounds on $J$ in (3) for Nesterov’s method as well.

We note that algorithms (2) are invariant under translation, i.e., if we let $\tilde{x}\mathrel{\mathop{:}}=x-\bar{x}$ and $g(\tilde{x})\mathrel{\mathop{:}}=f(\tilde{x}+\bar{x})$, then (2c), for example, satisfies

[TABLE]

Thus, in what follows, without loss of generality, we assume that $x^{\star}=0$ is the unique minimizer of (1).

IV-A An approach based on contraction mappings

Before we present our approach based on Linear Matrix Inequalities (LMIs), we provide a more intuitive approach that can be used to examine noise amplification of gradient descent. Let \varphi$$:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n} be a contraction mapping, i.e., there exists a positive scalar $\eta<1$ such that $\|\varphi(x)-\varphi(y)\|\leq\eta\|x-y\|$ for all $x$, $y\in\mathbb{R}^{n}$, and let $x^{\star}=0$ be the unique fixed point of $\varphi$, i.e, $\varphi(0)=0$. For the noisy recursion $x^{t+1}=\varphi(x^{t})+\sigma w^{t},$ where $w^{t}$ is a zero-mean white noise with identity covariance and $\operatorname{\mathbb{E}}((w^{t})^{T}\varphi(x^{t}))=0$, the contractiveness of $\varphi$ implies

[TABLE]

Since $\eta<1$, this relation yields

[TABLE]

If $\eta\mathrel{\mathop{:}}=\max\{|1-\alpha m|,|1-\alpha L|\}<1$, the map $\varphi(x)\mathrel{\mathop{:}}=x-\alpha\nabla f(x)$ is a contraction [56]. Thus, for the conventional stepsize $\alpha=1/L$ we have $\eta=1-1/\kappa$, and the bound becomes

[TABLE]

In the next section, we show that this upper bound is indeed tight for the class of functions $\mathcal{F}_{m}^{L}$. While this approach yields a tight upper bound for gradient descent, it cannot be used for Nesterov’s method (because it is not a contraction).

IV-B An approach based on linear matrix inequalities

For any function $f\in\mathcal{F}_{m}^{L}$, the nonlinear mapping \Delta$$:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}

[TABLE]

satisfies the quadratic inequality [49, Lemma 6]

[TABLE]

for all $y$, $y_{0}\in\mathbb{R}^{n}$, where the matrix $\Pi$ is given by

[TABLE]

We can bring algorithms (2) with constant parameters into a time-invariant state-space form

[TABLE]

In what follows, we demonstrate how property (19) of the nonlinear mapping $\Delta$ allows us to obtain upper bounds on $J$ when system (23a) is driven by the white stochastic input $w^{t}$ with zero mean and identity covariance. Lemma 1 uses a quadratic Lyapunov function of the form $V(\psi)=\psi^{T}X\psi$ and provides upper bounds on the steady-state second-order moment of the performance output $z^{t}$ in terms of solutions to a certain LMI. This approach yields a tight upper bound for gradient descent.

Lemma 1

Let the nonlinear function $u=\Delta(y)$ satisfy the quadratic inequality

[TABLE]

for some matrix $\Pi$, let $X$ be a positive semidefinite matrix, and let $\lambda$ be a nonnegative scalar such that system (23a) satisfies

[TABLE]

Then the steady-state second-order moment $J$ of the performance output $z^{t}$ in (23a) is bounded by

[TABLE]

Proof:

See Appendix -B. ∎

For Nesterov’s accelerated method with the parameters provided in Table I, computational experiments show that LMI (25) becomes infeasible for large values of the condition number $\kappa$. Thus, Lemma 1 does not provide sensible upper bounds on $J$ for Nesterov’s algorithm. This observation is consistent with the results of [49], where it was suggested that analyzing the convergence rate requires the use of additional quadratic inequalities, apart from (19), to further tighten the constraints on the gradient $\nabla f$ and reduce conservativeness. In what follows, we build on the results of [52] and present an alternative LMI in Lemma 2 that is obtained using a Lyapunov function of the form V(\psi)$$=\psi^{T}X\psi+f([\,0~{}\,I\,]\psi), where $X$ is a positive semidefinite matrix and $f$ is the objective function in (1). Such Lyapunov functions have been used to study convergence of optimization algorithms in [55]. The resulting approach allows us to establish an order-wise tight analytical upper bound on $J$ for Nesterov’s accelerated method.

Lemma 2

Let the matrix $M(m,L;\alpha,\beta)$ be defined as

[TABLE]

where

[TABLE]

Consider state-space model (23a)-(LABEL:eq.ss-na) for algorithm (2c) and let $\Pi$ be given by (22). Then, for any positive semidefinite matrix $X$ and scalars $\lambda_{1}\geq 0$ and $\lambda_{2}\geq 0$ that satisfy

[TABLE]

the steady-state second-order moment $J$ of the performance output $z^{t}$ in (23a) is bounded by

[TABLE]

Proof:

See Appendix -B. ∎

Remark 3

Since LMI (26) simplifies to (25) by setting $\lambda_{2}=0$, Lemma 2 represents a relaxed version of Lemma 1. This modification is the key enabler to establishing tight upper bound on $J$ for Nesterov’s method.

The upper bounds provided in Lemmas 1 and 2 are proportional to $\sigma^{2}$. In what follows, to make a connection between these bounds and our analytical expressions for the variance amplification in the quadratic case (Section III), we again set $\sigma=1$. The best upper bound on $J$ that can be obtained using Lemma 2 is given by the optimal objective value of the semidefinite program

[TABLE]

For system matrices (LABEL:eq.ss-na), LMI (26) is of size $3n\times 3n$ where $x^{t}\in\mathbb{R}^{n}$. However, if we impose the additional constraint that the matrix $X$ has the same block structure as $A$,

[TABLE]

for some scalars $x_{1}$, $x_{2}$, and $x_{0}$, then using appropriate permutation matrices, we can simplify (25) into an LMI of size $3\times 3$. Furthermore, imposing this constraint comes without loss of generality. In particular, the optimal objective value of problem (28) does not change if we require $X$ to have this structure; see [49, Section 4.2] for a discussion of this lossless dimensionality reduction for LMI constraints with similar structure.

In Theorem 5, we use Lemmas 1 and 2 to establish tight upper bounds on $J_{\mathrm{gd}}$ and $J_{\mathrm{na}}$ for all $f\in\mathcal{F}_{m}^{L}$.

Theorem 5

For gradient descent and Nesterov’s accelerated method with the parameters provided in Table I and $\sigma=1$, the performance measures $J_{\mathrm{gd}}$ and $J_{\mathrm{na}}$ of the error $x^{t}-x^{\star}\in\mathbb{R}^{n}$ satisfy

[TABLE]

where

[TABLE]

and $\kappa\mathrel{\mathop{:}}=L/m$ is the condition number of the set $\mathcal{F}_{m}^{L}$.

Proof:

See Appendix -B. ∎

The variance amplification of gradient descent and Nesterov’s method for $f(x)=\tfrac{m}{2}\,x^{T}x$ in $\mathcal{F}_{m}^{L}$ is determined by $q_{\mathrm{gd}}$ and $q_{\mathrm{na}}$, respectively, and these two quantities can be obtained using Theorem 1. In Theorem 5, we use this strongly convex quadratic objective function to certify the accuracy of the upper bounds on $\sup J$ for all $f\in\mathcal{F}_{m}^{L}$. In particular, we observe that the upper bound is exact for gradient descent and that it is within a $4.08$ factor of the optimal for Nesterov’s method.

For strongly convex objective functions with the condition number $\kappa$, Theorem 5 proves that gradient descent outperforms Nesterov’s accelerated method in terms of the largest noise amplification by a factor of $\sqrt{\kappa}$. This uncovers the fundamental performance limitation of Nesterov’s accelerated method when the gradient evaluation is subject to additive stochastic uncertainties.

V Tuning of algorithmic parameters

The parameters provided in Table II yield the optimal convergence rate for strongly convex quadratic problems. For these specific values, Theorem 4 establishes upper and lower bounds on the variance amplification that reveal the negative impact of acceleration. However, it is relevant to examine whether the parameters can be designed to provide acceleration while reducing the variance amplification.

While the convergence rate solely depends on the extreme eigenvalues $m=\lambda_{\min}(Q)$ and $L=\lambda_{\max}(Q)$ of the Hessian matrix $Q$, variance amplification is influenced by the entire spectrum of $Q$ and its minimization is challenging as it requires the use of all eigenvalues. In this section, we first consider the special case of eigenvalues being symmetrically distributed over the interval $[m,L]$ and demonstrate that for gradient descent and the heavy-ball method, the parameters provided in Table II yield a variance amplification that is within a constant factor of the optimal value. As we demonstrate in Section VI, symmetric distribution of the eigenvalues is encountered in distributed consensus over undirected torus networks. We also consider the problem of designing parameters for objective functions in which the problem size satisfies $n\ll\kappa$ and establish a trade-off between convergence rate and variance amplification. More specifically, we show that for any accelerating pair of parameters $\alpha$ and $\beta$ and bounded problem dimension $n$, the variance amplification of accelerated methods is larger than that of gradient descent by a factor of $\Omega(\sqrt{\kappa})$.

V-A Tuning of parameters using the whole spectrum

Let $L=\lambda_{1}\geq\lambda_{2}\,\geq\cdots\geq\lambda_{n}=m>0$ be the eigenvalues of the Hessian matrix $Q$ of the strongly convex quadratic objective function in (4). Algorithms (2) converge linearly in expected value to the optimizer $x^{\star}$ with the rate

[TABLE]

where $\hat{\rho}(\lambda_{i})$ is the spectral radius of the matrix $\hat{A}_{i}$ given by (8). For any scalar $c>0$ and fixed $\sigma$, let

[TABLE]

for gradient descent, where the expression for the variance amplification $J$ is provided in Theorem 1. Here, the constraints enforce a standard rate of linear convergence for gradient descent and an accelerated rate of linear convergence for the heavy-ball method parametrized with the constant $c$. Obtaining a closed form solution to (30) is challenging because $J$ depends on all eigenvalues of the Hessian matrix $Q$. Herein, we focus on objective functions for which the spectrum of $Q$ is symmetric, i.e., for any eigenvalue $\lambda$, the corresponding mirror image $\lambda^{\prime}:=L+m-\lambda$ with respect to $\frac{1}{2}(L+m)$ is also an eigenvalue with the same algebraic multiplicity. For this class of problems, Theorem 6 demonstrates that the parameters provided in Table II for gradient descent and the heavy-ball method yield variance amplification that is within a constant factor of the optimal.

Theorem 6

For any scalar $c>0$ and fixed $\sigma$, there exist constants $c_{1}\geq 1$ and $c_{2}>0$ such that for any strongly convex quadratic objective function in which the spectrum of the Hessian matrix $Q$ is symmetrically distributed over the interval $[m,L]$ with $\kappa\mathrel{\mathop{:}}=L/m>c_{1}$, we have

[TABLE]

where parameters $\alpha_{\mathrm{gd}}$ and ($\alpha_{\mathrm{hb}},\beta_{\mathrm{hb}}$) are provided in Table II, whereas $\alpha_{\mathrm{gd}}^{\star}(c)$ and ($\alpha^{\star}_{\mathrm{hb}}(c),\beta^{\star}_{\mathrm{hb}}(c)$) solve (30).

Proof:

See Appendix -C. ∎

For strongly convex quadratic objective functions with symmetric spectrum of the Hessian matrix over the interval $[m,L]$, Theorem 6 shows that the variance amplifications of gradient descent and the heavy-ball method with the parameters provided in Table II are within a constant factors of the optimal values. As we illustrate in Section VI, this class of problems is encountered in distributed averaging over noisy undirected networks. Combining this result with the lower bound on $J_{\mathrm{hb}}(\alpha_{\mathrm{hb}},\beta_{\mathrm{hb}})$ and the upper bound on $J_{\mathrm{gd}}(\alpha_{\mathrm{gd}})$ established in Theorem 4, we see that regardless of the choice of parameters, there is a fundamental gap of $\Omega(\sqrt{\kappa})$ between $J_{\mathrm{hb}}$ and $J_{\mathrm{gd}}$ as long as we require an accelerated rate of convergence.

V-B Fundamental lower bounds

We next establish lower bounds on the variance amplification of accelerated methods that hold for any pair of $\alpha$ and $\beta$ for strongly convex quadratic problems with $\kappa\gg 1$. In particular, we show that the variance amplification of accelerated algorithms is lower bounded by $\Omega(\kappa^{3/2})$ irrespective of the choice of $\alpha$ and $\beta$.

The next theorem establishes a fundamental tradeoff between the convergence rate and variance amplification for the heavy-ball method.

Theorem 7

For strongly convex quadratic problems with any stabilizing parameters $\alpha>0$ and $0<\beta<1$ and with a fixed noise magnitude $\sigma$, the heavy-ball method with the linear convergence rate $\rho$ satisfies

[TABLE]

Proof:

See Appendix -D. ∎

To gain additional insight, let us consider two special cases: (i) for $\alpha=1/L$ and $\beta\rightarrow 0^{+}$, we obtain gradient descent algorithm for which $1-\rho=\Theta({1}/{\kappa})$ and $J=\Theta(\kappa)$; (ii) for the heavy-ball method with the parameters provided in Table II, we have $1-\rho=\Theta({1}/{\sqrt{\kappa}})$ and $J=\Theta(\kappa\sqrt{\kappa})$. Thus, in both cases, ${J_{\mathrm{hb}}}/(1-\rho)=\Omega(\kappa^{2})$. Theorem 7 shows that this lower bound is fundamental and it therefore quantifies the tradeoff between the convergence rate and the variance amplification of the heavy-ball method for any choice of parameters $\alpha$ and $\beta$. It is also worth noting that the lower bound for $\sigma=\alpha$ depends on the largest eigenvalue $L$ of the Hessian matrix $Q$. Thus, this bound is meaningful when the value of $L$ is uniformly upper bounded. This scenario occurs in many applications including consensus over undirected tori networks; see Section VI.

While we are not able to show a similar lower bound for Nesterov’s method, in the next theorem, we establish an asymptotic lower bound on the variance amplification that holds for any pair of accelerating parameters ($\alpha,\beta$) for both Nesterov’s and heavy-ball methods.

Theorem 8

For a strongly convex quadratic objective function with condition number $\kappa$, let $c>0$ be a constant such that either Nesterov’s algorithm or the heavy-ball method with some (possibly problem dependent) parameters $\alpha>0$ and $0<\beta<1$ converges linearly with a rate $\rho\leq 1-c/\sqrt{\kappa}$. Then, for any fixed noise magnitude $\sigma$, the variance amplification satisfies

[TABLE]

Proof:

For the heavy-ball method, the result follows from combining Theorem 7 with the inequality $1-\rho\geq{c}/{\sqrt{\kappa}}$. For Nesterov’s method, the proof is provided in Appendix -D. ∎

For problems with $n\ll\kappa$, we recall that the variance amplification of gradient descent with conventional values of parameters scales as $O(\kappa)$; see Theorem 5. Irrespective of the choice of parameters $\alpha$ and $\beta$, this result in conjunction with Theorem 8 demonstrates that acceleration cannot be achieved without increasing the variance amplification $J$ by a factor of $\Omega(\sqrt{\kappa})$.

VI Application to distributed computation over undirected networks

Distributed computation over networks has received significant attention in optimization, control systems, signal processing, communications, and machine learning communities. In this problem, the goal is to optimize an objective function (e.g., for the purpose of training a model) using multiple processing units that are connected over a network. Clearly, the structure of the network (e.g., node dynamics and network topology) may impact the performance (e.g., convergence rate and noise amplification) of any optimization algorithm. As a first step toward understanding the impact of the network structure on performance of noisy first-order optimization algorithms, in this section, we examine the standard distributed consensus problem.

The consensus problem arises in applications ranging from social networks, to distributed computing networks, to cooperative control in multi-agent systems. In the simplest setup, each node updates a scalar value using the values of its neighbors such that they all agree on a single consensus value. Simple updating strategies of this kind can be obtained by applying a first-order algorithm to the convex quadratic problem

[TABLE]

where $\mathbf{L}=\mathbf{L}^{T}\in\mathbb{R}^{n\times n}$ is the Laplacian matrix of the graph associated with the underlying undirected network and $x\in\mathbb{R}^{n}$ is the vector of node values.

The graph Laplacian matrix $\mathbf{L}\succeq 0$ has a nontrivial null space that consists of the minimizers of problem (33). In the absence of noise, for gradient descent and both of its accelerated variants, it is straightforward to verify that the projections $v^{t}$ of the iterates $x^{t}$ onto the null space of $\mathbf{L}$ remain constant ($v^{t}=v^{0}$, for all $t$) and also that $x^{t}$ converges linearly to $v^{0}$. In the presence of additive noise, however, $v^{t}$ experiences a random walk which leads to an unbounded variance of $x^{t}$ as $t\to\infty$. Instead, as described in [37], the performance of algorithms in this case can be quantified by examining $\bar{J}\mathrel{\mathop{:}}=\lim_{t\,\to\,\infty}\operatorname{\mathbb{E}}\left(\|x^{t}-v^{t}\|^{2}\right).$ For connected networks, the null space of $\mathbf{L}$ is given by $\operatorname*{{\cal N}}(\mathbf{L})=\{c\mathds{1}\,|\,c\in\mathbb{R}\}$ and

[TABLE]

quantifies the mean-squared deviation from the network average, where $\mathds{1}$ denotes the vector of all ones, i.e., $\mathds{1}\mathrel{\mathop{:}}=[\,1\;\cdots\;1\,]^{T}$. Finally, it is straightforward to show that $\bar{J}$ can also be computed using the formulae in Theorem 1 by summing over the non-zero eigenvalues of $\mathbf{L}$.

In what follows, we consider a class of networks whose structure allows for the explicit evaluation of the eigenvalues of the Laplacian matrix $\mathbf{L}$. For $d$-dimensional torus networks, fundamental performance limitations of standard consensus algorithms in continuous time were established in [38], but it remains an open question if gradient descent and its accelerated variants suffer from these limitations. We utilize such torus networks to demonstrate that standard gradient descent exhibits the same scaling trends as consensus algorithms studied in [38] and that, in lower spatial dimensions, acceleration always increases variance amplification.

VI-A Explicit formulae for $d$-dimensional torus networks

We next examine the asymptotic scaling trends of the performance metric $\bar{J}$ given by (34) for large problem dimensions $n\gg 1$ and highlight the subtle influence of the distribution of the eigenvalues of $\mathbf{L}$ on the variance amplification for $d$-dimensional torus networks. Tori with nearest neighbor interactions generalize one-dimensional rings to higher spatial dimensions. Let $\mathbb{Z}_{{n_{0}}}$ denote the group of integers modulo ${n_{0}}$. A $d$-dimensional torus $\mathbb{T}^{d}_{n_{0}}$ consists of $n\mathrel{\mathop{:}}=n_{0}^{d}$ nodes denoted by $v_{a}$ where $a\in\mathbb{Z}_{n_{0}}^{d}$ and the set of edges $\{\{v_{a}~{}v_{b}\}\,|\,\|a-b\|=1\!\!\mod{n_{0}}\};$ nodes $v_{a}$ and $v_{b}$ are neighbors if and only if $a$ and $b$ differ exactly at a single entry by one. For example, $\mathbb{T}_{n_{0}}^{1}$ denotes a ring with $n={n_{0}}$ nodes and $\mathbb{T}_{n_{0}}^{5}$ denotes a five dimensional torus with $n=n_{0}^{5}$ nodes.

The multidimensional discrete Fourier transform can be used to determine the eigenvalues of the Laplacian matrix $\mathbf{L}$ of a $d$-dimensional torus $\mathbb{T}^{d}_{n_{0}}$,

[TABLE]

where $i\mathrel{\mathop{:}}=(i_{1},\ldots,i_{d})\in\mathbb{Z}_{n_{0}}^{d}$. We note that $\lambda_{0}=0$ is the only zero eigenvalue of $\mathbf{L}$ with the eigenvector $\mathds{1}$ and that all other eigenvalues are positive. Let $\kappa\mathrel{\mathop{:}}=\lambda_{\max}/\lambda_{\min}$ be the ratio of the largest and smallest nonzero eigenvalues of $\mathbf{L}$. A key observation is that, for ${n_{0}}\gg 1$,

[TABLE]

This is because $\lambda_{\min}=2d\,(1-\cos\,({2\pi}/{{n_{0}}}))$ goes to zero as ${n_{0}}\to\infty$, and the largest eigenvalue of $\mathbf{L}$, $\lambda_{\max}=2d\,(1-\cos\,(2\pi\lfloor\tfrac{{n_{0}}}{2}\rfloor/{{n_{0}}})),$ is equal to $4\,d$ for even ${n_{0}}$ and it approaches $4\,d$ from below for odd ${n_{0}}$.

As aforementioned, the performance metric $\bar{J}$ can be obtained by

[TABLE]

where $\hat{J}(\lambda)$ for each algorithm is determined in Theorem 1 and $\lambda_{i}$ are the non-zero eigenvalues of $\mathbf{L}$. The next theorem characterizes the asymptotic value of the network-size normalized mean-squared deviation from the network average, $\bar{J}/n$, for a fixed spatial dimension $d$ and condition number $\kappa\gg 1$. This result is obtained using analytical expression (35) for the eigenvalues of the Laplacian matrix $\mathbf{L}$.

Theorem 9

Let $\mathbf{L}\in\mathbb{R}^{n\times n}$ be the graph Laplacian of the $d$-dimensional undirected torus $\mathbb{T}_{n_{0}}^{d}$ with $n=n_{0}^{d}\gg 1$ nodes. For convex quadratic optimization problem (33), the network-size normalized performance metric $\bar{J}/n$ of noisy first-order algorithms with the parameters provided in Table II and $\sigma=1$, is determined by

[TABLE]

where $\kappa=\Theta(n^{2/d})$ is the condition number of $\mathbf{L}$ given in (36).

Proof:

See Appendix -E. ∎

Theorem 9 demonstrates that the variance amplification of gradient descent is equivalent to that of the standard consensus algorithm studied in [38] and that, in lower spatial dimensions, acceleration always negatively impacts the performance of noisy algorithms. Our results also highlight the subtle influence of the distribution of the eigenvalues of $\mathbf{L}$ on the variance amplification. For rings (i.e., $d=1$), lower bounds provided in Theorem 4 capture the trends that our detailed analysis based on the distribution of the entire spectrum of $\mathbf{L}$ reveals. In higher spatial dimensions, however, the lower bounds that are obtained using only the extreme eigenvalues of $\mathbf{L}$ are conservative. Similar conclusion can be made about the upper bounds provided in Theorem 4. This observation demonstrates that the naïve bounds that result only from the use of the extreme eigenvalues can be overly conservative.

We also note that gradient descent significantly outperforms Nesterov’s accelerated algorithm in lower spatial dimensions. In particular, while $\bar{J}/n$ becomes network-size-independent for $d=3$ for gradient descent, Nesterov’s algorithm reaches “critical connectivity” only for $d=5$. On the other hand, in any spatial dimension, there is no network-size independent upper bound on $\bar{J}/n$ for the heavy-ball method. These conclusions could not have been reached without performing an in-depth analysis of the impact of all eigenvalues on performance of noisy networks with $n\gg 1$ and $\kappa\gg 1$.

VII Concluding remarks

We study the robustness of noisy first-order algorithms for smooth, unconstrained, strongly convex optimization problems. Even though the underlying dynamics of these algorithms are in general nonlinear, we establish upper bounds on noise amplification that are accurate up to constant factors. For quadratic objective functions, we provide analytical expressions that quantify the effect of all eigenvalues of the Hessian matrix on variance amplification. We use these expressions to establish lower bounds demonstrating that although the acceleration techniques improve the convergence rate they significantly amplify noise for problems with large condition numbers. In problems of bounded dimension $n\ll\kappa$, the noise amplification increases from $O(\kappa)$ to $\Omega(\kappa^{{3}/{2}})$ when moving from standard gradient descent to accelerated algorithms. We specialize our results to the problem of distributed averaging over noisy undirected networks and also study the role of network size and topology on robustness of accelerated algorithms. Future research directions include (i) extension of our analysis to multiplicative and correlated noise; and (ii) robustness analysis of broader classes of optimization algorithms.

-A Quadratic problems

Proof:

For gradient descent, $\hat{A}_{i}=1-\alpha\lambda_{i}$ and $\hat{B}_{i}=1$ are scalars and the solution to (9) is given by

[TABLE]

For the accelerated methods, we note that for any $\hat{A}_{i}$ and $\hat{B}_{i}$ of the form

[TABLE]

the solution $\hat{P}_{i}$ to Lyapunov equation (9) is given by

[TABLE]

where

[TABLE]

The parameters $a_{i}$ and $b_{i}$ for Nesterov’s algorithm are $\{a_{i}=-\beta(1-\alpha\lambda_{i})$; $b_{i}=(1+\beta)(1-\alpha\lambda_{i})\}$ and for the heavy-ball method we have $\{a_{i}=-\beta$; $b_{i}=1+\beta-\alpha\lambda_{i}\}.$ Now, since $\hat{C_{i}}=1$ for gradient descent and $\hat{C}_{i}=[\,1\;\,0\,]$ for the accelerated algorithms, it follows that for all three algorithms we have \hat{J}(\lambda_{i})\mathrel{\mathop{:}}=\mathrm{trace}\,(\hat{C}_{i}\hat{P}_{i}\hat{C}_{i}^{T})={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}p_{i}. Finally, if we use the expression for $p_{i}$ for gradient descent and substitute for $a_{i}$ and $b_{i}$ in (37) for the accelerated algorithms, we obtain the expressions for $\hat{J}$ in the statement of the theorem. ∎

Proof:

To show that ${\hat{J}_{\mathrm{na}}(\lambda)}/{\hat{J}_{\mathrm{gd}}(\lambda)}$ is a decreasing function of $\lambda\in[m,L]$, we split this ratio into the sum of two homographic functions ${\hat{J}_{\mathrm{na}}(\lambda)}/{\hat{J}_{\mathrm{gd}}(\lambda)}=\sigma_{1}(\lambda)+\sigma_{2}(\lambda),$ where

[TABLE]

Now, if we substitute the parameters provided in Table II into (38), it follows that the signs of the derivatives $\mathrm{d}\sigma_{1}/\mathrm{d}\lambda$ and $\mathrm{d}\sigma_{2}/\mathrm{d}\lambda$ satisfy

[TABLE]

Furthermore, since the critical points of the functions $\sigma_{1}(\lambda)$ and $\sigma_{2}(\lambda)$ are outside the interval $[m,L]$,

[TABLE]

we conclude that both $\sigma_{1}$ and $\sigma_{2}$ are decreasing functions over the interval $[m,L]$. We next prove (13a) and (13b).

It is straightforward to verify that both $\hat{J}_{\mathrm{gd}}(\lambda)$ and $\hat{J}_{\mathrm{na}}(\lambda)$ are quasi-convex functions over the interval $[m,L]$ and that the respective minima are attained at the critical point $\lambda=1/\alpha$. Quasi-convexity also implies

[TABLE]

Now, letting $\alpha=2/(L+m)$ in the expression for $\hat{J}_{\mathrm{gd}}$ gives $\hat{J}_{\mathrm{gd}}(m)=\hat{J}_{\mathrm{gd}}(L)=(\kappa+1)^{2}/(4\kappa)$ which in conjunction with (39) complete the proof for (13a). Finally, since the ratio $\hat{J}_{\mathrm{na}}(\lambda)/\hat{J}_{\mathrm{gd}}(\lambda)$ is decreasing, we have ${\hat{J}_{\mathrm{na}}(L)}/{\hat{J}_{\mathrm{gd}}(L)}\leq{\hat{J}_{\mathrm{na}}(m)}/{\hat{J}_{\mathrm{gd}}(m)}.$ Combining this inequality with $\hat{J}_{\mathrm{gd}}(m)=\hat{J}_{\mathrm{gd}}(L)$ and (39) completes the proof of (13b). ∎

Proof:

From Proposition 1, it follows that

[TABLE]

For the upper bound, we have

[TABLE]

where the first inequality follows from (40a). The second inequality can be verified by multiplying both sides with the product of the denominators and using $\hat{J}_{\mathrm{gd}}(m)=\hat{J}_{\mathrm{gd}}(L)$, $\hat{J}_{\mathrm{na}}(m)\geq\hat{J}_{\mathrm{na}}(L)$, and (40b). Similarly, for the lower bound we can write

[TABLE]

Again, the first inequality follows from (40a) and the second inequality can be verified by multiplying both sides with the product of the denominators and using $\hat{J}_{\mathrm{gd}}(m)=\hat{J}_{\mathrm{gd}}(L)$, $\hat{J}_{\mathrm{na}}(m)\geq\hat{J}_{\mathrm{na}}(L)$, and (40b). ∎

Proof:

From Proposition 1, we have

[TABLE]

where $b\mathrel{\mathop{:}}=\sqrt{3\,\kappa+1}>2$. The upper and lower bounds on $\hat{J}_{\mathrm{na}}(m)$ are obtained as follows

[TABLE]

where the positive quantities $c_{1}(b)\mathrel{\mathop{:}}=b-2$ and $c_{2}(b)\mathrel{\mathop{:}}=b^{2}-2b$ are added to yield a simple upper bound. Similarly, for $\hat{J}_{\mathrm{na}}(L)$ we have

[TABLE]

where the positive quantities $c_{3}(b)\mathrel{\mathop{:}}=3b-3$, $c_{4}(b)\mathrel{\mathop{:}}=b^{2}-2b$, $c_{5}(b)\mathrel{\mathop{:}}=b-1$, and $c_{6}(b)\mathrel{\mathop{:}}=(1/2)b^{2}-b+1$ are introduced to obtain tractable bounds. ∎

-B General strongly convex problems

Proof:

Let us define the positive semidefinite function $V(\psi)\mathrel{\mathop{:}}=\psi^{T}X\psi$ and let $\eta\mathrel{\mathop{:}}=[\,\psi^{T}\;u^{T}\,]^{T}$. Using LMI (25) and (24), we can write

[TABLE]

Since $w^{t}$ is a zero-mean white input with identity covariance which is independent of $u^{t}$ and $x^{t}$, if we take the average of the above inequality over $t$ and expectation over different realizations of $w^{t}$, we obtain

[TABLE]

Therefore, letting $\bar{T}\rightarrow\infty$ and using $X\succeq 0$ lead to J\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}\mathrm{trace}\,(B_{w}^{T}\,X\,B_{w}), which completes the proof. ∎

In order to prove Lemma 2, we present a technical lemma which along the lines of results of [52] provides us with an upper bound on the difference between the objective value at two consecutive iterations.

Lemma 3

Let $f\in\mathcal{F}_{m}^{L}$ and $\kappa\mathrel{\mathop{:}}=L/m.$ Then, Nesterov’s accelerated method, with the notation introduced in Section IV, satisfies

[TABLE]

where $N_{1}$ and $N_{2}$ are defined in Lemma 2.

Proof:

For any $f\in\mathcal{F}_{m}^{L}$, the Lipschitz continuity of $\nabla f$ implies

[TABLE]

and the strong convexity of $f$ yields

[TABLE]

Moreover, the state and output equations in (5) lead to

[TABLE]

Summing up inequalities (41) and (42) and substituting for \left[\begin{array}[]{c}x^{t+2}\,-\,y^{t}\\ \nabla f(y^{t})\end{array}\right] and \left[\begin{array}[]{c}x^{t+2}\,-\,y^{t}\\ \nabla f(y^{t})\end{array}\right] from (53) completes the proof. ∎

Proof:

Let us define the positive semidefinite function $V(\psi)\mathrel{\mathop{:}}=\psi^{T}X\psi$ and let $\eta\mathrel{\mathop{:}}=[\,\psi^{T}\;\,u^{T}\,]^{T}$. Similar to the first part of the proof of Lemma 1, we can use LMI (26) and inequality (19) to write

[TABLE]

From Lemma 3, it follows that

[TABLE]

Now, combining inequalities (54) and (55) yields

[TABLE]

Since $w^{t}$ is a zero-mean white input with identity covariance which is independent of $u^{t}$ and $x^{t}$, taking the expectation of the last inequality yields

[TABLE]

and taking the average over the first $\bar{T}$ iterations results in

[TABLE]

Finally, using positive definiteness of the function $V$, strong convexity of the function $f$, and letting $\bar{T}\rightarrow\infty$, it follows that J\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sigma^{2}}(nL\lambda_{2}+\mathrm{trace}\,(B_{w}^{T}X\,B_{w})) as required. ∎

Proof:

Using Theorem (1), it is straightforward to show that for gradient descent and Nesterov’s method with the parameters provided in Table I, the function $f(x)\mathrel{\mathop{:}}=\tfrac{m}{2}\|x\|^{2}$ leads to the largest variance amplification $J$ among the quadratic objective functions within $\mathcal{F}_{m}^{L}$. This yields the lower bounds

[TABLE]

with $J_{\mathrm{gd}}$ and $J_{\mathrm{na}}$ corresponding to $f(x)=\tfrac{m}{2}\|x\|^{2}$. We next show that $J_{\mathrm{gd}}\leq q_{\mathrm{gd}}$.

To obtain the best upper bound on $J_{\mathrm{gd}}$ using Lemma 1, we minimize $\mathrm{trace}\,(B_{w}^{T}XB_{w})$ subject to LMI (25), $X\succeq 0$, and $\lambda\geq 0$. For gradient descent, if we use representation (23c), then the negative definiteness of the $(1,1)$-block of LMI (25) implies that

[TABLE]

It is straightforward to show that the pair

[TABLE]

is feasible as the LMI (25) becomes

[TABLE]

Thus, $X$ and $\lambda$ given by (62) provide a solution to LMI (25). Therefore, inequality (61) is tight and it provides the best achievable upper bound

[TABLE]

Finally, we show $J_{\mathrm{na}}\leq 4.08q_{\mathrm{na}}$ by finding a sub-optimal feasible point for (28). Let X\mathrel{\mathop{:}}=\left[\begin{array}[]{cc}x_{1}I&x_{0}I\\ x_{0}I&x_{2}I\end{array}\right] with

[TABLE]

and let $\lambda_{1}\mathrel{\mathop{:}}=(\kappa/L)^{2}/(2\kappa-1)$ and $\lambda_{2}\mathrel{\mathop{:}}={-x_{0}}/(Ls(\kappa)).$ We first show that $(\lambda_{1},\lambda_{2},X)$ is feasible for problem (28). It is straightforward to verify that $s(\kappa)$, $x_{1}s(\kappa)$, $x_{2}s(\kappa)$, and $-x_{0}s(\kappa)$ (which are polynomials of degree less than $7$ in $\sqrt{\kappa}$) are all positive for any $\kappa\geq 1$. Hence, $x_{1}>0$, $x_{2}>0$ and $\lambda_{2}>0$. It is also easy to see that $\lambda_{1}>0$ and that the determinant of $X$ satisfies

[TABLE]

which yields $X\succeq 0$. Moreover, it can be shown that the left-hand-side of LMI (26) becomes

[TABLE]

Therefore, the point $(\lambda_{1},\lambda_{2},X)$ is feasible to problem (28) and

[TABLE]

Comparing $p$ with $q_{\mathrm{na}}$, it can be verified that, for all $\kappa\geq 1$, $4.08q_{\mathrm{na}}(\kappa)\geq p(\kappa),$ which completes the proof. ∎

-C Proof of Theorem 6

Without loss of generality, let $\sigma=1$ and

[TABLE]

where $\lambda_{i}$ are the eigenvalues of the Hessian of the objective function $f$ and $\lambda_{i}^{\prime}=m+L-\lambda_{i}$ is the mirror image of $\lambda_{i}$ with respect to $(m+L)/2$. Since $J=\sum_{i}\hat{J}(\lambda_{i})$, if $\lambda_{i}$ are symmetrically distributed over the interval $[m,L]$ i.e., $(\lambda_{1},\cdots,\lambda_{n})=(\lambda_{n}^{\prime},\cdots,\lambda_{1}^{\prime})$, then for any parameters $\alpha$ and $\beta$ we have

[TABLE]

Equation (64) implies that any bound on $G$ simply carries over to $J$ within an accuracy of constant factors. Thus, we focus on $G$ and establish one of its useful properties in the next lemma that allows us to prove Theorem 6.

Lemma 4

The heavy-ball method with any stabilizing parameter $\beta$ satisfies

[TABLE]

where $\rho$ is the rate of linear convergence. Furthermore, if the Hessian of the quadratic objective function $f$ has a symmetric spectrum over the interval $[\lambda_{1},\lambda_{n}]=[m,L]$, then

[TABLE]

Proof:

The linear convergence rate $\rho$ is given by $\rho=\max_{1\,\leq\,i\,\leq\,n}\hat{\rho}(\lambda_{i})$, where $\hat{\rho}(\lambda)$ is the largest absolute value of the roots of the characteristic polynomial

[TABLE]

associated with the heavy-ball method and the eigenvalue $\lambda$ of the Hessian of the objective function $f$; See (8) for the form of $\hat{A}$. Thus, we have

[TABLE]

where $\Delta\mathrel{\mathop{:}}=(1+\beta-\alpha\lambda)^{2}-4\beta.$ This can be simplified to

[TABLE]

It is straightforward to show that $\hat{\rho}$ and $\hat{J}$ with $\sigma=1$ are explicit quasi-convex functions of $\mu\mathrel{\mathop{:}}=\alpha\lambda$ which are symmetric with respect to $\mu=1+\beta$. Quasi-convexity of $\hat{\rho}$ yields

[TABLE]

Let $\alpha^{\sharp}(\beta)=2(1+\beta)/(L+m)$. For any eigenvalue $\lambda_{i}$, from the symmetry of the spectrum, we have

[TABLE]

meaning that $\alpha^{\sharp}(\beta)\lambda_{i}$ and $\alpha^{\sharp}(\beta)\lambda^{\prime}_{i}$ are the mirror images with respect to the middle point $1+\beta$. Thus, from the quasi-convexity and symmetry of the functions $\hat{\rho}$ and $\hat{J}$, it follows that $\alpha^{\sharp}(\beta)$ minimizes $\rho$ as well as $\max\,\{\hat{J}(\lambda_{i}),\hat{J}(\lambda_{i}^{\prime})\}$ for all $i$, which completes the proof. ∎

Since gradient descent is obtained from the heavy-ball method by letting $\beta=0$, from Lemma 4 it immediately follows that $\alpha_{\mathrm{gd}}=2/(L+m)$ given in Table II optimizes both $G_{\mathrm{gd}}$ and the convergence rate $\rho_{\mathrm{gd}}$. This fact combined with (64) yields

[TABLE]

where $\alpha^{\star}_{\mathrm{gd}}(c)$ is given by (30b). This completes the proof for gradient descent.

We next use Lemma 4 to establish a bound on the parameter $\beta^{\star}_{\mathrm{hb}}(c)$ that allows us to prove the result for the heavy-ball method as well.

Lemma 5

There exists a positive constant $a$ such that

[TABLE]

where $\beta^{\star}_{\mathrm{hb}}(c)$ is given by (30a).

Proof:

We first show that for any parameters $\alpha$ and $\beta$, the convergence rate $\rho$ of the heavy-ball method given by (29) is lower bounded by

[TABLE]

The convergence rate satisfies

[TABLE]

where the function $\hat{\rho}(\lambda)$ is given by (see proof of Lemma 4 for the proof of this statement)

[TABLE]

and $\Delta\mathrel{\mathop{:}}=(1+\beta-\alpha\lambda)^{2}-4\beta.$ According to Lemma 4, $\alpha=2(1+\beta)/(L+m)$ optimizes the rate $\rho$. This value of $\alpha$ yields

[TABLE]

or equivalently

[TABLE]

which completes the proof of inequality (70). Now, if $\beta\geq(\sqrt{\kappa}-1)^{2}/(\sqrt{\kappa}+1)^{2}$, then (67) with $a=2$ follows immediately. Otherwise, from (70) we obtain

[TABLE]

which yields

[TABLE]

The convergence rate $\rho$ satisfies $(\sqrt{\kappa}-1)^{2}/(\sqrt{\kappa}+1)^{2}\leq\rho\leq 1-c/\sqrt{\kappa}$, where the lower bound follows from the optimal rate provided in Table II and the upper bound follows from the definition in (30a). Moreover, the derivative $\frac{\mathrm{d}v}{\mathrm{d}\rho}=0$ vanishes only at $\rho=(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1)$. Thus, we obtain a lower bound on $\beta$ as

[TABLE]

A simple manipulation of (75) allows us to find a constant $a$ that satisfies (67), which completes the proof. ∎

Let $(\hat{\alpha},\hat{\beta})$ be the optimal solution of the optimization problem

[TABLE]

where $G$ is defined in (63). We next show that there exists a scalar $c^{\prime}>0$ such that

[TABLE]

where $\alpha_{\mathrm{hb}}$ and $\beta_{\mathrm{hb}}$ are provided in Table II. Let $\hat{\alpha}(\beta)\mathrel{\mathop{:}}=2(1+\beta)/(L+m)$. It is straightforward to verify that

[TABLE]

which allows us to write

[TABLE]

Here, (i) determines partial minimization with respect to $\alpha$ which follows from Lemma 4; (ii) follows from (64); (iii) follows from (77), and (iv) follows from Lemma 5. Furthermore, it is easy to show the existence of a constant scalar $c^{\prime}$ such that

[TABLE]

Inequality (76) follows from combining (79) and (78). Finally, we obtain that

[TABLE]

where the first inequality follows from (64), the second inequality follows from the definition of $(\hat{\alpha},\hat{\beta})$, and the last inequality is given by (76). This completes the proof of Theorem 6 for the heavy-ball method.

-D Fundamental lower bounds

Proof:

We first prove (31a). Without loss of generality, let the noise magnitude $\sigma=1$. We define the trivial lower bound

[TABLE]

and show that $\dfrac{\hat{J}^{\star}}{1-\rho}\;\geq\;(\dfrac{\kappa+1}{8})^{2}.$ Let $\tilde{f}(x_{1},x_{2})\mathrel{\mathop{:}}=\tfrac{1}{2}\,(m\,x_{1}^{2}+L\,x_{2}^{2})$. The eigenvalues of the Hessian matrix $\nabla^{2}\tilde{f}$ are given by $m$ and $L$ which are clearly symmetric over the interval $[m,L]$. Thus, for any given value of $\beta$, $m$, and $L$, we can use Lemma 4 with the objective function $\tilde{f}$ to obtain

[TABLE]

For the stepsize $\hat{\alpha}(\beta)$, the rate of convergence $\rho$ is given by (73), i.e.,

[TABLE]

and the lower bound $\hat{J}^{\star}$ is given by

[TABLE]

Therefore, we obtain a lower bound on $\hat{J}^{\star}/(1-\rho)$ as

[TABLE]

where the last equality follows from (83) and (84). It can be shown that $v(\beta)$ attains its minimum at $\beta=(\sqrt{\kappa}-1)^{2}/(\sqrt{\kappa}+1)^{2}$; see Figure 5 for an illustration.

Therefore,

[TABLE]

which completes the proof of (31a). We next prove (31b) for $\sigma=\alpha$.

We analyze the two cases $\alpha>1/L$ and $\alpha\leq 1/L$ separately. If $\alpha>1/L$, inequality (31b) directly follows from inequality (31a)

[TABLE]

Here, the first inequality is given by (31a) and the second inequality holds since $\alpha>1/L$.

Now suppose $\alpha\leq 1/L$. The convergence rate of Polyak’s method is given by $\max_{i}\hat{\rho}(\lambda_{i})$, where

[TABLE]

and $\Delta\mathrel{\mathop{:}}=(1+\beta-\alpha\lambda)^{2}-4\beta$ (see the proof of Lemma 4). Thus, for $\sigma=\alpha$, we have the trivial lower bound

[TABLE]

Here, the first inequality follows from combining $J=\sum_{i}\hat{J}(\lambda_{i})$ and $\max_{i}\hat{\rho}(\lambda_{i})$, and the second inequality follows from $\alpha m\leq\alpha L\leq 1$. We next show that for any fixed $\alpha$, the function $p(\alpha,\cdot)$ attains its minimum at $\beta=(1-\sqrt{\alpha m})^{2}$. Before we do so, note that this fact allows us to use partial minimization with respect to $\beta$ and obtain

[TABLE]

which completes the proof of (31b).

For any fixed $\alpha$, it is straightforward to verify that $p(\alpha,\beta)$ is increasing with respect to $\beta$ over $[(1-\sqrt{\alpha m})^{2},\;1)$. Thus, it suffices to show that $p(\alpha,\beta)$ is decreasing with respect to $\beta$ over $[0,(1-\sqrt{\alpha m})^{2})$. To simplify the presentation, let us define the new set of parameters

[TABLE]

It is now straightforward to verify that $p(\alpha,\beta)=\alpha/(mq)$ for $\beta\in[(1-\sqrt{\alpha m})^{2},\;1)$. It thus follows that $p(\alpha,\beta)$ is decreasing with respect to $\beta$ over $[0,(1-\sqrt{\alpha m})^{2})$ if and only if $q^{\prime}=\mathrm{d}q/\mathrm{d}s\leq 0$ for $s\in(\sqrt{x}(2-\sqrt{x}),1]$. The derivative is given by

[TABLE]

Thus, we have

[TABLE]

It is easy to verify that both sides of the inequality in (88), namely, $(2s+x)\delta$ and $2s^{2}+3sx+x^{2}-4x$ are positive for the specified range of $s\in(\sqrt{x}(2-\sqrt{x}),1]$. Thus, we can square both sides and obtain that

[TABLE]

where $(i)$ follows from the definition of $\delta$ and $(ii)$ is obtained by expanding both sides and rearranging the terms. Finally, the inequality $8s+4x\;\leq\;16$ clearly holds since $s\leq 1$ and $x\leq 1$. This proves that $p(\alpha,\cdot)$ attains its minimum at $\beta=(1-\sqrt{\alpha m})^{2}$.

∎

Proof:

For the heavy-ball method, the result follows from combining Theorem 7 and the inequality $1-\rho>{c}/{\sqrt{\kappa}}$. Next, we present three additional lemmas that allow us to prove the result for Nesterov’s method.

The following lemma provides a lower bound on the function $\hat{J}(m)$ associated with Nesterov’s method which depends on $\kappa$ and $\beta$.

Lemma 6

For any strongly convex quadratic problem with condition number $\kappa>2$ and the smallest eigenvalue of the Hessian $m$, the function $\hat{J}$ associated with Nesterov’s accelerated method with any stabilizing pair of parameters $0<\alpha$, $0<\beta<1$, and $\sigma=1$ satisfies

[TABLE]

Proof:

We first show that Nesterov’s method with $0<\alpha$ and $0<\beta<1$ is stable if and only if

[TABLE]

The rate of linear convergence is given by $\rho=\max_{1\leq i\leq n}\;\hat{\rho}(\lambda_{i})$, where $\hat{\rho}(\lambda)$ is the largest absolute value of the roots of the characteristic polynomial

[TABLE]

associated with Nesterov’s method and the eigenvalue $\lambda$ of the Hessian of the objective function $f$; See (8) for the form of $\hat{A}$. For $\alpha>0$ and $0<\beta<1$, it can be shown that

[TABLE]

The stability of the algorithm is equivalent to $\hat{\rho}(\lambda_{i})<1$ for all eigenvalues $\lambda_{i}$. For any positive stepsize $\alpha$ and parameter $\beta\in(0,1)$, it can be shown that the function $\hat{\rho}(\lambda)$ is quasi-convex and $\hat{\rho}(\lambda)=1$ if and only if $\lambda\in\{0,\frac{2\beta+2}{\alpha(2\beta+1)}\}$. This fact along with $0<m\leq\lambda_{i}\leq L=\kappa\,m$ imply that $\hat{\rho}(\lambda_{i})<1$ for all $\lambda_{i}\in[m,L]$ if and only if $\kappa\,m\leq\frac{2\beta+2}{\alpha(2\beta+1)}$ which completes the proof of (90).

For Nesterov’s method, it is straightforward to show that the function $\hat{J}(\lambda)$ is quasi-convex over the interval $[0,\frac{2\beta+2}{\alpha(2\beta+1)}]$ and that it attains its minimum at $\lambda=1/\alpha$. Also, from (90), for $\kappa>2$ we obtain

[TABLE]

and thus,

[TABLE]

where the last inequality follows from the fact that $\beta\in(0,1)$. ∎

The following lemma presents a lower bound on any accelerating parameter $\beta$ for Nesterov’s method.

Lemma 7

For Nesterov’s method, under the conditions of Theorem 8, there exist positive constants $c_{3}$ and $c_{4}$ such that for any $\kappa>c_{3}$,

[TABLE]

Proof:

For any $\alpha>0$ and $\beta\in(0,1)$, Nesterov’s method converges with the rate $\rho=\max_{1\,\leq\,i\,\leq\,n}\hat{\rho}(\lambda_{i})$, where $\hat{\rho}(\lambda)$ is given by (93). We treat the two cases $(1-\beta)/(1+\beta)^{2}<\alpha m$ and $(1-\beta)/(1+\beta)^{2}\geq\alpha m$ separately. For $(1-\beta)/(1+\beta)^{2}<\alpha m$, we have

[TABLE]

where the last inequality follows from (90). Therefore, we obtain $\beta\geq 1-\sqrt{8}/\sqrt{\kappa}$ as required. Now, suppose $(1-\beta)/(1+\beta)^{2}\geq\alpha m$. The convergence rate $\rho$ satisfies

[TABLE]

Thus,

[TABLE]

which yields a lower bound on $\beta$,

[TABLE]

In what follows, we establish a lower bound for $\nu$. For a fixed $\alpha m$, the critical point of $\nu(\rho)$ is given by $\rho_{1}\mathrel{\mathop{:}}=1-\sqrt{\alpha m}$, i.e., ${\partial\nu}/{\partial\rho}=0$ for $\rho=\rho_{1}$. Furthermore, the optimal rate from Table II and the condition on convergence rate in Theorem 8 for any $\kappa>c_{1}$ yield upper and lower bounds $\rho_{3}<\rho<\rho_{2}$, where $\rho_{2}\mathrel{\mathop{:}}=1-{c_{2}}/{\sqrt{\kappa}}$ and $\rho_{3}\mathrel{\mathop{:}}=1-{2}/\sqrt{3\kappa+1}$. Thus, the lower bound on $\nu$ is given by

[TABLE]

From the stability condition (90), we have

[TABLE]

Furthermore, it can be shown that for any given $\rho\in(0,1)$ the function $\nu(\rho,\alpha m)$ is decreasing with respect to $\alpha m$. This fact combined with (97) and (98) yield

[TABLE]

If we substitute for $\rho_{1}$. $\rho_{2}$, and $\rho_{3}$ their values as functions of $\kappa$ and use $\alpha m<{2}/{\kappa}$, then the result follows immediately. In particular,

[TABLE]

which completes the proof. ∎

The next lemma provides a lower bound on $J_{\mathrm{na}}/(1-\rho)$ for Nesterov’s method with $\sigma=\alpha\leq 1/L$.

Lemma 8

Nesterov’s accelerated method with any stabilizing pair of parameters $0<\alpha\leq 1/L$ and $0<\beta<1$, and $\sigma=\alpha$ satisfies

[TABLE]

Proof:

The convergence rate of Nesterov’s method is given by $\max_{i}\hat{\rho}(\lambda_{i})$, where

[TABLE]

and $\Delta\mathrel{\mathop{:}}=(1+\beta)^{2}(1-\alpha\lambda)^{2}-4\beta(1-\alpha\lambda)$; see equation (93). Thus, we have the trivial lower bound

[TABLE]

where $\gamma\mathrel{\mathop{:}}=\dfrac{1-\sqrt{\alpha m}}{1+\sqrt{\alpha m}}$. Here, the first inequality can be obtained by combining $J=\sum_{i}\hat{J}(\lambda_{i})$ and $\max_{i}\hat{\rho}(\lambda_{i})$, and the second inequality follows from the fact that $0<\alpha m\leq 1$ and $0\leq\beta<1$. We next show that for any fixed $\alpha$, the function $p(\alpha,\cdot)$ attains its minimum at $\beta=\gamma$. Before we do so, note that this fact allows us to do partial minimization with respect to $\beta$ and obtain

[TABLE]

For any fixed $\alpha$, it is straightforward to verify that $p(\alpha,\beta)$ is increasing with respect to $\beta$ over $[\gamma,\;1)$. Thus, it suffices to show that $p(\alpha,\beta)$ is decreasing with respect to $\beta$ over $[0,\gamma)$. To simplify the presentation, let us define

[TABLE]

It is now straightforward to verify that $p(\alpha,\beta)=\alpha/(2mq)$ for $\beta\in[0,\gamma)$. It thus follows that $p(\alpha,\beta)$ is decreasing with respect to $\beta$ over $[0,\gamma)$ if and only if $q^{\prime}=\mathrm{d}q/\mathrm{d}s\geq 0$ for $s\in[0,(1-\sqrt{1-x})^{2})$. The derivative is given by

[TABLE]

Thus, we have

[TABLE]

It is easy to verify that both sides of the inequality in (103), namely, $(1-s)(2-x-s)+\delta^{2}$ and $(3-x-2s)\delta$ are positive for the specified range of $s\in[0,(1-\sqrt{1-x})^{2})$. Thus, we can square both sides and obtain that

[TABLE]

where $(i)$ follows from the definition of $\delta$ and $(ii)$ is obtained by expanding both sides and rearranging the terms. Finally, the inequality $4(x-1)^{2}(2s+x+1)\;\geq\;0$ trivially holds which completes the proof. ∎

We are now ready to prove Theorem 8 for Nesterov’s method. Inequality (32a) directly follows from combining (89) in Lemma 6 and (94) in Lemma 7. To show inequality (32b), we treat the two cases $\alpha>1/L$ and $\alpha\leq 1/L$ separately. If $\alpha>1/L$, then (32b) directly follows from (32a)

[TABLE]

Now suppose $\alpha\leq 1/L$. We can use Lemma 8 to obtain

[TABLE]

Here, the first inequality follows from Lemma 8 and the second inequality follows from the acceleration assumption $\rho\leq 1-c/\sqrt{\kappa}$. This completes the proof. ∎

-E Consensus over $d$-dimensional torus networks

The proof of Theorem 9 uses the explicit expression for the eigenvalues of torus in (35) to compute the variance amplification $\bar{J}=\sum_{i\neq 0}\hat{J}(\lambda_{i})$ for all three algorithms. Several technical results that we use in the proof are presented next.

We borrow the following lemma, which provides tight bounds on the sum of reciprocals of the eigenvalues of a $d$-dimensional torus network, from [38, Appendix B].

Lemma 9

The eigenvalues $\lambda_{i}$ of the graph Laplacian of the $d$-dimensional torus $\mathbb{T}_{n_{0}}^{d}$ with ${n_{0}}\gg 1$ satisfy

[TABLE]

where the function $B$ is given by

[TABLE]

We next use Lemma 9 to establish an asymptotic expression for the variance amplification of the gradient descent algorithm for a $d$-dimensional torus.

Lemma 10

For the consensus problem over a $d$-dimensional torus $\mathbb{T}_{n_{0}}^{d}$ with ${n_{0}}\gg 1$, the performance metric $\bar{J}_{\mathrm{gd}}$ corresponding to gradient decent with the stepsize $\alpha=2/(L+m)$ satisfies

[TABLE]

where the function $B$ is given in Lemma 9.

Proof:

Using the expression for the noise amplification of gradient descent from Theorem 1, we have

[TABLE]

The first approximation follows from the facts that the eigenvalues satisfy

[TABLE]

and that their distribution is asymptotically symmetric with respect to $\lambda=2d$. The second approximation follows from

[TABLE]

The bounds for the sum of reciprocals of $\lambda_{i}$ provided in Lemma 9 can now be used to complete the proof. ∎

The following lemma establishes a relationship between the variance amplifications of Nesterov’s method and gradient descent. This relationship allows us to compute tight bounds on $J_{\mathrm{na}}$ by splitting it into the sum of two terms. The first term depends linearly on $J_{\mathrm{gd}}$ which is already computed in Lemma 10 and the second term can be evaluated separately using integral approximations for consensus problem on torus networks. This result holds in general for the scenarios in which the largest eigenvalue $L=\Theta(1)$ is bounded and the smallest eigenvalue $m$ goes to zero causing the condition number $\kappa$ to go to infinity.

Lemma 11

For a strongly convex quadratic problem with $mI\preceq Q\preceq LI$ and condition number $\kappa\mathrel{\mathop{:}}=L/m\geq\kappa_{0}$, the ratio between variance amplifications of Nesterov’s algorithm and gradient descent with the parameters given in Table II satisfies the asymptotic bounds

[TABLE]

where $\kappa_{0}$, $c_{1}$, and $c_{2}$ are positive constants. Furthermore, depending on the distribution of the eigenvalues of the Laplacian matrix, $D$ can take values between

[TABLE]

where $c_{3}$ and $c_{4}$ are positive constants.

Proof:

We can split $\hat{J}_{\mathrm{na}}(\lambda)/\hat{J}_{\mathrm{gd}}(\lambda)$ into the sum of two decreasing homographic functions $\sigma_{1}(\lambda)\,+\,\sigma_{2}(\lambda),$ where $\sigma_{1}$ and $\sigma_{2}$ are defined in (38); see the proof of Proposition 1. Furthermore, for $\kappa\gg 1$, these functions attain their extrema over the interval $[m,L]$ at

[TABLE]

where we have kept the leading terms. It is straightforward to verify that

[TABLE]

This equation in conjunction with (105), yield inequalities in (104). Moreover, we obtain that

[TABLE]

This also implies that, asymptotically,

[TABLE]

which completes the proof. ∎

The next two lemmas provide us with asymptotic bounds on summations of the form $\sum_{i}{1}/{(\lambda_{i}^{2}+\mu\lambda_{i})}$, where $\lambda_{i}$ are the eigenvalues of the graph Laplacian matrix of a torus network. These bounds allow us to combine Lemma 10 and Lemma 11 to evaluate the variance amplification of Nesterov’s accelerated algorithm.

Lemma 12

For an integer $q\gg 1$ and any positive $a=O(q^{3})$, we have

[TABLE]

where $\omega={a}/{q^{2}}$.

Proof:

The function $h(x)\mathrel{\mathop{:}}=\|x\|^{4}+\omega\|x\|^{2}$ is strictly increasing over the positive orthant ($x\succ 0$) and $h((1/q){\mathds{1}})$ goes to 0 as $q$ goes to infinity where $\mathds{1}\in\mathbb{R}^{d}$ is the vector of all ones. Therefore, using the lower and upper Riemann sum approximations, it is straightforward to show that

[TABLE]

where $\Delta={1}/{q}$ is the incremental step in the Riemann approximation. Therefore, since $\omega=a\Delta^{2}$, we can write

[TABLE]

Finally, we obtain the result by transforming the integral into a $d$-dimensional polar coordinate system, i.e.,

[TABLE]

∎

Lemma 13

Let $\lambda_{i}$ be the eigenvalues of the Laplacian matrix for the $d$-dimensional torus $\mathbb{T}_{{n_{0}}}^{d}$. In the limit of large ${n_{0}}$, for any $\mu=O({n_{0}})$, we have

[TABLE]

where $\omega=\Theta(\mu)$.

Proof:

Let $\zeta\,\mathrel{\mathop{:}}=\,\sum_{0\neq i\in\mathbb{Z}_{n_{0}}^{d}}\dfrac{1}{\lambda_{i}^{2}+\mu\lambda_{i}}$, where $\lambda_{i}=2\sum_{l=1}^{d}\left(1-\cos(i_{l}\dfrac{2\pi}{{n_{0}}})\right)$ are the eigenvalues of the graph Laplacian matrix. Since $1-\cos(\cdot-\pi)$ is an even function, for large ${n_{0}}$,

[TABLE]

where $q=\lfloor{{n_{0}}}/{2}\rfloor$. It is well-known that the function $1-\cos(x)$ can be bounded by quadratic functions as ${x^{2}}/{\pi^{2}}\leq 1-\cos(x)\leq x^{2}$ for any $x\in[-\pi,\pi]$. Now, since for any $i\in\mathbb{Z}_{q}^{d}$, $i_{l}\tfrac{2\pi}{{n_{0}}}\in[0,\pi]$ for all $l$, we can use these quadratic bounds to obtain

[TABLE]

where $c$ is a bounded constant. Finally, equation (106) follows from Lemma 12 where we let $a=c\mu n_{0}^{2}$ and $q\approx{{n_{0}}}/{2}$. ∎

The following proposition characterizes the network-size-normalized asymptotic variance amplification of noisy consensus algorithms for $d$-dimensional torus networks. This result is used to prove Theorem 9.

Proposition 2

Let $\mathbf{L}\in\mathbb{R}^{n\times n}$ be the graph Laplacian of the $d$-dimensional undirected torus $\mathbb{T}_{n_{0}}^{d}$ with $n=n_{0}^{d}\gg 1$ nodes. For convex quadratic optimization problem (33), the network-size-normalized asymptotic variance amplification $\bar{J}/n$ of the first-order algorithms on the subspace $\mathds{1}^{\perp}$ is determined by

[TABLE]

Proof:

We prove the result for the three algorithms separately.

For gradient descent, the result follows from dividing the asymptotic bounds established in Lemma 10 with the total number of nodes $n=n_{0}^{d}$. 2. 2.

For Nesterov’s algorithm, we use the relation established in Lemma 11 to write

[TABLE]

where $c=2/\left((3\beta+1)\alpha_{\mathrm{na}}^{2}\right)\approx 9d^{2}/2$ and $\mu=(1-\beta)/(\alpha_{\mathrm{na}}\beta)=\Theta({1}/{\sqrt{\kappa}})=\Theta(n_{0}^{-1})$; see equation (36). We can use Lemma 13 to compute the second term

[TABLE]

where $\omega=\Theta(\mu)=\Theta(n_{0}^{-1})$. Evaluating the above integral for different values of $d\in\mathbb{N}$ and letting $\omega=\Theta(n_{0}^{-1})$, it is straightforward to show that

[TABLE]

Finally, the result follows from the asymptotic values for $\bar{J}_{\mathrm{gd}}/n$ (shown in Part 1) and substituting for the second term on the left-hand-side of equation (108) from the above asymptotic values and using $n=n_{0}^{d}$. We note that we used the following integrals to evaluate $\bar{J}_{\mathrm{na}}$,

[TABLE] 3. 3.

The result for the heavy-ball method directly follows from the first part of the proof, the relationship between variance amplifications of gradient descent and the heavy-ball method in Theorem 2, and equation (36).

∎

We now use Proposition 2 to proof Theorem 9 as follows.

Proof of Theorem 9:

As stated in (36), the condition number satisfies $\kappa=\Theta(n^{2/d})$ and the result follows from combining this asymptotic relation with those provided in Proposition 2.

Computational experiments

To complement our asymptotic theoretical results, we compute the performance measure $\bar{J}$ in (34) for the consensus problem over $d$-dimensional torus $\mathbb{T}_{n_{0}}^{d}$ with $n=n_{0}^{d}$ nodes for different values of ${n_{0}}$ and $d$. We use expression (35) for the eigenvalues of the graph Laplacian $\mathbf{L}$ to evalute the formulae provided in Theorem 1 for each algorithm. Figure 6 illustrates network-size normalized variance amplification $\bar{J}/n$ vs. condition number $\kappa$ and verifies the asymptotic relations provided in Theorem 9. It is noteworthy that, even though our analysis is asymptotic in the condition number (i.e., it assumes that $\kappa\gg 1$), our computational experiments exhibit similar scaling trends for small values of $\kappa$ as well.