Robust Accelerated Gradient Methods for Smooth Strongly Convex Functions

Introduction

For many large-scale convex optimization and machine learning problems, first-order methods have been the leading computational approach for computing low-to-medium accuracy solutions because of their cheap iterations and mild dependence on the problem dimension and data size. The typical analysis of first-order methods assumes the availability of exact gradient information and provides statements on the rate of convergence to the optimal solution as the main performance criterion. However, in many applications, the gradient contains deterministic or stochastic errors either because the gradient is computed by inexactly solving an auxiliary problem , or the method itself involves errors with respect to the full gradient as in standard incremental gradient, stochastic gradient, and stochastic approximation methods . When there are persistent errors in gradients, the iterates do not converge and could oscillate in a neighborhood of the optimal solution or may even diverge . This makes robustness of the algorithms to gradient errors (in terms of solution accuracy) another important performance objective . In particular, even though accelerated gradient method proposed by Nesterov converges faster than gradient descent (GD) in the absence of noise for convex problems , it was shown that they are less robust to errors, i.e., accelerated methods require higher precision gradient information than GD to achieve the same solution accuracy .

In this paper, we study the trade-offs between convergence rate and robustness to gradient errors in designing a first-order algorithm. We focus on GD and Nesterov’s accelerated gradient (AG) method for minimizing strongly convex smooth functions when the gradient has stochastic errors and investigate how the parameters of each algorithm should be set to achieve a particular trade-off between these two performance objectives. To study this question systematically, we employ tools from control theory whereby we represent each of the algorithms as a dynamical system. This approach has attracted recent attention and has already led to a number of insights for the design and analysis of optimization algorithms . The novelty of our work is to use this approach to provide explicit characterizations of robustness, which can then be placed in a computationally tractable optimization problem for selecting the algorithm parameters to systematically achieve a desired trade-off.

We first focus on problems with a strongly convex quadratic objective function. For this case, the rate of convergence of any of the two algorithms we study is given by the spectral radius of the “state-transition" matrix in the dynamical system representation. To characterize robustness, we consider the asymptotic expected suboptimality for the centered iterate sequence (output vector of the dynamical system) per unit noise which is a measure of the asymptotic accuracy of the iterates. For the quadratic case we show that this limit exists and can be characterized using the $H_{2}$ norm of a transformed linear dynamical system. The $H_{2}$ norm is a fundamental measure for quantifying robustness of a linear system to noise and admits various definitions and characterizations. We focus on a particular representation of the $H_{2}$ norm that requires the solution of a discrete Lyapunov equation. This representation leads to explicit expressions for robustness of GD and AG.

Using this result, we study the rate and robustness trade-off of the GD method for minimizing quadratic strongly convex functions. The spectral radius of the state-transition matrix corresponding to GD dynamics, hence, the rate of convergence for GD, can be expressed in terms of the smallest and largest eigenvalues of the positive definite matrix $Q$ defining the Hessian of the strongly convex quadratic objective. We show that our robustness measure admits a tractable characterization for GD in terms of the spectrum of $Q$. We also show a fundamental lower bound on the robustness level of an algorithm for any achievable convergence rate.

We next consider the AG method defined by two parameters: stepsize $\alpha$ and momentum parameter $\beta$. Our first step is to characterize the stability region of the method, i.e., the set of nonnegative $(\alpha,\beta)$ for which the spectral radius of the state-transition matrix is less than or equal to one. Similar to GD, we then provide an explicit characterization of the $H_{2}$ norm of the dynamical system representation of AG. We use these explicit expressions for both GD and AG within an optimization problem for selecting the parameters to minimize the robustness measure subject to a given upper bound on the convergence rate. Our results show that AG with properly selected parameters is superior to GD in the sense that AG can achieve the same rate with GD while being more robust to noise; similarly, AG can be tuned to be faster than GD while achieving the same robustness level. This behavior contrasts with the comparison of GD and AG in the deterministic gradient error setting in , which shows GD performance degrades gracefully while AG may accumulate error. These results show the random and deterministic noise settings have different behavior.

In addition to the above cited papers, Devolder’s Ph.D. thesis is closely related to our paper. Chapters 4 and 6 of this thesis, considered smooth weakly convex functions under a deterministic oracle model whereas Chapter 7 focused on a stochastic oracle model; these general oracles can model inexactness in the gradients as well as function evaluations. In the deterministic oracle case, Devolder shows that primal gradient method (PGM) and the dual gradient method (DGM) on smooth weakly convex objectives exhibit slow convergence with a rate $\mathcal{O}(1/k)$ but without accumulation of errors (the total effect of errors after $k$ iterations is equal to the individual error $\delta$ of each first-order information); whereas accelerated gradient methods converge faster with rate $\mathcal{O}(\frac{1}{k^{2}})$ but suffers from accumulation of errors at a linear rate $\mathcal{O}(k\delta)$. Based on these observations, Devolder et al. design a novel family of first-order methods called intermediate gradient methods (IGM) for solving smooth weakly convex problems; these methods have an intermediate speed and intermediate sensitivity to gradient errors, i.e., faster than classical gradient methods and more robust to noise than the accelerated gradient methods. In the stochastic oracle case, Devolder developed a class of accelerated gradient methods for weakly convex functions with decaying stepsize rules and showed that the expected suboptimality admits the convergence rate $\mathcal{O}\left(\frac{LR^{2}}{k^{2}}+\frac{\sigma R}{\sqrt{k}}\right)$ as opposed to the $\mathcal{O}\left(\frac{LR^{2}}{k}+\frac{\sigma R}{\sqrt{k}}\right)$ rate of PGM and DGM, where $R$ is the distance of the initial point to the optimal solution, $L$ is the Lipschitz constant for the gradient of the objective $f(x)$ and $\sigma$ is the level of the stochastic noise [13, Ch. 7]. In his thesis, Devolder studied also smooth and strongly convex objectives under the same deterministic oracle model, showing that both PGM and DGM converge with a rate that is proportional to $\exp(-k\frac{\mu}{L})$ without accumulation of errors where $\mu$ is the strong convexity constant, whereas accelerated gradient converges faster proportional to $\exp(-k\sqrt{\frac{\mu}{L}})$ while the error accumulation behaves like $\sqrt{\frac{L}{\mu}}\delta$ up to a constant [13, Chapter 5]. On the other hand, the smooth and strong convex objectives subject to stochastic errors was left as future work [13, Ch. 8.1.1]; and this is the setting considered in our paper where we focus on stochastic additive gradient errors for strongly convex objectives, which arises in a number of problems in machine learning and large-scale optimization . In this setting, Ghadimi and Lan propose an accelerated method called AC-SA for solving strongly convex composite optimization problems obtaining an optimal rate matching the lower complexity bounds for stochastic optimization. Flammarion and Bach considered accelerated versions of gradient descent for quadratic optimization that attain the optimal rates for both the bias and variance terms, respectively, in the performance bounds. Michalowsky and Ebenbauer posed the design of deterministic gradient algorithms as a state feedback problem and used robust control theory and linear matrix inequalities to study them. Mohammadi et al. examined the sensitivity of accelerated algorithms to stochastic noise for strongly convex quadratic functions in terms of the steady-state variance of the optimization variable. Finally, Dvurechensky et al. consider composite convex optimization problems with inexact first-order oracles having both deterministic and stochastic errors; indeed, their inexact oracle is an extension of the one adopted in to include stochastic errors. For this setting, Dvurechensky et al. propose a stochastic version of the intermediate gradient method in and analyze the convergence rate in terms of expected suboptimality and error accumulation due to inexact oracle; the proposed algorithm in has complexity bounds matching the optimal lower complexity bounds for composite convex problems with stochastic inexact oracle as in . Finally, Hu et al. analyze the stochastic gradient method under deterministic noise and study the effect of the stepsize on the convergence rate and the asymptotic neighborhood of convergence. These papers focus on convergence rate of the algorithms, whereas our goal is to define robustness and design algorithms to successfully trade-off different objectives. Furthermore, we make some connections between the robustness of a first-order method and its behavior when perturbed from the optimal solution and show that AG is more resilient to perturbations in the sense that it recovers the optimal point with less energy compared to GD for sufficiently small stepsizes. We will also demonstrate in our numerical experiments that the framework we propose is competitive in practice with the existing state-of-the-art algorithms from the literature and can outperform them in some problems, illustrating the potential of the proposed framework in practice. In fact, in a companion paper, we use our framework to develop a universally optimal multi-stage stochastic gradient algorithm for stochastic optimization which achieves the lower bounds without assuming a known bound for suboptimality or the variance of the gradient noise.

(see e.g. ) where the gradient $\nabla f$ is represented as a column vector. The ratio $\kappa\triangleq\frac{L}{\mu}$ is called the condition number of $f$. In many places, we also use the following relation for strongly convex smooth functions.

For our subsequent analysis, we represent the preceding relation in matrix form:

Optimization Algorithms as Dynamical Systems

Our goal is to design first-order algorithms with certain rate-robustness balance to solve

when the gradient $\nabla f$ is corrupted by random errors in the form of additive white noise. We denote the unique optimal solution of problem (3) by $x^{*}$. We will focus on Gradient Descent (GD) and Accelerated Gradient Descent (AG) and show how the parameters of these algorithms can be tuned to optimize various performance metrics.

Our analysis builds on a dynamical system representation of these algorithms. A discrete-time dynamical system with a feedback rule $\phi$ can be expressed as

which can be cast as (4) by setting $\xi_{k}=x_{k}$, $\phi(\cdot)=\nabla f(\cdot)$ and letting

On the other hand, when implemented on (3), the AG method with constant stepsize $\alpha>0$ and momentum parameter $\beta>0$ generates the iterates as follows for $k\geq 0$:

Setting $\phi(\cdot)=\nabla f(\cdot)$ and defining the state vector $\xi_{k}=\begin{bmatrix}x_{k}^{\top}&x_{k-1}^{\top}\end{bmatrix}^{\top}$, AG iterations can be rewritten as in (4) for

For both algorithms, the iterates $x_{k}$ are captured by the state $\xi_{k}$ of the dynamical system representation.

where $A$, $B$, $C$, and $D$ are selected according to (6) for GD or (8) for AG.Although our focus in this paper will be primarily on GD and AG dynamics under noise, it will be clear from our discussion that our ideas naturally extend to many other algorithms that admit such a dynamical system representation including the heavy-ball and the robust momentum methods . Except for Section 5 where we study deterministic perturbations, we assume throughout this paper that the sequence $\{w_{k}\}_{k}$ of random variables satisfies the following assumption.

It is worth emphasizing that robustness can also be studied in the solution space. Indeed, let $\{x_{k}\}_{k\geq 0}$ be a random iterate sequence corresponding to (9) where $\{w_{k}\}_{k}$ models the additive noise sequence and satisfies Assumption 2.1. Due to the noise injected at each step, the sequence $\{x_{k}\}$ will oscillate around the optimal solution with a non-zero variance; therefore, another natural metric to measure robustness is the worst-case limiting distance to the optimal solution $x^{*}$ along all possible iterate subsequences, i.e.,

Quadratic Functions

where $x^{*}=Q^{-1}p$ is the optimal solution to problem (3). Plugging the formula for the gradient $\nabla f(y_{k})$ from (12) into (9), we obtain

where $A_{Q}$ is the state-transition matrix given by $A_{Q}=A+BQC$.

Recall the robustness definition given in (10). In the next lemma, we focus on the suboptimality sequence, $\{f(x_{k})-f^{*}\}_{k}$ for quadratic $f$ and we show that the limit,

exists; moreover, for some $\{\varepsilon_{k}\}_{k}\subset[0,\infty)$ such that $\lim_{k\to\infty}\varepsilon_{k}=0$, we have

The limit in (16) can be evaluated by using the tools from standard $H_{2}$ theory arising in robust control of dynamical systems (see e.g. ) as we shall explain below. The $H_{2}$-norm is a well-known fundamental metric for quantifying the robustness of a linear dynamical system to noise in control engineering and has been widely used in designing the parameters of control systems subject to noise. Given arbitrary matrices $(A,B,C)$ and $D=0_{d}$, consider a linear system as in (4) but without feedback $\phi$. Suppose there exists $\xi^{*}$ and $y^{*}$ such that $\xi^{*}=A\xi^{*}$ and $y^{*}=C\xi^{*}$. The $H_{2}$-norm of this linear system, denoted by $H_{2}(A,B,C)$, measures the stationary variance of the output response $\{y_{k}\}$ to unit white noise input , i.e.,

The $H_{2}$ norm admits alternative definitions, which are all equivalent for linear systems (see e.g. ). When it is clear from the context, we will remove the dependency of the $H_{2}$ norm to the system matrices $(A,B,C)$. The $H_{2}$-norm can be computed as

(see e.g. ). Moreover, if $BB^{\top}$ is positive definite and $A$ is discrete-time stable (i.e., $\rho(A)<1$), the solution admits the following formula:

(see e.g. ). We will show in the following lemma that the limit $\mathcal{J}$ in (16) exists for quadratic objectives. Our proof technique is based on relating $\mathcal{J}$ to the $H_{2}$ norm of a transformed linear system as follows: We first rewrite the suboptimality $f(x_{k})-f^{*}$ in terms of the iterates $x_{k}$:

where we used the fact that $x_{k}=T\xi_{k}$, and $\frac{1}{2}Q=R^{\top}R$ is the Cholesky decomposition of $\frac{1}{2}Q$. If we consider the system defined by matrices $(A_{Q},B,RT)$, it follows from the definition of the $H_{2}$ norm (18) and (22) that

holds for some explicitly given positive constant $\psi_{0}$ that depends on the initialization $x_{0}$. Furthermore, when $A_{Q}$ is symmetric, $\varepsilon_{k}=0$ for every $k\geq 0$.

where the last equality comes from recursively using the first equality. This implies

where $\|.\|$ is the spectral norm, and the first inequality in (28) follows from the Von Neumann’s trace inequality which states that for any two $m\times m$ matrices $U$ and $V$ with singular values $||U\|_{2}=u_{1}\geq...\geq u_{m}$ and $||V\|_{2}=v_{1}\geq...\geq v_{m}$, respectively, we have$|\operatorname{Tr}(UV)|\leq\sum_{i=1}^{m}u_{i}v_{i}$. Finally, it follows from the Gelfand’s formula that there exists a sequence of non-negative numbers $\{\varepsilon_{k}\}_{k}$ such that for every $k\geq 0$, $\|A_{Q}^{k}\|_{2}\leq\left(\rho(A_{Q})+\varepsilon_{k}\right)^{k}$ and $\lim_{k}\varepsilon_{k}=0$. Note that when $A_{Q}$ is symmetric, we have $\|A_{Q}^{k}\|_{2}=\rho(A_{Q})^{k}$ so that we can choose $\varepsilon_{k}=0$. Inserting this bound into (28), we obtain the desired result. ∎

2 Gradient descent (GD) method

The dynamical system representation of GD, choosing the $A,B,C$ as in (6) yields

As shown in Lemma 3.1, the convergence rate of GD is given by $\rho(A_{Q})^{2}$. For GD, we will suppress the dependence of $\rho(A_{Q})$ on $A_{Q}$ and use the notation $\rho(\alpha)$ to highlight the effect of the stepsize $\alpha$. Since $A_{Q}$ is symmetric, $\rho(\alpha)$ can be computed as

$\alpha\in(0,2/L)$ is a necessary condition for global linear convergence; otherwise, $\rho(\alpha)\geq 1$. In particular, it is well-known that the fastest rate is achieved for the stepsize

which leads to a convergence rate of $\bar{\rho}=1-\frac{2}{\kappa+1}$. The choice of the stepsize not only affects the rate (see (29)) but also the robustness of the GD algorithm to gradient noise. The following proposition provides an analytical characterization of the robustness $\mathcal{J}$ of the GD method as a function of the stepsize, which we denote by $\mathcal{J}(\alpha)$ to highlight its dependence on $\alpha$.

Let $f$ be a quadratic function of the form $f(x)=\tfrac{1}{2}x^{\top}Qx-p^{\top}x+r$. Consider the GD iterations given by (5) with constant stepsize $\alpha\in(0,2/L)$. Then the robustness of the GD method is given by

where $0<\mu=\lambda_{1}\leq\lambda_{2}\leq...\lambda_{d}=L$ are the eigenvalues of $Q$.

We first show that without loss of generality we can assume $Q$ is a diagonal matrix. Let $Q=U\Lambda U^{\top}$ be the eigenvalue decomposition of $Q$ where $U$ is a unitary matrix and $\Lambda=\mathbf{diag}(\lambda_{1},...,\lambda_{d})$ is a diagonal matrix containing the eigenvalues of $Q$. Multiplying $A_{Q}$ by $U^{\top}$ and $U$ from left and right leads to

where $A_{\Lambda}\triangleq I_{d}-\alpha\Lambda$ is a diagonal matrix. Similarly, we multiply the Lyapunov equation (24) from left and right by $U^{\top}$ and $U$, which yields $U^{\top}A_{Q}XA_{Q}^{\top}U-U^{\top}XU+\alpha^{2}I_{d}=0$, where we have used the fact that $B=-\alpha I_{d}$ for the dynamical system representation of the GD method (see (6)). It follows from (32) that $A_{Q}=U(I_{d}-\alpha\Lambda)U^{\top}$, which when plugged into the Lyapunov equation above, yields $(I_{d}-\alpha\Lambda)U^{\top}XU(I_{d}-\alpha\Lambda)-U^{\top}XU+\alpha^{2}I_{d}=0$. This means that the matrix $U^{\top}XU$ solves the Lyapunov equation obtained by replacing $A_{Q}$ by $A_{\Lambda}$ in (24). Furthermore, the Cholesky decomposition of $\frac{1}{2}\Lambda$ is equal to $\sqrt{\frac{1}{2}}\Lambda^{1/2}$; thus, the robustness $\mathcal{J}$, corresponding to $H_{2}^{2}(A_{\Lambda},B,\sqrt{\frac{1}{2}}\Lambda^{1/2}T)$, is equal to

where we used $T=I_{d}$ for GD for the first equality and the fact that the Cholesky decomposition of $\frac{1}{2}Q$ is $(\sqrt{\frac{1}{2}}\Lambda^{1/2}U^{\top})^{\top}(\sqrt{\frac{1}{2}}\Lambda^{1/2}U^{\top})$ to obtain the second equality. Therefore, robustness $\mathcal{J}$ would be invariant if we were to replace $Q$ by $\Lambda$ and solve the Lyapunov equation (24) for $(A_{\Lambda},B)$ instead of $(A_{Q},B)$. With this replacement, it is easy to verify that the solution of the Lyapunov equation is $X_{\Lambda}=\mathbf{diag}\left(\frac{\alpha^{2}}{1-(1-\alpha\lambda_{1})^{2}},...,\frac{\alpha^{2}}{1-(1-\alpha\lambda_{d})^{2}})\right)$ as $A_{\Lambda}$ and $B$ are both diagonal. Plugging this solution into $\frac{1}{2}\operatorname{Tr}(\Lambda^{1/2}X_{\Lambda}\Lambda^{1/2})$ implies $\mathcal{J}(\alpha)=\sum_{i=1}^{d}\frac{\alpha^{2}\lambda_{i}}{2(1-(1-\alpha\lambda_{i})^{2})}=\alpha\sum_{i=1}^{d}\frac{1}{2(2-\alpha\lambda_{i})}$ which completes the proof. ∎

Proposition 3.2 also shows that the robustness $\mathcal{J}(\alpha)$ for the GD method is an increasing function of $\alpha$. This means choosing a smaller stepsize leads to GD being more robust which has been previously observed in the literature for both additive and multiplicative deterministic noise .

Having explicit expressions for both convergence rate and robustness for GD (see (29) and (31)), given an allowable deviation $\epsilon>0$ from the optimal convergence rate $\bar{\rho}=1-\frac{2}{\kappa+1}$, a natural approach to account for the trade-off between these two measures is to choose the stepsize $\alpha$ that results in the most robust algorithm satisfying the rate constraints, i.e., optimizing

This problem is equivalent to the following convex problem for $\epsilon\in[0,\frac{2}{\kappa-1})$ (which ensures that the upper bound on the rate is less than one and the optimization problem (33) admits a solution):

Indeed, $1/(1-\rho^{2})$ is a nondecreasing convex function for $\rho\in(0,1)$ and $\rho(\alpha)$ is convex in $\alpha$; therefore, both $1/(1-\rho(\alpha)^{2})$ and $\mathcal{J}(\alpha)$ in (31) are convex for $\alpha\in(0,\frac{2}{L})$ and is increasing in $\alpha$. Moreover, (34) satisfies the Slater condition. Thus, strong duality implies that there exists $\tau$ (which is a function of $\epsilon$) such that the above minimization problem is equivalent to the following unconstrained problem:

The parameter $\tau>0$ determines the trade-off between rate and robustness. For small $\tau$, the dominant term in the cost would be $\mathcal{J}(\alpha)$ so that we expect the optimal stepsize to be small since $\mathcal{J}(\alpha)$ is an increasing function of $\alpha$. On the other hand, for large enough $\tau$, the convergence rate is the dominant term in the cost; therefore, one would expect the optimal stepsize (that solves the problem (35)) to be close to $\bar{\alpha}$ which corresponds to the fastest achievable rate $\bar{\rho}$ (see (30)). In order to get more intuition about the effect of the choice of the stepsize parameter, we next give an illustrative example in dimension $d=2$ to show the behavior of the optimal $\alpha_{*}(\tau)$ as the tradeoff parameter $\tau$ is varied from zero to infinity. For computational tractability, we consider the unconstrained version of the problem given in (35).In Proposition A.1 of the appendix, we derive the first-order conditions for $\alpha_{*}(\tau)$ that allows it to be computed up to an arbitrary accuracy.

In dimension $d=2$, let $\tau=2$ and consider the parameters

The first-order optimality conditions for (35) is derived in Proposition A.1 which is equivalent to a polynomial root finding problem in $\alpha$ for a polynomial of degree $4$. The roots of polynomials can be found up to arbitrary accuracy by calculating the eigenvalues of the corresponding companion matrix , for instance using the roots function in Matlab. After a careful examination of all the roots, we conclude that the optimal stepsize $\alpha_{*}$ that minimizes the cost $F_{\tau}(\alpha)$ is $\alpha_{*}\approx 1.5055$ which gives the rate $\rho(\alpha_{*})\approx 0.8494$ and robustness $\mathcal{J}(\alpha_{*})\approx 1.9294$. This point is marked on Figure 1 below which shows the robustness level as a function of the optimal convergence rate $\rho$ when we change $\tau$ from zero (corresponds to the rightmost point in the curve) to infinity (corresponds to the uppermost point in the curve) for the parameters in (36).

The left and the middle panels of Figure 1 show the convergence rate and robustness corresponding to the optimal stepsize $\alpha_{*}$ as a function of the trade-off parameter $\tau$. As $\tau$ goes to , the robustness term is more dominant which requires a smaller stepsize; therefore, $\alpha_{*}$ goes to and $\rho(\alpha_{*})$ thus goes to $1$. As $\tau$ becomes larger, convergence rate becomes more important, and the stepsize also becomes larger to ensure faster convergence. In particular, as $\tau$ goes to infinity, $\alpha_{*}$ goes to $\bar{\alpha}$ given in (30)leading to the fastest rate, $\bar{\rho}=\frac{\kappa-1}{\kappa+1}\approx 0.8182$.

Finally, the rightmost panel of Figure 1 illustrates the trade-off between the rate and robustness. We se that for small $\tau$, the optimal stepsize $\alpha_{*}$ is smaller which implies improved robustness but slower convergence. As $\tau$ grows, we achieve faster rate at the expense of being less robust to the additive gradient noise. In addition, the points corresponding to the fastest rate, i.e., $\alpha=2/(\mu+L)$, and standard parameter choice $\alpha=1/L$ for GD has been marked on this trade-off curve.

We see from Figure 1 that smaller values of $\rho$ (or equivalently smaller values of $\frac{1}{1-\rho^{2}}$) are accompanied by larger values of $\mathcal{J}$. This suggests that the product $\mathcal{J}\frac{1}{1-\rho^{2}}$ cannot be too small for any choice of the stepsize $\alpha$. The next lemma shows that there are some fundamental limits (lower bounds) on how robust the GD can be.

Let $\rho(\alpha)$ and $\mathcal{J}(\alpha)$ be given by (29) and (31), respectively. Then, the following inequality holds \mathcal{J}(\alpha)\geq\big{(}1-\rho^{2}(\alpha)\big{)}\sum_{i=1}^{d}\frac{1}{8\lambda_{i}} for any choice of the stepsize $\alpha>0$.

It follows from (29) that for every $i\in\{1,...,d\}$, we have $\rho(\alpha)\geq|1-\alpha\lambda_{i}|$. This implies that $\tfrac{1}{1-\rho(\alpha)^{2}}\geq\tfrac{1}{1-(1-\alpha\lambda_{i})^{2}}$. Multiplying both sides by $\frac{\alpha^{2}\lambda_{i}}{2(1-(1-\alpha\lambda_{i})^{2})}$ and summing over all $i$ yields $\tfrac{1}{1-\rho^{2}}\sum_{i=1}^{d}\tfrac{\alpha^{2}\lambda_{i}}{2(1-(1-\alpha\lambda_{i})^{2})}\geq\sum_{i=1}^{d}\tfrac{\alpha^{2}\lambda_{i}}{2(1-(1-\alpha\lambda_{i})^{2})^{2}}$. Given the explicit characterization of $\mathcal{J}(\alpha)$ in Proposition 3.2 (see (31)) we obtain

The right hand side of (37) admits a lower bound as follows:

where the last inequality follows from the fact that $|2-\alpha\lambda_{i}|\leq 2$. Using the lower bound (38) along with (37) completes the proof. ∎

3 Accelerated gradient (AG) method

The dynamical system representation of AG, given $A,B,C$ in (8) leads to

We will first formulate an analogous problem to (35) for the AG method to design the parameters $(\alpha,\beta)$ in a way to find a trade-off between the rate and the robustness. Because AG has the pair $(\alpha,\beta)$ as design parameters, the analogue of (35) is

where $\mathcal{J}(\alpha,\beta)$ is the robustness to the noise for the system (14), $\rho(\alpha,\beta)$ is the convergence rate of AG with parameters $(\alpha,\beta)$ and $\mathcal{S}$ is the set of all possible choices of the tuple $(\alpha,\beta)$ so that the AG iterations are globally convergent, i.e.,

We call $\mathcal{S}$, the stability region of $AG$, in analogy with the stability region of numerical methods that arise in the discretization of continuous-time differential equations.

We next provide an explicit characterization for the convergence rate and robustness of AG for any given parameters $(\alpha,\beta)\in\mathcal{S}$. The convergence rate $\rho$ of the AG method as a function of $\alpha$ and $\beta$ is well-known. Diagonalizing the $A_{Q}$ matrix using the eigenvalue decomposition of $Q$, it can be shown after some computations that the rate $\rho=\rho(\alpha,\beta)$ admits the following formula

where $A_{Q}$ is defined by (39) and $\rho_{\lambda}$ is defined for $\lambda\in\{\mu,L\}$ as follows:

(see e.g. [33, Appendix A], [38, Section 4.3]). The explicit expression (42) for the rate allows us to characterize the set $\mathcal{S}$ in the next proposition whose proof can be found in the appendix. We illustrate the set $\mathcal{S}$ in Figure 2 for different choices of the parameters $\mu$ and $L$.We note that the stability region of a second-order difference equation that arises in accelerated algorithms that are sublinearly convergent for weakly convex quadratic functions has been studied in , however these results do not apply to the set $\mathcal{S}$ as we do not require the rate to be accelerated (we consider not only accelerated rates but also any rate $\rho$ less than one) and we consider strongly convex functions instead of weakly convex functions.

Let $\mathcal{S}$ be the stability set of Nesterov’s accelerated method defined by (41). Then its closure is given by the union of the following three sets:

with the convention that $\mathcal{S}_{3}$ is the empty set if $\mu<\frac{L}{2}$.

The next proposition gives a characterization of the robustness $\mathcal{J}(\alpha,\beta)$ of AG whose proof can be found in the Appendix C.

Let $f$ be a quadratic function of the form $f(x)=\tfrac{1}{2}x^{\top}Qx-p^{\top}x+r$. Consider the AG iterations given by (7) with parameters $(\alpha,\beta)\in\mathcal{S}$ . Then the robustness of the AG method is given by

where $\mu=\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{d}=L$ are the eigenvalues of $Q$ and

In the special case, choosing $\beta=0$ reduces to the formula (31) derived for GD.

Since we have an exact characterization of $\mathcal{J}(\alpha,\beta)$, we can derive the optimality conditions for the problem (40) by an approach similar to Proposition A.1, where the optimizer can be characterized as a root of some polynomial. In dimension $d=2$, given parameters $\mu$ and $L$, the optimizer is easy to compute. However, in high dimensions, this is computationally expensive as it would require determining all the eigenvalues of $Q$ which can be as expensive as optimizing the objective function $f$. Nevertheless, exploiting the convexity properties of the function $u_{\alpha,\beta}(\lambda)$, we develop a tractable upper bound for $\mathcal{J}(\alpha,\beta)$ that only depends on $\mu$ and $L$, hence tractable. Moreover, in the numerical experiments section, we present experiments illustrating that this approach can lead to good performance in terms of trading the speed and the robustness of an algorithm.

To develop this upper bound, first, we show in Lemma D.1 that the function $u_{\alpha,\beta}(\lambda)$ defined in (46) is convex in $\lambda\in[\mu,L]$ for fixed $(\alpha,\beta)\in\mathcal{S}$. Therefore, its maximum is attained at one of the endpoints of this interval, i.e.,

Substituting this upper bound in (45) and (40) leads to

This objective only depends on $\mu$ and $L$ and is differentiable everywhere in the interior of the stability region $\mathcal{S}$ except when the first term is not differentiable, i.e., when $u_{\alpha,\beta}(\mu)=u_{\alpha,\beta}(L)$, or the second term is not differentiable, i.e., when $\rho_{\mu}=\rho_{L}$ or $\Delta_{\mu}=0$ or $\Delta_{L}=0$. Furthermore, following a similar approach as in Example 3.4, the first order optimality conditions with respect to $\alpha$ and $\beta$ results in low-order polynomials (that are independent of the dimension $d$) which can be solved efficiently up to any accuracy. Thus, other than checking the non-differentiable points of $\bar{F}_{\tau}$, the bottleneck in computational complexity is determined by computing the roots of a polynomial with a small degree (whose degree is independent from the dimension $d$), which is easy to compute even in high dimensions.

Strongly Convex Functions

The goal is to extend the definitions of rate and robustness from the quadratic case to general strongly convex functions. We will use

(provided also in (10)) to define the robustness of an algorithm and study the convergence rate of the expected suboptimality to an interval around zero with radius $\sigma^{2}\mathcal{J}$. For both GD and AG, our main results provide upper bounds of the form:

where $\psi_{0}$, $R$, and $0<\rho<1$ are non-negative numbers all of which depend on algorithm parameters and the initial point $x_{0}$. Clearly, $R$ is an upper bound on $\mathcal{J}$; we will show in this section that our bounds are tight. Moreover, we also recover the fastest known rates in the literature in the absence of noise ($\sigma$=0). Our upper bounds only depend on $\mu$ and $L$, and are computationally tractable and explicit in some cases. With these upper bounds, one can formulate an optimization problem similar to that of the previous section to find the algorithm parameters that can achieve a particular trade-off between rate and robustness.

2 Rate and robustness trade-off analysis using Lyapunov functions

We use a Lyapunov function approach to provide a bound as in (51) for both GD and AG methods. In particular, we consider a family of Lyapunov functions parameterized by a non-negative constant $c$ and a positive semidefinite matrix $P$ as

where $V_{P}(\xi)\triangleq(\xi-\xi^{*})^{\top}P(\xi-\xi^{*})$, and study the change in the Lyapunov function $V_{P,c}(\xi)$ along $\{\xi_{k}\}_{k}$ generated by the dynamical system representation (49).

Consider the Lyapunov function $V_{P}(\xi)=(\xi-\xi^{*})^{\top}P(\xi-\xi^{*})$ where $P\succeq 0$. Then, we have

for every $k\geq 0$; hence, it follows from (57) and (58) along with (55) that

This MI based approach has been used in the literature to study the convergence rate of first-order methods, e.g., . Here we use it to characterize their rate and robustness under additive gradient noise.

3 Gradient descent (GD) method for strongly convex functions

which admits the dynamical system representation in (49) with $A,B,C$ as in (6). The next theorem extends the result of Proposition 3.2 to general strongly convex functions and characterize the behavior of $\{x_{k}\}_{k}$ under additive gradient error.

holds where $X_{0}=\begin{bmatrix}2\mu LI_{d}&-(\mu+L)I_{d}\\ -(\mu+L)I_{d}&2I_{d}\end{bmatrix}$ and $P=pI_{d}$. Then for all $k\geq 0$:

Noting that $\xi_{k}=y_{k}$ for GD, it follows from (2) with $x=\xi_{k}$ and $y=x^{*}$ that (58) holds for $X=X_{0}$ and $c=0$. Moreover, (61) implies that (57) holds for $X=X_{0}$ and $P=p\otimes I_{d}$; therefore, (59) yields

With $B=-\alpha I_{d}$ for GD, we have $\operatorname{Tr}(B^{\top}PB)=\alpha^{2}pd$ and this completes the proof. ∎

Note that for a fixed $\alpha$, a smaller $\rho$ makes both terms of (62) smaller as $\tfrac{1-\rho^{2k}}{1-\rho^{2}}$ is an increasing function of $\rho$. If $\alpha\in(0,2/L)$, it was shown in that there exist $(p,\rho)$ such that the MI in (61) holds; moreover, for a given $\alpha$ fixed, the smallest $\rho\in(0,1)$ for which such a positive $p$ exists is equal to

as in (29) given for quadratic functions. Using $\rho=\rho_{GD}(\alpha)$ in (62) leads to the following upper bound for GD. Trivially, $\mathcal{J}^{\prime}\leq\frac{\alpha^{2}d}{1-\rho_{GD}(\alpha)^{2}}$. More details are provided in Appendix E as a supplementary material.

where $\psi_{0}=\frac{L}{2}\left\|x_{0}-x^{*}\right\|^{2}$ and $\rho_{GD}(\alpha)$ is given in (64). As a consequence,

Using the fact that $f(x_{k})-f(x^{*})\leq\frac{L}{2}\|x_{k}-x^{*}\|^{2}$ for $k\geq 0$ together with Proposition 4.3 yields the desired result. ∎

Note that by substituting $\rho_{GD}$ in (65), we obtain $R_{GD}=\mathcal{O}(\alpha d)$. This bound is tight, as Proposition 3.2 implies that for quadratic functions $\mathcal{J}=\Theta(\alpha d)$.

4 Accelerated Gradient (AG) method for strongly convex functions

We next consider the AG algorithm with gradient noise given by

As before, these iterations admit the dynamical system representation in (49) with $A,B,C$ as in (8). We use the following result which extends Lemma 3 in to the case with noisy gradient.

Setting $x=x_{k}$ and $y=y_{k}$ in the second inequality in (1) leads to

Similarly, setting $x=x_{k+1}=y_{k}-\alpha\nabla f(y_{k})-\alpha w_{k}$ and $y=y_{k}$ in (1) yields to

Note that $x_{k}-y_{k}=x_{k}-((1+\beta)x_{k}-\beta x_{k-1})=\beta(x_{k-1}-x_{k})$; hence, (69) implies

Next, in a similar way, setting $x=x^{*}$ and $y=y_{k}$ in (1), and summing the second inequality with (68) leads to

Multiplying (70) by $\rho^{2}$ and (71) by $1-\rho^{2}$, and summing them will lead to the desired result. ∎

for $X_{1}$ and $X_{2}$ defined in Lemma 4.5. Then the following bounds hold for all $k\geq 0$:

Using (2) for $x=y_{k}$ and $y=x^{*}$ along with the fact $y_{k}=C\xi_{k}$ yields

This inequality along with Lemma 4.5 implies that (58) holds for $X=c_{0}X_{0}+cX(\rho)$ and $\Gamma=\tfrac{1}{2}L\alpha^{2}d$. Moreover, (72) implies that (57) holds for this $X$. Therefore, (59) holds and $\operatorname{Tr}(B^{\top}PB)=\alpha^{2}\operatorname{Tr}(P_{11})$ completes the proof. ∎

where $\psi_{0}=\frac{1}{c}V_{P,c}(\xi_{0})$. As a consequence, $\mathcal{J}\leq R_{AG}(\alpha,\beta)$.

Using this result, the next corollary characterizes the rate and robustness of the AG method with a particular parameterization.

where $\psi_{0}=V_{P,1}(\xi_{0})$, $\rho_{AG}(\alpha)\triangleq\sqrt{1-\sqrt{\alpha\mu}}$ and $R_{AG}(\alpha)\triangleq\frac{\alpha d}{2(1-\rho_{AG}(\alpha)^{2})}(1+\alpha L)$; hence, $\mathcal{J}\leq R_{AG}(\alpha){=\frac{\sqrt{\alpha}d}{2\sqrt{\mu}}(1+\alpha L)}$.

Therefore, the desired result follows from Corollary 4.7. ∎

5 Approximating the rate and robustness trade-off curve

In particular, for GD, the best robustness level while asking for linear convergence with rate $\rho_{{\rm GD},\epsilon}$ or faster is obtained by solving

where $\rho_{GD}(\alpha)$ is given in (64). The function $\rho_{GD}(\alpha)$ is convex and piecewise linear in $\alpha$ over the interval $[0,2/L]$ with a unique minimum at $\bar{\rho}_{GD}$ and it satisfies $\rho_{GD}(0)=\rho_{GD}(2/L)=1$ on the boundary points. Therefore, it follows from this property that, given $\epsilon\in(0,\tfrac{2}{\kappa-1})$, there are exactly two $\alpha_{\epsilon}>0$ values such that $\rho_{GD}(\alpha_{\epsilon})=\rho_{{\rm GD},\epsilon}$ which we can explicitly compute as $\alpha_{\epsilon}=\frac{2-\epsilon(\kappa-1)}{L+\mu}$ or $\alpha_{\epsilon}=\frac{2+\epsilon(\kappa-1)/\kappa}{L+\mu}$. The former value is strictly smaller as $\varepsilon>0$ and $\kappa>1$ here. From the formula (65), we have $R_{GD}(\alpha_{\varepsilon})=\frac{L\alpha_{\varepsilon}^{2}d}{2(1-\rho^{2}_{{\rm GD},\epsilon})}$. Clearly one should select the smaller $\alpha_{\epsilon}$ value to minimize the robustness bound, i.e., a choice of $\alpha_{\epsilon}=\frac{2-\epsilon(\kappa-1)}{L+\mu}$ leads to $\rho_{{\rm GD},\epsilon}$ rate with a robustness bound $R_{GD}(\alpha_{\epsilon})$, i.e., $\mathcal{J}_{{\rm GD},\epsilon}\leq R_{GD}(\alpha_{\epsilon})$.

For AG, we can also write an analogous optimization problem in order to trade rate with robustness:

The first approach is similar to the one we used for GD. In particular, consider Corollary 4.9, for $\alpha\in(0,1/L]$, choosing $\beta=\frac{1-\sqrt{\alpha\mu}}{1+\sqrt{\alpha\mu}}$ implies that $\rho_{AG}(\alpha)=\sqrt{1-\sqrt{\alpha\mu}}$. We get $\rho_{AG}(\alpha_{\epsilon})=\rho_{{\rm AG},\epsilon}$ for

with \epsilon\in\big{[}0,\sqrt{\frac{\sqrt{\kappa}}{\sqrt{\kappa}-1}}-1\big{)} to make sure the rate is smaller than $1$. Thus, choosing $(\alpha,\beta)=(\alpha_{\epsilon},\beta_{\epsilon})$ with $\beta_{\epsilon}\triangleq\frac{1-\sqrt{\alpha_{\epsilon}\mu}}{1+\sqrt{\alpha_{\epsilon}\mu}}$ guarantees the rate $\rho_{{\rm AG},\epsilon}$. In addition, Corollary 4.9 implies the robustness bound $R_{AG}(\alpha_{\epsilon})=\mathcal{O}(\frac{\sqrt{\alpha_{\epsilon}}d}{\sqrt{\mu}})$ for this case.

with $X_{0}$ and $X(\rho)$ defined in Corollary 4.7, and $\bar{R}_{AG}$ as given above. In fact, for a fixed $(\alpha,\beta)$, this is a small dimensional convex SDP problem and can be solved easily.

Thus, we first grid the AG parameter space, i.e., $\{(\alpha_{i_{1}},\beta_{i_{2}})\}_{i_{1}\in\mathcal{I}_{1},i_{2}\in\mathcal{I}_{2}}$ and for given trade-off parameter $\epsilon$, we solve $|\mathcal{I}_{1}||\mathcal{I}_{2}|$ many 4-dimensional SDPs, i.e., for each $(i_{1},i_{2})\in\mathcal{I}_{1}\times\mathcal{I}_{2}$,

and this bound can be achieved for some choices of $f$. For large $d$, we have clearly $\mathcal{J}_{max}(\alpha,\beta)/d\approx\max\left[u_{\alpha,\beta}(\mu),u_{\alpha,\beta}(L)\right]$. In Figure 3, we plot the latter quantity versus the convergence rate (marked in purple color) to demonstrate the rate-robustness curve for AG in the case of quadratic objective functions. We observe from Figure 3 that our bounds for the quadratic case are tighter than those for general strongly convex functions as expected.

Asymptotic stability of the optimum with respect to perturbations

Our discussion so far has focused on the robustness of first-order methods with respect to random noise in the gradients, which we quantify by $\mathcal{J}$ defined in (10). Our robustness measure $\mathcal{J}$ is based on the $H_{2}$ norm of an associated linear dynamical system. It is well known that the $H_{2}$ norm of a dynamical system is closely related to the asymptotic stability of the equilibrium (which is characterized by the optimal solution $x^{*}$ to (3) in our setup) in the sense that it quantifies how quickly the system can converge back to the equilibrium if it is unsettled from its equilibrium in the direction of a coordinate . More specifically, for each $i\in\{1,\ldots,d\}$, let $\{x_{k}^{i}\}_{k\geq 0}$ be the iterate sequence corresponding to (49) whenever $\{w_{k}\}_{k}=\delta[k]e_{i}$ for $k\geq 0$ where $e_{i}$ is the $i$-th basis vector, i.e., we perturb the system from its equilibrium with an impulse input in the direction of $e_{i}$. Let

where $\|x^{i}-x^{*}\|_{2}$ is the $l_{2}$ norm of the sequence $\{x_{k}^{i}-x^{*}\}_{k}$. It is worth noting that $\{x_{k}^{i}\}_{k}$ is the same as the iterate sequence of the noiseless system (4) with initial state $\xi^{*}+Be_{i}$, $D=0_{d}$ and $\phi(\cdot)=\nabla f(\cdot)$.

For GD, the following bound holds for all $\alpha\in(0,2/L)$

where $\rho_{GD}(\cdot)$ is defined in (64). Moreover, for AG, given $\alpha\in(0,1/L]$, setting $\beta(\alpha)=\frac{1-\sqrt{\alpha\mu}}{1+\sqrt{\alpha\mu}}$, the perturbation stability can be bounded as $\mathcal{J}_{*}(\alpha)\leq\frac{\alpha^{2}d}{\sqrt{\alpha\mu}}(1+\kappa)$.

Recall that $\{x_{k}^{i}\}_{k}$ is the same as the iterate sequence of the noiseless system (4) with initial state $\xi^{*}+Be_{i}$. Hence, Proposition 4.3 with $\sigma=0$ implies that

for some $\rho\in(0,1)$ and for any $1\leq i\leq d$ and stepsize $\alpha\in(0,2/L)$. Thus, $\sum_{k=0}^{\infty}\left\|x^{i}_{k}-x^{*}\right\|^{2}\leq\frac{1}{1-\rho^{2}}\left\|x^{i}_{0}-x^{*}\right\|^{2}$, which implies $\left\|x^{i}-x^{*}\right\|^{2}\leq\frac{\alpha^{2}}{1-\rho^{2}}$ for all $i=1,\ldots,d$ since $B=-\alpha I_{d}$ for GD. Therefore, we have $\mathcal{J}_{*}(\alpha)\leq\frac{\alpha^{2}d}{1-\rho^{2}}$. Moreover, given any stepsize $\alpha\in(0,2/L)$ for GD, using (64), which is the smallest $\rho$ value for which (88) holds, we obtain (87). On the other hand, for AG, using Corollary 4.9 with $\sigma=0$ and the fact that $\frac{\mu}{2}\left\|x_{k}-x^{*}\right\|^{2}\leq f(x^{k})-f^{*}$, we get $\left\|x^{i}_{k}-x^{*}\right\|^{2}\leq\rho_{AG}^{2k}(\left\|x_{0}^{i}-x^{*}\right\|^{2}+\tfrac{2}{\mu}(f(x^{i}_{0})-f^{*}))$ for $k\geq 0$ and $i=1,\ldots,d$, where we used $x_{0}^{i}=x_{-1}^{i}=x^{*}+Be_{i}$ for $i=1,\ldots,d$. Thus,

Numerical Experiments

Our first set of experiments concern a further study of Example 3.4 for comparing AG and GD in terms of performance. In the leftmost plot of Figure 4, we vary the trade-off parameter from $\tau=0$ to $\tau=\infty$ for AG and plot the robustness level $\mathcal{J}(\alpha_{*}(\tau),\beta_{*}(\tau))$ versus the rate $\rho(\alpha_{*}(\tau),\beta_{*}(\tau))$ corresponding to the optimal parameters $(\alpha_{*}(\tau),\beta_{*}(\tau))$, we also plot the analogous curve for GD (the same curve from Figure 1) to compare. We observe that for the same achievable convergence rate, the optimized AG parameters lead to more robust algorithms compared to the optimized GD algorithms as AG has an additional parameter $\beta$ to optimize robustness over. This shows that AG can improve GD in terms of both convergence rate and robustness at the same time when gradients are subject to white noise. This result is in contrast with the deterministic gradient error setting in , which shows that GD performance degrades gracefully while AG may accumulate error. Therefore, our results suggest that AG algorithms can tolerate random noise better than deterministic noise to preserve their accelerated rates, which is also inline with the theoretical findings of . Also it is interesting to note that the popular choice of parameters (blue and red dots), as well as the parameters that lead to the optimal (fastest) rate (green and purple dots) lie on curves that trade robustness with rate in an optimal fashion.

Next, we illustrate the tightness of our upper bound $\bar{\mathcal{J}}(\alpha,\beta)$ provided in (47) to the (true) robustness level $\mathcal{J}(\alpha,\beta)$. This upper bound results in a small scale optimization problem (48) that allows trading-off robustness and the convergence rate in a way that computationally tractable, even in high dimensions. The middle plot of Figure 4 shows the convergence rate and robustness obtained by solving (40) versus solving (48). The objective is a random quadratic function in dimension $d=100$ with parameters $\mu=0.1,L=1$. Our results show that for any trade-off parameter $\tau$ our upper bound is within a factor of $1.2$ of true parameters, illustrating the accuracy of this approximation to the optimal parameters for different levels of robustness, especially the approximation is more accurate when the trade-off parameter is larger (in which case the convergence rate is closer to 1). We obtain quantitatively similar results repeating this experiment with other randomly generated quadratic functions.

Next, we illustrate our framework to trade-off robustness and convergence rate on a quadratic optimization problem, similar to the one considered in where it is shown that AG algorithms with standard choice of parameters have difficulty to handle random gradient noise. We consider the quadratic function $f(x)=\frac{1}{2}x^{\top}Qx+b^{T}x+\delta\|x\|^{2}$ in dimension $d=100$ where $Q$ is the Laplacian of a cyclic graph, $\delta=0.1$ is a regularization parameter to make the problem strongly convex and $b$ is a random vector. As it can be seen in the rightmost plot of Figure 4, we show that when properly modified, AG can be both faster and more robust in comparison with GD.

In the leftmost plot of Figure 5, we compare the tuned AG with other algorithms such as AC-SA and the Flammarion-Bach algorithm . For this purpose, we consider the same quadratic test problem from in dimension $d=20$, where the eigenvalues of its Hessian $Q$ are set equal to $\lambda_{i}=i^{2}$ for $i=1,2,...,20$. Our results show that modified AG can trade robustness with the convergence rate successfully and can improve upon AC-SA and Flammarion-Bach algorithm on this example.

Conclusion

We consider the gradient descent (GD) and accelerated gradient (AG) methods for optimizing strongly convex functions. We developed a computationally tractable framework to design their parameters in a way to trade between two conflicting performance measures: the convergence rate and the robustness to additive white noise in the gradient computations measured in terms of final asymptotic variance of the algorithm output. For strongly convex quadratics, we show that this robustness measure is equal to the $H_{2}$ norm of a dynamical system associated to the optimization algorithm and give an explicit characterization of this quantity. Our results show that for the same achievable rate, AG can always be tuned to be more robust. Similarly, for the same robustness level, we show that AG can be tuned to be always faster than GD. We also give fundamental lower bounds on the achievable robustness level for gradient descent for a given achievable rate. We show how our analysis can be extended to smooth strongly convex functions and we derive upper bounds on the robustness measures for GD and AG.

Acknowledgments

The work of Necdet Serhat Aybat is partially supported by NSF Grant CMMI-1635106. Alireza Fallah is partially supported by Siebel Scholarship. Mert Gürbüzbalaban acknowledges support from the grants NSF DMS-1723085 and NSF CCF-1814888.

References

There exists an optimizer $\alpha_{*}(\tau)$ to the minimization problem (35). Furthermore, any optimizer is either $\alpha_{*}(\tau)=2/(\mu+L)$ or it satisfies one of the following two conditions:

Therefore, by examining the values of $F$ at the points that satisfy this equality and inequality constraints, we can determine the optimal stepsize $\alpha^{*}(\tau)$.

The optimal $\alpha^{*}$ cannot be attained on the boundary points of the interval $[0,2/L]$ as $F$ is not finite at these points. Therefore, it suffices to solve the optimization problem over the open interval $(0,2/L)$ where $F$ is differentiable with respect to $\alpha$ except when $|1-\alpha\mu|=|1-\alpha L|$, i.e. when $\alpha=2/(\mu+L)$. For $\alpha^{*}\neq 2/(\mu+L)$, we can write-down the first-order conditions of optimality $\tfrac{\partial F}{\partial\alpha}=0$ which leads to (89) and (90). ∎

Appendix B Proof of Proposition 3.6

In the light of the formula (42) that characterizes $\rho(A_{Q})$, the closure of the stability set $\mathcal{S}$ admits the representation $\mathcal{S}=\mathcal{S}_{\mu}\cap\mathcal{S}_{L}$ where for $\lambda\in\{\mu,L\}$ we define

We first write $\mathcal{S}_{\lambda}$ as a union of two disjoint sets depending on the signature of $\Delta_{\lambda}$: $\mathcal{S}_{\lambda}=\mathcal{S}_{\lambda,1}\cup\mathcal{S}_{\lambda,2}$ where

It follows from the definition of $\Delta_{\lambda}$ in (43) that $\Delta_{\lambda}\leq 0$ if and only if $0\leq 1-\alpha\lambda\leq\frac{4\beta}{(1+\beta)^{2}}$; and when this condition holds, $\rho_{\lambda}=\sqrt{\beta(1-\alpha\lambda)}\leq 1$ if and only if $0\leq 1-\alpha\lambda\leq\frac{1}{\beta}.$ Therefore,

We next focus on $\mathcal{S}_{\lambda,2}$. Note that $\Delta_{\lambda}\geq 0$ if and only if

If (94) is satisfied, then $\rho_{\lambda}\leq 1$ if and only if $\frac{1}{2}(1+\beta)(1-\alpha\lambda)\mbox{sign}(1-\alpha\lambda)+\frac{1}{2}\sqrt{\Delta_{\lambda}}\leq 1$. There are two cases:

1) $\Delta_{\lambda}>0$ and $1-\alpha\lambda<0$: In this case, $\rho_{\lambda}\leq 1$ if and only if $\sqrt{\Delta_{\lambda}}\leq 2-(1+\beta)c_{\lambda}$, where $c_{\lambda}=-(1-\alpha\lambda)>0$. By squaring both sides, this is if and only if, $\Delta_{\lambda}\leq(2-(1+\beta)c_{\lambda})^{2}\quad\mbox{and}\quad 2-(1+\beta)c_{\lambda}\geq 0$ The first inequality holds if $c_{\lambda}=-(1-\alpha\lambda)\leq\tfrac{1}{2\beta+1}$ whereas the second inequality holds if $c_{\lambda}=-(1-\alpha\lambda)\leq\tfrac{2}{\beta+1}.$ The first inequality is more binding, if it holds the second inequality holds too. Therefore,

2) $\Delta_{\lambda}>0$ and $1-\alpha\lambda>0$: In this case, $\rho_{\lambda}\leq 1$ if and only if $\sqrt{\Delta_{\lambda}}\leq 2-(1+\beta)d_{\lambda}$ where $d_{\lambda}:=-c_{\lambda}=(1-\alpha\lambda)>0$. After squaring both sides, this is if and only if

where the first inequality simplifies to $1\geq d_{\lambda}$. (96) along with (94) means $\tfrac{4\beta}{(1+\beta)^{2}}\leq 1-\alpha\lambda\leq\min\{1,\tfrac{2}{1+\beta}\}$ which implies $\beta\leq 1$; therefore,

To complete the proof, due to the representation (91), it suffices to compute the intersection $\mathcal{S}_{\mu}\cap\mathcal{S}_{L}$. There are several cases to consider depending on the value of $\alpha$:

1) First, consider $\alpha\in[0,\frac{1}{L}]$. In this case $1-\alpha\mu\geq 1-\alpha L\geq 0$, and hence (98) implies $1-\alpha\mu\leq\frac{2}{1+\beta}$ if $\beta\leq 1$ whereas $1-\alpha\mu\leq\frac{1}{\beta}$ if $\beta\geq 1$. Nevertheless, if $\beta\leq 1$ then $\frac{2}{1+\beta}\geq 1$, so the first case always holds; hence, $(1-\alpha\mu)\beta\leq 1$.

2) Now, assume $\alpha\in[\frac{1}{L},\min\{\frac{2}{L},\frac{1}{\mu}\}]$. Then $1-\alpha\mu\geq 0\geq 1-\alpha L$, and thus (98) yields $1-\alpha L\geq-\frac{1}{1+2\beta},\quad 1-\alpha\mu\leq\min\{\frac{1}{\beta},\frac{2}{1+\beta}\}$ where the second inequality again simplifies to $(1-\alpha\mu)\beta\leq 1$.

3) The last possible case happens when $\mu\geq\frac{L}{2}$ , and so $\alpha\in[\frac{1}{\mu},\frac{2}{L}]$ is possible. In this case $1-\alpha L\leq 1-\alpha\mu\leq 0$, and so using (98), we just need to check $1-\alpha L\geq-\frac{1}{1+2\beta}$ Considering all these cases along with the fact that (98) shows $\alpha$ cannot be greater than $\frac{2}{L}$ completes the proof.

Appendix C Proof of Proposition 3.7

Similar to the analysis for GD, we can assume without loss of generality that $Q$ is diagonal. The proof is also similar. Consider $U\Lambda U^{\top}$ be the eigenvalue decomposition of $Q$. Then $A_{Q}$ in (39) can be written as

Replacing $A_{Q}$ from (99) in Lyapunov equation (20) implies

Let $P_{\pi}$ be the permutation matrix associated with the permutation $\pi$ over the set $\{1,2,...,2d\}$ that satisfies $\pi(i)=2i-1$ for $1\leq i\leq d$ and $\pi(i)=2(i-d)$ for $d+1\leq i\leq 2d$. It is well-known that permutation matrices satisfy $P_{\pi}^{-1}=P_{\pi}^{\top}=P_{\pi^{-1}}$; therefore, multiplying Lyapunov equation (20) by $P_{\pi}$ and $P_{\pi}^{\top}$ from left and right, respectively, leads to

where $Y=P_{\pi}XP_{\pi}^{\top}$. It follows from (39) that

and $0<\mu=\lambda_{1}\leq\lambda_{2}\leq...\leq\lambda_{d}=L$ are the eigenvalues of $Q$. Since $BB^{T}=\begin{bmatrix}\alpha^{2}I_{d}&0_{d}\\ 0_{d}&0_{d}\end{bmatrix}$,$P_{\pi}BB^{\top}P_{\pi}^{\top}$ is a $2d$ by $2d$ diagonal matrix with $\alpha^{2}$ on entries $(1,1),(3,3),...,$ $(2d-1,2d-1)$ and zero elsewhere. Hence, $Y$ that solves (103) is a block diagonal matrix in the form: $Y=\mathbf{diag}([Y_{i}]_{i=1}^{d})$, where $Y_{i}=\begin{bmatrix}y_{i}^{u}&y_{i}^{o}\\ y_{i}^{o}&y_{i}^{d}\end{bmatrix}$ satisfies the equality

for all $i=1,\dots,d$. This is equivalent to the linear system:

Solving this system of equations, we obtain:

The $\mathcal{J}(\alpha,\beta)$ can be computed using