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

Nicklas Werge
email: nicklas.werge@upmc.fr

Olivier Wintenberger
email: olivier.wintenberger@upmc.fr

Learning from time-dependent streaming data with online stochastic algorithms

Keywords: stochastic optimization, machine learning, stochastic algorithms, online learning, streaming, time-dependent data

We study stochastic algorithms in a streaming framework, trained on samples coming from a dependent data source. In this streaming framework, we analyze the convergence of Stochastic Gradient (SG) methods in a non-asymptotic manner; this includes various SG methods such as the well-known stochastic gradient descent (i.e., Robbins-Monro algorithm), mini-batch SG methods, together with their averaged estimates (i.e., Polyak-Ruppert averaged). Our results form a heuristic by linking the level of dependency and convexity to the rest of the model parameters. This heuristic provides new insights into choosing the optimal learning rate, which can help increase the stability of SGbased methods; these investigations suggest large streaming batches with slow decaying learning rates for highly dependent data sources.

Introduction

Over the past decade, machine learning and artificial intelligence have become mainstream in many parts of society; substantial improvements in the performance and cost of mass storage devices and network systems have contributed to this. Traditional machine learning methods often work in a batch or offline learning setting, where the model is re-trained from scratch when new data arrive. Such learning methods suffer some critical drawbacks, such as expensive re-training costs when dealing with new data and thus poor scalability for large-scale and real-world applications. At the same time, these intelligent systems generate a practically infinite amount of large-scale data sets, many of which come as a continuous data stream, so-called streaming data.

Streaming data arrives as an endless sequence of samples (data points), which means that at any given time, the model must be able to adapt to the samples observed (so far) to predict/label new samples accurately. Such (streaming) models can never be seen as complete but must be updated continuously as newer samples arrive. Methods that recalculate the model from scratch on the arrival of new samples are impractical due to their high computational cost. Therefore we need procedures that effectively update the model as more samples arrive. This computational efficiency should not be at the expense of accuracy; the model's accuracy should be close to that achieved if we built a model from scratch using all the samples [START_REF] Bottou | Large scale online learning[END_REF].

Stochastic approximation algorithms have proven effective in overcoming the drawbacks of traditional (batch/offline) machine learning methods as they only use samples one by one without knowing their number in advance, especially the Stochastic Gradient (SG) method [START_REF] Robbins | A stochastic approximation method[END_REF]. These SG methods have proven scalable and robust in many areas ranging from smooth and strongly convex problems to complex non-convex ones, which makes them applicable in many large-scale machine learning tasks for real-world applications where data are large in size (and dimension) and arrive at a high velocity. Such first-order methods have been intensively studied in theory and practice in recent years [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF].

The classical analyses for SG methods typically require unbiased gradients drawn independently and identically distributed (i.i.d.) from some underlying (and unknown) data generation process [START_REF] Cesa-Bianchi | On the generalization ability of on-line learning algorithms[END_REF]. However, in practice, learning often happens with non-i.i.d. (and biased) data, e.g., network traffic, meteorological, financial time series, or other sensor data. We go beyond these standard assumptions by allowing dependent and biased gradients. SG methods can converge even when they only have access to biased gradients, but most analysis has been developed with specific applications in mind [START_REF] Ajalloeian | On the convergence of sgd with biased gradients[END_REF][START_REF] Bertsekas | Nonlinear Programming[END_REF]d'Aspremont, 2008;[START_REF] Devolder | Stochastic first order methods in smooth convex optimization[END_REF][START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF]. Stochastic learning algorithms for non-i.i.d. data are not as well understood as for i.i.d. data; however, some researchers have examined the convergence of statistical learning algorithm in non-i.i.d. settings [START_REF] Agarwal | The generalization ability of online algorithms for dependent data[END_REF][START_REF] Mohri | Stability bounds for stationary ϕ-mixing and β-mixing processes[END_REF][START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF].

Solving the problem of stochastic approximations using streaming SGs methods means we must approach the objective using the gradually arriving samples drawn according to some unknown dependent process. This leads to some new challenges, e.g., this endless stream of samples (may) changes at each step (and arrives sequentially), meaning that streaming SGs must be able to adapt to varying arrival speeds without compromising accuracy. We present and analyze streaming SGs that overcome these challenges and achieve convergence in various settings with long-and short-range dependence, model misspecification, and changing data streams.

Contributions. In this paper, we investigate SG methods in a streaming framework [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF], where the data comes from dependent stochastic processes. We provide non-asymptotic analysis and quantify the magnitude of achievable convergence rates under various dependency structures and convexity levels. Our framework covers many applications with dependence and biased gradients under weak gradient assumptions. Our results construct a heuristic between the level of dependency, noise, and convexity and the achievable learning rate to obtain optimal convergence. Generally, SG methods can achieve convergence using non-decreasing (streaming) batch sizes, which counteract the long-range (and short-range) dependence and model misspecification. We show that biased SG methods converge with the same accuracy as unbiased SG methods if the bias is not too large. More surprisingly, these heuristics can be used in practice to help increase the stability of SG-based methods.

Organization. Section 2 presents the streaming framework on which the non-asymptotic analysis relies; we introduce some key concepts, definitions, and assumptions. In particular, Section 2.2 contains the assumptions about dependency structures and gradients, with some examples of how these could be verified using mixing conditions. Our convergence results are presented in Section 3, with and without averaging (Sections 3.1 and 3.2). Each result is followed by a thorough discussion that relates to other work. All our convergence analysis depends on the assumptions in Section 2 and some additional conditions for the averaged case (Section 3.2). At last, experimentations of our findings are illustrated in Section 4, with some final remarks in Section 5.

Problem Formulation

We consider the Stochastic Optimization (SO) problem min θ∈Θ L(θ) = E t [l t (θ)], where Θ ⊆ R d is a convex body1 and l t : R d → R is some differentiable random functions (possibly non-convex), e.g, see Nesterov et al. (2018). We solve the SO problem in a streaming framework, where a block l t = (l t,1 , . . . , l t,n t ) of n t ∈ N random functions arrives at any given time t ∈ N. In solving the SO problem, we use the Stochastic Streaming Gradient (SSG) estimate proposed by [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF], given as

θ t = θ t-1 - γ t n t n t i=1 ∇ θ l t,i (θ t-1 ) , θ 0 ∈ Θ, (1)
where γ t is the learning rate satisfying the conditions [START_REF] Robbins | A stochastic approximation method[END_REF]. Note that if ∀t, n t = 1, SSG becomes the well-known SG method, which has attracted a lot of attention [START_REF] Bousquet | Stability and generalization[END_REF][START_REF] Hardt | Train faster, generalize better: Stability of stochastic gradient descent[END_REF][START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF][START_REF] Xiao | Dual averaging method for regularized stochastic learning and online optimization[END_REF][START_REF] Zhang | Solving large scale linear prediction problems using stochastic gradient descent algorithms[END_REF]. Almost surely convergence of SO algorithms were shown in [START_REF] Pelletier | On the almost sure asymptotic behaviour of stochastic algorithms[END_REF]. In many models, there may be constraints on the parameter space, which would require a projection of the parameters; therefore, we also introduce the Projected Stochastic Streaming Gradient (PSSG) estimate, defined by

∞ i=1 γ i = ∞ and ∞ i=1 γ 2 i < ∞ (
θ t = P Θ        θ t-1 - γ t n t n t i=1 ∇ θ l t,i (θ t-1 )        , θ 0 ∈ Θ, (2)
where P Θ denotes the Euclidean projection onto Θ, i.e., P Θ (θ) = arg min θ ∈Θ θθ 2 . To shorten notation, we let θ). An essential extension is the Polyak-Ruppert averaging [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], which guarantees optimal statistical efficiency without jeopardizing the computational cost; the Averaged Stochastic Streaming Gradient (ASSG) is given by θt

∇ θ l t (θ) = n -1 t n t i=1 ∇ θ l t,i (
= 1 N t t-1 i=0 n i+1 θ i , θ0 = 0, (3)
where N t = t i=1 n i is the accumulated sum of observations. Likewise, let PASSG denote the (Polyak-Ruppert) averaged estimate of PSSG (2).

Quasi-strong Convex Objectives

Following [START_REF] Gower | Sgd: General analysis and improved rates[END_REF]; [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], we assume that L has a unique global minimizer θ * ∈ Θ such that ∇ θ L(θ * ) = 0, and it is µ-quasi-strongly convex [START_REF] Karimi | Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition[END_REF][START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF], i.e, there exists µ > 0 such that ∀θ ∈ Θ,

L(θ * ) ≥ L(θ) + ∇ θ L(θ), θ -θ * + µ 2 θ -θ * 2 . ( 4
)
The µ-quasi-strongly convexity assumption is a non-strongly convex relaxation of the SO problem, which is more conservative than µ-strongly convexity. Relaxations of convexity is crucial in practice to ensure robustness and adaptiveness of the algorithms, e.g., for non-strongly convex SO, see [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n)[END_REF]; [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF]; [START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF].

Stochastic

Streaming Gradient Assumptions: Dependence, Biased, Expected Smoothness, and Gradient Noise We go beyond the classical assumptions that require unbiased (uniformly bounded) gradients by allowing the gradients to be dependent and biased estimates. Our aim is to non-asymptotically bound the SSG estimates (1) to (3) explicitly using the SO problem parameters. In order to do this, we let the natural filtration of the SO problem F t = σ(l i : i ≤ t), and assume the following about the gradients (∇ θ l t ):

Assumption 1-p (D ν ν t -dependence and B ν ν t -bias). Let θ 0 be F 0 -measurable. For each t ≥ 1, the random function ∇ θ l t (θ) is square-integrable, F t -measurable, and there exists a positive integer p such that for all

F t-1 -measurable θ ∈ Θ, E[ E[∇ θ l t (θ)|F t-1 ] -∇ θ L(θ) p ] ≤ ν p t (D p ν E[ θ -θ * p ] + B p ν ), (5)
for some positive sequence (ν t ) t≥1 with D ν , B ν ≥ 0.

In the classical convergence analysis of SG methods, one assumes that the SGs are uniformly bounded. However, this assumption is too restrictive as it only may hold for some losses [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF][START_REF] Nguyen | Sgd and hogwild! convergence without the bounded gradients assumption[END_REF]. Instead, we follow the same ideas as in [START_REF] Gower | Sgd: General analysis and improved rates[END_REF]; [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], to make the following assumption about the expected smoothness of the stochastic gradients (∇ θ l t ).

Assumption 2-p (κ t -expected smoothness). There exists a positive integer p such that ∀θ,

θ ∈ Θ, E[ ∇ θ l t (θ) - ∇ θ l t (θ ) p ] ≤ κ p t E[ θ -θ p ] for some positive sequence (κ t ) t≥1 .
Assumption 2-p can be seen as an assumption about the smoothness properties of (l t ). The last fundamental assumption (Assumption 3-p) is a very weak assumption, and should be seen as an assumption on Θ rather than on (l t ): Assumption 3-p (σ t -gradient noise). There exists a positive integer p such that E[ ∇ θ l t (θ * ) p ] ≤ σ p t for some positive sequence (σ t ) t≥1 .

These assumptions (Assumptions 1-p to 3-p) are milder than the standard assumptions for stochastic approximations, e.g., see [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF]; [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF]; [START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF]; [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]. They include classic examples such as stochastic approximation and learning from dependent data, which we will demonstrate later in Section 4. Assumption 1-p is on the form of mixing conditions for weakly dependence sequences, implying that dependence dilutes with the rate of ν t . It is possible to verify Assumption 1-p by using moment inequalities for partial sums of strongly mixing sequences [START_REF] Rio | Asymptotic theory of weakly dependent random processes[END_REF]; we will refer to this as short-range dependence. Note that for any positive integer p, Assumption 1-p can be upper bounded by

E[ E[∇ θ l t (θ)|F t-1 ] -∇ θ L(θ) p ] ≤ E[ ∇ θ l t (θ) -∇ θ L(θ) p ] = n -p t E[ S t p ], (6)
using Jensen's inequality, where

S t = n t i=1 (∇ θ l t,i (θ) -∇ θ L(θ)) is a d-dimensional vector. Let (∇ θ l t,i
) be a strictly stationary sequence and assume that there exists some r > p such that sup x>0 (x r Q(x)) 1/r < ∞, where Q(x) denotes the quantile function of ∇ θ l t,i . Suppose that (∇ θ l t,i ) is strongly α-mixing in the sense of [START_REF] Rosenblatt | A central limit theorem and a strong mixing condition[END_REF], with strong mixing coefficients (α t ) t≥1 satisfying α t = O(t -pr/(2r-2p) ). Then by Rio (2017, Corollary 6.1) [START_REF] Bradley | Basic properties of strong mixing conditions. a survey and some open questions[END_REF]; [START_REF] Doukhan | Mixing: properties and examples[END_REF] for more examples, other mixing coefficients of weak dependence and the relations between them. In relation to the form of Assumption 1-p, this means that B ν 0 in this case. However, having B ν = 0 is possible in well-specified examples, which we will see later in Section 4. Note that Assumptions 2-p and 3-p can be verified using α-mixing conditions by analogues arguments as for Assumption 1-p such that κ p t and σ p t is O(n -p/2 t ).

Convergence Analysis

In this section, we consider the stochastic streaming estimates in (1) to (3) with streaming-batches (n t ) arriving in non-decreasing streams. We aim to non-asymptotically bound δ t = E[ θ t -θ * 2 ] and δt = E[ θtθ * 2 ], such that they only depend on the parameters of the problem.

Learning rate and function forms. Throughout this paper, we consider learning rates on the form γ t = C γ n β t t -α with C γ > 0, β ∈ [0, 1], and α chosen accordingly to the expected streaming-batches n t . Obviously, (ν t ), (κ t ), and (σ t ) may be considered as uncertain terms depending on the streaming-batch n t . Thus, let ν t = n -ν t , κ t = C κ n -κ t , and σ t = C σ n -σ t with ν ∈ (0, ∞), κ, σ ∈ [0, 1/2], and C κ , C σ > 0. Having, σ, κ ∈ [0, 1/2] follows directly from [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF], since σ = κ = 1/2 corresponds to the i.i.d. case 2 , whereas σ, κ < 1/2 allows noisier outputs. Similarly, v t = 0 corresponds to the classical i.i.d. setting. Having ν t = n -ν t means Assumption 1-p, allow so-called long-range dependence (also known as long memory or long-range persistence) when ν ∈ (0, 1/2) and short-range dependence when ν ∈ [1/2, ∞). Thus, the i.i.d. case is when ν → ∞.

For the sake of simplicity, we consider streaming-batches (n t ) on the form C ρ t ρ with C ρ ∈ N and ρ ∈ [0, 1) such that n t ∈ N. This form of streaming-batches means that we are considering everything from vanilla SG and mini-batch SG methods, to more exotic learning designs, e.g., C ρ > 1 and ρ = 0 correspond to mini-batch SG of size C ρ . We will refer to C ρ as the streaming constant size and ρ as the streaming rate.

Stochastic Streaming Gradients

Theorem 1. Denote δ t = E[ θ t -θ * 2 ]
for some δ 0 ≥ 0, where (θ t ) follows the recursion in (1) or (2). Assume that Assumptions 1-p to 3-p hold true for p = 2. Suppose n

t = C ρ t ρ with ρ ∈ [0, 1) and C ρ ∈ N, such that µ ν = µ -1 {ρ=0} 2D ν C -ν ρ > 0. For α -ρβ ∈ (1/2, 1), we have δ t ≤ π t + 2 2+6ρν 1+ρ B 2 ν µµ ν C 2ν 1+ρ ρ N 2ρν 1+ρ t + 2 7+6ρσ 1+ρ C 2 σ C γ µ ν C 2σ-β-α 1+ρ ρ N ρ(2σ-β)+α 1+ρ t , (7)
with π t given in (22) such that π t = O(exp(-N (1+ρβ-α)/(1+ρ) t

)).

Sketch of proof.

Under Assumptions 1-p to 3-p with p = 2, it can be shown that (δ t ) satisfies the recursive relation (20),

δ t ≤ [1 -(µ -2D ν ν t )γ t + 2κ 2 t γ 2 t ]δ t-1 + B 2 ν µ ν 2 t γ t + 2σ 2 t γ 2 t ,
for any γ t , ν t , κ t , σ t , and n t . This recursive relation can be explicitly upper bounded in a non-asymptotic way (by Proposition 1) using classical techniques from stochastic approximations [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF][START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF]. As mentioned in [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], bounding the projected estimate in (2) follows directly from that

E[ P Θ (θ)- θ * 2 ] ≤ E[ θ -θ * 2 ], ∀θ ∈ R d , ∀θ * ∈ Θ, as Θ is a convex body.
Related work. Theorem 1 replicate the results of the well-specified i.i.d. case (with B ν = 0 and κ = σ = 1/2) considered in [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF]. Our findings also reproduce the results of [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], where they considered the well-specified i.i.d. case (under slightly different assumptions) using the vanilla SG method, namely, when C ρ = 1 and ρ = 0. Moreover, if the function L has C ∇ -Lipschitz continuous gradients3 , then (7) implies the bound on the objective function values of L, E[L(θ t ) -L(θ * )] ≤ C ∇ δ t /2 by Cauchy-Schwarz's inequality.

Decay of the initial conditions. The initial conditions that π t contains will be forgotten sub-exponentially fast, since π t = O(exp(-N (1+ρβ-α)/(1+ρ) t )) as long as µ ν = µ -1 {ρ=0} 2D ν C -ν ρ > 0. Note that the positivity of the dependence penalised convexity constant µ ν is essential in all terms of (7). Having µ ν > 0 depends solely on the level of dependence D ν but it is scaled by C -ν ρ , meaning if D ν is so large that µ ν is no longer positive, then we should take C ρ large enough such that µ ν becomes positive again; this is illustrated in Sections 4.2 and 4.3. The streaming constant C ρ contributes positively to all terms in (7), either directly or though µ ν .

The last term of ( 7) can be seen as the noise term decaying with O(N -(ρ(2σ-β)+α)/(1+ρ) t

) for αρβ ∈ (1/2, 1), e.g., for any ρ ∈ [0, 1), δ t = O(N -2/3 t ) when α = 2/3, β = 1/3, and σ = 1/2. In addition, the noise term is positively affected by large streaming constants C ρ when α + β < 2σ, which will be expressed as a variance reduction, e.g., see Section 4. In well-specified cases (B ν = 0) the noise term would also be the asymptotic term.

Behavior for B ν . The second term of (7) can be seen as an dependency term as it is determined solely by the level of dependence ν, the bias (misspecification error) B ν , and the convexity constant µ ν ; It is remarkable that the dependence term is unconnected from the choice of the learning rate (γ t ) but instead by the streaming rate through C ρ and ρ. The dependence term decay with O(N -2ρν/(1+ρ) t ) which requires ρ positive to decay since ν ∈ (0, ∞), e.g., to obtain O(N -1/2 t ) we would need ρ = 1 and ν = 1/2. It is surprising that Theorem 1 allows both long-range and short-range dependence. Indeed, long-range dependence leads to slow convergence (slower than O(N -1/2 t )) but it will still converge. Obviously, this only matters if

B ν 0. Overall, δ t = O(max{1 {B ν 0} N -2ρν/(1+ρ) t , N -(ρ(2σ-β)+α)/(1+ρ) t }).

Averaged Stochastic Streaming Gradients

In what follows, we consider the averaging estimate θn given in (3) with (θ t ) following the SSG estimate in (1) or the PSSG estimate in (2). Some additional assumptions is needed for bounding the rest terms of the averaging estimate: let the function L have C ∇ -Lipschitz continuous gradients, i.e., there exists a constant

C ∇ > 0, ∀θ, θ ∈ Θ ⊆ R d , ∇ θ L(θ) -∇ θ L(θ ) ≤ C ∇ θ -θ . ( 8
)
As discussed in [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF], this assumption ensures that ∇ θ L does not vary arbitrarily, making the gradient ∇ θ L a useful indicator on how to decrease L. Next, assume that the Hessian of

L is C ∇ -Lipschitz-continuous, that is, there exists C ∇ > 0 such that ∀θ, θ ∈ Θ ⊆ R d , ∇ 2 θ L(θ) -∇ 2 θ L(θ ) ≤ C ∇ θ -θ . (9)
Note that ( 8) and ( 9) only needs to hold true for θ = θ * . Moreover, in continuation of Assumption 3-p with

σ t = C σ n -σ t for σ ∈ [0, 1/2],
we make the following assumption:

Assumption 4. There exists a non-negative self-adjoint operator Σ such that ∀t ≥ 1,

n 2σ t E[∇ θ l t (θ * )∇ θ l t (θ * ) ] Σ + Σ t , where Σ t is a positive symmetric matrix with Tr(Σ t ) = C σ n -2σ t , C σ ≥ 0, and σ ∈ (0, 1/2].
Remark that in the independent or some well specified cases such as in Section 4.1.1, Assumption 4 is verified with σ = 1/2 and C σ = 0 [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF]. The short-range dependence cases is when σ = 1/2, as in Section 4.1.2, whereas, the long-range dependence case is for σ < 1/2. Moreover, Assumption 4 allows us to obtain leading term Λ/N t with Λ = Tr(∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ), which attains the Cramer-Rao bound; we will see this in Theorem 2.

To consider the averaging estimate θn given in (3) but with the use of the projected estimate PSSG from (2), which we will denote PASSG. An additional assumption is needed in order to avoid calculating the six-order moment, we make the unnecessary assumption that (∇ θ l t ) is uniformly bounded; the derivation of the six-order moment can be found in [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 (θ) 2 ≤ G 2 Θ a.s. Theorem 2. Denote δt = E[ θt -θ * 2 ]
with θn given by (3), where (θ t ) follows the recursion in (1) or (2). Assume that Assumptions 1-p to 3-p for p = 4 and Assumption 4 hold true. In addition, Assumption 5 must hold true if

(θ t ) follows the recursion in (2). Suppose n t = C ρ t ρ with ρ ∈ [0, 1) and C ρ ∈ N, such that µ ν = µ -1 {ρ=0} 2D ν C -ν ρ > 0. For α -ρβ ∈ (1/2, 1), we have δ1/2 t ≤ Λ 1/2 N 1/2 t 1 {σ=1/2} + 2 1/2 Λ 1/2 C 1-2σ 2(1+ρ) ρ N 1+2ρσ 2(1+ρ) t 1 {σ 1/2} + 2 1/2 C 1/2 σ C 1-2(σ+σ ) 2(1+ρ) ρ µN 1+2ρ(σ+σ ) 2(1+ρ) t (10)
+ O max N -2+ρ(2σ+β)-α 2(1+ρ) t , N -ρ(2σ-β)+α 1+ρ t + Õ N -δ+ρν 2(1+ρ) t + 1 {B ν 0} Ψ t , (11)
with δ = 1 {B ν =0} (ρ(2σ -β) + α) + 1 {B ν 0} min{ρ(2σ -β) + α, 2ρν}
and Ψ t given in (36), such that

Ψ t = Õ max N -ρ(σ+ν) 2(1+ρ) t , N -1+ρ(β+ν)-α 1+ρ t , N -1+2ρν 2(1+ρ) t , N -δ/2+ρν 2(1+ρ) t , N -2ρν 1+ρ t .
An explicit version of the bound is given in (37).

Sketch of proof.

In Lemma 3, we conduct a general study of the Polyak-Ruppert averaging estimate ( θt ) defined in (3) for (γ t ), (ν t ), (κ t ), (σ t ) and (n t ) on any form. Thus, Theorem 2 follows by Lemma 3 using the (specific) bounds of

δ t = E[ θ t -θ * 2 ] and ∆ t = E[ θ t -θ * 4 ]
in Theorem 1 (eq. ( 21)) and Lemma 2.

Related work. Similarly to the well-specified i.i.d. case [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF], the leading term of (10) is Λ/N t , which obtain the (asymptotically optimal) Cramer-Rao bound [START_REF] Murata | Statistical analysis of learning dynamics[END_REF]. Each term of ( 10) is a direct consequence of Assumption 4 and they are all independent of the choice of learning rate (γ t ). Moreover, as discussed in [START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF], the bound of δt can be seen as a bias-variance decomposition between the leading terms (10) and the remaining terms in (11).

Accelerated decay. By averaging it is possible (in some specific cases) to achieve the desirable Cramer-Rao bound, namely, the leading term Λ/N t could obtain the optimal and incorrigible rate of O(N -1 t ). This is always achieved in the well-specified case with σ = 1/2, even under short-range dependence (i.e., when ν > 1/2).

As for Theorem 1, the positivity of µ ν is essential for all terms in (11) even if it does not appear directly. In case of lack of convexity µ or high levels of dependence constant D ν , we can only ensure convergence by increasing C ρ , i.e., ensuring positivity of µ ν ; this is illustrated in Sections 4.2 and 4.3 for ARCH models.

The first term of (11) decay at the rate

O(N -(2+ρ(β+2σ)-α)/(1+ρ) t ) or O(N -2(ρ(2σ-β)+α)/(1+ρ) t
), which suggests choosing α, β such that α + ρ(2σ/3 -β) = 2/3, e.g., α = 2/3, β = 1/3 and σ = 1/2 yields a decay of O(N -4/3 t ) for any ρ. Thus, we can robustly achieve O(N -4/3 t ) for any streaming rate ρ by setting α = 2/3 and β = 1/3 if σ = 1/2. In general, the convergence is resilient to any streaming rate ρ by having α = 2/3 and β = 2σ/3. But taking β > 0 would damage the variance reduction effect from having C ρ large (e.g., see discussion after Theorem 1). Thus, there is a trade-off between accelerating the convergence by taking β = 2σ/3 > 0 or taking β = 0 to favor from variance reduction. In practice, an immediate choice would be to take β = 0, but if the data or model contains a low amount of noise, it can be advantageous to raise β to improve convergence [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF].

Next, the decay of the second term in ( 11) is tricky to interpret in a simple manner as it is a mixture of the learning rate, streaming rate, dependence, and bias (misspecification error). Nevertheless, some observations can be made: first, having β = 0 is beneficial for the decay rate δ in all cases. Similarly, increasing streaming rate ρ would also increase the decay. The most important thing to mention is that if B ν 0 then we would at most have a decay of O(N -2/3 t ), which is the same as for the SSG in Theorem 1.

Behavior for B ν . The influence of B ν is exclusively contained in Ψ t , with the exception of the second term of (11). Also, increasing ρ will always diminish the bad influence of this bias term. Surprisingly, Ψ t → 0 as t → ∞ for any ν, but long-range dependence is excluded if we wish to obtain the desired rate of δt = O(N -1 ). However, it does not seem to have any major influence in our experiments in Section 4.

Experiments

A way to illustrate our findings is by use of classical methods that aims to construct models for time-series analysis, modeling, and prediction of the underlying sequences of real-valued signals (X t ). These methods have been successfully used in a wide range of applications such as statistics, econometrics, and signal processing because of their ability to describe or predict time-varying (dependent) processes, e.g., the AutoRegressive (AR), Moving-Average (MA), and AutoRegressive Moving-Average (ARMA) models are the most well-known models for timeseries [START_REF] Box | Time series analysis: forecasting and control[END_REF][START_REF] Brockwell | Time series: theory and methods[END_REF][START_REF] Hamilton | Time series analysis[END_REF]. The standard time-series analysis often relies on independence and constant noise, but it can be relaxed by, e.g., the AutoRegressive Conditional Heteroskedasticity (ARCH) model [START_REF] Engle | Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation[END_REF]. Online learning algorithms of (both stationary and non-stationary) dependent timeseries have been studied in [START_REF] Agarwal | The generalization ability of online algorithms for dependent data[END_REF]; [START_REF] Anava | Online learning for time series prediction[END_REF]; [START_REF] Wintenberger | Stochastic online convex optimization; application to probabilistic time series forecasting[END_REF].

Our experiment measures the performance by the quadratic mean error E[ θ N t -θ * 2 ] over one thousand replications with θ 0 and θ * drawn randomly according to the models' specifications. It should be noted that averaging over several replications gives a reduction in variability, that mainly benefits the SSG. The experiments will demonstrate how the choice of C ρ and ρ affects the dependence D ν , bias B ν , and the (dependence) penalised convexity constant µ ν . To compare different data streams n t = C ρ t ρ through the selection of C ρ and ρ, we fix the following parameters: C γ = 1, α = 2/3, and β = 0.

AR model

A process (X t ) is called a (zero-mean) AR(1) process, if there exists real-valued parameter θ such that X t = θX t-1 + t , where ( t ) is some noise process with zero mean and noise σ . To illustrate the versatility of our results, we construct some noisy (heavy-tailed) data with long-range dependence: the noisiness is integrated using a Student's t-distribution with degrees of freedom above four, denoted by (z t ). The long-range dependence is incorporated by multiplying (z t ) with the fractional Gaussian noise G t (H) = B t+1 (H) -B t (H), where (B t (H)) is a fractional Brownian motion with Hurst index H ∈ (0, 1). (B H t ) can also be seen as a (zero-mean) Gaussian process with stationary and self-similar increments. Thus, let the AR(1) process X t be constructed using the noise process

t = √ G t (3/4)z t , where a Hurst index of H = 3/4 corresponds to ν 2 t , κ 2 t , σ 2 t is O(n -1/2 t ) and ν 4 t , κ 4 t , σ 4 t is O(n -3/4 t
) in Assumptions 1-p to 4.

Well-specified case

Consider the well-specified case, in which, we estimate an AR(1) model X t = θX t-1 + t from the underlying stationary AR(1) process X t = θ * X t-1 + t with |θ * | < 1. We omit to project our estimates as this will hide the dependence coming from D ν , which is what we wish to explore. For constant streaming-batch sizes of one, the squared loss is l t (θ) = (X t -θX t-1 ) 2 with ∇ θ l t (θ) = -2X t-1 (X t -θX t-1 ). This gives a mean squared loss

L(θ) = E t [l t (θ)] = E[(X t -θX t-1 ) 2 ] = E[(θ * X t-1 + t -θX t-1 ) 2 ] = (θ * -θ) 2 E[X 2 t-1 ] + σ 2 , with ∇ θ L(θ) = -2(θ * -θ)E[X 2 t-1 ]. Thus, Assumption 1-p (for p = 2 with σ(X t-1 ) ⊆ F t-1 ) yields E[ E[∇ θ l t (θ)|F t-1 ] -∇ θ L(θ) 2 ] =E[(E[2X t-1 (θX t-1 -X t )|F t-1 ] -2(θ -θ * )E[X 2 t-1 ]) 2 ] = 4(θ -θ * ) 2 E[(X 2 t-1 -E[X 2 t-1 ]) 2 ],
meaning that Assumption 1-p is fulfilled if X t has bounded moments of order p. Moreover, from this we can directly deduce that B ν = 0. Likewise, the remaining assumptions can be verified, in particular Assumption 4 is satisfied with Σ t = 0.

Misspecified case

Next, assume that the underlying data generating process follows the MA(1)-process, X t = φ t-1 + t , with φ ∈ R. The misspecification error of fitting an AR(1) model to a MA(1) process can be found by minimizing L(θ),

θ * = arg min θ E[(X t -θX t-1 ) 2 ] = arg min θ E[( t + φ t-1 -θ( t-1 + φ t-2 )) 2 ] = arg min θ E[( t + (φ -θ) t-1 -θφ t-2 ) 2 ] = arg min θ σ 2 + (φ -θ) 2 σ 2 + θ 2 φ 2 σ 2 ≡ arg min θ (φ -θ) 2 + θ 2 φ 2 = arg min θ L(θ), where L(θ) = (φ -θ) 2 + θ 2 φ 2 is a strictly convex function in θ. Thus, ∇ θ L(θ) = 0 ⇔ 2(θ -φ) + 2θφ 2 = 0 ⇔ 2θ(1 + φ 2 ) = 2φ ⇔ θ = φ/(1 + φ 2 ), meaning for φ ∈ R we have θ ∈ (-1/2, 1/2).
With this in mind, we can conduct our study of fitting an AR(1) model to the MA(1) process with φ drawn randomly from R.

ARCH model

A key element of time series analysis is modeling heteroscedasticity of the conditional variance, e.g., volatility clustering in financial time-series; ARCH models are some of the most well-known models that incorporate this feature. A process ( t ) is called an ARCH(1) process with parameters α 0 and α 1 if it satisfies

       t = σ t z t , σ 2 t = α 0 + α 1 2 t-1 , (12)
where α 0 > 0 and α 1 ≥ 0 ensures the non-negativity of the conditional variance process (σ 2 t ), and the innovations (z t ) is white noise. The ARCH process parameters are known to be challenging to estimate in empirical applications as the optimization algorithms can quickly fail or converge to irregular solutions. Therefore, projecting the estimates is vital for the optimization procedure. A well-discussed problem for the ARCH models is that small values of α 0 are tricky to estimate. Stabilizing the estimation of α 0 would not only improve the α 0 estimate but also have a positive impact on the other model parameters. One way to deal with small values of α 0 is by the using the models homogeneity, i.e., scaling an ARCH process (X t ) with parameters (α 0 , α 1 ) gives us an ARCH process ( √ λX t ) with parameters (λα 0 , α 1 ) with same innovations. To simplify our analysis we consider a stationary ARCH(1) model, where we fix α 0 at 1 and initialize it at 1/2. We employ the quasi-maximum likelihood procedure for the statistical inference as outlined in [START_REF] Werge | Adavol: An adaptive recursive volatility prediction method[END_REF]; the quasi likelihood losses is given by l

t (θ) = 2 -1 (X 2 t /σ 2 t (θ) + log(σ 2 t (θ)) with first-order derivative ∇ θ l t (θ) = ∇ θ σ 2 t (θ) σ 2 t (θ) -X 2 t 2σ 4 t (θ) where ∇ θ σ 2 t (θ) = (1, X 2 t-1 ) T .
Observe that the loss function (l t ) itself is not strongly convex but only the objective function L may be strongly convex; convexity conditions of ARCH was investigated in [START_REF] Wintenberger | Stochastic online convex optimization; application to probabilistic time series forecasting[END_REF]. There are different ways to overcome lack of convexity: first, projecting the estimates such that the (conditional) variance process (σ 2 t ) stays away from zero (and close to the unconditional variance). Second, in the specific example of ARCH model, one could also recover convexity by implementing variance targeting techniques; an example using Generalized ARCH (GARCH) models can be found in [START_REF] Werge | Adavol: An adaptive recursive volatility prediction method[END_REF].

AutoRegressive (AR)-AutoRegressive Conditional Heteroskedasticity (ARCH) Model

We complete our experiments by considering an AR models with ARCH noise: the process (X t ) is called an AR(1)-ARCH(1) process with parameters θ, α 0 and α 1 if it satisfies

             X t = θX t-1 + t , t = σ t z t , σ 2 t = α 0 + α 1 2 t-1 . ( 13
)
where the innovations (z t ) is white noise. The statistical inference of this model is done using the squared loss for the AR-part and the QMLE for the ARCH part, e.g., see Sections 4.1.1 and 4.2.

Discussion of experiments

The experiments described earlier in Sections 4.1 to 4.3 can be found in Figure 1; here {C ρ = 1, ρ = 0} corresponds to the classical SG method and its (Polyak-Ruppert) average estimate, {C ρ = 64, ρ = 0} is a mini-batch SSG/ASSG, and {C ρ = 64, ρ = 1/2} is an increasing SSG/ASSG with initial batch size of C ρ = 64. First consider the AR illustration in Figures 1a and1b: each pair of data streams converges, but it is clear that the traditional SG method experiences a large amount of noise initially, particularly affecting the average estimate period but not its decay rate. 4 Both methods show a noticeable reduction in variance when C ρ increases, which is particularly beneficial in the beginning. Nevertheless, too large streaming batch sizes C ρ may hinder the convergence as this leads to too few iterations. Moreover, Figures 1a and1b indicates improving decay for SSG when increasing the streaming rate ρ. Conversely, ASSG does not see improvements in the same way, as we do not exploit the potential of using multiple observations through the β parameter, which could accelerate convergence, e.g., see [START_REF] Godichon-Baggioni | Non-asymptotic analysis of stochastic approximation algorithms for streaming data[END_REF] for a discussion in the (well-specified) i.i.d. case. It is surprising that we do not see any effect from Σ t in Assumption 4, but this seems to be an artifact effect in the proof as we need fourth-order moments.

In Figures 1c and1d, we have the experiments for the stationary ARCH(1) models, with and without the ARpart, respectively, as outlined in Sections 4.2 and 4.3. These figures illustrate the lack of convexity when using small streaming batch sizes, e.g., the traditional SG method and its average estimate, {C ρ = 1, ρ = 0} diverges. Remark that the lack of convexity is expressed through the lack positively of µ ν , which only larger streaming batch sizes C ρ can counteract. Moreover, the traditional SG method, {C ρ = 1, ρ = 0} is omitted in Figure 1d due to lack of convexity. Figure 1d shows that large (C ρ = 64) and non-decreasing (ρ ≥ 0) streaming batches can converge under difficult settings.

Conclusions

We studied the SO problem in a streaming framework using dependent and biased (gradient) estimates. In particular, we explored convergence rates of the SSG and ASSG algorithms in a non-asymptotic manner. The theoretical results formed heuristics that links the level of dependency and convexity to the rest of the model parameters. These heuristics provided new insights into determining optimal learning rates, which can help increase the stability of SG-based methods. Our experimentation verified these investigations suggesting large streaming batches with slow decaying learning rates for highly dependent data sources. Moreover, in large-scale learning problems with dependence, noisy variables, and lack of convexity, we know how to achieve (and accelerate) convergence and reduce noise through the learning rate and the treatment pattern of the data.

There are several ways to expand our work: first, we can extend our analysis to include streaming batches of any size (not in terms of streaming batch size and streaming rates). Second, an extension to non-strongly convex goals could be beneficial as it will provide more insight into how we can choose robust learning rates [START_REF] Bach | Non-strongly-convex smooth stochastic approximation with convergence rate o (1/n)[END_REF][START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF][START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]. At the same time, this learning rate could be made adaptive such that it is robust to poor initialization and requires less fine-tuning. This last objective is the most important for practitioners as it builds a universality across applications.

Proofs

Let us start by giving a short sketch of how our proofs section is structured: we start by deriving recursive relations to the desired quantities. Next, we derive a general bounds to the recursive relationship for any (γ t ), (ν t ), (κ t ), (σ t ), and (n t ). Finally, we insert the specific functions forms of (γ t ), (ν t ), (κ t ), (σ t ), and (n t ), which yield the results seen in Theorems 1 and 2. Before doing the proofs, we recall a repeating argument used to non-asymptotically bound recursive relations of form ( 14):

Proposition 1 (Godichon-Baggioni et al. ( 2021)). Suppose (ω t ), (α t ), (η t ), and (β t ) to be some non-negative sequences satisfying the recursive relation,

ω t ≤ [1 -2λα t + η t α t ]ω t-1 + β t α t , (14)
with ω 0 ≥ 0 and λ > 0. Let C ω ≥ 1 be such that λα t ≤ 1 for all t ≥ t ω with t ω = inf{t ≥ 1 : C ω η t ≤ λ}. Then, for (α t ) and (η t ) decreasing, we have the upper bound on (ω t ) given by

ω t ≤ τ t + 1 λ max t/2≤i≤t β i , with τ t = exp         -λ t i=t/2 α i                 exp        C ω t i=1 η i α i        ω 0 + 1 λ max 1≤i≤t β i + t/2-1 i=1 β i α i         . (15)
Proposition 1 shows a simple way to bound (ω t ) in ( 14); the bound in (15) consists of a sub-exponential term τ t and a noise term λ -1 max t/2≤i≤t β i . Thus, our attention is on reducing the noise term without damaging the natural decay of the sub-exponential term where τ t → 0 exponentially fast as t → ∞.

Later in the proofs, we will insert some specific types of the sequences above, resulting in different generalized harmonic numbers, which can be bounded with the integral test for convergence. Moreover, to present our results in terms of N t = t i=1 n i , we will use that (N t /2C ρ ) 1/(1+ρ) ≤ t ≤ (2N t /C ρ ) 1/(1+ρ) . To ease notation, we will make use of the functions ψ x (t), ψ y x (t) : R + \ {0} → R, given as

ψ x (t) =              t 1-x /(1 -x) if x < 1, 1 + log(t) if x = 1, x/(x -1) if x > 1,
and

ψ y x (t) =              t (1-x)/(1+y) /(1 -x) if x < 1, 1 + log(t 1/(1+y) ) if x = 1, x/(x -1) if x > 1, (16)
with y ∈ R + such that ψ y x (t) = ψ x (t 1/(1+y) ). Thus, t i=1 i -x ≤ ψ x (t) for any x ≥ 0. Furthermore, with this notation, we have that ψ y x (t)/t = O(t -(x+y)/(1+y) ) if x < 1, ψ y x (t)/t = O(log(t)t -1 ) if x = 1, and ψ y x (t)/t = O(t -1 ) if x > 1. Hence, for any x 0 , x 1 , x 2 , y ≥ 0, ψ y x 0 (t)/t = Õ(t -(x 0 +y)/(1+y) ) and ψ y x 1 (t)ψ y x 2 (t)/t = Õ(t -(x 1 +x 2 +y-1)/(1+y) ), where the Õ(•) notation suppress logarithmic factors.

Proofs for Section 3.1

In the following lemma, we derive an explicit recursive relation of δ t = E[ θ t -θ * 2 ] to non-asymptotically bound the t-th estimate of (1) for any (γ t ), (ν t ), (κ t ), (σ t ), and (n t ) using classical techniques from stochastic approximations [START_REF] Benveniste | Adaptive algorithms and stochastic approximations[END_REF][START_REF] Kushner | Stochastic Approximation and Recursive Algorithms and Applications[END_REF]. As mentioned in [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], bounding the projected estimate in (2) follows directly from that

E[ P Θ (θ) -θ * 2 ] ≤ E[ θ -θ * 2 ], ∀θ ∈ R d , ∀θ * ∈ Θ, as Θ is a convex body.
Lemma 1 (Second-order moment). Assume that Assumptions 1-p to 3-p hold true for p = 2. Suppose that

µ ν = µ -1 {ν t =C} 2D ν ν t > 0. Let 1 {ν t =C} and 1 {ν t ¬C} indicate whether (ν t ) is constant or not. Denote δ t = E[ θ t -θ * 2 ]
for some δ 0 ≥ 0, where (θ t ) follows the recursion in (1) or (2). For any learning rate (γ t ), we have

δ t ≤ π t + 2B 2 ν µµ ν max t/2≤i≤t ν 2 i + 4 µ ν max t/2≤i≤t σ 2 i γ i , with π t = exp         - µ ν 2 t i=t/2 γ i                exp        1 {ν t ¬C} 2C δ D ν t i=1 ν i γ i        exp        2C δ t i=1 κ 2 i γ 2 i        δ 0 + 2B 2 ν µµ ν max 1≤i≤t ν 2 i + 4 µ ν max 1≤i≤t σ 2 i γ i + B 2 ν µ t/2-1 i=1 ν 2 i γ i + 2 t/2-1 i=1 σ 2 i γ 2 i         .
Proof of Lemma 1. By taking the quadratic norm on (1), expanding the norm, and taking the expectation, we can derive the equation,

δ t = δ t-1 + γ 2 t E[ ∇ θ l t (θ t-1 ) 2 ] -2γ t E[ ∇ θ l t (θ t-1 ), θ t-1 -θ * ], (17)
where

δ t = E[ θ t -θ * 2 ] with δ 0 ≥ 0.
To bound the second term in (17), we use Assumptions 2-p and 3-p for p = 2, to obtain that

E[ ∇ θ l t (θ t-1 ) 2 ] =E[ ∇ θ l t (θ t-1 ) -∇ θ l t (θ * ) + ∇ θ l t (θ * ) 2 ] ≤2E[ ∇ θ l t (θ t-1 ) -∇ θ l t (θ * ) 2 ] + 2E[ ∇ θ l t (θ * ) 2 ] ≤ 2κ 2 t δ t-1 + 2σ 2 t , (18)
as x+y p ≤ 2 p-1 ( x p + y p ). As noted in [START_REF] Bottou | Optimization methods for large-scale machine learning[END_REF]; Nesterov et al. (2018), (4) implies that

∇ θ L(θ), θ-θ * ≥ µ θ -θ * 2 for all θ ∈ Θ ⊆ R d .
Next, since L is µ-strongly convex (4) and θ t-1 is F t-1 -measurable (Assumption 1-p), we can bound the third term on the right-hand side of ( 17) by

E[ ∇ θ l t (θ t-1 ), θ t-1 -θ * ] =E[ ∇ θ L(θ t-1 ), θ t-1 -θ * ] + E[ E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ), θ t-1 -θ * ] ≥µδ t-1 -D ν ν t δ t-1 -B ν ν t δ 1/2 t-1 , (19)
since

E[ E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ), θ t-1 -θ * ] ≥ -E[ E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ) θ t-1 -θ * ] ≥ -E[ E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ) 2 ] E[ θ t-1 -θ * 2 ] ≥ -ν 2 t (D 2 ν δ t-1 + B 2 ν ) δ t-1 ≥ -D ν ν t δ t-1 -B ν ν t δ t-1
, by Jensen's inequality, Cauchy-Schwarz inequality, Hölder's inequality, and Assumption 1-p with p = 2. Hence, applying the inequalities ( 18) and ( 19) to ( 17), yields

δ t ≤ [1 -2µγ t + 2D ν ν t γ t + 2κ 2 t γ 2 t ]δ t-1 + 2B ν ν t γ t δ 1/2 t-1 + 2σ 2 t γ 2 t ≤ [1 -(µ -2D ν ν t )γ t + 2κ 2 t γ 2 t ]δ t-1 + B 2 ν µ ν 2 t γ t + 2σ 2 t γ 2 t , using Young's inequality 5 ; 2B ν ν t γ t δ 1/2 t-1 ≤ µγ t δ t-1 + B 2 ν ν 2 t γ t /µ.
Next, we introduce the indicator function for whether (ν t ) is constant (= C) or not (¬C), such that

δ t ≤ [1 -(µ ν -1 {ν t ¬C} 2D ν ν t )γ t + 2κ 2 t γ 2 t ]δ t-1 + B 2 ν µ ν 2 t γ t + 2σ 2 t γ 2 t , (20)
with µ ν = µ -1 {ν t =C} 2D ν ν t > 0. Let C δ be the constant fulfilling the conditions of Proposition 1 such that C δ is chosen larger than 1 verifying C δ (1 {ν t ¬C} 2D ν ν t + 2κ 2 t γ t ) ≤ µ ν /2 such that it's imply µ ν γ t /2 ≤ 1, which is possible as the sequence (ν t ) is non-increasing, and (κ t ) and (γ t ) is decreasing. At last, bounding (20) by Proposition 1 concludes the proof.

Proof of Theorem 1. Inserting the functions γ

t = C γ n β t t -α , ν t = n -ν t , κ t = C κ n -κ t , σ t = C σ n -σ
t , and n t = C ρ t ρ into the bound of Lemma 1 yields

δ t ≤π t + 2 1+2ρν B 2 ν µµ ν C 2ν ρ t 2ρν + 2 2+ρ(2σ-β)+α C 2 σ C γ C β ρ µ ν C 2σ ρ t ρ(2σ-β)+α (21) ≤π t + 2 (2+6ρν)/(1+ρ) B 2 ν µµ ν C 2ν/(1+ρ) ρ N 2ρν/(1+ρ) t + 2 (7+6ρσ)/(1+ρ) C 2 σ C γ µ ν C (2σ-β-α)/(1+ρ) ρ N (ρ(2σ-β)+α)/(1+ρ) t , with µ ν = µ -1 {ρ=0} 2D ν C -ν
ρ > 0, and π t can be bounded by exp

        - µ ν C γ C β ρ 2 t i=t/2 i ρβ-α                 exp        1 {ρ 0} 2C δ D ν C γ C β ρ C ν ρ t i=1 i ρ(β-ν)-α        exp         2C δ C 2 κ C 2 γ C 2β ρ C 2κ ρ t i=1 i 2ρ(β-κ)-2α                δ 0 + 2B 2 ν µµ ν C 2ν ρ + 4C 2 σ C γ C β ρ µ ν C 2σ ρ        + B 2 ν C γ C β ρ µC 2ν ρ t/2-1 i=1 i ρ(β-2ν)-α + 2C 2 σ C 2 γ C 2β ρ C 2σ ρ t/2-1 i=1 i 2ρ(β-σ)-2α         ≤ exp        - µ ν C γ C β ρ t 1+ρβ-α 2 2                exp        1 {ρ 0} 2C δ D ν C γ C β ρ ψ α-ρ(β-ν) (t) C ν ρ        exp         4(α -ρ(β -κ))C δ C 2 κ C 2 γ C 2β ρ (2α -2ρ(β -κ) -1)C 2κ ρ                δ 0 + 2B 2 ν µµ ν C 2ν ρ + 4C 2 σ C γ C β ρ µ ν C 2σ ρ        + B 2 ν C γ C β ρ ψ α-ρ(β-2ν) (t/2) µC 2ν ρ + 4(α -ρ(β -σ))C 2 σ C 2 γ C 2β ρ (2α -2ρ(β -σ) -1)C 2σ ρ         ≤ exp        - µC γ N (1+ρβ-α)/(1+ρ) t 2 (3+ρ(2+β)-α)/(1+ρ) C (1-β-α)/(1+ρ) ρ                exp         1 {ρ 0} 2C δ D ν C γ C β ρ ψ ρ α-ρ(β-ν) (2N t /C ρ ) C ν ρ         exp         4(α -ρ(β -κ))C δ C 2 κ C 2 γ C 2β ρ (2α -2ρ(β -κ) -1)C 2κ ρ                δ 0 + 2B 2 ν µµ ν C 2ν ρ + 4C 2 σ C γ C β ρ µ ν C 2σ ρ        + B 2 ν C γ C β ρ ψ ρ α-ρ(β-2ν) (N t /C ρ ) µC 2ν ρ + 4(α -ρ(β -σ))C 2 σ C 2 γ C 2β ρ (2α -2ρ(β -σ) -1)C 2σ ρ         , (22)
with help of an integral test for convergence6 , ψ x (t) and ψ y x (t) from ( 16), and by use of (N t /2C ρ ) 1/(1+ρ) ≤ t ≤ (2N t /C ρ ) 1/(1+ρ) .

Proofs for Section 3.2

As in Section 6.1, we begin the following sections by conducting a general study for any (γ t ), (ν t ), (κ t ), (σ t ), and (n t ), when applying the Polyak-Ruppert averaging estimate ( θt ) from (3). Moreover, we need to study fourth-order rate ∆ t = E[ θ t -θ * 4 ] of the recursive estimates (1) and (2).

Lemma 2 (Fourth-order moment). Assume that Assumptions 1-p to 3-p hold true for p = 4. Suppose that

µ ν = µ -1 {ν t =C} 2D 4 ν ν 4 t /µ 3 > 0. Let 1 {ν t =C} and 1 {ν t ¬C} indicate whether (ν t ) is constant or not. Denote ∆ t = E[ θ t -θ * 4 ]
for some ∆ 0 ≥ 0, where (θ t ) follows the recursion in (1) or (2). For any learning rate (γ t ), we have

∆ t ≤Π t + 4B 4 ν µ 3 µ ν max t/2≤i≤t ν 4 i + 1024 µµ ν max t/2≤i≤t σ 4 i γ 2 i + 96 µ ν max t/2≤i≤t σ 4 i γ 3 i , with Π t = exp         - µ ν 4 t i=t/2 γ i                exp        1 {ν t ¬C} C ∆ D 4 ν µ 3 t i=1 ν 4 i γ i        exp        256C ∆ µ t i=1 κ 4 i γ 3 i        exp        24C ∆ t i=1 κ 4 i γ 4 i        ∆ 0 + 4B 4 ν µ 3 µ ν max 1≤i≤t ν 4 i + 1024 µµ ν max 1≤i≤t σ 4 i γ 2 i + 96 µ ν max 1≤i≤t σ 4 i γ 3 i + B 4 ν µ 3 t/2-1 i=1 ν 4 i γ i + 256 µ t/2-1 i=1 σ 4 i γ 3 i + 24 t/2-1 i=1 σ 4 i γ 4 i         .
Proof of Lemma 2. The derivation of the recursive step sequence for the fourth-order moment ∆ t of (1) follows the same methodology as for the second-order moment in Lemma 1. In the same way we deduced ( 17), we can take the quadratic norm on (1), expand the norm, and take the square on both sides, to derive the equation

θ t -θ * 4 =( θ t-1 -θ * 2 + γ 2 t ∇ θ l t (θ t-1 ) 2 -2γ t ∇ θ l t (θ t-1 ), θ t-1 -θ * ) 2 = θ t-1 -θ * 4 + γ 4 t ∇ θ l t (θ t-1 ) 4 + 4γ 2 t ∇ θ l t (θ t-1 ), θ t-1 -θ * 2 + 2γ 2 t θ t-1 -θ * 2 ∇ θ l t (θ t-1 ) 2 -4γ t θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * -4γ 3 t ∇ θ l t (θ t-1 ) 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * .
Taking the expectation on both sides of the equality above gives us

∆ t =∆ t-1 + γ 4 t E[ ∇ θ l t (θ t-1 ) 4 ] + 4γ 2 t E[ ∇ θ l t (θ t-1 ), θ t-1 -θ * 2 ] + 2γ 2 t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ) 2 ] -4γ t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ] -4γ 3 t E[ ∇ θ l t (θ t-1 ) 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ] ≤∆ t-1 + γ 4 t E[ ∇ θ l t (θ t-1 ) 4 ] + 6γ 2 t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ) 2 ] -4γ t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ] + 4γ 3 t E[ θ t-1 -θ * ∇ θ l t (θ t-1 ) 3 ],
by use of Cauchy-Schwarz inequality. Next, Young's inequality yields

4γ 3 t θ t-1 -θ * ∇ θ l t (θ t-1 ) 3 ≤2γ 4 t ∇ θ l t (θ t-1 ) 4 + 2γ 2 t θ t-1 -θ * 2 ∇ θ l t (θ t-1 ) 2 , 8γ 2 t θ t-1 -θ * 2 ∇ θ l t (θ t-1 ) 2 ≤(µγ t /2) θ t-1 -θ * 4 + 32µ -1 γ 3 t ∇ θ l t (θ t-1 ) 4 ,
which helps us to obtain the simplified expression,

∆ t ≤[1 + µγ t /2]∆ t-1 + 3γ 4 t E[ ∇ θ l t (θ t-1 ) 4 ] + 32µ -1 γ 3 t E[ ∇ θ l t (θ t-1 ) 4 ] -4γ t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ].
To bound the fourth-order term E[ ∇ θ l t (θ t-1 ) 4 ], we make use of the Lipschitz continuity of ∇ θ l t (Assumption 2-p), Assumption 3-p, and that θ t-1 is F t-1 -measurable (Assumption 1-p), to have that

E[ ∇ θ l t (θ t-1 ) 4 ] ≤ 8κ 4 t ∆ t-1 + 8σ 4 t , (23)
using that x + y p ≤ 2 p-1 ( x p + y p ) for any p ∈ N. Thus,

∆ t ≤[1 + µγ t /2 + 256µ -1 κ 4 t γ 3 t + 24κ 4 t γ 4 t ]∆ t-1 + 256µ -1 σ 4 t γ 3 t + 24σ 4 t γ 4 t -4γ t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ]. (24)
Next, using the same arguments as in the proof of Lemma 1, Young's inequality, and Assumption 1-p with p = 4, we have

4γ t E[ θ t-1 -θ * 2 E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ), θ t-1 -θ * ] ≥ -4γ t E[ θ t-1 -θ * 3 E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ) ] ≥ -3µγ t ∆ t-1 -µ -3 γ t E[ E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ) 4 ] ≥ -3µγ t ∆ t-1 -µ -3 γ t D 4 ν ν 4 t ∆ t-1 -µ -3 γ t B 4 ν ν 4
t , such that the last term of ( 24) can be bounded as follows,

4γ t E[ θ t-1 -θ * 2 ∇ θ l t (θ t-1 ), θ t-1 -θ * ] = 4γ t E[ θ t-1 -θ * 2 E[∇ θ l t (θ t-1 )|F t-1 ], θ t-1 -θ * ] = 4γ t E[ θ t-1 -θ * 2 ∇ θ L(θ t-1 ), θ t-1 -θ * ] + 4γ t E[ θ t-1 -θ * 2 E[∇ θ l t (θ t-1 )|F t-1 ] -∇ θ L(θ t-1 ), θ t-1 -θ * ] ≥ µγ t ∆ t-1 -µ -3 γ t D 4 ν ν 4 t ∆ t-1 -µ -3 γ t B 4 ν ν 4 t .
Indeed, inserting this into (24) gives us

∆ t ≤ 1 - µ 2 - D 4 ν ν 4 t µ 3 γ t + 256κ 4 t γ 3 t µ + 24κ 4 t γ 4 t ∆ t-1 + B 4 ν ν 4 t γ t µ 3 + 256σ 4 t γ 3 t µ + 24σ 4 t γ 4 t ,
which can be modified with use the indicator function that determines whether (ν t ) is constant (= C) or not (¬C), such that

∆ t ≤ 1 - µ ν 2 - 1 {ν t ¬C} D 4 ν ν 4 t µ 3 γ t + 256κ 4 t γ 3 t µ + 24κ 4 t γ 4 t ∆ t-1 + B 4 ν ν 4 t γ t µ 3 + 256σ 4 t γ 3 t µ + 24σ 4 t γ 4 t , (25)
with

µ ν = µ -1 {ν t =C} 2D 4 ν ν 4 t /µ 3 > 0.
Note that µ ν from Lemma 1 is lower bounded by µ ν , and strictly lower bounded for (ν t ) constant, i.e., µ ν > µ ν > 0. Let C ∆ ≥ 1 fulfill the conditions of Proposition 1; the

C ∆ constant is chosen such that C ∆ (1 {ν t ¬C} D 4 ν ν 4 t /µ 3 + 256κ 4 t γ 2 t /µ + 24κ 4 t γ 3 t ) ≤ µ ν /2 implying µ ν γ t /2 ≤ 1,
which is possible as the sequence (ν t ) is non-increasing, and (κ t ) and (γ t ) decrease. Hence, by applying Proposition 1 on (25), we obtain the desired bound for ∆ t .

Corollary 1 (Fourth-order moment). Assume that Assumptions 1-p to 3-p hold true for p = 4.

Let γ t = C γ n β t t -α , ν t = n -ν t , κ t = C κ n -κ t , and σ t = C σ n -σ t with ν ∈ (0, ∞), β ∈ [0, 1], κ, σ ∈ [0, 1/2], and C γ , C κ , C σ > 0. Suppose n t = C ρ t ρ with ρ ∈ [0, 1) and C ρ ∈ N, such that µ ν = µ -1 {ρ=0} 2D 4 ν /µ 3 C 4ν ρ > 0. Denote ∆ t = E[ θ t -θ * 4
] for some ∆ 0 ≥ 0, where (θ t ) follows the recursion in (1) or (2). For αρβ ∈ (1/2, 1), we have

∆ t ≤ Π t + 2 2+4ρν B 4 ν µ 3 µ ν C 4ν ρ t 4ρν + 2 2ρ(2σ-β)+2α (2 10 µ -1 + 2 7 C γ C β ρ )C 4 σ C 2 γ C 2β ρ µ ν C 4σ ρ t 2ρ(2σ-β)+2α , ( 26
)
with Π t given in (27) such that Π t = O(exp(-N (1+ρβ-α)/(1+ρ) t )).
Proof of Corollary 1. Inserting the functions

γ t = C γ n β t t -α , ν t = n -ν t , κ t = C κ n -κ t , σ t = C σ n -σ
t , and n t = C ρ t ρ into the bound of Lemma 2 and using

γ 3 t ≤ C γ C β ρ γ 2 t as α -ρβ ∈ (1/2, 1), yields (26) with µ ν = µ -1 {ρ=0} 2D 4 ν /µ 3 C 4ν ρ > 0, where Π t can be bounded as follows, exp         - µ ν C γ C β ρ 4 t i=t/2 i ρβ-α                 exp        1 {ρ 0} C ∆ D 4 ν C γ C β ρ µ 3 C 4ν ρ t i=1 i ρ(β-4ν)-α        exp         2 8 C ∆ C 4 κ C 3 γ C 3β ρ µC 4κ ρ t i=1 i ρ(3β-4κ)-3α         exp         24C ∆ C 4 κ C 4 γ C 4β ρ C 4κ ρ t i=1 i 4ρ(β-κ)-4α                 ∆ 0 + 4B 4 ν µ 3 µ ν C 4ν ρ + 1024C 4 σ C 2 γ C 2β ρ µµ ν C 4σ ρ + 96C 4 σ C 3 γ C 3β ρ µ ν C 4σ ρ         + B 4 ν C γ C β ρ µ 3 C 4ν ρ t/2-1 i=1 i ρ(β-4ν)-α + 256C 4 σ C 3 γ C 3β ρ µC 4σ ρ t/2-1 i=1 i ρ(3β-4σ)-3α + 24C 4 σ C 4 γ C 4β ρ C 4σ ρ t/2-1 i=1 i 4ρ(β-σ)-4α         ≤ exp        - µ ν C γ C β ρ t 1+ρβ-α 2 3                exp         1 {ρ 0} C ∆ D 4 ν C γ C β ρ ψ 0 α-ρ(β-4ν) (t) µ 3 C 4ν ρ         exp         2 10 C ∆ C 4 κ C 3 γ C 3β ρ µC 4κ ρ         exp         2 6 C ∆ C 4 κ C 4 γ C 4β ρ C 4κ ρ                 ∆ 0 + 2 2 B 4 ν µ 3 µ ν C 4ν ρ + 2 10 C 4 σ C 2 γ C 2β ρ µµ ν C 4σ ρ + 2 7 C 4 σ C 3 γ C 3β ρ µ ν C 4σ ρ         + B 4 ν C γ C β ρ ψ 0 α-ρ(β-4ν) (t/2) µ 3 C 4ν ρ + 2 10 C 4 σ C 3 γ C 3β ρ µC 4σ ρ + 2 6 C 4 σ C 4 γ C 4β ρ C 4σ ρ         , ( 27
)
with help of the integral test for convergence; t i=1 i ρ(3β-4x)-3α ≤ 3 < 2 2 and t i=1 i 4ρ(β-x)-4α ≤ 2 for any x ≥ 0 as αρβ ∈ (1/2, 1). Lemma 3. Assume that Assumptions 1-p to 3-p for p = 4 and Assumption 4 hold true. Denote δt = E[ θtθ * 2 ] with θn given by (3), where (θ t ) follows the recursion in (1) or (2). In addition, Assumption 5 must hold true if (θ t ) follows the recursion in (2), which is indicated by 1 {D Θ <∞} . For any learning rate (γ t ), we have

δ1/2 t ≤ Λ 1/2 N t        t i=1 n 2(1-σ) i        1/2 + C 1/2 σ µN t        t i=1 n 2(1-σ-σ ) i        1/2 + 2 1/2 B 1/2 ν µN t         t j=2         n j ν j j-1 i=1 n i σ i                 1/2 + 1 µN t t-1 i=1 δ 1/2 i n i+1 γ i+1 - n i γ i + n t µγ t N t δ 1/2 t + n 1 µN t 1 γ 1 + 2 1/2 (C ∇ + κ 1 ) δ 1/2 0 + 2 1/2 µN t         t-1 i=1 n 2 i+1 (C 2 ∇ + κ 2 i+1 )δ i         1/2 + 2 3/4 µN t         t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν )n j+1 ν j+1 j-1 i=0 (C ∇ + κ i+1 )n i+1 δ 1/2 i                 1/2 + C ∇ µN t t-1 i=0 n i+1 ∆ 1/2 i , with Λ = Tr(∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ) and C ∇ = C ∇ /2 + 1 {D Θ <∞} 2G Θ /D 2
Θ . Proof of Lemma 3. The proof is divided into two parts; in the first part, (θ t ) follows (1), and the second part considers (2). Assume that (θ t ) is derived from the recursion in (1): following [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], we rewrite (1) to

1 γ t (θ t-1 -θ t ) = ∇ θ l t (θ t-1 ), (28)
where

∇ θ l t (θ t-1 ) = n -1 t n t i=1 ∇ θ l t,i (θ t-1 ). Observe that ∇ 2 θ L(θ * )(θ t-1 -θ * ) = -∇ θ l t (θ * ) + ∇ θ l t (θ t-1 ) -[∇ θ l t (θ t-1 ) -∇ θ l t (θ * ) -∇ θ L(θ t-1 )] -[∇ θ L(θ t-1 ) -∇ 2 θ L(θ * )(θ t-1 -θ * )],
where ∇ 2 θ L(θ * ) is invertible with lowest eigenvalue greater than µ, i.e., ∇ 2 θ L(θ * ) ≥ µI d . Thus, summing the parts, taking the quadratic norm and expectation, and using Minkowski's inequality, gives us 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 . ( 29
)
As (∇ θ l t (θ * )) is a square-integrable sequences on R d (Assumption 1-p), 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 2 i E ∇ 2 θ L (θ * ) -1 ∇ θ l i (θ * ) 2 + 2 N 2 t 1≤i< j≤t n i n j E ∇ 2 θ L (θ * ) -1 ∇ θ l i (θ * ) , ∇ 2 θ L (θ * ) -1 ∇ θ l j (θ * ) ,
where the first term can be bounded by Assumption 4,

1 N 2 t t i=1 n 2 i E ∇ 2 θ L (θ * ) -1 ∇ θ l i (θ * ) 2 ≤ 1 N 2 t t i=1 n 2(1-σ) i       Tr ∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 + C σ µ 2 n 2σ i       = Λ N 2 t t i=1 n 2(1-σ) i + C σ µ 2 N 2 t t i=1 n 2(1-σ-σ ) i , where Λ denotes Tr[∇ 2 θ L(θ * ) -1 Σ∇ 2 θ L(θ * ) -1 ]. For the next term, 2 N 2 t 1≤i< j≤t n i n j E ∇ 2 θ L (θ * ) -1 ∇ θ l i (θ * ) , ∇ 2 θ L (θ * ) -1 ∇ θ l j (θ * ) ≤ 2 µ 2 N 2 t 1≤i< j≤t n i n j E ∇ θ l i (θ * ) , ∇ θ l j (θ * ) = 2 µ 2 N 2 t 1≤i< j≤t n i n j E ∇ θ l i (θ * ) , E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * ) ≤ 2 µ 2 N 2 t 1≤i< j≤t n i n j E ∇ θ l i (θ * ) [E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * )] ≤ 2 µ 2 N 2 t 1≤i< j≤t n i n j E ∇ θ l i (θ * ) 2 E [E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * )] 2 ≤ 2B ν µ 2 N 2 t 1≤i< j≤t n i n j σ i ν j = 2B ν µ 2 N 2 t t j=2         n j ν j j-1 i=1 n i σ i         ,
by Cauchy-Schwarz inequality, Hölder's inequality, and Assumptions 1-p and 3-p. Thus,

         E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ * ) 2                   1/2 ≤ Λ 1/2 N t        t i=1 n 2(1-σ) i        1/2 + C 1/2 σ µN 1/2 t        t i=1 n 2(1-σ-σ ) i        1/2 + 2 1/2 B 1/2 ν µN t         t j=2         n j ν j j-1 i=1 n i σ i                 1/2 . (30)
Next, by the relation in ( 28), we have

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 notation of δ t = E[ θ t -θ * 2 ], the second term can be bounded by

         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 µγ t N t δ 1 2 t + n 1 µγ 1 N t δ 1 2 0 . (31)
For the third term, we can derive it as

E          ∇ 2 θ L(θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L(θ i-1 ) 2          = 1 µ 2 N 2 t t i=1 n 2 i E ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L(θ i-1 ) 2 + 2 µ 2 N 2 t t i< j n i n j E ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L(θ i-1 ), ∇ θ l j (θ j-1 ) -∇ θ l j (θ * ) -∇ θ L(θ j-1 ) , where t i=1 n 2 i E ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L (θ i-1 ) 2 ≤2 t i=1 n 2 i E ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) 2 + 2 t i=1 n 2 i E ∇ θ L (θ i-1 ) 2 ≤2 t i=1 n 2 i κ 2 i δ i-1 + 2C 2 ∇ t i=1 n 2 i δ i-1 ,
by the Cauchy-Schwarz inequality, Assumption 2-p and (8). For the other term, we note that

E[ ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L(θ i-1 ), ∇ θ l j (θ j-1 ) -∇ θ l j (θ * ) -∇ θ L(θ j-1 ) ] =E[ ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -[∇ θ L(θ i-1 ) -∇ θ L(θ * )], E[∇ θ l j (θ j-1 )|F j-1 ] -∇ θ L(θ j-1 ) -[E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * )] ] ≤ E[ ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -[∇ θ L(θ i-1 ) -∇ θ L(θ * )] 2 ] E[ E[∇ θ l j (θ j-1 )|F j-1 ] -∇ θ L(θ j-1 ) -[E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * )] 2 ] ≤ 2E[ ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) 2 ] + 2E[ ∇ θ L(θ i-1 ) -∇ θ L(θ * ) 2 ] 2E[ E[∇ θ l j (θ j-1 )|F j-1 ] -∇ θ L(θ j-1 ) 2 ] + 2E[ E[∇ θ l j (θ * )|F j-1 ] -∇ θ L(θ * ) 2 ] ≤ 2κ 2 i δ i-1 + 2C 2 ∇ δ i-1 2D 2 ν ν 2 j δ j-1 + 4B 2 ν ν 2 j ≤ 2 1/2 (κ i δ 1/2 i-1 + C ∇ δ 1/2 i-1 )(D ν ν j δ 1/2 j-1 + 2 1/2 B ν ν j ),
using F i-1 ⊂ F j-1 since i < j, Cauchy-Schwarz inequality, Hölder's inequality, a + b p ≤ 2 p-1 ( a p + b p ) with p ∈ N, Assumptions 1-p and 2-p, and (8). Thus,

         E          ∇ 2 θ L(θ * ) -1 1 N t t i=1 n i ∇ θ l i (θ i-1 ) -∇ θ l i (θ * ) -∇ θ L(θ i-1 ) 2                   1/2 ≤ 2 1/2 µN t        t i=1 n 2 i κ 2 i δ i-1        1/2 + 2 1/2 C ∇ µN t        t i=1 n 2 i δ i-1        1/2 + 2 3/4 µN t         t j=2         (D ν δ 1/2 j-1 + 2 1/2 B ν )n j ν j j-1 i=1 (C ∇ + κ i )n i δ 1/2 i-1                 1/2 . ( 32
)
The last term is directly bounded by ( 9), using that (9 [START_REF] Nesterov | Cubic regularization of newton method and its global performance[END_REF], giving us

) implies ∀θ, ∇ θ L(θ) -∇ 2 θ L(θ * )(θ -θ * ) ≤ C ∇ θ -θ * 2 /2 (Nesterov
         E          ∇ 2 θ L (θ * ) -1 1 N t t i=1 n i ∇ θ L (θ i-1 ) -∇ 2 θ L (θ * ) (θ i-1 -θ * ) 2                   1 2 ≤ C ∇ 2µN t t i=1 n i ∆ 1/2 i-1 , (33)
with the notion ∆ t = E[ θ t -θ * 4 ]. Combining the terms (30) to (33) into (29), gives us

δ1/2 t ≤ Λ 1/2 N t        t i=1 n 2(1-σ) i        1/2 + C 1/2 σ µN t        t i=1 n 2(1-σ-σ ) i        1/2 + 2 1/2 B 1/2 ν µN t         t j=2         n j ν j j-1 i=1 n i σ i                 1/2 + 1 µN t t-1 i=1 δ 1/2 i n i+1 γ i+1 - n i γ i + n t µγ t N t δ 1/2 t + n 1 µγ 1 N t δ 1/2 0 + 2 1/2 µN t             1/2 + C ∇ 2µN t t i=1 n i ∆ 1/2 i-1 ,
which gives the desired by shifting the indices and collecting the δ 0 terms,

δ1/2 t ≤ Λ 1/2 N t        t i=1 n 2(1-σ) i        1/2 + C 1/2 σ µN t        t i=1 n 2(1-σ-σ ) i        1/2 + 2 1/2 B 1/2 ν µN t         t j=2         n j ν j j-1 i=1 n i σ i                 1/2 + 1 µN t t-1 i=1 δ 1/2 i n i+1 γ i+1 - n i γ i + n t µγ t N t δ 1/2 t + n 1 µN t 1 γ 1 + 2 1/2 (C ∇ + κ 1 ) δ 1/2 0 + 2 1/2 µN t         t-1 i=1 n 2 i+1 (C 2 ∇ + κ 2 i+1 )δ i         1/2 + 2 3/4 µN t         t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν )n j+1 ν j+1 j-1 i=0 (C ∇ + κ i+1 )n i+1 δ 1/2 i                 1/2 + C ∇ 2µN t t-1 i=0 n i+1 ∆ 1/2 i . (34)
Now, assume that (θ t ) is derived from the recursion in (2): as above, we follow the steps of [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], in which, we can rewrite (2) to

1 γ t (θ t-1 -θ t ) = ∇ θ l t (θ t-1 ) - 1 γ t Ω t ,
where

Ω t = P Θ (θ t-1 -γ t ∇ θ l t (θ t-1 )) -(θ t-1 -γ t ∇ θ l t (θ t-1 )
). Thus, summing the parts, taking the norm and expectation, and using the Minkowski's inequality, yields the same terms as in ( 29), 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 ) Θ} , (35)
using Godichon-Baggioni (2016, Lemma 4.3). Next, we note that E[ Ω t 2 1 {θ t-1 -γ t ∇ θ l t (θ t-1 ) Θ} ] = 4γ 2 t G 2 Θ P[θ t-1γ t ∇ θ l t (θ t-1 ) Θ], since Ω t 2 = P Θ (θ t-1 -γ t ∇ θ l t (θ t-1 )) -θ t-1 + γ t ∇ θ l t (θ t-1 ) 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 (θ) 2 ≤ G 2 Θ for any θ ∈ Θ. Moreover, as in Godichon-Baggioni and Portier (2017, Theorem 4.2) with use of Lemma 2, we know that P[θ t-1 -γ t ∇ θ l t (θ t-1 ) Θ] ≤ ∆ t /D 4 Θ , where D Θ = inf θ∈∂Θ θθ * with ∂Θ denoting the frontier of Θ. Thus, (35) 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 ≤ 2G Θ µD 2 Θ N t t i=1 n i+1 ∆ 1/2 i ,
since the sequence (n t ) is either constant or increasing, meaning ∀t, n t /n t+1 ≤ 1. At last, this term can be combined into (34) with use of C ∇ = C ∇ /2 + 1 {D Θ <∞} 2G Θ /D 2 Θ , which indicates whether (θ t ) follows (2) or not. Proof of Theorem 2. The result follows by simplifying and bounding each term of Lemma 3, with use of Theorem 1 and Lemma 2. Thus, by inserting γ t = C γ n β t t -α , ν t = n -ν t , κ t = C κ n -κ t , σ t = C σ n -σ t , and n t = C ρ t ρ into the bound of Lemma 3, we obtain

δ1/2 t ≤ Λ 1/2 N 1/2 t 1 {σ=1/2} + Λ 1/2 C 1-σ ρ N t        t i=1 i 2ρ(1-σ)        1/2 1 {σ 1/2} + C 1/2 σ C 1-σ-σ ρ µN t        t i=1 i 2ρ(1-σ-σ )        1/2 + 2 1/2 B 1/2 ν C 1/2 σ C ρ µC (σ+ν)/2 ρ N t         t j=2         j ρ(1-ν) j-1 i=1 i ρ(1-σ)                 1/2 + (ρ(1 -β) + α)C ρ µC γ C β ρ N t t-1 i=1 i ρ(1-β)+α-1 δ 1/2 i + C ρ t ρ(1-β)+α µC γ C β ρ N t δ 1/2 t + C ρ µN t        1 C γ C β ρ + 2 1/2 C κ C κ ρ + C ∇        δ 1/2 0 + 2 1/2+ρ(1-κ) C κ C ρ µC κ ρ N t         t-1 i=1 i 2ρ(1-κ) δ i         1/2 + 2 1/2+ρ C ∇ C ρ µN t         t-1 i=1 i 2ρ δ i         1/2 + 2 3/4+ρ(2-ν)/2 C ρ µC ν/2 ρ N t         t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν ) j ρ(1-ν) j-1 i=1 C ∇ + 2 ρκ C κ C κ ρ i ρκ i ρ δ 1/2 i                 1/2 + 2 ρ C ∇ C ρ µN t t-1 i=0 i ρ ∆ 1/2 i ,
using n i+1 /n i ≤ 2 ρ and that |n i+1 /γ i+1n i /γ i | ≤ (ρ(1 -β) + α)C 1-β ρ /C γ i 1-ρ(1-β)-α as ρ(1 -β) + α ≤ 1 -ρ with ρ ∈ [0, 1). Next, as σ ∈ [0, 1/2] and σ ∈ (0, 1/2], we have t i=1 i 2ρ(1-σ-σ ) ≤ t 1+2ρ(1-σ-σ ) /(1 + 2ρ(1 -σσ )), where t ≤ (2N t /C ρ ) 1/(1+ρ) . Similarly, as ν ∈ (0, ∞), we have that

t-1 j=2         j ρ(1-ν) j-1 i=1 i ρ(1-σ)         ≤ t-1 j=1 j ρ(1-ν)
t-1 i=1 i ρ(1-σ) ≤ ψ ρ(ν-1) (t)ψ ρ(σ-1) (t) ≤ ψ ρ ρ(ν-1) (2N t /C ρ )ψ ρ ρ(σ-1) (2N t /C ρ ), using the ψ-function defined in ( 16), such that ψ ρ ρ(σ-1) (2N t /C ρ )ψ ρ ρ(ν-1) (2N t /C ρ )/N t ≤ Õ(N -ρ(σ+ν)/2(1+ρ) t

). Let D κ ∇ denote C ∇ + 2 ρκ C κ /C κ ρ with κ ∈ [0, 1/2], such that

2 1/2+ρ(1-κ) C κ C ρ µC κ ρ N t         t-1 i=1 i 2ρ(1-κ) δ i         1/2 + 2 1/2+ρ C ∇ C ρ µN t         t-1 i=1 i 2ρ δ i         1/2 ≤ 2 1/2+ρ D κ ∇ C ρ µN t         t-1 i=1 i 2ρ δ i         1/2
, and, likewise, we have that

t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν ) j ρ(1-ν) j-1 i=1 C ∇ + 2 ρκ C κ C κ ρ i ρκ i ρ δ 1/2 i         ≤ D κ ∇ t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν ) j ρ(1-ν) j-1 i=1 i ρ δ 1/2 i         .
From ( 21) we know that δ t ≤ D δ /t δ with

D δ = sup t∈N π t t δ + 2 1+2ρν B 2 ν µµ ν C 2ν ρ + 2 2+ρ(2σ-β)+α C 2 σ C γ C β ρ µ ν C 2σ ρ ,
and δ = 1 {B ν =0} (ρ(2σ -β) + α) + 1 {B ν 0} min{ρ(2σ -β) + α, 2ρν}, yielding

t-1 j=1         (D ν δ 1/2 j + 2 1/2 B ν ) j ρ(1-ν) j-1 i=1 i ρ δ 1/2 i         ≤ D 1/2 δ t-1 j=1         (D ν D 1/2 δ j -δ/2 + 2 1/2 B ν ) j ρ(1-ν) j-1 i=1 i ρ-δ/2         ≤ D 1/2 δ t-1 j=1 (D ν D 1/2 δ j -δ/2 + 2 1/2 B ν ) j ρ(1-ν) ψ δ/2-ρ (t) ≤ D ν D δ ψ δ/2-ρ (t)ψ δ/2+ρ(ν-1) (t) + 2 1/2 B ν D 1/2 δ ψ δ/2-ρ (t)ψ ρ(ν-1) (t) ≤ D ν D δ ψ ρ δ/2-ρ (2N t /C ρ )ψ ρ δ/2+ρ(ν-1) (2N t /C ρ ) + 2 1/2 B ν D 1/2 δ ψ ρ δ/2-ρ (2N t /C ρ )ψ ρ ρ(ν-1) (2N t /C ρ ),
if δ/2-ρ ≥ 0. Hence, ψ ρ δ/2-ρ (2N t /C ρ )ψ ρ δ/2+ρ(ν-1) (2N t /C ρ )/N t = Õ(N -(δ+ρν)/2(1+ρ) t

), and ψ ρ δ/2-ρ (2N t /C ρ )ψ ρ ρ(ν-1) (2N t /C ρ )/N t = Õ(N -(δ/2+ρν)/2(1+ρ) t ). Next, we define πt = t i=1 i 2 π i ≥ t i=1 π i such that π t ≤ t -1 t i=1 π i ≤ t -1 πt ≤ t -1 π∞ since π t is decreasing. Similarly, let Πt = t i=1 i ρ Π i . Both πt and Πt convergences to some finite constant depending on the

Figure 1 :

 1 Figure 1: Simulation of various data streams n t = C ρ t ρ . See Section 4 for details.

  (a) AR(1): well-specified case. See Section 4.1.1 for details. (b) AR(1): misspecified case. See Section 4.1.2 for details. (c) ARCH(1). See Section 4.2 for details. (d) AR(1)-ARCH(1). See Section 4.3 for details.

A convex body in R d is a compact convex set with non-empty interior.

You can't beat the system.

Later, in Section 3.2 for the averaged estimate (3), we assume in (8) that the function L has C ∇ -Lipschitz continuous gradients.

A modification of our average estimate to a weighted average version could improve convergence as it could limit the effect of poor initializations[START_REF] Boyer | On the asymptotic rate of convergence of stochastic newton algorithms and their weighted averaged versions[END_REF][START_REF] Mokkadem | A generalization of the averaging procedure: The use of two-time-scale algorithms[END_REF]. But despite this, we still achieve better convergence for the ASSG method.

If a, b, c > 0, p, q > 1 such that 1/p + 1/q = 1, then ab ≤ a p c p /p + b q /qc q .

Note that t i=1 i 2ρ(β-κ)-2α ≤ (2α -2ρ(β -κ))/(2α -2ρ(β -κ) -1) and t i=1 i 2ρ(β-σ)-2α ≤ (2α -2ρ(β -σ))/(2α -2ρ(β -σ) -1) as ν > 0, σ, κ ∈ [0, 1/2], ρ ∈ [0, 1), β ∈ [0, 1],and αρβ ∈ (1/2, 1).

model's parameters. With use of these notions, one can show that δ1/2

, and Ψ t given as

Furthermore, with Õ-notation one can yield,

ρ µµ 1/2 ν N (1+ρ(2σ-β)+α)/(2(1+ρ)) ), implying that ν > 1/2 to obtain the desired rate δt = O(N -1 ) if B ν = 0.