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

Wei Lu

Bruno Portier

Online estimation of the inverse of the Hessian for stochastic optimization with application to universal stochastic Newton algorithms

Keywords: Stochastic Newton algorithm, Stochastic Optimization, Robbins-Monro algorithm, online estimation

come

Online estimation of the inverse of the Hessian for stochastic optimization with application to universal stochastic Newton algorithms

Introduction

In this paper, we consider the usual stochastic optimization problem, which consists of estimating the parameter θ ∈ R d defined by θ = arg min

h∈R d G(h)
where the function G is defined for all h ∈ R d by : G(h) = E [g(X, h)] and X is a random vector of R p . The function g : R p × R d -→ R is assumed to be twice continuously differentiable. This problem arises in various contexts such as estimating the parameters of logistic regressions [START_REF] Bach | Adaptivity of averaged stochastic gradient descent to local strong convexity for logistic regression[END_REF][START_REF] Cohen | On projected stochastic gradient descent algorithm with weighted averaging for least squares regression[END_REF], geometric median and quantiles [START_REF] Cardot | Efficient and fast estimation of the geometric median in Hilbert spaces with an averaged stochastic gradient algorithm[END_REF](Cardot et al., , 2017)), or superquantiles (Bercu et al., 2020a;[START_REF] Costa | Non asymptotic controls on a recursive superquantile approximation[END_REF]. We denote ∇ h g and ∇ 2 h g as the gradient and Hessian matrix of g with respect to the second variable h, and ∇G and ∇ 2 G as the gradient and Hessian matrix of G. It is assumed that the matrix ∇ 2 G(θ) is positive definite.

Starting from a sequence of independent random vectors (X n ) n≥1 with the same distribution as X, we aim to online estimate the parameter θ. One of the most well-known methods in this context is certainly the stochastic gradient algorithm, recursively defined for all n ≥ 1 by:

θ SG n = θ SG n-1 -ν n ∇ h g X n , θ SG n-1
where θ SG 0 is an arbitrarily chosen initial value and (ν n ) n≥1 is a sequence of positive real numbers decreasing towards 0. These algorithms have been extensively studied, with asymptotic results found by [START_REF] Pelletier | On the almost sure asymptotic behaviour of stochastic algorithms[END_REF][START_REF] Pelletier | Asymptotic almost sure efficiency of averaged stochastic algorithms[END_REF] and non-asymptotic results, such as uniform bounds of the quadratic mean error, presented by [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF]; [START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF]; [START_REF] Godichon-Baggioni | Convergence in quadratic mean of averaged stochastic gradient algorithms without strong convexity nor bounded gradient[END_REF] to name a few. To ensure asymptotic efficiency, an additional step consists of considering an averaged version of the estimates [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF].

Despite their known efficiency, these methods can be very sensitive to ill-conditioned problems, where the Hessian has eigenvalues at different scales [START_REF] Leluc | Asymptotic optimality of conditioned stochastic gradient descent[END_REF][START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF]. To overcome this problem, second-order stochastic algorithms of the form

θ n = θ n-1 -ν n A n ∇ h g(X n , θ n-1 )
have been proposed and recently studied. Here, (ν n ) n≥1 is a sequence of positive real numbers decreasing towards 0 and the matrix A n is a recursive estimate of the inverse of the Hessian matrix of G at θ, i.e a recursive estimate of H -1 with H = ∇ 2 G(θ). The challenge lies in constructing the recursive estimate A n .

Several recursive second-order algorithms have been proposed and studied. For example, [START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF] propose an efficient stochastic Newton algorithm for estimating the parameters of a logistic regression model. In a recent work, [START_REF] Bercu | A stochastic gauss-newton algorithm for regularized semi-discrete optimal transport[END_REF] propose a stochastic Gauss-Newton algorithm to estimate the entropically regularized Optimal Transport cost between two discrete probability measures. [START_REF] Cénac | An efficient averaged stochastic gauss-newton algorithm for estimating parameters of non linear regressions models[END_REF] study the asymptotic properties of a stochastic Gauss-Newton algorithm for estimating the parameters of a non-linear regression model. [START_REF] Godichon-Baggioni | Recursive ridge regression using second-order stochastic algorithms[END_REF] propose second-order algorithms to solve the Ridge regression problem in the linear and logistic framework, while the case of the geometric median is introduced and studied by [START_REF] Godichon-Baggioni | Online stochastic newton methods for estimating the geometric median and applications[END_REF]. In all these algorithms, the estimate of the inverse of the Hessian matrix is recursively computed using the Riccati inversion formula (also called Sherman-Morrison formula, see e.g. Duflo (1990) p. 96). This calculation is made possible thanks to the particular form of the estimate of the Hessian matrix H, presented as (1/n) n k=1 a k ϕ k ϕ T k , where (a n ) n≥1 is a sequence of positive real random variables and (ϕ n ) n≥1 is a sequence of random vectors in R d .

However, it is not always possible to obtain such an estimate of the Hessian matrix. In this work, we propose to construct a direct recursive estimate of H -1 without first attempting to construct an estimate of H. This approach is based on the fact that we have HH -1 = H -1 H = I d and, consequently, the following relation:

E H -1 ∇ 2 h g(X, θ) + ∇ 2 h g(X, θ)H -1 -2I d = 0 (1)
where I d denotes the identity matrix of order d. Using a Robbins-Monro type algorithm, we propose a recursive estimate of the matrix H -1 defined for all n ≥ 1 by:

A n = A n-1 -γ n A n-1 ∇ 2 h g(X n , θ n-1 ) + ∇ 2 h g(X n , θ n-1 )A n-1 -2 I d
where (γ n ) n≥1 is a sequence of positive real numbers, decreasing towards 0 and θ n-1 is an estimate of θ.

However, the complexity of computing this estimate is of order O(d 3 ), which is the same as directly calculating the inverse of an estimate of matrix H. Nevertheless, we can introduce an algorithm with complexity of order O(d 2 ) based on the following observation: let Z be a centered random vector in R d with variance-covariance matrix I d , independent of the vector X. Then,

E H -1 Z Z T ∇ 2 h g(X, θ) + ∇ 2 h g(X, θ)Z Z T H -1 -2I d = 0. (2)
Therefore, considering a sequence (Z n ) n≥1 of random vectors in R d independent of the sequence (X n ) n≥1 leads to an estimate of the form:

Framework

We consider the problem of minimizing the convex function G : R d -→ R defined for all h ∈ R d by:

G(h) := E [g(X, h)] ,
where the loss g(X, •) is a convex, twice-differentiable function and X is a random vector of R p . We assume that there exists a unique value θ ∈ R d such that

∇G(θ) = 0.
This assumption, couple with strict convexity, ensures the existence of a minimizer for G and provides a welldefined optimization problem. Now, let's introduce the assumptions that underlie the parameter estimation framework for θ:

(A1) There exists C > 0 such that for all h ∈ R d , E ∥∇ h g(X, h)∥ 2 ≤ C (1 + G(h) -G(θ)) .
(A2) The functional G is twice continuously differentiable and ∇ 2 G(θ) is positive. In addition, the Hessian is uniformly bounded, i.e there exists a positive constant L ∇G such that for all h ∈ R d ,

∇ 2 G(h) op ≤ L ∇G .
(A3) The function ∇ 2 G is Lipschitz on a neighborhood of θ, i.e. there exist positive constants r > 0 and L r such that for all h ∈ B(θ, r)

∇ 2 G(h) -∇ 2 G(θ) op ≤ L r ∥θ -h∥ ,
where B(θ, r) denotes a ball of radius r centered at θ. (A4) There exists q > 2 and C q such that for all h ∈ R d ,

E ∇ 2 h g(X, h) q F ≤ C q .
These assumptions are very close to those presented in the literature [START_REF] Pelletier | Asymptotic almost sure efficiency of averaged stochastic algorithms[END_REF][START_REF] Gadat | Optimal non-asymptotic bound of the ruppert-polyak averaging without strong convexity[END_REF][START_REF] Godichon-Baggioni | Online estimation of the asymptotic variance for averaged stochastic gradient algorithms[END_REF]. Assumption (A1) controls the growth of the gradient and guarantees the stability of the estimation process. It ensures that the gradient remains bounded as the estimation progresses. Assumption (A2) ensures that the curvature of G at θ is well-behaved, allowing the estimation algorithm to reliably exploit the local structure of G. This assumption also guarantees that the gradient of G is Lipschitz continuous with a constant L ∇G . This Lipschitz continuity is crucial as convergence results are obtained using a Taylor's expansion of G up to the second order. Assumption (A3) indicates that the Hessian matrix does not exhibit abrupt changes within a neighborhood around θ. This assumption ensures the stability of the Hessian estimates during the estimation process. Assumption (A4) guarantees that the second-order derivative of G does not exhibit excessive fluctuations. It imposes a bound on the variation of the Hessian matrix, providing further stability to the estimation algorithm. It's worth noting that Hölder's inequality leads to

∇ 2 G(h) F ≤ E ∇ 2 h g(X, h) q F
1/q ≤ C 1/q q . Of course, this inequality intertwines with the bound in Assumption (A2), but we keep the notation L ∇G for the sake of clarity.

These assumptions, along with the differentiability properties of g and G, provide a solid foundation for developing efficient second order methods to solve the minimization problem and obtain reliable estimates of θ.

Estimation of the Hessian inverse

In this section, our focus is solely on estimating the inverse Hessian of function G in θ, denoted as H -1 with H = ∇ 2 G(θ). Even if our motivation is to estimate H -1 for proposing a second-order algorithm, the estimation of H -1 can also be valuable for recursively constructing confidence intervals or significance statistical tests for a component of parameter θ when the parameter θ is estimated using another asymptotically efficient algorithm like the averaged stochastic gradient algorithm. Indeed, in most cases, the asymptotic variance involving in the central limit theorem, generally depends on matrix H -1 and its estimation is then required.

Let (X n ) n≥1 be a sequence of independent random vectors in R p with the same distribution as X. Assume first that θ is known. From equality (1), the matrix H -1 satisfies an equation of the form Φ(H -1 ) = 0. We can then employ the Robbins-Monro procedure [START_REF] Robbins | A stochastic approximation method[END_REF] to recursively estimate the parameter H -1 . Denoting this estimator as A n , for any n ≥ 1, we have:

A n = A n-1 -γ n A n-1 ∇ 2 h g(X n , θ) + ∇ 2 h g(X n , θ) A n-1 -2 I d ,
where A 0 is an arbitrary symmetric positive definite matrix, and γ n = c γ n -γ with 1 2 < γ < 1 and c γ > 0. It is important to note that A n is symmetric for any n ≥ 1 due to its construction. However, since θ is unknown, we need to estimate it. Assuming we have an efficient recursive estimator of θ (e.g., a stochastic gradient estimator), we can easily derive an estimator of H -1 using a plug-in procedure:

A n = A n-1 -γ n A n-1 ∇ 2 h g(X n , θ n-1 ) + ∇ 2 h g(X n , θ n-1 ) A n-1 -2 I d .
This estimator is always symmetric but not necessarily positive definite. To ensure positive definiteness, we introduce a truncation based on the norm of ∇ 2 h g(X n , θ n-1 ), leading to the following estimator of H -1 :

A n = A n-1 -γ n A n-1 ∇ 2 h g(X n , θ n-1 ) + ∇ 2 h g(X n , θ n-1 ) A n-1 -2I d 1 ∥∇ 2 h g(Xn, θ n-1 )∥ op ≤βn , (3)
where β n = c β n β with 1-γ q-1 < β < γ -1 2 and 0 < c γ c β < 1 2 . Additionally, this truncation enables control over the smallest eigenvalue of A n , which is useful for studying an estimator of the parameter θ involving the matrix A n . This is particularly important in establishing the consistency of the Stochastic Newton algorithm presented in Section 4.

However, although this estimator is efficient, each update requires matrix multiplications, resulting in a computational complexity of order O(d 3 ), which is the same as matrix inversion. Hence, it is necessary to improve the complexity of each update of A n .

Building on equality (2), considering a sequence (Z n ) n≥1 of independent and identically distributed bounded random vectors of R d such that E [Z n ] = 0 and E Z n Z T n = I d , and independent of (X n ) n≥1 , we can propose another estimate of H -1 defined for any n ≥ 1 as follows:

P n = A n-1 Z n Q n = ∇ 2 h g(X n , θ n-1 )Z n A n = A n-1 -γ n P n Q T n + Q n P T n -2 I d 1 {∥Qn∥∥Zn∥≤βn} (4)
where A 0 is an arbitrary symmetric and positive definite matrix. We can observe that in this algorithm the truncation is only based on ∥Q n ∥ ∥Z n ∥, and not on

Q n Z T n op
as expected, because we have Q n Z T n op ≤ ∥Q n ∥ ∥Z n ∥. Notably, the computational complexity of each update of A n is now reduced to O(d 2 ). Moreover, following [START_REF] Mokkadem | A generalization of the averaging procedure: The use of two-timescale algorithms[END_REF]; [START_REF] Boyer | On the asymptotic rate of convergence of stochastic newton algorithms and their weighted averaged versions[END_REF], we can propose a weighted averaged estimate A n,τ , which performs better in practice when the initialization is poor. It is given by

A n,τ = 1 - ln(n + 1) τ n k=0 ln(k + 1) τ A n-1,τ + ln(n + 1) τ n k=0 ln(k + 1) τ A n .
(5)

This estimator can be recursively computed since for any n ≥ 1, n k=0 ln(k + 1) τ = ln(n + 1) τ + n-1 k=0 ln(k + 1) τ . The following theorem establishes the consistency of the estimators A n and A n,τ for the parameter H -1 in the context of estimating the inverse of the Hessian matrix. It states that the convergence rate depends on multiple factors, including the step sequence γ n , the truncation parameter β n , the regularization parameter τ , and the convergence rate of the estimate θ n .

Theorem 3.1. Assume that Assumptions (A2) to (A4) hold, and that there is an estimate θn satisfying for all δ > 0 θn -

θ 2 = o ln n 1+δ n a a.s.,
with a > 0. Then A n and A n,τ defined by (4) and (5) satisfy

A n -H -1 2 = o ln n 1+δ n min{γ,a,2β(q-1)} a.s. and

A n,τ -H -1 2 = o ln n 1+δ n min{1,a,2β(q-1)} a.s.

The proof is given in Section 6. Observe that if a = 1, one can achieve the usual rate of convergence taking β > 1 2(q-1) , which is only possible if q > 1 + 1 2γ-1 . For instance, taking the usual parametrization γ = 3/4 or 2/3, q must satisfy q > 3 or q > 4.

Universal Weighted Averaged Stochastic Newton Algorithm

In this section, we introduce the Universal Weighted Averaged Stochastic Newton algorithm and discuss its main properties. As mentioned in Theorem 3.1, using a Weighted Averaged version of the inverse Hessian estimate can yield improved theoretical results. Therefore, we incorporate this choice into the Stochastic Newton algorithm. Furthermore, we have observed that the convergence rate of the estimate for θ significantly influences the theoretical behavior of the estimates of H -1 . Consequently, we incorporate the best possible estimate for parameter θ, namely the Weighted Averaged Stochastic Newton estimates, into the latter. This reasoning leads to the following Weighted Averaged Stochastic Newton algorithm defined for all n ≥ 1 by

P n = A n-1 Z n Q n = ∇ 2 h g(X n , θ n-1,τ ′ )Z n θ n = θ n-1 -ν n A n-1,τ ∇ h g(X n , θ n-1 ) (6) θ n,τ ′ = 1 - ln(n + 1) τ ′ n k=0 ln(k + 1) τ ′ θ n-1,τ ′ + ln(n + 1) τ ′ n k=0 ln(k + 1) τ ′ θ n (7) A n = A n-1 -γ n P n Q T n + Q n P T n -2 I d 1 {∥Qn∥∥Zn∥≤βn} (8) A n,τ = (1 - ln(n + 1) τ n k=0 ln(k + 1) τ )A n-1,τ + ln(n + 1) τ n k=0 ln(k + 1) τ A n (9)
where (ν n ) n≥1 is a sequence of positive real numbers defined for any n ≥ 1 by ν n = c ν n ν with c ν > 0 and ν ∈ (1/2, 1 -β) satisfying γ + ν > 3/2. In addition, τ, τ ′ ≥ 0. The following theorem gives the consistency of the estimates defined by ( 6) and (7).

Theorem 4.1. Assume that Assumptions (A1) to (A4) hold. Let θ n and θ n,τ ′ be defined as in (6) and (7). Then, θ n a.s.

---→ n→∞ θ and θ n,τ ′ a.s.

---→ n→∞ θ.

The proof is given in Section 6. Note that the constraint γ + ν > 3/2 is of a purely technical nature and is crucial for the application of the Robbins-Siegmund Theorem and so that to get the consistency of the estimates. However, we believe this condition might not be necessary in practical applications. We can now give the almost sure rate of convergence of the estimates.

Theorem 4.2. Assume that Assumptions (A1) to (A4) hold. Then θ n and θ n,τ ′ defined by (6) and (7) satisfy for all δ > 0

∥θ n -θ∥ 2 = o ln n 1+δ n ν a.s. and θ n,τ ′ -θ 2 = o ln n 1+δ n a.s.
In addition, A n and A n,τ defined by (8) and (9) satisfy

A n -H -1 2 = o ln n 1+δ n γ a.s. and A n,τ -H -1 2 = o ln n 1+δ n a.s.
Moreover, the estimates θ n,τ ′ defined by (7) satisfy

√ n(θ n,τ ′ -θ) L ---→ n→∞ N (0, H -1 ΣH -1 ),
where

Σ = E ∇ h g(X, θ)∇ h g(X, θ) T .
The proof is given in Section 6. The Universal Weighted Averaged Stochastic Newton estimates so achieve the asymptotic efficiency, and so, under very weak assumptions.

Remark 4.1. Mention that for estimating parameter θ, it is also possible to consider the following simpler algorithm, which we refer to as the Universal Stochastic Newton Algorithm, and only relies on A n :

P n = A n-1 Z n Q n = ∇ 2 h g(X n , θ n-1 )Z n A n = A n-1 -γ n P n Q T n + Q n P T n -2I d 1 {∥ Qn∥∥Zn∥}≤βn θ n = θ n-1 -ν n A n-1 ∇ h g(X n , θ n-1 ).
By following the same scheme of proof as for Theorem 4.2, one could check that:

θ n -θ 2 = O ln n n ν a.s.
However, it can be observed that the convergence rate of θ n is not optimal. Nevertheless, this algorithm has the merit of being much simpler. In addition, mention that following the reasoning presented by [START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF], one could take a step sequence of the form ν n = 1 n leading to the Stochastic Newton algorithm. However, we are unfortunately not able to obtain the consistency of the estimates in this context.

Applications

The simulation section of this paper focuses on evaluating the performance of our novel methods, called Universal Stochastic Newton Algorithm (USNA) and Universal Weighted Averaged Stochastic Newton Algorithms (UWASNA). We begin by analyzing the performance of USNA and UWASNA in the context of logistic regression and the geometric median estimation. In both contexts, it was already feasible to employ second-order algorithms such as the Stochastic Newton Algorithms (SNA) and its Weighted Averaged version (WASNA). These two algorithms use the Riccati formula to recursively compute the inverse of the Hessian estimator. By demonstrating comparable results to SNA and WASNA, we establish USNA and UWASNA as viable alternatives with efficient performance. Additionally, we investigate the applicability of our method in challenging scenarios where using the Riccati formula in SNA is not feasible, particularly in estimating p-means and parameters of a spherical distribution. In these two cases, we compare the performances of USNA and UWASNA with the one of the Averaged Stochastic Gradient Descent (ASGD) introduced by [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF]. Through comprehensive simulations, our goal is to verify that USNA and UWASNA consistently exhibit favourable performances even in these contexts. To close this section, we extend our evaluation to real-world datasets to showcase the practicality of deploying USNA and UWASNA algorithm in applications.

Choice of the hyperparameters

In our experiments, the choice of hyperparameters involved in the different algorithms, plays a significant role in achieving desired outcomes. The decision to set these values is based on both theoretical justifications and empirical observations. 1. Setting of ν n for USNA: Despite the lack of a theoretical proof demonstrating the convergence rate of USNA when ν n = 1/n, this setting is adopted in our experiments for a direct comparison with SNA. Empirical results, as presented later, validate that this choice is effective in practice.

2. Conditions on β and γ: Condition β < γ -1/2 is only used to apply the Robbins-Siegmund theorem. It serves a theoretical purpose, and we advocate for its removal. On the contrary, condition β n+1 γ n+1 ≤ 1/2 is essential. It ensures the positiveness of the estimate of the inverse of the Hessian.

For this reason, in all simulations, we set

β n+1 = 1 2 γ -1 n .
3. Initialization of Estimators: For initializing the estimators of the Hessian inverse, we consistently use S -1 0 = I d for both SNA and WASNA. Similarly, A 0 = I d is chosen for USNA and UWASNA.

Comparison with Riccati Newton

The objective of this section is to demonstrate the comparable performance of USNA and UWASNA when contrasted with SNA and WASNA, particularly in scenarios where the use of the Riccati formula is applicable to recursively compute the inverse of an estimator of the Hessian. Let us begin by revisiting this context.

If we can estimate the Hessian matrix H = ∇ 2 G(θ) by an estimator of the form S n /n with S n defined by S n = n j=1 φ j φ T j where (φ n ) n≥1 is a sequence of random vectors of R d , then, thanks to the Riccati formula, we can recursively calculated matrix S -1 n for any n ≥ 1:

S -1 n = S -1 n-1 - 1 1 + φ T n S -1 n-1 φ n S -1 n-1 φ n φ T n S -1 n-1 , (10)
with S 0 = I d to avoid the invertibility problem. This formula finds application in various scenarios, as previously demonstrated by [START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF], for instance, to obtain efficient stochastic Newton algorithms.

In light of this, we can define the stochastic Newton algorithm (SN) for estimating parameter θ as followed:

U n = S -1 n-1 φ n (11) S -1 n = S -1 n-1 -(1 + φ T n U n ) -1 U n U T n (12) θ SN n = θ SN n-1 -S -1 n ∇ h g(X n , θ SN n-1 ) (13)
where S -1 0 = I d and θ SN 0 is arbitrarily chosen. Note that the random vector φ n is dependent on the current observation X n and the previous estimation θ SN n-1 . The Weighted Averaged Stochastic Newton Algorithm is defined by:

U n = S -1 n-1 φ n S -1 n = S -1 n-1 -(1 + φ T n U n ) -1 U n U T n θ n = θ n-1 -γ n+1 S -1 n ∇ h g(X n , θ n-1 ) θ W ASN n = 1 - ln(n + 1) τ ′ n k=0 ln(k + 1) τ ′ θ W ASN n-1 + ln(n + 1) τ ′ n k=0 ln(k + 1) τ ′ θ n
where S -1 0 = I d , θ 0 and θ W ASN 0 are arbitrarily chosen, the random vector φ n is dependent on the current observation X n and the previous estimation θ W ASN n-1 .

Logistic regression

Let (X, Y ) be a random vector taking values in R p × {0, 1} and let us set ϕ = (1, X T ) T . In the binary logistic regression framework, function G to minimize is defined for any h ∈ R p+1 by:

G(h) = E log(1 + exp(h T Φ)) -h T ΦY = E [g(X, Y, h)]
where the conditional distribution of the binary response Y knowing ϕ is a Bernoulli distribution of parameter π(θ T ϕ) with for any x ∈ R, π(x) = exp(x)/(1 + exp(x)) and θ ∈ R p+1 is the unknown parameter to be estimated. It is easy to show that θ = arg min

h∈R p+1 G(h)
and θ is the unique solution of equation ∇G(h) = 0. In addition, we have

H = E a(θ T ϕ)ϕϕ T with a(z) = π(z)(1 -π(z)).
Let (ϕ n , Y n ) n≥1 be a sequence of independent random vectors in R p+1 × {0, 1} with the same distribution as (ϕ, Y ). In this context, [START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF] uses algorithm ( 11)-( 13) to estimate θ with

φ n = a(ϕ T n θ SN n-1 )ϕ n and ∇ h g(X n , θ SN n-1 ) = -Φ n Y n -π(ϕ T n θ SN n-1
) . We conduct extensive simulations to evaluate the performance of our novel methods, USNA and UWASNA, in comparison to SNA and WASNA. For this purpose, we consider the logistic regression model introduced by [START_REF] Bercu | An efficient stochastic newton algorithm for parameter estimation in logistic regressions[END_REF] where p = 10 and the true coefficients are set as follows: θ = (0, 3, -9, 4, -9, 15, 0, -7, 1, 0) T .

We compare USNA and UWASNA against SNA and WASNA in terms of their ability to accurately estimate the true coefficients. The evaluation of the algorithms' performance is carried out using the mean squared error (MSE) metric. We simulate N = 100 independent sample of size n = 10 000. The results are averaged to mitigate the effects of sampling fluctuations. For all algorithms, we initialize the estimate of the parameter with θ init = θ + eϵ, where ϵ ∼ N (0, I p+1 ) and e = 1 or 2. As shown in the figure 1, the weighted averaged estimators converge more rapidly than the other two methods. Both USNA and UWASNA perform comparably to SNA and WASNA in accurately estimating the true coefficients of the logistic regression model. Remarkably, USNA and UWASNA achieve this without using the Riccati formula.

Geometric Median

Next, we conduct simulations in the context of geometric median estimation for a multivariate distribution. We focus on the model introduced by Godichon-Baggioni and Lu (2023). We generate n = 10 000 copies of the random vector X of R p with p = 10, where X ∼ N (0, Σ) with Σ i,j = |i -j| 0.5 . Recall that the geometric median is defined by: m = arg min

h∈R p E [∥X -h∥ -∥X∥] .
In this model, the result leads to m = (0, . . . , 0) T in this model. In addition, the Hessian matrix H is defined by:

H = E 1 ∥X -m∥ I p - (X -m)(X -m) T ∥X -m∥ 2 .
In this context, algorithm ( 11)-( 13) can be used to estimate

m taking ∇ h g(X n , m SN n-1 ) = -Xn-m n-1 ∥Xn-m n-1 ∥ and φ n = ∇ h g(X n , m SN n-1 + α n Z n ) -∇ h (X n , m SN n-1 ), where α n = 1 n ln (n+1) and (Z n ) n≥1 is a sequence of independent standard Gaussian vectors.
For a comprehensive evaluation of our methods' effectiveness in geometric median estimation, we also compare them against two baselines : SNA and WASNA. For all four algorithms, we initialize the estimate of the geometric median with m init = eϵ, where ϵ ∼ N (0, I p+1 ) and e = 1 or 2. Throughout the simulations, we recorded the MSE of the estimated medians for the four algorithms. The simulation results were averaged over multiple iterations (N = 100). In Figure 2, the results consistently indicate that when the sample size is relatively small (around 100), USNA and UWASNA converge slightly slower. However, beyond that point, they achieve performance on par with both WASNA and SNA in estimating the geometric median.

Cases where the Riccati formula cannot be used

In this section, we focus on the scenarios where using the Riccati formula in SNA or WASNA is not feasible. We explore the performance of USNA and UWASNA in contrast to an alternative averaged stochastic gradient-based method ("ASGD") proposed by [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF], which is defined as:

θ SGD n = θ SGD n-1 -η n ∇ h g(X n , θ SGD n-1 ) (14) θ ASGD n = θ ASGD n-1 + 1/n(θ SGD n -θ ASGD n-1 ) (15)
where (η n ) n≥1 is a sequence of learning rates and θ ASGD 0 = θ SGD 0 .

Spherical Distribution

In this paragraph, we focus on the estimation of the parameters of a spherical distribution (Godichon-Baggioni and Portier, 2017). The aim of the task is to fit a sphere onto a 3D point cloud with noise. In this context, we assume that the observations represent independent realizations of a random vector X, which is defined as

X = µ + rW U,
where U is uniformly distributed on the unit sphere of R 3 , W ∼ U ([1 -δ, 1 + δ]) with δ > 0, W and U are independent. The radius r > 0 and the center µ ∈ R 3 are parameters to be estimated. The unknown parameter θ = (µ, r) T is a local minimizer of the function

G : R 3 × R * + -→ R defined for all h = (a, b) ∈ R 3 × R * + by: G(h) := E [g(X, h)] = 1 2 E (∥X -a∥ -b) 2 .
In this scenario, one can use algorithm ( 14)-( 15) to estimate θ with

∇ h g(X, h) = a -X + b(X-a)
∥X-a∥ b -∥X -a∥ .

Second-order methods USNA and UWASNA can also be used with

∇ 2 h g(X, h) =   1 -b ∥X-a∥ I 3 + b(X-a)(X-a) T ∥X-a∥ 3 X-a ∥X-a∥ (X-a) T ∥X-a∥ 1   .
We emphasize that in this specific case, neither the conventional Stochastic Newton Algorithms (SNA) nor the Weighted Averaged version (WASNA) using the Riccati formula are applicable due to the nature of the problem. Thus, we conduct simulations to compare the performance of the USNA, UWASNA, and ASGD. The synthetic datasets were generated with a sample size of 10 000 and the true parameters of the spherical distribution were set as follows: µ = (0, 0, 0) and r = 2. In addition, we set δ = 0.2, which results in a Hessian matrix with eigenvalues of different order sizes. For all three algorithms, we initialize the estimate of the parameter by θ init = (µ init , r init ) T = (0, 0, 0, 2) T + eϵ where ϵ ∼ N (0, I 4 ) and e = 0.5 or 1. Multiple iterations (N = 100) were performed to reduce the impact of sampling variations. Figure 3 illustrates the comparison between the performances of USNA, UWASNA and ASGD in terms of the mean squared error (MSE) of the estimated parameters. Throughout the simulations, UWASNA amd USNA demonstrate superior performance in accurately estimating the parameters of the spherical distribution when compared to the gradient-based method. Additionally, the Hessian matrix H of the model can be explicitly calculated (Godichon-Baggioni and Portier, 2017), one has

H = I 3 -2 3 I 3 E [W ] W -1 0 0 1 .
Note that this matrix is diagonal, making its inverse computation straightforward. Therefore, we investigate the Frobenius norm of the difference between the estimated matrix A n and the true matrix H -1 . From Figure 4, it's evident that our methods provide a good estimation of the inverse of the Hessian matrix. Moreover, UWASNA offers a better estimation compared to USNA.

p-means

Now we focus on the estimation of p-means of a multivariate distribution [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF]. We consider a random vector X of R d with d = 40, where X ∼ N (0, Σ) with Σ i,j = |i -j| 0.5 . The p-mean m of X is defined as the minimizer of the functional G p : R d → R given for all h ∈ R d by: where 1 ≤ p < +∞. Note that in our model m = (0, . . . , 0) T . We can easily verify that the gradient and the Hessian of G p are given by:

G p (h) = 1 p E [∥X -h∥ p ] ,
∇G p (h) = -E (X -h) ∥X -h∥ p-2 , ∇ 2 G p (h) = E ∥X -h∥ p-2 I d -(2 -p) (X -h)(X -h) T ∥X -h∥ 2 .
We aim to compare the performance of USNA, UWASNA, and ASGD for estimating the p-mean of X through simulations. We consider the case p = 1.5 and we simulate N = 100 independent samples of size n = 10 000. For all three algorithms, we initialize the estimate with m init = eϵ where ϵ ∼ N (0, I d ) and e = 1 or 2.

As depicted in the Figure 5, we plotted the MSE versus sample size. It is evident that USNA and UWASNA consistently outperforms ASGD in estimating the p-means. Their superior performance can be attributed to the incorporation of information from the Hessian matrix. These results further emphasize the advantage of using USNA and UWASNA as alternative methods when it is impossible to use the Riccati formula.

Application to real data

We apply the algorithms to "COVTYPE" dataset, a well-known dataset used for classification tasks [START_REF] Blackard | Comparative accuracies of artificial neural networks and discriminant analysis in predicting forest cover types from cartographic variables[END_REF][START_REF] Schmidt | Minimizing finite sums with the stochastic average gradient[END_REF][START_REF] Toulis | Towards stability and optimality in stochastic gradient descent[END_REF]. The original study was based on a dataset comprising 581,011 observations and 54 predictors. The primary objective was to predict the cover type of forests within Roosevelt National Park. In our current investigation, we will narrow our focus to the most prevalent category of the target variable, which is "Spruce/Fir," accounting for 49% of the observations. As a result, we convert the "covertype" variable into a binary form, with the "fir" category assigned a value of 1, while all other categories are assigned a value of 0. The objective is to use logistic regression to predict the variable "covertype". We split the dataset into training (80%) and test (20%) sets. We implement the UWASNA, USNA, WASNA and SNA on the training set. As a baseline, we also apply a first-order algorithm ASGD on it. We evaluate the performance of each algorithm by calculating the accuracy of each one on both training set and test set. The results are summarized in Table 1 We observe that UWASNA, USNA, WASNA and SNA demonstrate similar performances, and they achieve higher accuracies on both training set and test set compared to ASGD. By successfully applying USNA and UWASNA on this dataset, we illustrate their practicality in real-world applications.

Conclusion

In this study, we thoroughly examined the stochastic optimization problem, primarily focusing on accurately estimating the unknown parameter. Our significant contribution is the introduction of a direct method to estimate the inverse of the Hessian matrix. Instead of the traditional approach, which estimates the Hessian matrix first, we directly addressed its inverse, using the Robbins-Monro procedure.

This approach led us to develop the Universal Weighted Averaged Stochastic Newton Algorithm. Through our extensive testing, we found that our newly proposed methods are both efficient and robust. When compared with standard first and second-order algorithms, our method consistently performed well. In certain scenarios, it even outperformed these conventional algorithms.

In summary, our findings emphasize the potential of second-order methods in optimization. Our approach to directly estimating the Hessian matrix's inverse represents a significant advancement, combining simplicity with efficiency.

Proofs

For the sake of simplicity, in the following we denote T n := ∇ 2 h g(X n , θn-1 )Z n Z T n . Recall that ∥Z n ∥ is bounded, i.e there is M such that ∥Z∥ ≤ M . In addition, let us recall that if ∥Q n+1 ∥ ∥Z n+1 ∥ ≤ β n+1 , it then follows that ∥T n+1 ∥ op ≤ β n+1 .

Study on the largest eigenvalue of A n

The following proposition provides an initial asymptotic bound of the largest eigenvalue of A n without requiring knowledge on the behavior of the estimate θn . This result is crucial to prove Theorems 3.1 and 4.1. Proposition 6.1. Under Assumptions (A3) and (A5), the largest eigenvalue of A n denoted by λ max (A n ) satisfies for all δ > 0

λ max (A n ) = o n 1-γ ln n 1+δ a.s. and λ max (A n , τ ) = o n 1-γ ln n 1+δ a.s. Proof of Proposition 6.1. Define W n := A n-1 T T n + T n A n-1 . By definition of A n , ∥A n+1 ∥ 2 F = ∥A n ∥ 2 F -2γ n+1 ⟨A n , W n+1 -2I d ⟩ F 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 + γ 2 n+1 ∥W n+1 -2I d ∥ 2 F 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 .
Therefore,

E ∥A n+1 ∥ 2 F |F n = ∥A n ∥ 2 F -2γ n+1 ⟨A n , ∇ 2 G( θn )A n + A n ∇ 2 G( θn ) -2I d ⟩ F + γ 2 n+1 E ∥W n+1 -2I d ∥ 2 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 |F n + 2γ n+1 E ⟨A n , W n+1 -2I d ⟩ F 1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 |F n . Assumption (A5) ensures that E ∥W n+1 -2I d ∥ 2 F 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 |F n ≤ 8 E ∇ 2 h g(X n+1 , θn ) 2 F E ∥Z n+1 ∥ 4 ∥A n ∥ 2 F + d ≤ 8 C 2/q q ∥A n ∥ 2 F + d .
Since ∥Q n+1 ∥ = ∥Q n+1 ∥ q ∥Q n+1 ∥ 1-q with q > 1, for any ζ n > 0,

γ n+1 E ⟨A n , W n+1 -2I d ⟩ F 1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 |F n ≤ 2γ n+1 ∥A n ∥ 2 F E ∥Q n+1 ∥ ∥Z n+1 ∥ 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 |F n + γ n+1 ζ n+1 ∥A n ∥ 2 F + ζ n+1 d ≤ 2C q M 2q γ n+1 β 1-q n+1 + γ n+1 ζ n+1 ∥A n ∥ 2 F + γ n+1 ζ n+1 d.
Since ∇ 2 G( θn ) and A n are positive, one has for any

ζ n > 0, -γ n+1 ⟨A n , ∇ 2 G( θn )A n + A n ∇ 2 G( θn ) -2I d ⟩ F ≤ γ n+1 ⟨A n , 2I d ⟩ ≤ γ n+1 ζ n+1 ∥A n ∥ 2 F + γ n+1 ζ n+1 d.
Finally,

E ∥A n+1 ∥ 2 F |F n ≤ 1 + 4γ n+1 ζ n+1 + 4C q γ n+1 β 1-q n+1 + 8γ 2 n+1 C 2/q q ∥A n ∥ 2 F + 4γ n+1 ζ n+1 d + 8γ 2 n+1 d.
Setting ζ n = n 1-γ ln n 1+δ with δ > 0, we can apply Lemma 6.2 with V n = ∥A n ∥ 2 F , and

a n = ζ 2 n . One then obtains ∥A n ∥ 2 F = o ζ 2 n = o n 2-2γ ln n 2+2δ , so that λ max (A n ) = o n 1-γ ln n 1+δ a.s.

Proof of Theorem 3.1

The aim is to provide an initial rate of convergence for A n . A more refined or faster rate will be established later. First, note that

A n+1 -H -1 = A n -H -1 -γ n+1 (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 = A n -H -1 -γ n+1 (W n+1 -2I d ) + γ n+1 (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 = A n -H -1 -γ n+1 (E [W n+1 |F n ] -2I d ) + γ n+1 ξ n+1 + γ n+1 (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 = A n -H -1 -γ n+1 ∇ 2 G( θn )A n + A n ∇ 2 G( θn ) -2I d + γ n+1 ξ n+1 + γ n+1 (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1
where ξ n+1 := -W n+1 + ∇ 2 G( θn )A n + A n ∇ 2 G( θn ). Let α * k be the function defined for all h ∈ M p (R) by

α * k (h) = (I d -γ k+1 H) h (I d -γ k+1 H) ,
we then have

A n+1 -A = α * n A n -H -1 -γ 2 n+1 H(A n -H -1 )H + γ n+1 r 1,n + γ n+1 r 2,n + γ n+1 s n + γ n+1 s ′ n + γ n+1 ξ n+1 (16)
where

r 1,n = (∇ 2 G( θn ) -H)H -1 + H -1 (∇ 2 G( θn ) -H), r 2,n = (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 -E (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 |F n , s n = E T n+1 A n -H -1 + A n -H -1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 |F n + ∇ 2 G( θn ) -H A n -H -1 + A n -H -1 ∇ 2 G( θn ) -H , s ′ n = E T n+1 H -1 + H -1 T n+1 -2I d 1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 |F n By induction, A n -H -1 = Ψ * n,0 A 0 -H -1 + n-1 k=0 Ψ * n,k+1 γ k+1 ξ k+1 + n-1 k=0 Ψ * n,k+1 γ k+1 r 1,k + n-1 k=0 Ψ * n,k+1 γ k+1 r 2,k + n-1 k=0 Ψ * n,k+1 γ k+1 s k + n-1 k=0 Ψ * n,k+1 γ k+1 s ′ k , (17)
where Ψ * n,k is the function defined for all h ∈ M p (R) by :

Ψ * n,k (h) =   n j=k+1 α * j   (h) =   n j=k+1 (I d -γ j H)   h   n j=k+1 (I d -γ j H)   .
Next, we will determine the rate of convergence for each term on the right-hand side of equation ( 17).

Rate of convergence of M n := n-1 k=0 Ψ * n,k+1 γ k+1 ξ k+1 : Recall that for all δ > 0, λ max (A n ) = o n 1-γ ln n 1+δ a.s. Thus, there exists a constant c ′ > 0 such that

E ∥ξ n+1 ∥ 2 F |F n ≤ c ′ 1 + ∥A n ∥ 2 F = o n 2-2γ
ln n 2+2δ a.s., and according to Lemma 6.3, it follows that

n-1 k=0 Ψ * n,k+1 γ k+1 ξ k+1 2 F = O n 2-3γ ln n 2+2δ a.s. ( 18
) Rate of convergence of R 1,n := n-1 k=0 Ψ * n,k+1 γ k+1 r 1,k : we have R 1,n+1 = (I d -γ n+1 H) R 1,n (I d -γ n+1 H) + γ n+1 r 1,n+1 .
Therefore, for n large enough

∥R 1,n+1 ∥ op ≤ ∥I d -γ n+1 H∥ 2 op ∥R 1,n ∥ op + γ n+1 ∥r 1,n+1 ∥ op ≤ (1 -λ min (H)γ n+1 ) 2 ∥R 1,n ∥ op + γ n+1 ∥r 1,n+1 ∥ op a.s.
By Assumption (A4), ∥r 1,n+1 ∥ op = o n -a/2 ln n (1+δ)/2 a.s. According to Lemma 6.4,

∥R 1,n+1 ∥ op = o n -a/2 ln n (1+δ)/2 a.s. ( 19
)
Rate of convergence of R 2,n := n-1 k=0 Ψ * n,k+1 γ k+1 r 2,k :

E (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 2 op |F n ≤12E ∥Q n+1 ∥ 2 ∥Z n+1 ∥ 2 ∥A n ∥ 2 + d 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 |F n .
Applying Markov's inequality, and since ∥Z n ∥ is supposed to be bounded by M ,

E (W n+1 -2I d )1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 2 |F n ≤ 12dE ∇ 2 g(X n+1 , θn ) q op Z n Z T n q |F n β -q n + 12 ∥A n ∥ 2 E ∥Q n ∥ 2 M 2 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 ≤ 12dM 2q E ∇ 2 g(X n+1 , θn ) q op |F n β -q n + 12C q ∥A n ∥ 2 β 2-q n+1
Thus, by Lemma 6.4,

∥R 2,n ∥ 2 op = O n 2-3γ ln n 2+2δ n 2β-qβ a.s. ( 20
)
Note that the rate of convergence for R 2,n is faster than that of M n .

Rate of convergence of S

′ n := n-1 k=0 Ψ * n,k+1 γ k+1 s ′ k : We have similarly s ′ n ≤ E T n+1 H -1 + H -1 T n+1 -2I d F 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 |F n ≤ 2 H -1 F E ∥Z n+1 ∥ 2 ∇ 2 h g X n+1 , θn F 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 |F n + 2qE 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 |F n = O β 1-q n .
Therefore, thanks to Lemma 6.4,

S ′ n = O β 1-q n . (21)
A first result for A n : Thanks to equalities ( 18) to ( 21), one can rewrite A n -H -1 as

A n -H -1 = n-1 k=0 Ψ * n,k+1 γ k+1 s k + S n with S n F = o n max(1-3 2 γ,-a/2,β(1-q)) ln n 1+δ .
In addition, we have

∥s n ∥ op ≤ A n -H -1 op E ∥Q n+1 ∥ ∥Z n+1 ∥ 1 ∥Q n+1 ∥∥Z n+1 ∥≥βn ∥F n + ∇ 2 G( θn ) -H A n -H -1 op = o A n -H -1 op a.s. Define S n := n-1 k=0 Ψ * n,k+1 γ k+1 s k .
There exists a positive sequence ( r n ) n≥0 with r n a.s.

---→ n→∞ 0 such that

∥S n+1 ∥ op ≤ (1 -λ min (H)γ n+1 ) 2 ∥S n ∥ op + γ n+1 r n A n -H -1 op ≤ (1 -λ min (H)γ n+1 ) 2 ∥S n ∥ op + γ n+1 r n ∥S n ∥ op + S n op .
Applying Lemma 6.4,

∥S n ∥ op = o n max(1-3 2 γ,-a/2,β(1-q)) ln n 1+δ , which implies that A n -H -1 F = o n max(1-3 2 γ,-a/2,β(1-q)) ln n 1+δ a.s.
Final rate of convergence of A n : The aim here is, with the help of this first result on A n , to give better rates of convergence for M n and R 2,n . First, note that if γ > 2/3, then 1 -3 2 γ < 0 and we have directly A n a.s.

---→ n→∞ H -1 . If γ ≤ 2/3, we have 1-3 2 γ > 0 , so that max{1-3 2 γ, -a/2, β(1-q)} = max{1-3 2 γ, β(1-q)} and ∥A n ∥ = o n max{1-3 2 γ,β(1-q)} ln n 1+δ a.s.. Thus, applying Lemma 6.3 in the worse case (i.e when γ > 2/3), one now has

∥M n ∥ 2 F = o n max{2-4γ,β(2-2q)-γ} ln n 2+2δ a.s. and ∥R 2,n ∥ 2 op = O n 2-4γ ln n 2+2δ n 2β-qβ a.s.
Following the same process as before, we now have that for any γ ∈ (1/2, 1), A n converges almost surely to H -1 , and in particular, one has ∥A n ∥ = O (1). Applying another time Lemma 6.3, we have

∥M n ∥ 2 F = O n -γ ln n 1+δ a.s. and ∥R 2,n ∥ 2 op = O n -γ ln n 2+2δ n 2β-qβ a.s.
Finally, we have

A n -H -1 2 F = o ln n 1+δ n min{γ,a,2β(q-1)} a.s. ( 22
)
Rate of convergence of A n,τ . Decomposition ( 16) can be written as

H(A n -H -1 ) + (A n -H -1 )H = A n -H -1 -(A n+1 -H -1 ) γ n+1 + r 1,n + r 2,n + s n + s ′ n + ξ n+1 , so that H(A n,τ -H -1 ) + (A n,τ -H -1 )H =t n n k=0 u k (A k -H -1 ) -(A k+1 -H -1 ) + t n n k=0 u k γ k+1 (r 1,n + r 2,n + s n + s ′ n + ξ n+1 ) (23)
where t n = 1 n k=0 log(k+1) τ and u k+1 = ln(k+1) τ γ k+1 . Our objective is to determine the rate of convergence for each term on the right-hand side of the decomposition given in (23).

Rate of convergence of t

n n k=0 u k+1 (A k -H -1 ) -(A k+1 -H -1
) : with the help of an Abel's transform,

t n n k=0 u k+1 (A k -H -1 ) -(A k+1 -H -1 ) = -t n A n+1 -H -1 u n+1 + t n A 0 -H -1 u 1 + t n n k=1 (A k -H -1 )(u k+1 -u k ).
It is obvious that t n A 0 -H -1 u 1 is negligible while thanks to equality (22),

t n A n+1 -H -1 u n+1 = o ln n 1+δ n 1-γ+ 1 2 min{γ,a,2β(q-1)} a.s
In addition, by the mean value theorem,

|u k+1 -u k | ≤ 2k γ-1 max{1, τ } ln(k + 1) τ , which implies that n k=1 (A k -H -1 )(u k+1 -u k ) F ≤ n k=1 A k -H -1 F 2k γ-1 max{1, τ } ln(k + 1) τ
and with the help of equation ( 22)

t n n k=0 u k (A k -H -1 ) -(A k+1 -H -1 ) F = o ln n 1+δ n 1-γ+ 1 2 min{γ,a,2β(q-1)} a.s
Rate of convergence of t n n k=0 u k+1 γ k+1 ξ k+1 . Remark that there exists c ′ > 0 such that for all δ > 0

E ∥ξ n+1 ∥ 2 F ∥F n ≤ c ′ 1 + ∥A n ∥ 2 F ≤ c ′ 1 + 2 A n -H -1 2 F + 2 H -1 2 F .
By applying a law of large numbers for martingales, one obtains for all δ > 0

t n n k=0 u k+1 γ k+1 ξ k+1 2 F = o ln n 1+δ n a.s.
Rate of convergence of t n n k=0 u k+1 γ k+1 r 1,k . First, let us give a useful equality: consider n k=0 k -ã/2 ln k (1+δ)/2 with ã > 0, since

n 1 x -ã/2 ln x (1+δ)/2 dx ≤ c(ln n) (1+δ)/2 x 1-ã if ã ̸ = 2 and n 1 x -ã/2 ln x (1+δ)/2 dx ≤ c ′ (ln n) (3+δ)/2 if ã = 2 then one can verify that t n n k=0 u k γ k+1 k -ã/2 ln k (1+δ)/2 2 F = o ln n 1+δ+21 ã=2 n min{2,ã} a.s. (24)
Then, recalling that ∥r 1,k+1

∥ op = o n -a/2 ln k (1+δ)/2 t n n k=0 u k γ k+1 r 1,k 2 F = o ln n 1+δ+21 a=2 n min{2,a} a.s. Rate of convergence of t n n k=0 u k+1 γ k+1 s ′ k . Recalling that s ′ k = O β 1-q k , it comes t n n k=0 u k γ k+1 s ′ k 2 F = o ln n 21 2β(q-1)=2 n min{2,2β(q-1)} a.s. Rate of convergence of t n n k=0 u k+1 γ k+1 s k . Recall that ∥s n ∥ op ≤ A n -H -1 op E ∥Q n+1 ∥ ∥Z n+1 ∥ 1 ∥Q n+1 ∥∥Z n+1 ∥≥βn ∥F n + ∇ 2 G( θn ) -H A n -H -1 op , Recalling that E ∥T n+1 ∥ op 1 ∥T n+1 ∥≥βn |F n = O β 1-q n
, and with the help of equality ( 22)

A n -H -1 op E ∥Q n+1 ∥ ∥Z n+1 ∥ 1 ∥Q n+1 ∥∥Z n+1 ∥≥βn ∥F n = o ln n (1+δ)/2 n β(q-1)+ 1 2 min{a,γ,2β(q-1)} a.s
In addition, one has

∇ 2 G( θn ) -H A n -H -1 op = o ln n 1+δ n a/2+ 1 2 min{a,γ,2β(q-1)} a.s.
leading to

∥s n ∥ = o ln n 1+δ
n min{a/2,β(q-1)}+ 1 2 min{a,γ,2β(q-1)} a.s.

Then, since min {a/2, β(q -1)} + 1 2 min {a, γ, 2β(q -1)} ≥ min {a, γ, 2β(q -1)},

t n n k=0 u k γ k+1 s k 2 F
= o ln n 1+δ n min{2,2 min{a,γ,2β(q-1)}} a.s.

Rate of convergence of t n n k=0 u k+1 γ k+1 r 2,k . Let us recall that r 2,n is a martingale difference, and with the help of a law of large numbers for martingales it comes that

t n n k=0 u k+1 γ k+1 r 2,k 2 F = o ln n 1+δ n a.s.
Rate of convergence of A n,τ . We so have

H(A n,τ -H -1 ) + (A n,τ -H -1 )H 2 = o ln n 1+δ n min{1,a,2β(q-1)} a.s.
Note that H is diagonalisable, and its eigenvalues are denoted by λ 1 , ..., λ d . Thus, for any matrix B ∈ M d , we have

∥HB + BH∥ 2 F = n i=1 n j=1 (HB + BH) 2 i,j = n i=1 n j=1 (λ i + λ j ) 2 B 2 i,j ≥ 4λ min (H) 2 n i=1 n j=1 B 2 i,j = 4λ min (H) 2 ∥B∥ 2 F . A n,τ -H -1 2 ≤ 1 4λ min (H) 2 H(A n,τ -H -1 ) + (A n,τ -H -1 )H 2 = o
ln n 1+δ n min{1,a,2β(q-1)} a.s.

Proof of Theorem 4.1

In order to prove this theorem, we first show that the eigenvalues of A n and A n,τ are well controlled, which implies the consistency of θ n . Then, we give an initial rate of convergence of θ n , which allows us to obtain the rate of convergence of A n and A n,τ . The last step is to deduce the rate of convergence of θ n and θ n,τ . We first study the smallest eigenvalue of A n . Proposition 6.2. Under Assumptions(A3) and (A5), the smallest eigenvalue of A n denoted by λ min (A n ) satisfies λ min (A n ) -1 = O n β a.s. and λ min (A n,τ ) -1 = O n β a.s.

Proof of Proposition 6.2. First of all, note that for all n ≥ 0, λ min (A n ) > 0. According to the definition of A n , one has

λ min (A n+1 ) = λ min A n -γ n+1 T n+1 A n + A n T T n+1 -2I d 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 ≥ λ min I d -γ n+1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 A n × I d -γ n+1 T T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 + γ n+1 I d -γ 2 n+1 T n+1 A n T T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 ≥ λ min I d -γ n+1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 A n × I d -γ n+1 T T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 + γ n+1 -γ 2 n+1 β 2 n+1 λ max (A n ) -γ n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1
For all h ∈ R p , one has

h t I d -γ n+1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 A n I d -γ n+1 T T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 h ≥ A 1/2 n I d -γ n+1 T n+1 1 ∥∥Q n+1 ∥∥Z n+1 ∥≤β n+1 ∥ h 2 ≥ λ min (A n ) I d -γ n+1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 h 2 ≥ λ min (A n ) (1 -γ n+1 β n+1 ) 2 ∥h∥ 2
Thus,

λ min I d -γ n+1 T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 A n I d -γ n+1 T T n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≤β n+1 ≥ λ min (A n ) (1 -γ n+1 β n+1 ) 2 . Therefore, λ min (A n+1 ) ≥ (1 -γ n+1 β n+1 ) 2 λ min (A n ) + γ n+1 -γ 2 n+1 β 2 n+1 λ max (A n ) -γ n+1 1 ∥Q n+1 ∥∥Z n+1 ∥>β n+1 .
Note that (β n γ n ) n is a decreasing sequence and β n γ n ≤ 1/2. Let U n and V n be sequences defined as

U n := n k=1 n j=k+1 (1 -γ k β k ) 2 γ 2 k+1 β 2 k+1 λ max (A n ),
and

V n := n k=0+1 n j=k+1 (1 -γ k β k ) 2 γ k+1 1 ∥Q n+1 ∥∥Z n+1 ∥≥β k+1 .
One can prove by induction that

λ min (A n ) ≥ λ min (A 0 ) n k=1 (1 -γ k β k ) 2 + n k=1 n j=k+1 (1 -γ j β j ) 2 γ k -U n -V n ,
with the convention n j=n+1 (1 -γ j β j ) 2 = 1. One has

U n+1 = (1 -γ n+1 β n+1 ) 2 U n + γ 2 n+1 β 2 n+1 λ max (A n ), since λ max (A n ) = o ln n 1+δ n 1-γ a.s., U n = o ln n 1+δ β n+1 n 1-γ γ n+1 = o 1 β n
since β < γ -1/2. One also has

V n+1 = (1 -γ n+1 β n+1 ) 2 V n + γ n+1 1 ∥Q n+1 ∥∥Z n+1 ∥≥β n+1 .
We define V ′ n := n v ln(n + 1) -(1+δ) V n for some v > 0, then, for n large enough,

E V ′ n+1 |F n ≤ ln n 1+δ ln(n + 1) 1+δ (n + 1) v n v (1 -γ n+1 β n+1 ) 2 V ′ n + γ n+1 (n + 1) v C q ln(n + 1) 1+δ β q n M 2q ≤ V ′ n + γ n+1 (n + 1) v C q ln(n + 1) 1+δ β q n M 2q .
In order to apply the Robbins-Siegmund theorem, let us take v = γ + qβ -1, then V ′ n converges almost surely to a finite random variable for all δ, which can be translated by

V n = o ln(n + 1) 1+δ n γ+qβ-1 a.s which is negligible as soon as β > 1-γ q-1 . It is obvious that λ min (A 0 ) n k=1 (1 -γ k β k ) 2 ≥ 0. Finally n k=1 n j=k+1 (1 -γ j β j ) 2 γ k ≥ n k=1 1 2β k n j=k+1 (1 -γ j β j ) 2 2γ k β k ≥ n k=1 1 2β k n j=k+1 (1 -γ j β j ) 2 2γ k β k -β 2 k γ 2 k = n k=1 1 2β k   n j=k+1 (1 -γ j β j ) 2 - n j=k (1 -γ j β j ) 2   Since 1
βn n is a decreasing sequence, one has

n k=1 n j=k+1 (1 -γ j β j ) 2 γ k ≥ 1 2β n n k=1   n j=k+1 (1 -γ j β j ) 2 - n j=k (1 -γ j β j ) 2   ≥ 1 2β n   1 - n j=1 (1 -γ j β j ) 2   ≥ 1 -(1 -γ 1 β 1 ) 2 2β n . Thus, λ min (A n ) ≥ c 1 βn + o 1 βn a.s with c 1 = 1 -(1 -γ 1 β 1 ) 2 /2, which implies that 1 λ min (A n ) = O (β n ) a.s.
We now prove the convergence of θ n . According to Taylor's theorem, and using Assumption (A2), we obtain

G(θ n+1 ) = G(θ n ) + ∇G(θ n ) T (θ n+1 -θ n ) + 1 2 (θ n+1 -θ n ) T 1 0 ∇ 2 G(θ n+1 + t(θ n -θ n+1 ))dt(θ n+1 -θ n ) ≤ G(θ n ) + ∇G(θ n ) T (θ n+1 -θ n ) + L ∇G 2 ∥θ n+1 -θ n ∥ 2 .
From the definition of θ n (6),

G(θ n+1 ) ≤ G(θ n ) -ν n+1 ∇G(θ n ) T A n,τ ∇g(X n+1 , θ n ) + L ∇G 2 ν 2 n+1 ∥A n,τ ∥ 2 op ∥∇g(X n+1 , θ n )∥ 2 .
Letting K n := G(θ n ) -G(θ), we can rewrite the inequality as

K n+1 ≤ K n -ν n+1 ∇G(θ n ) T A n,τ ∇ 2 g(X n+1 , θ n ) + L ∇G 2 ν 2 n+1 ∥A n,τ ∥ 2 op ∥∇g(X n+1 , θ n )∥ 2 ,
which implies

E [K n+1 |F n ] ≤ K n -ν n+1 ∇G(θ n ) T A n,τ ∇G (θ n ) + L ∇G 2 ν 2 n+1 ∥A n,τ ∥ 2 op E ∥∇g(X n+1 , θ n )∥ 2 |F n .
According to Assumption (A1), we have

E [K n+1 |F n ] ≤ 1 + CL ∇G 2 ν 2 n+1 ∥A n,τ ∥ 2 op K n -ν n+1 λ min (A n,τ ) ∥∇G(θ n )∥ 2 + CL ∇G 2 ν 2 n+1 ∥A n,τ ∥ 2 op .
With the help of Proposition 6.1,

∞ n=0 ν 2 n+1 ∥A n,τ ∥ 2 op < +∞ a.s.
Subsequently, according to Robbins-Siegmund Theorem, K n is guaranteed to converge almost surely to a finite random variable, and

∞ n=0 ν n+1 λ min (A n,τ ) ∥∇G(θ n )∥ 2 < +∞ a.s.
With the help of Proposition 6.2, ---→ n→∞ θ and θ n,τ ′ a.s.

---→ n→∞ θ.

Proof of Theorem 4.2

Rate of convergence of θ n . Recall that

E [G(θ n+1 ) -G(θ)|F n ] ≤ 1 + L ∇G C 2 ν 2 n+1 λ max (A n,τ ) 2 (G(θ n ) -G(θ)) -ν n+1 ∇G(θ n ) T A n,τ ∇G(θ n ) + L ∇G C 2 ν 2 n+1 λ max (A n,τ ) 2 .
According to equation ( 6) in Godichon-Baggioni (2021), there is a positive constant c 0 such that

∇G(θ n ) T A n,τ ∇G(θ n ) ≥ λ min (A n,τ ) ∥∇G(θ n )∥ 2 ≥ c 0 λ min (A n,τ ) (G(θ n ) -G(θ)) .
Recall that 2γ + 2ν -2 > 1, so that there exists η > 0 such that η < 2γ + 2ν -3. We define

V n := n η (G(θ n ) -G(θ)), so that E V n+1 |F n = n + 1 n η 1 + L ∇GC 2 ν 2 n+1 λ max (A n,τ ) 2 -c 0 ν n+1 λ min (A n,τ ) V n + L ∇GC 2 ν 2 n+1 λ max (A n,τ ) 2 (n + 1) η . Let U n := n+1 n η 1 + L ∇GC 2 ν 2 n+1 λ max (A n,τ ) 2 -c 0 ν 2 n+1 λ min (A n,τ ) , then E V n+1 |F n ≤ V n + L ∇GC 2 ν 2 n+1 λ max (A n,τ ) 2 (n + 1) η + V n 1 Un>1 .
As ν + β < 1, one can easily verify that 1 Un>1 a.s.

--→ 0. We can now apply the Robbins-Siegmund theorem. We therefore have ∥G(θ n+1 ) -G(θ)∥ = o (n -η ) for all η < 2γ + 2ν -3. Thanks to the local strong convexity of G, this leads to ∥θ n -θ∥ 2 = o (n -η ) a.s. We are now able to apply Theorem 3.1 and one obtains ∥A n -A∥ 2 = o ln n 1+δ n min{γ,η,2β(q-1)} a.s. and ∥A n,τ -A∥ 2 = o ln n 1+δ n min{η,2β(q-1)} a.s.

According to Theorem 4.2 and Theorem 4.3 in Boyer and Godichon-Baggioni (2022), one so has

∥θ n -θ∥ 2 = o ln n 1+δ n ν a.s. (25)
Rate of convergence of θ n,τ ′ . For all non-negative integers n,

θ n+1 -θ n = -ν n+1 A n,τ ∇ h g(X n+1 , θ n ), so that θ n -θ -(θ n+1 -θ) ν n+1 = θ n -θ n+1 ν n+1 = A n,τ ∇ h g(X n+1 , θ n ).
In addition,

θ n -θ = H -1 ∇ h g(X n+1 , θ n ) + H -1 ∇G(θ n ) -H -1 ∇ h g(X n+1 , θ n ) -H -1 ∇G(θ n ) -(θ n -θ) .
Finally,

θ n -θ = H -1 A -1 n,τ θ n -θ -(θ n+1 -θ) ν n+1 + H -1 (∇G(θ n ) -∇ h g(X n+1 , θ n )) -H -1 (∇G(θ n ) -H(θ n -θ)) , which implies that θ n,τ -θ = t ′ n n k=0 ln(k + 1) τ ′ H -1 A -1 k,τ θ k -θ -(θ k+1 -θ) ν k+1 + t ′ n n k=0 ln(k + 1) τ ′ H -1 (∇G(θ k ) -∇ h g(X k+1 , θ k )) -t ′ n n k=0 ln(k + 1) τ ′ H -1 (∇G(θ k ) -H(θ k -θ)) • (R n ) is a sequence of matrices lying in H such that ∥R n ∥ F = o (v n ) a.s where v n = (ln n) a n b ,
with a, b ∈ R;

• For all n ≥ 1 and 1 ≤ k ≤ n, and for all A ∈ H,

β n,k A = n j=k+1 (I q -γ j Γ) A n j=k+1 (I q -γ j Γ) and β n,n A = A,
where Γ ∈ H is symmetric and satisfies 0 < λ min (Γ) ≤ λ max (Γ) < +∞.

Then,

∥M n+1 ∥ 2 F = O γ n v 2 n ln n a.s.
Proof of Lemma 6.3. Let us now consider the events

A n = {∥R n ∥ F > v n or R 2,n > C} B n+1 = {∥R n ∥ F ≤ v n , R 2,n ≤ C, ∥ξ n+1 ∥ F ≤ δ n } C n+1 = {∥R n ∥ F ≤ v n , R 2,n ≤ C, ∥ξ n+1 ∥ F > δ n } with δ n = γ -1/2 n (ln n) -1/2 . One can check that A c n = B n+1 ⊔ C n+1 . Then, one can write M n+1 as M n+1 = n k=1 β n,k γ k R k ξ k+1 1 A k + n k=1 β n,k γ k R k ξ k+1 1 A c k = n k=1 β n,k γ k R k ξ k+1 1 A k + n k=1 β n,k γ k R k ξ k+1 1 B k+1 -E ξ k+1 1 B k+1 |F k + n k=1 β n,k γ k R k ξ k+1 1 C k+1 -E ξ k+1 1 C k+1 |F k .
Let us now give the rates of convergence of these three terms.

Bounding M 1,n+1 := n k=1 β n,k γ k R k ξ k+1 1 A k . There exists a rank n 0 such that for all n ≥ n 0 , ∥I q - γ n Γ∥ op ≤ (1 -λ min γ n ). Furthermore, M 1,n+1 = (I q -γ n Γ)M 1,n (I q -γ n Γ) + γ n R n ξ n+1 1 An . Then, for all n ≥ n 0 , E ∥M 1,n+1 ∥ 2 F |F n ≤ (1 -λ min γ n ) 4 ∥M 1,n ∥ 2 F + γ 2 n ∥R n ∥ 2 F (C + R 2,n ) 1 An . Considering V n+1 = n k=1 (1 + λ min γ k ) 4 ∥M 1,n+1 ∥ 2 F , it follows that E [V n+1 |F n ] ≤ 1 -λ 2 min γ 2 n 4 V n + n k=1 (1 + λ min γ k ) 4 γ 2 n ∥R n ∥ 2 F (C + R 2,n ) 1 An
Moreover, 1 An converges almost surely to 0 implying that

n≥1 n k=1 (1 + λ min γ k ) 4 γ 2 n ∥R n ∥ 2 F (C + R 2,n ) 1 An < +∞ a.s
and applying Robbins-Siegmund Theorem, V n converges almost surely to a finite random variable, i.e

∥M 1,n+1 ∥ 2 F = O n k=1 (1 + λ min γ k ) -4 a.s
and converges exponentially fast.

Bounding M 2,n+1 := n k=1 β n,k γ k R k (ξ k+1 1 B k+1 -E[ξ k+1 1 B k+1 |F k ]). Let us denote Ξ k+1 = R k (ξ k+1 1 B k+1 - E[ξ k+1 1 B k+1 |F k ]). Remark that (Ξ n
) is a sequence of martingale differences adapted to the filtration (F n ). As in [START_REF] Pinelis | Optimum bounds for the distributions of martingales in banach spaces[END_REF] (proofs of Theorems 3.1 and 3.2), let λ > 0 and consider for all t ∈ [0, 1] and j ≤ n,

φ(t) = E cosh λ j-1 k=1 β n,k γ k Ξ k+1 + tβ n,j γ j Ξ j+1 F |F j .
One can check that φ ′ (0) = 0 and (see Pinelis for more details) φ ′′ (t) ≤ λ 2 E ∥β n,j γ j Ξ j+1 ∥ 2 F e λt∥β n,j γ j Ξ j+1 ∥ F cosh λ .

Furthermore, let ϵ j+1 = ξ j+1 1 B j -E[ξ j+1 1 B j |F j ] and note that E[∥ϵ j+1 ∥ 2 F |F j ] ≤ 2C. Then, recalling that δ n = γ -1/2 n (ln n) -1/2 , and since for all k ≥ 2,

E ∥ϵ j+1 ∥ k F |F j ≤ 2 k-2 δ k-2 j E ∥ξ j+1 ∥ 2 F 1 B j |F j ≤ 2 k-1 Cδ k-2 j ,
ans since for any A ∈ H one has ∥β n,k A∥ F ≤ ∥β n,j ∥ op ∥A∥ F , ∥β n,j ∥ 2 op γ 2 j v 2 j exp 2λ ∥β n,j ∥ op v j γ j ln j   .

e j,n ≤ ∞ k=2 λ k ∥β n,j ∥ k op γ k j E ∥Ξ j+1 ∥ k F |F j ≤ ∞ k=2 λ k ∥β n,j ∥ k op γ k j v k j E ∥ϵ j+1 ∥ k F |F k ≤ ∞ k=2 λ k ∥β n,j ∥ k op γ k j v k j 2 k-1 Cδ k-2 j ≤ 2Cλ 2 ∥β n,j ∥ 2 op γ 2 j v 2 j ∞ k=2 (
Take λ = γ -1/2 n v -1 n √ ln n. Let C 0 = ∥β n 0 ,0 ∥ op and remark that for n ≥ 2n 0 (i.e such that γ n/2 λ max (Γ) ≤ 1), and for all j ≤ n/2, ∥β n,j ∥ op ≤ C 0 exp -cλ min (n/2) 1-α , so that for all j ≤ n/2, λ ∥β n,j ∥ op γ j v j a.s -----→ n→+∞ 0.

Furthermore, for all n ≥ 2n 0 , and for all j ≥ n/2, λ ∥β n,j ∥ op √ γ j v j √ ln j ≤ C 0 2 2b+2a+α+1 .

Then, there is a positive constant C ′′ such that for all n ≥ 1 and j ≤ n,

exp λ ∥β n,j ∥ op √ γ j v j ≤ C ′′
ϵ k+1 = ξ k+1 1 C k+1 -E[ξ k+1 1 C k+1 |F k ]
and remark that for n ≥ n 0 ,

E ∥M 3,n+1 ∥ 2 F |F n ≤ (1 -λ min γ n ) 4 ∥M 3,n ∥ 2 F + γ 2 n v 2 n E ∥ϵ n+1 ∥ 2 F |F n ≤ (1 -λ min γ n ) 4 ∥M 3,n ∥ 2 F + γ 2 n v 2 n E ∥ξ n+1 ∥ 2 F 1 ∥ξ n+1 ∥ 2 F ≥γ -1 n |F n Let V ′ n = γ -1 n v -2 n ∥M 3,n ∥ 2 F .
There are a rank n 1 and a positive constant c such that for all n ≥ n 1

E [V n+1 |F n ] ≤ (1 -cγ n ) V n + O γ n E ∥ξ n+1 ∥ 2 F 1 ∥ξ n+1 ∥ 2 F ≥γ -1
n |F n a.s.

Applying Robbins-Siegmund Theorem as well as equation ( 26), it comes

∥M 3,n+1 ∥ 2 F = O γ n v 2 n a.s.
Lemma 6.4. Let J n , K n , r n be sequences of positive random variables such that r n converges almost surely to 0 and J n+1 = (1 -c γ n+1 )J n + γ n+1 r n (J n + K n )

where γ n = c γ n γ with 1/2 < γ < 1 and c γ > 0. In addition, it is assumed that

K n = O(v n ) a.s.
where v n = c v n v (ln n) b with v ∈ R and b ≥ 0. Then

J n = O(v n ) a.s.
Proof of Lemma 6.4. For the sake of simplicity, let us assume that for every n ≥ 0, c γ n+1 ≤ 1 (up to take n large enough). Now, consider the event E n,c = {|r n | ≤ c/2}, and therefore 1 E C n,c converges almost surely to 0. Hence, J n+1 can be rewritten as:

J n+1 ≤ (1 -c γ n+1 ) J n + c 2 γ n+1 (J n + K n ) + =:δn γ n+1 r n (J n + K n ) 1 E C n,c ≤ 1 - c 2 γ n+1 J n + c 2 γ n+1 K n + δ n 1 E C n,c
By induction, one can check that for all n ≥ 0:

J n ≤ βn,0 J 0 + c 2 n-1 k=0 βn,k+1 γ k+1 K k =:J 1,n + n-1 k=0 βn,k+1 δ k 1 E C k,c =:J 2,n
with βn,k := n j=k+1 1 -c 2 γ j and βn,n := 1. Using standard calculations, we can easily show that βn,0 converges at an exponential rate. Furthermore, J 2,n can be written as βn,0

n-1 k=0 β-1 k,0 δ k 1 E C k,c
and since

1 E C n,c
converges almost surely to 0, the sum is almost surely finite, leading to J 2,n = O βn,0 a.s. and this term thus converges at an exponential rate. Finally, there exists a random variable K such that for every n ≥ 1, K n ≤ Kv n almost surely, leading to the induction relation:

Figure 1 :

 1 Figure 1: Evolution of the mean squared error with respect to the sample size for logistic regression.

Figure 2 :

 2 Figure 2: Evolution of the mean squared error with respect to the sample size for geometric median estimation.

Figure 3 :Figure 4 :

 34 Figure 3: Evolution of the mean squared error with respect to the sample size for parameters estimation in a spherical Gaussian distribution.

Figure 5 :

 5 Figure 5: Evolution of the mean squared error with respect to the sample size for p-means estimation.

ν

  n+1 λ min (A n,τ ) = +∞ a.s.It suggests that lim inf n ∥∇G(θ n )∥ = 0 almost surely and lim inf n K n = 0 almost surely. As K n converges almost surely to a random variable, G(θ n ) converges almost surely to G(θ). By local strict convexity,

  t)φ ′′ (t)dt ≤ (1 + e j,n ) cosh λ j-1 k=1 β n,k γ k Ξ k+1 F with e j,n = E[e λ∥β n,j γ j Ξ j+1 ∥ F -1 -λ∥β n,j γ j Ξ j+1 ∥ F |F j ], which is well defined since Ξ j+1 is a.s. finite. Additionally, considering G n+1 = cosh λ ∥ n k=1 β n,k γ k Ξ k+1 ∥ F n j=1 (1 + e j,n ) and G 0 = 1 and since E[G n+1 |F n ] = G n , it comes E[G n+1 ] = 1. For all r > 0, P ∥M 2,n+1 ∥ F ≥ r = P G n+1 ≥

  2λ) k-2 ∥β n,j ∥ k∥β n,j ∥ 2 op γ 2 j v 2 j exp 2λ ∥β n,j ∥ op √ γ j v jThen,P ∥M 2,n+1 ∥ F ≥ r ≤P   2G n+1 ≥ e λr n j=1 1 + 2Cλ 2 ∥β n,j ∥ 2 op γ 2 j v 2 j exp 2λ ∥β n,j ∥ op v j γ j ln j

Finally

  , one can easily check that (see Lemma E.2 in Cardot and Godichon-Baggioni (2017))There is a positive constant C ′′′ such thatP ∥M 2,n+1 ∥ F ≥ r ≤ exp -rv -1 n γ -1/2 n √ ln n + C ′′′ ln n Then , taking r = (2 + C ′′′ )v n √ γ n ln n, it comes P ∥M 2,n+1 ∥ F ≥ 2 + C ′′′ v n γ n ln n ≤ exp (-2 ln n) = 1 n 2and applying Borell Cantelli's lemma,∥M 2,n+1 ∥ F = O v n γ n ln n a.s. Bounding M 3,n+1 := n k=1 β n,k γ k R n (ξ k+1 1 C k+1 -E[ξ k+1 1 C k+1 |F k ]). Let us denote

Table 1 :

 1 . Accuracy of UWASNA, USNA, WASNA, SNA and ASGD algorithms on "COVTYPE" dataset.

		UWASNA USNA WASNA SNA ASGD
	Training Accuracy(%)	75.63	75.31	75.52	75.38 74.54
	Test Accuracy(%)	75.61	75.34	75.50	75.33 74.64

where t ′ n = 1 n k=0 log(k+1) τ ′ . Rate of convergence of t ′ n n k=0 log(k + 1) τ ′ H -1 (∇G(θ k ) -∇ h g(X k+1 , θ k )). Analogous to the proof of Theorem 3.1, one can check with the help of a law of large numbers for martingales that

With the help of Assumption (A3) and since θ n converges almost surely to θ, ∥G(θ n ) -H(θ n -θ)∥ = O ∥θ n -θ∥ 2 a.s. Then, with the help of Thanks to equality (25) and since ν < 1, we can prove with the help of the equation ( 24) that

and this term is negligible as soon as

. One can easily check that

.

With the help of an Abel's transform and since A n,τ converges almost surely to H -1 , one has

Thus, as ∥θ n -θ∥ 2 = o ln n 1+δ n ν a.s. we have for all δ > 0

which is negligible as soon as 2 -ν > 1. Finally, we can conclude that ∥θ n,τ -θ∥ 2 = o ln n 1+δ n a.s.

Useful lemmas

The following lemma is a corollary of the Robbins-Siegmund theorem [START_REF] Robbins | A convergence theorem for non negative almost supermartingales and some applications[END_REF].

Lemma 6.1. Let (V n ), (B n ), (D n ) and (a n ) be positive sequences adapted to F = (F n ). Assume that V 0 is integrable and, for all n ≥ 0,

s. The proof of this lemma is given in Chapter 1.III in Duflo (1990), and here we give a generalized version of it. Lemma 6.2. Let (V n ), (B n ), (D n ), E n and (a n ) be positive sequences adapted to F = (F n ). Assume that V 0 is integrable and, for all n ≥ 0,

s. Proof of Lemma 6.2. Note that the case where E n = 0 is exactly the case of Lemma 6.1. Therefore, we are going to study the case where E n ̸ = 0. We define α n := n k=0 (1 + E k ). Note that since ∞ n=0 E n converges almost surely, α n converges almost surely to a finite random variable α ∞ . Moreover, noting

B n < +∞ a.s.

According to Lemma 6.1, we have V ′ n = o(a n ) a.s., and therefore

We now give two lemmas which will be tools for the study of the rate of convergence associated to the estimates A n . Lemma 6.3. Let us denote by H = M q (R) the set of squared matrices of size q × q. Let us consider

where • (ξ n ) is a H-valued martingale differences sequence adapted to a filtration (F n ) such that

where C is a non-negative random variable and (R 2,n ) n converges almost surely to 0;

• γ n = cn -γ with c > 0 and γ ∈ (1/2, 1);