# Complexity reduction in online stochastic Newton methods with potential O(N d) total cost

Antoine Godichon-Baggioni, Bruno Portier, Guillaume Sallé

<figure><img src="image_1.png" /></figure>

## To cite this version:

Antoine Godichon-Baggioni, Bruno Portier, Guillaume Sallé. Complexity reduction in online stochastic Newton methods with potential O(N d) total cost. 2026. hal-05493562

HAL Id: hal-05493562

https://hal.science/hal-05493562v1

Preprint submitted on 4 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Complexity reduction in online stochastic Newton methods with potential $O(Nd)$ total cost

Antoine Godichon-Baggioni$^{(1)}$, Bruno Portier$^{(2)}$ and Guillaume Sallé$^{(2)}$

(1) Laboratoire de Probabilités, Statistique et Modélisation, Sorbonne Université, 75005 Paris, France, antoine.godichon_baggioni@upmc.fr

(2) Laboratoire de mathématiques de l'INSA, INSA Rouen Normandie, 76800 Saint-Etienne du Rouvray, France

## Abstract

Optimizing smooth convex functions in stochastic settings, where only noisy estimates of gradients and Hessians are available, is a fundamental problem in optimization. While first-order methods possess a low per-iteration cost, their convergence is slow for ill-conditioned problems. Stochastic Newton methods utilize second-order information to correct for local curvature, but the $O(d^3)$ per-iteration cost of computing and inverting a full Hessian, where $d$ is the problem dimension, is prohibitive in high dimensions. This paper introduces an online mini-batch stochastic Newton algorithm. The method employs a random masking strategy that selects a subset of Hessian columns at each iteration, substantially reducing the per-step computational cost. This approach allows the algorithm, in the mini-batch setting, to achieve a total computational cost for a single pass over $N$ data points of $O(Nd)$, which is comparable to first-order methods while retaining the advantages of second-order information. We establish the almost sure convergence and asymptotic efficiency of the resulting estimator. This property is obtained without requiring iterate averaging, which distinguishes this work from prior analyses.

**Keywords:** Convex optimization; Online algorithm; Stochastic Newton method; Mini-batch; Computational complexity; Asymptotic efficiency

## 1 Introduction

Estimating the minimizer $\theta^* \in \mathbb{R}^d$ of a smooth convex function $F: \mathbb{R}^d \to \mathbb{R}$ is a fundamental problem in applied mathematics and machine learning. This paper studies this problem within a stochastic online framework, wherein the function $F$ is not directly accessible. Instead, we assume access to stochastic oracles that provide noisy estimates of the gradient $\nabla F(\theta)$ and Hessian $\nabla^2 F(\theta)$. This setting is common in large-scale learning, where the objective $F$ is defined as an expectation of a known function $f$, $F(\theta) = \mathbb{E}_{\xi}[f(\xi, \theta)]$, with respect to a random variable $\xi$. The minimizer $\theta^*$ is then estimated from sequential samples of $\xi$.
For large-scale problems, particularly when datasets exceed memory capacity or arrive sequentially, online algorithms are essential. First-order methods, notably stochastic gradient descent (SGD), are widely employed due to their simplicity and low per-iteration computational cost, typically $O(d)$. Their theoretical properties are well-understood, with foundational asymptotic results in Pelletier (2000); Duflo (1996) and non-asymptotic analyses in Bach and Moulines (2011). We refer to Bottou et al. (2018) for a comprehensive survey. SGD employs a single scalar step-size, which can lead to slow convergence for ill-conditioned problems. To address this limitation, a broad class of algorithms, termed *conditioned SGD*, incorporates a pre-conditioning matrix $C_{n-1}$ into the update:

$$ \theta_n = \theta_{n-1} - \alpha_n C_{n-1} g_n(\theta_{n-1}), $$

where $\alpha_n > 0$ is the step-size and $g_n(\theta_{n-1})$ is a stochastic estimate of the gradient $\nabla F(\theta_{n-1})$. The choice of $C_{n-1}$ dictates the trade-off between convergence performance and computational cost. Well-known adaptive methods, such as AdaGrad Duchi et al. (2011), RMSProp Hinton (2012), and Adam Kingma and Ba (2017), utilize diagonal matrices for $C_{n-1}$. This strategy adapts the step-size for each coordinate but fails to capture the off-diagonal curvature information.

The ideal conditioning, inspired by Newton's method, sets $C_{n-1}$ to approximate the inverse Hessian, $(\nabla^2 F(\theta_{n-1}))^{-1}$. Methods that explicitly estimate this matrix are known as *stochastic Newton methods*. This pre-conditioning corrects for the local curvature of the objective function, rendering the algorithm robust to ill-conditioned problems. Stochastic Newton methods are capable of achieving *asymptotic efficiency*, a fundamental benchmark for asymptotic optimality introduced in Pelletier (2000); Duflo (1996). This property signifies that the estimator's asymptotic covariance matrix attains the Cramér-Rao lower bound in relevant statistical settings; see, e.g., Leluc and Portier (2023). Asymptotic efficiency can be achieved if the pre-conditioning matrix $C_n$ converges almost surely to the inverse Hessian at the optimum, $H^{-1} := (\nabla^2 F(\theta^*))^{-1}$, as $n \to \infty$ Leluc and Portier (2023). This result relaxes the stronger convergence rate conditions required by previous analyses Boyer and Godichon-Baggioni (2023).

The principal obstacle to their practical application, however, remains the computational cost. A straightforward implementation, such as that in Leluc and Portier (2023), requires computing and inverting an empirical Hessian at each iteration. This is an $O(d^3)$ operation that is prohibitive in high-dimensional settings ($d \gg 1$).

**Related works.** Significant research has focused on reducing the $O(d^3)$ computational burden of stochastic Newton methods. Stochastic Quasi-Newton methods, such as (L)-BFGS variants for the online setting Schraudolph et al. (2007); Moritz et al. (2016), circumvent direct Hessian calculations by using gradient differences to build an approximation of the inverse Hessian, typically at a cost of $O(md)$. Other approaches avoid forming or inverting the full Hessian by using iterative linear solvers that only require Hessian-vector products, such as stochastic Newton-CG methods Byrd et al. (2011).

While these methods reduce the per-iteration cost, they generally lack guarantees of asymptotic efficiency, as they do not ensure that the pre-conditioning matrix converges
to the true inverse Hessian.

In the particular case where the Hessian estimates are rank-one matrices, as in generalized linear models, algorithms with a reduced $O(d^2)$ cost have been developed.

For instance, Bercu et al. (2020) introduced an $O(d^2)$ algorithm that leverages the Sherman-Morrison formula for efficient rank-1 updates Dennis and Moré (1977). Boyer and Godichon-Baggioni (2023) extended this approach to a broader class of models, achieving asymptotic efficiency with $O(d^2)$ complexity per iteration.

Recently, an online algorithm was introduced by Godichon-Baggioni et al. (2025), which estimates the inverse Hessian recursively with a $O(d^2)$ cost per iteration, applies to general Hessian structures and possesses established convergence guarantees. However, this algorithm required an additional averaging step to achieve asymptotic efficiency.

Randomized iterative algorithms for matrix inversion have been proposed to reduce computational cost. These methods, which include stochastic variants of quasi-Newton updates, compute an approximation of a fixed matrix inverse by repeatedly solving a "sketched," or projected, version of the linear matrix equation Gower and Richtárik (2017). Such algorithms necessarily operate in an offline setting, where the target matrix is known and available at each iteration. This framework is distinct from the online estimation setting considered here, where the target Hessian is unknown and must be estimated sequentially from a stream of noisy data. Consequently, methods designed for offline iteration are not directly applicable to our problem.

An alternative path to achieving asymptotic efficiency, distinct from second-order conditioning, is the averaging method. Introduced independently by Polyak (1990) and Ruppert (1988), and further developed by Polyak and Juditsky (1992), this method, commonly known as Polyak-Ruppert averaging, allows first-order SGD variants to attain the optimal asymptotic covariance. While powerful, averaging does not inherently adapt to the local geometry in the same manner as methods employing Hessian information.

Mini-batching is another important consideration for large-scale applications, as it can stabilize estimates, leverage parallel architectures, and reduce update overhead Bottou et al. (2018). In the particular case where the Hessian estimates are rank-one matrices, Godichon-Baggioni and Werge (2025) proposed an asymptotically efficient mini-batch stochastic Newton algorithm with reduced cost by updating the pre-conditioning matrix less frequently. However, this approach is not directly applicable to general Hessian structures and incorporates less information from the Hessian estimates, as it utilizes only a single sample from the mini-batch for the pre-conditioner update rather than aggregating information from the entire batch. For the general case, the algorithm of Godichon-Baggioni et al. (2025) is asymptotically efficient in the purely online setting (batch size $b = 1$), but its extension to the mini-batch setting ($b > 1$) case remains a challenge.

**Contributions.** This paper introduces an asymptotically efficient online mini-batch (or streaming) algorithm, the *masked Stochastic Newton Algorithm* (mSNA), which extends and improves upon the method proposed in Godichon-Baggioni et al. (2025). Our work presents two principal advancements. First, we establish the almost sure convergence of
the mSNA estimator and its asymptotic efficiency *without* requiring iterate averaging, theoretical guarantees not provided in the prior work. Second, we develop a practical method with reduced computational cost for processing mini-batched Hessian information. Our algorithm operates in a fully online setting (unlike offline iterative solvers) and applies to general Hessian structures (rather than only rank-one). Its computational gains are achieved by deriving the inverse Hessian estimator itself from a stochastic gradient descent (SGD) on a specific matrix functional. We then employ a random masking strategy for this underlying SGD, which selects a subset of Hessian columns at each iteration; this technique can be seen as an adaptation of SGD with coordinate sampling Leluc and Portier (2022). This approach allows the algorithm to process mini-batches of size *b* with a per-iteration computational cost of $O(\ell bd + \ell d^2)$, where $\ell \ll d$ is a rank parameter. On a dataset of size *N*, the algorithm performs $N/b$ iterations. By setting the batch size $b=d$ and $\ell=1$, the total complexity for one pass over the data is of order $O(Nd)$, matching the complexity scaling of first-order methods while retaining the benefits of second-order information. We also provide a study of an averaged version of the algorithm, which can improve practical performance at a minimal additional computational cost.

The objective of our work is therefore to propose: (i) online stochastic Newton estimators, (ii) developed within a setting as general as possible, (iii) with reduced computational complexity and mini-batch processing (possibly $O(Nd)$ operations for one pass over the data), and (iv) with strong theoretical guarantees, in particular asymptotic efficiency.

To the best of our knowledge, no previous work has successfully combined these four aspects simultaneously.

**Paper Organization.** Section 2 introduces the notations, the optimization framework, and the underlying assumptions. Section 3 presents the algorithm for recursively estimating the inverse Hessian. This estimator is then incorporated into the proposed mini-batch stochastic Newton algorithm in Section 4, where we establish its theoretical convergence analysis. Section 5 details the implementation of this algorithm in the streaming (mini-batch) setting and analyzes its computational complexity. Finally, Section 6 provides numerical experiments to illustrate the algorithm’s performance. Proofs are deferred to the Appendix.

## 2 Framework

**Notations.** The gradient and Hessian operators are denoted by $\nabla$ and $\nabla^2$, respectively. We denote by $\mathcal{M}_d(\mathbb{R})$ the set of square matrices of order $d$ with real coefficients, and by $S_d^{++}(\mathbb{R})$ the set of symmetric positive definite matrices of order $d$. Let $I_d$ be the identity matrix of order $d$. On $\mathcal{M}_d(\mathbb{R})$, we denote by $\|\cdot\|_F$ (resp. $\|\cdot\|_{op}$) the Frobenius norm (resp. the operator norm), and by $\langle\cdot, \cdot\rangle_F$ the Frobenius inner product. For a vector $v \in \mathbb{R}^d$, we denote by $\|v\|$ its Euclidean norm. For a matrix $A \in \mathcal{M}_d(\mathbb{R})$, $A^T$ denotes its transpose and $\text{Tr}(A)$ its trace. We denote the convergence in distribution by $\xrightarrow{\mathcal{L}}$.

Instead of studying general step-sizes $(\alpha_n)_{n \ge 1}$ satisfying the Robbins-Monro conditions $\sum_n \alpha_n = +\infty$ and $\sum_n \alpha_n^2 < +\infty$, we restrict ourselves for simplification to step sizes of
the form $\alpha_n := c_\alpha/(n^\alpha + c'_\alpha)$, with $c_\alpha > 0, c'_\alpha \ge 0, 1/2 < \alpha \le 1$, and following Carpentier et al. (2015) we call such a sequence a $\sigma(\alpha, c_\alpha, c'_\alpha)$-sequence.

**Problem setting** Let $F: \mathbb{R}^d \to \mathbb{R}$ be a convex and twice differentiable function. Assuming its existence and uniqueness, we consider the problem of estimating the minimizer $\theta^*$ of $F$, that is,

$$ \theta^* := \arg \min_{\theta \in \mathbb{R}^d} F(\theta), $$

with a sequence of estimators $(\theta_n)_{n \ge 0}$. At each iteration $n \ge 1$, we assume access to stochastic oracles of the gradient and the Hessian. More precisely, given a point $\theta \in \mathbb{R}^d$, the stochastic oracles return random estimates $g_n(\theta)$ (resp. $h_n(\theta)$) of the gradient $\nabla F(\theta)$ (resp. the Hessian $\nabla^2 F(\theta)$). In addition, we denote by $H := \nabla^2 F(\theta^*)$ the Hessian at the minimizer. Let $(\mathcal{F}_n)_{n \ge 0}$ be the natural filtration generated by the iterates $(\theta_n)_{n \ge 0}$ and the oracles $(g_n, h_n)_{n \ge 1}$.

(A1) (Unbiased Estimates). For all $n \ge 1$ and $\theta \in \mathbb{R}^d$:

$$ \text{a)} \quad \mathbb{E}[g_n(\theta) | \mathcal{F}_{n-1}] = \nabla F(\theta) \quad \text{a.s.} $$

$$ \text{b)} \quad \mathbb{E}[h_n(\theta) | \mathcal{F}_{n-1}] = \nabla^2 F(\theta) \quad \text{a.s.} $$

(A2) (Growth Condition). There exists $\mathcal{L}_g, \mathcal{L}_h, \sigma^2 \ge 0$ such that for all $n \ge 1$ and $\theta \in \mathbb{R}^d$:

$$ \text{a)} \quad \mathbb{E} \left[ \|g_n(\theta)\|^2 | \mathcal{F}_{n-1} \right] \le \mathcal{L}_g (F(\theta) - F(\theta^*)) + \sigma^2 \quad \text{a.s.} $$

$$ \text{b)} \quad \mathbb{E} \left[ \|h_n(\theta)\|_{op}^2 | \mathcal{F}_{n-1} \right] \le \mathcal{L}_h \quad \text{a.s.} $$

(A3) (Hessian at Minimizer). The Hessian matrix $H := \nabla^2 F(\theta^*)$ at the minimizer is positive definite, and the mapping $\theta \mapsto \nabla^2 F(\theta)$ is continuous at $\theta^*$.

(A4) (Lyapunov Conditions). There exist $q, q' > 2$ and $M, M' > 0$ such that:

$$ \text{a)} \quad \sup_{\substack{n \ge 1 \\ \|\theta - \theta^*\| < M}} \mathbb{E} \left[ \|g_n(\theta) - \nabla F(\theta)\|^q | \mathcal{F}_{n-1} \right] < +\infty \quad \text{a.s.} $$

$$ \text{b)} \quad \sup_{\substack{n \ge 1 \\ \|\theta - \theta^*\| < M'}} \mathbb{E} \left[ \|h_n(\theta) - \nabla^2 F(\theta)\|_{F}^{q'} | \mathcal{F}_{n-1} \right] < +\infty \quad \text{a.s.} $$

(A5) (Covariance Limit). If $(\theta_n)_{n \ge 0}$ converges almost surely to $\theta^*$, there exists $\Gamma \in S_d^{++}$ such that $\text{Cov}(g_n(\theta_{n-1}) | \mathcal{F}_{n-1}) \xrightarrow{n \to \infty} \Gamma$.

These assumptions are standard in stochastic approximation and follow those in Pelletier (1998b) and Pelletier (1998a) (see also Leluc and Portier (2023)). Assumption (A1) could be relaxed to accommodate a vanishing bias, but for simplicity we do not consider this case; we refer to Surendran et al. (2024) for an analysis thereof. Assumption (A2) a) is crucial for obtaining the strong consistency of stochastic gradient estimators. It is related to the *expected smoothness* introduced by Gower et al. (2019) and, as in
Leluc and Portier (2023), permits changes in the distribution of the gradient oracles. Assumption (A2) b) is more restrictive and implies in particular that the Hessian of F is uniformly bounded. Assumption (A3) implies that the function F is locally strongly convex around $\theta^*$ and justifies the use of Newton methods. This assumption is crucial for attaining asymptotic efficiency with any gradient-based method (see, e.g., Pelletier (2000); Godichon-Baggioni (2019) for averaged stochastic gradient algorithms and Leluc and Portier (2023) for conditioned gradient algorithms). Assumption (A4) is a standard known Lyapunov condition, required to hold only locally around the minimizer $\theta^*$. Finally, Assumption (A5) is a standard condition required to establish the asymptotic normality of gradient-based estimators (Pelletier (1998b)).

**In practice** This framework, which makes assumptions on generic oracles, is more general than the standard setting where F is defined as the expectation of a known function f, i.e., $F(\theta) = \mathbb{E}_{\xi}[f(\xi, \theta)]$. This generality permits the oracle distributions to change over time (a setting also considered in Leluc and Portier (2023)). Furthermore, this oracle-based framework facilitates the study of mini-batching strategies, including those with varying batch-sizes or non-uniform sampling.

In the standard setting where $F(\theta) = \mathbb{E}_{\xi}[f(\xi, \theta)]$ and i.i.d. samples $(\xi_n)_{n \ge 1}$ are available, the natural oracles are

$$g_n(\theta) = \nabla f(\xi_n, \theta) \quad \text{and} \quad h_n(\theta) = \nabla^2 f(\xi_n, \theta).$$

In the online mini-batch (streaming) setting, data are processed in blocks of size $b > 1$. Letting $\{\xi_{n,i}\}_{i=1}^b$ denote the samples in the block at iteration $n$, the oracles are the empirical averages:

$$g_n(\theta) = \frac{1}{b} \sum_{i=1}^{b} \nabla f(\xi_{n,i}, \theta) \quad \text{and} \quad h_n(\theta) = \frac{1}{b} \sum_{i=1}^{b} \nabla^2 f(\xi_{n,i}, \theta).$$

The analysis of this mini-batch setting is detailed in Section 5.

These oracle assumptions are satisfied if corresponding assumptions hold for the function f. These assumptions are given in Appendix C. For classical problems such as logistic regression and p-means, a verification of the assumptions can be found in the appendix A of Godichon-Baggioni et al. (2025).

## 3 Online Estimation of the Inverse of a Positive Definite Matrix

In this section, we aim to estimate recursively the inverse of the matrix $H := \nabla^2 F(\theta^*)$ assumed to be positive definite, with the help of a sequence of estimates $(H_n)_{n \ge 1}$ of H adapted to the filtration $(\mathcal{F}_n)_{n \ge 0}$. Within our framework, these estimates are obtained by evaluating the Hessian oracle $h_n$ at an estimator $\hat{\theta}_{n-1}$ that converges to $\theta^*$. This sequence $(\hat{\theta}_n)_{n \ge 0}$ may represent the primary iterates $(\theta_n)_{n \ge 0}$ themselves, or an averaged version $(\hat{\theta}_n)_{n \ge 0}$ derived from them.
## 3.1 Stochastic Gradient Estimation of the Inverse

We first observe that $H^{-1}$ is the unique minimizer of the quadratic functional $J: \mathcal{M}_d(\mathbb{R}) \to \mathbb{R}_+$ defined for all $A \in \mathcal{M}_d(\mathbb{R})$ by:

$$J(A) = \|H^{1/2}(A - H^{-1})\|_F^2,$$

where $H^{1/2}$ denotes the unique symmetric positive definite square root of the matrix $H$. The function $J$ is twice differentiable, and standard matrix calculus (see Appendix D) yields the gradient:

$$\nabla J(A) = 2(H^{1/2})^T H^{1/2}(A - H^{-1}) = 2(HA - I_d).$$

Thus, the functional $J$ is $2\lambda_{\min}(H)$-strongly convex and $2\lambda_{\max}(H)$-smooth.

The main interest of the functional $J$ is that its gradient only depends on $H$ through the product $HA$. One could then recursively estimate $H^{-1}$ using a SGD on the function $J$, by replacing $H$ by its estimates $H_n$. The algorithm we will introduce in the next section is derived from a symmetric, positive definite factorization. To aid in its analysis, we also define the related symmetrized functional:

$$J_{sym}(A) := \frac{1}{2} (J(A) + J(A^T)),$$

which has gradient

$$\nabla J_{sym}(A) = \frac{1}{2} (\nabla J(A) + (\nabla J(A^T))^T) = HA + AH - 2I_d,$$

which is symmetric when $A$ is symmetric. A naive SGD estimator $(A_n)_{n \ge 0}$ based on either $J$ or $J_{sym}$ faces two significant limitations. First, a primary drawback of such an estimator is the computational cost. Indeed, at each step, the dense matrix multiplication requires $O(d^3)$ operations using standard algorithms. To reduce this cost, we will consider a sketched version of the gradient. A second limitation is that these updates do not guarantee the quadratic form associated to $A_n$ remains positive definite, i.e., $v^T A_n v > 0$ for all $v \in \mathbb{R}^d \setminus \{0\}$, which is crucial to ensure descent directions when $A_n$ is used as a preconditioner. We will thus propose a modification of the SGD update to ensure the positive definiteness of the estimator.

## 3.2 A Reduced-Cost Positive Estimator of the Inverse

Let $\ell \in \{1, \dots, d\}$ be an integer. Let $(M_{\ell,n})_{n \ge 1}$ be a sequence of random diagonal projection matrices, adapted to the filtration $(\mathcal{F}_n)_{n \ge 0}$. At each step $n \ge 1$, $M_{\ell,n}$ is constructed by selecting a subset $I_n \subset \{1, \dots, d\}$ of size $\ell$ uniformly at random, and setting $(M_{\ell,n})_{ii} = 1$ if $i \in I_n$ and $0$ otherwise. We assume this random selection at step $n$ is independent of $\mathcal{F}_{n-1}$ and of the stochastic Hessian $H_n$. It follows directly that $\mathbb{E}[M_{\ell,n} | \mathcal{F}_{n-1}] = \ell I_d$. Let $\tilde{H}_n := M_{\ell,n}H_n$ be the sketched random Hessian. Starting from
a positive definite initial estimator $A_0 := I_d$, we propose to estimate $H^{-1}$ by a sequence of estimators $(A_n)_{n \ge 0}$ defined recursively for all $n \ge 1$ by:

$$A_n := A_{n-1} - \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \left( \gamma_n \left( \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{\ell,n} \right) - \gamma_n^2 \tilde{H}_n A_{n-1} \tilde{H}_n^T \right), \quad (1)$$

where $(\gamma_n)_{n \ge 1}$ is a $\sigma(\gamma, c_\gamma, c'_\gamma)$-sequence with $\gamma \in (\frac{1}{2}, 1)$.

The update (1) warrants several remarks. First, the estimator $A_n$ is symmetric and positive definite by construction. As $A_0$ is positive definite, the sequence $(A_n)_{n \ge 0}$ remains positive definite by induction, as shown by the following factorization:

$$A_n = \left( I_d - \gamma_n \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \tilde{H}_n \right) A_{n-1} \left( I_d - \gamma_n \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \tilde{H}_n \right)^T + 2\gamma_n \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} M_{\ell,n}. \quad (2)$$

Indeed, $(I_d - \gamma_n \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \tilde{H}_n)$ is invertible, so by Sylvester's law of inertia, the first term is positive definite if $A_{n-1}$ is positive definite, and the second term is positive semi-definite since $M_{\ell,n}$ is an orthogonal projection matrix. Second, the $O(\gamma_n)$ term in the update, $(\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{\ell,n})$, is a symmetrisation of the projected stochastic gradient $2M_{\ell,n} (H_n A_{n-1} - I_d)$ of $J$ at point $A_{n-1}$, using the fact that $A_{n-1}$ is symmetric. Its expectation conditionally to $F_{n-1}$ is $\frac{\ell}{d} \nabla J_{sym}(A_{n-1})$. The sketching introduces an expected scaling factor of $\ell/d$, which slows down the convergence of the estimator to a neighborhood of $H^{-1}$. Third, the $O(\gamma_n^2)$ term in the update, $-\gamma_n^2 \tilde{H}_n A_{n-1} \tilde{H}_n^T$, is precisely the correction required to complete the quadratic form in the factorization above, which is what guarantees positive definiteness. Finally, the indicator $\mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}}$ is a standard truncation to ensure the stability of the recursion by bounding the effect of large stochastic samples.

### 3.3 Convergence Results

We now state the convergence properties of the estimator $(A_n)_{n \ge 0}$.

**Proposition 3.1.** Let $(\hat{\theta}_n)_{n \ge 0}$ be a sequence of estimators of $\theta^*$ adapted to the filtration $(F_n)_{n \ge 0}$, and let $(A_n)_{n \ge 0}$ be the sequence defined by (1) using $H_n := h_n(\hat{\theta}_{n-1})$. Suppose Assumptions (A1), (A2), b), (A3), and (A4)b) hold with $q' > \frac{2}{\gamma}$.

If $\hat{\theta}_n \to \theta^*$ almost surely, then $A_n \to H^{-1}$ almost surely.

In that case, let $\delta_n := \sup_{k \ge n} \|\hat{\theta}_{k-1} - \theta^*\|$. Let $\mu : \mathbb{R}_+ \to [0; +\infty]$ be the local modulus of continuity of $\nabla^2 F$ at $\theta^*$ defined by $\mu(r) := \sup_{\|\theta - \theta^*\| \le r} \|\nabla^2 F(\theta) - \nabla^2 F(\theta^*)\|_F$. By Assumption (A3), the mapping $\mu$ is finite in a neighborhood of 0 and satisfies $\lim_{r \to 0} \mu(r) = 0$. If $\gamma < 1$, then for any $\eta > 0$:

$$\|A_n - H^{-1}\|_F^2 = O\left(\gamma_n(\ln n)^{1+\eta} + \mu\left(\delta_{\left[\frac{n}{2}\right]}\right)^2\right) \quad \text{a.s.}$$
Proposition 3.1 establishes that the estimator sequence $(A_n)_{n \ge 0}$ is strongly consistent, provided that $\hat{\theta}_n \to \theta^*$ a.s. Furthermore, it quantifies the asymptotic error by decomposing it into two components: a variance term, $O(\gamma_n(\ln n)^{1+\eta})$, which is inherent to the stochastic approximation, and a bias term, $O(\mu(\delta_{\lceil \frac{n}{2} \rceil})^2)$, which is controlled by the convergence rate of the estimator $\hat{\theta}_n$. This sequence $(A_n)_{n \ge 0}$ will serve as the conditioning matrix in the stochastic Newton algorithm developed in the next section.

## 3.4 Relation to Previous Works

In Godichon-Baggioni et al. (2025), the authors proposed a similar online algorithm to estimate $H^{-1}$ based on a Robbins-Monro procedure to find a zero of the mapping $A \mapsto HA + AH - 2I_d$, which is the gradient of the functional $J_{sym}$ defined above, and incorporating sketching. Their algorithm ensures the positive definiteness of the estimator at each step by introducing a more strict truncation than ours, and a projection onto a ball of matrices with slowly growing radius.

In the offline setting, Agarwal et al. (2017) proposed an algorithm to estimate the inverse of a positive definite matrix based on a Taylor expansion of the inverse around a known matrix. The recursion they propose is exactly a SGD on the quadratic functional $J$ defined above with fixed step-size $\gamma_n = 1$, without sketching and without the positive definite correction term, and after a few recursion they compute only the product of this inverse hessian estimate with the stochastic gradient, resulting in a $O(Sd^2)$ cost per iteration, with $S$ the number of recursions. The inverse hessian estimate starts from scratch at each iteration of the main algorithm.

## 3.5 Remark on a Weighted Averaged Version

We conclude this section with a remark on weighted averaging. To obtain a faster convergence rate, one can also consider an averaged estimator (see, e.g. Polyak and Juditsky (1992), Pelletier (2000)). Since standard averaging can be sensitive to initialization, we consider a weighted averaged version (see Mokkadem and Pelletier (2011) or Boyer and Godichon-Baggioni (2023)) recursively defined by $\bar{A}_0 := A_0$ and, for all $n \ge 1$:

$$ \bar{A}_n = \left( 1 - \frac{\omega_n}{\sum_{k=0}^{n} \omega_k} \right) \bar{A}_{n-1} + \frac{\omega_n}{\sum_{k=0}^{n} \omega_k} A_n, \quad \text{where} \quad \omega_n = (\ln(n+1))^\tau, \quad (3) $$

where $\tau \ge 0$ is a weighting parameter; $\tau = 2$ is a common choice. As $\bar{A}_n$ is a convex combination of the positive definite estimators $(A_k)_{0 \le k \le n}$, it remains positive definite.

Despite the known benefits of averaging for optimizing asymptotic variance, the estimator $(\bar{A}_n)_{n \ge 0}$ will not be employed in the sequel. The primary challenge for the estimator $(A_n)_{n \ge 0}$ is not the stochastic noise around the optimum, but rather the slow convergence during the transient phase since the algorithm modifies only $\ell$ rows and columns at each iteration. Averaging introduces inertia, which would further impede this initial convergence. Furthermore, the averaging update (3) is a dense operation, requiring $O(d^2)$ memory access to read and write all components of the estimator. This dense
access requirements contrasts with the non-averaged update (1); while the calculation of that update term requires $O(\ell d^2)$ operations, the update itself only modifies the $O(\ell d)$ components of $A_{n-1}$ corresponding to the selected $\ell$ rows and columns. Given that the averaging step introduces both undesirable inertia and a dense memory access requirement, we will retain the non-averaged estimator $(A_n)_{n \ge 0}$ for the stochastic Newton algorithm in the next section. For a convergence analysis of a related weighted-averaging scheme for estimating the inverse of a positive definite matrix, one can refer to Godichon-Baggioni et al. (2025).

# 4 Algorithms

## 4.1 A New Stochastic Newton Algorithm

We propose here an algorithm named "masked Stochastic Newton algorithm" (mSNA) defined recursively for $n \ge 1$ by

$$ \theta_n := \theta_{n-1} - \alpha_n (A_{n-1} + \nu_n I_d) g_n(\theta_{n-1}) \quad (4) $$

$$ \tilde{H}_n := M_{\ell,n} h_n(\theta_{n-1}) \quad (5) $$

$$ A_n := A_{n-1} - \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \left( \gamma_n \left( \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{\ell,n} \right) - \gamma_n^2 \tilde{H}_n A_{n-1} \tilde{H}_n^T \right), \quad (6) $$

where $\theta_0 \in \mathbb{R}^d$ and $A_0 \in \mathcal{M}_d(\mathbb{R})$ are chosen arbitrarily, $(M_{\ell,n})_{n \ge 1}$ is the sequence of random projection matrices defined in Section 3, $(\alpha_n)_{n \ge 1}$ is a $\sigma(\alpha, c_\alpha, c'_\alpha)$-sequence, $(\gamma_n)_{n \ge 1}$ is a $\sigma(\gamma, c_\gamma, c'_\gamma)$-sequence with $\gamma \in (\frac{1}{2}, 1)$, and $(\nu_n)_{n \ge 1}$ decreases to 0 and satisfies $\sum_{n \ge 1} \alpha_n \nu_n = +\infty$. For that purpose, with any $\nu > 0$, we take $\nu_n = \frac{\nu}{n^{1-\alpha}}$ if $\alpha < 1$ and $\nu_n = \frac{\nu}{\ln n}$ if $\alpha = 1$. The term $\nu_n$ is crucial theoretically to lower bound the smallest eigenvalue of the conditioning matrix, but in practice $\nu$ can be chosen arbitrarily small. The following theorem gives the asymptotic properties of this stochastic Newton algorithm.

**Theorem 4.1.** Under Assumptions (A1) and (A2), if $\alpha + \gamma > \frac{3}{2}$, the estimator $(\theta_n)_{n \ge 0}$ defined by (4) converge almost surely to $\theta^*$. Assume also (A3), (A4) with $q' > \frac{2}{\gamma}$, and (A5). Then, denoting $\zeta := \mathbf{1}_{\alpha=1}/2c_\alpha$, if $1 - \zeta > 0$ it holds that

$$ A_n \xrightarrow[n \to +\infty]{a.s.} H^{-1} \quad \text{and} \quad \frac{1}{\sqrt{\alpha_n}} (\theta_n - \theta^*) \xrightarrow[n \to +\infty]{\mathcal{L}} \mathcal{N} \left( 0, \frac{1}{2(1-\zeta)} H^{-1} \Gamma H^{-1} \right). $$

The proof is given in Appendix, and relies on a recent result of Leluc and Portier (2023) that ensures the asymptotic normality of conditioned SGD as soon as the conditioning matrix is strongly consistent. The results of Theorem 4.1 present a main difference with Godichon-Baggioni et al. (2025) since we are able to obtain the convergence and asymptotic efficiency of the algorithm without an averaging step on $\theta_n$. This is given by the following corollary:

**Corollary 4.2.** Under the assumptions of Theorem 4.1, if $\alpha_n = 1/(n + c'_\alpha)$, then

$$ \sqrt{n} (\theta_n - \theta^*) \xrightarrow[n \to +\infty]{\mathcal{L}} \mathcal{N} (0, H^{-1} \Gamma H^{-1}). $$
In order to establish almost sure asymptotic rates of convergence using the results of Boyer and Godichon-Baggioni (2023) and Godichon-Baggioni and Werge (2025), we have to make a stronger assumption than (A3).

(A6) (Lipschitz at $\theta^*$). The mapping $\theta \mapsto \nabla^2 F(\theta)$ is Lipschitz continuous at $\theta^*$, i.e. there exists $r > 0$ and $L$ such that for all $\theta \in \mathcal{B}(\theta^*, r)$, $\|\nabla^2 F(\theta) - H\|_F \le L \|\theta - \theta^*\|$.

Observe that, as far as we know, averaged conditioned algorithms have only been studied in the literature with the additional Assumption (A6) that the Hessian is Lipschitz continuous at the minimizer, contrary to the non-averaged mSNA with step size $\alpha_n = 1/(n + c'_\alpha)$.

**Theorem 4.3.** Under the assumptions of Theorem 4.1 and Assumption (A6), if the exponent $q$ in Assumption (A4) satisfies $q > \frac{2}{\alpha}$, then for any $\eta > 0$:

$$ \|\theta_n - \theta^*\|^2 = \mathcal{O}\left(\frac{\ln n}{n^\alpha}\right) \quad \text{a.s. and} \quad \|A_n - H^{-1}\|_F^2 = \mathcal{O}\left(\frac{(\ln n)^{1+\eta}}{n^{\min\{\alpha, \gamma\}}}\right) \quad \text{a.s.} $$

## 4.2 A Weighted Averaged Version

We have seen that asymptotic efficiency can only be directly achieved by choosing $\alpha = 1$ and $c_\alpha = 1$, which corresponds to using steps of size $1/n$. However, this leads to very small updates, which can be problematic in the case of poor initialization (see Cénac et al. (2025); Boyer and Godichon-Baggioni (2023)). To address this issue, one can choose $\alpha < 1$ to allow larger step sizes, and then apply an averaging step to recover asymptotic efficiency. This averaging step can be applied to a wide range of conditioned SGD algorithms, including the mSNA algorithm. For any estimator $(\theta_n)_{n \ge 0}$, we define the averaged estimator $(\bar{\theta}_n)_{n \ge 0}$ as a weighted average of the iterates $\theta_k$ for $k \le n$, with weights of the form $\omega_k := (\ln(k+1))^\tau$ for a chosen $\tau \ge 0$ Boyer and Godichon-Baggioni (2023):

$$ \bar{\theta}_0 := \theta_0, \quad \text{and for } n \ge 1, \quad \bar{\theta}_n := \frac{\sum_{k=0}^{n} \omega_k \theta_k}{\sum_{k=0}^{n} \omega_k}. \qquad (7) $$

With the convention $0^0 = 1$, taking $\tau = 0$ corresponds to the simple average of the iterates, while taking $\tau > 0$ enables to give more weights to the last estimates $\theta_n$. This averaging step can be computed in a recursive manner:

$$ \bar{\theta}_n = \left(1 - \frac{\omega_n}{\sum_{k=0}^n \omega_k}\right) \bar{\theta}_{n-1} + \frac{\omega_n}{\sum_{k=0}^n \omega_k} \theta_n. $$

**Remark 4.4.** Observe that since the averaged estimator $\bar{\theta}_n$ is expected to converge faster than the iterates $\theta_n$, it is possible to query the Hessian oracle in $\bar{\theta}_n$ instead of $\theta_n$, which modifies equation (5) into

$$ \tilde{H}_n := M_{\ell,n} h_n(\bar{\theta}_{n-1}). \qquad (8) $$
However, the computations done to obtain an Hessian-vector product often require the gradient as intermediate step, and for the update of $\theta_n$ the stochastic gradient must be computed at $\theta_n$. Therefore, querying the stochastic oracle of the Hessian at $\bar{\theta}_n$ instead of $\theta_n$ can even be less efficient.

The following Theorem gives the strong consistency of the averaged estimates (with either update (5) or (8)), their almost sure rate of convergence as well as their asymptotic efficiency.

**Theorem 4.5.** Under Assumptions (A1) and (A2), if $\alpha + \gamma > \frac{3}{2}$, the estimates $\theta_n$ and $\bar{\theta}_n$ defined by (4), (5) (or (8)), (6) and (7) satisfy

$$ \theta_n \xrightarrow{n \to +\infty} \theta^* \quad \text{and} \quad \bar{\theta}_n \xrightarrow{n \to +\infty} \bar{\theta}^*. $$

Assume also (A3), (A4) with $q > \frac{2}{\alpha}$ and $q' > \frac{2}{\gamma}$, (A6), and $\alpha < 1$. Then for any $\eta > 0$:

$$ \|\theta_n - \theta^*\|^2 = \mathcal{O}\left(\frac{\ln n}{n^\alpha}\right) \quad \text{and} \quad \|\bar{\theta}_n - \theta^*\|^2 = \mathcal{O}\left(\frac{\ln n}{n}\right) \quad \text{a.s.,} $$

$$ \|A_n - H^{-1}\|_F^2 = \mathcal{O}\left(\frac{(\ln n)^{1+\eta}}{n^{\min\{\gamma, \alpha\}}}\right) \quad \text{a.s. and} \quad \|\bar{A}_n - H^{-1}\|_F^2 = \mathcal{O}\left(\frac{(\ln n)^{1+\eta}}{n}\right) \quad \text{a.s.} $$

Moreover, with Assumption (A5), we have

$$ \sqrt{n}(\bar{\theta}_n - \theta^*) \xrightarrow{n \to +\infty} \mathcal{N}(0, H^{-1}\Gamma H^{-1}). $$

Notably, if $\alpha < 1$, the averaged estimator attains the asymptotic efficiency for any $c_\alpha > 0$. It is then possible to slow down the decrease of the step-size $\alpha_n$, by taking a larger $c_\alpha$ balanced by a larger $c'_\alpha$, and still have an asymptotically efficient algorithm.

# 5 Streaming Newton algorithms

## 5.1 The Algorithm

In this section, we focus on a streaming version of the mSNA algorithm, with potentially $\mathcal{O}(Nd)$ total operations for one pass over the data. To this end, we consider the case where the objective function is the expectation of a known function $f$ twice differentiable in its second argument $\theta$:

$$ F(\theta) = \mathbb{E}[f(\xi, \theta)], $$

where $\xi$ is a random variable with values in a measured space $\mathcal{X}$. Following the idea presented in (Godichon-Baggioni and Werge, 2025), we assume from now that at each iteration $n$, we have access to $b$ new i.i.d copies of $\xi$ arriving in a block $\{\xi_{n,1}, \dots, \xi_{n,b}\}$. The streaming mSNA algorithm (and its weighted averaged version) are then defined by using the following stochastic oracles within the updates (4) and (6):

$$ g_n(\theta) = \frac{1}{b} \sum_{i=1}^{b} \nabla f(\xi_{n,i}, \theta) \quad \text{and} \quad h_n(\theta) = \frac{1}{b} \sum_{i=1}^{b} \nabla^2 f(\xi_{n,i}, \theta). \qquad (9) $$
Considering $N = nb$ as the total number of samples processed after $n$ iterations, the number of iterations required to process $N$ samples is $N/b$. This is a factor $b$ fewer than in the purely online ($b = 1$) setting, leading to a reduction in total computational complexity, as discussed below. Although the number of iterations is reduced by a factor of $b$, the covariance of the mini-batch oracles (9) is also reduced by a factor of $b$. Consequently, the limiting covariance matrix $\Gamma$ from Assumption (A5) becomes $\Sigma/b$, where $\Sigma := \text{Cov}[\nabla f(\xi, \theta^*)]$ is the covariance matrix for a single sample defined in Assumption (B5). The estimates are still asymptotically efficient, which is given by the following corollary of Theorems 4.1 and 4.5.

**Corollary 5.1.** Under Assumptions (B1) to (B5) (given in Appendix C) made on the oracles defined by (9), and assuming that (B4)b) holds with $q' > \frac{2}{\gamma}$, the estimator $\theta_n$ of the streaming mSNA algorithm with step-size $\alpha_n = 1/(n + c'_\alpha)$ and mini-batch size $b$ is asymptotically efficient:

$$ \sqrt{N}(\theta_n - \theta^*) \xrightarrow[n \to +\infty]{\mathcal{L}} \mathcal{N}(0, H^{-1}\Sigma H^{-1}), $$

where $\Sigma := \text{Cov}[\nabla f(\xi, \theta^*)]$. Moreover, with Assumption (A6), the estimator $\bar{\theta}_n$ of the streaming averaged mSNA algorithm with step-size exponent $\alpha < 1$ and with mini-batch size $b$ and assuming that (B4) holds with $q > \frac{2}{\alpha}$ and $q' > \frac{2}{\gamma}$ is asymptotically efficient:

$$ \sqrt{N}(\bar{\theta}_n - \theta^*) \xrightarrow[n \to +\infty]{\mathcal{L}} \mathcal{N}(0, H^{-1}\Sigma H^{-1}). $$

Notably, the streaming mSNA algorithm with $\alpha_n := 1/(n + c'_\alpha)$, as well as the streaming averaged mSNA algorithm with any $(\alpha, c_\alpha, c'_\alpha)$-sequence step-size with $\alpha < 1$, attain both the asymptotic efficiency and possess the same asymptotic distribution as the purely online versions with batch-size $b = 1$. Observe that in the case of the streaming averaged mSNA, one can take $c_\alpha = b^a$ with $a \in [0, 1]$ (usually $a = 1 - \alpha$). This enables larger steps when the number of iterations decreases (see (Godichon-Baggioni et al., 2023) for more details).

## 5.2 Discussion on the Computational Complexity

The per-iteration workload of our algorithm can be summarized as follows:

**Computing $\tilde{H}_n$.** The complexity of computing the selected columns of the Hessian oracle depends on the function $f$. We assume that computing the product of the Hessian $\nabla^2 f(\xi_{n,i}, \theta_{n-1})$ with a single vector (a Hessian-vector product) can be done efficiently, without forming the full matrix (see, e.g., (Pearlmutter, 1994)). Computing the $\ell$ columns selected by $M_{\ell,n}$ for a single sample $\xi_{n,i}$ thus requires $O(\ell d)$ operations in many standard cases. Computing $\tilde{H}_n$ involves computing these $\ell$ columns for each of the $b$ samples and averaging them, resulting in a complexity of order $O(\ell bd)$.
**Calculating** $\|\tilde{H}_n\|_{op}$. The calculation of the operator norm of $\tilde{H}_n$ requires $O(\ell^2d)$ operations, as only $\ell$ columns of $\tilde{H}_n$ are non-zero. As long as $\ell \le \sqrt{d}$, this is of order $O(d^2)$.

**Updating** $A_n$. Calculating the products $A_{n-1}\tilde{H}_n^T$ and $\tilde{H}_nA_{n-1}\tilde{H}_n^T$ requires $O(\ell d^2)$ operations. The overall update for $A_n$ is thus dominated by these matrix multiplications, requiring $O(\ell d^2)$ operations.

**Updating** $\theta_n$. This first requires $O(bd)$ operations to calculate the gradient oracle $g_n(\theta_{n-1})$, which is an average over the mini-batch. Computing the product $(A_{n-1} + \nu_n I_d) g_n(\theta_{n-1})$ then requires $O(d^2)$ operations. The total cost is $O(bd + d^2)$.

**Updating** $\bar{\theta}_n$. This operations requires $O(d)$ operations (vector scaling and addition).

**Total cost of an iteration of the algorithm** Summing these costs, the total cost per iteration is:

$$ \underbrace{O(\ell bd + \ell d^2)}_{\text{Update } A_n} + \underbrace{O(bd + d^2)}_{\text{Update } \theta_n} + \underbrace{O(d)}_{\text{Update } \bar{\theta}_n \text{ (if used)}} = O(\ell (bd + d^2)) $$

Therefore, as long as the batch-size $b$ is not larger than the dimension $d$, the cost of an iteration is of order $O(\ell d^2)$. In a streaming setting with a single pass over $N$ data points, the algorithm performs $N/b$ iterations, leading to a total complexity of $O(N\bar{\ell}d^2)$. Hence, by choosing $b = d$ and $\ell = 1$, the algorithm achieves a total complexity for one pass of order $O(Nd)$ operations.

# 6 Numerical Experiments

In this section, we empirically evaluate the performance of the proposed masked Stochastic Newton Algorithm (mSNA) and its weighted averaged variant, particularly focusing on the streaming setting. We compare our methods against Stochastic Gradient Descent (SGD) and its averaged variant on both synthetic and real-world datasets. Our implementation uses the PyTorch framework and the code is available at https://github.com/guillaume-salle/SNA. All experiments were conducted on a commodity laptop, with computations performed on a multi-core CPU using 6 threads for parallelism; no GPU acceleration was utilized. Although our simulations have focused on classical examples of linear and logistic regression, there are many other applications, see for instance Godichon-Baggioni et al. (2025).

## 6.1 Experimental Setup

We compare the performance of the following algorithms:
* **SGD**: The standard stochastic gradient descent update defined by

$$ \theta_n = \theta_{n-1} - \alpha_n g_n(\theta_{n-1}). \quad (10) $$

* **Averaged SGD**: The weighted-average variant with estimator $\bar{\theta}_n$ defined by (7).

* **mSNA**: Our proposed masked Stochastic Newton Algorithm defined by (4), (5), and (6).

* **Averaged mSNA**: The weighted-average variant of our method with estimator $\bar{\theta}_n$ defined by (7).

For all optimizers, we use a batch size equal $b$ to the dimension $d$ of $\theta$.

For non-averaged methods SGD and mSNA, we use the standard step size $\alpha_n := 1/(n + n_0)$. For the averaged SGD and averaged mSNA methods, we use the larger step size $\alpha_n := d^{0.25}/(n^{0.75} + d^{0.25} * n_0)$, and $\tau = 2$ for the averaging step (7). The hyper-parameter $n_0$ determines the initial step size $\alpha_1$, and is linked to the expected smoothness constant $\mathcal{L}$ of the objective function. It is set to $n_0 = d$ for synthetic data, and tuned by line search for real data.

In addition, in the case of the mSNA and averaged mSNA algorithms, the parameter $\ell$ determines the number of columns and lines of the conditioning matrix $A_n$ that are updated at each iteration. Three choices of $\ell$ are considered: $\ell = 1, d^{0.25}, \sqrt{d}$. We use the step size $\gamma_n = 1/(n^{0.75} + n_0)$ for the conditioning matrix update (6), and the regularization parameter $\nu = 0$ for the $\theta_n$ update (4). Note that for averaged mSNA, this choice of $\gamma = 0.75$ is the edge case of the condition $\alpha + \gamma > 3/2$ of Theorem 4.5 for the convergence of $\theta_n$.

## 6.2 Experiments with Synthetic Data

To assess the performance of the different methods, we analyze the evolution of the quadratic error $\|\theta_n - \theta^*\|^2$ with respect to the sample size. We also examine the evolution of the squared Frobenius norm error $\|A_n - H^{-1}\|_F^2$ for the inverse Hessian estimates, as well as the total computational time of the different algorithms. Observe that for both methods, the initialization $\theta_0$ is randomly chosen on the unit sphere centered at $\theta^*$ and $A_0$ is set to $I_d$.

### 6.2.1 Linear Regression

We consider a linear regression problem where $F(\theta) = \frac{1}{2}\mathbb{E}[(Y - X^T\theta)^2]$, with $\theta \in \mathbb{R}^d$. It is an interesting case since the Hessian is known and constant, given by $\nabla^2 F(\theta) = \mathbb{E}[XX^T] =: \Sigma_X$ for all $\theta \in \mathbb{R}^d$, and we can so calculate explicitly $H^{-1}$. To create a challenging scenario, we consider the following model:

* $X$ follows a multivariate Gaussian distribution $N(0, \Sigma_X)$. The covariance matrix $\Sigma_X$ is constructed to be ill-conditioned. Specifically, $\Sigma_X = U\Lambda U^T$, where $\Lambda$ is a diagonal matrix with eigenvalues spaced evenly between $10^{-2}$ and 1, and $U$ is
a random orthogonal matrix, yielding a condition number of $10^2$ for $\Sigma_X$. This construction ensures an ill-conditioned covariance matrix with non-axis-aligned principal components, making the optimization problem challenging even for conditioned SGD methods that use a diagonal conditioning matrix.

* $Y = X^T\theta^* + \epsilon$, where the parameter $\theta^*$ is randomly chosen from a standardized Gaussian distribution and $\epsilon \sim N(0, 1)$.

We take $d = 1000$ and $N = 10^7$, in order to have a moderately large dimension while still having a reasonable number of iterations of the algorithms with a batch size equal to $d$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Linear regression with $d = 1000$. From the left to the right, comparison of: evolution of the quadratic errors $\|\theta_n - \theta^*\|^2$, quadratic errors $\|A_n - H^{-1}\|_F^2$, computational time for the different methods.</figcaption></figure>

In Figure 1, we can see that, as expected in an ill-conditioned setting, standard gradient-based methods fail to converge effectively. In contrast, Newton-type methods demonstrate significantly better performance, with only a slightly higher computational cost. For all algorithmic variants considered, we observe that the averaged versions consistently outperform their non-averaged counterparts. As anticipated, increasing the mask size $\ell$ leads to improved performance, albeit at the cost of increased computation time.

### 6.2.2 Logistic Regression

We consider here a logistic regression problem:

$$F(\theta) = \mathbb{E}[\log(1 + \exp(X^T\theta)) - X^T\theta Y],$$

with $X \sim N(0, \Sigma_X)$ as in Section 6.2.1 and $Y | X \sim B(\sigma(X^T\theta^*))$, where $\sigma(x) = \frac{e^x}{1+e^x}$ is the sigmoid function.

For this case, we do not have a closed-form expression for the Hessian at the minimizer: $H = \nabla^2 F(\theta^*) = \mathbb{E}[\sigma(X^T\theta^*)(1 - \sigma(X^T\theta^*))XX^T]$. Using a Monte-Carlo approach, we compute an empirical estimate of $H$ on a large sample, and use it to compute an estimate of the inverse $H^{-1}$.
<figure><img src="image_4.png" /><figcaption>Figure 2: Logistic regression with $d = 1000$. From the left to the right, comparison of: evolution of the quadratic errors $||\theta_n - \theta^*||^2$, quadratic errors $||A_n - H^{-1}||_F^2$, computational time for the different methods.</figcaption></figure>

In Figure 2, we observe similar trends as in the linear regression setting. In the case of logistic regression, the Hessian $H$ exhibits very small eigenvalues and a large conditioning number, which further diminish the performance of the SGD optimizer.

## 6.3 Ridge logistic regression on Real Data

We evaluate the algorithms on real-world datasets for binary classification using Ridge logistic regression. The characteristics of the studied datasets are summarized in Table 1.

<table><caption>Table 1: Key characteristics of the datasets used in this study.</caption><thead><tr><th></th><th>Dimension</th><th>Training Set Size</th><th>Init Set Size</th><th>Testing Set Size</th></tr></thead><tbody><tr><td>Adult</td><td>98</td><td>38,682</td><td>390</td><td>9,769</td></tr><tr><td>Connect-4</td><td>85</td><td>53,504</td><td>540</td><td>13,512</td></tr><tr><td>Covtype</td><td>55</td><td>460,160</td><td>4,648</td><td>116,203</td></tr><tr><td>Mnist</td><td>785</td><td>54,430</td><td>1,570</td><td>14,000</td></tr><tr><td>Mushrooms</td><td>95</td><td>6,309</td><td>190</td><td>1,625</td></tr><tr><td>Phishing</td><td>39</td><td>8,755</td><td>88</td><td>2,211</td></tr></tbody></table>

Each data set is split into training and test sets. To evaluate the performance of the different methods, we consider both the prediction accuracy and the empirical loss function evaluated both on the train set and the test set.

Adult (Becker and Kohavi, 1996) contains census data for income prediction and includes multiple categorical variables converted into binary indicators. Connect-4 (Tromp, 1995) is a dataset for the Connect-4 game, where the goal is to predict the next move based on the current board state. COVTYPE (Blackard, 1998) originally includes multiple forest cover types; in this study, we focus on distinguishing “Spruce/Fir” (labeled as 1) from all other categories (labeled as 0). MNIST (Lecun et al., 1998) is a dataset for handwritten digit recognition, where we focus on distinguishing the digit '0' (labeled as 1) from all other digits (labeled as 0). Mushrooms (UCI Machine Learning Repository,
1981) consists of morphological attributes of mushrooms used to determine their edibility. Phishing (Mohammad and McCluskey, 2012) is designed to detect malicious websites and consists entirely of categorical features, which we encode in binary form. These datasets are frequently adopted to benchmark binary classifiers (Toulis and Airoldi, 2017).

Our objective is to apply a ridge regression model to predict the binary response variable for each dataset. In the ridge regression, the objective function is defined by

$$F(\theta) = \mathbb{E} \left[ \log(1 + \exp(X^T \theta)) - X^T \theta Y + \frac{\lambda}{2} \|\theta\|^2 \right]$$

where $Y \in \{0, 1\}$, $X \in \mathbb{R}^d$ and $\lambda > 0$ is an $L_2$-regularization parameter. For all optimizers, we perform an initialization for the initial estimator $\theta_0$ by doing 100 steps of gradient descent over an initial batch of $\max\{N/100, 2d\}$ data, starting at 0 and using a constant learning rate found by line search. We do not include this initialization of $\theta_0$ in the computational time of the algorithms, since it is shared by all optimizers.

For the mSNA and averaged mSNA optimizers, we also perform an initialization for the conditioning matrix $A_0$ by computing on the same batch an estimator of the hessian of $F$ at $\theta_0$ obtained before, and we inverse it. This initialization has a computational cost of $O(d^3)$, and is included in the computational time.

In Table 2, we observe that the mSNA algorithms, along with their averaged versions have analogous behavior or can outperform gradient-based methods, with higher but comparable computation times. These results highlight not only the superior performance of our algorithms but also their practical readiness for deployment.

## Conclusion

This paper introduced an efficient online mini-batch stochastic Newton algorithm (mSNA) designed for smooth convex optimization problems in stochastic settings. Our work advances the Universal Stochastic Newton Algorithm (USNA) proposed in Godichon-Baggioni et al. (2025) by establishing, in the case of local strong convexity around the minimizer, the asymptotic efficiency of our algorithm without requiring iterate averaging, a notable theoretical improvement. This was made possible by proving a general theorem on the convergence of Robbins-Monro algorithms for linear functions, which is of independent interest.

The proposed mSNA algorithm operates in a streaming fashion, processing data in mini-batches, making it particularly well-suited for large-scale datasets. This enables our algorithm to achieve a total computational complexity of $O(Nd)$ for a single pass over $N$ data points. Numerical experiments on both synthetic and real-world datasets for linear and logistic regression, with or without ridge penalization, demonstrate that mSNA and its averaged variant consistently achieve competitive performance compared to SGD and averaged SGD. These results highlight not only the strong theoretical guarantees of our method but also its practical readiness for real-world problems.
<table><caption>Table 2: Performance of streaming optimizers on various datasets.</caption><thead><tr><th>Dataset</th><th>Optimizer</th><th>Train Acc</th><th>Test Acc</th><th>Train Loss</th><th>Test Loss</th><th>Time</th></tr></thead><tbody><tr><td rowspan="8">adult</td><td>SGD-Avg</td><td><strong>76.86</strong></td><td><strong>76.87</strong></td><td>5.72e-01</td><td>5.75e-01</td><td>31.8 ms</td></tr><tr><td>SGD</td><td>76.75</td><td>76.72</td><td><strong>5.67e-01</strong></td><td><strong>5.70e-01</strong></td><td><strong>27.0 ms</strong></td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>151.4 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>149.8 ms</td></tr><tr><td>mSNA-Avg</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>78.6 ms</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>131.3 ms</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>146.7 ms</td></tr><tr><td>mSNA</td><td>76.17</td><td>76.17</td><td>5.74e-01</td><td>5.77e-01</td><td>78.1 ms</td></tr><tr><td rowspan="8">connect-4</td><td>SGD-Avg</td><td>77.04</td><td>76.72</td><td>4.68e-01</td><td>4.69e-01</td><td><strong>42.2 ms</strong></td></tr><tr><td>SGD</td><td>77.00</td><td>76.61</td><td>4.69e-01</td><td>4.70e-01</td><td>44.8 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td><strong>79.26</strong></td><td><strong>79.18</strong></td><td>4.38e-01</td><td>4.39e-01</td><td>205.0 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>79.26</td><td>79.18</td><td><strong>4.38e-01</strong></td><td><strong>4.39e-01</strong></td><td>217.4 ms</td></tr><tr><td>mSNA-Avg</td><td>79.25</td><td>79.17</td><td>4.38e-01</td><td>4.39e-01</td><td>123.5 ms</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>79.16</td><td>79.08</td><td>4.39e-01</td><td>4.41e-01</td><td>201.5 ms</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td>79.14</td><td>79.10</td><td>4.39e-01</td><td>4.41e-01</td><td>217.8 ms</td></tr><tr><td>mSNA</td><td>79.15</td><td>79.05</td><td>4.39e-01</td><td>4.41e-01</td><td>119.8 ms</td></tr><tr><td rowspan="8">covtype</td><td>SGD-Avg</td><td><strong>58.36</strong></td><td><strong>58.23</strong></td><td><strong>6.66e-01</strong></td><td><strong>6.66e-01</strong></td><td>527.6 ms</td></tr><tr><td>SGD</td><td>57.56</td><td>57.35</td><td>6.70e-01</td><td>6.71e-01</td><td><strong>494.9 ms</strong></td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td>55.64</td><td>55.68</td><td>6.79e-01</td><td>6.79e-01</td><td>2.35 s</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>55.64</td><td>55.68</td><td>6.79e-01</td><td>6.79e-01</td><td>2.22 s</td></tr><tr><td>mSNA-Avg</td><td>55.64</td><td>55.68</td><td>6.79e-01</td><td>6.79e-01</td><td>1.51 s</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>55.62</td><td>55.65</td><td>6.79e-01</td><td>6.79e-01</td><td>2.31 s</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td>55.62</td><td>55.65</td><td>6.79e-01</td><td>6.79e-01</td><td>2.19 s</td></tr><tr><td>mSNA</td><td>55.62</td><td>55.65</td><td>6.79e-01</td><td>6.79e-01</td><td>1.47 s</td></tr><tr><td rowspan="8">mnist</td><td>SGD-Avg</td><td>86.76</td><td>86.66</td><td>5.05e-01</td><td>5.21e-01</td><td>23.8 ms</td></tr><tr><td>SGD</td><td><strong>86.95</strong></td><td><strong>87.11</strong></td><td><strong>3.09e-01</strong></td><td><strong>3.14e-01</strong></td><td><strong>19.0 ms</strong></td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td>86.08</td><td>85.76</td><td>3.28e-01</td><td>3.35e-01</td><td>102.3 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>86.08</td><td>85.76</td><td>3.28e-01</td><td>3.35e-01</td><td>98.8 ms</td></tr><tr><td>mSNA-Avg</td><td>86.08</td><td>85.80</td><td>3.28e-01</td><td>3.35e-01</td><td>65.0 ms</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>85.75</td><td>85.49</td><td>3.34e-01</td><td>3.41e-01</td><td>243.2 ms</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td>85.75</td><td>85.49</td><td>3.34e-01</td><td>3.41e-01</td><td>133.7 ms</td></tr><tr><td>mSNA</td><td>85.76</td><td>85.49</td><td>3.35e-01</td><td>3.41e-01</td><td>50.2 ms</td></tr><tr><td rowspan="8">mushrooms</td><td>SGD-Avg</td><td>94.37</td><td>95.08</td><td>2.43e-01</td><td>2.19e-01</td><td>9.3 ms</td></tr><tr><td>SGD</td><td>94.37</td><td>95.08</td><td>2.41e-01</td><td>2.17e-01</td><td><strong>4.6 ms</strong></td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td>99.18</td><td>98.95</td><td>3.50e-02</td><td>4.38e-02</td><td>47.1 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>99.18</td><td>98.95</td><td>3.50e-02</td><td>4.37e-02</td><td>48.6 ms</td></tr><tr><td>mSNA-Avg</td><td>99.18</td><td>98.95</td><td>3.50e-02</td><td>4.38e-02</td><td>27.8 ms</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>99.69</td><td>99.08</td><td><strong>2.98e-02</strong></td><td>3.83e-02</td><td>47.4 ms</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td><strong>99.70</strong></td><td><strong>99.20</strong></td><td>2.98e-02</td><td><strong>3.81e-02</strong></td><td>49.3 ms</td></tr><tr><td>mSNA</td><td>99.66</td><td>99.08</td><td>2.98e-02</td><td>3.84e-02</td><td>30.3 ms</td></tr><tr><td rowspan="8">phishing</td><td>SGD-Avg</td><td>89.45</td><td>89.69</td><td>2.56e-01</td><td>2.46e-01</td><td>20.2 ms</td></tr><tr><td>SGD</td><td>89.26</td><td>89.51</td><td>2.59e-01</td><td>2.49e-01</td><td><strong>13.3 ms</strong></td></tr><tr><td>mSNA-Avg (ℓ = 0.25)</td><td>93.79</td><td>94.08</td><td>1.61e-01</td><td>1.63e-01</td><td>95.0 ms</td></tr><tr><td>mSNA-Avg (ℓ = 0.5)</td><td>93.77</td><td>94.08</td><td><strong>1.61e-01</strong></td><td><strong>1.62e-01</strong></td><td>124.2 ms</td></tr><tr><td>mSNA-Avg</td><td>93.75</td><td>94.08</td><td>1.62e-01</td><td>1.63e-01</td><td>52.5 ms</td></tr><tr><td>mSNA (ℓ = 0.25)</td><td>93.77</td><td>94.12</td><td>1.66e-01</td><td>1.66e-01</td><td>78.8 ms</td></tr><tr><td>mSNA (ℓ = 0.5)</td><td><strong>93.80</strong></td><td><strong>94.17</strong></td><td>1.65e-01</td><td>1.65e-01</td><td>95.4 ms</td></tr><tr><td>mSNA</td><td>93.77</td><td>94.08</td><td>1.67e-01</td><td>1.66e-01</td><td>62.7 ms</td></tr></tbody></table>
**Limitations and Future Work:** The current work focuses on the convex setting, and several avenues for refinement and extension within this scope exist.

First, formulating the problem of matrix inverse estimation as an optimization problem opens the door to the possibility of estimating the inverse under structural constraints, such as Ridge or Lasso regularization. This could be valuable when the matrix is sparse or when its smallest eigenvalue is too close to zero.

Second, exploring the robustness of the algorithm under weaker conditions, such as relaxing the smoothness assumptions (e.g., (A2) b)), would be a valuable direction for future work. Establishing the asymptotic efficiency of the averaged mSNA without relying on the Hessian Lipschitz continuity assumption (A6) remains an open and interesting challenge.

Finally, while the random mask strategy used for Hessian approximation is both simple and effective, investigating alternative sketching techniques or adaptive mechanisms for selecting Hessian information could lead to further improvements.

## References

Agarwal, N., Bullins, B., and Hazan, E. (2017). Second-order stochastic optimization for machine learning in linear time. *Journal of Machine Learning Research*, 18(116):1–40.

Bach, F. and Moulines, E. (2011). Non-asymptotic analysis of stochastic approximation algorithms for machine learning. In *Advances in Neural Information Processing Systems*.

Becker, B. and Kohavi, R. (1996). Adult. UCI Machine Learning Repository.

Bercu, B., Godichon, A., and Portier, B. (2020). An efficient stochastic newton algorithm for parameter estimation in logistic regressions. *SIAM Journal on Control and Optimization*, 58(1):348–367.

Blackard, J. (1998). Covertype. UCI Machine Learning Repository.

Bottou, L., Curtis, F. E., and Nocedal, J. (2018). Optimization methods for large-scale machine learning. *SIAM Review*, 60(2):223–311.

Boyer, C. and Godichon-Baggioni, A. (2023). On the asymptotic rate of convergence of stochastic newton algorithms and their weighted averaged versions. *Computational Optimization and Applications*, 84(3):921–972.

Byrd, R. H., Chin, G. M., Neveitt, W., and Nocedal, J. (2011). On the use of stochastic hessian information in optimization methods for machine learning. *SIAM Journal on Optimization*, 21(3):977–995.

Carpentier, P., Chancelier, J.-P., Cohen, G., and De Lara, M. (2015). *Stochastic Multi-Stage Optimization: At the Crossroads between Discrete Time Stochastic Control and Stochastic Programming*, volume 75 of *Probability Theory and Stochastic Modelling*. Springer International Publishing.
Cénac, P., Godichon-Baggioni, A., and Portier, B. (2025). An efficient averaged stochastic Gauss-Newton algorithm for estimating parameters of nonlinear regressions models. *Bernoulli*, 31(1):1–29.

Dennis, Jr., J. E. and Moré, J. J. (1977). Quasi-newton methods, motivation and theory. *SIAM Review*, 19(1):46–89.

Duchi, J., Hazan, E., and Singer, Y. (2011). Adaptive subgradient methods for online learning and stochastic optimization. *Journal of Machine Learning Research*, 12(61):2121–2159.

Duflo, M. (1996). *Algorithmes stochastiques*. Mathématiques et Applications. Springer Berlin, Heidelberg.

Godichon-Baggioni, A. (2019). Online estimation of the asymptotic variance for averaged stochastic gradient algorithms. *Journal of Statistical Planning and Inference*, 203:1–19.

Godichon-Baggioni, A., Lu, W., and Portier, B. (2025). Online estimation of the inverse of the hessian for stochastic optimization with application to universal stochastic newton algorithms. Preprint, available at https://hal.science/hal-04391570.

Godichon-Baggioni, A. and Werge, N. (2025). On adaptive stochastic optimization for streaming data: A newton's method with o(dN) operations. *Journal of Machine Learning Research*, 26(59):1–49.

Godichon-Baggioni, A., Werge, N., and Wintenberger, O. (2023). Non-asymptotic analysis of stochastic approximation algorithms for streaming data. *ESAIM: Probability and Statistics*, 27:482–514.

Gower, R. M., Loizou, N., Qian, X., Sailanbayev, A., Shulgin, E., and Richtárik, P. (2019). SGD: General analysis and improved rates. In Chaudhuri, K. and Salakhutdinov, R., editors, *Proceedings of the 36th International Conference on Machine Learning*, volume 97 of *Proceedings of Machine Learning Research*, pages 5200–5209. PMLR.

Gower, R. M. and Richtárik, P. (2017). Randomized quasi-newton updates are linearly convergent matrix inversion algorithms. *SIAM Journal on Matrix Analysis and Applications*, 38(4):1380–1409. Publisher: Society for Industrial and Applied Mathematics.

Hinton, G. (2012). Neural networks for machine learning, lecture 6. Online Course, Coursera.

Kingma, D. P. and Ba, J. (2017). Adam: A method for stochastic optimization.

Lecun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998). Gradient-based learning applied to document recognition. *Proceedings of the IEEE*, 86(11):2278–2324.

Leluc, R. and Portier, F. (2022). Sgd with coordinate sampling: Theory and practice. *Journal of Machine Learning Research*, 23(342):1–47.
Leluc, R. and Portier, F. (2023). Asymptotic analysis of conditioned stochastic gradient descent. *Transactions on Machine Learning Research*.

Mohammad, R. and McCluskey, L. (2012). Phishing Websites. UCI Machine Learning Repository.

Mokkadem, A. and Pelletier, M. (2011). A generalization of the averaging procedure: The use of two-time-scale algorithms. *SIAM Journal on Control and Optimization*, 49(4):1523-1543.

Moritz, P., Nishihara, R., and Jordan, M. (2016). A linearly-convergent stochastic l-BFGS algorithm. In *Proceedings of the 19th International Conference on Artificial Intelligence and Statistics*, pages 249-258. PMLR.

Pearlmutter, B. A. (1994). Fast exact multiplication by the hessian. *Neural Computation*, 6(1):147-160.

Pelletier, M. (1998a). On the almost sure asymptotic behaviour of stochastic algorithms. *Stochastic Processes and their Applications*, 78(2):217-244.

Pelletier, M. (1998b). Weak convergence rates for stochastic approximation with application to multiple targets and simulated annealing. *The Annals of Applied Probability*, 8(1):10-44.

Pelletier, M. (2000). Asymptotic almost sure efficiency of averaged stochastic algorithms. *SIAM Journal on Control and Optimization*, 39(1):49-72.

Polyak, B. T. (1990). A new method of stochastic approximation type. *Avtomatika i Telemekhanika*, (7):98-107. In Russian. English translation in: Automation and Remote Control, 51(7):937-946, 1990.

Polyak, B. T. and Juditsky, A. B. (1992). Acceleration of stochastic approximation by averaging. *SIAM Journal on Control and Optimization*, 30(4):838-855.

Ruppert, D. (1988). Efficient estimations from a slowly convergent robbins-monro process. Technical report, Cornell University Operations Research and Industrial Engineering.

Schraudolph, N. N., Yu, J., and Günter, S. (2007). A stochastic quasi-newton method for online convex optimization. In *Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics*, pages 436-443. PMLR.

Surendran, S., Fermanian, A., Godichon-Baggioni, A., and Le Corff, S. (2024). Non-asymptotic analysis of biased adaptive stochastic approximation. *Advances in Neural Information Processing Systems*, 37:12897-12943.

Toulis, P. and Airoldi, E. M. (2017). Asymptotic and finite-sample properties of estimators based on stochastic gradients. *The Annals of Statistics*, 45(4):1694-1727.

Tromp, J. (1995). Connect-4. UCI Machine Learning Repository.
# A Proofs

## A.1 Proof of Proposition 3.1

*Proof of Proposition 3.1.* The following lemma bounds almost surely and asymptotically the operator norm of $A_n$, without relying on information on the behavior of $\theta_n$.

**Lemma A.1.** Under Assumptions (A1) and (A2)b), for any $\eta > 0$, we have:

$$ \lambda_{\max}(A_n) = \mathcal{O}(n^{1-\gamma}(\ln n)^{1+\eta}) \quad \text{a.s.} $$

The proof of this lemma is technical and is deferred to the next section A.3.

We have

$$ A_n = A_{n-1} - \gamma_n (\nabla J_{\text{sym}}(A_{n-1}) + \epsilon_n + r_n), $$

with a noise

$$ \epsilon_n := (\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} - \gamma_n \tilde{H}_n A_{n-1} \tilde{H}_n^T) \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} \\ - \mathbb{E} \left[ (\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} - \gamma_n \tilde{H}_n A_{n-1} \tilde{H}_n^T) \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} | \mathcal{F}_{n-1} \right], $$

and residual terms

$$ r_n := r_n^{(1)} + r_n^{(2)} + r_n^{(3)} \\ r_n^{(1)} := A_{n-1} (\nabla^2 F(\theta_{n-1}) - H) + (\nabla^2 F(\theta_{n-1}) - H) A_{n-1} $$

$$ r_n^{(2)} := \mathbb{E} \left[ \gamma_n \tilde{H}_n A_{n-1} \tilde{H}_n^T \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} | \mathcal{F}_{n-1} \right] $$

$$ r_n^{(3)} := -\mathbb{E} \left[ (\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T + 2M_{n,\ell}) \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} > \frac{1}{2}} | \mathcal{F}_{n-1} \right]. $$

The term $r^{(1)}$ accounts for the error of the estimator $\theta_n$, leading to a bias in the estimation of $H$. The term $r^{(2)}$ accounts for the term added in order to ensure the positivity of $A_n$, and the term $r^{(3)}$ accounts for the truncation. By construction, the noise $\epsilon_n$ is such that $\mathbb{E}[\epsilon_n | \mathcal{F}_{n-1}] = 0$. Moreover, using inequality $(a+b)^p \le 2^{p-1}(a^p + b^p)$ for $p \ge 2$, we have

$$ \begin{aligned} \mathbb{E} \left[ \|\epsilon_n\|_F^2 | \mathcal{F}_{n-1} \right] & \le \mathbb{E} \left[ \left\| \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T + 2M_{n,\ell} + \gamma_n \tilde{H}_n A_{n-1} \tilde{H}_n^T \right\|_F^2 \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} | \mathcal{F}_{n-1} \right] \\ & \le \mathbb{E} \left[ \left( 2 \|A_{n-1}\|_F \|\tilde{H}_n\|_F + 2 \|M_{n,\ell}\|_F + \frac{1}{2} \|A_{n-1}\|_F \|\tilde{H}_n\|_F \right)^2 | \mathcal{F}_{n-1} \right] \\ & \le \mathbb{E} \left[ 2 \times \left( \frac{5}{2} \|A_{n-1}\|_F \|\tilde{H}_n\|_F \right)^2 + 2 \times (2 \|M_{n,\ell}\|_F)^2 | \mathcal{F}_{n-1} \right] \\ & = \mathcal{O} \left( 1 + \|A_{n-1}\|_F^2 \right). \end{aligned} $$
Since $\theta_n \xrightarrow{n \to \infty} \theta^*$ and using Assumption (A4)b), for $n$ large enough $\|\epsilon_n\|_F$ has a moment of order $q'$ and we have with the same reasoning as order 2:

$$ \mathbb{E}[\|\epsilon_n\|_F^{q'} | \mathcal{F}_{n-1}] = \mathcal{O}(1 + \|A_{n-1}\|_F^{q'}) . $$

For the residual term $r_n$, taking a non-increasing modulus of continuity $\mu$ of $\nabla^2 F$ in $\theta^*$, we have

$$ \|\nabla^2 F(\theta_{n-1}) - H\| \le \mu(\|\theta_{n-1} - \theta^*\|) \le \mu(\delta_n). $$

We get

$$ \|r_n^{(1)}\|_F \le 2 \|A_{n-1}\|_F \|\nabla^2 F(\theta_{n-1}) - H\|_F = \mathcal{O}(\|A_{n-1}\|_F \mu(\delta_n)), $$

$$ \|r_n^{(2)}\|_F \le \gamma_n \|A_{n-1}\|_F \mathbb{E} \left[ \|\tilde{H}_n\|_F^2 \mid \mathcal{F}_{n-1} \right] = \mathcal{O}(\gamma_n \|A_{n-1}\|_F) \text{ a.s.,} $$

$$ \begin{aligned} \|r_n^{(3)}\|_F &\le \mathbb{E} \left[ \left( 2 \|A_{n-1}\|_F \|\tilde{H}_n\|_F + 2 \|M_{n,l}\|_F \right) \times 2\gamma_n \|\tilde{H}_n\|_F \mid \mathcal{F}_{n-1} \right] \\ &= \mathcal{O}(\gamma_n (1 + \|A_{n-1}\|_F)) \text{ a.s.} \end{aligned} $$

Therefore,

$$ \|r_n\|_F = \mathcal{O}((1 + \|A_{n-1}\|_F)(\mu(\delta_n) + \gamma_n)) \text{ a.s.} $$

By applying the Proposition B.1 with $a_n := \mu(\delta_n) + \gamma_n$ which is non-increasing, we obtain for any $\eta > 0$:

$$ \begin{aligned} \|A_n - H^{-1}\|_F^2 &= \mathcal{O}\left(\gamma_n(\ln n)^{1+\eta} + \mu\left(\delta_{\lceil \frac{n}{2} \rceil}\right)^2 + \gamma_{\lceil \frac{n}{2} \rceil}^2\right) \quad \text{a.s.} \\ &= \mathcal{O}\left(\gamma_n(\ln n)^{1+\eta} + \mu\left(\delta_{\lceil \frac{n}{2} \rceil}\right)^2\right) \quad \text{a.s.} \end{aligned} $$

since $\gamma_{\lceil \frac{n}{2} \rceil}^2 = o(\gamma_n)$, which concludes this proof. $\square$

## A.2 Proof of Theorem 4.1

*Proof of Theorem 4.1.* We first prove the strong consistency of $(\theta_n)$. Recall that by induction on the factorization (2), $A_n$ is positive definite for any $n \ge 0$. With the help of Lemma A.1, we have for any $\eta > 0$:

$$ \lambda_{\min}(A_n + \nu_n I_d) \ge \nu_n. \quad \text{and} \quad \lambda_{\max}(A_n + \nu_n I_d) = \mathcal{O}(n^{1-\gamma}(\ln n)^{1+\eta}) \quad \text{a.s.} $$

By definition of $(\nu_n)$, and if $\alpha + \gamma > 3/2$, we have

$$ \sum_{n \ge 1} \alpha_n \nu_n = +\infty \quad \text{and} \quad \sum_{n \ge 1} \alpha_n n^{1-\gamma} (\ln n)^{1+\eta} < +\infty. $$

Therefore, we can apply Theorem 1 from Godichon-Baggioni and Werge (2025) to obtain the strong consistency of $(\theta_n)$ towards $\theta^*$. Then, by applying Theorem B.1, we obtain
an almost sure rate of convergence of $(A_n)$ toward $H^{-1}$, and in particular the strong consistency of $(A_n)$. Hence, if $1 - \zeta > 0$ we can apply Theorem 2 in Leluc and Portier (2023) to obtain the asymptotic normality of $\theta_n$:

$$ \frac{1}{\sqrt{\alpha_n}} (\theta_n - \theta^*) \xrightarrow{n \to +\infty} \mathcal{N} \left( 0, \frac{1}{2(1 - \zeta)} H^{-1} \Gamma H^{-1} \right). $$

To obtain the rate of almost sure convergence with Assumption (A6), from Boyer and Godichon-Baggioni (2023), if $\alpha = 1$ we apply Theorem 3.2, and if $\alpha \in (2/b, 1)$ we apply Theorem 4.2. $\square$

*Proof of Theorem 4.5.* The norm of the averaged matrices $\bar{A}_n$ still verify the upper bound of Lemma A.1, since $\bar{A}_n$ is a convex combination of $A_j$ for $j \le n$ and the upper bound is increasing. Hence, we can follow the same reasoning as the proof of Theorem 4.1 to obtain the results on $\theta_n$ and $A_n$.

By application of the Toeplitz lemma, the strong consistency of $\theta_n$ (resp. $A_n$) gives the strong consistency of $\bar{\theta}_n$ (resp. $\bar{A}_n$). To obtain the convergence rate and asymptotic normality of $\bar{\theta}_n$, we can apply Theorem 3 in Godichon-Baggioni and Werge (2025). Their Assumption 5 is implied by Assumptions (A6) and (A2)b). Then, Theorem B.1 with Assumption (A6) gives us the convergence rate of $\bar{A}_n$. $\square$

## A.3 Control on the largest eigenvalue of $A_n$

This proof of Lemma A.1 follows closely the proof of Proposition 6.1 in Godichon-Baggioni et al. (2025).

*Proof of Lemma A.1.* Let $n \ge 1$. We have:

$$ \mathbb{E} [\|A_n\|_F^2 | \mathcal{F}_{n-1}] = \|A_{n-1}\|_F^2 + 2\langle A_{n-1}, \mathbb{E}[A_n - A_{n-1} | \mathcal{F}_{n-1}] \rangle_F + \mathbb{E} [\|A_n - A_{n-1}\|_F^2 | \mathcal{F}_{n-1}]. \quad (11) $$

For the last term of the right-hand side of (11), we have using the truncation and then the inequality $(a+b)^2 \le 2a^2 + 2b^2$:

$$ \begin{aligned} & \mathbb{E} [\|A_n - A_{n-1}\|_F^2 | \mathcal{F}_{n-1}] \\ &= \mathbb{E} \left[ \left\| \gamma_n \left( \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} \right) - \gamma_n^2 \tilde{H}_n A_{n-1} \tilde{H}_n^T \right\|_F^2 \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} | \mathcal{F}_{n-1} \right] \\ &\le \gamma_n^2 \mathbb{E} \left[ \left( 2 \|\tilde{H}_n\|_F \|A_{n-1}\|_F + 2 \|M_{n,\ell}\|_F + \gamma_n \|\tilde{H}_n\|_F^2 \|A_{n-1}\|_F \right)^2 \mathbf{1}_{\gamma_n \|\tilde{H}_n\|_{op} \le \frac{1}{2}} | \mathcal{F}_{n-1} \right] \\ &\le \gamma_n^2 \mathbb{E} \left[ \left( 2 \|\tilde{H}_n\|_F \|A_{n-1}\|_F + 2 \|M_{n,\ell}\|_F + \frac{1}{2} \|\tilde{H}_n\|_F \|A_{n-1}\|_F \right)^2 | \mathcal{F}_{n-1} \right] \\ &\le \gamma_n^2 \left( 2 \times \frac{25}{4} \mathbb{E} [\|\tilde{H}_n\|_F^2 | \mathcal{F}_{n-1}] \|A_{n-1}\|_F^2 + 2 \times 4 \mathbb{E} [\|M_{n,\ell}\|_F^2 | \mathcal{F}_{n-1}] \right) \\ &= \mathcal{O} \left( \gamma_n^2 (\|A_{n-1}\|_F^2 + 1) \right) \text{ a.s.} \end{aligned} $$
Regarding the inner product term of (11), in order to form the expectation of $\tilde{H}_{n-1}$, we
can write

$$
\begin{align}
& \mathbb{E}[A_n - A_{n-1} \mid \mathcal{F}_{n-1}] \nonumber \\
&= \mathbb{E}\left[\left(-\gamma_n (\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell}) + \gamma_n^2 \tilde{H}_n A_{n-1} \tilde{H}_n^T\right) \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} \le \frac{1}{2}} \mid \mathcal{F}_{n-1}\right] \nonumber \\
&= -\gamma_n \mathbb{E}\left[\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} \mid \mathcal{F}_{n-1}\right] + \gamma_n^2 \mathbb{E}\left[\tilde{H}_n A_{n-1} \tilde{H}_n^T \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} \le \frac{1}{2}} \mid \mathcal{F}_{n-1}\right] \nonumber \\
&\quad + \gamma_n \mathbb{E}\left[\left(\tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell}\right) \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} > \frac{1}{2}} \mid \mathcal{F}_{n-1}\right]. \tag{12}
\end{align}
$$

We now bound the inner product of $A_{n-1}$ with each of the three terms of (12).

By Assumption (A1), the conditional expectation of $\tilde{H}_n$ is the hessian $\nabla^2 F(\theta_{n-1})$, which
is positive since $F$ is convex. We get:

$$
\begin{align*}
\langle A_{n-1}, \nabla^2 F(\theta_{n-1}) A_{n-1} \rangle_F &= \operatorname{Tr} (A_{n-1}^T \nabla^2 F(\theta_{n-1}) A_{n-1}) \ge 0 \\
\langle A_{n-1}, A_{n-1} \nabla^2 F(\theta_{n-1}) \rangle_F &= \operatorname{Tr} (A_{n-1}^T A_{n-1} \nabla^2 F(\theta_{n-1})) = \operatorname{Tr} (A_{n-1} \nabla^2 F(\theta_{n-1}) A_{n-1}^T) \ge 0.
\end{align*}
$$

Thus, for any $\zeta_n > 0$ which we will choose later in order to balance each component, we
get the following bound for the inner product of $A_{n-1}$ with the first term of (12):

$$
\begin{align*}
\langle A_{n-1}, -\gamma_n \mathbb{E} \left[ \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} \mid \mathcal{F}_{n-1} \right] \rangle_F &\le 2\gamma_n \langle A_{n-1}, \mathbb{E} [M_{n,\ell} \mid \mathcal{F}_{n-1}] \rangle_F \\
&\le \frac{\gamma_n}{\zeta_n} \|A_{n-1}\|_F^2 + \frac{\ell^2}{d^2} \gamma_n \zeta_n.
\end{align*}
$$

For the inner product of $A_{n-1}$ with the second term of (12), we get

$$
\langle A_{n-1}, \gamma_n^2 \mathbb{E} \left[ \tilde{H}_n A_{n-1} \tilde{H}_n^T \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} \le \frac{1}{2}} \mid \mathcal{F}_{n-1} \right] \rangle_F \le \gamma_n^2 \|A_{n-1}\|^2 \mathbb{E} \left[ \| \tilde{H}_n \|_F^2 \mid \mathcal{F}_{n-1} \right].
$$

For the inner product of $A_{n-1}$ with the third term of (12), we get:

$$
\begin{align*}
\langle A_{n-1}, \gamma_n \mathbb{E} \left[ \left( \tilde{H}_n A_{n-1} + A_{n-1} \tilde{H}_n^T - 2M_{n,\ell} \right) \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} > \frac{1}{2}} \mid \mathcal{F}_{n-1} \right] \rangle_F & \\
&\le 2\gamma_n \|A_{n-1}\|_F^2 \mathbb{E} \left[ \|\tilde{H}_n\|_F \mathbf{1}_{\gamma_n \| \tilde{H}_n \|_{op} > \frac{1}{2}} \mid \mathcal{F}_{n-1} \right] + 2\gamma_n \langle A_{n-1}, \mathbb{E}[M_{n,\ell} \mid \mathcal{F}_{n-1}] \rangle_F \\
&\le 2\gamma_n \|A_{n-1}\|_F^2 \mathbb{E} \left[ 2\gamma_n \|\tilde{H}_n\|_F^2 \mid \mathcal{F}_{n-1} \right] + \frac{\gamma_n}{\zeta_n} \|A_{n-1}\|_F^2 + \gamma_n \zeta_n d \\
&= \mathcal{O} \left( \left( \gamma_n^2 + \frac{\gamma_n}{\zeta_n} \right) \|A_{n-1}\|_F^2 + \gamma_n \zeta_n \right) \quad \text{a.s.}
\end{align*}
$$

Finally, we obtain

$$
\mathbb{E} [\|A_n\|_F^2 | \mathcal{F}_{n-1}] \le (1 + E_n) \|A_{n-1}\|_F^2 + B_n,
$$

with $E_n = O(\gamma_n^2 + \gamma_n/\zeta_n)$ and $B_n = O(\gamma_n^2 + \gamma_n\zeta_n)$. We can now set $\zeta_n := n^{1-\gamma}(\ln n)^{1+\eta}$ with $\eta > 0$, so that $(\gamma_n/\zeta_n)$ is summable. The Lemma B.4 applied with $a_n := \zeta_n^2$ and $D_n := 0$ gives us:

$$
\|A_{n-1}\|_F^2 = o(\zeta_{n-1}^2) \quad \text{a.s.}
$$

□
B Auxiliary results

**Theorem B.1.** Let $H \in \mathcal{M}_d(\mathbb{R})$ be a positive definite matrix and $z^* \in \mathbb{R}^d$. Consider the linear Robbins-Monro algorithm

$$z_n = z_{n-1} - \gamma_n (h(z_{n-1}) + \epsilon_n + r_n),$$

where $h(z) = H(z - z^*)$, $(\epsilon_n)$ and $(r_n)$ are random sequences adapted to a filtration $(\mathcal{F}_n)$, $z_0$ is an $\mathcal{F}_0$-measurable $\mathbb{R}^d$-valued random variable, and $\gamma_n = c/(n^\gamma + c')$ with $\gamma \in (\frac{1}{2}, 1)$, $c > 0$, and $c' \ge 0$. Assume there exists an exponent $p \ge 0$ such that $\|z_n\| = \mathcal{O}(n^p)$ a.s. Suppose further that for all $n \ge 1$,

$$\begin{align*}
\mathbb{E}[\epsilon_n | \mathcal{F}_{n-1}] &= 0 \text{ a.s.,} \\
\mathbb{E}[\|\epsilon_n\|^b | \mathcal{F}_{n-1}] &= \mathcal{O}((1 + \|z_{n-1}\|)^b) \text{ a.s.,} \\
\|r_n\| &= \mathcal{O}(a_n(1 + \|z_{n-1}\|)) \text{ a.s.,}
\end{align*}$$

for some $b > \frac{2}{\gamma}$, and a random sequence $(a_n)$ that almost surely is non-increasing and converging to 0. Then, for any $\eta > 0$,

$$\|z_n - z^*\|^2 = \mathcal{O}\left(\gamma_n(\ln n)^{1+\eta} + a_n^2 \left[\frac{1}{n/2}\right]\right) \text{ a.s.}$$

*Proof.* Define $\beta_{k,n} := \prod_{j=k}^n (I_d - \gamma_j H)$ for all $k \le n$, and $\beta_{k,n} := I_d$ for $k > n$. By induction, we express

$$z_n - z^* = \underbrace{\beta_{1,n}(z_0 - z^*)}_{Z_n} + \underbrace{\sum_{k=1}^{n} \beta_{k+1,n}\gamma_k \epsilon_k}_{E_n} + \underbrace{\sum_{k=1}^{n} \beta_{k+1,n}\gamma_k r_k}_{R_n} \quad (13)$$

Choose $n_0 \in \mathbb{N}$ such that for all $n \ge n_0$, $\lambda_{\max}(H)\gamma_n \le 1$. Since $H$ is positive, for $n \ge n_0$,

$$\|I_d - \gamma_n H\|_{op} \le \max\{1 - \lambda_{\min}(H)\gamma_n, 1 - \lambda_{\max}(H)\gamma_n\} = 1 - \lambda_{\min}(H)\gamma_n.$$

**Bounding $Z_n$.** The term $Z_n$ tends exponentially fast to 0. Indeed, for $n \ge n_0$,

$$\|\beta_{1,n}\|_{op} \le \|\beta_{1,n_0-1}\|_{op} \prod_{j=n_0}^{n} (1 - \lambda_{\min}(H)\gamma_j) \le \|\beta_{1,n_0-1}\|_{op} \exp \left( -\lambda_{\min}(H) \sum_{j=n_0}^{n} \gamma_j \right),$$

and $\sum_{j=n_0}^{n} \gamma_j \sim cn^{1-\gamma}/(1-\gamma)$.

**Bounding $E_n$.** Let $p \ge 0$ be such that $\|z_{n-1}\| = \mathcal{O}(n^p)$ a.s. Since $b > \frac{2}{\gamma}$, apply Theorem 6.1 and Remark 6.1 from Cénac et al. (2025), setting $T_n := (1 + \|z_{n-1}\|)I_d$ and $\xi_n := \epsilon_n/(1 + \|z_{n-1}\|)$. Then, $\|T_n\|_{op} = o(n^p(\ln n)^{n/2})$ a.s. for any $\eta > 0$. Although the theorem assumes non-negative exponents $a, b$ such that $\|T_n\|_{op} = o((\ln n)^b/n^a)$ a.s., it extends to negative values with the same proof, and we take $a := -p$. This yields

$$\|E_n\|^2 = \mathcal{O}(n^{2p}\gamma_n(\ln n)^{1+\eta}) \text{ a.s.}$$
**Bounding $R_n$.** We have for $n \ge n_0$,

$$ \|R_n\| \le \|I_d - \gamma_n H\|_{op} \|R_{n-1}\| + \gamma_n \|r_n\| \le (1 - \lambda_{\min}(H)\gamma_n) \|R_{n-1}\| + \gamma_n \|r_n\|. $$

Using equation (13) and the hypothesis on $(r_n)$, we have

$$ \|r_n\| = \mathcal{O} \left( a_n \left( 1 + \|Z_{n-1}\| + \|E_{n-1}\| + \|R_{n-1}\| \right) \right) \text{ a.s.} $$

We showed that $1 + \|Z_n\| + \|E_n\| = \mathcal{O}(\max\{1, n^{p-\gamma/2}(\ln n)^{(1+\eta)/2}\})$ a.s. Using the hypothesis on $(a_n)$, we can apply a.e. the deterministic result of Lemma B.2. This yields with some $\kappa > 0$:

$$ \|R_n\| = \mathcal{O} \left( \left( \exp(-\kappa n^{1-\gamma}) + a_{\lfloor n/2 \rfloor} \right) \max \left\{ 1, n^{p-\gamma/2} (\ln n)^{(1+\eta)/2} \right\} \right) \text{ a.s.} $$

We obtain with (13), after neglecting the exponential terms and using that $a_n = o(1)$ a.s.:

$$ \|z_n - z^*\| = \mathcal{O} \left( a_{\lfloor n/2 \rfloor} + n^{p-\gamma/2} (\ln n)^{(1+\eta)/2} \right) \text{ a.s.} \quad (14) $$

Let $k \in \mathbb{N}$ be the smallest integer such that $\|z_{n-1}\| = \mathcal{O}(n^{k\gamma/4})$ a.s. Applying Equation (14) with $p := k\gamma/4$, then using that $(\ln n)^{(1+\eta)/2} = o(n^{k\gamma/4})$ and $a_n = o(1)$ a.s., we get

$$ \|z_n\| \le \|z_n - z^*\| + \|z^*\| = \mathcal{O} \left( a_{\lfloor n/2 \rfloor} + n^{(k-2)\gamma/4} (\ln n)^{(1+\eta)/2} + 1 \right) \subset \mathcal{O} \left( n^{(k-1)\gamma/4} + 1 \right) \text{ a.s.} $$

If $k > 0$, it leads to $\|z_n\| = \mathcal{O}(n^{(k-1)\gamma/4})$ a.s., which is a contradiction with the definition of $k$. Therefore $k = 0$, and taking $p = 0$ in (14) gives us the desired rate of convergence. $\square$

**Lemma B.2.** Let $(r_n)_{n \ge 0}$, $(s_n)_{n \ge 1}$, $(a_n)_{n \ge 1}$ be three non-negative sequences verifying for $n \ge n_0$

$$ r_n \le (1 - \gamma_n) r_{n-1} + \gamma_n a_n (s_n + r_n), $$

where $\gamma_n := c_\gamma n^{-\gamma}$ with $\gamma \in (1/2, 1)$ and $c_\gamma > 0$. Assume that $(a_n)$ is non-increasing and converging to 0, and that $s_n = \mathcal{O}(v_n)$ with $(v_n)$ positive and non-decreasing. Then, for any $0 < \kappa < c_\gamma(1 - 2^{\gamma-1})/(1 - \gamma)$, we have

$$ r_n = \mathcal{O} \left( (\exp(-\kappa n^{1-\gamma}) + a_{\lfloor n/2 \rfloor}) v_n \right). $$

*Proof.* Let $n_1 \ge n_0$ be such that $a_{n_1} \le 1/2$, and $M \ge 0$ such that for any $n \in \mathbb{N}$, $s_n \le Mv_n$. For any $n \ge n_1$:

$$ r_n \le (1 - \gamma_n) r_{n-1} + M \gamma_n a_n v_n + \frac{1}{2} \gamma_n r_n \le \left(1 - \frac{1}{2} \gamma_n\right) r_{n-1} + M \gamma_n a_n v_n. $$

Since $(v_n)$ is positive and non-decreasing, we have for $n \ge n_1$:

$$ \frac{r_n}{v_n} \le \left(1 - \frac{1}{2}\gamma_n\right) \frac{r_{n-1}}{v_n} + M\gamma_n a_n \le \left(1 - \frac{1}{2}\gamma_n\right) \frac{r_{n-1}}{v_{n-1}} + \frac{1}{2}\gamma_n (2Ma_n). $$

We can apply Proposition B.3 to $(r_n/v_n)$, and we obtain for any $\kappa < c_\gamma(1 - 2^{\gamma-1})(1 - \gamma)^{-1}$:

$$ \frac{r_n}{v_n} = \mathcal{O} \left( \exp \left( -\kappa n^{1-\gamma} \right) + a_{\lfloor n/2 \rfloor} \right). $$

$\square$
**Proposition B.3.** Let $(r_n)_{n \ge 0}$, $(\gamma_n)_{n \ge 1}$, and $(a_n)_{n \ge 1}$ be three non-negative sequences verifying for $n \ge n_0$:

$$r_n \le (1 - \gamma_n) r_{n-1} + \gamma_n a_n. \qquad (15)$$

Assume that $(a_n)$ is non increasing. Then, for all $n, m$ with $n_0 < m \le n$, we have the upper bound:

$$r_n \le \exp \left( - \sum_{k=m}^{n} \gamma_k \right) \left( r_{n_0} + \sum_{k=n_0+1}^{m-1} \gamma_k a_k \right) + a_m.$$

In particular, if $\gamma_n = c/(n^\gamma + c')$ with $c > 0$, $\gamma \in (0, 1)$, $c' \ge 0$, and if $\lim_{n \to \infty} a_n = 0$, then for any $\kappa$ such that $0 < \kappa < c(1 - 2^{\gamma-1})(1 - \gamma)^{-1}$, we have:

$$r_n = \mathcal{O} \left( \exp(-\kappa n^{1-\gamma}) + a_{\lceil \frac{n}{2} \rceil} \right)$$

*Proof.* By induction on $n$, and by splitting the sum at any rank $n_0 < m \le n$, one has

$$r_n \le \prod_{j=n_0+1}^{n} (1 - \gamma_j) r_{n_0} + \sum_{k=n_0+1}^{m-1} \left( \prod_{j=k+1}^{n} (1 - \gamma_j) \right) \gamma_k a_k + \sum_{k=m}^{n} \left( \prod_{j=k+1}^{n} (1 - \gamma_j) \right) \gamma_k a_k$$

Using that $1 - x \le e^{-x}$, that $0 \le 1 - \gamma_j \le 1$ and that $(a_n)$ is non increasing, we get

$$r_n \le \exp \left( - \sum_{j=n_0+1}^{n} \gamma_j \right) r_{n_0} + \exp \left( - \sum_{j=m}^{n} \gamma_j \right) \underbrace{\sum_{k=n_0+1}^{m-1} \gamma_k a_k + a_m \sum_{k=m}^{n} \left( \prod_{j=k+1}^{n} (1 - \gamma_j) \right) \gamma_k}_{=:S}$$

By identifying a telescopic sum, we get:

$$S = a_m \sum_{k=m}^{n} \left( \prod_{j=k+1}^{n} (1 - \gamma_j) - \prod_{j=k}^{n} (1 - \gamma_j) \right) = a_m \left( 1 - \prod_{j=m}^{n} (1 - \gamma_j) \right) \le a_m,$$

which proves the first claim.

For the particular case, assume $\gamma_n = c/(n^\gamma + c')$ with $c > 0$, $c' \ge 0$, $\gamma \in (0, 1)$, and $\lim_{n \to \infty} a_n = 0$. Then, $\sum_{k=m}^{n} \gamma_k \sim c(1 - \gamma)^{-1}(n^{1-\gamma} - m^{1-\gamma})$, and choosing $m = \lceil n/2 \rceil$, the sum grows as $\kappa n^{1-\gamma}$ with $\kappa := c(1 - \gamma)^{-1}(1 - 2^{\gamma-1})$ positive. The exponential $\exp(-\sum_{k=m}^{n} \gamma_k)$ decays as $\exp(-\kappa n^{1-\gamma} + \text{const})$. The bound becomes

$$r_n = \mathcal{O} \left( \exp(-\kappa n^{1-\gamma}) (1 + n^{1-\gamma}) + a_{\lceil n/2 \rceil} \right) \subset \mathcal{O} \left( \exp(-\kappa' n^{1-\gamma}) + a_{\lceil n/2 \rceil} \right)$$

for $\kappa' < \kappa$, completing the proof. $\square$

<figure><img src="image_5.png" /></figure>

The following Lemma B.4 is a simple corollary of the Robbins-Siegmund theorem, and is stated and proved in Godichon-Baggioni et al. (2025).

**Lemma B.4.** Let $(V_n)$, $(B_n)$, $(E_n)$, $(D_n)$ and $(a_n)$ be five positive sequences adapted to $(\mathcal{F}_n)_{n \in \mathbb{N}}$ such that

$$\mathbb{E}[V_n | \mathcal{F}_{n-1}] \le (1 + E_n) V_{n-1} + B_n - D_n \quad \text{a.s.}$$

Assume also that $\sum_{n=0}^{\infty} \frac{B_n}{a_n} < +\infty$ a.s. If $a_n \to \infty$, then $V_n = o(a_n)$ a.s.
## C Assumptions for the streaming case

The objective of this section is to give the assumptions for the Corollary 5.1, directly stated on the function $f$. More precisely, in the case where $F(\theta) = E_{\xi}[f(\xi, \theta)]$ and considering the i.i.d. or mini-batch setting, the assumptions on the oracles are satisfied as soon as corresponding assumptions on the function $f$ are satisfied:

(B1) (Differentiability of $f$). For all $\theta \in \mathbb{R}^d$, the random function $\theta \mapsto f(\xi, \theta)$ is almost surely twice differentiable at $\theta$.

(B2) (Growth Condition). There exists $\mathcal{L}_g, \mathcal{L}_h, \sigma^2 \ge 0$ such that for all $\theta \in \mathbb{R}^d$:

$$ \text{a)} \quad \mathbb{E}[\|\nabla f(\xi, \theta)\|^2] \le 2\mathcal{L}_g(F(\theta) - F(\theta^*)) + \sigma^2 $$

$$ \text{b)} \quad \mathbb{E} \left[ \| \nabla^2 f(\xi, \theta) \|_{op}^2 \right] \le \mathcal{L}_h. $$

(B3) (Hessian at minimizer). The Hessian matrix $H := \nabla^2 F(\theta^*)$ at the minimizer is positive definite, and the mapping $\theta \mapsto \nabla^2 F(\theta)$ is continuous at $\theta^*$.

(B4) (Lyapunov Conditions). There exist $q, q' > 2$ and $M, M' > 0$ such that:

$$ \text{a)} \quad \sup_{\|\theta - \theta^*\| < M} \mathbb{E} [\|\nabla f(\xi, \theta) - \nabla F(\theta)\|^q] < +\infty $$

$$ \text{b)} \quad \sup_{\|\theta - \theta^*\| < M'} \mathbb{E} \left[ \| \nabla^2 f(\xi, \theta) - \nabla^2 F(\theta) \|_F^{q'} \right] < +\infty. $$

(B5) (Continuity of Covariance). The covariance function $\theta \mapsto \text{Cov}(\nabla f(\xi, \theta))$ is continuous at $\theta^*$, and we define $\Sigma := \text{Cov}(\nabla f(\xi, \theta^*))$, which, as $\nabla F(\theta^*) = 0$, is equivalent to $\Sigma = \mathbb{E}[\nabla f(\xi, \theta^*)\nabla f(\xi, \theta^*)^T]$.

## D Matrix gradient

### Gradient and adjoint in a Hilbert space

Let $(\mathcal{H}, \langle \cdot, \cdot \rangle)$ be a Hilbert space, $\| \cdot \|$ the norm induced by the inner product, $U \subset \mathcal{H}$ an open set, and $f: U \to \mathbb{R}$ a function differentiable at a point $a \in U$. The differential $df(a)$ being by definition a continuous linear operator from $\mathcal{H}$ to $\mathbb{R}$, the Riesz representation theorem ensures the existence of a unique element $\nabla f(a) \in \mathcal{H}$ such that

$$ \forall h \in \mathcal{H}, \quad df(a)(h) = \langle \nabla f(a), h \rangle. $$

The element $\nabla f(a)$ is called the gradient of $f$ at $a$.

Let $L: \mathcal{H} \to \mathcal{H}$ be a continuous linear operator. For any $y \in \mathcal{H}$, the mapping $h \mapsto \langle y, Lh \rangle$ is a bounded linear functional on $\mathcal{H}$, so by the Riesz representation theorem, there exists a unique element $L^*y \in \mathcal{H}$ such that

$$ \forall h \in \mathcal{H}, \quad \langle y, Lh \rangle = \langle L^*y, h \rangle. $$

The operator $L^*: \mathcal{H} \to \mathcal{H}$ is called the adjoint operator of $L$.
## Gradient of a quadratic functional

Let $L: \mathcal{H} \to \mathcal{H}$ be a continuous linear operator, and $b \in \mathcal{H}$. Consider the quadratic functional $f: \mathcal{H} \to \mathbb{R}$ defined as

$$f(x) := \|Lx + b\|^2.$$

Let $x, h \in \mathcal{H}$. We have

$$\begin{aligned}
f(x+h) &= \|L(x+h) + b\|^2 \\
&= \|Lx + b\|^2 + 2\langle Lx + b, Lh \rangle + \|Lh\|^2 \\
&= f(x) + 2\langle L^*(Lx + b), h \rangle + \|Lh\|^2.
\end{aligned}$$

The gradient of $f$ at $x$ is therefore given by

$$\nabla f(x) = 2L^*(Lx + b).$$

## Application to matrix spaces

We now apply the previous result to the Hilbert space $\mathcal{H} = \mathcal{M}_d(\mathbb{R})$, the space of real $d \times d$ matrices, endowed with the Frobenius inner product:

$$\langle M, N \rangle_F := \text{tr}(M^T N).$$

The space $\mathcal{H}$ is finite-dimensional, so all linear operators on $\mathcal{H}$ are continuous. Let $A, B \in \mathcal{H}$. The adjoint of the linear operator $M \mapsto AM$ is $N \mapsto A^T N$. Indeed, for any $M, N \in \mathcal{H}$,

$$\langle N, AM \rangle_F = \text{tr}(N^T AM) = \text{tr}((A^T N)^T M) = \langle A^T N, M \rangle_F.$$

Consider the quadratic functional $f: \mathcal{H} \to \mathbb{R}$ defined as

$$f(M) := \|AM + B\|_F^2.$$

The gradient of $f$ at point $M \in \mathcal{H}$ is given by

$$\nabla f(M) = 2A^T (AM + B).$$

## Gradient of composition with transpose

Let $f: \mathcal{H} \to \mathbb{R}$ be a differentiable function, and define $g: \mathcal{H} \to \mathbb{R}$ as

$$g(M) := f(M^T).$$

Denoting by $t: \mathcal{H} \to \mathcal{H}$ the transpose operator for matrices, we have $g = f \circ t$. The transpose operator is linear and self-adjoint since for any $M, N \in \mathcal{H}$,

$$\langle N, t(M) \rangle_F = \text{tr}(N^T M^T) = \text{tr}((N^T)^T M) = \langle t(N), M \rangle_F.$$
Using the chain rule for differentiation, we obtain that the differential of $g$ at point $M \in \mathcal{H}$ is given by

$$dg(M) = df(M^T) \circ t.$$

Therefore, the gradient of $g$ at point $M \in \mathcal{H}$ is

$$\nabla g(M) = t^*(\nabla f(M^T)) = t(\nabla f(M^T)) = (\nabla f(M^T))^T.$$