Keywords:

Dans de nombreux domaines, comme par exemple l'optimisation, la performance d'une méthode est souvent caractérisée par son taux de convergence. Cependant, accélérer un algorithme requiert une certaine connaissance de la structure du problème et de telles améliorations sont le fruit d'une étude au cas-par-cas. De nombreuses techniques d'accélération ont été développées ces dernières décennies et sont maintenant massivement utilisées. En dépit de leur simplicité, ces méthodes sont souvent basées sur des arguments purement algébriques et n'ont généralement pas d'explications intuitives.

Récemment, de nombreux travaux ont été menés pour faire des liens entre les algorithmes accélérés et d'autres domaines scientifiques, comme par exemple avec la théorie du contrôle ou des équations différentielles. Cependant, ces explications reposent souvent sur des arguments complexes et la plupart utilise des outils non-conventionnels dans leur analyse.

Une des contributions de cette thèse est une tentative d'explication des algorithmes accélérés en utilisant la théorie des méthodes d'intégration, qui a été très étudiée et jouit d'une analyse théorique solide. En particulier, nous montrons que les méthodes d'optimisation sont en réalité des instances de méthode d'intégration, lorsqu'on intègre l'équation du flot de gradient. Avec des arguments standards, nous expliquons intuitivement l'origine de l'accélération.

De l'autre côté, les méthodes accélérées ont besoin de paramètres supplémentaires, en comparaison d'autres méthodes plus lentes, qui sont généralement difficiles à estimer. De plus, ces schémas sont construits pour une configuration particulière et ne peuvent pas être utilisés autre part.

Ici, nous explorons une autre approche pour accélérer les algorithmes d'optimisation, qui utilise des arguments d'accélération générique. En analyse numérique, ces outils ont été développés pour accélérer des séquences de scalaires ou de vecteurs, en construisant parallèlement une autre séquence avec un meilleur taux de convergence. Ces méthodes peuvent être combinées avec un algorithme itératif, l'accélérant dans la plupart des cas. En pratique, ces méthodes d'extrapolation ne sont pas tellement utilisées, notamment dû à leur manque de garanties de convergence et leur instabilité. Nous étendons ces méthodes en les régularisant, ce qui permettra une analyse théorique plus profonde et des résultats de convergence plus fort, en particulier lorsqu'elles sont appliquées à des méthodes d'optimisation. v vi

Contributions

This thesis is divided into two main parts. Chapter 2 discusses the links between integration theory in numerical analysis and optimization methods. Chapter 3 develops a generic acceleration method for iterative algorithm, and Chapter 4 applies the results of Chapter 3 to optimization algorithms.

Chapter 1: This chapter is an introduction to optimization, integration theory and sequence acceleration in numerical analysis, which are the main topics of this thesis.

Chapter 2: This chapter introduces basics of integration theory and present standard integration methods. Then, it makes the link between these methods and optimization algorithms.

These links allow an unified analysis of many optimization schemes and give new intuitions about their theoretical properties, such as improved rate of convergence. This chapter is based on the results of [START_REF] Scieur | Integration methods and optimization algorithms[END_REF].

Chapter 3: This chapter introduces the notion of acceleration in numerical analysis, then proposes a novel algorithm for accelerating vector sequences. In comparison with previous work in this vein, this new scheme is more stable and has a proven rate of convergence. The analysis is then pushed to more complex vector sequences, thus making the acceleration method more versatile. The results of this chapter have been partially published in [START_REF] Scieur | Regularized nonlinear acceleration[END_REF][START_REF] Scieur | Nonlinear acceleration of stochastic algorithms[END_REF][START_REF] Scieur | Nonlinear acceleration of deep neural networks[END_REF], and the other have been revised.

Chapter 4: This chapter applies the results of Chapter 3 to standard optimization methods, in the deterministic or stochastic setting. It also analyses the effects of acceleration to the optimization of neural networks. Then, we present some numerical experiments on different data sets using various, standard optimization algorithms to show the effectiveness of the acceleration method. The results of this chapter have been published in [START_REF] Scieur | Regularized nonlinear acceleration[END_REF][START_REF] Scieur | Nonlinear acceleration of stochastic algorithms[END_REF][START_REF] Scieur | Nonlinear acceleration of deep neural networks[END_REF][START_REF] Scieur | Nonlinear acceleration of CNNs[END_REF].

Chapter 1

Introduction

In this thesis, we develop tools to analyse and improve optimization methods, especially for convex functions. We will thus recall some key results from optimization, but also from acceleration methods and integration schemes.

Convex Optimization

One of the main goals of optimization is the design of efficient algorithms for solving

min x∈C f (x)
in the variable x ∈ C ⊂ R d , where C is the feasible set and f is the objective function. In general, this is an impossible problem, unless we add some assumptions.

A common structural hypothesis is convexity. This makes the minimization problem tractable by generic algorithms. If an optimization problem is convex, this means that • The set C is convex, i.e., ∀x, y ∈ C, we also have αx + (1α)y ∈ C, α ∈ [0, 1].

• The function f is convex, i.e., its epigraph {(x, z) : x ∈ C, z ≥ f (x)} is a convex set. This is a sufficient condition to ensure that there exist methods, potentially slow, to solve this minimization problem. With additional structure it is possible to design more efficient algorithms with better convergence properties. For simplicity, in the rest of this thesis, we consider the case of unconstrained optimization, which means C = R d .

Smoothness and Strong Convexity

Classical assumptions in optimization are smoothness (requires that f is differentiable) and (strong) convexity. We assume thorough the thesis that ∇f exists.

• Smoothness. The function f is smooth with constant L if, for all x, y ∈ R d , f (y) ≤ f (x) + ∇f (x) T (yx) + L 2 yx 2 .

• Strong convexity. The function f is strongly convex with constant µ if, for all x, y ∈ R d ,

f (y) ≥ f (x) + ∇f (x) T (y -x) + µ 2 y -x 2 .
These definitions are valid for any norm, but we restrict here to the standard case where • stands for the Euclidean norm. These assumptions upper and lower bound the function by quadratics. With smoothness and strong convexity, it is possible to design efficient algorithms, with good theoretical guarantees.

Deterministic Optimization Algorithms

One of the oldest scheme in optimization is the so-called steepest descent, or gradient descent.

It consists in computing the gradient at one point, then move in its opposite direction. More formally,

x i+1 = x i -h∇f (x i ) (Gradient Step)
where h is the step size parameter. The convergence of this scheme is dependent of h and f : it is possible to show that gradient descent converges with a certain rate, function of the regularity of f , if h is well-chosen.

Proposition 1.1.1. (Nesterov [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]) Let the sequence {x i } i=1...N be generated by (4.3). Let f * be the minimum of f . If f is smooth and convex, and if 0 < h < 2 L then

f (x i ) -f * ≤ 2L(f (x 0 ) -f * ) x 0 -x * 2 2L x 0 -x * 2 + i • (f (x 0 ) -f * ) .
If, in addition, f is strongly convex and 0 < h < 2 µ+L , then the minimizer x * : f (x * ) = f * is unique and

x i -x * 2 ≤ 1 - 2hµL µ + L i x 0 -x * 2 .
Despite its simplicity, the rate of convergence of gradient descent is not so bad. However, it is possible to use an improved version which uses past iterates, originally found in [START_REF] Nemirovskii | Optimal methods of smooth convex minimization[END_REF]. Its latest version is written

   x i+1 = y i -1 L ∇f (y i ) y i+1 = (1 + β i )x i -β i x i-1 (Fast Gradient)
where β i = i-1 i-2 when the function is convex, or

β i = √ L- √ µ √ L+
√ µ when we also have strong convexity. In comparison with the gradient method, this method roughly square-roots the number of required iterations to reach -accuracy, f (x)f * ≤ .

More formally, the next proposition gives the rate of convergence of the fast gradient method. Proposition 1.1.2. (Nesterov [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]) Let the sequence {x i } i=1...N be generated by (Fast Gradient).

Let f * be the minimum of f . If f is smooth and convex, and if

β i = i-1 i-2 , then f (x i ) -f * ≤ constant • 1 1 + i 2 f (x 0 ) -f * + L 2 x 0 -x * 2 .
If, in addition, f is strongly convex and

β i = √ L- √ µ √ L+ √ µ , then f (x i ) -f * ≤ constant • 1 - µ L i f (x 0 ) -f * + L 2 x 0 -x * 2 .
In the second case, the choice of β is crucial but implies the knowledge of constants µ and L. Estimating the smoothness constant is actually easy. It suffices to check if the estimate of L, says L, satisfies the descent condition

f x - 1 L ∇f (x) ≤ f (x) - 1 2 L ∇f (x) 2
when we perform a gradient step. If not, we increase L until the condition is met. However, the strong convexity parameter is much harder in practice, and its estimation remains a challenge [START_REF] Fercoq | Restarting accelerated gradient methods with a rough strong convexity estimate[END_REF].

The algorithm (Fast Gradient) remains mysterious from an intuitive point of view. Of course, if we use more memory then we can only have a faster algorithm, but it remains unclear why it makes the method faster than gradient method. The proof relies on purely algebraic arguments and the key mechanism behind acceleration remains elusive, which led to various interpretations of it [START_REF] Bubeck | A geometric alternative to Nesterov's accelerated gradient descent[END_REF][START_REF] Zhu | Linear coupling: An ultimate unification of gradient and mirror descent[END_REF][START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF][START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF][START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF][START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF][START_REF] Wilson | A Lyapunov analysis of momentum methods in optimization[END_REF].

Stochastic Optimization

In the case of stochastic optimization, we do not have an explicit access to ∇f (x). Instead, we can only compute the noisy version

∇f (x) = ∇f (x) + ε.
This is a typical setting in machine learning, where the objective function f is a sum of error function f i , i.e.,

f (x) = m i=1 f i (x) + R(x)
where R(x) is a regularizer, for example, R(x) = x 2 . Instead of the full gradient, we only compute the gradient for function f i , where i is sampled uniformly in {1, . . . , m}.

Stochastic optimization is typically challenging as classical algorithms are not convergent (for example, gradient descent or Nesterov's accelerated gradient). Even the averaged version of stochastic gradient descent with constant step size does not converge.

Algorithms such as SAG [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF], SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], SDCA [START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss minimization[END_REF] and SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] accelerate convergence using a variance reduction technique akin to control variate in Monte-Carlo methods. Their rate of convergence depends on 1µ/L and thus does not exhibit an accelerated rate on par with the deterministic setting (in 1µ/L). Recently a generic acceleration algorithm called Catalyst [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], based on the proximal point methods improved this rate of convergence. On the other hand, recent papers, for example [START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF] (Accelerated SDCA) and [START_REF] Allen-Zhu | Katyusha: The first direct acceleration of stochastic gradient methods[END_REF] (Katyusha), propose algorithms with accelerated convergence rates, if the strong convexity parameter is given.

Links Between Optimization and Numerical Integration

We have seen that Nesterov's method is an efficient algorithm for minimizing a smooth, convex function f . However, its principal arguments are mostly algebraic, so it is hard to explain it with simple, intuitive arguments. On the other hand, there exist many different interpretations of gradient descent. One of the most popular way to explain it is the continuous time interpretation.

Gradient Flow and Euler's method

Consider the gradient descent algorithm,

x i+1 = x i -h∇f (x i ).
We can write it so that we isolate the x's from the gradient,

x i+1 -x i h = -∇f (x i ).
If we take the limit h → 0, we end up with a differential equation, called the gradient flow,

x(0) = x 0 dx dt ẋ = -∇f (x(t)).

(Gradient Flow)

From a geometrical point of view, the gradient flow equation tells us that the trajectory of x(t)

follows always the steepest descent, similarly to the trajectory of a water drop in a valley, and this corresponds quite accurately to the discrete trajectory of x i .

In fact, the discrete gradient descent corresponds exactly to the integration of (Gradient Flow) when using the well-known Euler's method. When we face an ordinary differential equation of the type dx dt = g(x(t)),

Euler's method estimates dx dt using the forward difference

dx dt ≈ x i+1 -x i t i+1 -t i , x(t) ≈ x(t i ).
If we assume t i+1t i = h constant, we end up with

x i+1 = x i + hg(x i ).
Here, we see a perfect match between Euler's method and gradient scheme when g = -∇f .

However, when moving to more complex algorithm such as Nesterov's fast gradient, things get more complicated.

Interpretation of Fast Gradient Algorithm

The main issue with the explanation of fast gradient descent using continuous time interpretation is the following: we can either refine the differential equation, or the integration method.

For example, on the basis of the work done by Su et al. [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF], [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF] proposed to analyse a complex differential equation, based on the Bregman Lagrangian function, written (in the Euclidean case)

L(x, v, t) = e αt+γt e -αt v 2 2e -βt f (x) , where α t , β t and γ t follows some "ideal scaling" conditions. Then, the differential equation is written d dt x(t) + e -αt dx(t) dt = -e αt+βt ∇f (x(t)).

In their work, the authors proved an accelerated rate in continuous time for this differential equation. Then they propose a smart, but unusual, discretization based on Euler's scheme which preserves the rate of convergence in the discrete case.

In our work [START_REF] Scieur | Integration methods and optimization algorithms[END_REF], we propose to keep the simple and intuitive differential equation (Gradient Flow). Instead, we will use a different integration rule to approximate the curve x(t). There exists two big families of integration methods, named Runge-Kutta and Linear Multistep methods. They all are well-studied, and used in many applications. For example in the Matlab software the function ode45 is a specific instance of Runge-Kutta method.

In this thesis, we will introduce the theory of integration methods and make the parallel with existing optimization scheme, and show that acceleration can be explained by the usage of larger step size. We focus on the family of (explicit) Linear Multistep methods for integrating (1.1), written s i=1

ρ i x i = h s-1 i=1 σ i g(x i )
for some coefficients ρ, σ and step-size h (later, ρ and σ will represent polynomials). The previous method is explicit because we can compute x s using the past iterates only, since

x s = 1 ρ s   s-1 i=1 ρ i x i -h s-1 i=1 σ i g(x i )   .
It happens there is a perfect match between many optimization algorithms and this family of methods, when g is actually -∇f . For example, the Nesterov's method can be written

x i+2 = 1 β = 1 ρs       (1 + β)x i+1 -x i+2 = s-1 i=1 ρixi + 1 L -β -∇f (x i+1 + (1 + β) -∇f (x i+1 = 1 h s-1 i=1 σig(xi)       where β = √ L- √ µ √ L+ √ µ .
We clearly see that Nesterov's method can also be an instance of a linear multi-step method. In addition, if we assume that the "Nesterov integration scheme" is stable and consistent (two important conditions that almost all integration methods should meet, they basically ask the method to integrate perfectly first order polynomials), we can retrieve its step size h. Surprisingly, for smooth and strongly convex problem, the step size is L µ larger than the step size of gradient descent, which means the method integrates the differential equation L µ times faster. Again, this corresponds to the theoretical speedup of Nesterov's method over the gradient scheme.

Acceleration Methods and Optimization Algorithms

We have seen that accelerated methods in optimization exist, and present an improved rate of convergence. However, their main weakness is their lack of adaptivity to the problem constants.

For example, in the case of Nesterov's method, we have seen it is possible to estimate the Lipchitz constant while the algorithm is running, but not the strong convexity parameter. If we have acess to a good estimate of the strong convexity constant, this is not a problem. However, if it is not the case, we have to choose between

• the strong-convex variant of Nesterov's scheme (β i = √ L- √ µ √ L+ √ µ
) with a bad estimate of µ, which leads to a bad convergence rate or

• the convex variant of Nesterov's scheme (β i = i-1 i-2 ), which does not have a linear rate of convergence.

On the other hand, there exist generic acceleration methods such as Catalyst [START_REF] Lin | A universal catalyst for first-order optimization[END_REF] or LBFGS [START_REF] Boyd | Convex optimization[END_REF]. However, the numerical performance of Catalyst is not as good as the Fast gradient method, and LBFGS lacks theoretical guarantees.

Anderson Acceleration

In this thesis, we propose a generic acceleration method based on Anderson Acceleration [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]. This method is capable to improve the rate of an iterative algorithm of the form

x i+1 = g(x i )
where the Jacobian of g is symmetric. For example, this corresponds to the gradient descent where g is

g(x) = x -h∇f (x).
The main idea of Anderson acceleration is to combine the previous iterates x i with coefficient c i to form an extrapolation x extr . This extrapolation aims to approximate the solution x * . In other words, Anderson acceleration computes c i such that

x extr = N i=0 c i x i ≈ x * .
As opposed to Nesterov's method, where the coefficient β i are known in advance, the coefficients c i are computed using only the iterates x i . However, the result has major stability issues (a small noise in the iterates x i produces large perturbation on the extrapolation x extr ), and thus Anderson acceleration is not suitable in most cases.

Regularized Nonlinear Acceleration

In Chapter 3, we introduce the Regularized Nonlinear Acceleration technique, which is basically a regularized version of Anderson acceleration. This makes the algorithm more stable, with proven convergence bound (Theorem 3.3.9). In addition regularization makes the acceleration scheme robust to any kind of noise, so we can derive a similar bound for stochastic algorithms (Proposition 3.3.11). Finally, we identify a larger class of algorithm, written

x i+1 = g(y i ), y i+1 = i+1 j=1 α (i) j x j + i j=0 β (i) j y j .
where the coefficients α and β are arbitrary. This class contains many different optimization schemes, such as gradient method but more importantly Nesterov's algorithm. We then propose a novel method which "accelerates" this family of methods, while keeping their theoretical properties (Proposition 3.5.3 and Algorithm 5).

In Chapter 4, we apply the results of Chapter 3 to optimization methods for deterministic and stochastic algorithms, and also on neural networks. In particular, we show this acceleration technique improves the rate of convergence of many algorithms, and usually outperforms optimal methods such as the fast gradient in the deterministic case, or Katyusha [START_REF] Allen-Zhu | Katyusha: The first direct acceleration of stochastic gradient methods[END_REF] for the stochastic setting.

Introduction

The gradient descent algorithm used to minimize a function f has a well-known simple numerical interpretation as the integration of the gradient flow equation, written

x(0) = x 0 ẋ(t) = -∇f (x(t)),
(Gradient Flow) using Euler's method. This appears to be a somewhat unique connection between optimization and numerical methods, since these two fields have inherently different goals. On one hand, numerical methods aim to get a precise discrete approximation of the solution x(t) on a finite time interval. More sophisticated methods than Euler's were developed to get better consistency with the continuous time solution but still focus on a finite time horizon [see for example 75].

On the other hand, optimization algorithms seek to find the minimizer of a function, which corresponds to the infinite time horizon of the gradient flow equation. Structural assumptions on f led to more sophisticated algorithms than the gradient method, such as the mirror gradient method [see for example 8, 6], proximal gradient method [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] or a combination thereof [START_REF] Duchi | Composite objective mirror descent[END_REF][START_REF] Nesterov | Universal gradient methods for convex optimization problems[END_REF].

Among them Nesterov's accelerated gradient algorithm [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o (1/k2)[END_REF] is proven to be optimal on the class of smooth convex or strongly convex functions. This last method was designed with the lower complexity bounds in mind, but the proof relies on purely algebraic arguments and the key mechanism behind acceleration remains elusive, which led to various interpretations of it [START_REF] Bubeck | A geometric alternative to Nesterov's accelerated gradient descent[END_REF][START_REF] Zhu | Linear coupling: An ultimate unification of gradient and mirror descent[END_REF][START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF].

A recent stream of papers used differential equations to model the acceleration behavior and offer a better interpretation of Nesterov's algorithm [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF][START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF][START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF][START_REF] Wilson | A Lyapunov analysis of momentum methods in optimization[END_REF]. However, the differential equation is often quite complex, being reverse-engineered from Nesterov's method itself, thus losing the intuition. Moreover, integration methods for these differential equations are often ignored or are not derived from standard numerical integration schemes, because the convergence proof of the algorithm does not require the continuous time interpretation.

Here, we take another approach. Rather than using a complicated differential equation, we use advanced multi-step methods to discretize the basic gradient flow equation in (Gradient Flow).

These lesser known methods, developed decades ago by numerical analysts, directly correspond to various well-known optimization algorithms. In particular, Nesterov's method can be seen as a stable (a small noise in the initial conditions produces bounded perturbations) and consistent (the limit of the numerical approximation when h → 0 is the differential equation) gradient flow discretization scheme that allows bigger step sizes in integration, leading to faster convergence.

The chapter is organized as follows. In Section 2.2 we present our setting and recall classical results on differential equations. We then review definitions and theoretical properties of integration methods called linear multi-step methods in Section 2.3. Linear one-step and twostep methods are detailed in Section 2.4 and linked to optimization algorithms in Section 2.5 (strongly convex case) and Section 2.6 (convex case). Finally, we propose several extensions for (Gradient Flow), which can be used to explain proximal methods (Section 2.7), non-euclidean method (Section 2.8) or a combination thereof (Section 2.9).

Gradient flow

We seek to minimize a L-smooth and µ-strongly convex function f defined on R d and look in that purpose at the discretization the gradient flow equation (Gradient Flow), given by the following ordinary differential equation (ODE)

d dt x(t) ẋ(t) = g(x(t)) x(0) = x 0 , (ODE)
where g comes from a potential -f , meaning g = -∇f . Smoothness of f means Lipschitz continuity of g, i.e.

g(x) -g(y) ≤ L x -y , for every x, y ∈ R d ,
where . is the Euclidean norm. This property ensures the existence and uniqueness of the solution of (ODE) (see [START_REF] Süli | An introduction to numerical analysis[END_REF]Theorem 12.1]). Strong convexity of f means strong monotonicity of -g, i.e., µ xy 2 ≤xy, g(x)g(y) , for every x, y ∈ R d , and ensures that (ODE) has a unique point x * such that g(x * ) = 0, called the equilibrium. This is the minimizer of f and the limit point of the solution, i.e., x(∞) = x * . In the numerical analysis literature, strong monotonicity of -g is referred to as one-sided Lipschitz-continuity of g. Whatever its name, this property essentially contracts solutions of (ODE) with different initial points as showed in the following proposition, where µ controls the rate of contraction. Proposition 2.2.1. Assume that x 1 (t), x 2 (t) are solutions of ODE with different initial conditions, where -g is strongly monotone. Then, x 1 (t) and x 2 (t) get closer as t increases, precisely,

x 1 (t) -x 2 (t) 2 ≤ e -2µt x 1 (0) -x 2 (0) 2 .
Proof. Let L(t) = x 1 (t)x 2 (t) 2 If we derive L(t) over time,

d dt L(t) = 2 x 1 (t) -x 2 (t), ẋ1 (t) -ẋ2 (t) = 2 x 1 (t) -x 2 (t), g(x 1 (t)) -g(x 2 (t))
≤ -2µL(t).

We thus have L(t) ≤ e -2µt L(0), which is the desired result.

Strong convexity allows us to control the convergence rate of the potential f and the solution x(t) as recalled in the following proposition.

Proposition 2.2.2. Let f be a L-smooth and µ-strongly convex function and x 0 ∈ dom(f ).

Writing x * the minimizer of f , the solution x(t) of (Gradient Flow) satisfies

f (x(t)) -f (x * ) ≤ (f (x 0 ) -f (x * ))e -2µt
(2.1)

x(t) -x * ≤ x 0 -x * e -µt . (2.2)
Proof. Indeed, if we derive the left-hand-side of (2.1),

d dt [f (x(t)) -f (x * )] = ∇f (x(t)), ẋ(t) = -∇f (x(t)) 2 .
Using that f is strongly convex, we have (see [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF])

f (x) -f (x * ) ≤ 1 2µ ∇f (x) 2 ,
and therefore

d dt [f (x(t)) -f (x * )] ≤ -2µ[f (x(t)) -f (x * )].
Solving this differential equation leads to the desired result. We can apply a similar technique for the proof of (2.2), using that

µ x -y 2 ≤ ∇f (x) -∇f (y), x -y ,
for strongly convex functions (see again [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]).

We now focus on numerical methods to integrate (ODE).

Numerical integration of differential equations

In general, we do not have access to an explicit solution x(t) of (ODE). We thus use integration algorithms to approximate the curve (t, x(t)) by a grid

(t k , x k ) ≈ (t k , x(t k )) on a finite interval [0, t max ].
For simplicity here, we assume the step size h k = t kt k-1 constant, i.e., h k = h and t k = kh. The goal is then to minimize the approximation error

x k -x(t k ) for k ∈ [0, t max /h].
We first introduce Euler's explicit method to illustrate this on a basic example.

Euler's explicit method

Euler's explicit method is one of the oldest and simplest schemes for integrating the curve x(t).

The idea stems from the Taylor expansion of x(t) which reads

x(t + h) = x(t) + h ẋ(t) + O(h 2 ).
When t = kh, Euler's explicit method approximates x(t + h) by x k+1 by neglecting the second order term,

x k+1 = x k + hg(x k ).
In optimization terms, we recognize the gradient descent algorithm used to minimize f . Approximation errors in an integration method accumulate with iterations, and as Euler's method uses only the last point to compute the next one, it has only limited control over the accumulated error.

Linear multi-step methods

Multi-step methods use a combination of several past iterates to improve convergence. Throughout the paper, we focus on linear s-step methods whose recurrence can be written

x k+s = - s-1 i=0 ρ i x k+i + h s i=0 σ i g(x k+i ), for k ≥ 0, (2.3)
where ρ i , σ i ∈ R are coefficients of polynomials ρ, σ and h is again the step size. Each new point

x k+s is a function of the information given by the s previous points. If σ s = 0, each new point is given explicitly by the s previous points and the method is then called explicit. Otherwise each new point requires solving an implicit equation and the method is then called implicit.

To simplify notation we use the shift operator E, which maps Ex k → x k+1 . Moreover, if we write g k = g(x k ), then the shift operator also maps Eg k → g k+1 . Recall that a univariate polynomial is called monic if its leading coefficient is equal to 1. We now give the following concise definition of s-step linear methods.

Definition 2.3.1. A linear s-step method integrates an (ODE) defined by g, x 0 , using a step size h and x 1 , . . . , x s-1 starting points by generating a sequence x k that satisfies

ρ(E)x k = hσ(E)g k , for every k ≥ 0,
where ρ is a monic polynomial (i.e. its leading coefficient is equal to 1) of degree s with coefficients ρ i , and σ a polynomial of degree s with coefficients σ i .

A linear s-step method is uniquely defined by the polynomials (ρ, σ). The sequence generated by the method then depends on the starting points and the step size. Each linear multi step method has a twin sister, called one-leg method, which generates the sequence xk following

ρ(E)x k = σ(1)hg σ(E) σ(1) xk .
It is possible to show a bijection between x k and xk [START_REF] Dahlquist | On one-leg multistep methods[END_REF]. We quickly mention one-leg methods in Section 2.7 but we will not go into more details here. We now recall a few results describing the performance of multi-step methods.

Stability

Stability is a key concept for integration methods. First of all, consider two curves x(t) and y(t), both solutions of an (ODE) defined by g, but starting from different points x(0) and y(0). If the function g is Lipchitz-continuous (but not especially monotone), we have show in Proposition 2.2.1 that the distance between x(t) and y(t) is bounded on a finite interval, i.e.

x(t) -y(t) ≤ C x(0) -y(0) ∀t ∈ [0, t max ],
where C may depend exponentially on t max . We would like to have a similar behaviour for our sequences x k and y k , approximating x(t k ) and y(t k ), i.e.

x k -y k ≈ x(t k ) -y(t k ) ≤ C x(0) -y(0) ∀k ∈ [0, t max /h], (2.4)
when h → 0, so k → ∞. Two issues quickly arise, namely:

• For a linear s-step method, we need s starting values x 0 , ..., x s-1 . Condition (2.4) will therefore depend on all these starting values and not only x 0 .

• Any discretization scheme introduces at each step an approximation error, called local error, which is accumulated over time. We denote this error by loc (x k+s ) and define it as

loc (x k+s ) x k+s -x(t k+s ) if x k+s is computed using x k = x(t k ), ..., x k+s-1 = x(t k+s-1 )
, where x(t) is the solution of the differential equation.

In other words, the difference between x k and y k can be described as follows

x ky k ≤ Error in the initial condition + Accumulation of local errors.

We now write a complete definition of stability, inspired from Definition 6.3.1 from Gautschi [START_REF] Gautschi | Numerical analysis[END_REF].

Definition 2.3.2. A linear multi-step method is stable if, for two sequences x k , y k generated by (ρ, σ) using any sufficiently small step size h > 0, from the starting values x 0 , ..., x s-1 , and y 0 , ..., y s-1 , we have

x k -y k ≤ C   max i∈{0,...,s-1} x i -y i + tmax/h i=1 loc (x i+s ) -loc (y i+s )   , (2.5
)

for any k ∈ [0, t max /h].
Here, the constant C may depend on t max but is independent of h.

When h tends to zero, we may recover equation (2.4) only if the accumulated local error tends also to zero. We thus need

lim h→0 1 h loc (x i+s ) -loc (y i+s ) = 0 ∀i ∈ [0, t max /h].
This condition is called consistency. Once this condition satisfied, we still need to ensure

x k -y k ≤ C max i∈{0,...,s-1} x i -y i , (2.6)
and this condition is called zero-stability (because, surprisingly, it suffices to check this condition for the special case g(x) = 0).

Truncation error and consistency

The truncation error of a linear multi-step method is a measure of the local error loc (x k ) made by the method, normalized by h. More precisely it is defined using the difference between the step performed by the algorithm and the step which reaches exactly x(t k+s ), with

T (h) x(t k+s ) -x k+s h assuming x k+i = x(t k+i ), i = 0, . . . , s -1. (2.7)
This definition does not depend on k but on the recurrence of the linear s-step method and on the (ODE) defined by g and x 0 . We can use this truncation error to define consistency.

Definition 2.3.3. An integration method for an (ODE) defined by g, x 0 is consistent if and only if, for any initial condition x 0 ,

lim h→0 T (h) = 0.
The following proposition shows there exist simple conditions to check consistency, which rely on comparing a Taylor expansion of the solution with the coefficients of the method. (2.8)

Proof. Assume t k = kh. If we expand g(x(t k )) we have,

g(x(t k )) = g(x 0 ) + O(h).
If we do the same thing with x(t k ), we have given (2.3),

x(t k ) = x 0 + kh ẋ(t 0 ) + O(h 2 ) = x 0 + khg(x 0 ) + O(h 2 ).
If we plug these results in the linear multi-step method,

T (h) = 1 h   x(t k+s ) + s-1 i=0 ρ i x(t k+i ) -h s i=0 σ i g(x(t k+i ))   = 1 h ρ(1)x 0 + (ρ (1) -σ(1))g(x 0 ) + O(h).
The limit is equal to zero if and only if we satisfy (2.8).

Consistency is crucial for integration methods, we give some intuition about how important are conditions defined in (2.8).

First condition, ρ(1) = 0. If the condition is not satisfied, then the method exhibits an artificial gain or damping. Assume we start at some equilibrium x * of ODE (i.e. ∇f (x * ) = 0), so x i = x * for the first s -1 steps. The next iterate becomes

x s = - s-1 i=0 ρ i x * + hσ(E) g(x * ) =0 ,
and if 1 +

s i=0 ρ i = ρ(1) = 0, we have that the next iterate x s is different from x * .
Second condition, ρ (1) = σ(1). If this relation is not satisfied, we actually are integrating another equation than (ODE) as shown in the following reasoning. Assume the first condition satisfied, 1 is a root of ρ. Consider then the factorization

ρ(z) = (z -1)ρ(z)
where ρ is a polynomial of degree s -1, and ρ (1) = ρ(1). The linear multi-step method becomes ρ(E)(y k+1y k ) = hσ(E)g(y k ).

If we sum up the above equation from the initial point, we get

ρ(E)(y k -y 0 ) = σ(E)G k ,
where G k = k i=0 hg(y i ). If h goes to zero, our iterates y k converge to a continuous curve y(t), and

G k → t 0 g(y(τ )) dτ , s i=0 ρi (y(t) -x(0)) = s-1 i=0 σ i t 0 g(y(τ )) dτ.
If we take the derivative over time, we get

ρ(1) ẏ(t) = σ(1)g(y(t)) ⇔ ρ (1) ẏ(t) = σ(1)g(y(t)).
which is different from the ODE we wanted to discretize, unless ρ (1) = σ(1).

Zero-stability and root condition

To get stability, assuming consistency holds as above, we also need to satisfy the zero-stability condition (2.6), which characterizes the sensitivity of a method to initial conditions. Actually, the name comes from an interesting fact: analyzing the special case where g = 0 is completely equivalent to the general case, as stated in the root condition theorem.

Theorem 2.3.5 (Root condition). Consider a linear multi-step method (ρ, σ). The method is zero-stable if and only if all roots of ρ are in the unit disk, and the roots on the unit circle are simple.

The proof of this theorem is technical and can be found as Theorem 6.3.4 of Gautschi [START_REF] Gautschi | Numerical analysis[END_REF]. The analysis becomes simplified, and reduces to standard linear algebra results because we only need to look at the solution of the homogeneous difference equation ρ(E)x k = 0.

Convergence of the global error and Dahlquist's equivalence theorem

Numerical analysis focuses on integrating an ODE on a finite interval of time [0, t max ]. It studies the behavior of the global error defined by the difference between the scheme and the solution of the ODE, i.e., x(t k )x k , as a function of the step size h. If the global error converges to 0 with the step size, the method is guaranteed to approximate correctly the ODE on the time interval, for h small enough. The faster it converges with h, the less points we need to guarantee a global error within a given accuracy.

We now state Dahlquist's equivalence theorem, which shows that the global error converges to zero when h does if the method is stable, i.e., when the method is consistent and zerostable. This naturally needs the additional assumption that the starting values x 0 , . . . , x s-1 are computed such that they converge to the solution (x(0), . . . , x(t s-1 )).

Theorem 2.3.6 (Dahlquist's equivalence theorem). Given an (ODE) defined by g and x 0 and a consistent linear multi-step method (ρ, σ), whose starting values are computed such that lim h→0 x i = x(t i ) for any i ∈ {0, . . . , s -1}, zero-stability is necessary and sufficient for being convergent, i.e., x(t k )x k tends to zero for any k when the step size h tends to zero.

Again, the proof of the theorem can be obtained from Gautschi [START_REF] Gautschi | Numerical analysis[END_REF]. Notice that this theorem is fundamental in numerical analysis. For example, it says that if a consistent method is not zero-stable, then the global error may be arbitrary large when the step size goes to zero, even if the local error decreases. In fact, if zero-stability is not satisfied, there exists a sequence generated by the linear multi-step method which grows with arbitrarily large factors.

Region of absolute stability

Stability and global error are ensured on finite time intervals, however solving optimization problems requires us to look at the infinite time horizon. We thus need more refined results and start by finding conditions ensuring that the numerical solution does not diverge when the time interval increases, i.e. that the numerical solution is stable with a constant C which does not depend of t max . Formally, for a fixed step-size h, we want to ensure

x k ≤ C max i∈{0,...,s-1}
x i for all k ∈ [0, t max /h] and t max > 0.

(2.9) This is not possible without assumptions on the function g as in the general case the solution x(t) itself may diverge. We begin with the simple scalar linear case which, given λ > 0, reads

ẋ(t) = -λx(t)
x(0) = x 0 .

(Scalar Linear ODE)

The recurrence of a linear multi-step methods with parameters (ρ, σ) applied to (Scalar Linear ODE) then reads

ρ(E)x k = -λhσ(E)x k ⇔ [ρ + λhσ](E)x k = 0,
where we recognize an homogeneous recurrence equation. Condition (2.9) is then controlled by the step size h and the constant λ, ensuring that this homogeneous recurrent equation produces bounded solutions. This leads us to the definition of the region of absolute stability.

Definition 2.3.7. The region of absolute stability of a linear multi-step method defined by (ρ, σ)

is the set of values λh such that the characteristic polynomial

π λh ρ + λhσ (2.10)
of the homogeneous recurrent equation π λh (E)x k = 0 produces bounded solutions.

Standard linear algebra links this condition to the roots of the characteristic polynomial as recalled in the next proposition (see Lemma 12.1 of Süli and Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF]).

Proposition 2.3.8. Let π be a polynomial and write x k a solution of the homogeneous recurrent equation π(E)x k = 0 with arbitrary initial values. If all roots of π are inside the unit disk and the ones on the unit circle have a multiplicity exactly equal to one, then x k ≤ ∞.

Absolute stability of a linear multi-step method determines its ability to integrate a linear ODE defined by

ẋ(t) = -Ax(t) x(0) = x 0 , (Linear ODE)
where A is a positive definite matrix whose eigenvalues belong to [µ, L] for 0 < µ ≤ L. In this case the step size h must indeed be chosen such that for any λ ∈ [µ, L], λh belongs to the region of absolute stability of the method. This (Linear ODE) is a special instance of (Gradient Flow)

where f is a quadratic function. Therefore absolute stability gives necessary (but not sufficient) condition to integrate (Gradient Flow) of L-smooth µ-strongly convex functions.

Convergence analysis in the linear case

By construction, absolute stability also gives us hints on the convergence of x k to the equilibrium.

More precisiely, it allows us to control the rate of convergence of x k , approximating the solution x(t) of (Linear ODE) (see Lemma 12.1 of Süli and Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF]).

Proposition 2.3.9. Given a (Linear ODE) defined by x 0 and a positive definite matrix A whose eigenvalues belong to [µ, L] for 0 < µ ≤ L, for a fixed step size h and a linear multi-step method defined by (ρ, σ), let r max be defined as

r max = max λ∈[µ,L] max r∈roots(π λh (z)) |r|,
where π λh is defined in (2.10). If r max < 1 and its multiplicity is equal to m, then the speed of convergence of the sequence x k produced by the linear multi-step method to the equilibrium x * of the differential equation is given by

x k -x * = O(k m-1 r k max ). (2.11)
We can now use these properties to analyze and design multi-step methods.

Analysis and design of multi-step methods

As shown before, we want to integrate (Gradient Flow) and Proposition 2.2.2 gives us a rate of convergence in the continuous case. If the method tracks x(t) with sufficient accuracy, then the rate of the method will be close to the rate of convergence of x(kh). So, larger values of h yield faster convergence of x(t) to the equilibrium x * . However h cannot be too large, as the method may be too inaccurate and/or unstable as h increases. Convergence rates of optimization algorithms are thus controlled by our ability to discretize the gradient flow equation using large step sizes. We recall the different conditions that proper linear multi-step methods should follow.

• Monic polynomial (Section 2.3.2). This is a convention, otherwise dividing both sides of the difference equation of the multi-step method by ρ s does not change the method.

• Explicit method (Section 2.3.2). We assume that the scheme is explicit in order to avoid solving a non-linear system at each step (Section 2.7 shows that implicit methods are linked to proximal methods).

• Consistency (Section 2.3.3). If the method is not consistent, then the local error does not converge when the step size goes to zero.

• Zero-stability (Section 2.3.3). Zero-stability ensures convergence of the global error (Section 2.3.4) when the method is also consistent.

• Region of absolute stability (Section 2.3.5). If λh is not inside the region of absolute stability for any λ ∈ [µ, L], then the method is divergent when t max increases.

Using the remaining degrees of freedom, we will tune the algorithm to have the best rate of convergence on (Linear ODE), which corresponds to the optimization of a quadratic function.

Indeed, as showed in Proposition 2.3.9, the largest root of π λh gives us the rate of convergence on quadratic functions (when λ ∈ [µ, L]). Since smooth and strongly convex functions are close to be quadratic (they are in fact sandwiched between two quadratics), this will also give us a good idea of the rate of convergence on these functions.

We do not derive a proof of convergence of the sequence for a general smooth and (strongly) convex function (in fact, it is already proved by [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] or using Lyapunov techniques by Wilson et al. [START_REF] Wilson | A Lyapunov analysis of momentum methods in optimization[END_REF]). But our results provide intuition on why accelerated methods converge faster.

Analysis and design of explicit Euler's method (s = 1)

In Section 2.3.1 we introduced Euler's method. In fact, we can view it as an explicit linear "multi-step" method with s = 1 defined by the polynomials

ρ(z) = -1 + z, σ(z) = 1.
We can check easily that it is consistent (using Proposition 2.3.4) and zero-stable since ρ has only one root which lies on the unit circle (Theorems 2.3.5 and 2.3.6). We need to determine the region of absolute stability in order to have an idea about the maximum value that h > 0 can take before the method becomes unstable. Assume we want to integrate any µ-strongly convex and L-smooth function f , with 0 < µ < L with any starting value x 0 . Then, we need to find the set of values of h such that the roots of the polynomial

π λh (z) = [ρ + λhσ](z) = -1 + λh + z, λ ∈ [µ, L]
are small. The unique root is 1λh and we need to solve the following minimax problem

min h max λ∈[µ,L] |1 -λh| ,
in the variable h > 0. The solution of this optimization problem is h * = 2 L+µ , its optimal value is (Lµ)/(L + µ) and its rate of convergence is then

x k -x * = O 1 -µ/L 1 + µ/L k .
We recover the optimal step size and the rate of convergence of the gradient method for a general smooth and strongly convex function [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF].

2.4.2 Analysis of two-step methods (s = 2)

We will now analyze two-step methods. First we write the conditions of a good linear multi step method, introduced at the beginning of this section, into constraints on the coefficients.

ρ 2 = 1 (Monic polynomial) σ 2 = 0 (Explicit method) ρ 0 + ρ 1 + ρ 2 = 0 (Consistency) σ 0 + σ 1 + σ 2 = ρ 1 + 2ρ 2 (Consistency) |Roots(ρ)| ≤ 1 (Zero-stability).
If we use all equalities, we finally have three linear constraints, defined by

L = ρ 0 , ρ 1 , σ 1 : ρ 1 = -(1 + ρ 0 ); σ 1 = 1 -ρ 0 -σ 0 ; |ρ 0 | < 1 . (2.12)
We will now try to find some condition on the remaining parameters in order to have a stable method. At first, let us analyze a condition on the roots of second order equations. Absolute stability requires that all roots of the polynomial π λh are inside the unit circle. The following proposition gives us the values of the roots of π λh as a function of the parameters ρ i and σ i .

Proposition 2.4.1. Given constants 0 < µ ≤ L, a step size h > 0 and a linear two-step method defined by (ρ, σ), under the conditions

(ρ 1 + µhσ 1 ) 2 ≤ 4(ρ 0 + µhσ 0 ), (ρ 1 + Lhσ 1 ) 2 ≤ 4(ρ 0 + Lhσ 0 ),
the roots r ± (λ) of π λh , defined in (2.10), are complex for any λ ∈ [µ, L]. Moreover, the largest modulus root is equal to

max λ∈[µ,L] |r ± (λ)| 2 = max {ρ 0 + µhσ 0 , ρ 0 + Lhσ 0 } . (2.13)
Proof. We begin by analyzing the roots r ± of the generic polynomial

z 2 + bz + c,
where b and c are real numbers, corresponding to the coefficients of π λh , i.e. b = ρ 1 + λhσ 1 and c = ρ 0 + λhσ 0 . For a fixed λ roots are complex if and only if

b 2 ≤ 4c ⇔ (ρ 1 + λhσ 1 ) 2 -4(ρ 0 + λhσ 0 ) ≤ 0.
Since the left side is a convex function in λ, it is equivalent to check only for the extreme values

(ρ 1 + µhσ 1 ) 2 ≤ 4(ρ 0 + µhσ 0 ), (ρ 1 + Lhσ 1 ) 2 ≤ 4(ρ 0 + Lhσ 0 ).
As roots are complex conjugates,

|r ± (λ)| 2 = |c| = |ρ 0 + λhσ 0 |.
Because this one-dimensional function is convex in λ, the maximum is attained for an extreme value of λ,

max λ∈[µ,L] |r ± (λ)| 2 = max {ρ 0 + µhσ 0 , ρ 0 + Lhσ 0 } , which is the desired result.
The next step is to minimize the largest modulus defined in (2.13) in the coefficients ρ i and σ i to get the best rate of convergence, assuming the roots are complex. We will not develop the case where the roots are real because this leads to weaker results.

Design of optimal two-step method for quadratics

We now have all ingredients to build a two-step method for which the sequence x k converges quickly to x * for quadratic functions. We need to solve the following problem,

minimize max {ρ 0 + µhσ 0 , ρ 0 + Lhσ 0 } s.t. (ρ 0 , ρ 1 , σ 1 ) ∈ L (ρ 1 + µhσ 1 ) 2 ≤ 4(ρ 0 + µhσ 0 ) (ρ 1 + Lhσ 1 ) 2 ≤ 4(ρ 0 + Lhσ 0 ), in the variables ρ 0 , ρ 1 , σ 0 , σ 1 , h > 0,
where L is defined in (2.12). If we use the equality constraints in (2.12) and make the following change of variables,

         ĥ = h(1 -ρ 0 ), c µ = ρ 0 + µhσ 0 , c L = ρ 0 + Lhσ 0 , (2.14)
the problem becomes, for fixed ĥ,

minimize max c µ , c L s.t. (-1 -c µ + µ ĥ) 2 ≤ 4c µ (-1 -c L + L ĥ) 2 ≤ 4c L |Lc µ -µc L | < |L -µ|,
in the variables c µ , c L . In that case, the optimal solution is given by

c * µ = 1 -µ ĥ 2 , c * L = 1 -L ĥ 2 , (2.15)
obtained by tightening the two first inequalities, for ĥ ∈]0, (1+µ/L) 2

L

[ such that the last inequality is satisfied. Now if we fix ĥ we can recover an optimal two step linear method defined by (ρ, σ) and an optimal step size h by using the equations in (2.14). We will use the following quantity

β 1 -µ/L 1 + µ/L . ( 2

.16)

A suboptimal two-step method. We can fix ĥ = 1/L for example. All computations done, the parameters of this two-step method, called method M 1 , are

M 1 =          ρ(z) = β -(1 + β)z + z 2 , σ(z) = -β(1 -β) + (1 -β 2 )z, h = 1 L(1-β) , (2.17)
and its largest modulus root (2.13) is given by rate

(M 1 ) = max{c µ , c L } = √ c µ = 1 -µ/L.
Optimal two-step method for quadratics. We can compute the optimal ĥ which minimizes the maximum of the two roots c * µ and c * L defined in (2.15). The solution is simply the one which balances the two terms in the maximum:

ĥ * = (1 + β) 2 L ⇒ c * µ = c * L .
This choice of ĥ leads to the method M 2 , described by

M 2 =            ρ(z) = β 2 -(1 + β 2 )z + z 2 , σ(z) = (1 -β 2 )z, h = 1 √ µL , (2.18)
with the rate of convergence

rate(M 2 ) = √ c µ = √ c L = β < rate(M 1 ).
We will now see that methods M 1 and M 2 are actually related to Nesterov's method and Polyak's heavy ball algorithms.

On the link between integration and optimization

In the previous section, we derived a family of linear multi-step methods, parametrized by ĥ. We will now compare these methods to common optimization algorithms used to minimize L-smooth, µ-strongly convex functions.

Polyak's heavy ball method

The heavy ball method was proposed by Polyak [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF]. It adds a momentum term to the gradient step

x k+2 = x k+1 -c 1 ∇f (x k+1 ) + c 2 (x k+1 -x k ),
where

c 1 = (1 -β 2 )/ √ µL and c 2 = β 2
where β is defined in (2.16). We can organize the terms in the sequence to match the general structure of linear multi-step methods, to get

β 2 x k -(1 + β 2 )x k+1 + x k+2 = (1 -β 2 )/ µL -∇f (x k+1 ) . We easily identify ρ(z) = β 2 -(1 + β 2 )z + z 2 and hσ(z) = (1 -β 2 )/ √ µLz.
To extract h, we will assume that the method is consistent (see conditions (2.8)), which means

ρ(1) = 0 Always satisfied hρ (1) = hσ(1) ⇒ h = 1 √ µL .
All computations done, we can identify the "hidden" linear multi-step method as

M Polyak =            ρ(z) = β 2 -(1 + β 2 )z + 1 σ(z) = (1 -β 2 )z h = 1 √ µL .
(2. [START_REF] Deng | Recent advances in deep learning for speech research at Microsoft[END_REF] This shows that M Polyak = M 2 . In fact, this result was expected since Polyak's method is known to be optimal for quadratic functions. However, it is also known that Polyak's algorithm does not converge for a general smooth and strongly convex function [START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF].

Nesterov's accelerated gradient

Nesterov's accelerated method in its simplest form is described by two sequences x k and y k , with

y k+1 = x k - 1 L ∇f (x k ), x k+1 = y k+1 + β(y k+1 -y k ).
As above, we will write Nesterov's accelerated gradient as a linear multi-step method by expanding y k in the definition of x k , to get

βx k -(1 + β)x k+1 + x k+2 = 1 L -β(-∇f (x k )) + (1 + β)(-∇f (x k+1 )) .
Consistency of the method is then ensured by

ρ(1) = 0 Always satisfied hρ (1) = hσ(1) ⇒ h = 1 L(1 -β) .
After identification,

M Nest =          ρ(z) = β -(1 + β)z + z 2 , σ(z) = -β(1 -β) + (1 -β 2 )z, h = 1 L(1-β) ,
which means that M 1 = M Nest .

Nesterov's method interpretation as a faster stable integration method

Pushing the analysis a little bit further, we can show why Nesterov's algorithm is faster than gradient method. There is of course a complete proof of its rate of convergence [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF], even using the argument of differential equations [START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF][START_REF] Wilson | A Lyapunov analysis of momentum methods in optimization[END_REF], but we take a more intuitive approach here. The key parameter is the step size h. If we compare it with the one in classical gradient method, Nesterov's method uses a step size which is

(1 -β) -1 ≈ L/µ larger.
Recall that, in continuous time, we have seen the rate of convergence of x(t) to x * given by

f (x(t)) -f (x * ) ≤ e -2µt (f (x 0 ) -f (x * )).
The gradient method tries to approximate x(t) using Euler approximation with step size h = 1/L, which means x

(grad) k ≈ x(k/L), so f (x (grad) k ) -f (x * ) ≈ f (x(k/L)) -f (x * ) ≤ (f (x 0 ) -f (x *
))e -2k µ L . On the left, the optimal step size is used. Because Polyak's algorithm is the one with the biggest step size, it needs less iterations than Nesterov or Euler to approximate x(t max ). On the right side, we reduced the step size to 1/L for all methods. We clearly observe that they all track the same ODE: the gradient flow.

However, Nesterov's method has the step size

h Nest = 1 L(1 -β) = 1 + µ/L 2 √ µL ≈ 1 √ 4µL which means x nest k ≈ x k/ 4µL .
In that case, the estimated rate of convergence becomes

f (x nest k ) -f (x * ) ≈ f (x(k/ 4µL)) -f (x * ) ≤ (f (x 0 ) -f (x * ))e -k √ µ/L ,
which is approximatively the rate of convergence of Nesterov's algorithm in discrete time and we recover the accelerated rate in µ/L versus µ/L for gradient descent. The accelerated method is more efficient because it integrates the gradient flow faster than simple gradient descent, making longer steps. A numerical simulation in Figure 2.1 makes this argument more visual. This intuitive argument is still valid for the convex counterpart of Nesterov's accelerated gradient.

Acceleration for convex functions

By matching the coefficients of Nesterov's method, we deduced the value of the step-size used for the integration of (Gradient Flow). Then, using the rate of convergence of x(t) to x * , we estimated the rate of convergence of Nesterov's method assuming x k ≈ x(t k ). Here, we will do the same but without assuming strong convexity. However, the estimation of the rate of convergence in discrete time needs the one in continuous time, described by the following proposition.

Proposition 2.6.1. Let f be L-smooth convex function, x * one of its minimizers and x(t) be the solution of (Gradient Flow). Then

f (x(t)) -f (x * ) ≤ x 0 -x * 2 t + (2/L) . (2.20) Proof. Let L(x(t)) = f (x(t))-f (x * ). We notice that ∇L(x(t)) = ∇f (x(t)) and dL(x(t))/ dt = -∇f (x(t)) 2 . Since f is convex, L(x(t)) ≤ ∇f (x(t)), x(t) -x * ≤ ∇f (x(t)) x(t) -x * .
By consequence,

-∇f (x(t)) 2 ≤ - L(x(t)) 2 x(t) -x * 2 ≤ - L(x(t)) 2 x 0 -x * 2 .
(2.21)

The last inequality comes from the fact that x(t)x * decreases over time,

d dt x(t) -x * 2 = 2 ẋ(t), x(t) -x * , = -2 ∇f (x(t)), x(t) -x * , ≤ 0 since f is convex.
From (2.21), we deduce the differential inequality

d dt L(x(t)) ≤ - L(x(t)) 2 x 0 -x * 2 .
The solution is obtained by integration,

t 0 dL(x(τ ))/ dτ L(x(τ )) 2 dτ ≤ t 0 -1 x 0 -x * 2 .
The general solution is thus

L(x(t)) ≤ 1 t x0-x * 2 + C ,
for some constant C > 0. Since the inequality is valid for all time t ≥ 0, the following condition on C,

L(x(t)) ≤ 1 t x0-x * 2 + C ≤ 1 C for t ≥ 0, is sufficient. Since L(x(t)) = f (x(t)) -f (x * ), setting C = 1 f (x0)-f (x * ) satisfies the above inequal- ity. Using smoothness of f , f (x 0 ) -f (x * ) ≤ L 2 x 0 -x * 2 ,
we get the desired result.

Assume we use Euler's method with step size h = 1 L , the estimated rate of convergence will be

f (x k ) -f (x * ) ≈ f (x(kh)) -f (x * ) ≤ L x 0 -x * 2 k + 2 ,
which is close to the rate of convergence of the classical gradient method for convex function.

Now, consider Nesterov's method for minimizing a smooth and convex function f :

y k+1 = x k - 1 L ∇f (x k ) x k+1 = -β k x k + (1 + β k )x k+1 ,
where

β k ≈ k-2 k+1 .
If we expand everything, we get after rearrangement,

β k x k-1 -(1 + β k )x k + x k+1 = 1 L β k (-∇f (x k-1 )) -(1 + β k )(-∇f (x k )) .
In other terms, we have an expression of the form

ρ k (E)x k = h k σ k (E)(-∇f (x k )).
We can identify h if we assume the method consistent, which means

ρ(1) = 0 Always satisfied h k ρ k (1) = h k σ k (1) ⇒ h k = 1 L(1 -β k+1 ) = (k + 2) 3L .
We can estimate, the rate of convergence of Nesterov's method using (2.20). Since x k ≈ x(t k ),

x k ≈ x k i=0 h i ≈ x k 2 6L .
In terms of convergence to the optimal value,

f (x k ) -f (x * ) ≈ f (x(t k ) -f (x * ) ≤ 6L x 0 -x * 2 k 2 + 12 ,
which is close to the bound from Nesterov [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. Again, because the step-size of Nesterov's algorithm is larger (while keeping a stable sequence), we converge faster than the Euler's method.

Proximal algorithms and implicit integration methods

We present here links between proximal algorithms and implicit numerical methods that integrate the gradient flow equation. We begin with Euler's implicit method that corresponds to the proximal point algorithm.

Euler's implicit method and proximal point algorithm

We saw in Section 2.3.1 that Euler's explicit method used the Taylor expansion of the solution x(t) of the (ODE) at the current point. The implicit version uses the Taylor expansion at the next point which reads

x(t) = x(t + h) -h ẋ(t + h) + O(h 2 ).
If t = kh, by neglecting the second order term we get implicit Euler's method,

x k+1 = x k + hg(x k+1 ). (2.22)
This recurrent equation requires to solve an implicit equation at each step that may be costly.

However it provides in some cases better stability than the explicit version, especially when the differential. This is generally the case for implicit methods (see Süli and Mayers [START_REF] Süli | An introduction to numerical analysis[END_REF] for further details on implicit methods).

Now assume that g comes from a potential -f such that we are integrating (Gradient Flow).

Solving the implicit equation (2.22) is equivalent to compute the proximal operator of f defined

as prox f,h (x) = argmin z 1 2 z -x 2 2 + hf (z). (2.23)
This can be easily verified by checking the first-order optimality conditions of the minimization problem. Euler's implicit method applied to (Gradient Flow) reads then

x k+1 = prox f,h (x k ),
where we recognize the proximal point algorithm [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF].

We present now Mixed ODE that corresponds to composite optimization problems.

Implicit Explicit methods and proximal gradient descent

In numerical analysis, it is common to consider the differential equation

ẋ = g(x) + ω(x), (Mixed ODE)
where g(x) is considered as the "non-stiff" part of the problem and ω the stiff one, where stiffness may be assimilated to bad conditioning [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Frank | On the stability of implicit-explicit linear multistep methods[END_REF]. Usually, we assume ω integrable using an implicit method.

If ω derives from a potential -Ω (meaning ω = -∇Ω), this is equivalent to assume that the proximal operator of Ω defined in (2.23) can be computed exactly.

We approximate the solution of (Mixed ODE) using IMplicit-EXplicit schemes (IMEX). In our case, we will focus on the following multi-step based IMEX scheme,

ρ(E)x k = h σ(E)g(x k ) + γ(E)ω(x k ) ,
where ρ, σ and γ are polynomials of degrees s, s -1 (the explicit part) and s respectively and ρ is monic. It means that, at each iteration, we need to solve, in x k+s ,

x k+s = s-1 i=0 -ρ i x k+i + σ i hg(x k+i ) + γ i hω(x k+i ) known +γ s ω(x k+s ).
In terms of optimization the mixed ODE corresponds to composite minimization problems of the form

minimize f (x) + Ω(x), (2.24)
where f, Ω are convex and Ω has a computable proximal operator. We can link IMEX schemes with many optimization algorithms which use the proximal operator, such as proximal gradient method, FISTA or Nesterov's method. For example, proximal gradient is written

y k+1 = x k -h∇f (x k )
x k+1 = prox hΩ (y k+1 ).

After expansion, we get

x k+1 = y k+1 -h∇Ω(x k+1 ) = x k + hg(x k ) + hω(x k+1 ),
which corresponds to the IMEX method with polynomials

ρ(z) = -1 + z, σ(z) = 1, γ(z) = z.
However, for Fista and Nesterov's method, we need to use a variant of linear multi-step algorithms, called one leg methods [START_REF] Dahlquist | On one-leg multistep methods[END_REF][START_REF] Zhang | Stability and convergence analysis of implicit-explicit one-leg methods for stiff delay differential equations[END_REF]. Instead of combining the gradients, the idea is to compute g at a linear combination of the previous points, i.e.

ρ(E)x k = h g(σ(E)x k ) + ω(γ(E)x k ) .
Their analysis (convergence, consistency, interpretation of h, etc...) is slightly different from linear multi-step method, so we will not go into details in this chapter, but the correspondence still holds.

Non-smooth gradient flow

In the last subsection we assumed that ω comes from a potential. However in the optimization literature, composite problems have a smooth convex part and a non-smooth sub-differentiable convex part which prevents us from interpreting the problem with the gradient flow ODE. Nonsmooth convex optimization problems can be treated with differential inclusions (see [START_REF] Bolte | The Łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems[END_REF] for recent results on it)

ẋ(t) + ∂f (x(t)) 0,
where f is a sub-differentiable function whose sub-differential at x is written ∂f (x). Composite problems (2.24) can then be seen as the discretization of the differential inclusion ẋ(t) + ∇f (x(t)) + ∂Ω(x(t)) 0.

Mirror gradient descent and non-Euclidean gradient flow

In many optimization problems, it is common to replace the Euclidean geometry with a distancegenerating function called d(x), with the associated Bregman divergence

B d (x, y) = d(x) -d(y) -∇d(y), x -y ,
with d strongly-convex and lower semi-continuous. To take into account this geometry we consider the Non-Euclidean Gradient Flow [START_REF] Krichene | Accelerated mirror descent in continuous and discrete time[END_REF]

ẏ(t) = -∇f x(t) x(t) = ∇d * (y(t))
x(0) = x 0 , y(0) = ∇d(x 0 ).

(NEGF)

Here ∇d maps primal variables to dual ones and, as d is strongly convex, (∇d) -1 = ∇d * , where d * is the Fenchel conjugate of d. In fact, we can write (NEGF) using only one variable y, but this formulation has the advantage to exhibit both primal and dual variables x(t) and y(t). Applying the forward Euler's explicit method we get the following recurrent equation

y k+1 -y k = -h∇f (x k ), x k+1 = ∇d * y k+1 .
Now consider the mirror gradient scheme :

x k+1 = argmin x h ∇f (x k ), x + B d (x, x k ).
First optimality condition reads

∇ x h ∇f (x k ), x + B d (x, x k ) x=x k+1 = h∇f (x k ) + ∇d(x k+1 ) -∇d(x k ) = 0 Using that (∇d) -1 = ∇d * we get h∇f (x k ) + y k+1 -y k = 0, x k+1 = ∇d * y k+1 ,
which is exactly Euler's explicit method defined in (NEGF).

Universal gradient descent and generalized gradient flow

Consider the Generalized Gradient Flow, which combines the ideas of (Mixed ODE) and (NEGF),

ẏ(t) = -∇f (x(t)) -∇Ω(x(t))
x(t) = ∇d * (y(t))

x(0) = x 0 , y(0) = ∇d(x 0 ).

(GGF)

We can write its ODE counterpart, called the "Generalized ODE",

ẏ(t) = g(x(t)) + ω(x(t))
x(t) = ∇d * (y(t)),

x(0) = x 0 , y(0) = ∇d(x 0 ).

(GODE)

where g = -∇f , with f a smooth convex function, d a strongly convex and semi-continuous distance generating function and ω = -∂Ω, where Ω is a simple convex function. If Ω is not differentiable we can consider the corresponding differential inclusion as presented in Section 2.7.3. Here we focus on (GODE) and (GGF) to highlight the links with integration methods.

The discretization of this ODE is able to generate many algorithms in many different settings.

For example, consider the IMEX explicit-implicit Euler's method,

y k+1 -y k h = g(x k ) + ω(∇d * (y k+1 )), x k+1 = ∇d * (y k+1 ),
which can be decomposed into three steps,

z k+1 = y k + hg (x k ) (Gradient step in dual space), y k+1 = prox h(Ω•∇d * ) (z k+1 ) (Projection step in dual space), (2.25)
x k+1 = ∇d * (y k+1 ) (Mapping back in primal space). Now consider the universal gradient method scheme presented by Nesterov [START_REF] Nesterov | Universal gradient methods for convex optimization problems[END_REF]:

x k+1 = arg min x ∇f (x k ), x -x k + Ω(x) + B d (x, x k ).
Again we can show that both recursions are the same: if we write the first optimality condition,

0 = ∇ x h ∇f (x k ), x -x k + hΩ(x) + B(x, x k ) x=x k+1 = hg(x k ) + h∂Ω(x k+1 ) + ∇d(x k+1 ) -∇d(x k ) = hg(x k ) -y k =z k+1 +h∂Ω(∇d * (y k+1 )) -y k+1 .
We thus need to solve the non-linear system of equations

y k+1 = z k+1 + h∂Ω(∇d * (y k+1 )),
which is equivalent to the projection step (2.25). Then we simply recover x k+1 by applying ∇d * on y k+1 .

Conclusion and future works

We connected several optimization algorithms to multi-step integration methods in numerical analysis. By using the theory of linear multi-step methods on the basic gradient flow equation, we recover Polyak's and Nesterov's method using some optimality arguments. This provides an intuitive interpretation for the design of these methods, with optimal step sizes giving a direct explanation of the acceleration phenomenon. Our approach generalizes to more structured problems by introducing the appropriate integration method and/or looking at a generalized gradient flow equation that takes into account the geometry of the problem.

We described a simple interpretation of the acceleration phenomenon, but our analysis is still restricted to quadratic problems. The study of G-stability [START_REF] Dahlquist | G-stability is equivalent toa-stability[END_REF][START_REF] Butcher | Thirty years of g-stability[END_REF] may generalize our approach to smooth strongly convex functions. For the non-strongly convex case, the dependence in k of Nesterov's algorithm makes its links with integration methods less clear.

Finally the study of the links between optimization and numerical methods may provide new algorithms for both fields. Conversely, Nesterov's algorithm may lead to new integration methods to integrate the gradient flow equation of a convex function.

Chapter 3

Regularized Nonlinear Acceleration

Abstract

We describe a generic convergence acceleration technique for iterative algorithms. Our scheme computes estimates of fixed-points from a nonlinear average of the iterates produced by the iterative method. The weights in this average are computed via a simple linear system, whose solution can be updated online. This acceleration scheme runs in parallel to the base algorithm, providing improved estimates of the solution on the fly, while the original method is running.

In this chapter, we study different settings such as linear algorithms with arbitrary, potentially random, perturbations. We also propose an acceleration method when another linear combination with arbitrary coefficients is maintained at each iteration, which usually corresponds to algorithms with momentum or averaging.

Main references: The results in this chapter have been partially published in [START_REF] Scieur | Regularized nonlinear acceleration[END_REF], [START_REF] Scieur | Nonlinear acceleration of stochastic algorithms[END_REF] (NIPS 2016 and 2017) and [START_REF] Scieur | Nonlinear acceleration of deep neural networks[END_REF] (under review NIPS 2018). A revised version of [START_REF] Scieur | Regularized nonlinear acceleration[END_REF] is under review for the journal Mathematical Programming. In this thesis, in comparison with the published papers, we improved the theoretical analysis of the proposed algorithm, especially the convergence bounds.

Introduction

Suppose we run an iterative algorithm of the form

x i+1 = g(x i )
which converges to a solution x * . For example, such an algorithm is used in optimization for solving the minimization problem

min x∈R d f (x). (3.1)
Here, we will focus on improving our estimates of the solution to problem (3.1) by tracking only the iterate sequence {x i } produced by an optimization algorithm, without any further calls to oracles on g(x).

We will study methods for accelerating convergence of iterative algorithms, originally studied in from numerical analysis. In particular we propose a regularized version of Anderson's extrapolation [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]. This method uses the iterates produced by any (converging) iterative algorithm, and estimates the solution directly from this sequence, up to some regularity conditions. Our scheme is based on the idea behind Aitken's ∆ 2 -algorithm [START_REF] Aitken | On bernoulli's numerical solution of algebraic equations[END_REF], generalized as the Shanks transform [START_REF] Shanks | Non-linear transformations of divergent and slowly convergent sequences[END_REF],

whose recursive formulation is known as the ε-algorithm [START_REF] Wynn | On a device for computing the e m (s n) transformation[END_REF] (see e.g. [START_REF] Brezinski | Accélération de la convergence en analyse numérique[END_REF][START_REF] Sidi | Acceleration of convergence of vector sequences[END_REF] for a survey). In a nutshell, these methods fit geometrical models to linearly converging scalar sequences, then extrapolate their limit from the fitted model. Our method can be viewed as an extension to vector sequences.

We start from a formulation of these techniques known as anderson Acceleration [4], Mešina's algorithm [START_REF] Mešina | Convergence acceleration for the iterative solution of the equations x= ax+ f[END_REF] or minimal polynomial extrapolation (MPE) [START_REF] Sidi | Acceleration of convergence of vector sequences[END_REF][START_REF] Smith | Extrapolation methods for vector sequences[END_REF]. They use the minimal polynomial of the linear operator driving iterations to estimate the optimum by a nonlinear average of the iterates (i.e. computing a weighted average using weights which are nonlinear functions of the iterates).

Our contribution here is to regularize this procedure and produce explicit bounds on the distance to optimality by controlling stability, thus explicitly quantifying acceleration. We show that these extrapolation algorithms reach optimal performance (asymptotically) and describe

several numerical examples where these stabilized estimates often speed up convergence by an order of magnitude. So far, for all the techniques cited above, no proofs of convergence of the estimates were given when the estimation process became unstable (for example when extrapolating iterates produced by nonlinear functions). Furthermore, the acceleration scheme can be run in parallel with the original algorithm, providing improved estimates of the solution on the fly, while the original method is progressing, so its numerical complexity is marginal.

Notations For the sake of simplicity, throughout the thesis we use • for the 2 norm of a vector, and its corresponding induced matrix norm

A = max v: v =1
Av .

In this chapter, we distinguish perturbed iteration or operators with their "noiseless" counterpart using the tilde ˜notation. In particular, we distinguish two sequences of points generated by

x i+1 = g(x i ), xi+1 = g(x i ).
where g and g share the same fixed-point, x * . On the left, the function g is written

g(x) = G(x -x * ) + x *
where G is a matrix. Its perturbed version reads

g(x i ) = g(x i ) + e i
for some perturbation e i .

Convergence Acceleration for Linear Algorithms

We begin by recalling the core arguments behind convergence acceleration. These ideas have taken various forms over time, known for example as Anderson acceleration [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF], the Eddy-Mesina method [START_REF] Mešina | Convergence acceleration for the iterative solution of the equations x= ax+ f[END_REF][START_REF] Eddy | Extrapolating to the limit of a vector sequence[END_REF] and minimal polynomial extrapolation [START_REF] Cabay | A polynomial extrapolation method for finding limits and antilimits of vector sequences[END_REF][START_REF] Smith | Extrapolation methods for vector sequences[END_REF]. The core idea behind these methods is to use a Taylor expansion of

xi+1 = g(x i ) = G(x i -x * ) + e i (3.2)
to approximate the fixed point iterations by a vector autoregressive model, then compute a weighted mean of the iterates x i to produce a better estimate of the (unique) limit x * . In this section, we will study the rate of convergence of these algorithms when g is linear, i.e., when

x i+1 is generated with the linear fixed point iteration

x i+1 = g(x i ) = G(x i -x * ) + x * , (LFPI)
where we assume G to be symmetric, positive semi-definite and G < (1κ), where κ < 1.

When optimizing a quadratic function using the gradient method, κ is often refereed to be the inverse of the condition number.

This model corresponds, for example, to a linearization of (3.2). Suppose g(x) is differentiable and let G be the Jacobian of g evaluated at x * . If we linearize (3.2) around x * we obtain

x i+1 = g(x * ) + G(x i -x * ) + O( x i -x * 2 ), for i = 1, . . . , k.
By neglecting the second order term, and because g(x * ) = x * , we obtain the linear fixed-point iteration (LFPI).

Because G 2 ≤ (1κ) < 1, the iterates x i converge to x * at a linear rate, with

x i -x * ≤ (1 -κ) x i-1 -x * ≤ (1 -κ) i x 0 -x * ,
where • stands for the Euclidean norm. Similarly, for linear g, the residual

r i g(x i ) -x i (= x i+1 -x i ) (3.3)
converges with the same rate since it also follow the autoregressive process

r i = (x i+1 -x * ) -(x i -x * ), = G (x i -x * ) -(x i-1 -x * ) , = Gr i-1 , (3.4)
so we easily deduce that

r i ≤ G r i-1 ≤ (1 -κ) i r 0 .
We will now see how to improve convergence rates using a linear combination of the previous iterates. Suppose we run N iterations of (LFPI). We will see that a good linear combination called extrapolation,

x extr = N i=0 c i x i , (3.5)
can minimize the residual (3.3). In other words, we have to find c * which satisfies

c * = argmin c g (x extr ) -x extr . (3.6)
The next proposition shows that minimizing the combination of residuals (3.3) leads to the same result, if the coefficients c sum to one. In particular, this implies

g (x extr ) -x extr = N i=0 c * i r i . (3.8)
Proof. Equation (3.6) can be expanded using (3.5) and (LFPI),

g (x extr ) -x extr = G   N i=0 c i x i -x *   + x * - N i=0 c i x i .
Since N i=0 c i = 1, we can factor the sum,

g (x extr ) -x extr = N i=0 c i G (x i -x * ) + x * =g(xi) -x i .
By definition of the residual (3.3), we obtain

g (x extr ) -x extr = N i=0 c i r i .
This conclude the proof, since the objective functions of (3.6) and (3.7) are equal.

We can now minimize the residual of the extrapolation x extr . It suffices to minimize the combination of the residuals r i with coefficients c i , then use these coefficients to combine the points x i of the sequence. The next theorem gives an upper bound on the residual norm of extrapolation given some assumptions on the spectrum of the matrix G in (LFPI).

Theorem 3.2.2. (Scieur, d'Aspremont and Bach [START_REF] Scieur | Regularized nonlinear acceleration[END_REF]). Suppose the iterates x i for i = 0, . . . , N are computed using (LFPI), with G symmetric positive semi-definite, satisfying G ≤ (1κ) for κ < 1. Let x * be the unique fixed point of g and m be the number of distinct eigenvalues of G. Then, the 2 norm of the residual of the extrapolation x extr (3.5) is bounded by

g(x extr ) -x extr 2 ≤      2β N 1 + β 2N r 0 2 if N < m 0 otherwise (3.9)
where

β = 1 - √ κ 1 + √ κ < 1. (3.10)
Proof. First, we rewrite the residual of x extr using polynomials. In view of Proposition 3.2.1 and equation (3.4),

g(x extr ) -x extr = N i=0 c i r i = N i=0 c i G i =p(G) r 0 .
The term p(G) is a polynomial in the matrix G, with coefficients c i . In this proof, we bound the value

min p∈R N [x]:p(1)=1 p(G)r 0 = min c:1 T c=1 g(x extr ) -x extr , where R N [x]
is the linear space of polynomials of degree at most N .

Because G is symmetric, it admits the eigenvalue decomposition

G = Q * ΛQ
where Λ is the diagonal matrix of eigenvalues {λ i , i = 1, . . . , m}, and Q is an orthogonal matrix.

For any polynomial p ∈ R N [x] (the space of polynomials of degree at most N ), we have

p(G)r 0 = Q * p(Λ)Qr 0 ≤ p(Λ) r 0 = max i=1,...,m |p(λ i )| r 0 .
There are two cases to analyze here, which are N < m and N ≥ m.

Case N ≥ m. Using Cayley-Hamilton theorem, if p(x) is the characteristic polynomial of G,

then p(G) = p(Λ) = 0.
By assumption none of the λ i is equal to 1 (we assumed

λ i ∈ [0, 1 -κ]
with 0 < κ < 1) so we can normalize p so that p(1) = 1 and p(λ i ) = 0 for all i = 1, . . . , m.

Case N < m. We have, for any polynomial p,

max i=1,...,m |p(λ i )| ≤ max λ∈[λmin,λmax] |p(λ)| Because 0 G 1 -κ, we have 0 ≤ λ i ≤ 1 -κ, so the error bound becomes min {p∈R N [x]: p(1)=1} p(G)(x 0 -x * ) 2 ≤ min {p∈R N [x]: p(1)=1} max λ∈[0,1-κ] |p(λ)| (x 0 -x * ) 2 (3.11)
where the right hand side involves a minmax problem, explicitly solved using Chebyshev's polynomials. Let C N (x) be the Chebyshev polynomial of degree N . By definition, C N is a monic polynomial (i.e. a polynomial whose leading coefficient is one) solving

C N (x) argmin {p∈R N [x]: p=1} max x∈[-1,1] |p(x)|.
Golub and Varga [START_REF] Golub | Chebyshev semi-iterative methods, successive overrelaxation iterative methods, and second order Richardson iterative methods[END_REF] use a variant of C N (x) to solve the problem in (3.11), whose solution is a rescaled Chebyshev polynomial given by

T N (x, κ) = C N (t(x, κ)) C N (t(1, κ)) , where t(x, κ) = 2x -(1 -κ) 1 -κ , (3.12)
where t(x, κ) is simply a linear mapping from interval [0, 1κ] to [-1, 1]. Moreover, they show that the maximal value is reached at the largest value of the interval,

min {p∈R N [x]: p(1)=1} max λ∈[0,1-κ] |p(λ)| = max λ∈[0,1-κ] |T N (λ, κ)| = |T N (1 -κ, κ)| = 2β N 1 + β 2N , (3.13)
where β is given by

β = 1 - √ κ 1 + √ κ < 1 -κ < 1.
Injecting the result of (3.13) in (3.11) yields the desired result.

Since in (3.7) we involve r N in the computation of the coefficients c * , it implies we have computed x N +1 because r N = x i+1x i . However, x extr does not contain the information of the last point of the sequence because x extr = N i=0 c * i x i . This can be considered as a waste of information, since x N +1 is more likely closer to the optimum x * . The following corollary shows that using the same coefficients c i we can incorporate x N +1 in the extrapolation process, and recover g(x extr ) instead of x extr .

Corollary 3.2.3. If g is a linear function, we have

g(x extr ) = N +1 i=1 c * i x i .
If we use this point instead of x extr , we have the following bound,

g g(x extr ) -g(x extr ) ≤ (1 -κ) g(x extr ) -x extr .
Proof. By linearity of g, and because 1 T c * = 1, we have

N +1 i=1 c * i x i = N i=0 c * i g(x i ) = g   N i=0 c * i x i   = g(x extr ).
The proof of Theorem 3.2.2 suggests T N (x, κ) is a good universal solution for the convergence acceleration problem. For illustration, we plot both T 3 (x, κ) and T 5 (x, κ) for x ∈ [0, 1] and κ = 0.15 in Figure 3.1. This solution is called the Chebyshev semi-iterative method in [START_REF] Golub | Chebyshev semi-iterative methods, successive overrelaxation iterative methods, and second order Richardson iterative methods[END_REF] and was further studied by e.g. [START_REF] Nemirovskii | Iterative methods for solving linear ill-posed problems under precise information[END_REF]. Combining the N iterates x i using the coefficients c i in T N (x, κ), we ensure

g(x extr ) -x extr ≤ (1 - √ κ) N r 0 (1 -κ) N r 0 ,
which means convergence is indeed accelerated. However, this method has some key drawbacks.

First, we need to know κ to form T N (x, κ), which is not always the case. For example, when g is in fact nonlinear, κ depends on the spectrum of the Jacobian at the optimum, which is clearly not observed. Second, the algorithm does not allow us to control the magnitude of the coefficients in the polynomial T (x, κ), which has a strong impact on the stability of the algorithm in the presence of numerical errors, or when the iterates are generated by a non-linear function g.

Because of stability issues with Chebyshev acceleration, we focus instead on the method (3.7) which computes the c i adaptively (which also has the same bound of convergence). We describe it formally in Algorithm 1. This acceleration algorithm is called nonlinear because the coefficients c i vary with x i . This method is also known as Anderson acceleration [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF], the Eddy-Mesina algorithm [START_REF] Mešina | Convergence acceleration for the iterative solution of the equations x= ax+ f[END_REF][START_REF] Eddy | Extrapolating to the limit of a vector sequence[END_REF], minimal polynomial extrapolation [START_REF] Cabay | A polynomial extrapolation method for finding limits and antilimits of vector sequences[END_REF], or reduced rank extrapolation

Algorithm 1 Nonlinear Acceleration of Convergence

Require:

Iterates x 0 , x 1 , . . . , x N +1 ∈ R d . 1: Form R = [r 0 , ..., r N ],
where r i is defined in (3.3). [ [START_REF] Sidi | Acceleration of convergence of vector sequences[END_REF][START_REF] Smith | Extrapolation methods for vector sequences[END_REF]. There are small variations between all these methods, which lie in the way they solve the minimization problem in (3.7).

In a sense, theses approaches are more statistical in nature. They assume an approximately linear model holds for iterations near the optimum, and estimate this model using the iterates.

When working with scalar sequences, it has strong links with Wynn's algorithm [START_REF] Wynn | On a device for computing the e m (s n) transformation[END_REF], directly connected to the Levinson-Durbin algorithm [START_REF] Levinson | The wiener rms error criterion in filter design and prediction, appendix b of wiener[END_REF][START_REF] Durbin | The fitting of time-series models[END_REF] used to solve Toeplitz systems recursively and fit autoregressive models (the Shanks transform solves Hankel rather than Toeplitz systems, but this is essentially the same problem [START_REF] Heinig | Fast algorithms for Toeplitz and Hankel matrices[END_REF]). The key difference in these extrapolation techniques is that estimating the autocovariance operator G is not required, as we only focus on the limit.

We see that Algorithm in the variable c ∈ R N , where R is a d × N matrix assumed to be of rank N , is given by

c * = (R T R) -1 1 1 T (R T R) -1 1 . (3.15)
Proof. Let µ be the dual variable of the equality constraint. Both c * and µ * should satisfy the KKT system

2R T R 1 1 T 0 c * µ * = 0 1 (3.16)
This block matrix can be inverted explicitly, with

2R T R 1 1 T 0 -1 = 1 1 T (R T R) -1 1 1 2 (R T R) -1 1 T (R T R) -1 1I -11 T (R T R) -1 (R T R) -1 1 1 T (R T R) -1 -2 .
Using this inverse we easily solve the linear system, which gives the result in (3.15).

In practice, instead of computing the inverse of the matrix R T R, we solve the linear system

R T Rz = 1, then set c * = z/(1 T z).
Of course, we need to assume the matrix invertible, and this likely holds until N > m in Theorem 3.2.2. In the case where the matrix is singular, it is possible to show, using the pseudo-inverse of R T R, that any solution of the system recovers exactly x * .

The formula (3.14) is used in Anderson acceleration algorithm and in the Mesina method.

Other algorithms usually force the last coefficient c N to be equal to one, solve the remaining linear system, then normalize the vector. However, these alternative strategies are harder to analyze when the iterates are generated by a non-linear function g.

Regularized Nonlinear Acceleration of Convergence

So far, we have only considered linear functions g in (LFPI), without perturbations, when computing the iterates x i . In general, the fixed-point iteration (3.2) is usually generated by a nonlinear function g, thus inducing a second order error term compared to the dynamics in (LFPI). We can also assume each iteration is corrupted by a random noise. In this section, we will study the impact of perturbations on the acceleration algorithm and propose a new variant which is robust to noise.

Sensitivity Analysis

We now study the sensitivity of the acceleration algorithm to perturbations. Consider the following perturbed linear fixed point iteration

xi+1 -x * g(x i ) -x * = G(x i -x * ) + e i (Pert. LFPI)
where e i is the noise injected at iteration i. For now, we do not assume any structure on the noise, so e i may come from nonlinearities of g, stochastic noise, roundoff error, etc. The iterates of this process will be compared to their noiseless counterpart,

x i+1 -x * = g(x i ) -x * = G(x i -x * ),
with x 0 = x0 . We now apply our acceleration algorithm on the sequences x i and xi and compare the results. We first form the residuals,

r i = g(x i ) -x i = x i+1 -x i and ri = g(x i ) -xi = xi+1 -xi .
Consider the matrices of residuals R = [r 0 , . . . , r N ] and R = [r 0 , . . . , rN ]. We write P the perturbation matrix defined as

P RT R -R T R. (3.17)
The next proposition describes the sensitivity of Algorithm 1 using R and P . Then the norm of ∆c * is bounded by

∆c * ≤ P ( RT R) -1 c * . (3.19)
Proof. We start with the sequence xi . Let μ * be the dual variable of the equality constraint of (3.14). Both c * = c * + ∆c * and μ * = µ * + ∆µ * should satisfy the KKT system

2 RT R 1 1 T 0 c * μ * = 0 1 ⇔ 2(R T R + P ) 1 1 T 0 c * + ∆c * µ * + ∆µ * = 0 1 .
Indeed, using the definition of c * and µ * in (3.16),

2(R T R + P ) 1 1 T 0 c * + ∆c * µ * + ∆µ * = 2R T R 1 1 T 0 c * µ * + 2R T R 1 1 T 0 ∆c * ∆µ * + 2P 0 0 0 c * + ∆c * µ * + ∆µ * , = 0 1 + 2R T R 1 1 T 0 ∆c * ∆µ * + 2P 0 0 0 c * + ∆c * µ * + ∆µ * .
With this simplification, the system becomes

2R T R 1 1 T 0 ∆c * ∆µ * + 2P 0 0 0 c * + ∆c * µ * + ∆µ * = 0 0 . It remains to isolate c * , 2(R T R + P ) 1 1 T 0 ∆c * ∆µ * = 2P c * 0 .
The explicit solution in obtained by inverting the block matrix, and is written

∆c * = I - (R T R + P ) -1 11 T 1 T (R T R + P ) -1 1 (R T R + P ) -1 P c * .
We can bound the norm of ∆c * by

∆c * = I - ( RT R) -1 11 T 1 T ( RT R) -1 1 ( RT R) -1 P c * .
Because the first factor is the norm of a projector of rank N , its value is bounded by 1, so we get the desired result.

This proposition bounds the relative error on c * in comparison with c * . We will see that the perturbation magnitude can be arbitrarily large, which is the key issue with the convergence results in [73, §7]. Even when P is small, the term (R T R + P ) -1 is problematic. Our problem here is the structure of the residuals matrix R, R = [r 0 , Gr 0 , G 2 r 0 , ..., G N r 0 ],

which matches exactly the structure of Krylov matrices, i.e. square matrices K formed using a matrix M and a vector v, and computed as

K = [v, M v, M 2 v, ..., M N v].
Tyrtyshnikov [START_REF] Tyrtyshnikov | How bad are Hankel matrices[END_REF] showed that the condition number of Krylov matrices is lower bounded by a function which grows exponentially with N . Now, the error bound (3.19) contains the norm of the inverse of a perturbed squared Krylov matrix, which makes the situation even worse. In other words, even if the perturbations are small, their impact on the solution can be arbitrarily large. Even in practical cases where N is small (for example, N = 5), RT R is usually a singular or nearlysingular matrix. This particular issue means the linear system (R T R) -1 1 in (3.15) needs to be regularized.

Regularized Nonlinear Acceleration

In this section, we will analyze the following acceleration algorithm, which uses Tikhonov regularization to solve the linear system in (3.15).

c T 1=1 Rc 2 + λ c 2 ,
or equivalently solve

( RT R + λI)z = 1 then set c * λ = z/1 T z. Ensure: Approximation of x * computed as N i=0 (c * λ ) i xi , or N +1 i=1 (c * λ ) i xi
Regularization controls the norm of the coefficients produced by the algorithm and reduces the impact of perturbations, as shown in the following proposition.

c * λ = ( RT R + λI) -1 1 1 T ( RT R + λI) -1 1 , (3.22) ∆c * λ = I - ( RT R + λI) -1 11 T 1 T ( RT R + λI) -1 1 ( RT R + λI) -1 P c * λ . (3.23)
We begin by the bound on c * λ . Indeed, with (3.22),

c * λ 2 = 1 T ( RT R + λI) -2 1 (1 T ( RT R + λI) -1 1) 2 , ≤ 1 N + 1 max v =1 v T ( RT R + λI) -2 v (v T ( RT R + λI) -1 v) 2 , = 1 N + 1 max v =1 ( RT R + λI) -1 2 ( RT R + λI) -1 2 v 2 ( RT R + λI) -1 2 v 4 , ≤ 1 N + 1 ( RT R + λI) -1 2 2 max v =1 1 ( RT R + λI) -1 2 v 2 , = 1 N + 1 ( RT R + λI) -1 2 2 ( RT R + λI) 1 2 2 .
The norm of the coefficients c * λ are thus bounded by

c * λ ≤ 1 N + 1 RT R + λ λ = R 2 + λ (N + 1)λ .
We will now bound ∆c λ . With equation (3.23),

∆c λ = I - ( RT R + λI) -1 11 T 1 T ( RT R + λI) -1 1 ( RT R + λI) -1 P c * λ , ≤ I - ( RT R + λI) -1 11 T 1 T ( RT R + λI) -1 1 ( RT R + λI) -1 P c * λ , ≤ ( RT R + λI) -1 P c * λ ,
where the last inequality is obtained by bounding the norm of a projector. Since RT R 0, we have ( RT R + λI) λI, we get

∆c λ ≤ P λ c * λ ,
which is the desired result.

Regularization allows a better control of the impact of perturbations, but also changes the solution c * into c * λ . The next section analyses the impact of regularization on the extrapolated solution when there are no perturbations.

Constrained Chebyshev Polynomial

The previous section shows that regularization is important to control the perturbations present in (Pert. LFPI). However, convergence analysis becomes more complicated in the perturbationfree case, and we introduce constrained Chebyshev polynomials. Definition 3.3.3. The constrained Chebyshev polynomial C * κ (x, N, γ) of degree N , range (1κ) and norm constraint γ > 1 is defined as the solution of

C * κ (x, N, λ) = argmin C∈R N [x] : C(1)=1 max x∈[0,1-κ] C(x), s.t. C ≤ γ √ 1 + N
where C corresponds to the 2 norm of the coefficients of polynomial C. We write the maximum value as

S κ (N, γ) max x∈[0,1-κ] C * κ (x, N, κ). (3.24)
Unfortunately, as far as we know, there exists no explicit bound for S κ (N, γ). However, it is still possible to compute this value using sum-of-squares.

Constrained Chebyshev Polynomials and SOS

We briefly recall basic results on Sum of Squares (SOS) polynomials and moment problems [START_REF] Nesterov | Squared functional systems and optimization problems[END_REF][START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Parrilo | Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization[END_REF], which will allow us to formulate problem (3.24) as a (tractable) semidefinite program.

A univariate polynomial is positive if and only if it is a sum of squares. Furthermore, if we let m(x) = (1, x, . . . , x N ) T we have, for any q(x) ∈ R [2N ] , the space of real polynomials of degree at most 2N , q(x) ≥ 0, for all x ∈ R q(x) = m(x) T Cm(x), for some C 0, which means that checking if a polynomial is non-negative on the real line is equivalent to solving a linear matrix inequality (see e.g. [8, §4.2] for details). We can thus write the problem of computing the maximum of a polynomial over the real line as minimize t

subject to t -p(x) = m(x) T Cm(x), for all x ∈ R C 0, (3.25)
which is a semidefinite program in the variables p ∈ R N +1 , C ∈ S N +1 (the space of real positivesemidefinite matrices) and t ∈ R, because the first contraint is equivalent to a set of linear equality constraints. Then, showing that p(x) ≥ 0 on the segment [0, 1κ] is equivalent to showing that the rational fraction

p (1 -κ)x 2 1 + x 2
is non-negative on the real line, or equivalently, that the polynomial

(1 + x 2 ) N p (1 -κ)x 2 1 + x 2
is non-negative on the real line. Overall, this implies that problem (3.24) can be written

S κ (N, γ) = min t s.t. t -(1 + x 2 ) N +1 1 -(1-κ)x 2 1+x 2 q (1-κ)x 2 1+x 2 = m(x) T Cm(x), for all x ∈ R 1 T q = 1, q 2 ≤ γ, C 0, (3.26)
which is a semidefinite program in the variables q ∈ R N +1 , C ∈ S N +2 and t ∈ R.

Constrained Chebyshev Polynomials and Acceleration

Using S κ (N, γ) we can now bound the accuracy of the extrapolated point using the regularized algorithm. For simplicity, we consider here the constrained version of Algorithm 2, where the coefficients are computed using the norm constraint

c * γ = argmin c:1 T c=1 Rc s.t. c ≤ γ √ 1 + N .
Equivalently, c * γ is computed using R instead of R. This assumption is not restrictive, because we control the norm of the coefficients c * λ (Proposition 3.

c * λ ≤ 1 √ N + 1 1 + R 2 λ Since we want c * λ ≤ γ √ 1+N
, in the worst case,

1 √ N + 1 1 + R 2 λ ≤ γ √ N + 1 ⇔ λ ≥ R 2 (γ 2 -1)
.

The previous proposition also works when we apply the acceleration algorithm on (LFPI).

We can now show the link between constrained Chebyshev polynomials and the accuracy of the extrapolation in the perturbation-free case.

Proposition 3.3.5. Assume we have a sequence of points x i generated by (LFPI). Assume we use Algorithm 2, where the norm of the coefficients is bounded by γ √ 1+N . Then, the residual norm of the extrapolation in the perturbation-free case is bounded by

g   N i=0 (c * γ ) i x i   - N i=0 (c * γ ) i x i ≤ r 0 S κ (N, γ) . (3.28)
Proof. In equation (3.8), we have seen that if the coefficients c i sum to one then the residual of the extrapolation is equal to the combination of residuals, so

g N i=0 (c * γ ) i x i - N i=0 (c * γ ) i x i = N i=0 (c * γ ) i r i . By definition of c * γ , N i=0 (c * γ ) i r i = min c:1 T c=1 Rc s.t. c ≤ γ √ 1 + N .
Using the Krylov structure of the matrix R, for any c we have

Rc ≤ r 0 N i=0 (c) i G i = r 0 N i=0 p(G) .
where p is a polynomial of degree N with coefficient c. We thus have as minimization problem,

N i=0 (c * γ ) i r i ≤ min p∈R [N ] :p(1)=1 r 0 N i=0 p(G) s.t. p ≤ γ √ 1 + N .
Since 0 G 1-κ, after maximization over all matrices G with spectrum contained in [0, (1-κ)],

N i=0 (c * γ ) i r i ≤ r 0 min p∈R [N ] :p(1)=1 max x∈[0,(1-κ)] N i=0 p(x) s.t. p ≤ γ √ 1 + N .
The right-hand-side is exactly the definition of the constrained Chebyshev polynomial, thus

N i=0 (c * γ ) i r i ≤ r 0 S κ (N, γ).
This result will be useful in Proposition 3.3.6 for bounding the norm g(x extr )x extr .

Convergence rate

We will now prove global accuracy bounds. Let the extrapolation be

x extr = N i=0 (c * γ ) i xi .
where c * γ is the solution of min

c:1 T c=1 Rc , c ≤ γ √ 1 + N . (3.29)
We recall that this way of computing the coefficients is not that different from Algorithm 2. In fact, it suffices to set some λ, compute c * λ using Algorithm 2 then set γ = c * λ √ 1 + N .

We will link the norm of the residual of the extrapolation,

g(x extr ) -x extr ,
with the optimal value of the optimization problem (3.29). We base the proof the the following decomposition of the residual norm,

g(x extr ) -x extr ≤ g(x extr ) -x extr + e extr , = g(x extr ) + N i=0 (c * γ ) i e i - N i=0 (c * γ ) i e i -x extr + e extr , = N i=0 (c * γ ) i xi+1 - N i=0 (c * γ ) i e i -x extr + e extr , ≤ N i=0 (c * γ ) i ri Acceleration + N i=0 (c * γ ) i e i Stability + e extr ,
Local Error (3.30) where e i is defined in equation (Pert. LFPI) and e extr = g(x extr )g(x extr ). The bound will be local because we have no control on the perturbation e extr . The next theorem bounds separately the two first terms using constrained Chebyshev polynomials.

Proposition 3.3.6. Let x extr be computed using (2), and the matrix E = [e 0 , e 1 , . . .] be the matrix of noises. Then, the residual of the extrapolation x extr computed by Algorithm 2 satisfies

g(x extr ) -x extr ≤ r 0 S κ (N, γ) + γ √ 1 + N R -R + E + e extr . ( 3

.31)

Proof. We first bound the term Stability. Since we have the constraint

c γ ≤ γ √ 1+N , N i=0 (c * γ ) i e i ≤ E γ √ 1 + N .
Now we bound the term Acceleration. By definition of c * γ ,

Rc * γ = min c:1 T c=1 Rc s.t. c ≤ γ √ 1 + N .
In addition,

Rc ≤ Rc + ( R -R)c .
Again, due to the constraint on c ,

Rc * γ ≤ γ R -R √ 1 + N + min c:1 T c=1 Rc s.t. c ≤ γ √ 1 + N .
The right-hand-side is bounded in Proposition 3.3.5, thus

Rc * γ ≤ γ R -R √ 1 + N + r 0 S κ (N, γ).
Finally, putting all together,

g(x extr ) -x extr ≤ r 0 S κ (N, γ) + γ √ 1 + N R -R + E + e extr .
In these bound, we clearly see that the constraint on the norm limits the impact of perturbations. The next proposition can get rid of the term e extr , but at a cost of a factor 1 κ . The idea consists in bounding x 0x * rather than the residual. This bound better controls the impact of perturbation, because they all are multiplied by γ. Proposition 3.3.7. The distance of the extrapolated point x extr to the limit x * is bounded by

x extr -x * ≤ 1 κ x 0 -x * S κ (N, γ) + γ √ 1 + N R -R + E . (3.32)
Proof. Indeed,

x extr -x * = (G -I) -1 (G -I)(x extr -x * ) ≤ (G -I) -1 g(x extr ) -x extr
Since -I G -I -κI, we have

x extr -x * ≤ 1 κ g(x extr ) -x extr .
The term g(x extr )x extr corresponds to the sum of Acceleration and Stability in equation (3.30), bounded in Theorem 3.3.6. Finally,

r 0 = (G -I)(x 0 -x * ) ≤ x 0 -x * .
Because -I G -I -κI, thus G -I ≤ 1. Combining the two bounds together leads to the desired result.

We will now link R -R with E, if we assume all the perturbation are bounded by a maximum perturbation .

Proposition 3.3.8. Assume e i ≤ ∀i = 0 . . . N . Then,

1 √ 1 + N R -R + E ≤ 1 + 2 κ
Proof. We begin by the bound on E . Since e i ≤ ,

E ≤ N i=0 e i 2 ≤ √ N + 1.
The error R -R can be computed column-wise, where

R -R ≤ N i=1 ri -r i ≤ 2 N i=1 xi -x i However, xi -x i = g(x i-1 ) + e i-1 -g(x i-1 ), = g(x i-1 -x i-1 ) + e i-1 , = e i-1 + Ge i-2 + G 2 e i-3 . . . = i j=1 G i-j e j-1 ,
whose norm is bounded by

xi -x i ≤ i j=1 (1 -κ) i-j e j-1 ≤ i-1 k=0 (1 -κ) k = 1 -(1 -κ) i-1 κ .
To simplify further results, we will use the coarse bound

1 -(1 -κ) i-1 κ ≤ κ .
This bounds R -R by

2 √ 1 + N κ
The sum of the two bounds leads to the desired result.

We now have the final expression of the accuracy bound for the RNA algorithm.

Theorem 3.3.9. Assume we apply the RNA algorithm on a sequence of vector generated by (Pert. LFPI), where G ≤ 1κ and the norm of the perturbation e i is bounded by . Let the norm of the coefficients c * λ computed by Algorithm 2 be bounded by γ √ N +1 with γ ≥ 1. Then, the accuracy x extrx * is bounded by

x extr -x * ≤ x 0 -x * 1 κ S κ (N, γ) Accelerated rate + γ x 0 -x * 1 + 1 κ Impact of perturbation . (3.33)
The next proposition gives us a condition on γ to recover an optimal rate of convergence when the perturbation becomes smaller. , where γ satisfies

γ = γ( ), γ ≥ 1, lim →0 γ √ 1 + N = ∞ and lim →0 γ x 0 -x * = 0,
then we asymptotically recover the optimal rate (3.9) with a factor 1 κ , i.e.,

x extr -x * ≤ x 0 -x * κ 2β N 1 + β 2N .
Proof. The proof is straightforward using Theorem 3.2.2. It suffices to remark that S κ (N, ∞) corresponds to an unconstrained Chebychev polynomial.

In short, this bounds tells us that once the perturbations become small, the norm of the coefficients should be large, but not too large in comparison with the ratio between the perturbation and the distance to the optimum x 0x * . This is an expected result, since larger coefficients lead to a better acceleration term, while keeping the norm not too large avoid stability issues.

Regularized Acceleration for Stochastic Algorithms

If we assume the perturbations e i are sampled from some probability distribution whose first and second moments are bounded, then we also recover an accelerated rate of convergence in expectation.

Proposition 3.3.11. Assume we generate a sequence {x i } from (Pert. LFPI). Let e i be a random noise where

E[e i ] ≤ ν, E[e i -E[e i ]] 2 ≤ σ 2 . Assume c * λ ≤ γ √
1+N . In this case, the bound of Theorem 3.3.9 holds in expectation with = ν + σ.

Proof. Indeed, from Theorem 3.3.9,

E x extr -x * ≤ x 0 -x * 1 κ   S κ (N, γ) + E γ x 0 -x * 1 + 1 κ   .
In fact, there is no source of noise in S κ (N, γ), since it corresponds to the rate of convergence in the noiseless case. In the right hand side, the only source of randomness is , which bounds xix i ≤ e i ≤ ∀i.

Let ν i = E[e i ].
If we take the expectation,

E[ e i ] = E e i -ν i + ν i ≤ E e i -ν i + ν i .
By assumption, the right hand side is bounded by ν. The left hand side becomes

E e i -ν i = E e i -ν i 2 ≤ E e i -ν i 2 ,
which, by assumption, is bounded by σ.

In short, this proposition ensures acceleration even in the presence of random perturbation, if the ratio ν + σ x 0x * is small. In practice, if we have a noise independent of the iterates, then this is not always true, especially when we get closer to the optimum x * . Intuitively, this means the magnitude of the noise becomes too large in comparison with the residues. In other words, the "signal to noise" ratio becomes worse as we approach the limit x * . We will see that this will not be the case with variance-reduced algorithms.

Computational Complexity of Convergence Acceleration

In Algorithm 2, computing the coefficients c * λ means solving the N × N system ( RT R + λI)z = 1. We then get c * λ = z/(1 T z). This can be done in both batch and online mode. We will see that, in any case, we end up with a complexity of O(dN 2 + N 3 ), for a small value of N (usually,

N ∈ [START_REF] Ascher | Implicit-explicit methods for time-dependent partial differential equations[END_REF][START_REF] Drori | Performance of first-order methods for smooth convex minimization: a novel approach[END_REF]). The complexity of the acceleration algorithm is linear in the dimension, thus adding a negligible additional computation cost to the original procedure.

Online updates

Here, we receive the vectors ri one by one from the optimization algorithm, and we would like to solve the linear system in parallel of the optimization algorithm. In this case, we perform low-rank updates on the Cholesky factorization of the system matrix [START_REF] Golub | Matrix computations[END_REF]. At iteration i, we have the Cholesky factorization LL T = RT R + λI, where L is a triangular matrix. We receive a new vector r + and we want

L + L T + = L 0 a T b L T a 0 b = RT R + λI RT r + ( RT r + ) T r T + r + + λ .
We can explicitly solve this system in variables a and b, and the solutions are

a = L -1 RT r + , b = a T a + λ.
The complexity of this update is thus O(i d + i 2 ), i.e. the matrix-vector multiplication of RT r + with cost O(i d) and solving a i × i triangular system with cost O(i 2 ). Since we need to do it N times, the final complexity is thus O(dN 2 + N 3 ).

Batch mode

The complexity is divided in two parts: First, we need to build the linear system itself. Since R ∈ R d×N , it takes O(dN 2 ) flops to perform the multiplication. Then we need to solve the linear system ( RT R + λI)z = 1 which can be done by Gaussian elimination (in particular when N is small), by Cholesky factorization or by using an iterative method like conjugate gradient. It takes O(N 3 ) flops to solve the linear system in the worst case, meaning that the overall complexity is O(dN 2 + N 3 ).

Online Regularized Nonlinear Acceleration and Adaptive Momentum-based Algorithms

So far, we analyzed an algorithm which extrapolates the limit point x * using the iterates from (Pert. LFPI). However, these results are limited to single-step algorithms such as gradient descent (due to the assumption that the matrix G in (Pert. LFPI) has to be symmetric), thus excluding much faster momentum-based methods such as Nesterov's fast gradient algorithm.

Accelerating a Class of Non-symmetric Algorithms

Here, we will add an additional structural assumption on the dynamic of the points x i , allowing (variable) linear combinations between two calls of the function g,

   x i = G(y i-1 -x * ) + x * (⇔ Ax + b), y i = i j=1 α (i) j x j + i-1 j=0 β (i) j y j , (Comb. LFPI)
and the perturbed version,

   xi = G(ỹ i-1 -x * ) + x * + e i , (⇔ xi = g(y i-1 ) + i ), ỹi = i j=1 α (i) j xj + i-1 j=0 β (i) j ỹj , (Comb. Pert. LFPI)
where the coefficient α and β are arbitrary. We will see in the next chapter that this class of algorithms includes, for example, momentum-based optimization methods. We wrote x * explicitly in the algorithm to simplify proofs and notations, but in practice this quantity is of course not known, and does not need to be known when using the acceleration method.

It is not possible to apply RNA directly on this sequence since it cannot be written as (Pert. LFPI). This is mainly due to the dependence of x i to the previous x i-1 , x i-2 , . . . , x 1 .

Even if we augment artificially the states of the problem, for example with the concatenation

z i = x i ; x i-1 ; . . . ; y i-1 ; . . .

T

then the underlying matrix G describing the algorithm is not symmetric, so our convergence analysis does not apply. In addition, the coefficients α

c T 1=1 Rc 2 + λ c 2 ,
or equivalently solve

( RT R + λI)z = 1 then set c * λ = z/1 T z. Ensure: Approximation of x * computed as N +1 i=1 (c * λ ) i xi
The next Theorem shows that when there are no perturbations, under mild assumptions on the coefficients we recover the optimal rate in Theorem 3.2.2 if we do not use regularization.

Theorem 3.5.1. Let the sequences of points {y 0 . . . y N -1 } and {x 1 . . . x N } be generated by

(Comb. LFPI) with 0 G (1 -κ)I. Let i j=1 α (i) j + i-1 j=0 β (i) j = 1, α (i) i = 0.
If we apply Algorithm 3, the rate of convergence is bounded by

g(x extr ) -x extr ≤ 2β N 1 + β 2N r 0 where β = 1 - √ κ 1 + √ κ < 1. (3.34)
Proof. This proof extend the one of Theorem 3.2.2. We have to show that minimizing the residues

min c:1 T c=1 N i=0 c i r i
is equivalent to finding the minimal polynomial

min p∈R [N ] :p(1)=1 p(G)r 0 .
It is equivalent to prove that r i can be written as a polynomial in G of degree exactly i, whose coefficients sums to one. Indeed, assume we have

r i = p i (G)r 0 , p i (1) = 1.
In this case, the family of polynomial {p i } i=0...N generates the space of polynomials of degree N . We can thus create p(G)r 0 for any p from a linear combination of r i . In addition, if the coefficients of this linear combination sum to one, we create a polynomial p(G) whose coefficients sums to one,

N i=0 c i p i (1) = N i=0 c i = 1.
Now, we have proved that

r i = p i (G)r 0 , p(1) = 1 ⇒ min c:1 T c=1 N i=0 c i r i = min p∈R N [x]:p(1)=1 p(G)r 0 .
The right hand side is already bounded by Theorem 3.2.2. It remains to show the equality, which is true if r i = p i (G)r 0 , where p(1) = 1 and deg(p i ) = i. For simpler notations, we remove the superscript (i) for coefficients α (i) j and β (i) j without loss of generality. We will prove the claim by recursion. We start with the simple case i = 0, i.e., r 0 = Ir 0 and the identity matrix I indeed corresponds to the polynomial p(G) = 1 • G 0 of degree 0, whose coefficients sums to one. Now, assume the claim true for r i-1 , and we will show it is also true for r i . Indeed,

r i = x i+1 -y i , = G(y i -x * ) + x * -y i , (3.35)
= (G -I)

  i j=1 α j x j + i-1 j=0 β j y j -x *   . Since i j=1 α j + i-1 j=0 β j = 1, (G -I)   i j=1 α j x j + i-1 j=0 β j y j -x *   = (G -I)   i j=1 α j (x j -x * ) + i-1 j=0 β j (y j -x * )   .
In addition,

x j -x * = G(y j-1 -x * ) + x * -x * = G(y j-1 -x * ).
Finally, using equation (3.35),

(G -I)(y j -x * ) = r j .
Bringing it all together leads to

(G -I)   i j=1 α j x j + i-1 j=0 β j y j -x *   , = (G -I)   i j=1 α j (x j -x * ) + i-1 j=0 β j (y j -x * )   , = (G -I)   i-1 j=0 (α j G + β j I)(y j -x * )   , = i-1 j=0 (α j+1 G + β j I)(G -I)(y j -x * ), = i-1 j=0 (α j+1 G + β j I)r j .
We assumed by recursion that r j = p j (G)r 0 , where p j is a polynomial of degree exactly j, where p j (1) = 1. The previous line is thus written

r i = i-1 j=0 α j+1 Gp j (G) + β j p j (G) =pi(G) r 0 .
Let us check if p i (1) = 1. Indeed,

p i (1) = i-1 j=0 α j+1 p j (1) + β j p j (1) = i-1 j=0 α j+1 p j (1) + β j p j (1) = 1.
Now, we check if the degree of the polynomial is indeed equal to i. We have

deg(p i (x)) = deg   i-1 j=0 (α j x + β j )p j (x)   .
Since by assumption deg(p

j (x)) < deg(p i-1 (x)) for j < i -1, deg(p i (x)) = deg (α i x + β i-1 )p i-1 (x) .
Because α i = 0 by assumption,

deg(p i (x)) = deg α i xp i-1 (x) = 1 + deg(p i-1 (x)).
Finally, because the degree of p i-1 is exactly equal to i -1,

deg(p i ) = 1.
We proved by recursion that deg(p i ) = i and p i (1) = 1, which concludes the proof.

The previous result also shows that our method is adaptive. Whatever the algorithm, if the iteration g(x) converges, the coefficients α i and β i sum to one and we have α N non-zero, we obtain an optimal rate of convergence. In comparison, in many case coefficients α and β are function of κ and aim to improve the rate of convergence of x ix * . In practice, it may be difficult to estimate κ, and thus the coefficients may not be optimal and this lead to sub-optimal rates.

We can now reproduce the bound of Theorem 3.3.9 using similar tricks.

x extr -x * ≤ x 0 -x * 1 κ S κ (N, γ) + Constant • γ x 0 -x * .
(3.36)

Proof. Indeed,

x extr = N i=0 (c * λ ) i xi+1 = N i=0 (c * λ ) i g(ỹ i ) = N i=0 (c * λ ) i g(ỹ i ) + e i
Thus, the distance to the optimum is bounded by

x extr -x * ≤ N i=0 (c * λ ) i g(ỹ i ) -x * + N i=0 (c * λ ) i e i . Because g(y) -x * = G(y -x * ), N i=0 (c * λ ) i g(ỹ i ) -x * + N i=0 (c * λ ) i e i = G N i=0 (c * λ ) i (ỹ i -x * ) + N i=0 (c * λ ) i e i .
In addition, because (G -I)(y

-x * ) = g(y) -y, G N i=0 (c * λ ) i ỹi ) -x * + N i=0 (c * λ ) i e i ≤ (G-I) -1 G N i=0 (c * λ ) i g(ỹ i ) -ỹi + (G -I) N i=0 (c * λ ) i e i .
Finally, because g(ỹ i ) + e i = ri , (G -

I) 1-≤ 1/κ and G ≤ 1 -κ, (G-I) -1 G N i=0 (c * λ ) i g(ỹ i ) -ỹi + (G -I) N i=0 (c * λ ) i e i ≤ 1 -κ κ    N i=0 (c * λ ) i ri + N i=0 (c * λ ) i e i    .
Here, we recover exactly the structure in equation (3.30) (without the term e extr ), which was the basis of the proof in Theorem 3.3.9.

Theses results allows us to introduce the online-RNA algorithm, which consists in using RNA at every iterations.

Online RNA

In the RNA algorithm, the output corresponds to a linear combination of the points x i , whose coefficients sums to one. It corresponds, in fact, to what we described in equation (Comb. Pert. LFPI)

where, at iteration N ,

α (N ) = c * λ , β (N ) = 0.
This describes the Online RNA algorithm.

Algorithm 4 Online Regularized Nonlinear Acceleration

Require: Function g(x) with fixed-point x * , starting point ỹ0 , regularization λ.

1: i = 1.
2: while the norm of residual g(ỹ i-1 )ỹi-1 ≥ Tolerance do 3:

Compute the next iterate xi = g(ỹ i-1 ).

4:

Compute the extrapolation ỹi = RNA({x j }, {ỹ j }, λ) (See Algorithm 3).

5:

i = i + 1.

6: end while

Since this algorithm matches the description of (Comb. Pert. LFPI), it satisfies the assumptions of Theorem 3.5.1, and thus presents an optimal rate of convergence on linear fixed-points iterations. Contrary to many acceleration techniques, this can be done without any knowledge on the problem's constants.

Since the complexity grows quadratically with N (the number of points in the sequences {x i } and {y i }), Algorithm 4 is not convenient because because N continuously increases over time.

In practice, it is better to use a fixed window where N is bounded. This strategy is similar to the limited version of BFGS, where a fixed number of gradients is stored in memory [START_REF] Boyd | Convex optimization[END_REF].

Adaptive Momentum-based Algorithms

Previously, we showed two different aspects of the acceleration of scheme (Comb. Pert. LFPI).

One is capable to extrapolates x extr , an approximation of the solution x * , from sequences {x i } and {ỹ i }. The other takes as input an iterative method g and accelerate it using RNA at each step. Ideally, we would like to combine the two approaches: using RNA at each iteration on algorithm (Comb. Pert. LFPI). In this section, we will try to unify these two strategies.

Assume that scheme (Comb. Pert. LFPI) has a proven convergence rate, which bounds the error at each iteration. This rate at iteration N is ensured if the condition h N (x N , ỹN-1 ) ≤ 0 is satisfied, for some function h N . For example, in the context of optimization, this may correspond to a sufficient descent condition for the point xN and g(x) is a simple gradient step with a certain step length. Assume this condition is automatically satisfied if

xN = g(ỹ N -1 ).
Again, in optimization, this is true when the step length is equal to 1/L. We will describe a way to combine the RNA algorithm without losing the theoretical guarantee of the original algorithm.

For simplicity, we consider β i = 0 but it is possible to extend to the case where β is non-zero.

Ideally, we would like to use acceleration as much as we can, but without losing the potential convergence guarantees of the algorithm. We will see in the next proposition that is some cases, it is possible to test a condition which decides if we should use acceleration or not. This test is based on the application of the function h to a "conditional point" z i . Proposition 3.5.3. Assume a convergence rate is proved for sequences {x i } and {ỹ i } if

h N (x N , ỹN-1 ) ≤ 0, y N = N i=1 α i xi , α N = 0, N i=1 α i = 1
where the first condition is ensured when xN = g(ỹ N -1 ) (we omitted the superscript (N ) for clarity). Let c * λ be the coefficients computed by RNA with sequences {x i } and {ỹ i }. Consider the "conditional point"

zN = N j=1 ρ i xi where ρ N = (c * λ ) N α N , ρ j<N = (c * λ ) j -α j α N .
Then, the rate of convergence is preserved if

h N (z N , ỹN-1 ) ≤ 0, y N = N j=1 (c * λ ) j xj .
This means we can use acceleration if h N (z N , ỹN-1 ) ≤ 0. Otherwise, we should perform a step of the original algorithm.

Proof. Since the convergence rate is proved for {x i } and {ỹ i }, it is also proved for any z,

ỹN-1 s.t. h(z, ỹN-1 ) ≤ 0, y N = α N z + N -1 i=1 α i xi .
In particular, this is true for z = zN . In this case, ỹN becomes

ỹN = α N zN + N -1 i=1 α i xi .
If we replace zN by its expression,

ỹN = α N N i=1 ρ i xi + N -1 i=1 α i xi , = (c * λ ) N xN N -1 i=1 (α i + (c * λ ) j -α j )x i , = N i=1 (c * λ ) i xi .
This theorem shows that the conditional point zN is crucial the design of an adaptive momentum-based algorithm, and determines if we have to use the adaptive coefficient c * λ or the user defined coefficients α and β.

Since the coefficients c *

λ are adaptive, we can usually expect a better rate of convergence when using the RNA algorithm. However, we have seen in Theorem 3.3.9 that the upper bound for the accuracy is local and may be very bad in some situation.

Combining the two situations, i.e., the adaptivity of RNA and the convergence guarantees of algorithm (Comb. Pert. LFPI) takes the best of the two worlds, at the cost of one evaluation of h(•, ỹN-1 ). We formalize the Online RNA applied to momentum-based methods in Algorithm 5.

Algorithm 5 Online Regularized Nonlinear Acceleration for Momentum-Based Methods

Require: Algorithm (Comb. Pert. LFPI), condition function h(•, •), starting point ỹ0 , regularization λ.

1: i = 1.
2: while the norm of residual g(ỹ i-1 )ỹi-1 ≥ Tolerance do

3:

Compute the next iterate xi = g(ỹ i-1 ).

4:

Compute the extrapolation x extr ; c * λ = RNA {x j }, {ỹ j }, λ .

5:

Compute the conditional point zi = (c * λ

) i xN + i-1 j=0 (c * λ ) j -α j x (i) j 1/α (i) i . 6:
if the condition h(z i , ỹi-1 ) ≤ 0 is satisfied then

7:
Use the extrapolation ỹi = x extr .

i = i + 1.
12: end while This algorithm, capable to combine robustness, adaptive acceleration and theoretical guarantees, concludes this chapter on generic acceleration.

Conclusion and Perspectives

In this chapter, we introduced the Regularized Nonlinear Acceleration technique, which is able to accelerate the perturbed fixed-point operation with linear combination described in equation (Comb. Pert. LFPI), if the underlying matrix G is symmetric. We deduced a rate of convergence, whose expression depends on the maximum value of a constrained Chebyshev Polynomial (3.24) and the ratio between the (maximum) noise level and the distance to the optimum x * . Its computational complexity is small in comparison with the computation of g(x), since it is linear in the dimension d when the length of the sequence to extrapolate is small. This chapter raises two main questions. We have seen that we can accelerate a fixed-point iteration when G is symmetric, but the case where G is not symmetric anymore is less clear, because we have to compute minimal polynomial on the complex plane instead of the real line (due to the presence of complex eigenvalues). In addition, working with non symmetric matrices makes the spectral radius and the norm two different quantities, which makes the analysis even more difficult.

The second question is about the constrained Chebyshev polynomial. Even if we have rates of convergence, computing a constrained Chebyshev polynomial is hard because it involves modeling a sum-of-squares problem, which is then lifted to an SDP formulation. This is known to be not scalable for large dimensions and solving such problem is computationally hard. In addition, we do not have access to an explicit formula, thus preventing further theoretical analysis. Working on a good approximation of constrained Chebyshev polynomial may lead to better understanding of nonlinear acceleration algorithms.

67

For example, the gradient method, written,

x i+1 = x i -h∇f (x i )
with step length h > 0 matches exactly the definition, where g(x) = x -h∇f (x), α

i = 1 and β i = 0. Here, we will study the impact of RNA on optimization methods for accelerating their convergence.

Since the publication of Nesterov's optimal first-order smooth convex minimization algorithm [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o (1/k2)[END_REF], significant efforts have been focused on either providing more interpretable views on current acceleration techniques, or on replicating these complexity gains using different, more intuitive schemes. Early efforts sought to directly extend the original acceleration result in [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate o (1/k2)[END_REF] to broader function classes [START_REF] Nemirovskii | Optimal methods of smooth convex minimization[END_REF], allow for generic metrics, line searches, produce simpler proofs [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] or adaptive accelerated algorithms [START_REF] Nesterov | Universal gradient methods for convex optimization problems[END_REF], etc. More recently however, several authors [START_REF] Drori | Performance of first-order methods for smooth convex minimization: a novel approach[END_REF][START_REF] Lessard | Analysis and design of optimization algorithms via integral quadratic constraints[END_REF] have started using classical results from control theory to obtain numerical bounds on convergence rates that match the optimal rates. Others have studied the second order ODEs obtained as the limit for small step sizes of classical accelerated schemes, to better understand their convergence [START_REF] Su | A differential equation for modeling nesterov's accelerated gradient method: Theory and insights[END_REF][START_REF] Wibisono | A variational perspective on accelerated methods in optimization[END_REF]. Finally, recent results have also shown how to wrap classical algorithms in an outer optimization loop, to accelerate convergence and reach optimal complexity bounds [START_REF] Lin | A universal catalyst for first-order optimization[END_REF] on certain structured problems.

The approach we propose here is completely different as our method is built on top of existing optimization algorithms. Classical algorithms typically retain only the last iterate or the average [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF] of iterates as their best estimate of the optimum, throwing away all the information contained in the converging sequence of iterates. As it is highly wasteful from a statistical perspective, extrapolation schemes estimate the optimum of an optimization problem using a weighted average of the last iterates produced by an algorithm, where the weights depend on the iterates. Overall, nonlinear acceleration has marginal computational complexity. On convex problems, the online version of RNA is competitive with L-BFGS in our experiments and is robust to misspecified strong convexity parameters.

In the case where we can only access a stochastic gradient, we have seen in the previous chapter that all results on acceleration still hold. Such setting is typical in machine learning, for example when f is a sum of convex functions, and we only have access to the gradient of one randomly selected function [START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF]. Stochastic optimization is typically challenging as classical algorithms are not convergent (for example, gradient descent or Nesterov's accelerated gradient).

Even the averaged version of stochastic gradient descent with constant step size does not converge to the solution but to another point whose proximity to the real minimizer depends of the step size [START_REF] Nedić | Convergence rate of incremental subgradient algorithms[END_REF][START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF].

In this chapter, we will explore optimization algorithms which can be written as

x i = g(x i-1 ) = G(y i-1 -x * ) + x * + e i , (4.1)
y i = i j=1 α (i) j x j + i j=0 β (i) j y j ,
where

0 G (1 -κ)I, 0 < κ < 1, i j=1 α (i) j + i j=0 β (i) j = 1, α (i) j = 0.
(We recall we dropped the tilde ˜notation). Then, we will study the impact of acceleration Algorithm 3 and 5 on these optimization methods. We begin with deterministic scheme, then move to stochastic one and finish with the optimization of neural networks.

Acceleration of Deterministic Algorithms

Introduction

In this section, we consider the minimization of a smooth and strongly convex function f with constants 0 < µ ≤ L, i.e.,

µ y -x ≤ ∇f (x) ≤ L y -x ∀ x, y ∈ R d
where • stands for the 2 norm. We also assume the function two times differentiable, so

µI ∇ 2 f (x) LI , ∀x ∈ R d .
In addition, we also assume the Hessian is Lipchitz-continuous with constant M , which means

∇ 2 f (y) -∇ 2 f (x) ≤ M y -x .
This last assumption is necessary to easily bound the error e i in (4.1). Before going further, we will show that we can reduce any convex function f to a strongly convex one, whose constant depends on the desired accuracy.

Reduction to Strongly Convex Functions

Using a simple regularization trick which we trace back at least to [START_REF] Hazan | [END_REF], we can extend our results

to smooth functions that are not strongly convex. Suppose we seek to solve

min x∈R n f (x)
in the variable x ∈ R n , where f (x) has a Lipschitz continuous gradient with parameter L with respect to the Euclidean norm, but is not strongly convex. Assume for simplicity that the initial iterate x 0 is close enough to the optimum so that D ≥ x kx * for any k ≥ 0. We can approximate the above problem by

min x∈R n f ε (x) f (x) + ε 2D 2 x 0 -x 2 2 (4.2)
in the variable x ∈ R n , where f ε (x) has a Lipschitz continuous gradient with parameter L + ε/D 2 with respect to the Euclidean norm, is strongly convex with parameter ε/D 2 with respect to the same norm. Furthermore f ε (x) is an ε approximation of f (x) near the optimum and we get

f (x k ) -f (x * ) = f ε (x k ) - ε 2D 2 x 0 -x k 2 2 -f ε (x * ) + ε 2D 2 x 0 -x * 2 2 , ≤ f ε (x k ) -f ε (x * ) + ε 2 , ≤ f ε (x k ) -f ε (x * ε ) + ε 2 ,
using the smoothness of f ε (x) and writing x * ε the optimum of problem (4.2). It suffices to optimize f ε up to ε/2 to find an ε-solution for the original problem. The linear convergence of gradient [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] algorithms guarantees

f ε (x k ) -f ε (x * ε ) = (L + ε)D 2 2 r k , r = 1 - 2ε LD 2 + 2ε
.

The number of iterations required to reach a target precision ε/2 is thus bounded by

k = O log((L + ε)D 2 /ε) log(1/r) .
By replacing the value of r, we have

log(1/r) ∼ 1 -r = LD 2 ε , while accelerated algorithms have r = 1 -ε/(LD 2 + ε) which yields log(1/r) ∼ LD 2 ε .
Up to a logarithmic constant, these upper bounds match the complexity of gradient and accelerated gradient methods. Overall, an algorithm for strongly convex function used with this regularization trick recovers an ε-approximated solution. This means we can always reduce a not strongly convex problem to strongly convex one. Of course, in practice one does not wants to use this king of "hack", but strong convexity, for now, is crucial in the analysis of generic acceleration methods. The extension to convex functions still remains an open problem.

Rate of Convergence

Local Rate of Convergence

We will study several classical optimization algorithm and show they can be written into the form (4.1), but before we will analyze the potential gain from the RNA algorithm. In the case of deterministic optimization, all the algorithms we will list have the same nonlinear function g for generating x i , which is the simple gradient step

g(x) = x - 1 L ∇f (x). (4.3)
If we linearize the algorithm around x * , we get

g(x) = x * + I - 1 L ∇ 2 f (x * ) (x -x * ) + O x -x * 2 .
We identify here the matrix

G = I -1 L ∇ 2 f (x * ) , whose eigenvalues lies in [0, 1 -κ] where κ = µ
L is the inverse of the condition number. We will now show that the error is bounded in function of x 0x * if the hessian is Lipchitz continuous. Proposition 4.2.1. (Nesterov and Polyak [START_REF] Nesterov | Cubic regularization of Newton method and its global performance[END_REF]) Let f be two times differentiable, whose Hessian is Lipchitz-continuous with constant M . Then,

g(x) -x * -I - 1 L ∇ 2 f (x * ) (x -x * ) = e(x) ≤ M 2L x -x * 2 . ( 4.4)
Proof. Indeed,

g(x) -x * -I - 1 L ∇ 2 f (x * ) (x -x * ) = x - 1 L ∇f (x) -x * -I - 1 L ∇ 2 f (x * ) (x -x * ) , = 1 L ∇f (x) - 1 L ∇ 2 f (x * )(x -x * ) , = 1 L ∇f (x) -∇f (x * ) -∇ 2 f (x * )(x -x * ) ,
where ∇f (x * ) = 0 since x * is a critical point. Because

∇f (x) -∇f (x * ) = 1 0 ∇ 2 f (x + τ (x * -x)) (x -x * ) dτ,
we can bound the norm by

1 L ∇f (x) -∇f (x * ) -∇ 2 f (x * )(x -x * ) ≤ 1 L 1 0 ∇ 2 f (x + τ (x * -x)) -∇f (x * ) (x -x * ) dτ , ≤ M L x -x * 1 0 τ dτ,
which is the desired result.

If we assume that all the points x i are going to be inside some region around x * , ∀i,

x i -x * ≤ D,
then the bound of Theorem 3.3.9 is written

x extr -x * ≤ D 1 κ S κ (N, γ) + γ 2M D L 1 + 1 κ . (4.5)
We can now derive an upper bound for γ, and thus β, to recover the asymptotic rate of convergence in Proposition 3. This means RNA recovers the rate of convergence of Proposition 3.3.10.

Proof. Indeed, for 0 < a < 1, we have

lim D→0 D -a = ∞ and lim D→0 D 1-a = 0.
We have to satisfy the following relationship between λ and γ (Proposition 3.3.4)

R 2 (γ 2 -1) ≤ λ to ensure c * λ < γ √ 1+N . Since R 2 = O(D 2
) and (γ 2 -1) ≥ (γ -1) 2 for γ > 1, then the following value of λ satisfies the conditions,

λ = O(D 2 ) (γ -1) 2 .
Replacing γ by 1 + O(D -a ) leads to the desired result.

This result is more permissive that in our associated publications [START_REF] Scieur | Regularized nonlinear acceleration[END_REF] and [START_REF] Scieur | Nonlinear acceleration of stochastic algorithms[END_REF], where the range of λ was D 2 < λ < D 8 3 , meanwhile we can now go up to D 4 .

Global Rate of Convergence

However, in the theoretical bounds, the interactions between the perturbation and the algorithm are not taken into account, because we only study it from a worst-case point of view. However, we can show that there exists λ for which we can ensure global convergence, if xix * is converging monotonically. In fact, it suffices to take λ = ∞ (or equivalently, γ = 1) to see that

c * ∞ = 1 N [1, 1, . . . , 1] T .
This corresponds to take the average of previous points. If the distance to the optimum converges at rate (1r) i , the rate of convergence for the averaged of N point becomes

(1/N ) N i=0 xi -x * ≤ (1/N ) N i=0 xi -x * ≤ (1/N ) N i=0 (1 -r) i x 0 -x * ≤ 1 N κ x 0 -x * .
Indeed, this slows down the rate of the original algorithm, but one can simply perform a comparison between the extrapolation step and the vanilla step and take the best between the two.

First Order Optimization Algorithms

We will now list several standard optimization algorithms and write them into the form of (Comb. Pert. LFPI). In this section, g(x) corresponds to a gradient step with step-size h = 1/L, as described in equation (4.3). This makes the expression of the residual very simple, since it corresponds to the gradient, i.e., r i = 1 L ∇f (x i ).

Gradient Method with Line-search

As we just described, the gradient method can be written RT W R ). The window size (maximum sequence length) is set to 10. We run the experiment on Madelon dataset, whole regularization is set such that the condition number κ = 10 -6 . We clearly see that regularization is crucial for the stability of the algorithm. (V) (vanilla) corresponds to the original methods, and P (parallel) corresponds to the application of RNA without interaction with the original algorithm.

x i = g(y i-1 ), y i = x i .
L(w) with 2 penalty τ is defined as

L(w) = m i=1 log(1 + exp(-ϑ i Θ T i w)) + τ 2 w 2 .
Its smoothness constant is bounded above by Θ 2 2 /4 + τ and its strong convexity parameter is bounded below by τ .

We will test the RNA algorithm on previous methods. We compare three strategies:

• (V) Vanilla: we run the original algorithm with no acceleration.

• (P) Offline version, in parallel: we run the RNA algorithm on the side (Algorithm 3), with no interaction on the original optimization method.

• (O) Online: we use the online version of RNA (Algorithm 5).

Comparison between RNA and Anderson Acceleration

Here, we compare the RNA algorithm with memory size N = 10 and regularization parameter λ = 10 -8 RT R (where R corresponds to the matrix of gradients ∇f (x i ) when g(x) is a gradient step) and the Anderson Acceleration, which is basically RNA without regularization. In Figure 4.1, we clearly see that the regularization in crucial for acceleration because it plays an important role in the algorithm stability.

Observations

In the previous experiments, we can observe that online acceleration (as well as BFGS) is always beneficial in comparison with the original optimization method, and globally, BFGS and RNA are competitive with each other, with a slight advantage to BFGS in the previous experiments.

However, it is complicated to compare the two since there is a line-search subroutine in BFGS to ensure its convergence, while in RNA there is only one comparison between function values. In addition, there is no (good) global convergence rate for -BFGS, while Algorithm (5) preserves the rate.

Surprisingly, accelerating gradient descent with RNA looks to be a better strategy than accelerating Nesterov's method. In our experiments, in all cases, RNA combined with gradient descent always performs better than RNA with Nesterov, but the difference is quite small. On the other hand, accelerating an accelerated method guarantees an optimal rate of convergence. So acceleration with good global convergence rate comes at the cost of (sometimes) worse numerical convergence.

The claim that RNA makes momentum methods adaptive to the strong convexity parameter is verified in practice. Indeed, if we compare the two variants of Nesterov's method without acceleration, we have a big gap in the convergence rate (both numerically and theoretically).

However, once they are accelerated with RNA, they present a similar rate of convergence.

Even if we have some extra computation effort, the difference in time between the original and accelerated version of an algorithm is negligible. It can be explained that we have a sum of function in our experiments, making the oracle call expensive. In comparison, RNA only needs to solve an N × N (where N = 10) system of linear equations. We can even speed up RNA using a recursive algorithm for solving this system, if needed.

In summary, online acceleration never hurts and often improves performance. Even from the theoretical point of view, in the worst case we multiply the number of call to f (x) by two for getting the same result. It presents a similar numerical performance as BFGS despite the absence of line-search, but in addition a rate of convergence is guaranteed. Using online acceleration on accelerated methods ensures optimal rates, but come at the cost of slightly worse numerical performance.

Acceleration of Stochastic Algorithms

Introduction

We focus on the problem

min x∈R d f (x) (4.6) Proposition 3.3.11 becomes E x extr -x * ≤ x 0 -x * 1 κ S κ (N, γ) + γσ x 0 -x * 1 + 1 κ .
The rate of convergence of acceleration thus depends on the ratio σ x0-x * , close to the concept of signal to noise ratio. In particular, if σ is dependent of x 0x * , a typical behavior when using SGD on quadratics, then acceleration does no converge when this ratio becomes too large. However, we can expect a better convergence for early iterations (which is the part where SGD generalizes the best [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]) because x 0x * will be large in comparison with σ.

On the other hand with variance-reduction techniques, such as SAGA and SVRG, σ asymptotically goes to zero. To reduce the ratio / x 0x * , it suffices to sufficiently wait between two "iterations" to reduce this ratio, i.e., considering x i+1 to be x i+1 = g(g(. . . g(x i ) . . .)), for a large number of calls of g.

We will make a coarse, non totally rigorous, analysis to give an intuition why waiting several iteration leads to smaller variance. For example, consider the case where the standard deviation of x i is proportional to the distance to the optimum, i.e.,

E[ x

i -E[x i ] ] = σ E[x i ] -x * .
In addition, assume the method has a certain rate of convergence 1r,

E[ x i+1 -x * ] ≤ (1 -r)E[ x i -x * ]
Then, if we wait s iterations,

E[ x s -x * ] ≤ (1 -r) s E[ x 0 -x * ]
which means that the variance of x s is bounded by

E[ x s -E[x s ] ] ≤ σE[ x i -x * ] ≤ (1 -r) s σ x 0 -x * ,
which decreases over time. The upper bound for the convergence rate becomes, in that case,

E x extr -x * ≤ x 0 -x * 1 κ S κ (N, γ) + γ(1 -r) s σ 1 + 1 κ .
With s sufficiently large, we can make the right hand side arbitrary small. Of course, we have to find the best trade-off between accelerating as much as we can, and waiting s non-accelerated iterations.

For simplicity, we consider in the rest of this section s = m (i.e., we wait one epoch), the number of data points in the data set, because it seems to be sufficiently large in our numerical experiments. Again, for simplicity, we consider x i+1 to be actually x i+1 = g m (x i ), which means that "one iteration" for RNA corresponds to m stochastic iterations.

Numerical Experiments (Stochastic)

We focus on the composite problem

min x∈R d f (x) = m i=1 1 m f i (x) + τ 2 x 2 ,
where f i are convex and L-smooth functions and τ is the regularization parameter. We will use classical methods for minimizing f (x) such as SGD (with fixed step size), SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF], SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF], and also the accelerated algorithm Katyusha [START_REF] Allen-Zhu | Katyusha: The first direct acceleration of stochastic gradient methods[END_REF]. We will compare their performances with and without the (potential) acceleration provided by Algorithm 3. We set λ = 10 -8 RT R and N = 10 for all the experiments.

Comparison Between Acceleration Strategies

In order to balance the complexity of the extrapolation algorithm and the optimization method we wait several data queries before adding the current point (the "snapshot") of the method to the sequence. Indeed, the extrapolation algorithm has a complexity of O(N 2 d), where N is small (10 in our experiments). If we wait at least O(d) updates, then the extrapolation method is of the same order of complexity as the optimization algorithm. A similar technique is used in [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF].

• SGD. We add the current point after m data queries (i.e. one epoch). The residual is considered to be x i+1x i where x i is the variable before the epoch and x i+1 the variable after the epoch.

• SVRG. We compute the gradient exactly, then perform m queries (the inner-loop of SVRG). The iterate x i corresponds to the point before computing the full gradient, and

x i+1 is the point computed after one iteration of the outer-loop. This means one iteration of SVRG is twice the cost of one iteration of SGD.

• SAGA. The main difference between SAGA and SGD is the presence of a table of computed gradients, but this does not matter when using RNA algorithm. We accelerate it exactly as SGD, and we do not touch the table of gradient between two accelerations.

• Katyusha. We compute the gradient exactly, then perform 4m gradient calls (the innerloop of Katyusha). We will not accelerate Katyusha with RNA, as its expression does not Because the iterates of SAGA and SVRG have low variance, their accelerated version converges faster to the optimum. Even if we can see that, in our experiments, SAGA performs usually better than SVRG, their online accelerated version is comparable.

When we accelerate SAGA and SVRG using the offline version of RNA, we obtain a method whose performance is better than Katyusha, but without the specification of the strong convexity parameter.

Like in the deterministic case, the online acceleration works better than its offline counterpart, itself improving the original optimization method, when we have a variance reduction technique.

Optimizing Convolutional Neural Networks

Introduction

Stochastic gradient descent is a popular and effective method to train neural networks [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF][START_REF] Deng | Recent advances in deep learning for speech research at Microsoft[END_REF].

A lot of efforts have been invested in accelerating stochastic algorithms, in particular by deriving methods that adapt to the structure of the problem. Algorithms such as RMSProp [START_REF] Tieleman | Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude[END_REF] or Adam [START_REF] Kingma | A method for stochastic optimization[END_REF] are examples of direct modifications of gradient descent, and estimate some statistical momentum during optimization to speed up convergence. Unfortunately, these methods can fail to converge on some simple problems [START_REF] Reddi | On the convergence of Adam and beyond[END_REF], and may fail to achieve state-of-the-art test accuracy on image classification problems. Other techniques have been developed to improve convergence speed or accuracy, such as adaptive batch-size for distributed SGD [START_REF] Goyal | Accurate, large minibatch sgd: training imagenet in 1 hour[END_REF], or quasi-second-order methods which improve the rate of convergence of stochastic algorithms [START_REF] Bollapragada | A progressive batching L-BFGS method for machine learning[END_REF]. However, only a limited number of settings are covered by such techniques, which are not compatible with state-of-the-art architectures.

The approach we propose here is completely different as RNA algorithm is built on top of existing optimization algorithms. Classical algorithms typically retain only the last iterate or the average [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF] of iterates as their best estimate of the optimum, throwing away all the information contained in the converging sequence of iterates. As it is highly wasteful from a statistical perspective, extrapolation schemes estimate the optimum of an optimization problem using a weighted average of the last iterates produced by an algorithm, where the weights depend on the iterates.

Network ensembling [START_REF] Zhou | Ensembling neural networks: many could be better than all[END_REF] can also be seen as combining several neural networks to improve convergence. However, contrary to our strategy, ensembling consists in training different networks from different starting points and then averaging the predictions or the parameters. Averaging the weights of successive different neural networks from SGD iterations has been studied in [START_REF] Izmailov | Averaging weights leads to wider optima and better generalization[END_REF],

and RNA can also be seen as an extension of this averaging idea for SGD, but with non-uniform adaptive weights.

Overall, nonlinear acceleration has marginal computational complexity. On neural network training problems, we will see the offline version improves both the test accuracy of early iterations, as well as the final accuracy with only minor modifications of existing learning pipelines. It never hurts performance as it runs on top of the algorithm and does not affect iterations. Finally, the scheme produces smoother learning curves. When training neural networks, we observe in the experiments that the convergence speedup produced by acceleration is much more significant in early iterations, which means it could serve for "rapid prototyping" of network architectures, with significant computational savings.

Numerical Experiments

The following numerical experiments seek to highlight the potential benefits of RNA in its offline and online versions when applied to the gradient method (with or without momentum term).

Since the complexity grows quadratically with the number N of points in the sequences {x i } and {y i }, we will use RNA with a fixed window size (N = 5 instead of N = 10 used previously, because neural networks has millions of parameters) in all these experiments. These values are sufficiently large to show a significant improvement in the rate of convergence, but can of course be fine-tuned. For simplicity, like in the previous experiments we fix λ = 10 be the parameters after one stochastic gradient step.

Then, for a data set of size m,

x i -y i-1 ≈ 1 m m t=1 (x (t+1) i -y (t) i ) = -h 1 m m t=1 ∇f (y (t) i ),
where h is the stepsize, also called learning rate in the community of deep neural networks. This means the matrix R in Algorithm 3 will be the matrix of (estimated) gradients.

Because the learning curve is highly dependent on the learning rate schedule, we decided to use a linearly decaying learning rate to better illustrate the benefits of acceleration, even if acceleration also works with a constant learning rate schedule (see [START_REF] Scieur | Nonlinear acceleration of CNNs[END_REF] and We then continue the optimization during 10 additional epochs using h T to stabilize the curve. We summarize the parameters used for the optimization in of SGD with and without momentum. In addition, we show a substantial improvement of the accuracy over the non-accelerated baseline. The improvement in the accuracy is reported in Figure 4.5. Interestingly, the resulting training loss is smoother than its non accelerated counterpart, which indicates a noise reduction.

Conclusion and Perspectives

In this chapter, we analyzed the effect of acceleration on several optimization algorithms, in the deterministic and stochastic setting and even for non-convex neural networks. In the deterministic case, we studied the local rate of convergence, and we deduced that when λ = O(D 2(1+a) )

with a ∈]0, 1[ where D is an upper bound for x 0x * , then we recover asymptotically an Chapter 5

Conclusion and Perspectives

In this thesis, we studied acceleration under different forms. First, we analyzed existing accelerated algorithms under another point of view, linking them with integration methods for the gradient flow. In particular, we showed that many existing optimization methods are special instances of standard integration schemes. If we assume classical conditions, such as consistency or stability, we showed that accelerated algorithms present larger step size, thus integrating (and converging) faster. Second, we modified and studied an extrapolation algorithm for vector sequences (Anderson Acceleration and Minimal Polynomial Method) and extended the analysis to optimization algorithms. We showed that this scheme, when regularized, improves significantly optimization methods performance, even when they are already accelerated. This is explained by its adaptivity to the local structure of the problem.

In the future, we hope that our links with integration methods can be extended. In particular, we are able to design methods with optimal rates, but only for quadratics. The extension to nonlinear function can be eventually done by studying G-stability [START_REF] Dahlquist | G-stability is equivalent toa-stability[END_REF][START_REF] Butcher | Thirty years of g-stability[END_REF], which generalizes the study of integration schemes to smooth and strongly convex functions. On the other hand, strong convexity is a key component in our approach, and moving to convex function can be an interesting extension.

The analysis of the regularized nonlinear acceleration algorithm, even generic, does not depend on the problem's structure and the convergence bound is very conservative. If we introduce strong convexity and smoothness sooner in the analysis, it may be possible to derive more accurate convergence bounds. In addition, we need to use a regularized version of Chebyshev polynomial, whose exact formula is unknown. We also assumed the operator which generates the sequences to be symmetric. Considering a non-symmetric operator (for example when using primal-dual algorithms) can be interesting. Finally, in other experiments, the efficiency of RNA on non-smooth problem looks promising. However our analysis does not hold because by essence non-smooth problems are not differentiable. Using other arguments, it might be possible to show that generic acceleration is also achievable for sharp problems. 97 Résumé Dans de nombreux domaines, comme par exemple l'optimisation, la performance d'une méthode est souvent caractérisée par son taux de convergence. Cependant, accélérer un algorithme requiert une certaine connaissance de la structure du problème et de telles améliorations sont le fruit d'une étude au cas-par-cas. De nombreuses techniques d'accélération ont été développées ces dernières décennies et sont maintenant massivement utilisées. En dépit de leur simplicité, ces méthodes sont souvent basées sur des arguments purement algébriques et n'ont généralement pas d'explications intuitives. Récemment, de nombreux travaux ont été menés pour faire des liens entre les algorithmes accélérés et d'autres domaines scientifiques, comme par exemple avec la théorie du contrôle ou des équations différentielles. Cependant, ces explications reposent souvent sur des arguments complexes et la plupart utilise des outils nonconventionnels dans leur analyse. Une des contributions de cette thèse est une tentative d'explication des algorithmes accélérés en utilisant la théorie des méthodes d'intégration, qui a été très étudiée et jouit d'une analyse théorique solide. En particulier, nous montrons que les méthodes d'optimisation sont en réalité des instances de méthode d'intégration, lorsqu'on intègre l'équation du flot de gradient. Avec des arguments standards, nous expliquons intuitivement l'origine de l'accélération. De l'autre côté, les méthodes accélérées ont besoin de paramètres supplémentaires, en comparaison d'autres méthodes plus lentes, qui sont généralement difficiles à estimer. De plus, ces schémas sont construits pour une configuration particulière et ne peuvent pas être utilisés autre part. Ici, nous explorons une autre approche pour accélérer les algorithmes d'optimisation, qui utilise des arguments d'accélération générique. En analyse numérique, ces outils ont été développés pour accélérer des séquences de scalaires ou de vecteurs, en construisant parallèlement une autre séquence avec un meilleur taux de convergence. Ces méthodes peuvent être combinées avec un algorithme itératif, l'accélérant dans la plupart des cas. En pratique, ces méthodes d'extrapolation ne sont pas tellement utilisées, notamment dû à leur manque de garanties de convergence et leur instabilité. Nous étendons ces méthodes en les régularisant, ce qui permettra une analyse théorique plus profonde et des résultats de convergence plus fort, en particulier lorsqu'elles sont appliquées à des méthodes d'optimisation.

Abstract

In many different fields such as optimization, the performance of a method is often characterized by its rate of convergence. However, accelerating an algorithm requires a lot of knowledge about the problem's structure, and such improvement is done on a case-by-case basis. Many accelerated schemes have been developed in the past few decades and are massively used in practice. Despite their simplicity, such methods are usually based on purely algebraic arguments and often do not have an intuitive explanation. Recently, heavy work has been done to link accelerated algorithms with other fields of science, such as control theory or differential equations. However, these explanations often rely on complex arguments, usually using non-conventional tools in their analysis. One of the contributions of this thesis is a tentative explanation of optimization algorithms using the theory of integration methods, which has been well studied and enjoys a solid theoretical analysis. In particular, we will show that optimization scheme are special instance of integration methods when integrating the classical gradient flow. With standard arguments, we intuitively explain the origin of acceleration. On the other hand, accelerated methods usually need additional parameters in comparison with slower one, which are in most cases difficult to estimate. In addition, these schemes are designed for one particular setting and cannot be used elsewhere. In this thesis, we explore a new approach for accelerating optimization algorithms, which uses generic acceleration arguments. In numerical analysis, these tools have been developed for accelerating sequences of scalars or vectors, by building on the side another sequence with a better convergence rate. These methods can be combined with an iterative algorithm, speeding it up in most cases. In practice, extrapolation schemes are not widely used due to their lack of theoretical guarantees and their instability. We will extend these methods by regularizing them, allowing a deeper theoretical analysis and stronger convergence results, especially when applied to optimization methods.

Proposition 2 . 3 . 4 .

 234 A linear multi-step method defined by polynomials (ρ, σ) is consistent if and only if ρ(1) = 0 and ρ (1) = σ(1).

Figure 2 . 1 :

 21 Figure2.1: Integration of a (Linear ODE) (which corresponds to the minimization of a quadratic function) using Euler's, Nesterov's and Polyak's methods between [0, t max ]. On the left, the optimal step size is used. Because Polyak's algorithm is the one with the biggest step size, it needs less iterations than Nesterov or Euler to approximate x(t max ). On the right side, we reduced the step size to 1/L for all methods. We clearly observe that they all track the same ODE: the gradient flow.

Proposition 3 . 2 . 1 .

 321 Let c * be the solution of (3.6) when we ensure 1 T c = 1. Then c * can be also recovered by minimizing the combination of residuals (3.3), i.e., c * = argmin c:c T 1=1 N i=0 c i r i .(3.7)

Figure 3 . 1 :

 31 Figure 3.1: We plot both T 3 (x, κ) (blue) and T 5 (x, κ) (black) for x ∈ [0, 1] and σ = 0.85. The maximum value of the image of [0, 1κ] by T N is clearly smaller when N grows, implying a better rate of convergence.

2 :Rc 3 :

 23 Solve c * = argmin {c∈R N +1 : c T 1=1} Compute x extr = N i=0 c i x i and g(x extr ) = N i=0 c i x i+1 . Ensure: Residual g(x extr )x extr small.

( 1 )

 1 needs to solve an optimization problem to find the c i . Hopefully, there exists a simpler way to perform the computation. The next proposition gives us an explicit solution to step 2 in Algorithm 1, involving the inversion of a N -by-N matrix. Proposition 3.2.4. The explicit solution of the problem c * = argmin c T 1=1 Rc (3.14)

Proposition 3 . 3 . 1 .

 331 Let the sequences x i be generated by (LFPI) and xi by (Pert. LFPI), with x 0 = x0 , with R and R the residual matrices defined above with RT R invertible and P the perturbation matrix in (3.17). Assume c * and c * are computed using formula (3.15) with matrices R and R respectively. Let ∆c * c *c * . (3.18)

Algorithm 2 1 :

 21 Regularized Nonlinear Acceleration (RNA) Require: Iterates x0 , x1 , ..., xN+1 ∈ R d produced by (Pert. LFPI), and a regularization parameter λ > 0. Compute R = [r 0 , ..., rN ], where ri = xi+1xi 2: Solve c * λ = argmin

Proposition 3 . 3 . 2 .

 332 Consider the sequences x i satisfying (LFPI) and xi satisfying (Pert. LFPI) with x 0 = x0 . Let c * λ and c * λ the output of Algorithm 2 with parameter λ applied to x i and xi respectively. Let R and R the matrices of residuals and P be defined in (3.17). Define ∆c * λ = c * λc * λ . Then, we have the following bounds, which control the stability of the solution c * λ . Proof. Using the same proof technique of Propositions 3.2.4 and 3.3.1, we have

3 . 2 )R 2 (γ 2

 3222 thanks to the parameter λ. The next proposition shows the relation between γ and λ when we use Algorithm 2 on (Pert. LFPI). Proposition 3.3.4. Let the sequence xi generated from (Pert. LFPI). Assume we want to use Algorithm 2 on this sequence, with c * λ ≤ γ √ 1+N . In this case, λ should satisfies λ ≥ -1) . (3.27) Proof. Using the result from Proposition 3.3.2,

Algorithm 3 1 :

 31 resulting in a time-varying matrix G i . Hopefully, it is possible to slightly modify RNA to overcome this difficulty. It suffices to use the residues r i = x i+1y i instead of x i+1x i as described in Algorithm 3. Regularized Nonlinear Acceleration (Improved version) Require: Iterates ỹ0 , ỹ1 , ..., ỹN ; x1 , x2 , ..., xN+1 from (Comb. Pert. LFPI), regularization parameter λ > 0. Compute R = [r 0 , ..., rN ], where ri = xi+1ỹi 2: Solve c * λ = argmin

Proposition 3 . 5 . 2 .

 352 Assume we apply Algorithm 3 on a sequence of vectors generated by (Comb. Pert. LFPI), where the norm of the perturbation e i is bounded by . Let the norm of the coefficients c * λ computed by Algorithm 3 be bounded by γ √ N +1 with γ ≥ 1. Then, the accuracy x extrx * is bounded by

3 . 10 .Dγ = 0 .Proposition 4 . 2 . 2 .

 3100422 It suffices to set γ in function of D such that lim D→0 γ = ∞, lim D→0 The next proposition gives the relation between D and the value of γ (and λ) to satisfy the above condition. If we have, for 0 < a < 1, γ = 1 + O D -a or λ = O D 2(1+a) ,

Figure 4 . 1 :

 41 Figure 4.1: Comparison between non-regularized RNA (or Anderson acceleration) and regularized nonlinear acceleration algorithm in the offline fashion (Left: λ = 0 and Right:λ = 10 -8 RT W R ).The window size (maximum sequence length) is set to 10. We run the experiment on Madelon dataset, whole regularization is set such that the condition number κ = 10 -6 . We clearly see that regularization is crucial for the stability of the algorithm. (V) (vanilla) corresponds to the original methods, and P (parallel) corresponds to the application of RNA without interaction with the original algorithm.

Figure 4 . 6 )

 46 . In all our experiments, until epoch T , the learning rate decreases linearly from an initial value h 0 to a final value h T , withh k = h 0 + (k/T )(h Th 0 ).(4.8)

Figure 4 . 5 :

 45 Figure 4.5: Training a Resnet-18 (top), VGG19 (middle) and a DenseNet-121 (bottom) on Cifar10for 100 epochs. SGD with and without momentum, and their off-line accelerated versions with a window size of 5. On the left, the evolution of the loss value. On the right, the error percentage on the testing set.

Figure

  Figure 4.6: Prototyping networks: acceleration of SGD (bottom) gives a smoother convergence, producing a clearer ranking of architectures, earlier (flat learning rate). Right plot zooms on left one. The thin line corresponds to SGD, the bold to RNA+SGD.

Table 4

 4 .1. and 4.4. This suggests that the offline version of RNA is a good candidate for training neural networks, as it converges faster while guaranteeing performance at least as good as the reference algorithm.

Table 4 .

 4 4.6: Prototyping networks: acceleration of SGD (bottom) gives a smoother convergence, producing a clearer ranking of architectures, earlier (flat learning rate). Right plot zooms on left one. The thin line corresponds to SGD, the bold to RNA+SGD. 2: Number of epochs required to reach the best test accuracy + Tolerance% on CIFAR10 with a Resnet18, using several algorithms (best accuracy is 5% here). The speed-up compared to the SGD+momentum baseline is in parenthesis.

	Tolerance	SGD	SGD+momentum SGD+RNA SGD+momentum+RNA
	5.0%	68 (0.87×)	59	21 (2.81×)	16 (3.69×)
	2.0%	78 (0.99×)	77	47 (1.64×)	40 (1.93×)
	1.0%	82 (1.00×)	82	67 (1.22×)	59 (1.39×)
	0.5%	84 (1.02×)	86	75 (1.15×)	63 (1.37×)
	0.2%	86 (1.13×)	97	84 (1.15×)	85 (1.14×)
	Tolerance	SGD	SGD+momentum SGD+RNA SGD+momentum+RNA
	5.0%	69 (0.87×)	60	26 (2.31×)	24 (2.50×)
	2.0%	83 (0.99×)	82	52 (1.58×)	45 (1.82×)
	1.0%	84 (1.02×)	86	71 (1.21×)	60 (1.43×)
	0.5%	89 (0.98×)	87	73 (1.19×)	62 (1.40×)
	0.2%	N/A	90	99 (0.90×)	63 (1.43×)

Table 4 .

 4 3: Number of epochs required to reach the best test accuracy + Tolerance% on CIFAR10 with VGG19. Here, the best accuracy is 6.54%ImagenetHere, we apply the RNA algorithm to the standard ImageNet dataset. We trained the networks during 90 epochs (T = 80) with a weight decay of 10 -4 (it corresponds to add a 2 regularization in the objective). We reported the test accuracy on Figure4.7 for the networks ResNet-50 and ResNet-152. We only tested the offline version of RNA here, because in previous experiments it gives better result than its online counterpart.We again observe that the offline version of Algorithm 3 improves the convergence speed

	Tolerance	SGD	SGD+momentum SGD+RNA SGD+momentum+RNA
	5.0%	65 (0.86×)	56	22 (2.55×)	13 (4.31×)
	2.0%	80 (0.98×)	78	45 (1.73×)	38 (2.05×)
	1.0%	83 (1.00×)	83	60 (1.38×)	56 (1.48×)
	0.5%	87 (0.99×)	86	80 (1.08×)	66 (1.30×)
	0.2%	92 (1.01×)	93	86 (1.08×)	75 (1.24×)

Table 4 .

 4 4: Number of epochs required to reach the best test accuracy + Tolerance% on CIFAR10 with DenseNet121. Here, the best accuracy is 4.62%.Figure 4.7: Training a Resnet-52 (left) and ResNet-152 (right) on validation ImageNet for 90 epochs using SGD with and without momentum, and their off-line accelerated versions.

		50					50
	Percentage validation error	25 30 35 40 45	SGD SGD + momentum SGD + RNA SGD + momentum + RNA	Percentage validation error	25 30 35 40 45	SGD + momentum SGD + momentum + RNA
		20	1	30	60	90	20	1	30	60	90
				Epoch						Epoch
				Pytorch	SGD	SGD+mom.		SGD+RNA	SGD+mom.+RNA
	Resnet-50	23.85	23.808	23.346	23.412 (-0.396%) 22.914 (-0.432%)
	Resnet-152	21.69	N/A	21.294		N/A		20.884 (-0.410%)

Table 4 .

 4 5: Best validation top-1 error percentage on ImageNet. In parenthesis the improvement due to RNA. The first column corresponds to the performance of Pytorch pre-trained models.

Chapter 4

Nonlinear Acceleration of Optimization Algorithms

Abstract

In this chapter we will apply the generic acceleration technique RNA to several optimization methods in different settings, such as deterministic or stochastic optimization. In the deterministic setting, we will explicitly write the convergence bounds using standard assumption on the objective function, such as smoothness and strong convexity. For the stochastic one, we will analyze the accuracy bounds in function of the noise level and deduce that variance reduction techniques are crucial for acceleration. We provide numerical experiments in both cases. In addition, we use it to train deep neural networks on CIFAR10 and ImageNet data sets, and empirically observe that a straightforward application of RNA makes networks more accurate.

Main references: The results in this chapter have been partially published in [START_REF] Scieur | Regularized nonlinear acceleration[END_REF], [START_REF] Scieur | Nonlinear acceleration of stochastic algorithms[END_REF] (NIPS 2016 and 2017), [START_REF] Scieur | Nonlinear acceleration of CNNs[END_REF] (ICLR workshop 2018) and [START_REF] Scieur | Nonlinear acceleration of deep neural networks[END_REF] (under review NIPS 2018). A revised version of [START_REF] Scieur | Regularized nonlinear acceleration[END_REF] is under review for the journal Mathematical Programming.

Introduction

In this chapter, we are interested in minimizing a smooth convex function f (x), min

in the variable x ∈ R n . To do so, we will use a first-order algorithm which can be written under the form of (Comb. Pert. LFPI) (we will drop the tilde ˜notation in this chapter for clarity).

However, we often equip it with a (backtracking) line-search to improve its performance. This means we have a variable step-size h i , and the method is written

This can be captured by coefficients α (i) and β (i) as following,

In this case, we can verify that

Fast Gradient Method

We now write the Nesterov's method [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF] into the form of (Comb. Pert. LFPI),

κ if we know the strong convexity constant µ. In this case, the algorithm converges with rate (1 -√ κ) i . Otherwise, we can use θ i = i 1+i , which has a rate of convergence of the order of 1 1+i 2 . If we want to use the adaptive momentum Algorithm 5, the temporary quantity z i is written

Then, it suffices to check if

In that way, we preserve the original global rate of Nesterov's method.

Numerical Experiments (Deterministic)

We will compare different methods for the optimization of a logistic loss on three datasets, named Sonar [START_REF] Gorman | Analysis of hidden units in a layered network trained to classify sonar targets[END_REF], Madelon [START_REF] Guyon | Feature extraction: foundations and applications[END_REF] and sido0 [START_REF] Guyon | Sido: A phamacology dataset[END_REF]. The logistic loss is defined as follows. Let a dataset with feature matrix Θ ∈ R m×d which contains m vector of d features each. Consider the {-1, 1}

vector ϑ which assigns each vectors of Θ to a positive or negative label. Then the logistic loss

Comparison Between Acceleration Strategies

In addition, we run the -BFGS algorithm. We use three different regularization parameters to tune the condition number κ, which are well conditioned (κ = 10 -3 ), normally conditioned

) and badly conditioned (κ = 10 -9 ). In these experiments, for the RNA algorithm, we use as regularization parameter λ = 10 -8 RT R , and a window size of 10. Of course, these parameters can be fine-tuned to enjoy the best possible performance, but for simplicity we chose to keep these values fixed. For the -BFGS algorithm we also have a window size of 10 gradients, and a Wolfe-line-search.

We used an unified color and linestyle specification, summarized in where f (x) is a smooth and strongly convex function with respect to the Euclidean norm. We consider a stochastic first-order oracle, which gives a noisy estimate of the gradient of f (x), with

where ε is a noise term with bounded variance. This is the case for example when f is a sum of strongly convex functions, and we only have access to the gradient of one randomly selected function. Stochastic optimization (4.7) is typically challenging as classical algorithms are not convergent for a constant step size (for example, gradient descent or Nesterov's accelerated gradient). Even the averaged version of stochastic gradient descent with constant step size does not converge to the solution of (4.6), but to another point whose proximity to the real minimizer depends of the step size [START_REF] Nedić | Convergence rate of incremental subgradient algorithms[END_REF][START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF].

When f is a finite sum of N functions, then algorithms such as SAG [START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF], SAGA [START_REF] Defazio | Saga: A fast incremental gradient method with support for non-strongly convex composite objectives[END_REF],

SDCA [START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss minimization[END_REF] and SVRG [START_REF] Johnson | Accelerating stochastic gradient descent using predictive variance reduction[END_REF] accelerate convergence using a variance reduction technique akin to control variate in Monte-Carlo methods. Their rate of convergence depends on 1µ/L and thus does exhibit an accelerated rate on par with the deterministic setting (in 1µ/L), where L is the smoothness constant and µ the strong convexity constant of f . Recently a generic acceleration algorithm called Catalyst [START_REF] Lin | A universal catalyst for first-order optimization[END_REF], based on the proximal point methods improved this rate of convergence, at least in theory. Unfortunately, numerical experiments show this algorithm to be conservative, thus limiting practical performances. On the other hand, recent papers, for example [START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF] (Accelerated SDCA) and [START_REF] Allen-Zhu | Katyusha: The first direct acceleration of stochastic gradient methods[END_REF] (Katyusha), propose algorithms with accelerated convergence rates, if the strong convexity parameter is given.

When f is a quadratic function then averaged SGD converges, but the rate of decay of initial conditions is very slow. Recently, some results have focused on accelerated versions of SGD for quadratic optimization, showing that with a two step recursion it is possible to enjoy both the optimal rate for the bias term and the variance [START_REF] Flammarion | From averaging to acceleration, there is only a step-size[END_REF], given an estimate of the ratio between the distance to the solution and the variance of ε.

In this section, we will analyse the impact of nonlinear acceleration to stochastic algorithms.

We hope this strategy will still be adaptive to the strong convexity constant, whose online estimation is still a challenge, even in the deterministic case [START_REF] Fercoq | Restarting accelerated gradient methods with a rough strong convexity estimate[END_REF].

Before going in the numerical part, we analyse quickly the impact of stochastic noise on the extrapolation.

Acceleration with Noisy Iterates

We have seen here that in addition with non-linearities, a stochastic noise is present in (Comb. Pert. LFPI). We already analyzed the impact of non-linearities in the previous section, so here we focus only on stochastic noises with zero-mean and variance σ 2 . In this case, match perfectly (Comb. Pert. LFPI).

We compare these various methods for minimizing a logistic loss on Sonar [START_REF] Gorman | Analysis of hidden units in a layered network trained to classify sonar targets[END_REF], Madelon [START_REF] Guyon | Feature extraction: foundations and applications[END_REF] and sido0 [START_REF] Guyon | Sido: A phamacology dataset[END_REF], with several condition numbers κ: well (κ = 100/m), moderately (κ = 1/m) and badly (κ = 1/100m) conditioned.

Like before, we compare the original optimization algorithm with its offline and online acceleration using RNA:

• (V) Vanilla: we run the original algorithm with no acceleration.

• (P) Offline version, in parallel: we run the RNA algorithm on the side (Algorithm 3), with no interaction on the original optimization method.

• (O) Online: we use the online version of RNA (Algorithm 4). We do not use Algorithm 5, because a comparison between function values (i.e., computing h(x)) is costly in stochastic optimization as it requires one data pass.

Observations

In the previous figures, we clearly see that both SGD and its accelerated versions do not converge. This is mainly due to the fact that the variance is contant over iterations as we do not perform averaging. In addition, excepted for quadratic problems, the averaged version of SGD with constant step size does not converge with arbitrary precision.

h 0 h T momentum SGD and Online RNA 1.0 0.01 0 SGD + momentum 0.1 0.001 0.9 Table 4.1: Parameters used in (4.8) to generate the learning rate for optimizers. We used the same setting for their accelerated version with RNA. Algorithm 4). On the left, the evolution of the loss value. On the right, the error percentage on the testing set. We see the online version does not improve the performance of SGD. This may be explained by the presence of the batchnorm module. The presence of a statistical estimation while the optimization is running is not taken into account in the iteration model (Comb. Pert. LFPI), and thus makes theoretical analysis more complex. In addition, because we are optimizing a nonconvex function, and the online version interferes with the optimization scheme, the final performance is not guaranteed to be as good as the non-accelerated version of the optimizer. optimal rate of convergence when D → 0.

In the case of stochastic algorithm, we noticed that the rate mostly depends on the ratio between the noise variance and the distance to the optimum, which is sort of the inverse of the signal to noise ratio. This means acceleration cannot work on algorithm with large variance, such as SGD. However, if the distance between two stochastic iteration is big enough, then we can accelerate an optimization algorithm with variance-reduction technique such as SAGA or SVRG. Waiting several stochastic iteration is also beneficial since it allows to balance the time taken by the algorithm's iterations and RNA computation.

In the two previous cases, offline and online acceleration improve the numerical performance of optimization algorithms when minimizing a logistic loss, whatever the condition number of the problem is. Even compared with optimal methods, such as Nesterov's accelerated gradient (for deterministic optimization) and Katyusha (Stochastic optimization), the online RNA algorithm is much faster, due to its adaptivity to the problem's structure.

In the case of neural network, the conclusion is more mitigated, as the online version does not works as well as expected. In particular due to non-convexity, high dimension and statistical estimations during the optimization process (batchnorm modules), convonlutional neural networks are extremely hard to optimize. This may explain why online acceleration does not perform as well as in the convex case.

On the contrary to its online counterpart, offline acceleration preserves the performance as is does not interact with the optimization scheme. Surprisingly, this gives good acceleration results on neural networks. Even after a few number of iteration, the prediction accuracy of the extrapolated model is close to be optimal.

Here, we did not explore the limits of extrapolation. For example, is it possible to design an accelerated method for non-smooth problem? Even if some numerical experiments have been conducted in [START_REF] Scieur | Regularized nonlinear acceleration[END_REF] on Lasso or dual-SVM, also in [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF] on dual-Lasso, the theoretical background is still incomplete, since non-smooth problems are inherently non-differentiable.

Also, the similarity between BFGS and RNA performance may mean there is a link between the two methods. In fact, in [START_REF] Fang | Two classes of multisecant methods for nonlinear acceleration[END_REF] is studied several classes of non-linear acceleration methods, such as Anderson Acceleration (highly related to RNA) and Broyden's family (linked with BFGS). A deeper analysis may lead to new convergence results.