Nicolas Flammarion
email: nicolas.flammarion@ens.fr

Francis Bach
email: francis.bach@ens.fr

From Averaging to Acceleration, There is Only a Step-size

We show that accelerated gradient descent, averaged gradient descent and the heavyball method for non-strongly-convex problems may be reformulated as constant parameter second-order difference equation algorithms, where stability of the system is equivalent to convergence at rate O(1/n 2 ), where n is the number of iterations. We provide a detailed analysis of the eigenvalues of the corresponding linear dynamical system, showing various oscillatory and non-oscillatory behaviors, together with a sharp stability result with explicit constants. We also consider the situation where noisy gradients are available, where we extend our general convergence result, which suggests an alternative algorithm (i.e., with different step sizes) that exhibits the good aspects of both averaging and acceleration.

Introduction

Many problems in machine learning are naturally cast as convex optimization problems over a Euclidean space; for supervised learning this includes least-squares regression, logistic regression, and the support vector machine. Faced with large amounts of data, practitioners often favor first-order techniques based on gradient descent, leading to algorithms with many cheap iterations. For smooth problems, two extensions of gradient descent have had important theoretical and practical impacts: acceleration and averaging.

Acceleration techniques date back to [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF] and have their roots in momentum techniques and conjugate gradient [START_REF] Polyak | Introduction to Optimization[END_REF]. For convex problems, with an appropriately weighted momentum term which requires to store two iterates, [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF] showed that the traditional convergence rate of O(1/n) for the function values after n iterations of gradient descent goes down to O(1/n 2 ) for accelerated gradient descent, such a rate being optimal among first-order techniques that can access only sequences of gradients [START_REF] Nesterov | Introductory Lectures on Convex Optimization[END_REF]. Like conjugate gradient methods for solving linear systems, these methods are however more sensitive to noise in the gradients; that is, to preserve their improved convergence rates, significantly less noise may be tolerated [START_REF]Smooth optimization with approximate gradient[END_REF][START_REF] Schmidt | Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization[END_REF][START_REF] Devolder | First-order methods of smooth convex optimization with inexact oracle[END_REF].

1 Averaging techniques which consist in replacing the iterates by the average of all iterates have also been thoroughly considered, either because they sometimes lead to simpler proofs, or because they lead to improved behavior. In the noiseless case where gradients are exactly available, they do not improve the convergence rate in the convex case; worse, for stronglyconvex problems, they are not linearly convergent while regular gradient descent is. Their main advantage comes with random unbiased gradients, where it has been shown that they lead to better convergence rates than the unaveraged counterparts, in particular because they allow larger step-sizes [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Bach | Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Machine Learning[END_REF]. For example, for least-squares regression with stochastic gradients, they lead to convergence rates of O(1/n), even in the non-strongly convex case [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n)[END_REF].

In this paper, we show that for quadratic problems, both averaging and acceleration are two instances of the same second-order finite difference equation, with different step-sizes. They may thus be analyzed jointly, together with a non-strongly convex version of the heavyball method (Polyak, 1987, Section 3.2). In presence of random zero-mean noise on the gradients, this joint analysis allows to design a novel intermediate algorithm that exhibits the good aspects of both acceleration (quick forgetting of initial conditions) and averaging (robustness to noise).

In this paper, we make the following contributions:

-We show in Section 2 that accelerated gradient descent, averaged gradient descent and the heavy-ball method for non-strongly-convex problems may be reformulated as constant parameter second-order difference equation algorithms, where stability of the system is equivalent to convergence at rate O(1/n 2 ).

-In Section 3, we provide a detailed analysis of the eigenvalues of the corresponding linear dynamical system, showing various oscillatory and non-oscillatory behaviors, together with a sharp stability result with explicit constants.

-In Section 4, we consider the situation where noisy gradients are available, where we extend our general convergence result, which suggests an alternative algorithm (i.e., with different step sizes) that exhibits the good aspects of both averaging and acceleration.

-In Section 5, we illustrate our results with simulations on synthetic examples.

2 Second-Order Iterative Algorithms for Quadratic Functions

Throughout this paper, we consider minimizing a convex quadratic function f : R d → R defined as:

f (θ) = 1 2 θ, Hθ -q, θ , (1)
with H ∈ R d×d a symmetric positive semi-definite matrix and q ∈ R d . Without loss of generality, H is assumed invertible (by projecting onto the orthogonal of its null space), though its eigenvalues could be arbitrarily small. The solution is known to be θ * = H -1 q, but the inverse of the Hessian is often too expensive to compute when d is large. The excess cost function may be simply expressed as f (θ n ) -f (θ * ) = 1 2 θ n -θ * , H(θ n -θ * ) .

Second-order algorithms

In this paper we study second-order iterative algorithms of the form:

θ n+1 = A n θ n + B n θ n-1 + c n , (2)
started with θ 1 = θ 0 in R d , with A n ∈ R d×d , B n ∈ R d×d and c n ∈ R d for all n ∈ N * . We impose the natural restriction that the optimum θ * is a stationary point of this recursion, that is, for all n ∈ N * :

θ * = A n θ * + B n θ * + c n .
(θ * -stationarity)

By letting φ n = θ n -θ * we then have φ n+1 = A n φ n +B n φ n-1 , started from φ 0 = φ 1 = θ 0 -θ * . Thus, we restrict our problem to the study of the convergence of an iterative system to 0.

In connection with accelerated methods, we are interested in algorithms for which f (θ n )f (θ * ) = 1 2 φ n , Hφ n converges to 0 at a speed of O 1/n 2 . Within this context we impose that A n and B n have the form :

A n = n n + 1 A and B n = n -1 n + 1 B ∀n ∈ N with A, B ∈ R d×d . (n-scalability)
By letting η n = nφ n = n(θ n -θ * ), we can now study the simple iterative system with constant terms η n+1 = Aη n + Bη n-1 , started at η 0 = 0 and η 1 = θ 0 -θ * . Showing that the function values remain bounded, we directly have the convergence of f (θ n ) to f (θ * ) at the speed O 1/n 2 . Thus the n-scalability property allows to switch from a convergence problem to a stability problem.

For feasibility concerns the method can only access H through matrix-vector products. Therefore A and B should be polynomials in H and c a polynomial in H times q, if possible of low degree. The following theorem clarifies the general form of iterative systems which share these three properties (see proof in Appendix B).

Theorem 1. Let (P n , Q n , R n ) ∈ (R[X]
) 3 for all n ∈ N, be a sequence of polynomials. If the iterative algorithm defined by Eq. ( 2) with A n = P n (H), B n = Q n (H) and c n = R(H)q satisfies the θ * -stationarity and n-scalability properties, there are polynomials ( Ā, B) ∈ (R[X]) 2 such that:

A n = 2 n n + 1 I - Ā(H) + B(H) 2 H , B n = - n -1 n + 1 I -B(H)H and c n = n Ā(H) + B(H) n + 1 q.
Note that our result prevents A n and B n from being zero, thus requiring the algorithm to strictly be of second order. This illustrates the fact that first-order algorithms as gradient descent do not have the convergence rate in O(1/n 2 ).

We now restrict our class of algorithms to lowest possible order polynomials, that is, Ā = αI and B = βI with (α, β) ∈ R 2 , which correspond to the fewest matrix-vector products per iteration, leading to the constant-coefficient recursion for η n = nφ n = n(θ n -θ * ):

η n+1 = (I -αH) η n + (I -βH) (η n -η n-1 ) . (3)
Expression with gradients of f . The recursion in Eq. ( 3) may be written with gradients of f in multiple ways. In order to preserve the parallel with accelerated techniques, we rewrite it as:

θ n+1 = 2n n + 1 θ n - n -1 n + 1 θ n-1 - nα + β n + 1 f n(α + β) nα + β θ n - (n -1)β nα + β θ n-1 . (4)
It may be interpreted as a modified gradient recursion with two potentially different affine (i.e., with coefficients that sum to one) combinations of the two past iterates. This reformulation will also be crucial when using noisy gradients. The allowed values for (α, β) ∈ R 2 will be determined in the following sections.

Examples

Averaged gradient descent. We consider averaged gradient descent (referred to from now on as "Av-GD") [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF] with step-size γ ∈ R defined by:

ψ n+1 = ψ n -γf (ψ n ), θ n+1 = 1 n + 1 n+1 i=1 ψ i .
When computing the average online as θ n+1 = θ n + 1 n+1 (ψ n+1 -θ n ) and seeing the average as the main iterate, the algorithm becomes (see proof in Appendix B.2):

θ n+1 = 2n n + 1 θ n - n -1 n + 1 θ n-1 - γ n + 1 f nθ n -(n -1)θ n-1 .
This corresponds to Eq. ( 4) with α = 0 and β = γ.

Accelerated gradient descent. We consider the accelerated gradient descent (referred to from now on as "Acc-GD") [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O(1/k 2 )[END_REF] with step-sizes (γ, δ n ) ∈ R 2 :

θ n+1 = ω n -γf (ω n ), ω n = θ n + δ n (θ n -θ n-1 ).
For smooth optimization the accelerated literature [START_REF] Nesterov | Introductory Lectures on Convex Optimization[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] uses the step-size δ n = 1 -3 n+1 and their results are not valid for bigger step-size δ n . However δ n = 1 -2 n+1 is compatible with the framework of [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF] and is more convenient for our set-up. This corresponds to Eq. ( 4) with α = γ and β = γ. Note that accelerated techniques are more generally applicable, e.g., to composite optimization with smooth functions [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF][START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

Heavy ball. We consider the heavy-ball algorithm (referred to from now on as "HB") [START_REF] Polyak | Some methods of speeding up the convergence of iteration methods[END_REF] with step-sizes (γ, δ n ) ∈ R 2 :

θ n+1 = θ n -γf (θ n ) + δ n (θ n -θ n-1 ), when δ n = 1 -2 n+1 .
We note that typically δ n is constant for strongly-convex problems. This corresponds to Eq. ( 4) with α = γ and β = 0.

Convergence with Noiseless Gradients

We study the convergence of the iterates defined by: η n+1 = (I -αH) η n +(I -βH) (η n -η n-1 ). This is a second-order iterative system with constant coefficients that it is standard to cast in a linear framework (see, e.g., [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF]. We may rewrite it as:

Θ n = F Θ n-1 , with Θ n = η n η n-1 and F = 2I -(α + β) H βH -I I 0 ∈ R 2d×2d .
Thus Θ n = F n Θ 0 . Following O'Donoghue and Candes ( 2013), if we consider an eigenvalue decomposition of H, i.e., H = P Diag(h)P with P an orthogonal matrix and (h i ) the eigenvalues of H, sorted in decreasing order:

h d = L ≥ h d-1 ≥ • • • ≥ h 2 ≥ h 1 = µ > 0, then Eq.
(3) may be rewritten as:

P η n+1 = (I -αDiag (h)) P η n + (I -βDiag (h)) P η n -P η n-1 . (5)
Thus there is no interaction between the different eigenspaces and we may consider, for the analysis only, d different recursions with η i n = p i η n , i ∈ {1, ..., d}, where p i ∈ R d is the i-th column of P :

η i n+1 = (1 -αh i ) η i n + (1 -βh i ) η i n -η i n-1 . (6)

Characteristic polynomial and eigenvalues

In this section, we consider a fixed i ∈ {1, . . . , d} and study the stability in the corresponding eigenspace. This linear dynamical system may be analyzed by studying the eigenvalues of the 2 × 2-matrix

F i = 2 -(α + β)h i βh i -1 1 0
. These eigenvalues are the roots of its characteristic polynomial which is:

det(XI -F i ) = det (X (X -2 + (α + β)h i ) + 1 -βh i ) = X 2 -2X 1- α + β 2 h i +1-βh i .
To compute the roots of the second-order polynomial, we compute its reduced discriminant:

∆ i = 1 - α + β 2 h i 2 -1 + βh i = h i α + β 2 2 h i -α .
Depending on the sign of the discriminant ∆ i , there will be two real distinct eigenvalues (∆ i > 0), two complex conjugate eigenvalues (∆ i < 0) or a single real eigenvalue (∆ i = 0). In the interior of the triangle, the convergence is linear.

We will now study the sign of ∆ i . In each different case, we will determine under what conditions on α and β the modulus of the eigenvalues is less than one, which means that the iterates (η i n ) n remain bounded and the iterates (θ n ) n converge to θ * . We may then compute function values as

f (θ n ) -f (θ * ) = 1 2n 2 d i=1 (η i n ) 2 h i = 1 2 d i=1 (φ i n ) 2 h i .
The various regimes are summarized in Figure 1: there is a triangle of values of (αh i , βh i ) for which the algorithm remains stable (i.e., the iterates (η n ) n do not diverge), with either complex or real eigenvalues. In the following lemmas (see proof in Appendix C), we provide a detailed analysis that leads to Figure 1.

Lemma 1 (Real eigenvalues). The discriminant ∆ i is strictly positive and the algorithm is stable if and only if

α ≥ 0, α + 2β ≤ 4/h i , α + β > 2 α/h i .
We then have two real roots

r ± i = r i ± √ ∆ i , with r i = 1 -( α+β 2 )h i .
Moreover, we have:

(φ i n ) 2 h i = (φ i 1 ) 2 h i 4n 2 (r i + √ ∆ i ) n -(r i - √ ∆ i ) n 2 ∆ i . (7)
Therefore, for real eigenvalues, ((φ i n ) 2 h i ) n will converge to 0 at a speed of O(1/n 2 ) however the constant ∆ i may be arbitrarily small (and thus the scaling factor arbitrarily large). Furthermore we have linear convergence if the inequalities in the lemmas are strict.

Lemma 2 (Complex eigenvalues). The discriminant ∆ i is stricly negative and the algorithm is stable if and only if α ≥ 0, β ≥ 0, α + β < α/h i .

We then have two complex conjugate eigenvalues:

r ± i = r i ± √ -1 √ -∆ i .
Moreover, we have:

(φ i n ) 2 h i = (φ i 1 ) 2 n 2 sin 2 (ω i n) α -( α+β 2 ) 2 h i ρ 2n . ( 8
)
with ρ i = √ 1 -βh i , and ω i defined through sin(ω i ) = √ -∆ i /ρ i and cos(ω i ) = r i /ρ i .

Therefore, for complex eigenvalues, there is a linear convergence if the inequalities in the lemma are strict. Moreover, ((φ i n ) 2 h i ) n oscillates to 0 at a speed of O(1/n 2 ) even if h i is arbitrarily small. Coalescing eigenvalues. When the discriminant goes to zero in the explicit formulas of the real and complex cases, both the denominator and numerator of ((φ i n ) 2 h i ) n will go to zero. In the limit case, when the discriminant is equal to zero, we will have a double real eigenvalue. This happens for β = 2 α/h i -α. Then the eigenvalue is r i = 1-√ αh i , and the algorithm is stable for 0 < α < 4/h i , we then have ( -1) . This can be obtained by letting ∆ i goes to 0 in the real and complex cases (see also Appendix C.3).

φ i n ) 2 h i = h i (φ i 1 ) 2 (1 - √ αh i ) 2(n

Summary. To conclude the iterate (η

i n ) n = (n(θ i n -θ i * )) n will be stable for α ∈ [0, 4/h i ] and β ∈ [0, 2/h i -α/2].
According to the values of α and β this iterate will have a different behavior. In the complex case, the roots are complex conjugate with magnitude √ 1 -βh i . Thus, when β > 0, (η i n ) n will converge to 0, oscillating, at rate √ 1 -βh i . In the real case, the two roots are real and distinct. However the product of the two roots is equal to √ 1 -βh i , thus one will have a higher magnitude and (η i n ) n will converges to 0 at rate higher than in the complex case (as long as α and β belong to the interior of the stability region).

Finally, for a given quadratic function f , all the d iterates (η i n ) n should be bounded, therefore we must have α ∈ [0, 4/L] and β ∈ [0, 2/L -α/2]. Then, depending on the value of h i , some eigenvalues may be complex or real.

Classical examples

For particular choices of α and β, displayed in Figure 1, the eigenvalues are either all real or all complex, as shown in the table below.

Av-GD

Acc-GD Heavy ball

α 0 γ γ β γ γ 0 ∆ i (γh i ) 2 -γh i (1 -γh i ) -γh i (1 -γh i 4 ) r ± i 1, 1 -γh i √ 1 -γh i e ±iω i e ±iω i cos(ω i ) √ 1 -γh i 1 -γ 2 h i ρ i √ 1 -γh i 1
Averaged gradient descent loses linear convergence for strongly-convex problems, because r + i = 1 for all eigensubspaces. Similarly, the heavy ball method is not adaptive to strong convexity because ρ i = 1. However, accelerated gradient descent, although designed for non-strongly-convex problems, is adaptive because ρ i = √ 1 -γh i depends on h i while α and β do not. These last two algorithms have an oscillatory behavior which can be observed in practice and has been already studied [START_REF] Su | A Differential Equation for Modeling Nesterov's Accelerated Gradient Method: Theory and Insights[END_REF].

Note that all the classical methods choose step-sizes α and β either having all the eigenvalues real either complex; whereas we will see in Section 4, that it is significant to combine both behaviors in presence of noise.

General bound

Even if the exact formulas in Lemmas 1 and 2 are computable, they are not easily interpretable. In particular when the two roots become close, the denominator will go to zero, which prevents from bounding them easily. When we further restrict the domain of (α, β), we can always bound the iterate by the general bound (see proof in Appendix D):

Theorem 2. For α ≤ 1/h i and 0 ≤ β ≤ 2/h i -α, we have (η i n ) 2 ≤ min 2(η i 1 ) 2 αh i , 8(η i 1 ) 2 n (α + β)h i , 16(η i 1 ) 2 (α + β) 2 h 2 i . (9)
These bounds are shown by dividing the set of (α, β) in three regions where we obtain specific bounds. They do not depend on the regime of the eigenvalues (complex or real); this enables us to get the following general bound on the function values, our main result for the deterministic case.

Corollary 1. For α ≤ 1/L and 0 ≤ β ≤ 2/L -α: f (θ n ) -f (θ * ) ≤ min θ 0 -θ * 2 αn 2 , 4 θ 0 -θ * 2 (α + β)n . ( 10
)
We can make the following observations:

-The first bound θ 0 -θ * 2 αn 2 corresponds to the traditional acceleration result, and is only relevant for α > 0 (that is, for Nesterov acceleration and the heavy-ball method, but not for averaging). We recover the traditional convergence rate of second-order methods for quadratic functions in the singular case, such as conjugate gradient (Polyak, 1987, Section 6.1).

-While the result above focuses on function values, like most results in the non-strongly convex case, the distance to optimum θ n -θ * 2 typically does not go to zero (although it remains bounded in our situation).

-When α = 0 (averaged gradient descent), then the second bound 4 θ 0 -θ * 2 (α+β)n provides a convergence rate of O(1/n) if no assumption is made regarding the starting point θ 0 , while the last bound of Theorem 2 would lead to a bound 8

H -1/2 (θ 0 -θ * ) 2 (α+β) 2 n 2
, that is a rate of O(1/n 2 ), only for some starting points.

-As shown in Appendix E by exhibiting explicit sequences of quadratic functions, the inverse dependence in αn 2 and (α + β)n in Eq. ( 10) is not improvable.

Quadratic Optimization with Additive Noise

In many practical situations, the gradient of f is not available for the recursion in Eq. ( 4), but only a noisy version. In this paper, we only consider additive uncorrelated noise with finite variance.

Stochastic difference equation

We now assume that the true gradient is not available and we rather have access to a noisy oracle for the gradient of f . In Eq. ( 4), we assume that the oracle outputs a noisy gradient

f n(α+β) nα+β θ n -(n-1)β nα+β θ n-1 -ε n+1 . The noise (ε n ) is assumed to be uncorrelated zero-mean with bounded covariance, i.e., E[ε n ⊗ ε m ] = 0 for all n = m and E[ε n ⊗ ε n ] C, where A B means that B -A is positive semi-definite.
For quadratic functions, for the reduced variable η n = nφ n = n(θ n -θ * ), we get:

η n+1 = (I -αH)η n + (I -βH)(η n -η n-1 ) + [nα + β]ε n+1 . (11)
Note that algorithms with α = 0 will have an important level of noise because of the term nαε n+1 . We denote by ξ n+1 = [nα + β]ε n+1 0 and we now have the recursion:

Θ n+1 = F Θ n + ξ n+1 , (12)
which is a standard noisy linear dynamical system (see, e.g., Arnold, 1998) with uncorrelated noise process (ξ n ). We may thus express Θ n directly as Θ n = F n Θ 0 + n k=1 F n-k ξ k , and its expected second-order moment as,

E Θ n Θ n = F n Θ 0 Θ 0 (F n ) + n k=1 F n-k E ξ k ξ k (F n-k ) .
In order to obtain the expected excess cost function, we simply need to compute tr 0

H 0 0 E Θ n Θ n ,
which thus decomposes as a term that only depends on initial conditions (which is exactly the one computed and studied in Section 3.3), and a new term that depends on the noise.

Convergence result

For a quadratic function f with arbitrarily small eigenvalues and uncorrelated noise with finite covariance, we obtain the following convergence result (see proof in Appendix F); since we will allow the parameters α and β to depend on the time we stop the algorithm, we introduce the horizon N :

Theorem 3 (Convergence rates with noisy gradients).

With E[ε n ⊗ ε n ] = C for all n ∈ N, for α ≤ 1 L and 0 ≤ β ≤ 2 L -α.
Then for any N ∈ N, we have:

Ef (θ N ) -f (θ * ) ≤ min θ 0 -θ * 2 αN 2 + (αN + β) 2 αN tr(C), 4 θ 0 -θ * 2 (α + β)N + 4(αN + β) 2 α + β tr(C) . ( 13
)
We can make the following observations:

-Although we only provide an upper-bound, the proof technique relies on direct moment computations in each eigensubspace with few inequalities, and we conjecture that the scalings with respect to n are tight.

-For α = 0 and β = 1/L (which corresponds to averaged gradient descent), the second bound leads to 4L θ 0 -θ * 2 N + 4 tr(C) L , which is bounded but not converging to zero. We recover a result from Bach and Moulines (2011, Theorem 1).

-For α = β = 1/L (which corresponds to Nesterov's acceleration), the first bound leads

to L θ 0 -θ * 2 N 2 + (N +1) tr(C)

L

, and our bound suggests that the algorithm diverges, which we have observed in our experiments in Appendix A.

-For α = 0 and β = 1/L √ N , the second bound leads to

4L θ 0 -θ * 2 √ N + 4 tr(C) L √
N , and we recover the traditional rate of 1/ √ N for stochastic gradient in the non-strongly-convex case.

-When the values of the bias and the variance are known we can choose α and β such that the trade-off between the bias and the variance is optimal in our bound, as the following corrollary shows. Note that in the bound below, taking a non zero β enables the bias term to be adaptive to hidden strong-convexity.

Corollary 2. For α = min

θ 0 -θ * 2 √ tr CN 3/2 , 1/L and β ∈ [0, min{N α, 1/L}],
we have:

Ef (θ N ) -f (θ * ) ≤ 2L θ 0 -θ * 2 N 2 + 4 √ tr C θ 0 -θ * √ N .

Structured noise and least-square regression

When only the noise total variance tr(C) is considered, as shown in Section 4.4, Corollary 2 recover existing (more general) results. Our framework however leads to improved result for structured noise processes frequent in machine learning, in particular in least-squares regression which we now consider but this goes beyond (see, e.g. [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate O(1/n)[END_REF].

Assume we observe independent and identically distributed pairs (x n , y n ) ∈ R d × R and we want to minimize the expected loss f

(θ) = 1 2 E[(y n -θ, x n ) 2 ].
We denote by H = E(x n ⊗ x n ) the covariance matrix which is assumed invertible. The global minimum of f is attained at θ * ∈ R d defined as before and we denote by r n = y n -θ * , x n the statistical noise, which we assume bounded by σ. We have E[r n x n ] = 0. In an online setting, we observe the gradient (x n ⊗ x n )(θ -θ * ) -r n x n , whose expectation is the gradient f (θ). This corresponds to a noise in the gradient of

ε n = (H -x n ⊗ x n )(θ -θ * ) + r n x n .
Given θ, if the data (x n , y n ) are almost surely bounded, the covariance matrix of this noise is bounded by a constant times H. This suggests to characterize the noise convergence by tr(CH -1 ), which is bounded even though H has arbitrarily small eigenvalues. However, our result will not apply to stochastic gradient descent (SGD) for least-squares, because of the term (H -x n ⊗ x n )(θ -θ * ) which depends on θ, but to a "semi-stochastic" recursion where the noisy gradient is

H(θ -θ * ) -r n x n , with a noise process ε n = r n x n , which is such that E[ε n ⊗ ε n ] σ 2 H,
and has been used by [START_REF] Bach | Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Machine Learning[END_REF] and [START_REF] Dieuleveut | Non-parametric Stochastic Approximation with Large Step sizes[END_REF] to prove results on regular stochastic gradient descent. We conjecture that our algorithm (and results) also applies in the regular SGD case, and we provide encouraging experiments in Section 5.

For this particular structured noise we can take advantage of a large β:

Theorem 4 (Convergence rates with structured noisy gradients).

Let α ≤ 1 L and 0 ≤ β ≤ 3 2L -α 2 . For any N ∈ N, Ef (θ N ) -f (θ * ) is upper-bounded by: min θ 0 -θ * 2 N 2 α + (αN + β) 2 αβN 2 tr(CH -1 ), 4L θ 0 -θ * 2 (α + β)N + 8(αN + β) 2 tr(CH -1 ) (α + β) 2 N . ( 14
)
We can make the following observations:

-For α = 0 and β = 1/L (which corresponds to averaged gradient descent), the second bound leads to 4L θ 0 -θ * 2 N + 8 tr(CH -1 )

N

. We recover a result from Bach and Moulines (2013, Theorem 1). Note that when C σ 2 H, tr(CH -1 ) σ 2 d.

-For α = β = 1/L (which corresponds to Nesterov's acceleration), the first bound leads

to L θ 0 -θ * 2 N 2
+ tr(CH -1 ), which is bounded but not converging to zero (as opposed to the the unstructured noise where the algorithm may diverge).

-For α = 1/(LN a ) with 0 ≤ a ≤ 1 and β = 1/L, the first bound leads to

L θ 0 -θ * 2 N 2-a + tr(CH -1 ) N a
. We thus obtain an explicit bias-variance trade-off by changing the value of a.

-When the values of the bias and the variance are known we can choose α and β with an optimized trade-off, as the following corrollary shows:

Corollary 3. For α = min θ 0 -θ * √ L tr(CH -1 )N
, 1/L and β = min {N α, 1/L} we have:

Ef (θ N ) -f (θ * ) ≤ max 5 tr(CH -1 ) N , 5 tr(CH -1 )L θ 0 -θ * N , 2 θ 0 -θ * 2 L N 2 . ( 15
)

Related work

Acceleration and noisy gradients. Several authors [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF][START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF][START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF] have shown that using a step-size proportional to 1/N 3/2 accelerated methods with noisy gradients lead to the same convergence rate of

O L θ 0 -θ * 2 N 2 + θ 0 -θ * √ tr(C) √ N
than in Corollary 2, for smooth functions. Thus, for unstructured noise, our analysis provides insights in the behavior of second-order algorithms, without improving bounds. We get significant improvements for structured noises.

Least-squares regression. When the noise is structured as in least-square regression and more generally in linear supervised learning, [START_REF] Bach | Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Machine Learning[END_REF] have shown that using averaged stochastic gradient descent with constant step-size leads to the convergence rate of

O L θ 0 -θ 0 2 N + σ 2 d N .
It has been highlighted by [START_REF] Défossez | Constant step size least-mean-square: Bias-variance trade-offs and optimal sampling distributions[END_REF] that the bias term L θ 0 -θ * 2 N may often be the dominant one in practice. Our result in Corollary 3 leads to an improved bias term in O(1/N 2 ) with the price of a potentially slightly worse constant in the variance term. However, with optimal constants in Corollary 3, the new algorithm is always an improvement over averaged stochastic gradient descent in all situations. If constants are unknown, we may use α = 1/(LN a ) with 0 ≤ a ≤ 1 and β = 1/L and we choose a depending on the emphasis we want to put on bias or variance.

Minimax convergence rates. For noisy quadratic problems, the convergence rate nicely decomposes into two terms, a bias term which corresponds to the noiseless problem and the variance term which corresponds to a problem started at θ * . For each of these two terms, lower bounds are known. For the bias term, if N ≤ d, then the lower bound is, up to constants, L θ 0 -θ * 2 /N 2 (Nesterov, 2004, Theorem 2.1.7). For the variance term, for the general noisy gradient situation, we show in Appendix H that for N ≤ d, it is (tr C)/(L √ N ), while for least-squares regression, it is σ 2 d/N [START_REF] Tsybakov | Optimal rates of aggregation[END_REF]. Thus, for the two situations, we attain the two lower bounds simultaneously for situations where respectively L θ 0 -θ * 2 ≤ (tr C)/L and L θ 0 -θ * 2 ≤ dσ 2 . It remains an open problem to achieve the two minimax terms in all situations.

Other algorithms as special cases. We also note as shown in Appendix G that in the special case of quadratic functions, the algorithms of Lan (2012); [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF]; [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF] could be unified into our framework (although they have significantly different formulations and justifications in the smooth case).

Experiments

In this section, we illustrate our theoretical results on synthetic examples. We consider a matrix H that has random eigenvectors and eigenvalues 1/k m , for k = 1, . . . , d and m ∈ N. We take a random optimum θ * and a random starting point θ 0 such that r = θ 0 -θ * = 1 (unless otherwise specified). In Appendix A, we illustrate the noiseless results of Section 3, in particular the oscillatory behaviors and the influence of all eigenvalues, as well as unstructured noisy gradients. In this section, we focus on noisy gradients with structured noise (as described in Section 4.3), where our new algorithms show significant improvements.

We compare our algorithm to other stochastic accelerated algorithms, that is, AC-SA [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF], SAGE [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF] and Acc-RDA [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF] which are presented in Appendix G. For all these algorithms (and ours) we take the optimal step-sizes defined in these papers. We show results averaged over 10 replications.

Homoscedastic noise. We first consider an i.i.d. zero mean noise whose covariance matrix is proportional to H. We also consider a variant of our algorithm with an any-time step-size function of n rather than N (for which we currently have no proof of convergence).

In Figure 3, we take into account two different set-ups. In the left plot, the variance dominates the bias (with r = θ 0 -θ * = σ). We see that (a) Acc-GD does not converge to the optimum but does not diverge either, (b) Av-GD and our algorithms achieve the optimal rate of convergence of O(σ 2 d/n), whereas (c) other accelerated algorithms only converge at rate O(1/ √ n). In the right plot, the bias dominates the variance (r = 10 and σ = 0.1). In this situation our algorithm outperforms all others.

Application to least-squares regression. We now see how these algorithms behave for least-squares regressions and the regular (non-homoscedastic) stochastic gradients described in Section 4.3. We consider normally distributed inputs. The covariance matrix H is the same as before. The outputs are generated from a linear function with homoscedatic noise with a signal-to-noise ratio of σ. We consider d = 20. We show results averaged over 10 replications. In Figure 4, we consider again a situation where the bias dominates (left) and vice versa (right). We see that our algorithm has the same good behavior than in the homoscedastic noise case and we conjecture that our bounds also hold in this situation.

Conclusion

We have provided a joint analysis of averaging and acceleration for non-strongly-convex quadratic functions in a single framework, both with noiseless and noisy gradients. This allows to define a class of algorithms that can benefit simultaneously of the known improvements of averaging and accelerations: faster forgetting of initial conditions (for acceleration), and better robustness to noise when the noise covariance is proportional to the Hessian (for averaging).

Our current analysis of our class of algorithms in Eq. ( 4), that considers two different affine combinations of previous iterates (instead of one for traditional acceleration), is limited to quadratic functions; an extension of its analysis to all smooth or self-concordant-like functions would widen its applicability. Similarly, an extension to least-squares regression with natural heteroscedastic stochastic gradient, as suggested by our simulations, would be an interesting development.

A Additional experimental results

In this appendix, we provide additional experimental results to illustrate our theoretical results.

A.1 Deterministic convergence

Comparaison for d = 1. In Figure 5, we minimize a one-dimensional quadratic function f (θ) = 1 2 θ 2 for a fixed step-size α = 1/10 and different step-sizes β. In the left plot, we compare Acc-GD, HB and Av-GD. We see that HB and Acc-GD both oscillate and that Acc-GD leverages strong convexity to converge faster. In the right plot, we compare the behavior of the algorithm for different values of β. We see that the optimal rate is achieved for β = β * defined to be the one for which there is a double coalescent eigenvalue, where the convergence is linear at speed O(1 -√ αL) n . When β > β * , we are in the real case and when β < β * the algorithm oscillates to the solution.

Comparison between the different eigenspaces. Figure 6 shows interactions between different eigenspaces. In the left plot, we optimize a quadratic function of dimension d = 2. The first eigenvalue is L = 1 and the second is µ = 2 -8 . For Av-GD the convergence is of order O(1/n) since the problem is "not" strongly convex (i.e., not appearing as strongly convex since nµ remains small). The convergence is at the beginning the same for HB and Acc-GD, with oscillation at speed O(1/n 2 ), since the small eigenvalue prevents Acc-GD from having a linear convergence. Then for large n, the convergence becomes linear for Acc-GD, since µn becomes large. In the right plot, we optimize a quadratic function in dimension d = 5 with eigenvalues from 1 to 0.1. We show the function values of the projections of the iterates η n on the different eigenspaces. We see that high eigenvalues first dominate, but converge quickly to zero, whereas small ones keep oscillating, and converge more slowly. Comparison for d = 20. In Figure 7, we optimize two 20-dimensional quadratic functions with different eigenvalues with Av-GD, HB and Acc-GD for a fixed step-size γ = 1/10. In the left plot, the eigenvalues are 1/k 2 and in the right one, they are 1/k 8 , for k = 1, . . . , d. We see that in both cases, Av-GD converges at a rate of O(1/n) and HB at a rate of O(1/n 2 ). For Acc-GD the convergence is linear when µ is large (left plot) and becomes sublinear at a rate of O(1/n 2 ) when µ becomes small (right plot).

A.2 Noisy convergence with unstructured additive noise

We optimize the same quadratic function, but now with noisy gradients. We compare our algorithm to other stochastic accelerated algorithms, that is, AC-SA [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF], SAGE [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF] and Acc-RDA [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF], which are presented in Appendix G. For all these algorithms (and ours) we take the optimal step-sizes defined in these papers. We plot the results averaged over 10 replications.

We consider in Figure 8 an i.i.d. zero mean noise of variance C = I. We see that all the accelerated algorithms achieve the same precision whereas Av-GD with constant step-size does not converge and Acc-Gd diverges. However SAGE and AC-SA are anytime algorithms and are faster at the beginning since their step-sizes are decreasing and not a constant (with respect to n) function of the horizon N .

Let (P n , Q n , R n ) ∈ (R[X]
) 3 for all n ∈ N be a sequence of polynomials. We consider the iterates defined for all n ∈ N * by

θ n+1 = P n (H)θ n + Q n (H)θ n-1 + R(H)q, started from θ 0 = θ 1 ∈ R d .
The θ * -stationarity property gives for n ∈ N * :

θ * = P n (H)θ * + Q n (H)θ * + R n (H)q.
Since θ * = H -1 q we get for all q ∈ R d

H -1 q = P n (H)H -1 q + Q n (H)H -1 q + R n (H)q.
For all q ∈ R d we apply this relation to vectors q = H q: q = P n (H)q + Q n (H)q + R n (H)H q ∀q ∈ R d , and we get

I = P n (H) + Q n (H) + R n (H)H ∀n ∈ N * .
Therefore there are polynomials ( Pn , Qn ) ∈ (R[X]) 2 and q n ∈ R for all n ∈ N * such that we have for all n ∈ N:

P n (X) = (1 -q n )I + X Pn (X) Q n (X) = q n I + X Qn (X) R n (X) = -( Pn (X) + Qn (X)). ( 16
)
The n-scalability property means that there are polynomials (P, Q) ∈ (R[X]) 2 independent of n such that:

P n (X) = n n + 1 P (X), Q n (X) = n -1 n + 1 Q(X).
And in connection with Eq. ( 16) we can rewrite P and Q as:

P (X) = p + X P (X), Q(X) = q + X Q(X),
with (p, q) ∈ R 2 and ( P , Q) ∈ (R[X]) 2 . Thus for all n ∈ N:

q n = n -1 n + 1 q (17) Qn (X) = n -1 n Q(X) n n + 1 p = (1 -q n ) (18) Pn (X) = n n + 1 P (X).
Eq. ( 17) and Eq. ( 18) give:

n n + 1 p = 1 - n -1 n + 1 q .
Thus for n = 1, we have p = 2. Then -n-1 n+1 q = 2n n+1 -1 = n-1 n+1 and q = -1. Therefore

P n (H) = 2n n + 1 I + n n + 1 P (H)H Q n (H) = - n -1 n I + Q(H)H R n (H) = - n P (H) + (n -1) Q(H) n + 1 .
We let Ā = -( P + Q) and B = Q so that we have:

P n (H) = 2n n + 1 I - Ā(H) + B(H) 2 H Q n (H) = - n -1 n I -B(H)H R n (H) = n Ā(H) + B(H) n + 1 ,
and with φ n = θ n -θ * for all n ∈ N, the algorithm can be written under the form:

φ n+1 = I - n n + 1 Ā(H) + 1 n + 1 B(H) H φ n + 1 - 2 n + 1 I -B(H)H (φ n -φ n-1 ).

B.2 Av-GD as two steps-algorithm

We show now that when the averaged iterate of Av-GD is seen as the main iterate we have that Av-GD with step-size γ ∈ R is equivalent to:

θ n+1 = 2n n + 1 θ n - n -1 n + 1 θ n-1 - γ n + 1 f nθ n -(n -1)θ n-1 .
We remind

ψ n+1 = ψ n -γf (ψ n ), θ n+1 = θ n + 1 n + 1 (ψ n+1 -θ n ).
Thus, we have:

θ n+1 = θ n + 1 n + 1 (ψ n+1 -θ n ) = θ n + 1 n + 1 (ψ n -γf (ψ n ) -θ n ) = θ n + 1 n + 1 (θ n + (n -1)(θ n -θ n-1 ) -γf (θ n + (n -1)(θ n -θ n-1 )) -θ n ) = 2n n + 1 θ n - n -1 n + 1 θ n-1 - γ n + 1 f (nθ n -(n -1)θ n-1 ).

C Proof of Section 3 C.1 Proof of Lemma 1

The discriminant ∆ i is strictly positive when α+β 2 2 h i -α > 0. This is always true for α strictly negative. For α positive and for h i = 0, this is true for Then we determine when the modulus of the eigenvalues is less than one (which corresponds to -1 ≤ r - i ≤ r + i ≤ 1).

| α+β 2 | > α/h i . Thus the discriminant ∆ i is strictly positive for α < 0 or α ≥ 0 and β < -α -2 α/h i or β > -α + 2 α/h i . 20 0 1 2 3 4 0 1 2 βh i αh i r + i ≤ 1 r - i ≥ -1 ∆ i > 0 r + i = 1 r - i = -1 ∆ i = 0
r + i ≤ 1 ⇔ h i α + β 2 2 h i -α ≤ β + α 2 h i ⇔ h i β + α 2 2 h i -α ≤ β + α 2 h i 2 and α + β 2 ≥ 0 ⇔ h i α ≥ 0 and α + β 2 ≥ 0 ⇔ α ≥ 0 and α + β ≥ 0.
Moreover, we have :

r - i ≥ -1 ⇔ h i β + α 2 2 h i -α ≤ 2 - β + α 2 h i ⇔ h i β + α 2 2 h i -α ≤ 2 - β + α 2 h i 2 and 2 - β + α 2 h i ≥ 0 ⇔ h i β + α 2 2 h i -α ≤ 4 -4 β + α 2 h i + β + α 2 h i 2 and β + α 2 ≤ 2/h i ⇔ -h i α ≤ 4 -4 β + α 2 h i and β 2 ≤ 2/h i - α 2 ⇔ β ≤ 2/h i -α/2 and β ≤ 4/h i -α.
Figure 9 (where we plot all the constraints we have so far) enables to conclude that the discriminant ∆ i is strictly positive and the algorithm is stable when the following three conditions are satisfied:

α ≥ 0 α + 2β ≤ 4/h i α + β ≥ 2 α/h 1 .
For any of those α et β we will have:

η i n = c 1 (r - i ) n + c 2 (r + i ) n .
Since η i 0 = 0, c 1 + c 2 = 0 and for n = 1, c 1 = η i 1 /(r - i -r + i ); we thus have:

η i n = η i 1 2 (r + i ) n -(r - i ) n √ ∆ i .
Thus, we get the final expression: The discriminant ∆ i is strictly negative if and only if α+β 2 2 h i -α < 0. This implies

(φ i n ) 2 h i = (φ i 1 ) 2 4n 2 r i + √ ∆ i n -r i - √ ∆ i n 2 ∆ i /h i . C.2 Proof of Lemma 2 0 1 2 3 4 0 1 2 βh i αh i ∆ 1 < 0 |r ± i | ≤ 1 |r ± i | = 1 ∆ i = 0
| α+β 2 | < α/h i . The modulus of the eigenvalues is |r ± i | 2 = 1 -βh i .
Thus the discriminant ∆ i is strictly negative and the algorithm is stable for

α, β ≥ 0 α + β < α/h i ,
as shown in Figure 10.

For any of those α et β we have:

η i n = [c 1 cos(ω i n) + c 2 sin(ω i n)]ρ n i , with ρ i = √ 1 -βh i , sin(ω i ) = √ -∆ i /ρ i and cos(ω i ) = r i /ρ i . Since η i 0 = 0, c 1 = 0 and we have for n = 1, c 2 = η i 1 /(sin(ω i )ρ i ).
Therefore

η i n = η i 1 sin(ω i n) √ -∆ i (1 -βh i ) n/2 ,
and

(φ i n ) 2 h i = (φ i 1 ) 2 n 2 sin 2 (ω i n) sin 2 (ω i )/h i (1 -βh i ) n-1 .

C.3 Coalescing eigenvalues

When β = 2 α/h i -α, the discriminant ∆ i is equal to zero and we have a double real eigenvalue:

r i = 1 -αh i .
Thus the algorithm is stable for α < 4 h i . For any of those α et β we have:

η i n = (c 1 + nc 2 )r n .
This gives with η i 0 = 0, c 1 = 0 and c 2 = η i 1 /r. Therefore

η i n = nη i (1 -αh i ) n-1 ,
and:

(φ i n ) 2 h i = h i (φ i 1 ) 2 (1 -αh i ) 2(n-1) .
In the presence of coalescing eigenvalues the convergence is linear if 0 < α < 4/h i and h i > 0, however one might worry about the behavior of ((φ i n ) 2 h i ) n when h i becomes small. Using the bound x 2 exp(-x) ≤ 1 for x ≤ 1, we have for α < 4/h i :

h i (1 -αh i ) 2n = h i exp(2n log(|1 -αh i |)) ≤ hi exp(-2n min{ αh i , 2 -αh i }) ≤ h i min{ √ αh i , 2 - √ αh i } 2 ≤ max 1 α , h i (2 - √ αh i ) 2 .
Therefore we always have the following bound for α < 4/h i :

(φ i n ) 2 h i ≤ (φ i 1 ) 2 4n 2 max 1 α , h i (2 - √ αh i ) 2 .
Thus for αh i ≤ 1 we get: We divide the domain of validity of Theorem 2 in three subdomains as explained in Figure 14. On the domain described in Figure 11 we have a first bound on the iterate η i n : Lemma 3. For 0 ≤ α ≤ 1/h i and 1 -√ 1 -αh i < βh i < 1 + √ 1 -αh i , we have:

(φ i n ) 2 h i ≤ (φ i 1 ) 2 4n 2 α .
(η i n ) 2 ≤ (η i 1 ) 2 αh i .
And on the domain described Figure 12 we also have:

Lemma 4. For 0 ≤ α ≤ 1/h i and β ≤ α we have:

(η i n ) 2 ≤ 2(η i 1 ) 2 αh i .
These two lemmas enable us to prove the first bound of Theorem 2 since the domain of this theorem is included in the intersection of the two domains of these lemmas as shown in Figure 14.

Then we have the following bound on domain described in Figure 13:

Lemma 5. For 0 ≤ α ≤ 2/h i and 0 ≤ β ≤ 2/h i -α, we have:

|η i n | ≤ min 2 √ 2n (α + β)h i , 4 (α + β)h i .
Since the domain of definition of Theorem 2 is included in the domain of definition of Lemma 5 (as shown in Figure 14), this lemma proves the last two bounds of the theorem.

D.2 Outline of the proofs of the Lemmas

-We find a Lyapunov function G from R 2 to R such that the sequence (G(η i n , η i n-1 )) decrease along the iterates.

-We also prove that G(η i n , η i n-1 ) dominates c η i n 2 when we want to have a bound on

η i n 2 of the form 1 c G(η i 1 , η i 0 ) = 1 c G(θ i 0 -θ i * , 0).
For readability, we remove the index i and take h i = 1 without loss of generality.

D.3 Proof of Lemma 3

We first consider a quadratic Lyapunov function

η n η n-1 G 1 η n η n-1 with G 1 = 1 α -1 α -1 1 -α .
We note that G 1 is symmetric positive semi-definite for α ≤ 1. We recall

F i = 2 -(α + β) β -1 1 0 .
For the result to be true we need for 0 ≤ α ≤ 1 and 1 -

√ 1 -α < β < 1 + √ 1 -α two properties: F i G 1 F i G 1 , (19)
and

α 1 0 0 0 G 1 . (20)
Proof of Eq. ( 20). We have:

G 1 -α 1 0 0 0 = (1 -α) 1 1 1 1 0 for α ≤ 1. Proof of Eq. (19). Since β → F i (β) G 1 F i (β) -G 1 is convex in β (G 1 is symmetric positive semi-definite)
, we only have to show Eq. ( 19) for the boundaries of the interval in β. For x ∈ R * + :

x

2 -x x 1 0 1 -x 2 -x 2 x 2 x 2 -x x 1 0 - 1 -x 2 -x 2 x 2 = -(1 -x 2 ) 2 1 0 0 0 0.
This especially shows Eq. ( 19) for the boundaries of the interval with x = ± √ 1 -α.

Bound. Thus, because η 0 = 0, we have

αη 2 n+1 ≤ Θ n G 1 Θ n ≤ Θ n-1 G 1 Θ n-1 ≤ Θ 0 G 1 Θ 0 ≤ η 2 1 .
This shows that for 0

≤ α ≤ 1/h i and 1 - √ 1 -αh i < βh i < 1 + √ 1 -αh i : (η i n ) 2 ≤ (η i 1 ) 2 αh i .

D.4 Proof of Lemma 4

We consider now a second Lyapunov function

G 2 (η n , η n-1 ) = (η n -rη n-1 ) 2 -∆(η n-1 ) 2 .
We have:

G 2 (η n , η n-1 ) = (η n -rη n-1 ) 2 -∆η 2 n-1 = (rη n-1 -(1 -β)η n-2 ) 2 -∆η 2 n-1 = (r 2 -∆)η 2 n-1 + (1 -β) 2 η 2 n-2 -2(1 -β)rη n-1 η n-2 = ((1 -β)η 2 n-1 + (1 -β)(r 2 -∆)η 2 n-2 -2(1 -β)rη n-1 η n-2 = (1 -β)[(η n-1 -rη n-2 ) 2 -∆(η n-2 ) 2 ]. = (1 -β)G 2 (η n-1 , η n-2 ).
Where we have used twice

r 2 -∆ = (1 -β) and η n = 2rη n-1 -(1 -β)η n-2 . Moreover G 2 (η n , η n-1
) can be rewritten as:

G 2 (η n , η n-1 ) = (1 - α + β 2 )(η n -η n-1 ) 2 + α -β 2 (η n-1 ) 2 + α + β 2 (η n ) 2 .
Thus for α + β ≤ 2 and β ≤ α we have:

α 2 (η n ) 2 ≤ G 2 (η n , η n-1 ) = (1 -β) n-1 G 2 (η 1 , η 0 ) = (1 -β) n-1 (η 1 ) 2 .
Therefore for α + β ≤ 2/h i and β ≤ α, we have:

(η i n ) 2 ≤ 2(η i 1 ) 2 αh i .

D.5 Proof of Lemma 5

We may write η n as

η n = rη n-1 + (r + ) n + (r -) n .
Moreover, we have:

|(r + ) n + (r -) n | ≤ 2, therefore for α + β ≤ 2, |η n | ≤ r|η n-1 | + 2 ≤ 2 1 -r n 1 -r ≤ 2 1 -(1 -( α+β 2 )) n ( α+β 2 )
.

Thus

|η n | ≤ 2 ( α+β 2 )h
.

Moreover for all u ∈ [0, 1] and n ≥ 1 we have 1

-(1 -u) n ≤ √ nu, since 1 -(1 -u) n ≤ 1 and 1 -(1 -u) n = u (1 -u) k ≤ nu. Thus |η n | ≤ 2 √ n ( α+β 2 )
.

Therefore for 0 ≤ α ≤ 2/h i and α + β ≤ 2/h i we have:

|η i n | ≤ min 2 √ 2n (α + β)h i , 4 (α + β)h i .

E Lower bound

We have the following lower-bound for the bound shown in Corollary 1, which shows that depending on which of the two terms dominates, we may always find a sequence of functions that makes it tight.

Proposition 1. Let L ≥ 0. For all sequences 0 ≤ α n ≤ 1/L and 0 ≤ β n ≤ 2/L -α n , such that α n + β n = o(nα n )
there exists a sequence of one-dimensional quadratic functions (f n ) n with second-derivative less than L such that:

lim α n n 2 (f n (θ n ) -f n (θ * )) = θ 0 -θ * 2 2 .
For all sequences 0 ≤ α n ≤ 1/L and 0 ≤ β n ≤ 2/L -α n , such that nα n = o(α n + β n ), there exists a sequence of one-dimensional quadratic functions (g n ) n with second-derivative less than L such that:

lim n(α n + β n )(g n (θ n ) -g n (θ * )) = (1 -exp(-2)) 2 θ 0 -θ * 2 4 .
Proof of the first lower-bound. For the first lower bound we consider 0 ≤ α n ≤ 1/L and 0 ≤ β n ≤ 2/L -α, such that α n + β n = o(nα n ). We define f n = π 2 /(4α n n 2 ) and we consider the sequence of quadratic functions f n (θ) = fnθ 2 2 . We consider the iterate (η n ) n defined by our algorithm. We will show that

lim α n f n (η n ) = η 2 1 2 .
We have, from Lemma 2,

f n (η n ) = η 2 n f n 2 = η 2 1 sin 2 (ω n n)ρ 2n n 2α n (1 -π 2 (αn+βn) 2 (4αnn) 2 )
.

Moreover,

ρ 2n n = 1 - β n π 2 4α n n 2 n = exp n log 1 - β n π 2 4α n n 2 = 1 + o(1), since βn αnn = o(1). Also, 1 -π 2 (αn+βn) 2 (4αnn) 2 = 1 + o(1), since α n + β n = o(nα n ). Moreover sin(ω n ) = √ -∆ n ρ n = √ f n α n -(αn+βn) 2 4 f n √ 1 -β n f n = π/(2n) + o(1/n), thus ω n = π/(2n) + o(1/n) and sin(nω n ) = 1 + o(1).
Proof of the second lower-bound. We consider now the situation where the second bound is active. Thus we take sequences (α n ) and (β n ), such that nα n = o(α n + β n ). We define g n = 2 n(αn+βn) + 4αn (αn+βn) 2 and consider the sequence of quadratic functions g n (θ) = gnθ 2

2 . We will show for the iterate (η n ) defined by our algorithm that:

lim n(α n + β n )(g n (θ n ) -g n (θ * )) = (1 -exp(-2)) 2 θ 0 -θ * 2 4 .
We will use Lemma 1. We first have

∆ n = α n + β n 2 2 g 2 n -α n g n = g n α n + β n 2 1 n . Thus (n∆ n )/g n = αn+βn 2 and ∆ n = 1 n 2 + 2α n n(α n + β n ) = 1 n 1 + 2α n n α n + β n = 1 n + α n α n + β n + o α n α n + β n . Moreover r n = 1 - α n + β n 2 g n = 1 - 1 n - 2α n α n + β n . Thus r + = 1 - α n α n+ β n + o α n α n + β n , and
r n + = exp(n log(r + )) = exp - nα n α n + β n + o nα n α n + β n = 1 + o(1). Furthermore r -= 1 - 2 n - 3α n α n+ β n + o α n α n + β n , and
r n -= exp(n log(r + )) = exp -2 - 3α n n α n + β n + o nα n α n + β n = exp(-2) + o(1)
.

Thus (r n + -r n -) 2 = (1 -exp(-2)) 2 + o(1)
. Finally, we have:

(α n + β n )n[g n (θ n ) -g n (θ * )] = α n + β n 2n θ 0 -θ * 2 [r n + -r n -] 2 4∆ n /g n = θ 0 -θ * 2 4 [r n + -r n -] 2 = θ 0 -θ * 2 4 (1 -exp(-2)) 2 + o(1).
F Proofs of Section 4

F.1 Proofs of Theorem 3 and Theorem 4

We decompose again vectors in an eigenvector basis of H with η i n = p i η n and ε i n = p i ε n :

η i n+1 = (1 -αh i )η i n + (1 -βh i )(η i n -η i n-1 ) + (nα + β)ε i n+1 .
We denote by ξ i n+1 =

[nα + β]ε i n+1 0 and we have the reduced equation:

Θ i n+1 = F i Θ i n + ξ i n+1 .
Unfortunately F i is not Hermitian and this formulation will not be convenient for calculus. Without loss of generality, we assume r - i = r + i even if it means having r - i -r + i goes to 0 in the final bound. Let

Q i = r - i r + i 1 1 be the transfer matrix of F i , i.e., F i = Q i D i Q -1 i with D i = r - i 0 0 r + i and Q -1 i = 1 r - i -r + i 1 -r + i -1 r - i
. We can reparametrize the problem in the following way:

Q -1 i Θ i n+1 = Q -1 i F i Θ i n + Q -1 i ξ i n+1 = Q -1 i F i Q i Q -1 i Θ i n + Q -1 i ξ i n+1 = D i (Q -1 i Θ i n ) + Q -1 i ξ i n+1 . With Θi n = Q -1 i Θ i n and ξi n = Q -1 i ξ i n we now have: Θi n+1 = D i Θi n + ξi n+1 , (21)
with now D i Hermitian (even diagonal).

Thus it is easier to tackle using standard techniques for stochastic approximation (see, e.g., [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Bach | Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Machine Learning[END_REF]:

Θi n = D n i Θi 0 + n k=1 D n-k i ξi k . Let M i = h 1/2 i h 1/2 i 0 0
, we then get using standard martingale square moment inequalities, since for n = m, ε i n and ε i m are uncorrelated (i.e., E[ε i n ε i m ] = 0):

E M i Θi n 2 = M i D n i Θi 0 2 + E n k=1 M i D n-k i ξi k 2 .
This is a bias-variance decomposition; the left term only depends on the initial condition and the right term only depends on the noise process.

We have with

M i = h 1/2 i h 1/2 i 0 0 , M i Q -1 i = 0 h 1/2 i 0 0 , and
M i Θi n = √ h i η i n 0
. Thus, we have access to the function values through:

M i Θi n 2 = h i (η i n ) 2 .

Moreover we have Θ

i 0 = φ i 1 /(r - i -r + i ) -φ i 1 /(r - i -r + i )
. Thus

M i D n i Θi 0 2 = (φ i 1 ) 2 h i (r + i ) n -(r - i ) n 2 (r + i -r - i ) 2
. This is the bias term we have studied in Section 3.3 which we bound with Theorem 2. The variance term is controlled by the next proposition.

Proposition 2. With E[(ε i n ) 2 ] = c i for all n ∈ N, for α ≤ 1/h i and 0 ≤ β ≤ 2/h i -α , we have 1 n 2 E n k=1 M i D n-k i ξi k 2 ≤ min 2(αn + β) 2 αβ(4 -(α + 2β)h i )n 2 c i h i , 16(nα + β) 2 n(α + β) 2 c i h i , 2 (αn + β) 2 nα c i , 8(nα + β) 2 α + β c i .
The last two bounds prove Theorem 3.

We note that if we restrict β to β ≤ 3/(2h i ) -α/2, then 4 -(α + 2β)h i ≥ 1 and the first bound of Proposition 2 is simplified to 2(αn+β) 2 αβn 2 c i h i . This allows to conclude to prove Theorem 4.

F.2 Proof of Corollary 3

We

let ν = θ 0 -θ * √ L tr(CH -1 )
and consider three different regimes depending on ν and L.

If ν < 1/L, we have ν/N < 1/L and thus α = ν/N and β = ν. Therefore

θ 0 -θ * 2 N 2 α + (αN + β) 2 αβN 2 tr(CH -1 ) = θ 0 -θ * 2 νN + 4 tr(CH -1 ) N ≤ L tr(CH -1 ) θ 0 -θ * N + 4 tr(CH -1 ) N ≤ 5 tr(CH -1 ) N ,
where we have used √ L θ 0 -θ * < tr(CH -1 ) since ν < 1/L.

If ν > 1/L and ν < N/L, we have α = ν/N and β = 1/L. Therefore

θ 0 -θ * 2 N 2 α + (αN + β) 2 αβN 2 tr(CH -1 ) ≤ θ 0 -θ * 2 νN + 4 tr(CH -1 ) LνN ≤ L tr(CH -1 ) θ 0 -θ * N + 4 tr(CH -1 ) N ≤ 5 L tr(CH -1 ) θ 0 -θ * N ,
where we have used

√ L θ 0 -θ * > tr(CH -1 ) since ν > 1/L.
If ν > N/L, we have α = 1/L and β = 1/L. Therefore

θ 0 -θ * 2 N 2 α + (α(N -1) + β) 2 αβN 2 tr(CH -1 ) = L θ 0 -θ * 2 N 2 + tr(CH -1 ) ≤ L θ 0 -θ * 2 N 2 + L θ 0 -θ * 2 N 2 ≤ 2L θ 0 -θ * 2 N 2 ,
where we have used that the real bound in Proposition 2 is in fact in (N -1)α + β, (see Lemma 6) and that tr(CH -1 ) <

L θ 0 -θ * 2 N 2 since ν > N/L.
F.3 Proof of Proposition 2

F.3.1 Proof outline

To prove Proposition 2 we will use Lemmas 6, 7 and 8, that are stated and proved in Section F.3.2.

We want to bound E

[ n k=1 M i D n-k i ξi k 2
] and according to Lemma 6, we have an explicit expansion using the roots of the characteristic polynomial:

E M i D k i ξi k 2 = h i ((k -1)α + β) 2 E[(ε i ) 2 ] [(r - i ) n-k -(r + i ) n-k ] 2 (r - i -r + i ) 2
. Thus, by bouding (k -1)α + β by (n -1)α + β, we get

E n k=1 M i D n-k i ξi k 2 ≤ h i ((n -1)α + β) 2 E[ε i 2 ] n k=1 [(r - i ) n-k -(r + i ) n-k ] 2 (r - i -r + i ) 2 . ( 22
)
Then, we have from Lemma 7 the inequality:

n k=1 [(r - i ) k -(r + i ) k ] 2 [(r - i ) -(r + i )] 2 ≤ 2 -βh i 4αβh 2 i (1 -( 1 4 α + 1 2 β)h i )
.

Therefore

E n k=1 M 1/2 i D n-k i ξi k 2 ≤ E[ε i 2 ] h i ((n -1)α + β) 2 4αβ 2 -βh i (1 -( 1 4 α + 1 2 β)h i )
.

This allows to prove the first part of the bound. The other parts are much simpler and are done in Lemma 8. Thus, adding these bounds gives for α ≤ 1/h i and 0 ≤ β ≤ 2/h i -α:

1 n 2 E n k=1 M i D n-k i ξi k 2 ≤ min 2(α(n -1) + β) 2 αβn 2 (4 -(α + 2β)h i ) c h i , 16((n -1)α + β) 2 n(α + β) 2 c h i , 2 (α(n -1) + β) 2 nα c i , 8((n -1)α + β) 2 α + β c i .

F.3.2 Some technical Lemmas

We first compute an explicit expansion of the noise term as a function of the eigenvalues of the dynamical system.

Lemma 6. For all α ≤ 1/h i and 0 ≤ β ≤ 2/h i -α we have

E M i D k i ξi k 2 = h i ((k -1)α + β) 2 E[(ε i ) 2 ] [(r - i ) n-k -(r + i ) n-k ] 2 (r - i -r + i ) 2 .
Proof. We first turn the Euclidean norm into a trace, using that tr[AB] = tr[BA] for two matrices A and B and that tr[x] = x for a real x.

E M i D n-k i ξi k 2 = TrD n-k i M i M i D n-k i E[ ξi k ( ξi k ) ], (23)
This enables us to separate the noise term from the rest of the formula. Then we compute the latter from the definition of ξi k in Eq. ( 21) :

E[ ξi k ( ξi k ) ] = ((k -1)α + β) 2 (r - i -r + i ) 2 E[(ε i ) 2 ] 1 -1 -1 1 .
And the first part of Eq. ( 23) is equal to:

D n-k i M i M i D n-k i = h i (r - i ) 2(n-k) (r - i ) (n-k) -(r + i ) (n-k) (r - i ) (n-k) -(r + i ) (n-k) (r + i ) 2(n-k) , because D i = r - i 0 0 r + i and M i = h 1/2 i h 1/2 i 0 0 . Therefore: E M i D n-k i ξi k 2 = h i ((k -1)α + β) 2 (r - i -r + i ) 2 E[ε i 2 ][(r - i ) n-k -(r + i ) n-k ] 2 .
In the following leamma, we bound a certain sum of powers of the roots.

Lemma 7. For all α ≤ 1/h i and 0

≤ β ≤ 2/h i -α we have n k=1 (r - i ) k -(r + i ) k ] 2 [(r - i ) -(r + i ) 2 ≤ 2 -βh i 4αβh 2 i (1 -( 1 4 α + 1 2 β)h i )
.

We first note that when the two roots become close, the denominator and the numerator will go to zero, which prevents from bounding the numerator easily. We also note that this bound is very tight since the difference between the two terms goes to zero when n goes to infinity.

Proof. We first expand the square of the difference of the powers of the roots and compute their sums.

n k=1 (r - i ) k -(r + i ) k 2 = n k=1 r + i 2k + r - i 2k -2(r + i r - i ) k = 1 -r + i 2n 1 -r + i 2 + 1 -r - i 2n 1 -r - i 2 -2 1 -(r + i r - i ) n 1 -(r + i r - i ) = 1 1 -r + i 2 + 1 1 -r - i 2 - 2 1 -(r + i r - i ) - r + i 2n 1 -r + i 2 + r - i 2n 1 -r - i 2 -2 (r + i r - i ) n 1 -(r + i r - i ) = 1 1 -r + i 2 + 1 1 -r - i 2 - 2 1 -(r + i r - i ) -I n , with I n = r + i 2n 1-r + i 2 + r - i 2n 1-r - i 2 -2 (r + i r - i ) n 1-(r + i r - i )
. This sum is therefore equal to the sum of one term we will compute explicitly and one other term which will go to zero. We have for the first term:

1 1 -r + i 2 + 1 1 -r - i 2 - 2 1 -(r + i r - i ) = (1 -r - i 2 )(1 -(r + i r - i )) -(1 -r - i 2 )(1 -r + i 2 ) (1 -r + i 2 )(1 -r - i 2 )(1 -(r + i r - i )) + (1 -r + i 2 )(1 -(r + i r - i )) -(1 -r - i 2 )(1 -r + i 2 ) (1 -r + i 2 )(1 -r - i 2 )(1 -(r + i r - i )) , with (1 -r - i 2 )(1 -(r + i r - i )) -(1 -r - i 2 )(1 -r + i 2 ) = (1 -r - i 2 )[(1 -(r + i r - i )) -(1 -r + i 2 )] = r + i (1 -r - i 2 )(r + i -r - i ),
and

(1 -r + i 2 )(1 -(r + i r - i )) -(1 -r - i 2 )(1 -r + i 2 ) = -r - i (1 -r - i
2 )(r + i -r - i ), and

r + i (1 -r - i 2 )(r + i -r - i ) -r - i (1 -r - i 2 )(r + i -r - i ) = (r + i -r - i )[r + i (1 -r - i 2 ) -r - i (1 -r - i 2 )] = (r + i -r - i )[r + i -r - i + r + i r - i (r + i -r - i )] = (r + i -r - i ) 2 [1 + r + i r - i ].
Therefore the first term is equal to:

1 1 -r + i 2 + 1 1 -r - i 2 - 2 1 -(r + i r - i ) = (r + i -r - i ) 2 [1 + r + i r - i ] (1 -r + i 2 )(1 -r - i 2 )(1 -(r + i r - i ))
, and the sum can be expanded as:

n k=1 [(r - i ) k -(r + i ) k ] 2 [r - i -r + i ] 2 = [1 + r + i r - i ] (1 -r + i 2 )(1 -r - i 2 )(1 -(r + i r - i )) -J n , with J n = In [(r - i )-(r + i )]
2 . Then we simplify the first term of this sum using the explicit values of the roots. We recall

r ± i = r i ± √ ∆ i = 1 -α+β 2 h i ± α+β 2 2 h 2 i -αh i , therefore r + i r - i = r 2 i -∆ 2 i = 1 - α + β 2 h i 2 - α + β 2 h i 2 + αh i = 1 -βh i ,
and

(1 -r + i 2 )(1 -r - i 2 ) = [(1 -r - i )(1 -r + i )][(1 + r + i )(1 + r + i )] = [(1 -r i + ∆ i )((1 -r i -∆ i )][(1 + r i + ∆ i )((1 + r i -∆ i )] = [(1 -r i ) 2 -∆ i ][(1 + r i ) 2 -∆ i ][(1 -r i ) 2 -∆ i ] = 4αh i 1 - 1 4 α + 1 2 β h i . Thus n k=1 [(r - i ) k -(r + i ) k ] 2 [(r - i ) -(r + i )] 2 = 2 -βh i 4αβh 2 i (1 -( 1 4 α + 1 2 β)h i ) -J n .
Even if J n will be asymptotically small, we want a non-asymptotic bound, thus we will show that J n is always positive.

In the real case [(r - i ) -(r + i )] 2 ≥ 0 and using a 2 + b 2 ≥ 2ab, for all (a, b) ∈ R 2 , we have

r + i 2n 1 -r + i 2 + r - i 2n 1 -r - i 2 ≥ 2 (r + i r - i ) n (1 -r + i 2 )(1 -r - i 2 )
, and using r

+ i 2 + r - i 2 ≥ 2r - i r - i we have (1 -r + i 2 )(1 -r - i 2 ) ≤ 1 -(r + i r - i ), since (1 -r + i 2 )(1 -r - i 2 ) -[1 -(r + i r - i )] 2 = 1 -r + i 2 -r - i 2 + (r + i r - i ) 2 -1 + 2r + i r - i -(r + i r - i ) 2 = 2r + i r - i -r + i 2 -r - i 2 ≤ 0. Thus r + i 2n 1 -r + i 2 + r - i 2n 1 -r - i 2 -2 (r + i r - i ) n 1 -(r + i r - i ) ≥ 0.
and J n ≥ 0 in the real case.

In the complex case, [(r - i ) -(r + i )] 2 ≤ 0, and using z 2 + z2 ≤ 2z z for all z ∈ C, we have

r + i 2n 1 -r + i 2 + r - i 2n 1 -r - i 2 ≤ 2 (r + i r - i ) n (1 -r + i 2 )(1 -r - i 2 )
, and using r

+ i 2 + r - i 2 ≤ 2r - i r - i we have (1 -r + i 2 )(1 -r - i 2 ) ≥ 1 -(r + i r - i ). Thus r + i 2n 1 -r + i 2 + r - i 2n 1 -r - i 2 -2 (r + i r - i ) n 1 -(r + i r - i ) ≤ 0.
and J n ≥ 0 in the complexe case.

Therefore we always have:

J n ≥ 0,
and

n k=1 [(r - i ) k -(r + i ) k ] 2 [(r - i ) -(r + i )] 2 ≤ 2 -βh i 4αβh 2 i (1 -( 1 4 α + 1 2 β)h i )
.

their notations. The only difference with their Theorem 1 in the different choice of two functions f + i and f - i , which we choose to be:

f ± i (x) = c i (x i ± r 2 ) 2 ,
with a non-increasing sequence (c i ) to be chosen later. The function g α that is optimized is thus:

g α (x) = 1 d d i=1 ( 1 2 + α i δ)f + i (x) + ( 1 2 -α i δ)f - i (x) .
This function is quadratic and its Hessian has eigenvalues equal to 2c i /d. Thus, its largest eigenvalue is 2c 1 /d, which we choose equal to L.

Noisy gradients are obtained by sampling d independent Bernoulli random variables b i , i = 1, . . . , d, with parameters ( 1 2 + α i δ) and using the gradient of the random function

1 d d i=1 b i f + i (x) + (1 -b i )f - i (x)
. The variance of the random gradient is equal to

V = d i=1 1 d 2 var b i c i (x i + r/2) -c i (x i -r/2) = 1 d 2 d i=1 c 2 i r 2 (1/4 -δ 2 ).
The function g α is minimized for x = -αδr, and the discrepancy measure between two functions g α and g β is greater than

1 d d i=1 inf x f + i (x)+f - i (x) -inf x f + i (x)-inf x f - i (x) 1 α i =β i 1 d d i=1 3c i r 2 δ 2 4 1 α i =β i 1 d 3c d r 2 δ 2 4 ∆(α, β).
Since the vectors α, β ∈ {-1, 1} d are so that their Hamming distance ∆(α, β) d/4 for α = β, we have a discrepancy measure greater than 3c d r 2 δ 2 16 . Thus, for a an approximate optimality of ε = c d r 2 δ 2 38 , we have, following the proof of Theorem 1 (equation ( 29)) from Agarwal et al. (2012), for N iterations of any method that accesses a random gradient, we have:

1/3 1 -2 16N dδ 2 + log 2 d log(2/ √ e) .
Thus, for d large, we get, up to constants, δ 2 1/N and thus ε r 2 c d N . For c 1 = 2Ld and c i = L √ d for the remaining ones, we get (up to constants):

ε V L √ d N .
This leads to the desired result for N d.

Figure 1 :

 1 Figure1: Area of stability of the algorithm, with the three traditional algorithms represented. In the interior of the triangle, the convergence is linear.

Figure 2 :

 2 Figure 2: Trade-off between averaged and accelerated methods for noisy gradients.

Figure 3 :

 3 Figure 3: Quadratic optimization with regression noise. Left σ = 1, r = 1. Right σ = 0.1, r = 10.

Figure 4 :

 4 Figure 4: Least-Square Regression. Left σ = 1, r = 1. Right σ = 0.1, r = 10.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Deterministic case for d = 1 and α = 1/10. Left: classical algorithms, right: different oscillatory behaviors.

  Figure 8: Quadratic optimization with additive noise.

Figure 9 :

 9 Figure 9: Stability in the real case, with all contraints plotted.

Figure 10 :

 10 Figure 10: Stability in the complex case, with all constraints plotted.

Figure 14 :

 14 Figure 11: Validity of Lemma 3

Acknowledgements

This work was partially supported by the MSR-Inria Joint Centre and a grant by the European Research Council (SIERRA project 239993). The authors would like to thank Aymeric Dieuleveut for helpful discussions.

A. Agarwal, P. L. Bartlett, P. Ravikumar, and M. J. Wainwright. Information-theoretic lower bounds on the oracle complexity of stochastic convex optimization. Information Theory, IEEE Transactions on, 58(5):3235-3249, 2012.

L. Arnold. Random dynamical systems. Springer Monographs in Mathematics. Springer-Verlag, 1998.

However we can also bound roughly Eq. ( 22) using Theorem 2 since we recall we have

. This gives us the following lemma which enables to prove the second part of Proposition 2.

Lemma 8. For all α ≤ 1/h i and 0 ≤ β ≤ 2/h i -α we have

Proof. From Lemma 6, we get

G Comparison with additional other algorithms G.1 Summary

When the objective function f is quadratic and for correct choices of step-sizes, the AC-SA algorithm of [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF], the SAGE algorithm of [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF] and the Accelerated RDA algorithm of [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF] are all equivalent to:

where we use H n θ + ε n as an unbiased estimate of the gradient and δ n as step-size which values will be specified later. [START_REF] Lan | An optimal method for stochastic composite optimization[END_REF] and [START_REF] Hu | Accelerated gradient methods for stochastic optimization and online learning[END_REF] only consider bounded cases by projecting their iterates on a bounded space. Xiao (2010) deals with the unbounded case and prove the following convergence result:

Theorem 5. (Xiao, 2010, Theorem 6). With

This result is significantly more general than ours since it is valid for composite optimization and general noise on the gradients.

We now present the different algorithms and show they all share the same form.

G.2 AC-SA

Lemma 9. AC-SA algorithm with step size γ n and β n and gradient estimate H n+1 θ n +ε n+1 is equivalent to:

Proof. We recall the general AC-SA algorithm:

• Let the initial points x ag 1 = x 1 , and the step-sizes {β n } n≤1 and {γ n } n≤1 be given. Set n = 1

• Step 3. Set n → n + 1 and go to step 1.

When f is quadratic we will have G(

n + γ n ε n+1 , and:

These give the result for θ n = x ag n .

G.3 SAGE

Lemma 10. The algorithm SAGE with step-sizes L n and α n is equivalent to:

Proof. We recall the general SAGE algorithm:

• Let the initial points x 0 = z 0 = 0, and the step-sizes {β n } n≤1 and {L n } n≤1 be given.

Set n = 1

• Step 3. Set n → n + 1 and go to step 1.

We have

and

Thus

These give the result for θ n = y n .

G.4 Accelerated RDA method

Lemma 11. The algorithm AccRDA with step-sizes β and α n is equivalent to:

Proof. We recall the general Accelerated RDA method:

• Let the initial points w 0 = v 0 , A 0 = 0, g0 = 0 and the step-sizes {α n } n≤1 and {β n } n≤1 be given.

Set n = 1

• Step 1. Set A n = A n-1 + α n and θ n = αn An .

• Step 2. Compute the query point

and update the weighted average gn gn = (1 -θ n )g n-1 + θ n g n .

•

Step 4. Set v n = v 0 -An L+βn gn .

• Step 6. Set n → n + 1 and go to step 1.

First we have

With β n = β we have v n = v n-1 -αn L+β (H n+1 u n + ε n+1 )] and

Since v n-1 = θ -1 n-1 w n-1 -θ -1 n-1 (1 -θ n-1 )w n-2 , then

and

H Lower bound for stochastic optimization for least-squares

In this section, we show a lower bound for optimization of quadratic functions with noisy access to gradients. We follow very closely the framework of Agarwal et al. (2012) and use


None; M A Muselier

M Souchier

B Mathieu

S Grégoire

F Nicot

F Baudouin

L Bretillon

A M Bron

C Creuzot-Garcher

None ; B Souchier

None; S Mathieu

None; F Grégoire

None; F Nicot

None; L Baudouin

None A M Bretillon

None; C Bron

None. Creuzot-Garcher

Ocular Surface Response in Contact Lens Wearers

Keywords: contact lens, cornea: tears/tear film/dry eye, lipids

Purpose:To evaluate the ocular surface changes induced by contact lens wear.

Methods:Data of 7 tolerant contact lens wearers, 7 non-tolerant contact lens wearers and 7 healthy subjects were collected. Every patient underwent a thorough ophthalmic examination and tear osmolarity test, conjunctival impression cytology and meibomian lipid sampling. Symptoms, Break-Up Time (BUT), tear osmolarity, conjunctival expression of human leucocyte antigen (HLA) DR and meibomian fatty acid composition were evaluated.

Results:Twenty one subjects (8 men and 13 women) were included (42 eyes). Symptoms reported by tolerant and non-tolerant contact lens wearers were as follows: superficial foreign body sensation (29% vs 21%), burning (14% vs 21%) and ocular dryness (28% vs 78%). No symptoms were reported by healthy subjects. There was no significant difference in tear osmolarity (p=0.63) and meibomian fatty acid composition (p=0.97) in both groups. Break-up time was significantly lower in tolerant (p<0.05) and non-tolerant (p<0.0001) contact lens wearers than in healthy subjects. The HLA-DR expression was significantly higher in non-tolerant (50%) contact lens wearers than in healthy subjects (16%) (p<0.01).

Conclusions:Contact lens wear is responsible for ocular surface alterations which lead to intolerance. These patterns are very similar to those reported in dry eye symptoms.