Acceleration by Random Stepsizes:
Hedging, Equalization, and the Arcsine Stepsize Schedule

The rate in Theorem 1.2 is optimal since there is an exactly matching lower bound: there is an explicit $\kappa$-conditioned quadratic function for which any (Krylov) first-order algorithm cannot contract to the optimum faster than this rate $\tfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$ [41, Theorem 2.1.12].

Theorem 1.2 first appeared as Theorem 6.1 of the 2018 MS thesis [2] of the first author (advised by the second author). This was the first convergence rate for any stepsize schedule (deterministic or randomized) that improved over the textbook rate for constant-stepsize GD, in any setting beyond quadratics. Indeed, even though many stepsize schedules have been investigated in both theory and practice, none of these alternative schedules had led to an analysis that outperforms (even by an arbitrarily small amount) the “unaccelerated” contraction factor $(\tfrac{\|x_{n}-x^{*}\|}{\|x_{0}-x^{*}\|})^{1/n}\leqslant\tfrac{\kappa-1}{\kappa+1}$ of constant stepsize GD, which corresponds to an iteration complexity of $\Theta(\kappa\log\tfrac{1}{\varepsilon})$. This comparison statement includes even adaptive stepsize schedules such as exact line search [42, 9, 45, 15], Armijo-Goldstein schedules [40], Polyak-type schedules [45], and Barzilai-Borwein-type schedules [5]. Some of these schedules fail to accelerate even on quadratics (and therefore also on separable convex problems), a classic example being exact line search [33, §3.2.2].

In the years since [2], a line of work has shown partially accelerated convergence rates for general convex optimization using deterministic stepsize schedules [2, 12, 23, 13, 26, 3, 4, 27, 28, 58, 62, 29, 8]. In particular, [3] showed the Silver Convergence Rate $O(\kappa^{\log_{1+\sqrt{2}}2}\log\tfrac{1}{\varepsilon})\approx O(\kappa^{0.78}\log\tfrac{1}{\varepsilon})$ for minimizing $\kappa$-conditioned functions by using a certain fractal-like deterministic stepsize schedule called the Silver Stepsize Schedule. This rate is conjecturally optimal among all possible deterministic stepsize schedules [3] and naturally extends to non-strongly convex optimization [4]. Theorem 1.2 is incomparable: it shows a fully accelerated rate $O(\kappa^{1/2}\log\tfrac{1}{\varepsilon})$, but under the additional assumption of separability. This incomparability raises several interesting questions about the possible benefit of randomization for GD. The future work section §7 collects several such conjectures.

Starting from the seminal work of Nesterov in 1983 [41], the traditional approach for obtaining accelerated convergence rates is to modify the GD algorithm beyond just changing stepsizes, for example by adding momentum, extra state, or other internal dynamics. These results hold for arbitrary convex optimization and do not require separability. See for example [11, 21, 52, 56] for some recent such algorithms, and see the recent survey [22] for a comprehensive account of this extensive line of work. The purpose of this paper is to investigate an alternative mechanism for acceleration: can GD be accelerated in its original form, by just choosing random stepsizes?

The assumptions in Theorem 1.2 can be relaxed. Separability can be generalized to radial-separability, and strong convexity and smoothness can be relaxed to less stringent assumptions on the spectrum of the Hessian of $f$; see §4 for details. Note also that since GD does not depend on the choice of the basis $U$ in the separable decomposition (1.5), Theorem 1.2 only requires the existence of some such $U$ (which is the definition of separability). In particular, the algorithm does not need to know $U$. Also, Theorem 1.2 is stated for convergence in terms of distance to optimality, but extends to all standard desiderata; see §4. Theorem 1.2 also extends to settings where gradient descent is performed inexactly, e.g., due to using inexact gradients; see §6.1.

Randomized algorithms have variability in their trajectories. This provides potential opportunities for running multiple times in parallel and selecting the best run, see §5.3. This best run can be significantly better than the typical run (whose rate converges almost surely to the accelerated rate, see §5.2) and the worst run (whose rate can diverge but only with exponentially small probability, see §5.1).

This paper considers random stepsizes in order to accelerate GD. This randomized perspective is also helpful for the dual problem of constructing lower bounds, i.e., exhibiting functions for which gradient descent converges slowly. In particular, lifting to distributions over functions leads to a game between an algorithm (who chooses a distribution over stepsizes) and an adversary (who chooses a distribution over functions) which is symmetric in that the Arcsine distribution is the optimal answer for both. This leads to new perspectives on classical lower bounds for GD; details in §6.2.

A distinguishing feature of our analysis is a connection to potential theory: the variational characterization for the optimal stepsize distribution mirrors the variational characterization for a certain equilibrium distribution which configures a unit of charged particles so as to minimize the logarithmic potential energy. The optimal distribution satisfies a remarkable “equalization property” which in the physical context amounts to a constant potential over space, and in the optimization context amounts to an identical convergence rate over all quadratic functions. A key technical insight is that martingale arguments extend this phenomenon to all separable convex functions. While links between Krylov-subspace algorithms and potential theory are classical (see the excellent survey [19]), these connections are restricted to the special case of quadratics and moreover, even for that setting, these connections have not been exploited for the analysis of GD with random stepsizes. See §2 for a conceptual overview of our analysis techniques and how we exploit this equalization property to prove Theorem 1.2.

Without loss of generality after translating,

$\displaystyle f(x)=\frac{1}{2}x^{T}Hx\,\qquad\text{ where }\qquad mI\preceq H\preceq MI\,.$

(2.1)

Since the function is quadratic, its gradient is linear $\nabla f(x)=Hx$, hence GD is a linear map

$\displaystyle x_{t+1}=x_{t}-\alpha_{t}\nabla f(x_{t})=(I-\alpha_{t}H)x_{t}\,,$

(2.2)

and therefore the final iterate is

$\displaystyle x_{n}=\prod_{t=0}^{n-1}(I-\alpha_{t}H)x_{0}\,.$

(2.3)

The convergence rate for a given stepsize schedule—in the worst case over objectives functions in $\mathcal{Q}$ and initializations $x_{0}$—is therefore

	$\displaystyle\sup_{f\in\mathcal{Q},\,x_{0}\neq x^{*}}\left(\frac{\|x_{n}-x^{*}\|}{\|x_{0}-x^{*}\|}\right)^{1/n}$	$\displaystyle=\sup_{mI\preceq H\preceq MI,\,x_{0}\neq x^{*}}\left(\frac{\|\prod_{t=0}^{n-1}(I-\alpha_{t}H)x_{0}\|}{\|x_{0}\|}\right)^{1/n}$
		$\displaystyle=\sup_{mI\preceq H\preceq MI}\|\prod_{t=0}^{n-1}(I-\alpha_{t}H)\|^{1/n}$		(2.4)
		$\displaystyle=\sup_{\lambda\in[m,M]}\prod_{t=0}^{n-1}|1-\alpha_{t}\lambda|^{1/n}\,.$		(2.5)

Above, the penultimate step is by definition of the spectral norm, and the final step is by diagonalizing. It is a simple and celebrated fact, dating back to Young in 1953 [61], that minimizing this convergence rate (2.5) over stepsizes is equivalent to a certain extremal polynomial problem, and that as a consequence, the optimal stepsizes minimizing (2.5) are the inverses of the roots of a (translated and scaled) Chebyshev polynomial, in any order. Said more precisely: for any horizon length $n$, the optimal $n$ stepsizes are any permutation of

$\displaystyle\left\{\left(\frac{M+m}{2}+\frac{M-m}{2}\cos\left(\frac{2t+1}{2n}\pi\right)\right)^{-1}\right\}_{t=0}^{n-1}$

(2.6)

and the corresponding $n$-step convergence rate in (2.5) is

$\displaystyle\sup_{f\in\mathcal{Q},\,x_{0}\neq x^{*}}\left(\frac{\|x_{n}-x^{*}\|}{\|x_{0}-x^{*}\|}\right)^{1/n}=\left[2\frac{(\sqrt{\kappa}+1)^{n}(\sqrt{\kappa}-1)^{n}}{(\sqrt{\kappa+1})^{2n}+(\sqrt{\kappa-1}^{2n})}\right]^{1/n}\overset{n\to\infty}{\longrightarrow}\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\,.$

(2.7)

Much less well-known is the fact that this accelerated rate can also be asymptotically achieved by simply using i.i.d. stepsizes from an appropriate distribution [32, 47, 63]. The intuition is that the convergence rate (2.5) is invariant in the order of the stepsizes, instead depending only on their empirical measure, and as $n\to\infty$ one can achieve the same empirical measure (and thus the same convergence rate) by replacing the carefully chosen deterministic Chebyshev stepsize schedule (2.6) with i.i.d. stepsizes from their limiting empirical measure. See [2, Chapter 3] for a discussion of this through the lens of exchangeability and De Finetti’s Theorem.

Concretely, if the inverse stepsizes $\beta_{t}:=1/\alpha_{t}$ are i.i.d. from some distribution $\mu$, then (forgoing measure-theoretic technicalities) the Law of Large Numbers ensures that the $n$-step convergence rate in (2.5) converges almost surely and in expectation to

	$\displaystyle\sup_{f\in\mathcal{Q},\,x_{0}\neq x^{*}}\left(\frac{\|x_{n}-x^{*}\|}{\|x_{0}-x^{*}\|}\right)^{1/n}$	$\displaystyle=\sup_{\lambda\in[m,M]}\prod_{t=0}^{n-1}|1-\alpha_{t}\lambda|^{1/n}$
		$\displaystyle=\exp\left(\sup_{\lambda\in[m,M]}\frac{1}{n}\sum_{t=0}^{n-1}\log|1-\beta_{t}^{-1}\lambda|\right)$
		$\displaystyle\overset{n\to\infty}{\longrightarrow}\exp\left(\sup_{\lambda\in[m,M]}\mathbb{E}_{\beta\sim\mu}\log|1-\beta^{-1}\lambda|\right)\,$		(2.8)

This informal derivation suggests that the distribution $\mu$ minimizing the asymptotic convergence rate (2.8) should be the appropriate limit as $n\to\infty$ of the empirical measure for the inverse Chebyshev stepsize schedule (this limit is the $\operatorname{{\mathrm{Arcsine}}(m,M)}$ distribution, since $a+b\cos(\Theta)$ has the Arcsine distribution on $[a-b,a+b]$ if $\Theta$ is uniformly distributed on $[0,\pi],$ see Figure 2), and moreover the optimal value should be the accelerated rate $\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$. This intuition is formalized by the following lemma.

This lemma was proven in the optimization literature via a long, tour-de-force integral calculation [32]. In Appendix A, we present an alternative, short proof from Altschuler’s thesis [2] that argues by connecting this extremal problem for stepsize distributions to a classical extremal problem in potential theory about electrostatic capacitors. The key feature of this connection is the following equalization property of the Arcsine distribution.

This equalization property uniquely defines the Arcsine distribution in that there is no other distribution satisfying (2.10); see Appendix A. In words, this property shows that the Arcsine distribution “equalizes” the convergence rate over all quadratic functions (parameterized by their curvatures $\lambda$). We emphasize that this is not just an algebraic identity: this equalization property also has a natural interpretation in terms of a game between an optimization algorithm (who chooses a stepsize to minimize the convergence rate) and an adversary (who chooses a quadratic function to maximize the convergence rate), see §6.2. In Appendix A we elaborate on the connections to potential theory, and leverage that connection to provide short proofs of these lemmas.

Clearly, if $f$ is not quadratic, then the gradient $\nabla f$ is not linear, and as a consequence the GD map is not linear. Assuming as before that $x^{*}=0$ without loss of generality, one can generalize the linear GD map (2.3) to the non-quadratic setting as

$\displaystyle x_{n}=\prod_{t=0}^{n-1}(I-\alpha_{t}H_{t})x_{0}\qquad\text{where}\qquad H_{t}:=\int_{0}^{1}\nabla^{2}f(ux_{t})du\,,$

(2.11)

since by the Fundamental Theorem of Calculus, $\nabla f(x_{t})=H_{t}x_{t}$ and thus $x_{t+1}=x_{t}-\alpha_{t}\nabla f(x_{t})=(I-\alpha_{t}H_{t})x_{t}$. The interpretation of $H_{t}$ is that it is the average Hessian between $x^{*}=0$ and $x_{t}$, i.e., it measures the relevant local curvature in the $t$-th iteration of GD. In the quadratic setting, $H_{t}=H$ for all $t$ since the Hessian $\nabla^{2}f\equiv H$ is constant, which recovers (2.3) from (2.11). However, for the general non-quadratic setting, the curvature $H_{t}$ varies in $t$. This is the overarching challenge mentioned above. For the purpose of accelerating GD, this challenge manifests in two interrelated ways: the lack of commutativity and predictability of the curvature $H_{t}$ along the trajectory of GD. Below we describe both challenges and how we circumvent them.

Here, we provide a short back-of-the-envelope calculation in order to demonstrate these tensions and the benefit of parallelization in the simple setting of univariate quadratic functions. For simplicity consider a single synchronization (i.e., $k=n$); the general case of multiple synchronizations ($k<n$) follows by replacing $n$ by $k$ in this calculation, and then repeating this convergence rate $n/k$ times. For shorthand, let $x_{t}^{(j)}$ denote the $t$-th iterate on the $j$-th machine, let $\alpha_{t}^{(j)}$ denote the corresponding stepsize, and let $f(x)=\frac{\lambda}{2}(x-x^{*})^{2}$ for $\lambda\in[m,M]$ (all univariate $\kappa$-conditioned quadratics are of this form).

Step 1: single machine. By the same analysis as in our main results (see §2.1 or §3.1) and then a Central Limit Theorem approximation, the convergence rate on the $j$-th machine is

$\displaystyle\left(\frac{\left\lvert x_{n}^{(j)}-x^{*}\right\rvert}{\left\lvert x_{0}^{(j)}-x^{*}\right\rvert}\right)^{1/n}=\exp\left(\frac{1}{n}\sum_{t=0}^{n-1}\log\left\lvert 1-\alpha_{t}^{(j)}\lambda\right\rvert\right)\approx\exp\left(\log\operatorname{\mathrm{R_{acc}}}+\frac{\sigma}{\sqrt{n}}\xi_{j}\right)=\operatorname{\mathrm{R_{acc}}}\cdot\exp\left(\frac{\sigma}{\sqrt{n}}\xi_{j}\right)\,,$

where $\xi_{j}\sim\mathcal{N}(0,1)$ and $\sigma=\operatorname{\mathrm{Var}}_{\alpha}(\log\left\lvert 1-\alpha\lambda\right\rvert)$. Note that the error of this approximation can be tightly controlled using standard techniques such as the Berry-Esseen CLT [57, Theorem 2.1.13].

Step 2: best machine. Taking the minimum over all $p$ processors gives the convergence rate of the parallelized scheme which takes the best final iterate:

$\displaystyle\min_{j\in[p]}\left(\frac{\left\lvert x_{n}^{(j)}-x^{*}\right\rvert}{\left\lvert x_{0}^{(j)}-x^{*}\right\rvert}\right)^{1/n}\approx\operatorname{\mathrm{R_{acc}}}\cdot\exp\left(\frac{\sigma}{\sqrt{n}}\min_{j\in[p]}\xi_{j}\right)\approx\operatorname{\mathrm{R_{acc}}}\cdot\exp\left(-\sigma\sqrt{\frac{2\log p}{n}}\right)\,.$

Above, the final step is due to the well-known fact that the minimum of $p$ i.i.d. standard Gaussians concentrates tightly around $-\sqrt{2\log p}$, see, e.g., [57, Chapter 2.5]. This factor $\exp(-\sigma\sqrt{(2\log p)/n})$ quantifies the speedup of parallelization since it is the size of the variation in Figure 3. By a Taylor approximation of $\log\operatorname{\mathrm{R_{acc}}}$, this parallelization provides asymptotic speedups when $\sigma\sqrt{(2\log p)/n}\gtrsim|\log\operatorname{\mathrm{R_{acc}}}|=|\log\tfrac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}|\approx\tfrac{1}{\sqrt{\kappa}}$, or equivalently when the number of processors $p\gtrsim\exp(\tfrac{n}{2\kappa\sigma^{2}})$.