Complexity Guarantees for Polyak Steps with Momentum

Introduction

We focus on unconstrained optimization problems of the form

where $f$ is strongly convex and has a Lipschitz continuous gradient with respect to the Euclidean norm. Very broadly speaking, the current numerical toolbox to solve these convex minimization problems contains two types of methods. On one hand, simple numerical schemes with explicit albeit conservative theoretical guarantees. These include gradient methods and their accelerated variants, and require knowing problem parameters, such as strong convexity parameters, or Hölderian error bounds (Bolte et al., 2007). On the other hand, adaptive methods, such as conjugate gradients or quasi-Newton, adapting much better to the objective function by estimating some of its regularity properties. For these methods, we typically have no theoretical justification for their improved performances or no computational complexity bounds at all.

Empirically, adaptive methods often perform significantly better than their parametric counterparts, and, by nature, require much less tuning. For example, roughly estimating regularity constants on-the-fly and plugging these estimates in parametric algorithms often produces fast algorithms with no theoretical guarantees. This phenomenon is illustrated in Figure 1 on logistic regression.

Although many advances have been made in designing optimization schemes adaptive to some types of parameters (e.g., Lipschitz constants, see discussions below), these results still leave a huge gap between theory and practice (as in Figure 1). In particular, estimating strong convexity coefficients while preserving convergence guarantees remains a challenging issue. Restart schemes are probably the most effective option among existing approaches for adapting to this type of parameters and do provide improved complexity estimates without any knowledge of strong convexity parameters, at the expense of a log scale grid search. However, while on paper the complexity of these schemes is nearly optimal, the presence of an outer loop clearly limits their practical effectiveness and their capacity to adapt to the function’s local regularity, which leaves a lot of margin for improvement, on the numerical front. Producing single loop algorithms adapting to local strong convexity (or Hölderian error bounds) and have nearly optimal complexity bounds is an important open problem which is the main focus of this work.

Here, we study the complexity of adaptive methods using Polyak steps, estimating the momentum term using information on the optimum objective value $f_{*}$ instead of the strong convexity constant. In some scenarios, such as “interpolation” in machine learning problems, the value of $f_{*}$ is known a priori (usually zero), and estimating it is much easier than estimating strong convexity, see e.g., (Asi and Duchi, 2019) for a recent discussion on these model assumptions.

The obvious next research question in this direction is to substitute knowledge on $f_{*}$ by weaker bounds. A first step in this direction is for example (Hazan and Kakade, 2019) which uses successive refinements of a lower bound on $f_{*}$. As it is, the proof in (Hazan and Kakade, 2019) contains several errors, but can be fixed. We hope, and believe, that such a mechanism could be used together with our momentum version of the Polyak steps.

For smooth optimization problems, simple line search strategies provide accelerated algorithms that adapt to the local gradient Lipschitz constant (Nesterov, 2013) and explicit adaptive complexity bounds can be derived for certain variants using the mean root Lipschitz constant (Scheinberg et al., 2014).

Restarts.

For smooth and strongly convex optimization problems (or more generally problems satisfying Hölderian error bounds), accelerated methods with optimal complexity bounds require knowledge of the strong convexity constant to compute iterates (Nesterov, 2013, 2018). In particular, Arjevani and Shamir (2016) show that this information is necessary when using oblivious steps. This quantity can be hard to estimate and a lot of effort has been put in the development of adaptive optimization methods preserving fast convergence rates (Lin and Xiao, 2014; Fercoq and Qu, 2016; Roulet and d’Aspremont, 2017). All these works are based on restart strategies (O’Donoghue and Candes, 2015; Nesterov, 2013) and although they exhibit fast theoretical convergence rates, they often contain parameters that have to be tuned in order to get good practical results, or require additional information on the function itself (e.g., its minimum $f^{*}$). Once again, while on paper the complexity of restart schemes is nearly optimal, the presence of an outer loop generally limits their capacity to adapt to the function’s local regularity and significantly affects empirical performance.

Quasi-Newton methods.

An important family of adaptive algorithms is composed with quasi-Newton methods. As the name suggests, these methods try to mimic the behavior of Newton schemes, by constructing an estimate of the hessian at the current point, using previous gradients. The most notable quasi-Newton method is certainly L-BFGS (Liu and Nocedal, 1989). These commonly used algorithms exhibit very fast empirical converge rates but only classical convergence rates comparable to that of gradient descent have been proven at this point (Byrd et al., 1987).

Conjugate gradient methods.

Conjugate gradient methods are probably among the most famous examples of adaptive algorithm. Firstly introduced for quadratic minimization (Hestenes and Stiefel, 1952), and motivated by nice theoretical guarantees (such as finite-time convergence), many variants have been introduced for going beyond quadratics (Fletcher and Reeves, 1964; Polyak, 1969; Fletcher, 1987)—see, for example, the nice survey (Hager and Zhang, 2006). Roughly speaking, at each iteration, the method constructs an update direction based on the gradient at the current iterate, and on the knowledge of the previous search directions. The next iterate is obtained by line-search in the update direction. Whereas conjugate gradient methods are widely used in practice (e.g Rodi and Mackie (2001); Volkwein (2004); Zhao et al. (2015)), and perform very well when they applies, there are barely any non-asymptotic convergence guarantees for those methods beyond unconstrained quadratic minimization.

Polyak step-sizes.

When the optimal value of the objective function value is known, a well-known adaptive strategy consists in using the so-called “Polyak step-sizes”—see e.g., (Polyak, 1987, Section 5.3.2) or (Nedic and Bertsekas, 2001; Boyd et al., 2003). The method consists in iterating gradient steps with step-sizes proportional to the primal gap at the current iterate. As opposed to most adaptive gradient methods mentioned above, this method comes with explicit theoretical properties, even beyond the quadratic optimization case.

Barzilai-Borwein step-sizes.

The Barzilai-Borwein (Barzilai and Borwein, 1988; Fletcher, 2005) method consists in gradient steps with adaptive step-sizes. It is another case with complete theory for quadratic optimization, but barely any performance guarantees in non-quadratic cases (it is even known to diverge on some problem instances).

Adaptive gradient steps

In (Malitsky and Mishchenko, 2019) the authors developed a step-size policy that adapts to the local geometry, together with nice theoretical guarantees.

2 Contributions

We develop and analyze an accelerated variant of the gradient method with Polyak steps that includes a momentum term and has better dependence on the condition number. We believe the Performance Estimation Program (PEP) technique used for obtaining the worst-case convergence guarantees is also of independent interest. As a byproduct, we also slightly improve convergence bounds for variants of the classical gradient method with Polyak steps (i.e. without momentum).

Classical Polyak Steps and Variants

Let us start with complexity bounds for gradient methods with Polyak steps for smooth and strongly convex optimization problems. Note that Polyak step sizes are usually discussed in the nondifferentiable setting—see (Polyak, 1987, Section 5.3.2) or (Nedic and Bertsekas, 2001; Boyd et al., 2003). We first recall the complexity of the gradient method with Polyak steps in the smooth strongly convex case, then derive similar bounds for two variants. For the first variant, we scale the steps by a factor two compared to standard Polyak steps, yielding a simple convergence proof with slightly improved theoretical guarantees. The second variant is a descent method, where the complexity bound is written in terms of the primal gap. We delay a full discussion of the proof mechanisms to Section 4, and the proofs themselves to the appendix.

The classical step size rule (Polyak) was mostly studied in the nonsmooth convex case (Polyak, 1987). For smooth strongly convex problems, it is known (see e.g., (Hazan and Kakade, 2019)) that

The two following propositions show that different step sizes policies (namely (Variant I) and (Variant II)) produce slightly improved convergence rates, matching the best known rates for gradient methods with known $\mu$ and $L$. The $\gamma_{k}$ are always well defined except when $x_{k}$ has a zero gradient, in this case we can simply stop the method since we have reached optimality. When it is well defined, $\gamma_{k}\in[\tfrac{1}{L},\tfrac{1}{\mu}]$ for (Variant I) and $\gamma_{k}\in[\tfrac{1}{L},\tfrac{2-\mu/L}{L}]$ for (Variant II). First, let us state that if we seek to decrease the distance to the optimal point, (Variant I) provides a rate that matches that of gradient descent with optimal (non-adaptive) step sizes (Nesterov, 2018).

where $\rho(\gamma)=\tfrac{(\gamma L-1)(1-\gamma\mu)}{\gamma(L+\mu)-1}$, and $\underset{\gamma\in[\tfrac{1}{L},\tfrac{1}{\mu}]}{\max}\;\rho(\gamma)=\tfrac{(L-\mu)^{2}}{(L+\mu)^{2}}$. Otherwise $\nabla f(x_{k})=0$ with $k\in[0,N]$.

If on the other hand we seek to decrease the primal gap, (Variant II) provides a rate that matches that of gradient descent with exact line search (de Klerk et al., 2017), at the expense of knowledge on L.

where $\rho(\gamma)=(L\gamma-1)\left(L\gamma(3-\gamma(L+\mu))-1\right)$, and $\underset{\gamma\in[\tfrac{1}{L},\tfrac{2L-\mu}{L^{2}}]}{\max}\;\rho(\gamma)=\tfrac{(L-\mu)^{2}}{(L+\mu)^{2}}$. Otherwise $\nabla f(x_{k})=0$ with $k\in[0,N]$.

In the following section, we study variants of those methods, where we aim to speed up convergence by incorporating a momentum term. Those methods follow in spirit the line of works on Nesterov’s acceleration (Nesterov, 2013), where we supersede knowledge of $\mu$ by that of $f_{*}$.

Acceleration with Polyak momentum

In the following, AGM refers to the Accelerated Gradient Method with momentum introduced by Nesterov (Nesterov, 1983, 2018). We are interested in optimizing a function $f\in\mathcal{F}_{\mu,L}$ without any information on the strong convexity constant $\mu$. However, as in the Polyak gradient method, we rely on the knowledge of $f^{*}$. We describe a single loop adaptive accelerated method (i.e. without restarts), with convergence rate of order $1-\left({\mu}/{L}\right)^{3/4}$, compared with $1-\mu/{L}$ for gradient descent, and $1-\left({\mu}/{L}\right)^{1/2}$ for its accelerated version with perfect knowledge of $\mu$.

Algorithm 2 is based on the AGM algorithm (Nesterov, 2018), in which the knowledge of $\mu$ is essential to set the constant momentum term $\beta_{k}=\beta_{*}=(\sqrt{L}-\sqrt{{\mu}})/(\sqrt{L}+\sqrt{{\mu}})$. Common convergence guarantees require a lower bound on the strong convexity. As a first step towards producing adaptive versions of AGM, Lemma 3.1 and Corollary 3.2 below guarantee that AGM with any momentum factor $\beta_{k}$ in $$ converges at least as fast as the classical gradient method.

where $V(x,y)=\frac{L-\mu}{2}\|x-y\|^{2}+f(y)-f_{*}$ and $\rho=1-\frac{\mu}{L}$.

We then get the following corollary on the primal gap.

Direct from Lemma 3.1 with $x_{0}=y_{0}$.

This result shows the robustness of AGM with respect to the momentum parameter. Adaptive strategies, that modify the momentum term in the algorithm automatically, thus at least enjoy the gradient method’s convergence rate when $\beta_{k}$ is kept within the interval $$—this is the case for both (Acc. Variant I) and (Acc. Variant II). To our knowledge, only non-blowup properties (Lin and Xiao, 2014, Lemma 1) were known when overestimating$\mu$.

The momentum term in (Acc. Variant I) was designed using the inverse of Polyak’s step as an estimate of the strong convexity parameter. The motivation for this choice of strong convexity estimate is the fact that under some mild assumptions on $f$ (i.e., for quadratic or self-concordant $f$), the quantity $\frac{\|\nabla f(z_{k})\|^{2}}{2(f(z_{k})-f_{*})}$ converges to the strong convexity constant at optimum when the $z_{k}$ are iterates of gradient descent algorithm with step-size 1/L.

Use Lemma 3.4 recursively and notice that $V(x_{0},y_{0})=f(x_{0})-f_{*}$.

where $V(x,y)=\frac{L}{2}\|\frac{1}{\sqrt{\rho}}(x-x_{*})-\sqrt{\rho}(y-x_{*})\|^{2}+f(y)-f_{*}$ and $\rho=\left(1+\left(\frac{\mu}{L}\right)^{\frac{3}{4}}\right)^{-1}$.

where $C=\left(\left(\tfrac{1}{\rho_{1}}-1\right)\left(1+\sqrt{\tfrac{L}{2\mu}}\right)^{2}+1\right)$, $\rho_{1}=\left(1+\left(\tfrac{\mu}{L}\right)^{\tfrac{3}{4}}\right)^{-1}$ and $\rho_{2}=\left(1+\sqrt{\tfrac{\mu}{L}}\right)^{-1}$. Otherwise $\nabla f(y_{k})=0$ with $k\in[0,N]$.

Proof mechanisms

Starting with the work of Drori and Teboulle (2014), computer-aided worst-case analyses of convex optimization methods have provided a generic technique producing convergence rates for many classical first-order algorithms. The results in (Drori and Teboulle, 2014; Taylor et al., 2017) use an interpolation argument to write the problem of finding the worst case behavior of an algorithm, given a convergence criterion, as a tractable semidefinite program—often referred to as a Performance Estimation Program (PEP). We adapted the technique for generating the complexity bounds on gradient methods with Polyak steps.

Our proofs were obtained by searching for Lyapunov (or potential) functions (see e.g. (Bansal and Gupta, 2019) for a recent survey). Due to space constraints, we do not detail how these potentials were obtained here, and refer the reader to the discussions on PEPs in (Taylor and Bach, 2019; Taylor et al., 2018) for more details. A related line of works (equivalent in many situations) is that of integral quadratic constraints (Lessard et al., 2016), which leverage results from control theory to perform worst-case complexity analysis. All these approaches were originally developed for non adaptive methods and in what follows, we show how we used the PEP approach for adaptive algorithms. A similar reasoning would allow adapting IQCs for adaptive methods as well.

To fix ideas and illustrate our procedure, we first analyze the worst case complexity of a variant of the classical gradient method with Polyak steps, and show improved convergence bounds compared to classical results (see Hazan and Kakade (2019) for a recent treatment). We consider the gradient method with Polyak steps described in Algorithm 1 with (Variant I) for $f\in\mathcal{F}_{\mu,L}$. Notice that there is a factor two in the step-size that is not present in the original Polyak step. This factor simplifies, and improves, the analysis for the convergence in terms of distance to the optimum.

To prove a linear convergence rate, we can focus on the improvement yielded by a single iteration of the form

We seek to bound the worst case (i.e., smallest) decrease in $\|x_{k+1}-x_{*}\|^{2}$ relative to $\|x_{k}-x_{*}\|^{2}$ when $x_{k+1}$ is obtained using the iteration in (7) for any function $f\in\mathcal{F}_{\mu,L}$ and any point $x_{k}$. In other words we seek to solve the following optimization problem

The key argument in (Drori and Teboulle, 2014; Taylor et al., 2017) is that the constraint on the regularity of the function $f$ in problem (8) can be replaced by a finite number of inequalities from Lemma 4.1. We get an upper bound on the optimum of problem (8) by relaxing the constraint $f\in\mathcal{F}_{\mu,L}$, keeping just two inequalities from Lemma 4.1 relating $x_{k}$ and $x_{*}$ to obtain the following relaxed problem

which is a semidefinite program. Given $\gamma$, $\rho(\gamma)$ can thus be computed efficiently and our relaxation upper bound on the convergence rate of the method is then given by the maximum value of $\rho(\gamma)$. Note that due to the definition of the step size, we only need to study $\rho(\gamma)$ on the interval $[\tfrac{1}{L},\tfrac{1}{\mu}]$. Figure 2 (left) plots $\rho(\gamma)$ for fixed values $\mu=0.1$ and $L=1$, and shows (right) the maximum value of $\rho(\gamma)$ for various condition numbers. In this experiment, the worst case convergence rates we obtained numerically appear to perfectly match the bound ${(L-\mu)^{2}}/{(L+\mu)^{2}}$.

These numerical observations can in fact be proven analytically as follows. Given a target convergence rate $\rho\in$, we need to show that

In practice, the numerical solution of the semidefinite program in (12) giving $\rho(\gamma)$ can be used to greedily narrow down the list of valid inequalities required by the proof.

Note that since (11) is a semialgebraic problem, we could have used sum-of-squares techniques to prove the convergence rate. However, the multipliers and the rates are fractions in $\gamma_{k}$. Since one usually doesn’t know in advance the form of the denominators, one needs relatively high degree polynomials in the SOS program. This means this approach suffers from the usual SOS issues of poor conditioning and scaling.

Numerical experiments

Numerical experiments with our algorithms are provided in Figure 3, respectively on least squares, regularized logistic regression and Lasso problems. For solving the Lasso problems, we used a proximal variant of Algorithm 2, whose details are provided in Appendix C.5. We respectively used the Sonar (Gorman and Sejnowski, 1988) and Musk (Dietterich et al., 1997) datasets.

In the experiments, when no analytical version of $f_{*}$ was available (for logistic regression and Lasso), we used ad hoc methods to obtain higher precision estimates of $f_{*}$. As previously discussed, a fundamental next step is to incorporate successive refinements of a lower bound on $f_{*}$ (a first step in this direction is for example (Hazan and Kakade, 2019)). One should notice that vanilla Polyak steps without momentum actually perform very well when they apply (see Appendix C.6 for a discussion on the performances of vanilla Polyak steps). We believe that modifying the accelerated Polyak so that it also adapts to the Lipschitz constant could make it more competitive, but the current state of the proofs does not allow it yet.

Conclusion and perspectives

We provided a momentum version of the Polyak steps, with an accelerated linear convergence rate. When $f_{*}$ is available, this method is easy to implement and requires no tuning at all. On the way, we illustrated the methodology that was used for obtaining those rates, for the special case of a gradient method with Polyak steps. This methodology relies on the recent developments on performance estimation problems (Drori and Teboulle, 2014; Taylor et al., 2017), which we adapted for studying our adaptive methods.

One of the main questions that remains open is to understand whether there exists a way to get the same convergence guarantees without using $f_{*}$. The robustness result of Lemma 3.1 is reassuring in the sense that a misspecified $f_{*}$ cannot break the algorithm (albeit worsening the convergence rate). We are confident that ideas introduced by Hazan and Kakade (2019) for Polyak steps could be used for our algorithm as well, and could potentially allow dealing with unknown $f_{*}$ at a reasonable cost. However it still appears as an unnatural trick that adds complexity to the method.

Let us mention that the problem of designing theoretically supported adaptive methods is an open question. We managed to design (Variant II), for which we used our methodology—to find a method that would use Polyak steps to make the primal gap decrease linearly at each iterations—, but designing adaptive accelerated methods appeared as much more daunting task.

Finally, we note that regular Polyak steps do not enjoy a known (working) proximal extension. On the contrary, our results suggest that its accelerated counterparts do work with proximal operators (for minimizing composite objective functions with a non-smooth term). Therefore, developing the theory in this direction is another natural next step.

The code used to obtain Figures 2-4-3 and to verify proofs is available at \urlhttps://github.com/mathbarre/PerformanceEstimationPolyakSteps.

The authors thank Konstantin Mishchenko and Yura Malitsky for insightful discussions on Polyak steps, and comments on a preliminary version of this work. The authors also thank three anonymous reviewers for their constructive feedbacks on the manuscript.

MB acknowledges support from an AMX fellowship. AT acknowledges support from the European Research Council (grant SEQUOIA 724063). AA is at CNRS & département d’informatique, École normale supérieure, UMR CNRS 8548, 45 rue d’Ulm 75005 Paris, France, INRIA and PSL Research University. AA acknowledges support from the French government under management of Agence Nationale de la Recherche as part of the ”Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), the ML & Optimisation joint research initiative with the fonds AXA pour la recherche and Kamet Ventures, as well as a Google focused award.

References

Appendix A Proof of Proposition 2.1

For proving the desired result, it is only necessary to consider a single iteration of Algorithm 1 with (Variant I). We use the following (in)equalities obtained from Lemma 4.1:

smoothness and strong convexity between $x_{k}$ and $x_{*}$, with multiplier $\lambda_{1}=\frac{2\gamma_{k}(\gamma_{k}L-1)}{\gamma_{k}(L+\mu)-1}$:

smoothness and strong convexity between $x_{*}$ and $x_{k}$, with multiplier $\lambda_{2}=\frac{2\gamma_{k}(1-\gamma_{k}\mu)}{\gamma_{k}(L+\mu)-1}$:

definition of the step-size policy, with multiplier $\lambda_{3}=\frac{\gamma_{k}(2-\gamma_{k}(L+\mu))}{\gamma_{k}(L+\mu)-1}$:

Given that $\lambda_{1},\lambda_{2}\geq 0$ (since $\tfrac{1}{L}\leq\gamma_{k}\leq\tfrac{1}{\mu}$), the following weighted sum is a valid inequality:

Using the fact that $x_{k+1}=x_{k}-\gamma_{k}\nabla f(x_{k})$, this weighted sum can be reformulated exactly as

(one can verify that both expressions are equal) with $\rho(\gamma)=\frac{(\gamma L-1)(1-\gamma\mu)}{\gamma(L+\mu)-1}$. Therefore, after $N$ iterations, we get

In addition, distance to optimality decreases, in the worst-case, with rate $\max_{\gamma}\rho(\gamma)$, with

because $\rho(\gamma)$ is a concave function of $\gamma$ on the interval $[\frac{1}{L},\frac{1}{\mu}]$, as $\rho^{\prime\prime}(\gamma)=-\frac{2L\mu}{(\gamma(L+\mu)-1)^{3}}\leq 0$, whose maximum is attained at $\gamma_{*}=\frac{2}{L+\mu}$. Note that substituting the expression of $\gamma_{k}$ inside the interpolation inequalities, instead of using it as an independent equality constraints, yields a considerably less tractable result.

Appendix B Proof of Proposition 2.2

Let us consider a single iteration of Algorithm 1, with step sizes (Variant II). The proof is a consequence of the following combination of inequalities obtained from Lemma 4.1:

smoothness and strong convexity between $x_{k}$ and $x_{*}$, with multiplier $\lambda_{1}=\gamma_{k}\mu(L\gamma_{k}-1)$:

smoothness and strong convexity between $x_{k+1}$ and $x_{*}$, with multiplier $\lambda_{2}=\gamma_{k}\mu$:

smoothness and strong convexity between $x_{k+1}$ and $x_{k}$, with multiplier $\lambda_{3}=1-\gamma_{k}\mu$:

definition of the step-size policy, with multiplier $\lambda_{4}=\frac{\gamma_{k}}{2}((L+\mu)\gamma_{k}-2)$:

Given that $\lambda_{1},\lambda_{2},\lambda_{3}\geq 0$ (due to $\tfrac{1}{L}\leq\gamma_{k}\leq\tfrac{2-\tfrac{\mu}{L}}{L}$), the following weighted sum is a valid inequality:

Using the expression $x_{k+1}=x_{k}-\gamma_{k}\nabla f(x_{k})$ (without substituting the expression of $\gamma_{k}$, whose value is encoded through the last equality of the list), this weighted sum can be rewritten exactly as

with $\rho(\gamma)=(L\gamma-1)\left(L\gamma(3-\gamma(L+\mu))-1\right)$ which, in turns, give

Finally, the worst-case convergence rate is $\max_{\gamma}\rho(\gamma)$ on the interval $[\tfrac{1}{L},\tfrac{2-{\mu}/{L}}{L}]$, for which

The proof follows from the following steps:

First, on the boundaries of the interval: (i) $\rho(\frac{1}{L})=0$ and (ii) $\rho(\frac{2-\frac{\mu}{L}}{L})=\frac{\left(L-\mu\right)^{4}}{L^{4}}\leq\frac{\left(L-\mu\right)^{2}}{\left(L+\mu\right)^{2}}$.

Secondly, in the interior of the interval: $\rho^{\prime}(\gamma)=L(3L\gamma-2)(2-(L+\mu)\gamma)$ is zero at $\gamma_{*}=\frac{2}{L+\mu}$ (inside the interval).

Therefore ${\rho(\gamma_{*})=\frac{\left(L-\mu\right)^{2}}{\left(L+\mu\right)^{2}}}$ and this is the maximum on the interval.

Appendix C Proof of § 3

In this section, we use $\rho=1-\mu/L$. The proof consists in combining the following inequalities obtained from Lemma 4.1:

smoothness and strong convexity between $x_{k}$ and $y_{k}$ with multiplier $\lambda_{1}=\rho$:

smoothness and strong convexity between $y_{k+1}$ and $x_{*}$ with multiplier $\lambda_{2}=1-\rho$:

smoothness and strong convexity between $y_{k+1}$ and $x_{k}$ with multiplier $\lambda_{3}=\rho$:

Given that $\lambda_{1},\lambda_{2},\lambda_{3}\geq 0$, the following weighted sum is a valid inequality

which can be reformulated exactly, using the notation

along with the expression of $\rho$, in the form

Therefore, using the assumption $\beta_{k}\in$, we finally arrive to the desired

C.2 Proof of Lemma 3.4

In this setting, we write $\rho(x)=\frac{1}{1+\frac{x}{L}}$. The proof consists in the following combination of inequalities obtained from Lemma 4.1:

The weighted sum is a valid inequality given that $\lambda_{1},\lambda_{2},\lambda_{3}\geq 0$:

which can be reformulated exactly, using the notation

along with the expression for $\rho(x)$, in the form

where the inequality follows from the sign of the term we removed, so it remains to show that

Indeed, evaluating the sign of the previous expression boils down to study that of $g(x)=4\sqrt{x}-\left(x-2x\sqrt{x}\right)-x^{2}$ on $$, which follows from:

C.3 Proof of Lemma 3.8

Our statement follows from a weighted sum of inequalities obtained from Lemma 4.1:

smoothness and strong convexity between $y_{k+1}$ and $x_{k}$, with multiplier $\lambda_{1}=1$:

smoothness and strong convexity between $x_{k}$ and $x_{*}$, with multiplier $\lambda_{2}=1-\rho$:

convexity between $x_{k}$ and $y_{k}$, with multiplier $\lambda_{3}=\rho$:

The weighted sum is a valid inequality given that $\lambda_{1},\lambda_{2},\lambda_{3}\geq 0$:

This inequality can be reformulated using the notations

where we used the facts that the following coefficients were nonnegative (proofs below) on the domain of interest:

$\tfrac{1}{2(L-\mu)}\geq 0$ (clear from the assumption $\mu\leq L$),

$\tfrac{1-\rho}{L}\geq 0$ (clear from $\rho\leq 1$),

$\tfrac{L\left(\rho^{3}-\beta^{2}\right)}{2\rho}\geq 0$ follows from $\left(\rho^{3}-\beta^{2}\right)\geq 0$, proved below,

$\tfrac{L^{2}(1-\rho)\left(\tfrac{\mu}{L}\rho(2\beta\rho-\beta(\beta+2)+\rho)+(\rho-1)(\beta-\rho)^{2}\right)}{2\left(\rho^{3}-\beta^{2}\right)(L-\mu)}\geq 0$ follows from previous points along with

The missing proofs are as follow. First, let us define $\kappa:=\tfrac{\mu}{L}\in$, the (inverse) condition number, and recall that we want to prove the expressions above to be nonnegative when $\rho=\tfrac{1}{1+\kappa^{3/4}}$ and $\beta_{-}\leq\beta\leq\beta_{+}$ with $\beta_{-}=\tfrac{\sqrt{1}-\sqrt{\kappa}}{\sqrt{1}+\sqrt{\kappa}}$ and $\beta_{+}=\tfrac{\sqrt{1}-\sqrt{\kappa}}{\sqrt{1}+\sqrt{\kappa}}$.

To show that $\rho^{3}-\beta^{2}\geq 0$, let us remark that the expression is a second order polynomial in the variable $\beta$ with negative curvature. Therefore, its minimum values are achieved on the boundary of the interval, and it is sufficient to show $\rho^{3}-\beta_{-}^{2}\geq 0$ and $\rho^{3}-\beta_{+}^{2}\geq 0$ for establishing our claim. For the case $\beta=\beta_{-}$, we get:

and we need to show that $\left(4-8\kappa^{1/4}+9\sqrt{\kappa}-4\kappa^{3/4}-4\kappa+9\kappa^{5/4}-8\kappa^{3/2}+4\kappa^{7/4}-\kappa^{2}\right)$ is non negative for all $\kappa\in$. For showing that, we perform the change of variable $x\leftarrow\kappa^{1/4}$ (which is invertible since $\kappa\in$), and study the polynomial

hence finally $\rho^{3}-\beta_{-}^{2}\geq 0$. For the case $\beta=\beta_{+}$, we obtain:

is nonnegative for all $\kappa\in$. After changing variable $x\leftarrow\kappa^{1/4}$ (which is invertible since $\kappa\in$), we study the polynomial

Similarly, the expression $p_{3}(\kappa)=\left(\kappa\rho(2\beta\rho-\beta(\beta+2)+\rho)+(\rho-1)(\beta-\rho)^{2}\right)$ is also a second order polynomial in $\beta$, with leading coefficient

Therefore, this quadratic function is also concave and we only need to verify the inequality on the boundary of the interval $[\beta_{-},\beta_{+}]$. In the case $\beta=\beta_{-}$, we get:

and we need to show that $\left(\kappa^{2}-\kappa^{7/4}+2\kappa^{3/2}+3\kappa^{5/4}-7\kappa+5\kappa^{3/4}+4\sqrt{\kappa}-7\sqrt{\kappa}+4\right)$ is nonnegative for $\kappa\in$. We change variables $x\leftarrow\kappa^{1/4}$ (which is invertible since $\kappa\in$), and study the polynomial

hence $p_{3}(\beta_{+})\geq 0$, which concludes the proof.

C.4 Proof of Proposition 3.9

The case $m=0$ results from Lemma 3.8 applied recursively and the case $m=\infty$ result from Proposition 3.5. In the following we consider that $m\in[1,N]$. Then for $(y_{k},x_{k})_{k\in[m+1,N]}$,

We can now apply Corollary 3.7. From the definition of $m$, we have

Therefore, by denoting $\rho_{2}=\left(1+\sqrt{\frac{\mu}{L}}\right)^{-1}$, we have the following inequalities

C.5 Proximal variants

A natural extension of smooth and strongly convex optimization is the case composite optimization

where $f\in\mathcal{F}_{\mu,L}$ and $h\in\mathcal{F}_{0,\infty}$ is a proper convex function with proximal operator available.

C.6 Study of standard Polyak steps

From numerical experiments, we noticed that (Variant I) was actually typically performing only slightly better than vanilla gradient descent. From a worst-case point of view, this is expected. However, our experiments (see Figure 1-3) suggest that regular Polyak steps (Polyak) actually perform much better than one could expect from its worst-case guarantees.

In this section, we provide a tentative explanation of this behavior, through experiments on a toy example. Figure 4 (top) was obtained by running the methods on a least squares problem (we used a rescaled version of the Sonar dataset, with regularity parameters $L=1$ and $\mu=0.01$).

Similar in spirit as in Figure 2 (left), we provide, in Figure 4, the worst-case ratio of $\|x_{k+1}-x_{*}\|^{2}/\|x_{k}-x_{*}\|^{2}$ (by solving (8) numerically for regular Polyak steps). One can observe that the worst case rate (using distances to optimum as the criterion) is slightly worse than that of (Variant I) (note that this rate can be improved through the use of refined Lyapunov functions).

In Figure 4, we provide the distributions of step size magnitudes observed through the optimization process on the toy example. One can notice that the distribution does not fully concentrate around the worst-case value (the value of $\gamma$ that achieves the worst-case) for (Polyak). A large proportion of effective step size values are even located in regions of fast convergence. On the contrary, for (Variant I), the distribution is much more concentrated around its worst-case. Those distributions strongly suggest that worst case analyses might not be the best way to explain the good practical behaviors of such adaptive methods.