Antoine Godichon-Baggioni
email: antoine.godichon_baggioni@upmc.fr

Nicklas Werge
email: nicklas.werge@upmc.fr

Olivier Wintenberger
email: olivier.wintenberger@upmc.fr

Non-Asymptotic Analysis of Stochastic Approximation Algorithms for Streaming Data

Keywords: stochastic optimization, machine learning, stochastic algorithms, online learning, streaming, mini-batch

We consider the stochastic approximation problem in a streaming framework where an objective is minimized through unbiased estimates of its gradients. In this streaming framework, we consider time-varying data streams that must be processed sequentially. Our methods are Stochastic Gradient (SG) based due to their applicability and computational advantages. We provide a non-asymptotic analysis of the convergence of various SG-based methods; this includes the famous SG descent (a.k.a. Robbins-Monro algorithm), constant and time-varying mini-batch SG methods, and their averaged estimates (a.k.a. Polyak-Ruppert averaging). Our analysis suggests choosing the learning rate according to the expected data streams, which can speed up the convergence. In addition, we show how the averaged estimate can achieve optimal convergence in terms of attaining Cramer-Rao's lower bound while being robust to any data stream rate. In particular, our analysis shows how Polyak-Ruppert averaging of time-varying mini-batches can provide variance reduction and accelerate convergence simultaneously, which is advantageous for large-scale learning problems. These theoretical results are illustrated for various data streams, showing the effectiveness of the proposed algorithms.

Introduction

Machine learning and artificial intelligence have become an integral part of modern society. This massive utilization of intelligent systems generates an endless sequence of data, many of which come as streaming data such as weather, traffic, stock trade, or other real-time sensor data. These continuously generated data should be processed sequentially with the property that the data stream may change over time. Such a streaming framework requires computationally efficient and robust algorithms that can quickly update the model as more data arrives.

Stochastic approximation algorithms have proven effective in handling large amounts of data; Bottou et al. [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF] reviews such stochastic algorithms for large-scale machine learning, including noise reduction and second-order methods, among others. Among these, the most well-known is presumably the Stochastic Gradient (SG) descent proposed by [START_REF] Robbins | A stochastic approximation method[END_REF], which is used for many models within machine learning. Since its introduction, much work has been spent on analyzing, developing, and improving various SG-based methods, e.g., see Kushner and Yin [START_REF] Kushner | approximation and recursive algorithms and applications[END_REF], Lan [START_REF] Lan | First-order and stochastic optimization methods for machine learning[END_REF], Shalev-Shwartz et al. [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF]. An essential extension is the Polyak-Ruppert procedure (ASG) proposed by Polyak and Juditsky [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], Ruppert [START_REF] Ruppert | Efficient estimations from a slowly convergent Robbins-Monro process[END_REF], which guarantees optimal statistical efficiency without jeopardizing the computational complexity; this average aggregates the estimates sequentially, which reduces the estimate variance while accelerating convergence [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF].

Contributions. A fundamental aspect of this paper is to explore how changing data streams affect these stochastic algorithms. Our analysis extends the work of Moulines and Bach [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF] to a streaming framework. We investigate several kinds of data streams, from vanilla SG descent and ASG to more exotic learning designs such as time-varying mini-batch SG and ASG. Our main theoretical contribution is the non-asymptotic analysis of SG-based methods in this streaming framework. Our results show a noticeable improvement in convergence rates by having learning rates that adapt to the expected data streams. In particular, we show how to obtain optimal convergence rates robust to any data streaming rate.

Organization. Section 2 presents the streaming framework on which the non-asymptotic analysis relies. Our convergence results are presented in Sections 3 and 4, with and without averaging. Both sections includes analysis of unbounded and uniform bounded gradients. These theoretical results are illustrated in Section 5 for a variety of data streams. At last, some final remarks are done in Section 6.

Problem Formulation

The analysis of statistical and machine learning models often involves some form of optimization [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Goodfellow | Deep learning[END_REF][START_REF] Lecun | Deep learning[END_REF]. Many of these optimization problems aim to minimize functions of the form,

L(θ) = E[l t (θ)], (1)
with respect to θ ∈ Θ, where Θ ⊆ R d is a convex body; a convex body in R d is a compact convex set with nonempty interior. For streaming (or large-scale) problems, it would be too expensive to compute the full gradient ∇ θ L(θ). Instead, the minimization of L in ( 1) is achieved without evaluating it directly but by unbiased estimates of its gradients, namely, through ∇ θ l t : Θ → R. Observe that the principles for biased functions (l t ) are rather different [START_REF] Aspremont | Smooth optimization with approximate gradient[END_REF][START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF]. Thus, we let (l t ) constitute a sequence of independent differentiable random (possibly non-convex) functions and their gradients unbiased estimates of ∇ θ L, e.g., see Nesterov et al. [24] for definitions and properties of such functions.

Many problems, from classification, and regression to ranking, can be written on this form (1), e.g., see Teo et al. [START_REF] Teo | A scalable modular convex solver for regularized risk minimization[END_REF] for examples of scalar and vectorial loss functions and their derivatives. For example, consider the simple case where we have some samples {(X t , Y t )}, t = 1, . . . , n. Our interest is to find predictor h θ over some parameterization {h θ } θ∈Θ , by minimizing (1) with l t (θ) = l(h θ (X t ), Y t )+λΩ(θ), where l is some loss function, λ > 0 a regulizer parameter, and Ω : Θ → R some regularizer, e.g., the l 1 or l 2 regularization; here the loss l could be the quadratic, logistic, (squared) hinge, or Huber's (robust) loss, but it depends on the experiments that one wants to perform [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Lan | First-order and stochastic optimization methods for machine learning[END_REF]24].

Streaming framework. Let us now describe our streaming framework in which we will solve our problem in (1): at each time t ∈ N, a block consisting of n t ∈ N random functions l t = (l t,1 , . . . , l t,n t ) arrive. To solve this, we introduce the Stochastic Streaming Gradient (SSG), defined as

θ t = θ t-1 - γ t n t n t i=1 ∇ θ l t,i (θ t-1 ), θ 0 ∈ Θ, (2)
where (γ t ) is a decreasing sequence of positive numbers also referred to as the learning rate, satisfying t i=1 γ i = ∞ and t i=1 γ 2 i < ∞ for t → ∞ [START_REF] Robbins | A stochastic approximation method[END_REF]. In the same way, we introduce the Projected SSG (PSSG), defined by

θ t = P Θ        θ t-1 - γ t n t n t i=1 ∇ θ l t,i (θ t-1 )        , θ 0 ∈ Θ, (3)
where P Θ denotes the Euclidean projection onto the convex body Θ ⊆ R d , i.e., P Θ (θ) = arg min θ ∈Θ θθ 2 . The PSSG estimate in (3) is very convenient for models with conditions on the parameters space, and thereby, requires a projection of the parameters. Next, to guarantee optimal convergence properties [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Ruppert | Efficient estimations from a slowly convergent Robbins-Monro process[END_REF], we introduce the Polyak-Ruppert average of (2), called Averaged SSG (ASSG), given as

θt = 1 N t t-1 i=0 n i+1 θ i , θ0 = 0, (4)
where N t = t i=1 n i denotes the accumulated sum of observations. Similarly, we define the (Polyak-Ruppert) Average PSSG (APSSG) estimate as when ( θt ) (in (4)) is derived using [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF]. These assumptions are modified versions of the standard assumptions for stochastic approximations [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF][START_REF] Kushner | approximation and recursive algorithms and applications[END_REF] as they hold for any i = 1, . . . , n t . Note that Assumption 2-p only needs to hold for θ = θ * . By the smoothness assumption (Assumption 2-p), we avoid the unfavorable uniformly bounded gradients assumption, which is too restrictive and only holds for a few losses. Assumption 3-p is a weak assumption that should be seen as an assumption on Θ rather than on (l t,i ). For SSG and PSSG, we only need Assumptions 2-p and 3-p to hold for p = 2, whereas, for ASSG and APSSG, we need p = 4 in order to bound the fourth-order moment.

Our streaming framework include classic examples: stochastic approximation (Robbins-Monro setting [START_REF] Robbins | A stochastic approximation method[END_REF]) and learning from i.i.d. data, such as linear regression, logistic regression, general ridge regressions and quantile regression, p-means, and softmax regression, under regularity conditions [START_REF] Cardot | Efficient and fast estimation of the geometric median in hilbert spaces with an averaged stochastic gradient algorithm[END_REF][START_REF] Steinwart | Estimating conditional quantiles with the help of the pinball loss[END_REF][START_REF] Teo | A scalable modular convex solver for regularized risk minimization[END_REF].

In the following theorem, we derive an explicit upper bound on the t-th estimate of ( 2) and ( 3) for any learning rate (γ t ) using classical techniques from stochastic approximations [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF][START_REF] Kushner | approximation and recursive algorithms and applications[END_REF].

Theorem 1 (SSG/PSSG). Denote δ t = E[ θ t -θ * 2 ]
for some δ 0 ≥ 0, where (θ t ) follows (2) or (3). Under Assumption 1, Assumptions 2-p and 3-p with p = 2, we have for any learning rate (γ t ) that

δ t ≤ exp         -µ t i=t/2 γ i         π δ t + 2σ 2 µ max t/2≤i≤t γ i n i , (8)
with π δ t = exp(4C 2 l t i=1 γ 2 i /n i ) exp(2C 2 ∇ t i=1 1 {ni>1} γ 2 i )(δ 0 + 2σ 2 /C 2 l ).
Sketch of proof. Under Assumption 1, Assumptions 2-p and 3-p with p = 2, we derive from (2) that (δ t ) satisfies the recursive relation

δ t ≤ [1 -2µγ t + (2C 2 l + (n t -1)C 2 ∇ )n -1 t γ 2 t ]δ t-1 + 2σ 2 n -1 t γ 2 t , (9)
for any (n t ) and (γ t ) fulfilling the conditions imposed on the learning rate [START_REF] Robbins | A stochastic approximation method[END_REF]. This recursive relation is then bounded in a non-asymptotic manner using Proposition 5. Bounding the projected estimate in (3) follows directly from the fact that

E[ P Θ (θ) -θ * 2 ] ≤ E[ θ -θ * 2 ]
, ∀θ ∈ Θ [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]. Alternatively, the projected estimate can also be shown without Assumptions 2-p and 3-p but instead with a bounded gradient assumption (Assumption 5), e.g., see Moulines and Bach [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]. Related work. When n t = 1 in (8), we obtain the usual SG descent studied in Moulines and Bach [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]. Similarly, Theorem 1 provides an upper bound on the function values, E[L(θ t )-L(θ * )] ≤ C l δ t /2; this follows by Cauchy-Schwarz inequality and Assumption 2-p.

Natural decay imposed by Robbins and Monro [START_REF] Robbins | A stochastic approximation method[END_REF]. The learning rate (γ t ) should satisfy the following requirements:

t i=1 γ i = ∞ and t i=1 γ 2 i /n i ≤ t i=1 γ 2 i < ∞ for t → ∞.
These conditions directly imply that π δ t < ∞ as t → ∞. Thus, our attention is on reducing the noise term max t/2≤i≤t γ i /n i without damaging the natural decay of the sub-exponential term exp(-µ t i=t/2 γ i ). In particular, this non-asymptotic bound shows convergence in quadratic mean for any learning rate, fulfilling these conditions. In addition, the scaling with (n t ) in the noise term shows an apparent variance reduction when we increase the streaming batches (n t ).

Throughout this paper, we will consider learning rates on the form γ t = C γ n β t t -α with hyper-parameters C γ > 0, β ∈ [0, 1], and α chosen accordingly to the expected streaming batches denoted by n t . We start by considering constant streaming batches (i.e., mini-batch SG), where n t follows the constant streaming batch size C ρ ∈ N:

Corollary 1 (SSG/PSSG, constant streaming batches). Denote δ t = E[ θ t -θ * 2 ], where (θ t ) follows (2) or (3). Suppose γ t = C γ n β t t -α with n t = C ρ for C ρ ∈ N, such that α ∈ (1/2, 1
). Under Assumption 1, Assumptions 2-p and 3-p with p = 2, we have

δ t ≤ exp        - µC γ N 1-α t 2 1-α C 1-α-β ρ        π c ∞ + 2 1+α σ 2 C γ µC 1-α-β ρ N α t , (10)
where

π c ∞ = exp(4αC 2 γ (2C 2 l + C ρ 1 {Cρ>1} C 2 ∇ )/(2α -1)C 1-2β ρ )(δ 0 + 2σ 2 /C 2 l ) is a finite constant.
Decay of the initial conditions. The bound in Corollary 1 depends on the initial condition δ 0 = θ 0 -θ * 2 and the variance σ 2 in the noise term. The initial condition δ 0 vanish sub-exponentially fast for α ∈ (1/2, 1). Thus, the asymptotic term is 2 1+α σ 2 C γ /µC 1-α-β ρ N α t , i.e., δ t = O(N -α t ). Moreover, the bound in [START_REF] Godichon-Baggioni | An averaged projected robbins-monro algorithm for estimating the parameters of a truncated spherical distribution[END_REF] is optimal (up to some constants) for quadratic functions (l t,i ), since the deterministic recursion in (9) would be with equality. It is worth noting that if C γ C l or C γ C ∇ is chosen too large, they may produce a large π c ∞ constant. In addition, π c ∞ is positively affected by C ρ when β < 1/2. Obviously, the hyper-parameter β only comes into play if the streaming batch size is larger than one, i.e., C ρ > 1. Nonetheless, the effect of π c ∞ will decrease exponentially fast due to the sub-exponentially decaying factor in front.

Variance reduction. The asymptotic term is divided by C 1-α-β ρ , implying we could achieve variance reduction by taking α + β ≤ 1 when C ρ is large. Taking a large streaming batch size, e.g., C ρ = t, one accelerates the vanilla SG descent convergence rate to O(N 1-β t ). However, this large streaming batch size would be unsuitable in practice, and it would mean that we would take few steps until convergence is achieved.

The safe choice of having β = 0 functions well for the SSG method for any streaming batch size C ρ , but fixed-sized streaming batches are not the most realistic streaming setting. These streaming batches are far more likely to vary in size depending on the data streams. For the sake of simplicity, we consider time-varying streaming batches where n t are on the form C ρ t ρ with C ρ ∈ N and ρ ∈ (-1, 1) such that n t ≥ 1 for all t. We will refer to ρ as the streaming rate. For the convenience of notation, let ρ = ρ1 {ρ≥0} .

Corollary 2 (SSG/PSSG, time-varying streaming batches). Denote

δ t = E[ θ t -θ * 2 ], where (θ t ) follows (2) or (3). Suppose γ t = C γ n β t t -α with n t = C ρ t ρ for C ρ ∈ N and ρ ∈ (-1, 1), such that α -β ρ ∈ (1/2, 1). Under Assumption 1,
Assumptions 2-p and 3-p with p = 2, we have

δ t ≤ exp        - µC γ N 1-φ t 2 (2+ρ)(1-φ) C 1-β-φ ρ        π v ∞ + 2 1+(2+ρ)φ σ 2 C γ µC (1-β)1 {ρ≥0} -φ ρ N φ t , (11)
where φ = ((1 -β) ρ + α)/(1 + ρ) and π v ∞ = exp(4(α -β ρ)C 2 γ C 2β ρ (2C 2 l + C 2 ∇ )/(2(α -β ρ) -1))(δ 0 + 2σ 2 /C 2 l ) is a finite constant.
Decay of the initial conditions. As mentioned for Corollary 1, the condition of having αβ ρ ∈ (1/2, 1) is a natural restriction from Robbins and Monro [START_REF] Robbins | A stochastic approximation method[END_REF], which relaxes the usual condition of having α ∈ (1/2, 1) for ρ nonnegative. For ρ ∈ (-1, 1/2), setting α = 2/3 and β = 1/3 would give same decay rate, δ t = O(N -2/3 t ) as we saw for Corollary 1 when α = 2/3. However, accelerated convergence could be achieved by, e.g., setting α = 1 and β = 1/2 for streaming rate ρ ∈ (0, 1), giving us

δ t = O(N -(1+ρ/2)/(1+ρ) t ).
Variance reduction. Similarly to Corollary 1, the sub-exponential and asymptotic term is scaled by C 1-β-φ ρ for ρ ≥ 0, implying we should take α + β ≤ 1 to obtain variance reduction. These conclusions will change when we consider the averaging estimate in Section 4.

The reasoning in Corollary 2 could be expanded to include random streaming batches where n t is given such that

C L t ρ L ≤ n t ≤ C H t ρ H with ρ L , ρ H ∈ (-1, 1) and C L , C H ≥ 1.
This yields the modified rate φ = ((1 -β)ρ L + α)/(1 + ρ H ); nevertheless, we will leave the proof to the reader.

Averaged Stochastic Streaming Gradients

In what follows, we consider the averaging estimate ( θn ) given in (4) derived with use of (θ t ) from (2) (Section 4.1) or (3) (Section 4.2). Besides having Assumptions 2-p and 3-p to hold for p = 4, an additional assumption is needed for bounding the rest term of the averaging estimate. Assumption 4. There exists a non-negative self-adjoint operator Σ such that E[∇ θ l t,i (θ * )∇ θ l t,i (θ * ) ] Σ.

Note that the operator Σ always exists when σ is finite for order p = 4 in Assumption 3-p.

Unbounded Gradients

As in Section 3, we conduct a general study for any learning rate (γ t ) when applying the Polyak-Ruppert averaging estimate from (4):

Theorem 2 (ASSG). Denote δt = E[ θt -θ * 2 ]
with ( θt ) given by (4), where (θ t ) follows (2). Under Assumption 1, Assumptions 2-p and 3-p with p = 4, and Assumption 4, we have for any learning rate (γ t ) that

δ1/2 t ≤ Λ 1/2 N 1/2 t + 1 µN t t-1 i=1 n i+1 γ i+1 - n i γ i δ 1/2 i + n t µγ t N t δ 1/2 t + n 1 µN t 1 γ 1 + C l δ 1/2 0 + C l µN t         t-1 i=1 n i+1 δ i         1/2 + C ∇ µN t t-1 i=0 n i+1 ∆ 1/2 i , (12)
where

Λ = Tr(∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ) and ∆ t = E[ θ t -θ * 4 ]
for some ∆ 0 ≥ 0. As noticed in Polyak and Juditsky [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], the leading term Λ/N t achieves the Cramer-Rao lower bound [START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF][START_REF] Murata | Statistical analysis of learning dynamics[END_REF]. Note that the leading term Λ/N t is invariant of the learning rate (γ t ). Moreover, this bound of O(N -1 t ) is achieved without inverting the Hessian. Next, the processes (δ t ) and (∆ t ) can be bounded by the recursive relations in [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF] and [START_REF] Murata | Statistical analysis of learning dynamics[END_REF]. There are no sub-exponential decaying terms for the initial conditions in Theorem 2, which is a common problem for averaging. However, as mentioned previously, we are more interested in advancing the decay of the asymptotic terms. To ease notation, we make use of the functions ψ y x (t) : R → R, given as

ψ y x (t) =              t (1-x)/(1+y) /(1 -x) if x < 1, (1 + y) log(t) if x = 1, x/(x -1) if x > 1, with y ∈ R + , such that t i=1 i -x ≤ ψ 0 x (t)
δ1/2 t ≤ Λ 1/2 N 1/2 t + 6σC (1-α-β)/2 ρ µ 3/2 C 1/2 γ N 1-α/2 t + 2 α 6C ∇ σ 2 C γ µ 2 C 1-α-β ρ N α t + 2C l σC 1/2 γ µ 3/2 C (1-α-β)/2 ρ N (1+α)/2 t + C ρ Γ c µN t + C 2-α-β ρ √ π c ∞ A c ∞ µC γ N 2-α t + (6 + 71 {Cρ>1} )2 3α/2 C ∇ σ 2 C 3/2 γ C 3β/2 ρ ψ 0 3α/2 (N t /C ρ ) µ 3/2 N t , with Γ c given by (1/C γ C β ρ + C l )δ 1/2 0 + C l π c ∞ A c ∞ /C ρ + √ π c ∞ A c ∞ /C γ C β ρ + C ∇ √ Π c ∞ A c ∞ , consisting of the finite constants π c ∞ , Π c ∞ and A c ∞ , that only depends on µ, δ 0 , ∆ 0 , C l , σ, C ∇ , C ∇ , C γ , C ρ , β and α.
Accelerated decay the initial conditions. By averaging, we have increased the rate of convergence from O(N -α t ) to the optimal rate O(N -1 t ). The two subsequent terms are the main remaining terms decaying at the rate O(N α-2

t

) and O(N -2α t ), which suggests setting α = 2/3 would be optimal. The remaining terms are negligible. Next, it is worth noting that having α + β = 1 in Corollary 3, we would give no impact in the main remaining terms from the streaming batch size C ρ . Moreover, taking α = 2/3 and β ≤ 1/3 would be an optimal choice of hyper-parameters such that the streaming batch size C ρ have a positive or no impact. At last, as we do not rely on sub-exponentially decaying terms, we need to be more careful when picking our hyper-parameters, e.g., taking C γ C l too large may cause Γ c to be significant. Nevertheless, the term consisting of Γ c decay at a rate of at least O(N -2 t ).

Corollary 4 (ASSG, time-varying streaming batches). Denote δt = E[ θtθ * 2 ] with ( θt ) given by (4), where (θ t ) follows (2). Suppose γ t = C γ n β t t -α with n t = C ρ t ρ for C ρ ∈ N and ρ ∈ (-1, 1), such that αβ ρ ∈ (1/2, 1). Under Assumption 1, Assumptions 2-p and 3-p with p = 4, and Assumption 4, we have

δ1/2 t ≤ Λ 1/2 N 1/2 t + 2 3+φ(1+ ρ) σC (1-φ-β)/21 {ρ≥0} ρ µ 3/2 C 1/2 γ N 1-φ/2 t + 2 (1+φ)(1+ ρ)-2 C ∇ σ 2 C γ µ 2 C 1-φ-β ρ N φ t + 2 φ(1+ ρ)/2 C l σC 1/2 γ µ 3/2 C (1-φ-β)/21 {ρ≥0} ρ N (1+φ)/2 t + C ρ Γ v µN t + C 2-φ-β ρ √ π v ∞ A v ∞ µC γ N 2-φ t + 2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ ψ ρ 3(α-β ρ)/2 (N t /C ρ ) µ 3/2 C 1 {ρ≥0} ρ N t , with Γ v given by (1/C γ C β ρ + C l )δ 1/2 0 + 2 ρC l π v ∞ A v ∞ /C ρ + 2 √ π v ∞ A v ∞ /C γ C β ρ + 2 ρC ∇ √ Π v ∞ A v ∞ , consisting of the finite constants π v ∞ , Π v ∞ and A v ∞ , that only depends on µ, δ 0 , ∆ 0 , C l , σ, C ∇ , C ∇ , C γ , C ρ , β and α.
Robustness towards streaming rates ρ: Following the arguments above, the two main remainder terms reveal that φ = 2/3 ⇔ αβ ρ = (2ρ)/3, e.g., by setting β = 0, we should pick α = (2ρ)/3. Likewise, if ρ = 0, we yield the same conclusion as in Corollary 3, namely α = 2/3. However, these hyper-parameter choices are not resilient against any arrival schedule ρ. Nonetheless, we can robustly achieve φ = 2/3 for any ρ ∈ (-1, 1) by setting α = 2/3 and β = 1/3. In other words, we can achieve optimal convergence for any data stream by having α = 2/3 and β = 1/3.

Bounded Gradients

In what follows, we consider the averaging estimate θn given in (4) but with the use of the projected estimate PSSG from (3). To avoid calculating the six-order moment, we make the unnecessary assumption that ∇ θ l t,i (θ) is uniformly bounded for any θ ∈ Θ; the derivation of the six-order moment can be found in Godichon-Baggioni [START_REF] Godichon-Baggioni | Estimating the geometric median in hilbert spaces with stochastic gradient algorithms: Lp and almost sure rates of convergence[END_REF].

Assumption 5. Let D θ = inf θ∈∂Θ θ -θ * > 0 with ∂Θ denoting the frontier of Θ. Moreover, there exists G Θ > 0 such that ∀t ≥ 1, sup θ∈Θ ∇ θ l t,i (θ) 2 ≤ G 2 Θ a.s., with i = 1, . . . , n t .
Corollary 5 (APSSG, constant streaming batches). Denote δt = E[ θtθ * 2 ] with ( θt ) given by (4), where (θ t ) follows (3). Suppose γ t = C γ n β t t -α with n t = C ρ for C ρ ∈ N, such that α ∈ (1/2, 1). Under Assumption 1, Assumptions 2-p and 3-p with p = 4, Assumptions 4 and 5, we have

δ1/2 t ≤ Λ 1/2 N 1/2 t + 6σC (1-α-β)/2 ρ µ 3/2 C 1/2 γ N 1-α/2 t + 2 α 6C ∇ σ 2 C γ µ 2 C 1-α-β ρ N α t + 2C l σC 1/2 γ µ 3/2 C (1-α-β)/2 ρ N (1+α)/2 t + C ρ Γ c µN t + C 2-α-β ρ √ π c ∞ A c ∞ µC γ N 2-α t + (6 + 71 {Cρ>1} )2 3α/2 C ∇ σ 2 C 3/2 γ C 3β/2 ρ ψ 0 3α/2 (N t /C ρ ) µ 3/2 N t , with C ∇ = C ∇ +2 2 G Θ /D 2 θ and Γ c given by (1/C γ C β ρ +C l )δ 1/2 0 +C l π c ∞ A c ∞ /C ρ + √ π c ∞ A c ∞ /C γ C β ρ +C ∇ √ Π c ∞ A c ∞ , consisting of the finite constants π c ∞ , Π c ∞ and A c ∞ , that only depends on µ, δ 0 , ∆ 0 , C l , σ, C ∇ , C ∇ , C γ , C ρ , β and α.
Corollary 6 (APSSG, time-varying streaming batches). Denote δt = E[ θtθ * 2 ] with ( θt ) given by (4), where

(θ t ) follows (3). Suppose γ t = C γ n β t t -α with n t = C ρ t ρ for C ρ ∈ N and ρ ∈ (-1, 1), such that α -β ρ ∈ (1/2, 1
). Under Assumption 1, Assumptions 2-p and 3-p with p = 4, Assumptions 4 and 5, we have

δ1/2 t ≤ Λ 1/2 N 1/2 t + 2 3+φ(1+ ρ) σC (1-φ-β)/21 {ρ≥0} ρ µ 3/2 C 1/2 γ N 1-φ/2 t + 2 (1+φ)(1+ ρ)-2 C ∇ σ 2 C γ µ 2 C 1-φ-β ρ N φ t + 2 φ(1+ ρ)/2 C l σC 1/2 γ µ 3/2 C (1-φ-β)/21 {ρ≥0} ρ N (1+φ)/2 t + C ρ Γ v µN t + C 2-φ-β ρ √ π v ∞ A v ∞ µC γ N 2-φ t + 2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ ψ ρ 3(α-β ρ)/2 (N t /C ρ ) µ 3/2 C 1 {ρ≥0} ρ N t , with C ∇ = C ∇ + 2 2 G Θ /D 2 θ and Γ v given by (1/C γ C β ρ + C l )δ 1/2 0 + 2 ρC l π v ∞ A v ∞ /C ρ + 2 √ π v ∞ A v ∞ /C γ C β ρ + 2 ρC ∇ √ Π v ∞ A v ∞ , consisting of the finite constants π v ∞ , Π v ∞ and A v ∞ , that only depends on µ, δ 0 , ∆ 0 , C l , σ, C ∇ , C ∇ , C γ , C ρ , β and α.

Experiments

In this section, we demonstrate the theoretical results presented in Sections 3 and 4 for various data streams. In Section 5.1, we illustrate the unbounded gradient case (Sections 3 and 4.1) using linear regression. Where in Section 5.2, we present the bounded gradient case (Sections 3 and 4.2) by considering the geometric median. To measure the performance, we use the quadratic mean error of the parameter estimates over one-hundred replications, given by (E[ θ N t -θ * 2 ]) t≥1 . Note that averaging over several iterations gives a reduction in variability, which mainly benefits the SSG and PSSG.

Linear Regression

Consider the linear regression defined by y t = X T t θ + t , where X t ∈ R d is a random features vector, θ ∈ R d is the parameters vector, and t is a random variable with zero mean, independent from X t . Moreover, (X t , t ) t≥1 are independent and identically distributed. Thus, θ * is the minimizer of

L(θ) = E[(y t -X T t θ) 2 ].
In this example, we fix d = 10, set θ = (-4, -3, 2, 1, 0, 1, 2, 3, 4, 5) T ∈ R 10 , and let (X t ) and ( t ) be standard Gaussian. It is well-known that C γ can substantially impact convergence; when C γ is too large, instability can occur, leading to an explosion during the first iterations. If C γ is too small, the convergence can become very slow and destroy the desired rate α. To focus on the various data streams, we set C γ = 1/2 and α = 2/3.

In Figure 1a, we consider constant data streams to illustrate the results in Corollaries 1 and 3. The figures show a solid decay rate proportional to α = 2/3 for any streaming batch size C ρ ∈ {1, 8, 64, 128} with β = 0, as shown in Corollary 1. In addition, we see an acceleration in decay by averaging, as explained in Corollary 3. Both methods show a noticeable reduction in variance when C ρ increases which are particularly beneficial in the beginning. Moreover, as mentioned in Remark 1, the stationary phase may also commence earlier when we raise the streaming batch size C ρ . Next, in Figures 1b to 1e, we vary the streaming rate ρ for streaming batch sizes C ρ = 1, 8, 64, and 128, respectively, with β = 0. These figures shows an increase in decay of the SSG when the streaming rate ρ increase. Despite this, we still achieve better convergence for the ASSG method, which seems more immune to the different choices of streaming rate ρ, e.g., see the discussion after Corollary 4. We know this from Corollary 2, as φ = ( ρ + α)/(1 + ρ) ≥ α for β = 0. In addition, we see that C ρ has a positive effect on the noise (i.e, variance reduction), but if C ρ becomes too large, it may slow down convergence (as seen in Figure 1e). Alternatively, we could think around the problem in another way: how can we choose α and β such that we have optimal decay of φ = 2/3 for any ρ. In other words, for any arrival schedule that may occur, how should we choose our hyper-parameters such that we achieve optimal decay of φ = 2/3. As discussed after Corollary 4, one example of this could be achieved by setting α = 2/3 and β = 1/3 such that φ = 2/3 for any ρ. Figure 1f shows an example of this where we (indeed) achieve the same decay rate for any streaming rate ρ.

Geometric Median

The geometric median is a generalization of the real median introduced by Haldane [START_REF] Haldane | Note on the median of a multivariate distribution[END_REF]. Robust estimators such as the geometric median may be preferred over the mean when the data is noisy. Moreover, in our streaming framework, stochastic algorithms are preferred as they efficiently handle large samples of high-dimensional data [START_REF] Cardot | Efficient and fast estimation of the geometric median in hilbert spaces with an averaged stochastic gradient algorithm[END_REF][START_REF] Godichon-Baggioni | Estimating the geometric median in hilbert spaces with stochastic gradient algorithms: Lp and almost sure rates of convergence[END_REF]. The geometric median of

X ∈ R d is defined by θ * ∈ R d which minimizes the convex function L(θ) = E[ X -θ -X ],
e.g., see Gervini [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF], Kemperman [START_REF] Kemperman | The median of a finite measure on a banach space[END_REF] for properties such as existence, uniqueness, and robustness (breakdown point). Thus, the gradient

∇ θ L(θ) = E[∇ θ l t (θ)] with ∇ θ l t (θ) = -(X t -θ)/ X t -θ is bounded as ∇ θ l t (θ) ≤ 1.
As discussed in Gadat and Panloup [START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF], this boundedness implies that it is not necessary to project our estimates explicitly. Similarly to above, we fix d = 10 and let (X t ) be standard Gaussian centered at θ = (-4, -3, 2, 1, 0, 1, 2, 3, 4, 5) T ∈ R 10 . Moreover, following the reasoning of Cardot et al. [START_REF] Cardot | Efficient and fast estimation of the geometric median in hilbert spaces with an averaged stochastic gradient algorithm[END_REF], we set

C γ = √ d = √ 10
, and let α = 2/3. Figure 2a shows the variance reduction effect for different constant streaming batches C ρ with β = 0. However, the robustness of the geometric median leaves only a small positive impact for further variance reduction. Thus, too large (constant) streaming batch sizes C ρ hinders the convergence as we make too few iterations. These findings can be extended to Figures 2b to 2e, where we vary the streaming rate ρ for streaming batch sizes C ρ = 1, 8, 64, and 128, respectively, with β = 0. The lack of convergence improvements comes from β = 0, which means we do not exploit the potential of using more observations to accelerate convergence. As shown in Figure 2f, we can achieve this acceleration by simply taking β = 1/3. In addition, β = 1/3 provides optimal convergence robust to any streaming rate ρ. Choosing a proper β > 0 is particularly important when C ρ is large, as robustness is an integral part of the geometric median method.

Conclusions

We considered the stochastic approximation problem in a streaming framework where we had to minimize convex objectives using only unbiased estimates of its gradients. We introduced and studied the convergence rates of the stochastic streaming algorithms in a non-asymptotic manner. This investigation was derived using learning rates of the form γ t = C γ n β t t -α under time-varying data streams (n t ). The theoretical results and our experiments showed a noticeable improvement in the convergence rate by choosing the learning rate (hyper-parameters) according to the expected data streams. For ASSG and APSSG, we showed that this choice of learning rate led to optimal convergence rates and was robust to any data stream rate we may encounter. Moreover, in large-scale learning problems, we know how to accelerate convergence and reduce variance through the learning rate and the treatment pattern of the data.

There are several ways to expand our work but let us give some examples: first, we can extend our analysis to include streaming batches of any size in the spirit of the discussion after Corollary 2. Second, many machine learning problems encounter correlated variables and high-dimensional data, making an extension to non-strongly convex objectives advantageous Bach and Moulines [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n)[END_REF], e.g, in Werge and Wintenberger [START_REF] Werge | Adavol: An adaptive recursive volatility prediction method[END_REF], they use SG-based optimization methods for volatility prediction through GARCH modeling. Third, Assumption 1 requires unbiased (and independent) gradient estimates, thus, an obvious extension could incorporate a more realistic dependency assumption, thereby increasing the applicability for more models. Moreover, studying dependence may give insight into how to process dependent information optimally. Next, a natural extension would be to modify our averaging estimate from (4) to a weighted averaged version (WASSG) proposed by Mokkadem and Pelletier [START_REF] Mokkadem | A generalization of the averaging procedure: The use of two-time-scale algorithms[END_REF] and Boyer and Godichon-Baggioni [START_REF] Boyer | On the asymptotic rate of convergence of stochastic newton algorithms and their weighted averaged versions[END_REF], given as for λ > 0 with (θ t ) following ( 2) or (3). We can limit the effect of bad initializations by placing more weight on the newest estimates. Following the demonstrations in Section 5, an example of this WASSG estimate ( θt,λ ) can be found in Figure 3 with use of λ = 2. Here we see that although the WASSG estimate in (13) may not achieve a better final error (compared to the ASSG and APSSG estimates in Figures 1f and2f), but it still achieves a better decay along the way, often referred to as parameter tracking.

θt,λ = 1 t i=1 n i log(1 + i) λ t i=1 n i log(1 + i) λ θ i-1 , θ0,λ = 0, (13)

Proofs

In this section, we provide detailed proofs of the results presented in the manuscript. Purely technical results used in the proofs can be found in Appendix A. Let (F t ) t≥1 be an increasing family of σ-fields, namely F t = σ(l 1 , . . . , l t ) with l t = (l t,1 , . . . , l t,n t ). Furthermore, we expand the notation with F t-1,i = σ(l 1,1 , . . . , l t-1,n t-1 , l t,1 , . . . , l t,i ) such that F t-1,0 = F t-1 . Meaning, ∀0 ≤ i < j, we have F t-1 ⊆ F t-1,i ⊂ F t-1, j . Thus, by the independence of the random (differentiable) functions (l t,i ), Assumption 1 yields that ∀t ≥ 1, E[∇ θ l t,i (θ t-1 )|F t-1,i-1 ] = ∇ θ L(θ t-1 ) with i = 1, . . . , n t .

Proofs for Section 3

The section is structured such that we start by analyzing the recursive relations and bounding them for every choice of learning rate. Next, we look at specific choices of learning rates.

Proof of Theorem 1. Taking the quadratic norm on both sides of (2), expanding it, and take the conditional expectation, yields

E[ θ t -θ * 2 |F t-1 ] = θ t-1 -θ * 2 + γ 2 t n 2 t E          n t i=1 ∇ θ l t,i (θ t-1 ) 2 F t-1          - 2γ t n t n t i=1 E[ ∇ θ l t,i (θ t-1 ), θ t-1 -θ * |F t-1 ]. ( 14
)
To bound the second term (on the right-hand side) of ( 14), we first expand it as follows,

n t i=1 E[ ∇ θ l t,i (θ t-1 ) 2 |F t-1 ] + n t i j E[ ∇ θ l t,i (θ t-1 ), ∇ θ l t, j (θ t-1 ) |F t-1 ]. ( 15
)
For first term of (15), we utilize the Lipschitz continuity of ∇ θ l t,i , together with Assumptions 1 to 3-p, to obtain

E[ ∇ θ l t,i (θ t-1 ) 2 |F t-1 ] ≤ 2E[ ∇ θ l t,i (θ t-1 ) -∇ θ l t,i (θ * ) 2 |F t-1 ] + 2E[ ∇ θ l t,i (θ * ) 2 |F t-1 ] ≤ 2C 2 l θ t-1 -θ * 2 + 2σ 2 , (16)
using x + y 2 ≤ 2( x 2 + y 2 ). Next, for the second term in [START_REF] Kurdyka | On gradients of functions definable in o-minimal structures[END_REF]: as F t-1 ⊆ F t-1,i ⊂ F t-1, j for all 0 ≤ i < j, we have

E[ ∇ θ l t,i (θ t-1 ), ∇ θ l t, j (θ t-1 ) |F t-1 ] = E[E[ ∇ θ l t,i (θ t-1 ), ∇ θ L(θ t-1 ) |F t-1,i-1 ]|F t-1 ],
since θ t-1 and l t,i are F t-1, j-1 -measurable for all 0 ≤ i < j, and similarly, as θ t-1 is F t-1 -measurable and F t-1,i-1measurable for all i ≥ 0, we also have

E[E[ ∇ θ l t,i (θ t-1 ), ∇ θ L(θ t-1 ) |F t-1,i-1 ]|F t-1 ] = E[ E[∇ θ l t,i (θ t-1 )|F t-1,i-1 ], ∇ θ L(θ t-1 ) |F t-1 ] = ∇ θ L(θ t-1 ) 2 ,
where

∇ θ L(θ t-1 ) 2 ≤ C 2 ∇ θ t-1 -θ * 2 as ∇ θ L is C ∇ -Lipschitz continuous and ∇ θ L(θ * ) = 0.
Thus, we obtained a bound for the second term (on the right-hand side) of ( 14) using the bounds of the two terms in [START_REF] Kurdyka | On gradients of functions definable in o-minimal structures[END_REF]:

n t i=1 (2C 2 l θ t-1 -θ * 2 + 2σ 2 ) + n t i j C 2 ∇ θ t-1 -θ * 2 = (2C 2 l n t + C 2 ∇ (n t -1)n t ) θ t-1 -θ * 2 + 2σ 2 n t . (17)
For the third term (on the right-hand side) of ( 14) we use that L is µ-quasi-strong convex and θ t-

1 is F t-1 -measurable, E[ ∇ θ l t,i (θ t-1 ), θ t-1 -θ * |F t-1 ] = E[∇ θ l t,i (θ t-1 )|F t-1 ], θ t-1 -θ * = ∇ θ L(θ t-1 ), θ t-1 -θ * ≥ µ θ t-1 -θ * 2 , (18)
by Assumption 1. Combining inequalities from ( 17) and ( 18) into ( 14) and taking the expectation on both sides of the inequality, yields the recursive relation ( 9):

δ t ≤ [1 -2µγ t + (2C 2 l + (n t -1)C 2 ∇ )n -1 t γ 2 t ]δ t-1 + 2σ 2 n -1 t γ 2 t , with δ t = E[ θ t -θ * 2 ]
with some δ 0 ≥ 0. At last, by Proposition 5, we obtain the desired inequality in (8), namely

δ t ≤ exp         -µ t i=t/2 γ i         exp        4C 2 l t i=1 γ 2 i n i        exp        2C 2 ∇ t i=1 1 {ni>1} γ 2 i              δ 0 + 2σ 2 C 2 l       + 2σ 2 µ max t/2≤i≤t γ i n i .
using that (n t -1)n -1

t ≤ 1 {nt >1} , n t ≥ 1, and that max 1≤i≤t 2σ 2 /(2C

2 l + (n i -1)C 2 ∇ ) ≤ max 1≤i≤t 2σ 2 /2C 2 l = σ 2 /C 2 l . Remark 1.
The decrease of (2C 2 l + (n t -1)C 2 ∇ )n -1 t γ t determines when the stationary phase occurs. This is more clearly seen in Proposition 4, where the inner terms directly depend on the inception of the stationary phase. Thus, by increasing n t , we decrease (2C 2 l + (n t -1)C 2 ∇ )n -1 t γ t , and especially it dominates the constant C l . Proof of Corollary 1. By Theorem 1, we have the upper bound giving as

δ t ≤ exp         -µ t i=t/2 γ i         π c t + 2σ 2 µC ρ max t/2≤i≤t γ i . (19)
as

n t = C ρ , with π c t = exp((4C 2 l /C ρ ) t i=1 γ 2 i ) exp(2C 2 ∇ 1 {Cρ>1} t i=1 γ 2 i )(δ 0 +σ 2 /C 2 l ). The sum term t i=1 γ 2 i = C 2 γ C 2β ρ t i=1 i -2α in π c
t can be bounded with the help of integral tests for convergence,

t i=1 i -2α = 1 + t i=2 i -2α ≤ 1 + t 1 x -2α dx ≤ 1 + 1/(2α -1) = 2α/(2α -1), as α ∈ (1/2, 1). Likewise, plugging γ t = C γ C β ρ t -α into the first term of (19), gives exp         -µ t i=t/2 γ i         = exp         -µC γ C β ρ t i=t/2 i -α         ≤ exp -µC γ C β ρ t t/2 x -α dx ≤ exp        - µC γ C β ρ t 1-α 2 1-α        ,
using the integral test for convergence. Next, as (γ t ) t≥1 is decreasing, then max t/2≤i≤t γ t = γ t/2 . Combining all these findings into [START_REF] Lojasiewicz | A topological property of real analytic subsets[END_REF], gives us

δ t ≤ exp        - µC γ C β ρ t 1-α 2 1-α        π c ∞ + 2 1+α σ 2 C γ µC 1-β ρ t α , (20)
with

π c ∞ = exp(4αC 2 γ (2C 2 l + C ρ 1 {Cρ>1} C 2 ∇ )/(2α -1)C 1-2β ρ )(δ 0 + 2σ 2 /C 2 l ).
At last, converting [START_REF] Mokkadem | A generalization of the averaging procedure: The use of two-time-scale algorithms[END_REF] into terms of N t using N t = C ρ t, yields the desired.

Proof of Corollary 2. For convenience, divided the proof into two cases to comprehend that n t ≥ 1 for all t: first, we bound each term of (8) (from Theorem 1) after inserting,

γ t = C γ n β t t -α = C γ C β ρ t βρ-α if ρ ≥ 0, or γ t ≥ C γ t -α
if ρ < 0 (using that β ≥ 0) into the inequality. If ρ ≥ 0, the first term of (8) can be bounded, as follows:

exp         -µ t i=t/2 γ i         = exp         -µC γ C β ρ t i=t/2 i βρ-α         ≤ exp        - µC γ C β ρ t 1+βρ-α 2 1+βρ-α        ,
using that αβρ ∈ (1/2, 1) and the integral test for convergence. In a same way, if ρ < 0, one has exp

        -µ t i=t/2 γ i         ≤ exp         -µC γ t i=t/2 i -α         ≤ exp - µC γ t 1-α 2 1-α .
Likewise, with the help of integral tests for convergence, we have for ρ ≥ 0, that t i=1

γ 2 i /n i ≤ t i=1 γ 2 i ≤ 2(α - βρ)C 2 γ C 2β ρ /(2(α -βρ) -1), as n t ≥ 1 and α -ρβ > 1/2. If ρ < 0, one has t i=1 γ 2 i /n i ≤ t i=1 γ 2 i ≤ 2αC 2 γ C 2β ρ /(2α -1) since C ρ ≥ n t ≥ 1.
Next, as (1 -β)ρ + α > 0 for ρ ≥ 0, then we can bound the last term of (8) by

2σ 2 µ max t/2≤i≤t γ i n i = 2σ 2 C γ µC 1-β ρ max t/2≤i≤t 1 i (1-β)ρ+α ≤ 2 1+(1-β)ρ+α σ 2 C γ µC 1-β ρ t (1-β)ρ+α
.

Likewise, if ρ < 0, we have

2σ 2 µ max t/2≤i≤t γ i n i = 2σ 2 C γ µ max t/2≤i≤t 1 n 1-β i i α ≤ 2 1+α σ 2 C γ µt α ,
since n t ≥ 1 and β ≤ 1. Combining all these findings gives

δ t ≤ exp         - µC γ C β1 {ρ≥0} ρ t (1-φ)(1+ ρ) 2 (1-φ)(1+ ρ)         π v ∞ + 2 1+φ(1+ ρ) σ 2 C γ µC (1-β)1 {ρ≥0} ρ t φ(1+ ρ) , (21)
where

π v ∞ = exp(4(α -β ρ)C 2 γ C 2β ρ (2C 2 l + C 2 ∇ )/2(α -β ρ) -1
) with ρ = ρ1 {ρ≥0} and φ = ((1 -β) ρ + α)/(1 + ρ). To write this in terms of N t , we use that

N t = t i=1 n i = C ρ t i=1 i ρ = C ρ (t ρ + t-1 i=1 i ρ ) ≤ C ρ (t ρ + t 1 x ρ dx) ≤ C ρ (t ρ + t ρ t 1 dx) = C ρ (t ρ + t 1+ρ ) ≤ 2C ρ t 1+ρ
, for ρ ≥ 0, thus, t ≥ (N t /2C ρ ) 1/(1+ρ) . For ρ < 0, we have N t ≤ C ρ t, i.e, t ≥ N t /C ρ .

Proofs for Section 4

Lemma 1 (ASSG/APSSG). Denote ∆ t = E[ θ t -θ * 4 ] for some ∆ 0 ≥ 0, where (θ t ) follows (2) or (3). Under Assumption 1, Assumptions 2-p and 3-p with p = 4 and Assumption 4, we have for any learning rate (γ t ) that

∆ t ≤ exp         -µ t i=t/2 γ i         Π ∆ t + 32σ 4 µ 2 max t/2≤i≤t γ 2 i n 2 i + 48σ 4 µ max t/2≤i≤t γ 3 i n 3 i + 114σ 4 µ max t/2≤i≤t γ 3 i 1 {ni>1} n 2 i , (22)
with Π ∆ t given in [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF]. Proof of Lemma 1. We will now derive the recursive step sequence for the fourth-order moment using the same arguments as in proof for Theorem 1. Thus, one can show that

E[ θ t -θ * 4 |F t-1 ] ≤ θ t-1 -θ * 4 + γ 4 t n 4 t E          n t i=1 ∇ θ l t,i (θ t-1 ) 4 F t-1          + 4γ 2 t n 2 t E         n t i=1 ∇ θ l t,i (θ t-1 ) , θ t-1 -θ * 2 F t-1         + 2γ 2 t n 2 t θ t-1 -θ * 2 E          n t i=1 ∇ θ l t,i (θ t-1 ) 2 F t-1          - 4γ t n t θ t-1 -θ * 2 n t i=1 E[∇ θ l t,i (θ t-1 )|F t-1 ], θ t-1 -θ * + 4γ 3 t n 3 t E          n t i=1 ∇ θ l t,i (θ t-1 ) 2 n t i=1 ∇ θ l t,i (θ t-1 ) , θ t-1 -θ * F t-1          , θ t-1 is F t-1 -measurable. Note, by Assumption 1, we have E[∇ θ l t,i (θ t-1 )|F t-1 ], θ t-1 -θ * = ∇ θ L(θ t-1 ), θ t-1 -θ * ≥ µ θ t-1 -θ * 2 ,
as L is µ-quasi-strong convex. Combining this with the Cauchy-Schwarz inequality x, y ≤ x y , we obtain the simplified expression:

E[ θ t -θ * 4 |F t-1 ] ≤ θ t-1 -θ * 4 + γ 4 t n 4 t E          n t i=1 ∇ θ l t,i (θ t-1 ) 4 F t-1          + 6γ 2 t n 2 t θ t-1 -θ * 2 E          n t i=1 ∇ θ l t,i (θ t-1 ) 2 F t-1          -4µγ t θ t-1 -θ * 4 + 4γ 3 t n 3 t θ t-1 -θ * E          n t i=1 ∇ θ l t,i (θ t-1 ) 3 F t-1          .
Next, recall Young's Inequality, i.e., for any a t , b t , c t > 0 we have

a t b t ≤ a 2 t c 2 t /2 + b 2 t /2c 2 t , n t i=1 ∇ θ l t,i (θ t-1 ) 3 ≤ γ t 2n t θ t-1 -θ * n t i=1 ∇ θ l t,i (θ t-1 ) 4 + 2n t θ t-1 -θ * γ t n t i=1 ∇ θ l t,i (θ t-1 ) 2 ,
giving us

E[ θ t -θ * 4 |F t-1 ] ≤(1 -4µγ t ) θ t-1 -θ * 4 + 3γ 4 t n 4 t E          n t i=1 ∇ θ l t,i (θ t-1 ) 4 F t-1          + 8γ 2 t n 2 t θ t-1 -θ * 2 E          n t i=1 ∇ θ l t,i (θ t-1 ) 2 F t-1          . (23)
To bound the second and fourth-order terms in [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF], we would need to study the recursive sequences: firstly, utilizing the Lipschitz continuity of ∇ θ l t,i , together with Assumptions 2-p and 3-p, and that θ t-1 is F t-1 -measurable (Assumption 1), we obtain

E[ ∇ θ l t,i (θ t-1 ) p |F t-1 ] ≤ 2 p-1 [E[ ∇ θ l t,i (θ t-1 ) -∇ θ l t,i (θ * ) p |F t-1 ] + E[ ∇ θ l t,i (θ * ) p |F t-1 ]] ≤ 2 p-1 [C p l θ t-1 -θ * p + σ p ], (24)
for any p ∈ [1, 4] using the bound x + y p ≤ 2 p-1 ( x p + y p ). Thus, we can bound the second-order term in ( 23) by

E          t i=1 ∇ θ l t,i (θ t-1 ) 2 F t-1          ≤[2C 2 l n t + C 2 ∇ (n t -1)n t ] θ t-1 -θ * 2 + 2σ 2 n t ≤ [2C 2 l n t + C 2 ∇ n 2 t 1 {nt >1} ] θ t-1 -θ * 2 + 2σ 2 n t , (25)
following the same steps in the proof of Theorem 1, but with use of (24). Bounding the fourth-order term is a bit heavier computationally, but let us recall that i x i 2 = i x i 2 + i j x i , x j = i x i 2 + 2 i≤ j x i , x j . Then, we have that

n t i=1 ∇ θ l t,i (θ t-1 ) 4 =         n t i=1 ∇ θ l t,i (θ t-1 ) 2 + n t i j ∇ θ l t,i (θ t-1 ), ∇ θ l t, j (θ t-1 )         2 ≤2        n t i=1 ∇ θ l t,i (θ t-1 ) 2        2 + 4         n t i< j ∇ θ l t,i (θ t-1 ), ∇ θ l t, j (θ t-1 )         2 , (26)
as (x + y) 2 ≤ 2x 2 + 2y 2 . For the first term of (26), we have

E                n 1 i=1 ∇ θ l t,i (θ t-1 ) 2        2 F t-1         = n t i=1 E[ ∇ θ l t,i (θ t-1 ) 4 |F t-1 ] + n t i j E[ ∇ θ l t,i (θ t-1 ) 2 ∇ θ l t, j (θ t-1 ) 2 |F t-1 ] ≤8n t [C 4 l θ t-1 -θ * 4 + σ 4 ] + 4n 2 t 1 {nt >1} [C 2 l θ t-1 -θ * 2 + σ 2 ] 2 ,
using the bound from (24), n (n t -1) ≤ n 2 t 1 {nt >1} , and that F t-1 ⊆ F t-1,i ⊂ F t-1, j for all 0 ≤ i < j. To bound the second term of ( 26), we ease notation by denoting ∇ θ l t,i (θ t-1 ) by υ i , giving us

        n t i< j υ i , υ j         2 = n t i< j υ i , υ j 2 + n t i< j,k<l (i, j) (k,l) υ i , υ j υ k , υ l = n t i< j υ i , υ j 2 A + n t i< j,k<l (i, j) (k,l), j=l υ i , υ j υ k , υ l B + n t i< j,k<l (i, j) (k,l), j l υ i , υ j υ k , υ l C .
By Cauchy-Schwarz inequality, we can bound the first term A, by

E[A|F t-1 ] ≤ n t i< j E[ υ i 2 υ j 2 |F t-1 ] ≤ 2n t (n t -1)[C 2 l θ t-1 -θ * 2 + σ 2 ] 2 ≤ 2n 2 t 1 {nt >1} [C 2 l θ t-1 -θ * 2 + σ 2 ] 2 ,
using that F t-1 ⊆ F t-1,i ⊂ F t-1, j for all 0 ≤ i < j. Next, since l = j implies i k, we have

E[B|F t-1 ] = n t i< j,k<l,i k, j=l E[ υ i , υ j υ k , υ l |F t-1 ] = n t i< j,k<l,i k, j=l E[E[ E[υ i |F t-1,i-1 ], υ j E[υ k |F t-1,k-1 ], υ l |F t-1,l-1 ]|F t-1 ] = n t i< j,k<l,i k, j=l E[E[ ∇ θ L(θ t-1 ), υ l 2 |F t-1,l-1 ]|F t-1 ] ≤ n t i< j,k<l,i k, j=l E[ ∇ θ L(θ t-1 ) 2 E[ υ l 2 |F t-1,l-1 ]|F t-1 ] ≤ n t i< j,k<l,i k, j=l 2C 2 ∇ θ t-1 -θ * 2 [C 2 l θ t-1 -θ * 2 + σ 2 ] = n t (n t -1)(n t -2)C 2 ∇ θ t-1 -θ * 2 [C 2 l θ t-1 -θ * 2 + σ 2 ] ≤n 3 t 1 {nt >1} C 2 ∇ θ t-1 -θ * 2 [C 2 l θ t-1 -θ * 2 + σ 2 ],
using the Cauchy-Schwarz inequality and the bound in (24). In the same way, as j l includes (i, j) (k, l), we can rewrite C as

C = n t i< j,k<l, j l υ i , υ j υ k , υ l = n t i< j,k<l,i=k, j l υ i , υ j υ k , υ l C 1 + n t i< j,k<l,i k, j l υ i , υ j υ k , υ l C 2
,

where

E[C 1 |F t-1 ] = E[B|F t-1 ]
. Finally, we can rewrite C 2 as

C 2 = n t i< j,k<l,i k, j l,i=l, j k υ i υ j υ k υ l C 2,1 + n t i< j,k<l,i k, j l,i l, j=k υ i υ j υ k υ l C 2,2 + n t i< j,k<l,i j k l υ i υ j υ k υ l C 2,3
,

where

E[C 2,1 |F t-1 ] = E[C 2,2 |F t-1 ] = E[B|F t-1 ],
and

E[C 2,3 |F t-1 ] = n t i< j,k<l,i j k l E[ ∇ θ L(θ t-1 ) 4 |F t-1 ] ≤ n t (n t -1)(n t -2)(n t -3)C 4 ∇ θ t-1 -θ * 4 ≤ n 4 t 1 {nt >1} C 4 ∇ θ t-1 -θ * 4 .
Thus, the fourth-order term of ( 23), is bounded by

E          n t i=1 ∇ θ l t,i (θ t-1 ) 4 F t-1          ≤16n t [C 4 l θ t-1 -θ * 4 + σ 4 ] + 16n 2 t 1 {nt >1} [C 2 l θ t-1 -θ * 2 + σ 2 ] 2 + 12n 3 t 1 {nt >1} C 2 ∇ θ t-1 -θ * 2 [C 2 l θ t-1 -θ * 2 + σ 2 ] + 4n 4 t 1 {nt >1} C 4 ∇ θ t-1 -θ * 4 ≤[16C 4 l n t + 16C 4 l n 2 t 1 {nt >1} + 12C 2 ∇ C 2 l n 3 t 1 {nt >1} + 4C 4 ∇ n 4 t 1 {nt >1} ] θ t-1 -θ * 4 + [32C 2 l σ 2 n 2 t 1 {nt >1} + 12C 2 ∇ σ 2 n 3 t 1 {nt >1} ] θ t-1 -θ * 2 + 16σ 4 n t + 16σ 4 n 2 t 1 {nt >1} . ( 27
)
the bound from ( 25) and ( 27) into [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF], we obtain the recursive relation for the fourth-order moment:

E[ θ t -θ * 4 |F t-1 ] ≤[1 -4µγ t + 8C 2 ∇ 1 {nt >1} γ 2 t + 16C 2 l n -1 t γ 2 t + 48C 4 l n -3 t γ 4 t + 48C 4 l n -2 t 1 {nt >1} γ 4 t + 36C 2 ∇ C 2 l n -1 t 1 {nt >1} γ 4 t + 12C 4 ∇ 1 {nt >1} γ 4 t ] θ t-1 -θ * 4 + [16σ 2 n -1 t γ 2 t + 96C 2 l σ 2 n -2 t 1 {nt >1} γ 4 t + 36C 2 ∇ σ 2 n -1 t 1 {nt >1} γ 4 t ] θ t-1 -θ * 2 + 48σ 4 n -3 t γ 4 t + 48σ 4 n -2 t 1 {nt >1} γ 4 t .
Using the Young's inequalities, 2C 2

∇ C 2 l ≤ n t C 4 ∇ + n -1 t C 4 l , 16σ 2 n -1 t γ 2 t θ t-1 -θ * 2 ≤ 2µγ t θ t -θ * 4 + 32σ 4 µ -1 n -2 t γ 3 t , 2C 2 l σ 2 n -2 t 1 {nt >1} γ 4 t θ t-1 -θ * 2 ≤ C 4 l n -2 t 1 {nt >1} γ 4 t θ t -θ * 4 +σ 4 n -2 t 1 {nt >1} γ 4 t , and 2C 2 ∇ σ 2 n -1 t 1 {nt >1} γ 4 t θ t-1 -θ * 2 ≤ C 4 ∇ 1 {nt >1} γ 4 t θ t - θ * 4 + σ 4 n -2 t 1 {nt >1} γ 4 t , yields, E[ θ t -θ * 4 |F t-1 ] ≤[1 -2µγ t + 8C 2 ∇ 1 {nt >1} γ 2 t + 16C 2 l n -1 t γ 2 t + 48C 4 l n -3 t γ 4 t + 114C 4 l n -2 t 1 {nt >1} γ 4 t + 48C 4 ∇ 1 {nt >1} γ 4 t ] θ t-1 -θ * 4 + 32µ -1 σ 4 n -2 t γ 3 t + 48σ 4 n -3 t γ 4 t + 114σ 4 n -2 t 1 {nt >1} γ 4 t . (28)
Taking, the expectation on both sides of the inequality in ( 28) yields the recursive relation for the fourth-order moment:

∆ t ≤[1 -2µγ t + 8C 2 ∇ 1 {nt >1} γ 2 t + 16C 2 l n -1 t γ 2 t + 48C 4 l n -3 t γ 4 t + 114C 4 l n -2 t 1 {nt >1} γ 4 t + 48C 4 ∇ 1 {nt >1} γ 4 t ]∆ t-1 + 32µ -1 σ 4 n -2 t γ 3 t + 48σ 4 n -3 t γ 4 t + 114σ 4 n -2 t 1 {nt >1} γ 4 t . ( 29
)
with

∆ t = E[ θ t -θ * 4
] for some ∆ 0 ≥ 0. By Proposition 5, we achieve the (upper) bound of ∆ t in ( 29), given as

∆ t ≤ exp         -µ t i=t/2 γ i         Π ∆ t + 32σ 4 µ 2 max t/2≤i≤t γ 2 i n 2 i + 48σ 4 µ max t/2≤i≤t γ 3 i n 3 i + 114σ 4 µ max t/2≤i≤t γ 3 i 1 {ni>1} n 2 i .
where

Π ∆ t is given by exp        32C 2 l t i=1 γ 2 i n i        exp        96C 4 l t i=1 γ 4 i n 3 i        exp        228C 4 l t i=1 1 {ni>1} γ 4 i n 2 i        exp        16C 2 ∇ t i=1 1 {ni>1} γ 2 i        exp        96C 4 ∇ t i=1 1 {ni>1} γ 4 i              ∆ 0 + 2σ 4 C 4 l + 4σ 4 γ 1 µC 2 l n 1       , (30)
with use of

max 1≤i≤t 32µ -1 σ 4 n -2 i γ i + 48σ 4 n -3 i γ 2 i + 114σ 4 n -2 i 1 {ni>1} γ 2 i 8C 2 ∇ 1 {ni>1} + 16C 2 l n -1 i + 48C 4 l n -3 i γ 2 i + 114C 4 l n -2 i 1 {ni>1} γ 2 i + 48C 4 ∇ 1 {ni>1} γ 2 i ≤ σ 4 C 4 l + 2σ 4 γ 1 µC 2 l n 1 .
At last, bounding the projected estimate (3) follows from that

E[ P Θ (θ) -θ * 2 ] ≤ E[ θ -θ * 2 ], ∀θ ∈ Θ.

Proofs for Section 4.1

Proof of Theorem 2. Following Polyak and Juditsky [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], we rewrite (2) to

θ t = θ t-1 - γ t n t n t i=1 ∇ θ l t,i (θ t-1 ) ⇐⇒ 1 γ t (θ t-1 -θ t ) = ∇ θ l t (θ t-1 ), ( 31
)
where ∇ θ l t (θ t-1 ) denotes n -1

t n t i=1 ∇ θ l t,i (θ t-1 ). Note ∇ θ l t (θ t-1 ) ≈ ∇ θ l t (θ * ) + ∇ 2 θ l t (θ * )(θ t-1 -θ * )
, and that ∇ θ l t (θ * ) and ∇ θ l t (θ)-∇ θ L(θ) behaves almost like an i.i.d. sequences with zero mean. Thus, θt -θ * behaves like

-∇ θ L(θ * ) -1 N -1 t t i=1 n i ∇ θ l i (θ * ) leading to a bound in O( √ N t ). Observe that ∇ 2 θ L(θ * )(θ t-1 -θ * ) =∇ θ l t (θ t-1 ) -∇ θ l t (θ * ) -[∇ θ l t (θ t-1 ) -∇ θ l t (θ * ) -∇ θ L(θ t-1 )] martingale term -[∇ θ L(θ t-1 ) -∇ 2 θ L(θ * )(θ t-1 -θ * )] rest term
, where ∇ 2 θ L(θ * is invertible with lowest eigenvalue greater than µ, i.e., ∇ 2 θ L(θ * ) ≥ µ. Thus, summing the parts and using the Minkowski's inequality, we obtain the inequality:

E θt -θ * 2 1 2 ≤          E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ * ) 2                   1 2 +          E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) 2                   1 2 +          E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L (θ i-1 ) 2                   1 2 +          E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ L (θ i-1 ) -∇ 2 θ L (θ * ) (θ i-1 -θ * ) 2                   1 2
.

As (∇ θ l t,i (θ * )) is a square-integrable martingale increment sequences on R d (Assumption 1), we have

E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ * ) 2          ≤ 1 N 2 t t i=1 n i j=1 E ∇ 2 θ L (θ * ) -1 ∇ θ l i, j (θ * ) 2 ≤ Tr ∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 N t , (32)
using Assumption 4. To ease notation, we denote Tr[∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ] by Λ. Next, note that for all t ≥ 1, we have the relation in [START_REF] Sridharan | Fast rates for regularized objectives[END_REF], giving us

1 N t t i=1 n i ∇ θ l i (θ i-1 ) = 1 N t t i=1 n i γ i (θ i-1 -θ i ) = 1 N t t-1 i=1 (θ i -θ * ) n i+1 γ i+1 - n i γ i - 1 N t (θ t -θ * ) n t γ t + 1 N t (θ 0 -θ * ) n 1 γ 1 , leading to ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) ≤ 1 N t µ t-1 i=1 θ i -θ * n i+1 γ i+1 - n i γ i + 1 N t µ θ t -θ * n t γ t + 1 N t µ θ 0 -θ * n 1 γ 1 .
Hence, with the notion of δ t = E[ θ t -θ * 2 ] this expression can be simplified to

         E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) 2                   1 2 ≤ 1 N t µ t-1 i=1 δ 1 2 i n i+1 γ i+1 - n i γ i + n t N t γ t µ δ 1 2 t + n 1 N t γ 1 µ δ 1 2 0 . (33)
For the martingale term, we have

E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L (θ i-1 ) 2          ≤ 1 N 2 t µ 2 t i=1 n 2 i E ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) 2 = 1 N 2 t µ 2 t i=1 E           n i j=1 ∇ θ l i, j (θ i-1 ) -∇ θ l i, j (θ * ) 2           ≤ 1 N 2 t µ 2 t i=1 n i j=1 E ∇ θ l i, j (θ i-1 ) -∇ θ l i, j (θ * ) 2 1 2 ≤ C 2 l N 2 t µ 2 t i=1 n i δ i-1 , (34)
by the Cauchy-Schwarz inequality and Assumption 2-p. For all t ≥ 1, the rest term is directly bounded by [START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF]:

         E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ L (θ i-1 ) -∇ 2 θ L (θ * ) (θ i-1 -θ * ) 2                   1 2 ≤ C ∇ N t µ t i=1 n i ∆ 1 2 i-1 , (35)
with the notion

∆ t = E[ θ t -θ * 4 ]
. Finally, combining the terms from (32) to [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], gives us

δ1/2 t ≤ Λ 1/2 N 1/2 t + 1 N t µ t-1 i=1 δ 1/2 i n i+1 γ i+1 - n i γ i + n t N t γ t µ δ 1/2 t + n 1 N t γ 1 µ δ 1/2 0 + C l N t µ        t i=1 n i δ i-1        1/2 + C ∇ N t µ t i=1 n i ∆ 1/2 i-1 , (36)
where δt = E[ θtθ * 2 ], which can be simplified into (12) by shifting the indices and collecting the δ 0 terms.

of Corollary 3. As n t = C ρ for all t ≥ 1, we simplify the bound for δt in [START_REF] Haldane | Note on the median of a multivariate distribution[END_REF] to

δ1/2 t ≤ Λ 1/2 N 1/2 t + C ρ N t µ t-1 i=1 δ 1/2 i 1 γ i+1 - 1 γ i + C ρ N t γ t µ δ 1/2 t + C ρ N t µ 1 γ 1 + C l δ 1/2 0 + C l C 1 2 ρ N t µ         t-1 i=1 δ i         1/2 + C ∇ C ρ N t µ t-1 i=0 ∆ 1/2 i . (37)
The second-order moment δ t is bounded by Corollary 1 but with use of (20) as we work in terms of t. The fourth-order moment ∆ t from Lemma 1 can be simplified to:

∆ t ≤ exp         -µ t i=t/2 γ i         Π c ∞ + 1 µ       32σ 4 µC 2 ρ max t/2≤i≤t γ 2 i + 48σ 4 C 3 ρ max t/2≤i≤t γ 3 i + 114σ 4 1 {Cρ>1} C 2 ρ max t/2≤i≤t γ 3 i       ≤ exp        - µC γ C β ρ t 1-α 2 1-α        Π c ∞ + 1 µ         2 2α 32σ 4 C 2 γ C 2β ρ µC 2 ρ t 2α + 2 3α 48σ 4 C 3 γ C 3β ρ C 3 ρ t 3α + 2 3α 114σ 4 C 3 γ C 3β ρ 1 {Cρ>1} C 2 ρ t 3α         , using that γ t = C γ C β ρ t -α is decreasing as α ∈ (1/2, 1
). Regarding Π ∆ t defined in [START_REF] Shalev-Shwartz | Online learning and online convex optimization[END_REF], we can bound it by

Π c ∞ = exp         64αC 2 l C 2 γ C 2β ρ (2α -1) C ρ         exp         (192 + 456C ρ 1 {Cρ>1} )C 4 l C 4 γ C 4β ρ C 3 ρ         exp         32αC 2 ∇ C 2 γ C 2β ρ 1 {Cρ>1} 2α -1         exp 192C 4 ∇ C 4 γ C 4β ρ 1 {Cρ>1}        ∆ 0 + 2σ 4 C 4 l + 4σ 4 C γ µC 2 l C 1-β ρ        , using t i=1 i -2α ≤ 2α/(2α -1) and t i=1 i -4α ≤ 2. Note that Π c ∞ is a finite constant, independent of t.
To bound the first term of (37), namely

C ρ N t µ t-1 i=1 δ 1/2 i |γ -1 i+1 -γ -1 i |, we remark that |γ -1 t+1 -γ -1 t | ≤ C -1 γ C -β ρ αt α-1 , one has (since √ a + b ≤ √ a + √ b), C ρ N t µ t-1 i=1 δ 1 2 i 1 γ i+1 - 1 γ i ≤ C 1-β ρ α C γ µN t t i=1 i α-1           exp        - µC γ C β ρ i 1-α 2 2-α        π c ∞ + 2 1+α 2 σ C γ √ µC 1-β 2 ρ i α/2           . ( 38
)
For simplicity, let us denote

A c ∞ = ∞ i=0 exp        - µC γ C β ρ i 1-α 2 2-α        ≥ ∞ i=0 i α-1 exp        - µC γ C β ρ i 1-α 2 2-α        ,
as α < 1. Thus, the first part of (38) is bounded as follows:

C 1-β ρ α √ π c ∞ C γ µN t t i=1 i α-1 exp        - µC γ C β ρ i 1-α 2 2-α        ≤ C 1-β ρ α √ π c ∞ A c ∞ C γ µN t .
Furthermore, with the help of an integral test for convergence, one has

t i=1 i α/2-1 ≤ 1 + t 1 s α/2-1 ds = 1 + (2/α)t α/2 - (2/α) ≤ (2/α)t α/2
, such that the second part of (38) can be bounded by

2 1+α 2 σC 1-β 2 ρ α C 1/2 γ µ 3/2 N t t i=1 i α/2-1 ≤ 2 3+α 2 σC 1-β 2 ρ t α/2 C 1/2 γ µ 3/2 N t = 2 3+α 2 σC 1-α-β 2 ρ C 1/2 γ µ 3/2 N 1-α/2 t .
By combining this, we get

C ρ N t µ t-1 i=1 δ 1 2 i 1 γ i+1 - 1 γ i ≤ C 1-β ρ α √ π c ∞ A c ∞ C γ µN t + 2 3+α 2 σC 1-α-β 2 ρ C γ µ 3/2 N 1-α/2 t . (39)
Similarly, second term of (37), can bounded by

C ρ N t γ t µ δ 1 2 t ≤ C 1-α-β ρ C γ µN 1-α t           exp        - µC γ C β ρ t 1-α 2 2-α        π c ∞ + 2 1+α 2 σ C γ √ µC 1-β 2 ρ t α/2           ≤ C 2-α-β ρ √ π c ∞ A c ∞ C γ µN 2-α t + 2 1+α 2 C 1-α-β 2 ρ σ C γ µ 3/2 N 1-α/2 t , using exp(-µC γ C β ρ t 1-α /2 2-α ) = A c t ≤ t -1 t i=1 A c i ≤ t -1 A c ∞ as A c
t is decreasing. In a same way, one has

C l C 1 2 ρ N t µ         t-1 i=1 δ i         1 2 ≤ C l C 1 2 ρ N t µ        A c ∞ π c ∞ + 2 1+α σ 2 C γ t 1-α (1 -α) µC 1-β ρ        1/2 ≤ C l C 1 2 ρ √ π c ∞ √ A c ∞ N t µ + 2 1+α 2 C l σ C γ C 1-α-β 2 ρ µ 3/2 N 1+α 2 t .
Bound the last term of (37), is done as follows,

C ∇ C ρ N t µ t-1 i=0 ∆ 1 2 i ≤ C ∇ C ρ N t µ t-1 i=0 exp        - µC γ C β ρ i 1-α 2 2-α        Π c ∞ + 2 α 6C ∇ σ 2 C γ C β ρ N t µ 2 t-1 i=1 i -α + (6 + 71 {Cρ>1} )2 3α/2 C ∇ σ 2 C 3/2 γ C 3β/2 ρ N t µ 3/2 t-1 i=1 i -3α/2 ≤ C ∇ C ρ √ Π c ∞ A c ∞ N t µ + 2 α 6C ∇ σ 2 C γ C 1-α-β ρ µ 2 N α t + (6 + 71 {Cρ>1} )2 3α/2 C ∇ σ 2 C 3/2 γ C 3β/2 ρ ψ 0 3α/2 (N t /C ρ ) µ 3/2 N t .
Thus, by collecting the terms above, we obtain:

δ1/2 t ≤ Λ 1/2 N 1/2 t + 6σC 1-α-β 2 ρ C γ µ 3/2 N 1-α/2 t + 2 α 6C ∇ σ 2 C γ C 1-α-β ρ µ 2 N α t + C 2-α-β ρ √ π c ∞ A c ∞ C γ µN 2-α t + 2 1+α 2 C l σ C γ C 1-α-β 2 ρ µ 3/2 N 1+α 2 t + C ρ Γ c µN t + (6 + 71 {Cρ>1} )2 3α/2 C ∇ σ 2 C 3/2 γ C 3β/2 ρ µ 3/2 ψ 0 3α/2 (N t /C ρ ) -1 N t ,
where

Γ c = (1/C γ C β ρ + C l )δ 1/2 0 + C l √ π c ∞ A c ∞ /C 1/2 ρ + √ π c ∞ A c ∞ /C γ C β ρ + C ∇ √ Π c ∞ A c ∞ .
Proof of Corollary 4. The steps of the proof follows the ones of Corollary 3 with the smart notation of φ and ρ: The bound for δt in [START_REF] Haldane | Note on the median of a multivariate distribution[END_REF] is given by

δ1/2 t ≤ Λ 1/2 N 1/2 t + 1 N t µ t-1 i=1 δ 1/2 i n i+1 γ i+1 - n i γ i + n t N t γ t µ δ 1/2 t + n 1 N t µ 1 γ 1 + C l δ 1/2 0 + C l N t µ         t-1 i=1 n i+1 δ i         1/2 + C ∇ N t µ t-1 i=0 n i+1 ∆ 1/2 i , (40)
where the learning rate is on the form γ t = C γ n β t t -α with n t = C ρ t ρ . The second-order moment δ t is upper bounded by ( 21) from Corollary 2. The fourth-order moment ∆ t from Lemma 1 can be simplified as follows,

∆ t ≤ exp         -µ t i=t/2 γ i         Π v ∞ + 32σ 4 µ 2 max t/2≤i≤t γ 2 i n 2 i + 162σ 4 µ max t/2≤i≤t γ 3 i n 2 i ,
as n t ≥ 1 for any t ≥ 1 and β ≤ 1, and

Π v ∞ = exp         32(α -β ρ)C 2 γ C 2β ρ (2C 2 l + C 2 ∇ ) 2(α -β ρ) -1         exp 192C 4 γ C 4β ρ (4C 4 l + C 4 ∇ )        ∆ 0 + 2σ 4 C 4 l + 4σ 4 C γ µC 2 l C 1-β ρ ,       
using that t i=1 i -a ≤ 2 for a ≥ 2. Next, for ρ ≥ 0, we have

∆ t ≤ exp        - µC γ C β ρ t 1+βρ-α 2 1+βρ-α        Π v ∞ + 2 2α-2βρ+2ρ 32σ 4 C 2 γ C 2β ρ µ 2 C 2 ρ t 2α-2βρ+2ρ + 2 3α-3βρ+2ρ 162σ 4 C 3 γ C 3β ρ µC 2 ρ t 3α-3βρ+2ρ
, that αβρ ∈ (1/2, 1). If ρ < 0, one directly have

∆ t ≤ exp        - µC γ C β ρ t 1-α 2 1-α        Π v ∞ + 2 2α 32σ 4 C 2 γ C 2β ρ µ 2 t 2α + 2 3α 162σ 4 C 3 γ C 3β ρ µt 3α .
With the notion of φ and ρ, we can combine the two ρ-cases as follows:

∆ t ≤ exp         - µC γ C β1 {ρ≥0} ρ t (1-φ)(1+ ρ) 2 (1-φ)(1+ ρ)         Π v ∞ + 2 2φ(1+ ρ) 32σ 4 C 2 γ C 2β ρ µ 2 C 21 {ρ≥0} ρ t 2φ(1+ ρ) + 2 3φ(1+ ρ)-ρ162σ 4 C 3 γ C 3β ρ µC 21 {ρ≥0} ρ t 3φ(1+ ρ)-ρ .
We will in the following bound the terms for t but afterwards we will translate it to terms in

N t . If ρ ≥ 0, the first relation is t ≥ (N t /2C ρ ) 1/(1+ρ) since N t = C ρ (t ρ + t-1 i=1 i ρ ) ≤ C ρ (t ρ + t 1 x ρ dx) ≤ C ρ (t ρ +t ρ t 1 dx) = C ρ (t ρ +t 1+ρ ) ≤ 2C ρ t 1+ρ by use the integral test for convergence. Similarly, N t = C ρ t i=1 i ρ ≥ C ρ t 0 x ρ dx = C ρ t ρ+1 , thus, t ≤ (N t /C ρ ) 1/(1+ρ) . If ρ < 0, one has t ≤ N t and N t ≤ C ρ t, i.e., t ≥ N t /C ρ .
Bounding 1

N t µ t-1 i=1 δ 1/2 i |n i+1 /γ i+1 -n i /γ i |, we first note n t /γ t = C -1 γ C 1-β ρ t (1-β
)ρ+α for ρ ≥ 0. Thus, by the mean value theorem, we obtain:

n i+1 γ i+1 - n i γ i ≤ ((1 -β)ρ + α) C 1-β ρ C γ sup ν∈(i,i+1) ν (1-β)ρ+α-1 ≤ ((1 -β)ρ + α) C 1-β ρ C γ i 1-(1-β)ρ-α , (41)
as α + (1 -β)ρ ≤ 1 -ρ since α -βρ ∈ (1/2, 1
). For ρ < 0, the mean value theorem gives us

n i+1 γ i+1 - n i γ i = 1 C γ n 1-β i+1 (i + 1) α -n 1-β i i α ≤ C 1-β ρ C γ |(i + 1) α -i α | ≤ αC 1-β ρ C γ sup ν∈(i,i+1) ν α-1 ≤ αC 1-β ρ C γ i 1-α ,
as (n t ) t≥1 is a decreasing sequence and β ≤ 1. Thus, for any ρ ∈ (-1, 1), we have

n i+1 γ i+1 - n i γ i ≤ φ(1 + ρ)C 1-β ρ C γ i 1-φ(1+ ρ) .
By using this, we obtain:

1 N t µ t-1 i=1 δ 1 2 i n i+1 γ i+1 - n i γ i ≤ φ(1 + ρ)C 1-β ρ N t µC γ t i=1 i φ(1+ ρ)-1           exp         - µC γ C β1 {ρ≥0} ρ i (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         π v ∞ + 2 1+φ(1+ ρ) 2 σ C γ √ µC (1-β) 2 1 {ρ≥0} ρ i φ(1+ ρ) 2           . Next, let us denote A v ∞ = ∞ i=0 i ρ exp         - µC γ C β1 {ρ≥0} ρ i (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         ≥ ∞ i=0 i φ(1+ ρ)-1 exp         - µC γ C β1 {ρ≥0} ρ i (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         , since φ(1 + ρ) -1 = α + (1 -β) ρ -1 ≤ ρ. Thus, φ(1 + ρ)C 1-β ρ √ π v ∞ N t µC γ t i=1 i φ(1+ ρ)-1 exp         - µC γ C β1 {ρ≥0} ρ i (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         ≤ φ(1 + ρ)C 1-β ρ √ π v ∞ A v ∞ N t µC γ .
Furthermore, with the help of an integral test for convergence, we have

φ(1 + ρ)2 1+φ(1+ ρ) 2 σC 1-β 2 1 {ρ≥0} ρ µ 3/2 C γ N t t i=1 i φ(1+ ρ) 2 -1 ≤ 2 3+φ(1+ ρ) 2 σC 1-β 2 1 {ρ≥0} ρ t φ(1+ ρ) 2 µ 3/2 C γ N t ≤ 2 3+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t .
Summarising, with use of φ(1 + < 2, we obtain

1 N t µ t-1 i=1 δ 1 2 i n i+1 γ i+1 - n i γ i ≤ φ(1 + ρ)C 1-β ρ √ π v ∞ A v ∞ N t µC γ + 2 3+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t ≤ 2C 1-β ρ √ π v ∞ A v ∞ µC γ N t + 2 3+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t .
Similarly, for n t N t γ t µ δ 1/2 t , one have

n t N t γ t µ δ 1 2 t ≤ C 1-β ρ √ π v ∞ t φ(1+ ρ) N t C γ µ exp         - µC γ C β1 {ρ≥0} ρ t (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         + 2 1+φ(1+ ρ) 2 σC 1-β 2 1 {ρ≥0} ρ t φ(1+ ρ) 2 µ 3/2 C γ N t ≤ C 2-φ-β ρ √ π v ∞ A v ∞ µC γ N 2-φ t + 2 1+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t .
For n 1 N t µ (γ -1 1 + C l )δ 1/2 0 , we insert the definition of our learning functions, giving us

n 1 N t µ 1 γ 1 + C l δ 1/2 0 = C ρ N t µ        1 C γ C β ρ + C l        δ 1/2 0 .
Bounding C l N t µ ( t-1 i=1 n i+1 δ i ) 1/2 , follows the ideas from above, using that n t+1 ≤ 2 ρn t , to obtain

C l N t µ         t-1 i=1 n i+1 δ i                 C ρ π v ∞ A v ∞ + 2 φ(1+ ρ) σ 2 C γ C β1 {ρ≥0} ρ t (1-φ)(1+ ρ) µ         1 2 ≤ 2 ρ/2 C l C ρ √ π v ∞ √ A v ∞ µN t + 2 φ(1+ ρ) 2 C l σ C γ C β/21 {ρ≥0} ρ t (1-φ)(1+ ρ) 2 µ 3/2 N t ≤ 2 ρ/2 C l C ρ √ π v ∞ √ A v ∞ µN t + 2 φ(1+ ρ) 2 C l σ C γ µ 3/2 C 1-φ-β 2 1 {ρ≥0} ρ N 1+φ 2 t .
Likewise, for

C ∇ N t µ t-1 i=0 n i+1 ∆ 1/2 i , we get C ∇ N t µ t-1 i=0 n i+1 ∆ 1 2      - µC γ C β1 {ρ≥0} ρ i (1-φ)(1+ ρ) 2 1+(1-φ)(1+ ρ)         Π v ∞ + 2 φ(1+ ρ) 6σ 2 C γ C β ρ µC 1 {ρ≥0} ρ i φ(1+ ρ) + 2 3φ(1+ ρ)/2-ρ/2 13σ 2 C 3/2 γ C 3β/2 ρ µ 1/2 C 1 {ρ≥0} ρ i 3φ(1+ ρ)/2         ≤ 2 ρC ∇ C ρ √ Π v ∞ A v ∞ µN t + 2 φ(1+ ρ)+ ρC ∇ σ 2 C γ C 1+β ρ µ 2 C 1 {ρ≥0} ρ N t t-1 i=1 i β ρ-α + 2 3φ(1+ ρ)/2+ ρ/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ µ 3/2 C 1 {ρ≥0} ρ N t t-1 i=1 i 3(β ρ-α)/2 ,
where the second term can be bounded as

2 (1+φ)(1+ ρ)-1 C ∇ σ 2 C γ C 1+β ρ µ 2 C 1 {ρ≥0} ρ N t t-1 i=1 i β ρ-α ≤ 2 (1+φ)(1+ ρ)-1 C ∇ σ 2 C γ C 1+β ρ t 1+β ρ-α (1 + β ρ -α)µ 2 C 1 {ρ≥0} ρ N t ≤ 2 (1+φ)(1+ ρ)-2 C ∇ σ 2 C γ µ 2 C 1-φ-β ρ N φ t ,
the third term by

2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ µ 3/2 C 1 {ρ≥0} ρ N t t-1 i=1 i 3(β ρ-α)/2 ≤ 2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ ψ ρ 3(α-β ρ)/2 (N t /C ρ ) µ 3/2 C 1 {ρ≥0} ρ N t .
By collecting these bounds, we get

C ∇ N t µ t-1 i=0 n i+1 ∆ 1 2 i ≤ 2 ρC ∇ C ρ √ Π v ∞ A v ∞ µN t + 2 (1+φ)(1+ ρ)-2 C ∇ σ 2 C γ µ 2 C 1-φ-β ρ N φ t + 2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ ψ ρ 3(α-β ρ)/2 (N t /C ρ ) µ 3/2 C 1 {ρ≥0} ρ N t .
Combining our findings from above, we have

δ1/2 t ≤ Λ 1/2 N 1/2 t + 2C 1-β ρ √ π v ∞ A v ∞ µC γ N t + 2 3+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t + C 2-φ-β ρ √ π v ∞ A v ∞ µC γ N 2-φ t + 2 1+φ(1+ ρ) 2 σC 1-φ-β 2 1 {ρ≥0} ρ µ 3/2 C γ N 1-φ/2 t + C ρ N t µ        1 C γ C β ρ + C l        δ 1 2 0 + 2 ρ/2 C l C ρ √ π v ∞ √ A v ∞ µN t + 2 φ(1+ ρ) 2 C l σ C γ µ 3/2 C 1-φ-β 2 1 {ρ≥0} ρ N 1+φ 2 t + 2 ρC ∇ C ρ √ Π v ∞ A v ∞ µN t + 2 (1+φ)(1+ ρ)-2 C ∇ σ 2 C γ µ 2 C 1-φ-β ρ N φ t + 2 3(1+φ)(1+ ρ)/2 C ∇ σ 2 C 3/2 γ C 1+3β/2 ρ ψ ρ 3(α-β ρ)/2 (N t /C ρ ) µ 3/2 C 1 {ρ≥0} ρ N t .
This can be simplified to the desired using Γ v given by ( 1

/C γ C β ρ + C l )δ 1/2 0 + 2 ρC l √ π v ∞ A v ∞ /C 1/2 ρ + 2 √ π v ∞ A v ∞ /C γ C β ρ + 2 ρC ∇ √ Π v ∞ A v ∞ , consisting of the finite constants π v ∞ , Π v ∞ and A v ∞ . 7
t ≤ Λ 1/2 N 1/2 t + 1 N t µ t-1 i=1 n i+1 γ i+1 - n i γ i δ 1/2 i + n t N t γ t µ δ 1/2 t + n 1 N t µ 1 γ 1 + C l δ 1/2 0 + C l N t µ         t-1 i=1 n i+1 δ i         1/2 + C ∇ N t µ t i=0 n i+1 ∆ 1/2 i where Λ = Tr(∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ) and C ∇ = C ∇ + 2 2 G Θ /D 2 θ . Proof of Theorem 3. Denote E[ θt -θ * 2 ]
by δt with ( θt ) given by ( 4) using (θ t ) from [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF]. As in the proof Theorem 2, we follow the steps of Polyak and Juditsky [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], in which, we can rewrite (3) to

1 γ t (θ t-1 -θ t ) = ∇ θ l t (θ t-1 ) - 1 γ t Ω t , where ∇ θ l t (θ t-1 ) = n -1 t n t i=1 ∇ θ l t,i (θ t-1 ) and Ω t = P Θ (θ t-1 -γ t ∇ θ l t (θ t-1 )) -(θ t-1 -γ t ∇ θ l t (θ t-1 )
). Thus, summing the parts, using the Minkowski's inequality, and bounding each term gives us the same bound as in Theorem 2, but with an additional term regarding Ω t , namely

         E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i γ i Ω i 2                   1 2 ≤ 1 µN t t i=1 n i γ i E Ω i 2 = 1 µN t t i=1 n i γ i E Ω i 2 1 {θ i-1 -γ i ∇ θ l i (θ i-1 ) Θ} , (42)
using Godichon-Baggioni [9, Lemma 4.3]. Next, we note that

E[ Ω t 2 1 {θ t-1 -γt ∇ θ lt (θ t-1 ) Θ} ] = 4γ 2 t G 2 Θ P[θ t-1 -γ t ∇ θ l t (θ t-1 ) Θ], since Ω t 2 ≤2 P Θ (θ t-1 -γ t ∇ θ l t (θ t-1 )) -θ t-1 2 + 2γ 2 t ∇ θ l t (θ t-1 ) 2 = 2 P Θ (θ t-1 -γ t ∇ θ l t (θ t-1 )) -P Θ (θ t-1 ) 2 + 2γ 2 t ∇ θ l t (θ t-1 ) 2 ≤2 θ t-1 -γ t ∇ θ l t (θ t-1 ) -θ t-1 2 + 2γ 2 t ∇ θ l t (θ t-1 ) 2 = 4γ 2 t ∇ θ l t (θ t-1 ) 2 ≤ 4γ 2 t G 2 Θ ,
as P Θ is Lipschitz and θ l t,i (θ) 2 ≤ G 2 Θ for any θ ∈ Θ. Moreover, as in Godichon-Baggioni and Portier [10, Theorem 4.2], we know that P[θ t-1 -γ t ∇ θ l t (θ t-1 ) Θ] ≤ ∆ t /D 4 θ , where D θ = inf θ∈∂Θ θθ * with ∂Θ denoting the frontier of Θ. Thus, (42) can then be bounded by

1 µN t t i=1 n i γ i E Ω i 2 1 {θ i-1 -γ i ∇ θ l i (θ i-1 ) Θ} ≤ 2G Θ µD 2 θ N t t i=1 n i ∆ 1/2 i ≤ 2 2 G Θ µD 2 θ N t t i=1 n i+1 ∆ 1/2 i ,
using that the sequence (n t ) is either constant or time-varying, meaning n t+1 /n t ≤ 2.

Proof of Corollary 5. The proof follows directly from Corollary 3 but with use of Theorem 3.

Proof of Corollary 6. The proof follows directly from Corollary 4 but with use of Theorem 3.

A. Technical Proofs

Appendix A contains purely technical results used in the proofs presented in Section 7. In what follows, we use the convention inf ∅ = 0, 0 t=1 = 0, and 0 t=1 = 1. Proposition 1. Let (γ t ) t≥1 be a positive sequence. For any k ≤ t, and ω > 0, we have

t i=k t j=i+1 1 + ωγ j γ i ≤ 1 ω t j=k 1 + ωγ j ≤ 1 ω exp         ω t j=k γ j         . (A.1)
Proof of Proposition 1. We begin with considering the first inequality in (A.1), which follows by expanding the sum of product: As the (positive) terms cancel out, we end up with the first inequality in (A.1):

1 ω t i=k         t j=i 1 + ωγ j - t j=i+1 1 + ωγ j         = 1 ω         t j=k 1 + ωγ j - t j=k+1 1 + ωγ j + • • • - t j=t+1 1 + ωγ j         = 1 ω         t j=k 1 + ωγ j - t j=t+1 1 + ωγ j         = 1 ω         t j=k 1 + ωγ j -1         ≤ 1 ω t j=k 1 + ωγ j ,
as t t+1 = 1 for all t ∈ N. Using the (simple) bound of 1 + t ≤ exp(t) for all t ∈ R, we obtain the second inequality of (A.1): 1 -ωγ j 1 -

ωγ i -1 = - 1 ω t i=k         t j=i+1 1 -ωγ j 1 -ωγ i - t j=i+1 1 -ωγ j         = - 1 ω t i=k         t j=i 1 -ωγ j - t j=i+1 1 -ωγ j         = 1 ω t i=k         t j=i+1 1 -ωγ j - t j=i 1 -ωγ j         .
As we only have positive terms, we can upper bound the term: with δ 0 ≥ 0 and ω > 0. Denote t 0 = inf {t ≥ 1 : η t ≤ ω}, and suppose that for all t ≥ t 0 + 1, one has ωγ t ≤ 1. Then, for γ t and η t decreasing, we have the upper bound on (δ t ):

1 ω t i=k         t j=i+1 1 -ωγ j - t j=i 1 -ωγ j         ≤ 1 ω         1 - t j=k 1 -ωγ j         ≤ 1
δ t ≤ exp         -ω t i=t/2 γ i                 exp        t 0 i=1 η i γ i        δ 0 + max 1≤i≤t 0 ν i η i + t/2-1 i=t 0 +1 ν i γ i         + 1 ω max t/2≤i≤t ν i , (A.6)
for all t ∈ N with the convention that t/2 t 0 = 0 if t/2 < t 0 . Proof of Proposition 4. Applying the recursive relation from (A.5) t times, we derive:

δ t ≤ t i=1
A t,2 = t 0 i=1 t j=i+1 1 -2ωγ j + η j γ j ν i γ i =         t j=t 0 +1 1 -2ωγ j + η j γ j         t 0 i=1 t 0 j=i+1 1 -2ωγ j + η j γ j ν i γ i ≤         t j=t 0 +1 1 -ωγ j        

Figure 1 :

 1 Figure 1: Linear regression for various data streams n t = C ρ t ρ . See Section 5.1 for details.

  (a) Constant streaming batches, ρ = 0, β = 0 (b) Time-varying streaming batches, C ρ = 1, β = 0 (c) Time-varying streaming batches, C ρ = 8, β = 0 (d) Time-varying streaming batches, C ρ = 64, β = 0 (e) Time-varying streaming batches, C ρ = 128, β = 0 (f) Time-varying streaming batches, C ρ = 8, β = 1/3

Figure 2 :

 2 Figure 2: Geometric median for various data streams n t = C ρ t ρ . See Section 5.2 for details. (a) Constant streaming batches, ρ = 0, β = 0 (b) Time-varying streaming batches, C ρ = 1, β = 0

Figure 3 :

 3 Figure 3: WASSG for various data streams n t = C ρ t ρ . See Section 6 for details.(a) Linear regression, time-varying streaming batches, C ρ = 8, β = 1/3 (b) Geometric median, time-varying streaming batches, C ρ = 8, β = 1/3

1 + 1 +

 11 ωγ j 1 + ωγ iωγ j 1 + ωγ i -

  Let (γ t ) t≥1 be a positive sequence. Let ω > 0 and k ≤ t such that for all i ≥ k, ωγ i ≤ 1, then We start with expanding the sums of products term in (A.2), given us

	of Proposition 2. t i=k t j=i+1 1 -ωγ j γ i = -	1 ω	t i=k	t j=i+1
			1 ω	t j=k	1 + ωγ j ≤	1 ω	t j=k	exp ωγ j =	1 ω	exp	        ω	t j=k	γ j	        .
	Proposition 2. t i=k	t j=i+1	1 -ωγ j γ i ≤	1 ω	.		(A.2)

  ω , using t j=k [1 -ωγ j ] ≥ 0, showing the inequality in (A.2). Proposition 3. Let (γ t ) t≥1 and (η t ) t≥1 be positive sequences. For any k ≤ t, we can obtain the (upper) bounds: Furthermore, suppose that for all i ≥ k, ωγ i ≤ 1, then Proposition 4. Let (δ t ) t≥0 , (γ t ) t≥1 , (η t ) t≥1 , and (ν t ) t≥1 be some positive sequences satisfying the recursive relation:δ t ≤ (1 -2ωγ t + η t γ t ) δ t-1 + ν t γ t , (A.5)

		t i=k	t j=i+1	1 + ωγ j η i γ i ≤	1 ω	max k≤i≤t	η i exp	        ω	t j=k	γ j	        ,	(A.3)
	with ω > 0. k≤i≤t	η i	t i=k	t j=i+1	1 + ωγ j γ i ≤	1 ω	max k≤i≤t	η i	t j=k	1 + ωγ j ≤	1 ω	max k≤i≤t	η i exp	        ω	t j=k	γ j	        .
	Similarly, for the inequality in (A.4), we have
	t i=k	t j=i+1	1 -ωγ j η i γ i ≤ max k≤i≤t	η i	t i=k	t j=i+1	1 -ωγ j γ i ≤	1 ω	max k≤i≤t	η i ,
	by Proposition 2.

t i=k t j=i+1 1 -ωγ j η i ≤ 1 ω max k≤i≤t η i . (A.4)

Proof of Proposition 3. We obtain the inequality in (A.3) directly by Proposition 1:

t i=k t j=i+1 1 + ωγ j η i γ i ≤ max

  -2ωγ j + η j γ j ν i γ i A t, where B t can be seen a transient term only depending on the initialisation δ 0 , and a stationary term A t . The transient term B t can be divided into two products, before and after t 0 ,  Next, the stationary term A t can (similarly) be divided into two sums (after and before t 0 ): -2ωγ j + η j γ j ν i γ i [1 -ωγ j ] ≤ exp(-ω t j=t/2 γ j ) as 1 + t ≤ exp(t) for all t ∈ R. Thus, for all t ∈ R, = 0 if t/2 < t 0 . The second stationary term A t,2 can be bounded, thanks to Proposition 1, as follows:

	1 -2ωγ i + η i γ i j=i+1 t 1 A t,1 1 -2ωγ t 0 +1≤i≤t δ 0 + t t     . i=t 0 +1 t + t 0 i=1 t j=i+1 ν i t i=t 0 +1 t j=i+1 1 -ωγ j γ i = 1 ω max t 0 +1≤i≤t ν i ≤ by Proposition 3. Furthermore, if t/2 > t 0 + 1, we get A t = A t,1 ≤ t i=t 0 +1 t j=i+1 1 -ωγ j ν i γ i = t/2-1 i=t 0 +1 t j=i+1 1 -ωγ j ν i γ i + t i=t/2 t j=i+1 1 -ωγ j ν i γ i 1 ω max t/2≤i≤t ν i , ≤ t/2-1 i=t 0 +1 t j=t/2 t/2≤i≤t ν i t i=t/2 t j=i+1 1 -ωγ j γ i = t j=t/2 1 -ωγ j t/2-1 i=t 0 +1 where t j=t/2 A t,1 ≤ exp         -ω t j=t/2 γ j         t/2-1 i=t 0 +1 ν i γ i + 1 ω max t/2≤i≤t ν i , 1  where t/2 t 0	ν i γ i +	1 ω	max t/2≤i≤t	ν i ,	(A.7) (A.8)
	i=1	j=i+1
	B t

j + η j γ j ν i γ i A t,2 .

The first stationary term A t,1 (with t > t 0 ) can be bounded as follows: if t/2 ≤ t 0 + 1, we have

A t,1 ≤ max 1 -ωγ j ν i γ i + max

Acknowledgement

This work was supported by grants from Région Île-de-France.

) and (γ t ) are positive sequences.