Stochastic Newton Proximal Extragradient Method

Authors: ruichen  jiang
michał  derezi
aryan  mokhtari
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 15:56:29

Abstract: Stochastic second-order methods are known to achieve fast local convergence in strongly convex optimization by relying on noisy Hessian estimates to precondition the gradient. Yet, most of these methods achieve superlinear convergence only when the stochastic Hessian noise diminishes, requiring an increase in the per-iteration cost as time progresses. Recent work in \cite{na2022hessian} addressed this issue via a Hessian averaging scheme that achieves a superlinear convergence rate without increasing the per-iteration cost. However, the considered method exhibits a slow global convergence rate, requiring up to $\tilde{\mathcal{O}}(\kappa^2)$ iterations to reach the superlinear rate of $\tilde{\mathcal{O}}((1/t)^{t/2})$, where $\kappa$ is the problem's condition number. In this paper, we propose a novel stochastic Newton proximal extragradient method that significantly improves these bounds, achieving a faster global linear rate and reaching the same fast superlinear rate in  $\tilde{\mathcal{O}}(\kappa)$ iterations. We achieve this by developing a novel extension of the Hybrid Proximal Extragradient (HPE) framework, which simultaneously achieves fast global and local convergence rates for strongly convex functions with access to a noisy Hessian oracle.

Main paper: 
 introduction: 
 In this paper, we focus on the use of second-order methods for solving the optimization problemEQUATIONwhere f : R d → R is strongly convex and twice differentiable. There is an extensive literature on second-order methods and their fast local convergence properties; e.g., [2] [3] [4] [5] . However, these results necessitate access to the exact Hessian, which can pose computational challenges. To address this issue, several studies have explored scenarios where only the exact gradient can be queried, while a stochastic estimate of the Hessian is available-similar to the setting we investigate in this paper. This oracle model is commonly encountered in large-scale machine learning problems, as computing the gradient is often much less expensive than computing the Hessian, and approximating the Hessian is a more affordable approach. Specifically, consider a finite-sum minimization problem min x∈R d n i=1 f i (x), where n denotes the number of data points and d denotes the dimension of the problem. To achieve a fast convergence rate, standard first-order methods need to compute one full gradient in each iteration, resulting in a per-iteration computational cost of O(nd). In contrast, implementing a second-order method such as damped Newton's method involves computing the full Hessian, which costs O(nd 2 ). An inexact Hessian estimate can be constructed efficiently at a cost of O(sd 2 ), where s is the sketch size or subsampling size [1, 6] . Hence, when the number of samples n significantly exceeds d, the per-iteration cost of stochastic second-order methods becomes comparable to that of first-order methods. Moreover, using second-order information often reduces the number of iterations needed to converge, thereby lowering overall computational complexity.A common template among stochastic second-order methods is to combine a deterministic secondorder method, such as Newton's method or cubic regularized Newton method, with techniques such as Hessian subsampling [7] [8] [9] [10] [11] [12] or Hessian sketching [4, 6, 13] that only require a noisy estimate Table 1 : Comparison between Algorithm 1 and the stochastic Newton method in [1] , in terms of how many iterations it takes to transition to each phase, and the convergence rates achieved. We drop constant factors as well as logarithmic dependence and 1/δ, and assume 1/poly(κ) ≤ Υ ≤ O(κ).Weights Linear phase Initial superlinear phase Final superlinear phase T1 rate ϕ T2 rate θ(1) t T3 rate θ(2) tNewton [1] Uniform Υ 2 1 -κ -2 κ 3κ 3 t κ 6 Υ 2 Υ log(t) t Non-uniform Υ 2 1 -κ -2 κ 2 κ 4 log(κ)+1 t log(t) κ 2 Υ log(t) √ t Stochastic NPE (Ours) Uniform Υ 2 κ 2 1 -κ -1 κ 2 κ 2 t κ 4 Υ 2 Υ log(t) t Non-uniform Υ 2 κ 2 1 -κ -1 Υ 2 + κ (Υ 2 +κ) log(Υ 2 +κ)+1t log(t)Υ 2 + κ Υ log(t) √ t of the Hessian. We refer the reader to [14, 15] for recent surveys and empirical comparisons. In terms of convergence guarantees, the majority of these works, including [4, 6, [8] [9] [10] [11] 13] , have shown that stochastic second-order methods exhibit a global linear convergence and a local linear-quadratic convergence, either with high probability or in expectation. The linear-quadratic behavior holds whenEQUATIONwhere x * denotes the optimal solution of Problem (1) and c 1 , c 2 are constants depending on the sample/sketch size at each step. In particular, the presence of the linear term in (2) implies that the algorithm can only achieve linear convergence when the iterate is sufficiently close to the optimal solution x * . Consequently, as discussed in [9, 10] , to achieve superlinear convergence, the coefficient c 1 = c 1,t needs to gradually decrease to zero as t increases. However, since c 1 is determined by the magnitude of the stochastic noise in the Hessian estimate, this in turn demands the sample/sketch size to increase across the iterations, leading to a blow-up of the per-iteration computational cost.The only prior work addressing this limitation and achieving a superlinear rate for a stochastic second-order method without requiring the stochastic Hessian noise to converge to zero is by [1] . It uses a weighted average of all past Hessian approximations as the current Hessian estimate. This approach reduces stochastic noise variance in the Hessian estimate, though it introduces bias to the Hessian approximation matrix. When combined with Newton's method, it was shown that the proposed method achieves local superlinear convergence with a non-asymptotic rate of (Υ log(t)/t) t with high probability, where Υ characterizes the noise level of the stochastic Hessian oracle (see Assumption 4) . However, the method may require many iterations to achieve superlinear convergence. Specifically, with the uniform averaging scheme, it takes Õ(κ 3 ) iterations before the method starts converging superlinearly and Õ(κ 6 /Υ 2 ) iterations before it reaches the final superlinear rate. Here, κ = L 1 /µ denotes the condition number of the function f , where L 1 is the Lipschitz constant of the gradient and µ is the strong convexity parameter. To address this, [1] proposed a weighted averaging scheme that assigns more weight to recent Hessian estimates, improving both transition points to Õ(Υ 2 + κ 2 ) while achieving a slightly slower superlinear rate of O(Υ log(t)/ √ t).Our contributions. In this paper, we improve the complexity of Stochastic Newton in [1] with a method that attains a superlinear rate in significantly fewer iterations. As shown in Table 1 , our method requires fewer iterations for linear convergence, denoted as T 1 , by a factor of κ 2 compared to [1] . Additionally, our method achieves a linear convergence rate of (1 -O(1/κ)) t , outperforming the (1 -O(1/κ 2 )) t rate in [1] . Thus, our method reaches the local neighborhood of the optimal solution x * and transitions from linear to superlinear convergence faster. Specifically, the second transition point, T 2 , is smaller by a factor of κ in both uniform and non-uniform averaging schemes when Υ = O( √ κ). Similarly, our method's initial superlinear rate has a better dependence on κ, leading to fewer iterations, T 3 , to enter the final superlinear phase. To achieve this result, we use the hybrid proximal extragradient (HPE) framework [16, 17] instead of Newton's method as the base algorithm. The HPE framework provides a principled approach for designing second-order methods with superior global convergence guarantees [3, [17] [18] [19] [20] . However, [16] and subsequent works focus on cases where f is merely convex, not leveraging strong convexity. Thus, we modify the HPE framework to suit our setting. Specifically, we relax the error condition for computing the proximal step in HPE, enabling a larger step size when the iterate is close to the optimal solution, crucial for achieving the final superlinear convergence rate. 
preliminaries: 
 In this section, we formally present our assumptions. Assumption 1. The function f is twice differentiable and µ-strongly convex. Assumption 2. The Hessian ∇ 2 f satisfies ∥∇ 2 f (x) -∇ 2 f (y)∥ ≤ M 1 . Assumption 3. The Hessian ∇ 2 f is L 2 -Lipschitz, i.e., ∥∇ 2 f (x) -∇ 2 f (y)∥ ≤ L 2 ∥x -y∥ 2 .Assumption 2 is more general than the assumption that ∇f is L 1 -Lipschitz. In particular, if the latter assumption holds, then M 1 ≤ L 1 . Moreover, we define κ ≜ M 1 /µ as the condition number.To simplify our notation, we denote the exact gradient ∇f (x) and the exact Hessian ∇ 2 f (x) of the objective function by g(x) and H(x), respectively. As mentioned earlier, we assume that we have access to the exact gradient, but we only have access to a noisy estimate of the Hessian denoted by Ĥ(x). In fact, we require a mild assumption on the Hessian noise. We define the stochastic Hessian noise as E(x) ≜ Ĥ(x) -H(x), where it is assumed to be mean zero and sub-exponential.Assumption 4. If we define E(x) ≜ Ĥ(x) -H(x), then E[E(x)] = 0 and E[∥E(x)∥ p ] ≤ p!Υ pE /2 for all integers p ≥ 2. Also, define Υ ≜ Υ E /µ to be the relative noise level.Assumption 5. The Hessian approximation matrix is positive semi-definite, i.e., Ĥ(x) ⪰ 0, ∀x ∈ R d .Stochastic Hessian construction. The two most popular approaches to construct stochastic Hessian approximations are "subsampling" and "sketching". Hessian subsampling is designed for a finite-sum objective of the form f(x) = 1 n n i=1 f i (x),where n is the number of samples. In each iteration, a subset S ⊂ {1, 2, . . . , n} is drawn uniformly at random, and then the subsampled Hessian at x is constructed as Ĥ(x) = 1 |S| i∈S ∇ 2 f i (x).In this case, if each f i is convex, then the condition in Assumption 5 is satisfied. Moreover, if we further assume that ∥∇ 2 f i (x)∥ ≤ cM 1 for some c > 0 and for all i, then Assumption 4 is satisfied with Example 1] ). The other approach is Hessian sketching, applicable when the Hessian H can be easily factorized as H = M ⊤ M, where M ∈ R n×d is the square-root Hessian matrix, and n is the number of samples. This is the case for generalized linear models; see [1] . To form the sketched Hessian, we draw a random sketch matrix S ∈ R s×n with sketch size s from a distribution D that satisfies E D [S ⊤ S] = I. The sketched Hessian is then Ĥ = M ⊤ S ⊤ SM. In this case, Assumption 5 is automatically satisfied. Moreover, for Gaussian sketch, Assumption 4 is satisfied with Υ = O(κ( d/s + d/s)) (see [1, Example 2] ). Remark 1. The above assumptions are common in the study of stochastic second-order methods, appearing in works on Subsampled Newton [7] [8] [9] [10] 12] , Newton Sketch [4, 6, 13] , and notably, [1] . The strong convexity requirement is crucial as stochastic second-order methods have a clear advantage over first-order methods like gradient descent when the function is strongly convex. Specifically, stochastic second-order methods attain a superlinear convergence rate, as shown in this paper, which is superior to the linear rate of first-order methods.Υ = O( cκ log(d)/|S| + cκ log(d)/|S|) (see [1, 
stochastic newton proximal extragradient: 
 Our approach involves developing a stochastic Newton-type method grounded in the Hybrid Proximal Extragradient (HPE) framework and its second-order variant. Therefore, before introducing our proposed algorithm, we will provide a brief overview of the core principles of the HPE framework. Following this, we will present our method as it applies to the specific setting addressed in this paper.Hybrid Proximal Extragradient. Next, we first present the Hybrid Proximal Extragradient (HPE) framework for strongly convex functions. To solve problem (1), the HPE algorithm consists of two steps. In the first step, given x t , we find a mid-point xt by applying an inexact proximal point update xt ≈ x t -η t ∇f (x t ), where η t is the step size. More precisely, we requireEQUATIONwhere γ t = 1+2η t µ, µ is the strong convexity parameter, and α ∈ (0, 1) is a user-specified parameter. Then, in the second step, we perform the extra-gradient update and compute x t+1 based onEQUATIONThe weights 1 γt in the above convex combination are chosen to optimize the convergence rate. Remark 2. When µ = 0, the algorithm outline above reduces to the original HPE framework studied in [16, 18] . Our modification in (3) is inspired by [21] and allows a larger error when performing the inexact proximal point update, which turns out to be crucial for achieving a fast superlinear convergence rate. Moreover, the modification in (4) has been adopted in [22] .Stochastic Newton Proximal Extragradient (SNPE). The HPE method described above provides a useful algorithmic framework, instead of a directly implementable method. The main challenge comes from implementing the first step in (3) , which involves an inexact proximal point update. Specifically, the naive approach is to solve the implicit nonlinear equation xx t + η t ∇f (x) = 0, which can be as costly as solving the original problem in (1) . To address this issue, [18] proposed to approximate the gradient operator ∇f (x) by its local linearization ∇f (x t ) + ∇ 2 f (x t )(x -x t ), and then compute xt by solving the linear system of equations xtx t + η t (∇f (x t ) + ∇ 2 f (x t )(x t -x t )) = 0. This leads to the Newton proximal extragradient method that was proposed and analyzed in [18] . However, in our setting, the exact Hessian ∇ 2 f (x t ) is not available. Thus, we construct a stochastic Hessian approximation Ht from our noisy Hessian oracle as a surrogate of ∇ 2 f (x t ). We will elaborate on the construction of Ht later, but for the present discussion assume that this stochastic Hessian approximation Ht is already provided. Then in the first step, we will compute xt byEQUATIONwhere we replace ∇f (x t ) by its local linear approximation ∇f (x t ) + Ht (x t -x t ). Moreover, ( 5) is equivalent to solving the following linear system of equations(I + η t Ht )(x -x t ) = -η t ∇f (x t ).For ease of presentation, we set xt as the exact solution of this system, leading toEQUATIONHowever, we note that an inexact solution to this linear system is also sufficient for our convergence guarantees so long as∥(I + η t Ht )(x t -x t ) + η t ∇f (x t )∥ ≤ α 2 ∥x t -x t ∥;We refer the reader to Appendix A.3 for details. Additionally, since we employed a linear approximation to determine the mid-point xt , the condition in (3) may no longer be satisfied. Consequently, it is crucial to verify the accuracy of our approximation after selecting xt . To achieve this, we implement a line-search scheme to ensure that the step size is not large and the linear approximation error is small. Next, we discuss constructing the stochastic Hessian approximation Ht . A simple strategy is using Ĥ(x t ) instead of ∇ 2 f (x t ), but the Hessian noise would lead to a highly inaccurate approximation of the prox operator, ruining the superlinear convergence rate. To reduce Hessian noise, we follow [1] and use an averaged Hessian estimate H(x t ). We consider two schemes: (i) uniform averaging; (ii) nonuniform averaging with general weights. In the first case, Ht = 1 t+1 t i=0 Ĥ(x t )uniformly averages past stochastic Hessian approximations. Motivated by the central limit theorem for martingale differences, we expect Ht to have smaller variance than Ĥ(x t ). It can be implemented online as Ht = t t+1 Ht-1 + 1 t+1 Ĥ(x t ), without storing past Hessian estimates. However, H(x t ) is a biased estimator of ∇ 2 f (x t ), since it incorporates stale Hessian information. To address the bias-variance trade-off, the second case uses non-uniform averaging to weight recent Hessian estimates more. Given an increasing non-negative weight sequence {w t } ∞ t=-1 with w -1 = 0, the running average is:EQUATIONEquivalently, with z i,t = wi-wi-1 wt , Ht can be written as t i=0 z i,t Ĥ(x i ). We discuss uniform averaging in Section 4 and non-uniform averaging in Section 5.Building on the discussion thus far, we are ready to integrate all the components and present our Stochastic Newton Proximal Extragradient (SNPE) method. The steps of SNPE are summarized in Algorithm 1. Each iteration of our SNPE method includes two stages. In the first stage, starting with the current point x t , we first query the noisy Hessian oracle and compute the averaged stochastic Hessian Ht from (7) , as stated in Step 4. Then given the gradient ∇f (x t ), the Hessian approximation Ht , and an initial trial step size σ t , we employ a backtracking line search to obtain η t and xt , as stated in Step 6 of Algorithm 1. Specifically, in this step, we set η t ← σ t and compute xt as suggestedAlgorithm 1 Stochastic NPE 1: Input: x0 ∈ R d , weights {wt} ∞ t=0, line-search parameters α, β ∈ (0, 1), initial step size σ0 > 0 2: Initialize: H-1 = 0 and w-1 = 0 3: for t = 0, 1, . . . do 4:Obtain a stochastic Hessian Ĥt = Ĥ(xt) 5:Compute Ht = w t-1w t Ht-1 + (1 - w t-1 w t ) Ĥt 6:(ηt, xt) = BLS(xt, ∇f (xt), Ht, α, β, σt) 7:Let γt = 1 + 2ηtµ and computext+1 = 1 γ t (xt -ηt∇f (xt)) + (1 -1 γ t )xt 8:Set σt+1 = ηt/β 9: end for Subroutine 1 (η, x) = BLS(x, g, H, α, β, σ)1: Input: current iterate x ∈ R d , gradient g ∈ R d ,Hessian approximation H ∈ R d×d , line-search parameters α, β ∈ (0,1), initial trial step size σ > 0 2: Set η ← σ and x ← x -η(I + η H) -1 g 3: Set γ ← 1 + 2ηµ 4: while ∥x -x + η∇f (x)∥ > α √ γ∥x -x∥ do 5:Set η ← βη and x ← x -η(I + η H) -1 g 6:Set γ ← 1 + 2ηµ Remark 3. Similar to the analysis in [22] , we can show that the total number of line search steps after t iterations can be bounded by 2t -1 + log( σ0 ηt-1 ). Moreover, when t is large enough, on average the line search requires 2 steps per iteration. We defer the details to Appendix A.4. Remark 4. Our motivation behind the choice σ t+1 = η t /β is to allow the step size to grow, which is necessary for achieving a superlinear convergence rate. Specifically, as shown in Proposition 1 below, we require the step size η t to go to infinity to ensure that lim t→∞ ∥xt+1-x * ∥ ∥xt-x * ∥ = 0. Note that this would not be possible if we simply set σ t+1 = η t , since it would automatically result in η t+1 ≤ σ t+1 ≤ η t . Moreover, this condition σ t+1 = η t /β is explicitly utilized in Lemmas 8 and 16 in the Appendix, where we demonstrate that η t can be lower bounded by the minimum of σ 0 /β t and another term. We should also note that this more aggressive choice of the initial step size at each round could potentially increase the number of backtracking steps. However, as mentioned above, this does not cause a significant issue, since the average number of backtracking steps per iteration can be bounded by a constant close to 2.This section outlines key properties of SNPE, applied in Sections 4 and 5 to determine its convergence rates. The first result reveals the connection between SNPE's convergence rate and the step size η t . Proposition 1. Let {x t } t≥0 and {x t } t≥0 be the iterates generated by Algorithm 1. Then for any t ≥ 0, we have∥x t+1 -x * ∥ 2 ≤ ∥x t -x * ∥ 2 (1 + 2η t µ) -1 .Proposition 1 guarantees that the distance to the optimal solution is monotonically decreasing, and it shows a larger step size implies faster convergence. Hence, we need to provide an explicit lower bound on the step size. This task is accomplished in the next lemma. For ease of notation, we let B be the set of iteration indices where the line search scheme backtracks, i.e., B ≜ {t : η t < σ t }. Moreover, we use g(x) and H(x) to denote the gradient ∇f (x) and the Hessian ∇ 2 f (x), respectively. Lemma 1. For t / ∈ B, we have η t = σ t . For t ∈ B, let ηt = η t /β and xt = x t -ηt (I + ηt Ht ) -1 ∇f (x t ). Then, ∥x t -x t ∥ ≤ 1 β ∥x t -x t ∥. Moreover,η t ≥ max αβ∥x t -x t ∥ ∥g(x t ) -g(x t ) -Ht (x t -x t )∥ , 2α 2 βµ∥x t -x t ∥ 2 ∥g(x t ) -g(x t ) -Ht (x t -x t )∥ 2 .As Lemma 1 demonstrates, in the first case where t / ∈ B, we have η t = σ t . Moreover, since we set σ t = η t-1 /β for t ≥ 1, in this case the step size will increase by a factor of 1/β. In the second case that t ∈ B, our lower bound on the step size η t depends inversely on the normalized approximation errorE t = ∥g(xt)-g(xt)-Ht( xt-xt)∥ ∥xt-xt∥. Also, note that E t involves an auxiliary iterate xt instead of the actual iterate xt accepted by our line search. We use the first result to relate ∥x t -x t ∥ to ∥x -x t ∥.To shed light on our analysis, we use the triangle inequality and decompose this error into two terms:EQUATIONThe first term in (8) represents the intrinsic error from the linear approximation in the inexact proximal update, while the second term arises from the Hessian approximation error. Using the smoothness properties of f , we can upper bound the first term, as shown in the following lemma. Lemma 2. Under Assumptions 2 and 3, we haveEQUATION)Lemma 2 shows that the linear approximation error is upper bounded by M 1 . Moreover, the second upper bound is O(∥x t -x * ∥). Thus, as Algorithm 1 converges to the optimal solution x * , the second bound in (9) will become tighter than the first one, and the right hand side approaches zero.To analyze the second term in (8), we isolate the noise component in our averaged Hessian estimate. Specifically, recall Ht = t i=0 z i,t Ĥi and Ĥi = H i + E i . Thus, we have Ht = Ht + Ēt , where Ht = t i=0 z i,t H i is the aggregated Hessian and Ēt = t i=0 z i,t E i is the aggregated Hessian noise, and it follows from the triangle inequality that ∥H t -Ht ∥ ≤ ∥H t -Ht ∥ + ∥ Ēt ∥. We refer to the first part, ∥H t -Ht ∥, as the bias of our Hessian estimate, and the second part, ∥ Ēt ∥, as the averaged stochastic error. There is an intrinsic trade-off between the two error terms. For the fastest error concentration, we assign equal weights to all past stochastic Hessian noises, i.e., z i,t = 1/(t + 1) for all 0 ≤ i ≤ t, corresponding to the uniform averaging scheme discussed in Section 4. To eliminate bias, we assign all weights to the most recent Hessian matrix H t , i.e., z t,t = 1 and z i,t = 0 for all i < t, but this incurs a large stochastic error. To balance these, we present a weighted averaging scheme in Section 5, gradually assigning more weight to recent stochastic Hessian approximations. 
analysis of uniform hessian averaging: 
 In this section, we present the convergence analysis of the uniform Hessian averaging scheme, where w t = t + 1. In this case, we have Ht = 1 t+1 t i=0 Ĥi . As discussed in Section 3.1, our main task is to lower bound the step size η t , which requires us to control the approximation error E t by analyzing the two error terms in (8) . The first term is bounded by Lemma 2, and the second term can be bounded as ∥H t -Ht ∥ ≤ ∥H t -Ht ∥ + ∥ Ēt ∥. Next, we establish a bound on ∥ Ēt ∥, referred to as the Averaged Stochastic Error, and a bound on ∥H t -Ht ∥, referred to as the Bias Term.Averaged stochastic error. To control the averaged Hessian noise ∥ Ēt ∥, we rely on the concentration of sub-exponential martingale difference, as shown in [1] . Lemma 3 ([1, Lemma 2] ). Let δ ∈ (0, 1) with d/δ ≥ e. Then with probability 1 -δπ 2 /6, we have∥ Ēt ∥ ≤ 8Υ E log(d(t+1)/δ) t+1for any t ≥ 4 log(d/δ).Lemma 3 shows that, with high probability, the norm of averaged Hessian noise ∥ Ēt ∥ approaches zero at the rate of Õ(Υ E / √ t). As discussed in Section 4.1, this error eventually becomes the dominant factor in the approximation error E t and determines the final superlinear rate of our algorithm. Remark 5. Our subsequent results are conditioned on the event that the bound on ∥ Ēt ∥ stated in Lemma 3 is satisfied for all t ≥ 4 log(d/δ). Thus, to avoid redundancy, we will omit the "with high probability" qualification in the following discussion.Bias. We proceed to establish an upper bound on ∥H t -Ht ∥. The proof can be found in Appendix B.1.Lemma 4. If Ht = 1 t+1 t i=0 Ĥi , then ∥H t -Ht ∥ ≤ 1 t+1 t i=0 ∥H t -H i ∥. Moreover, for any i ≥ 0, we have ∥H t -H i ∥ ≤ max{M 1 , 2L 2 ∥x i -x * ∥}.The analysis of the bias term is more complicated. Specifically, to obtain the best result, we break the sum in Lemma 4 into two parts,1 t I-1 i=0 ∥H t -H i ∥ and 1 t t i=I ∥H t -H i ∥,where I is an integer to be specified later. The first part corresponds to the bias from stale Hessian information and converges to zero at O(M 1 I/t), as shown by the first bound in Lemma 4. The second part is the bias from recent Hessian information when the iterates are near the optimal solution x * . Using the second bound in Lemma 4, we show this part contributes less to the total bias and is dominated by the first part. Thus, we can conclude that ∥H t -Ht ∥ = O( M1I t+1 ). Based on the previous discussions, it is evident that the terms contributing to the upper bound of E t all converge to zero, albeit at different rates. Furthermore, the linear approximation error and bias term display distinct global and local convergence patterns, depending on the distance ∥x t -x * ∥. Hence, this necessitates a multi-phase convergence analysis, which we undertake in the following section.Similar to [1] , we consider four convergence phases with three transitions points T 1 , T 2 , and T 3 , whose expressions will be specified later. Due to space limitations, in the following we provide an overview of the four phases and relegate the details to Appendix B.Warm-up phase 0 ≤ t < T 1 . At the beginning of the algorithm, the averaged Hessian estimate is dominated by stochastic noise and provides little useful information for convergence. Thus, there are generally no guarantees on the convergence rate for 0 ≤ t < T 1 . However, due to the line search scheme, Proposition 1 ensures that the distance to x * is non-increasing, i.e., ∥x t+1 -x * ∥ ≤ ∥x t -x * ∥ for all t ≥ 0. During the warm-up phase, the averaged Hessian noise ∥ Ēt ∥, which contributes most to the approximation error E t , is gradually suppressed. Once the averaged Hessian noise is sufficiently concentrated, Algorithm 1 transitions to the linear convergence phase, denoted by T 1 .Linear convergence phase T 1 ≤ t < T 2 . After T 1 iterations, Algorithm 1 starts converging linearly to the optimal solution x * . Moreover, during this phase, all the three errors discussed in Section 3.1 continue to decrease. Specifically, Lemma 2 shows the linear approximation error is bounded by O(∥x t -x * ∥), which converges to zero at a linear rate. Furthermore, Lemma 3 implies that the averaged Hessian error ∥ Ēt ∥ diminishes at a rate of Õ( Υ E √ t ). Finally, regarding the bias term, it can be shown ∥H t -Ht ∥ = O( M1I t ) following the discussions after Lemma 4. Thus, once all the three errors are sufficiently small, Algorithm 1 moves to the superlinear phase, denoted by T 2 .Superlinear phases T 2 ≤ t < T 3 and T 3 ≤ t < T 4 . After T 2 iterations, Algorithm 1 converges at a superlinear rate. Moreover, the superlinear rate is determined by the averaged noise ∥ Ēt ∥, which decays at the rate of Õ( Υ E √ t ), and the bias of our averaged Hessian estimate Ht , which decays at the rate of O( M1I t ). Hence, as the number of iterations t increases, the averaged noise will dominate and the algorithm transitions from the initial superlinear rate to the final superlinear rate.We summarize our convergence guarantees in the following theorem and the proofs are in Appendix B.Theorem 1. Suppose Assumptions 1-5 hold and the weights for Hessian averaging in SNPE are uniform, and define C :=1 2β √ 1-α 2 + 5.Then, the followings hold: (a) Warm-up phase:If 0 ≤ t < T 1 , then ∥x t+1 -x * ∥ ≤ ∥x t -x * ∥, where T 1 = Õ( Υ 2 κ 2 ). (b) Linear convergence phase: If T 1 ≤ t < T 2 , then ∥x t+1 -x * ∥ 2 ≤ ∥x t -x * ∥ 2 (1 + 2αβ 3κ ) -1 , where T 2 = Õ(max{ Υ 2 κ + κ 2 , Υ 2 }) = Õ(Υ 2 + κ 2 ). (c) Initial superlinear phase: For T 2 ≤ t < T 3 , we have ∥x t+1 -x * ∥ ≤ Cρ (1) t ∥x t -x * ∥, where ρ (1) t = 6κI α √ 2β(t+1) = Õ Υ 2 /κ+κ 2 t with I defined in (28) and T 3 = Õ( (Υ 2 /κ+κ 2 ) 2 Υ 2).(d) Final superlinear phase: Finally, for t ≥ T 3 , we have ∥x t+1 -x * ∥ ≤ Cρ t ∥x t -x * ∥, where ρ (2) t = 8 √ 2Υ α √ β log(d(t+1)/δ) t+1 = O Υ log(t) t .Comparison with [1] . As shown in Table 1 , our method in Algorithm 1 with uniform averaging achieves the same final superlinear convergence rate as the stochastic Newton method in [1] . However, it transitions to the linear and superlinear phases much earlier. Specifically, the initial transition point T 1 is improved by a factor of κ 2 , and our linear rate in Lemma 6 is faster. This reduces the iterations needed to reach the local neighborhood, cutting the time to reach T 2 and T 3 by factors of κ and κ 2 . 
analysis of weighted hessian averaging: 
 Previously, we showed Algorithm 1 with uniform averaging eventually achieves superlinear convergence. However, as per Theorem 1, it requires Õ( κ 4 Υ 2 ) iterations to reach this rate. To achieve a faster transition, we follow [1] and use Hessian averaging with a general weight sequence {w t }. We show this method also outperforms the stochastic Newton method in [1] . Specifically, we set w t = w(t) for all integer t ≥ 0, where w(•) : R → R satisfies certain regularity conditions as in [1, Assumption 3] . Assumption 6. (i) w(•) is twice differentiable; (ii) w(-1) = 0, w(t) > 0, ∀t ≥ 0; (iii) w ′ (-1) ≥ 0; (iv) w ′′ (t) ≥ 0, ∀t ≥ -1; (v) max w(t+1) w(t) , w ′ (t+1) w ′ (t) ≤ Ψ, ∀t ≥ 0 for some Ψ ≥ 1.Choosing w(t) = t p for any p ≥ 1 satisfies Assumption 6. Additionally, as discussed in [1] , a suitable choice is w(t) = (t + 1) log(t+4) , allowing us to achieve the optimal transition to the superlinear rate. Since the analysis in this section closely resembles that in Section 4 on uniform averaging, we will only present the final result here for brevity. The four stages of convergence are detailed in the following theorem, with intermediate lemmas and proofs in the appendix. To simplify our bounds, we report results for non-uniform averaging with w(t) = (t + 1) log(t+4) . Theorem 2. Suppose Assumptions 1-5 hold and the weights for Hessian averaging in SNPE are defined as w(t) = (t + 1) log(t+4) , and defineC ′ := ( 1 10β √ 2(1-α 2 ) + 1 √ 2 ). Then, the following hold: (a) Warm-up phase: If 0 ≤ t < U 1 , then ∥x t+1 -x * ∥ ≤ ∥x t -x * ∥, where U 1 = Õ( Υ 2 κ 2 ). (b) Linear convergence phase: If U 1 ≤ t < U 2 , then ∥x t+1 -x * ∥ 2 ≤ ∥x t -x * ∥ 2 (1+ 2αβ 3κ ) -1 , where U 2 = Õ(max{ Υ 2 κ 2 + κ, Υ 2 }) = Õ(Υ 2 + κ). (c) Initial superlinear phase: If U 2 ≤ t < U 3 , then ∥x t+1 -x * ∥ ≤ C ′ θ (1) t ∥x t -x * ∥, where θ (1) t = 5κw(J ) α √ 2βw(t) = Õ κ(Υ 2 +κ) log(Υ 2 +κ)/t log t with J defined in (49) and U 3 = Õ(Υ 2 +κ).(d) Final superlinear phase:Finally, if t ≥ U 3 , then ∥x t+1 -x * ∥ ≤ C ′ θ (2) t ∥x t -x * ∥, where θ (2) t = 8 √ 2Υ α √ β w ′ (t) log(d t+1 δ ) w(t) = O Υ log(t) √ t .In the weighted averaging case, similar to the uniform averaging scenario, we observe four distinct phases of convergence. The warm-up phase for SNPE, during which the distance to the optimal solution does not increase, has the same duration as in the uniform averaging case but is shorter than the warm-up phase for the stochastic Newton method in [1] by a factor of 1/κ 2 . The linear convergence rates of both uniform and weighted Hessian averaging methods are 1 -κ -1 , improving over the 1 -κ -2 rate achieved by the stochastic Newton method in [1] . The number of iterations to reach the initial superlinear phase is Õ(Υ 2 +κ), smaller than the Õ(κ 2 ) needed for uniform averaging in SNPE when we focus on the regime where Υ = O( √ κ). The non-uniform averaging method in [1] requires κ 2 iterations to achieve the initial superlinear phase, whereas the non-uniform SNPE achieves an initial superlinear rate of O(κ log(κ)+1 /t) t , improving over the rate of O(κ 4 log(κ)+1 /t log(t) ) t in [1] . Finally, while the ultimate superlinear rates in all cases are comparable at approximately Õ((1/ √ t) t ), the non-uniform version of SNPE requires Õ(Υ 2 + κ) iterations to attain this fast rate, whereas [1] 's non-uniform version requires Õ(κ 2 ) iterations, which is less favorable. Remark 6. As discussed in [1, Example 1] , with a subsampling size of s, we have Υ = Õ(κ/s) for subsampled Newton. This implies that when s = Ω( √ κ), we achieve Υ = O( √ κ).Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: The code and data are attached in the supplementary material.Guidelines:• The answer NA means that paper does not include experiments requiring code. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
discussion and complexity comparison: 
 In this section, we compare the complexity of our method with accelerated gradient descent (AGD), damped Newton, and stochastic Newton [1] methods. The iteration complexities of these methods are summarized in Table 2 , which we use to establish their overall complexities. Note that since both stochastic Newton and our proposed SNPE method achieve superlinear convergence, their Iteration O( √ κ log(1/ϵ)) O(κ 2 +log log(1/ϵ)) O(Υ 2 +κ 2 + log(1/ϵ) log log(1/ϵ) ) O(Υ 2 +κ+ log(1/ϵ) log log(1/ϵ) )complexities depend on the target accuracy ϵ in the form O( log(1/ϵ) log log(1/ϵ) ), which is provably better than the complexity of AGD by at least a factor of log log(ϵ -1 ). Further details are provided in Appendix D.1. To better compare them, we focus on the finite-sum problem with n functions:min x∈R d n i=1 f i (x).Let ϵ be the target accuracy, κ the condition number, and Υ the noise level in Assumption 4. In this case, computing the exact gradient and Hessian costs O(nd) and O(nd 2 ), respectively. Thus, the per-iteration cost for AGD is O(nd). Each iteration of damped Newton's method requires computing the full Hessian and solving a linear system, resulting in a total periteration cost of O(nd 2 + d 3 ). For both stochastic Newton in [1] and our SNPE method, the per-iteration cost depends on how the stochastic Hessian is constructed. For example, Subsampled Newton constructs the Hessian estimate with a cost of O(sd 2 ), where s denotes the sample size. Newton Sketch has a similar computation cost (see [6] ). Additionally, it takes O(nd) to compute the full gradient and O(d 3 ) to solve the linear system. Since the sample/sketch size s is typically chosen as s = O(d), the total per-iteration cost is O(nd + d 3 ).• Compared to AGD, SNPE achieves better iteration complexity. Specifically, when the noise level Υ and target accuracy ϵ are relatively small (Υ = O( √ κ) and log 1 ϵ = Ω( √ κ)), SNPE converges in fewer iterations. Additionally, when n ≥ d 2 (indicating many samples), the per-iteration costs of both methods are O(nd), giving our method a better overall complexity.• Compared to damped Newton's method, our method's iteration complexity depends better on the condition number κ, while damped Newton's depends better on ϵ. However, when n ≥ d 2 , the per-iteration cost of damped Newton is O(nd 2 ), significantly more than our method's O(nd).• Compared to the stochastic Newton method in [1] , the per-iteration costs of both methods are similar. However, Table 2 shows that our iteration complexity is strictly better. Specifically, when the noise level is relatively small compared to the condition number, i.e., Υ = O( √ κ), the complexity of SNPE improves by an additional factor of κ over the stochastic Newton method.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: We explained how we performed the experiments and included the implementation details in Section 7.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
numerical experiments: 
 While our focus is on improving theoretical guarantees, we also provide simple experiments to showcase our method's improvements. All simulations are implemented on a Windows PC with an AMD processor and 16GB Memory. We consider minimizing the regularized log-sum-exp objective f (x) = ρ log( n i=1 exp(a ⊤ i x-bi ρ)) + λ 2 ∥x∥ 2 , a common test function for second-order methods [23, 24] due to its high ill-conditioning. Here, λ > 0 is a regularization parameter, ρ > 0 is a smoothing parameter, and the entries of the vectors a 1 , . . . , a n ∈ R d and b ∈ R d are randomly generated from the standard normal distribution and the uniform distribution over [0, 1], respectively. In our experiments, the regularization parameter λ is 10 -3 , the dimension d is 500, and the number of samples n is chosen from 50,000, 10,000, and 150,000, respectively. We compare our SNPE method with the stochastic Newton method in [1] , using both uniform Hessian averaging (Section 4) and weighted averaging (Section 5). In addition, we evaluate it against accelerated gradient descent (AGD), damped Newton's method, and Newton Proximal Extragradient (NPE), which corresponds to our SNPE method with the exact Hessian. For the stochastic Hessian estimate, we use a subsampling strategy with a subsampling size of s = 500. Empirically, we found that the extragradient step (Line 7 in Algorithm 1) tends to slow down convergence. To address this, we consider a variant of our method without the extragradient step, modifying Line 7 to x k+1 = xk . In Figure 3 of the Appendix, we further explore the effect of the extragradient step. We also note that similar observations were made in [19] , where the simple iteration x t+1 = x t -η t (I + η t H t ) -1 g t outperformed "accelerated" second-order methods. From Figures 1 and 2 , we have the following observations: Comparison with stochastic Newton. In all cases, our SNPE method outperforms stochastic Newton in both the number of iterations and runtime, due to the problem's highly ill-conditioned nature and our method's better dependence on the condition number.Comparison with AGD. From Figure 1 , we observe that our SNPE method, with either uniform or weighted averaging, requires far fewer iterations to converge than AGD due to the use of second-order information. Consequently, while SNPE has a higher per-iteration cost than AGD, it converges faster overall in terms of runtime, as demonstrated in Figure 2 .Comparison with the damped Newton's method and NPE. As expected, since both damped Newton and NPE use exact Hessian, Figure 1 shows that they exhibit superlinear convergence and converge in fewer iterations than the other algorithms. However, since the exact Hessian matrix is expensive to compute, they incur a high per-iteration computational cost and overall take more time than our proposed SNPE method to converge (see Figure 2 ). Moreover, the gap between these two methods and SNPE widens as the number of samples n increases, demonstrating the advantage of our method in the large data regime.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [No] Justification: Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
conclusions and limitations: 
 We introduced a stochastic variant of the Newton Proximal Extragradient method (SNPE) for minimizing a strongly convex and smooth function with access to a noisy Hessian. Our contributions include establishing convergence guarantees under two Hessian averaging schemes: uniform and non-uniform. We characterized the computational complexity in both cases and demonstrated that SNPE outperforms the best-known results for the considered problem. A limitation of our theory is the assumption of strong convexity. Extending the theory to the convex setting would make it more general. This extension is left for future work due to space limitations.Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: The compute resources are reported in Section 7.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: The research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Answer: [NA] Justification: Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: A.1 Proof of Proposition 1In this section, we prove Proposition 1. In fact, using the same proof, we can also show an additional result that upper bounds ∥x t -xt ∥, which will useful in the proof of Lemma 2. Therefore, we present the full version below for completeness. Proposition 2 (Full version of Proposition 1). Let {x t } t≥0 and {x t } t≥0 be the iterates generated by Algorithm 1. Then for any t ≥ 0, we have∥x k+1 -x * ∥ 2 ≤ ∥x k -x * ∥ 2 (1 + 2η k µ) -1 and ∥x t -xt ∥ ≤ 1 √ 1-α 2 ∥x t -x * ∥.Proof. Our proof is inspired by the approach in [22, Proposition 1] . For any x ∈ R d , we first writeEQUATION)To begin with, we bound the first term in (10) byEQUATION)where the first inequality is due to Cauchy-Schwarz inequality, the second inequality is due to the condition in (3), and the last inequality is due to Young's inequality. Moreover, for the second term in (10) , we use the three-point equality to getEQUATION)By combining (10) , ( 11) and ( 12), we obtain thatEQUATION)Moreover, it follows from the update rule in (4) that η t ∇f (x t ) = x t -x t+1 + 2η t µ(x t -x t+1 ). This implies that, for anyEQUATIONwhere we applied the three-point equality twice in the last equality. Thus, by combining (13) with x = x t+1 and (15) with x = x * , we getEQUATION) Since ∇f (x * ) = 0 and f is µ-strongly convex, we further have ⟨∇f (x t ), xt -x * ⟩ = ⟨∇f (x t ) -∇f (x * ), xt -x * ⟩ ≥ µ∥x t -x * ∥ 2 .(17) Combining ( 16) and (17) and rearranging the terms, we obtain thatEQUATION)Since α < 1, the last term in ( 18) is negative and we immediately obtain that ∥x t+1 -x * ∥ 2 ≤ ∥x t - x * ∥ 2 (1 + 2η t µ) -1 . Moreover, since 1+2ηtµ 2 ∥x t+1 -x * ∥ 2 ≥ 0, it follows that 1-α 2 2 ∥x t -xt ∥ 2 ≤ 1 2 ∥x t -x * ∥ 2 , which leads to ∥x t -xt ∥ ≤ 1 √ 1-α 2 ∥x t -x * ∥.The proof is complete.Recall that in our backtracking line search scheme in Algorithm 1, the step size η t starts from σ t and keeps backtracking until the condition in (3) is satisfied. Hence, by the definition of B, it immediately follows that η t = σ t if t / ∈ B. Moreover, if t ∈ B, then the step size ηt = η t /β and the corresponding iterate xt must have failed the condition in (3) (Otherwise, our line search scheme would have accepted the step size ηt instead). This implies thatEQUATIONSince xt = x t -ηt (I + ηt Ht ) -1 ∇f (x t ), we havext -x t + ηt Ht (x t -x t ) = -η t ∇f (x t ) ⇔ xt -x t = -η t ∇f (x t ) + Ht (x t -x t ) .and hence the left-hand side in (19) equals to ηt ∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥. Thus, we obtain from (19) thatEQUATIONBy substituting η t = β ηt , we further haveη t > αβ 1 + 2η t µ/β∥x t -x t ∥ ∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥ .Using the fact that 1 + 2η t µ/β ≥ 1, we getEQUATIONMoreover, using the fact that 1 + 2η t µ/β ≥ 2η t µ/β, we can also conclude thatη t > αβ 2η t µ/β∥x t -x t ∥ ∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥ ⇒ η t > 2α 2 βµ∥x t -x t ∥ 2 ∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥ 2 .(22) By combining (21) and (22) , we obtain the lower bound in Lemma 1.Finally, when Ht ⪰ 0, it holds that I + ηt Ht ⪰ I + η t Ht ⪰ 0 and thus (I + η t Ht ) -1 ⪰ (I + ηt Ht ) -1 ⪰ 0. This further implies that ∥(I + η t Ht ) -1 ∇f (x t )∥ ≥ ∥(I + ηt Ht ) -1 ∇f (x t )∥. Hence, we can conclude that∥x t -x t ∥ = ηt ∥(I + ηt Ht ) -1 ∇f (x t )∥ ≤ η t β ∥(I + η t Ht ) -1 ∇f (x t )∥ ≤ 1 β ∥x t -x t ∥.This completes the proof.In this section, we extend our convergence results to the case where the linear system in (6) is solved inexactly, i.e., we find xt such thatEQUATIONIn this case, since the proof of Proposition 2 does not rely on the update rule in (6), Proposition 2 continues to hold. However, we need to modify the proof of Lemma 1 and replace it by the following lemma. We note that the two results differ only by an absolute constant. Lemma 5 (Extension to Lemma 1). For t / ∈ B, we have η t = σ t . For t ∈ B, let ηt = η t /β and xt be the corresponding iterate rejected by our line search scheme. Then, ∥x t -x t ∥ ≤ 3 β ∥x t -x t ∥. Moreover,η t ≥ max αβ∥x t -x t ∥ 2∥g(x t ) -g(x t ) -Ht (x t -x t )∥ , α 2 βµ∥x t -x t ∥ 2 2∥g(x t ) -g(x t ) -Ht (x t -x t )∥ 2 .Proof. We follow a similar argument as in the proof of Lemma 1. If t / ∈ B, it immediately follows that η t = σ t . Further, if t ∈ B, then ηt and the corresponding iterate xt must fail to satisfy the condition in (3). This implies thatEQUATIONMoreover, by our inexactness condition in (23) , xt satisfies∥(I + ηt Ht )(x t -x t ) + ηt ∇f (x t )∥ ≤ α 2 ∥x t -x t ∥.Hence, by using triangle inequality, we haveηt ∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥ = (x t -x t + ηt ∇f (x t )) -(I + ηt Ht )(x t -x t ) + ηt ∇f (x t ) ≥ ∥x t -x t + ηt ∇f (x t )∥ -∥(I + ηt Ht )(x t -x t ) + ηt ∇f (x t )∥ ≥ α 2 1 + 2η t µ∥x t -x t ∥.Thus, we obtain thatηt > α √ 1 + 2η t µ∥x t -x t ∥ 2∥∇f (x t ) -∇f (x t ) -Ht (x t -x t )∥ ,Combining this with (20), we observe that these two bounds differ only by a constant factor of 2.Hence, the rest of the proof for the lower bound follows similarly as in Lemma 1.Next, we prove the inequality ∥x t -x t ∥ ≤ 1 β ∥x t -x t ∥. Define x * t = x t -η t (I + η t Ht ) -1 ∇f (x t ) and x * t = x t -ηt (I + ηt Ht ) -1 ∇f (x t ), i.e, they are the exact solutions to the corresponding linear systems. By the argument in the proof of Lemma 1, we have ∥x *t -x t ∥ ≤ 1 β ∥x * t -x t ∥.In the following, we first prove thatEQUATION)It suffices to prove the first set of inequalities, since the second one follows similarly. To see this, note that the condition in ( 23) can be rewritten as∥(I + η t Ht )(x t -x * t )∥ ≤ α 2 ∥x t -x t ∥.Since Ht ⪰ 0, this further implies that ∥x t -x * t ∥ ≤ ∥(I + η t Ht )(x t -x * t )∥ ≤ α 2 ∥x t -x t ∥. Hence, by the triangle inequality, we obtain that∥x * t -x t ∥ ≤ ∥x * t -xt ∥ + ∥x t -x t ∥ ≤ 1 + α 2 ∥x t -x t ∥ ≤ 3 2 ∥x t -x t ∥, ∥x * t -x t ∥ ≥ ∥x * t -xt ∥ -∥x t -x t ∥ ≥ 1 - α 2 ∥x t -x t ∥ ≥ 1 2 ∥x t -x t ∥,which lead to (25) . Finally, we conclude that ∥x t -x t ∥ ≤ 2∥x * t -x t ∥ ≤ 2 β ∥x * t -x t ∥ ≤ 3 β ∥x t -x t ∥.This completes the proof.Let l t denote the number of line search steps in iteration t. We first note that η t = σ t β lt-1 by our line search subroutine, which implies l t = log 1/β (σ t /η t ) + 1. Moreover, recall that σ t = η t-1 /β for t ≥ 1. Hence, the total number of line search steps after t iterations can be bounded by (cf. [22, Lemma 22] ):t-1 i=0 l i = t-1 i=0 log 1/β σ i η i + 1 = 2t -1 + log σ 0 η t-1 .Moreover, it can be shown that η t-1 ≥ αβ/(3M 1 ) when t = Ω(Υ 2 /κ 2 ) in both the uniform averaging and the non-uniform averaging cases (cf. Corollaries 2 and 3). This implies that the total number of line search steps can be bounded by 2t -1 + log(3M 1 σ 0 /αβ).Recall that we use g(x) and H(x) to denote the gradient ∇f (x) and the Hessian ∇ 2 f (x), respectively. We first consider the inequality in (9) . By the fundamental theorem of calculus, we can write∇f (x t ) -∇f (x t ) = 1 0 ∇ 2 f (x t + τ (x t -x t ))(x t -x t ) dτ.Therefore, we can further use the triangle inequality to getEQUATIONMoreover, we have∇ 2 f (x t + τ (x t -x t )) -∇ 2 f (x t ) ≤ M 1 for any τ ∈ [0, 1] by Assumption 2.Together with (26), this further implies that ∥∇f(xt)-∇f (xt)-∇ 2 f (xt)(xt-xt)∥ ∥xt-xt∥≤ M 1 , which proves the first bound in (9) .Next, we consider the second bound in (9) . By Assumption 3, it follows from standard arguments that ∥∇f(x t ) -∇f (x t ) -∇ 2 f (x t )(x t -x t )∥ ≤ L2 2 ∥x t -x t ∥ 2 (e.g., see [25, Lemma 1.2.4]). This leads to ∥∇f (x t ) -∇f (x t ) -∇ 2 f (x t )(x t -x t )∥ ∥x t -x t ∥ ≤ L 2 ∥x t -x t ∥ 2 ≤ L 2 ∥x t -x t ∥ 2β ≤ L 2 ∥x t -x * ∥ 2β √ 1 -α 2 ,where we used ∥x t -x t ∥ ≤ i ∥ ≤ L 2 ∥x t -x i ∥ ≤ L 2 (∥x t -x * ∥ + ∥x i -x * ∥). Since i ≤ t and ∥x t -x * ∥ is non-increasing in t by Proposition 1, we further have ∥x t -x * ∥ ≤ ∥x i -x * ∥, which proves ∥H t -H i ∥ ≤ 2L 2 ∥x i -x * ∥.Warm-up phase. To determine the transition point T 1 , recall that both the linear approximation error and the bias term can be bounded by M 1 according to Lemmas 2 and 4. Since Lemma 3 shows that ∥ Ēt ∥ = Õ(Υ E / √ t), we will have∥ Ēt ∥ ≤ M 1 when t = Ω(Υ 2 E /M 2 1 ) = Ω(Υ 2 /κ 2 ). More specifically, the transition point is given byEQUATIONwhere δ ∈ (0, 1) satisfies d/δ ≥ e, α, β ∈ (0, 1) are line-search parameters, and σ 0 is the initial step size.Linear convergence phase. In the following lemma, we prove the linear convergence of Algorithm 1 with uniform averaging. Lemma 6. Assume that β ≤ 1/2 and recall the definition of T 1 in (27). For any t ≥ T 1 , we have∥x t+1 -x * ∥ 2 ≤ ∥x t -x * ∥ 2 1 + 2αβ 3κ -1, where κ ≜ M 1 /µ is the condition number.Proof. See Appendix B.3. Now we discuss the transition point T 2 . At a high level, the algorithm transitions to the superlinear phase if all three errors discussed in Section 3.1 are reduced from O(M ) to O(µ). For this to happen, first the iterate x t needs to reach a local neighborhood satisfying ∥x t -x * ∥ = O(µ/L 2 ).As a corollary of Lemma 6, this holds at most after an additional Õ(κ) iterations. Specifically, let ν ∈ (0, 1) be a parameter and defineEQUATIONwhere D = ∥x 0 -x * ∥ is the initial distance to the optimal solution. Then Lemma 6 implies that we have ∥x t -x * ∥ ≤ νµ/L 2 for all t ≥ I. Moreover, Lemma 3 implies that the averaged Hessian noise satisfies∥ Ēt ∥ = O(µ) when t = Ω(Υ 2 E /µ 2 ) = Ω(Υ 2). Finally, regarding the bias term, following the discussions after Lemma 4, it can be shown that ∥H t -Ht ∥ = O M1I t+1 . Thus, ∥H t -Ht ∥ = O(µ) when t = Ω(κI). Combining all pieces together, we formally define the second transition point byEQUATIONSinceI = Õ(T 1 +κ) = Õ(Υ 2 /κ 2 +κ), we note that T 2 = Õ(max{Υ 2 , Υ 2 /κ+κ 2 }) = Õ(Υ 2 +κ 2 ).Superlinear phase. In the following theorem, we show that after T 2 iterations, Algorithm 1 with uniform averaging converges at a superlinear rate. See Appendix B.4 for proof.Theorem 3. Let ν ∈ (0, 1) be a parameter satisfying5 2αβ √ (1-α 2 )β + 25 α √ 2βν ≤ 1. and recall the definition of T 2 in (29). Then for any t ≥ T 2 ,∥x t+1 -x * ∥ ≤ 1 2β √ 1 -α 2 + 5 ρ t ∥x t -x * ∥,whereEQUATION)In Theorem 3, we observe that the rate ρ t in (30) goes to zero as the number of iterations t increases, and thus it implies that the iterates converge to x * superlinearly. Moreover, the rate ρ t consists of two terms. The first term comes from the averaged noise ∥ Ēt ∥, which decays at the rate of Õ(Υ/ √ t). In addition, the second term is due to the bias of our averaged Hessian estimate Ht , which decays at the rate of O(κI/t). Hence, when t is sufficiently large, the averaged noise will dominate and the superlinear rate settles for the slower rate of Õ(Υ/ √ t). Specifically, the algorithm transitions from the initial superlinear rate to the final superlinear rate when the two terms in (30) are balanced. Hence, we define the third transition point T 3 as the root of64(T 3 + 1) log(d(T 3 + 1)/δ) = 9κ 2 I 2 Υ 2 . Since I = Õ(Υ 2 /κ 2 + κ), T 3 = Õ((Υ 2 /κ + κ 2 ) 2 /Υ 2 ). We summarize our discussions in the following corollary. Corollary 1. For T 2 ≤ t ≤ T 3 -1, we haveEQUATIONt ∥x t -x * ∥,where ρ(1)t = 6κI α √ 2β(t+1) . Moreover, for t ≥ T 3 , we have ∥x t+1 -x * ∥ ≤ 1 2β √ 1 -α 2 + 5 ρ (2) t ∥x t -x * ∥,where ρ(2)t = 8 √ 2Υ α √ β log(d(t+1)/δ) t+1.We divide the proof of Lemma 6 into the following three steps. First, in Lemma 7, we provide a lower bound on the the step size η t in those iterations where our line search scheme backtracks the step size, i.e., t ∈ B. Building on Lemma 7, we use induction in Lemma 8 to prove a lower bound for all t ≥ 0. This allows us to establish η t = Ω(1/M 1 ) for all t ≥ T 1 in Corollary 2, from which Lemma 6 immediately follows.To simplify the notation, we define the functionEQUATION)Then Lemma 3 can be equivalently written as ∥ Ēt ∥ ≤ ϕ(t) for all t ≥ 4 log(d/δ). We are now ready to state our first lemma. Lemma 7. If t ∈ B, then we have η t ≥ αβ/(2M 1 + ϕ(t)).Proof. If t ∈ B, by Lemma 1 we can lower bound the step size η t byEQUATIONwhereE t ≜ ∥g(xt)-g(xt)-Ht( xt-xt)∥ ∥xt-xt∥is the normalized approximation error. Moreover, as outlined in Section 3.1, we can apply the triangle inequality to upper bound E t . Specifically, we haveEQUATION)By (9) in Lemma 2, the first term in (33) is upper bounded by M 1 . Moreover, it also follows from Lemma 4 that ∥H t -Ht ∥ ≤1 t+1 t i=0 ∥H t -H i ∥ ≤ M 1 .Hence, we further obtainE t ≤ 2M 1 + ∥ Ēt ∥ ≤ 2M 1 + ϕ(t).Combining this with (32), we obtain the desired result.Lemma 7 provides a lower bound on the step size η t , but only for the case where t ∈ B. In the next lemma, we further use induction to show a lower bound for the step sizes in all iterations. Lemma 8. Assume that β ≤ 1 2 . For any t ≥ 0, we haveη t ≥ min αβ 2M 1 + ϕ(t) , σ 0 β t (34)Proof. We prove this lemma by induction. For the base case t = 0, we consider two subcases. If 0 ∈ B, then by Lemma 7 we obtain that η 0 ≥ αβ 2M1+ϕ(0) . Otherwise, if 0 / ∈ B, we have η 0 = σ 0 . In both cases, we observe that (34) is satisfied for the base case t = 0. Now assume that (34) is satisfied for t = s where s ≥ 0. For t = s + 1, we again distinguish two subcases. If s + 1 ∈ B, then by Lemma 7 we obtain that η s+1 ≥ αβ 2M1+ϕ(s+1) , which implies that (34) is satisfied. Otherwise, if s + 1 / ∈ B, then we haveEQUATION)where we used the induction hypothesis in the last inequality. Furthermore, note that ϕ(s)/ϕ(s + 1) ≤ s+2 s+1 ≤ 2 ≤ 1 β , which implies that ϕ(s) ≤ ϕ(s + 1)/β. Hence, we haveα 2M 1 + ϕ(s) ≥ αβ 2βM 1 + ϕ(s + 1) ≥ αβ 2M 1 + ϕ(s + 1).Therefore, (35) implies that η s+1 ≥ minαβ 2M1+ϕ(s+1) , σ0 β s+1and thus (34) also holds in this subcase. This completes the induction and we conclude that (34) holds for all t ≥ 0.As a corollary of Lemma 8, we obtain the following lower bound on η t for t ≥ T -x * ∥ 2 ≤ ∥x t -x * ∥ 2 (1 + 2αβµ/(3M 1 )) -1 = ∥x t -x * ∥ 2 (1 + 2αβ/(3κ)) -1 .In this section, we present the proof of Theorem 3. The proof consists of four steps. To begin with, in Lemma 9 we show that the iterate x t stays in a local neighborhood of x * when t ≥ I, where I is defined in (28). Next, we use this result in Lemma 10 to upper bound 1 √ µηt in those iterations where our line search scheme backtracks the step size, i.e., t ∈ B. Then we use induction in Lemma 11 to prove an upper bound for all t ≥ 0. Furthermore, we again use induction in Lemma 12 to establish an improved upper bound on 1 √ µηt when t ≥ T 2 , where T 2 is defined in (29). After proving Lemma 12, Theorem 3 then follows from Proposition 1. Lemma 9. We have ∥x t -x * ∥ ≤ νµ/L 2 for all t ≥ I, where I is given in (28).T1+u -x * ∥ 2 ≤ ∥x T1 -x * ∥ 2 1 + 2αβ 3κ -u. Moreover, since ∥x t -x * ∥ is non-increasing in t, we have ∥x T1 -x * ∥ ≤ ∥x 0 -x * ∥ = D. Thus, we have∥x T1+u -x * ∥ ≤ νµ L 2 ⇐ D 2 1 + 2αβ 3κ -u ≤ ν 2 µ 2 L 2 2 ⇐ u ≥ 2 1 + 3κ 2αβ log L 2 D νµ .This completes the proof.Note that by Proposition 1, we haveEQUATION)Hence, to characterize the convergence rate of our method, it is sufficient to upper bound the quantity 1/ √ 2η t µ. We achieve this goal in the subsequent lemmas. Lemma 10. If t ∈ B and t ≥ I, thenEQUATION)Moreover, it also holds thatEQUATION)Proof. By using the second bound in Lemma 1, we obtain thatEQUATION)Furthermore, by combining (33) and ( 9) in Lemma 2, we further haveEQUATION)It remains to bound the bias term ∥H t -Ht ∥. By Lemma 4, we have∥H t -Ht ∥ ≤ 1 t + 1 t i=0 ∥H t -H i ∥ = 1 t + 1 I-1 i=0 ∥H t -H i ∥ + 1 t + 1 t i=I ∥H t -H i ∥.For i = 0, 1 . . . , I -1, we use the first upper bound on ∥H t -H i ∥ in Lemma 4 to bound ∥H t -H i ∥ ≤ M 1 , and thusI-1 i=0 ∥H t -H i ∥ ≤ M 1 I.Moreover, for i = I, I + 1, . . . , t, we use the second upper bound in Lemma 4 to get ∥H t -H i ∥ ≤ 2L 2 ∥x i -x * ∥. Moreover, note that x i converges linearly to x * when i ≥ I by Lemma 6. Hence, we further haveEQUATIONIn the last inequality, we used the fact that ∥x I -x * ∥ ≤ νµ/L 2 and∞ i=0 (1 + ϕ) -i/2 = 1/(1 - (1 + ϕ) -1/2 ) = (1 + ϕ) 1/2 /((1 + ϕ) 1/2 -1) = (1 + ϕ) 1/2 ((1 + ϕ) 1/2 + 1)/ϕ ≤ 2(1 + 1/ϕ),where ϕ = 2αβ/(3κ). Moreover, since I ≥ 2 1 + 3κ 2αβ and ν ≤ 1, from (41) we further have1 t+1 t i=I ∥H t -H i ∥ ≤ 2µI t+1 .Combining the above inequalities, we arrive atEQUATIONCombining ( 39), (40), and (42) leads to the first result in (37). Finally, (38) follows from the fact that ∥ Ēt ∥ ≤ ϕ(t) and ∥x t -x * ∥ ≤ νµ/L 2 for all t ≥ I.Lemma 11. Assume that β ≤ 1/2. For any t ≥ I, we haveEQUATIONwhere ρ t is defined in (30).Proof. We shall use induction to prove Lemma 11. For t = I, note that by Corollary 2, we have η I ≥ αβ/(3M 1 ). Thus, this implies that1 √ 2η I µ ≤ √ 3κ √ 2αβ ≤ 3κI 2α √ β(I + 1) ,where we used the fact that κ ≥ 1, α < 1, β ≤ 1/2 and I ≥ 4. This proves the base case where t = I. Now assume that (43) holds for t = s, where s ≥ I. For t = s + 1, we distinguish two subcases. If s + 1 ∈ B, then by Lemma 10 we obtain that (43) is satisfied for t = s + 1. Otherwise, if s + 1 / ∈ B, then we have η s+1 = σ s+1 = η s /β. Hence, by using the induction hypothesis, we have1 √ 2η s+1 µ = √ β √ 2η s µ ≤ ν 4αβ (1 -α 2 )β + √ βϕ(s) 2α √ βµ + 3 √ βκI 2α √ β(s + 1) .Since β ≤ 1/2 and I ≥ 2, we have (s + 2)/(s + 1) ≤ (I + 2)/(I + 1) ≤ 1.4 ≤ 1/ √ β and ϕ(s) ≤ ϕ(s + 1)/ √ β. Thus, we further have1 √ 2η s+1 µ ≤ ν 4αβ (1 -α 2 )β + ϕ(s + 1) 2α √ βµ + 3κI 2α √ β(s + 2) .This shows that (43) also holds in this subcase. This completes the induction and we conclude that (43) holds for all t ≥ I.Before proving Lemma 12, recall the definition of ϕ in (31) and first defineEQUATIONNote that we have ϕ(t) ≤ νµ when t ≥ 256Υ 2 ν 2 log 8dΥ δν . Hence, by the definition of (29), it holds that T 2 ≥ T ′ 2 . Lemma 12. Recall the definition of ρ t in (30). For any t ≥ T ′ 2 , we haveEQUATIONProof. By Lemma 10, if t ∈ B, then1 √ µη t ≤ L 2 ∥x t -x * ∥ 2αβ 2(1 -α 2 )βµ + ρ t .We shall prove (45) by induction. First consider the base case where t = T ′ 2 , where T ′ 2 is defined in (44). To begin with, we will show that ν2α √ β ≤ ρ T ′ 2 ≤ 5ν 2α √ β . Since T ′ 2 is the maximum of I ′ and κI ν , we have either T ′ 2 = I ′ or T ′ 2 = κI ν -1.In the former case, we can lower boundρ T ′ 2 ≥ 1 2α √ βµ ϕ(I ′ ) ≥ ν 2α √ β .In the latter case, we can lower bound ρT ′ 2 ≥ 3κI 2α √ β(T ′ 2 +1) = 3ν 2α √ β .Combining both cases leads to the lower bound on ρ T ′ 2 . Furthermore, note that both two terms in ρ t are a decreasing function in terms of t, and hence we haveρ T ′ 2 ≤ 1 2α √ βµ ϕ(I ′ ) + 3κI 2α √ β(κI/ν) ≤ 2 2α √ βµ ϕ(I ′ + 1) + 3ν 2α √ β ≤ 5ν 2α √ β .This proves the upper bound on ρ T ′ 2 . Now we return to the proof in the base case where t = T ′ 2 . since ∥x T ′ 2 -x * ∥ ≤ νµ/L 2 by Lemma 9, we obtain thatL2∥x T ′ 2 -x * ∥ 2α √ βµ ≤ ν 2α √ β ≤ ρ T ′ 2 .Moreover, by Lemma 11, we have1 2η T ′ 2 µ ≤ ν 4αβ (1 -α 2 )β + ρ(T ′ 2 ) ≤ ν 4αβ (1 -α 2 )β + 5ν 2α √ β ≤ 1 2β √ 1 -α 2 + 5 ρ T ′ 2 .This shows that (45) holds for t = T ′ 2 . Now assume that (45) holds for t = s ≥ T ′ 2 . For t = s + 1, by using the induction hypothesis and (36), we obtain thatL 2 ∥x s+1 -x * ∥ 2α √ βµ ≤ L 2 ∥x s -x * ∥ 2α √ βµ √ 2η s µ ≤ 1 √ 2 1 2β √ 1 -α 2 + 5 ρ 2 s .Moreover, since ρ s /2 ≤ ρ s+1 , it suffices to show that 1√ 2 1 2β √ 1-α 2 + 5 ρ 2 s ≤ ρ s /2, which is equivalent to √ 2 1 2β √ 1-α 2 + 5 ρ s ≤ 1. Furthermore, since ρ s is non-increasing, we further have √ 2 1 2β √ 1 -α 2 + 5 ρ s ≤ √ 2 1 2β √ 1 -α 2 + 5 ρ T ′ 2 ≤ √ 2 1 2β √ 1 -α 2 + 5 5ν 2α √ β ≤ 1,where we used the condition on ν stated in Theorem 3 in the last inequality. This proves the first inequality in (45).To prove the second inequality in (45) for t = s + 1, we distinguish two subcases. If s + 1 ∈ B, then by Lemma 10, we have1 √ 2η s+1 µ ≤ L 2 ∥x s+1 -x * ∥ 4αβ (1 -α 2 )βµ + ∥ Ēs+1 ∥ 2α √ βµ + 3κI 2α √ β(s + 1) ≤ 1 2β √ 1 -α 2 ρ s+1 + ρ s+1 ≤ 1 2β √ 1 -α 2 + 5 ρ s+1 .Otherwise, if s + 1 / ∈ B, then we have η s+1 = η s /β and hence1 √ 2η s+1 µ = √ β √ 2η s µ ≤ 1 2β √ 1 -α 2 + 5 βρ s .Since T ′ 2 ≥ I ≥ 2 and β ≤ 1/2, we have ρ s /ρ s+1 ≤ (I + 2)/(I + 1) ≤ 1.4 ≤ 1/ √ β. Thus, we also proved that1 √ µηs+1 ≤ 1 2β √ 1-α 2 + 5 ρ s+1. This completes the induction.Proof of Theorem 3. It immediately follows from (36) and Lemma 12.In this section, we will present the formal version of Theorem 2. Our proof largely mirrors the developments in Section 4.Averaged stochastic error. Similar to Lemma 3, we can use tools from concentration inequalities to prove the following upper bound on the averaged stochastic error. Lemma 13 ([1, Lemma 6] ). Let δ ∈ (0, 1) with d/δ ≥ e. Then with probability 1 -δπ 2 /6, we have, for any t ≥ 0,EQUATIONWarm-up phase. Similar to the case of uniform averaging, we can only ensure that the distance to x * is monotonically non-increasing by Proposition 1 during this phase. Moreover, Algorithm 1 transitions to the linear phase when ∥ Ēt ∥ ≤ M 1 . Specifically, the transition point U 1 is given byEQUATIONWhen w(t) = (t + 1) log(t+4) , we have w ′ (t)/w(t) = O (log(t)/t). Thus, we conclude thatU 1 = Õ(Υ 2 /κ 2 ).Linear convergence phase. In the following lemma, we prove the linear convergence of Algorithm 1 with weighted averaging. Lemma 14. Assume that β ≤ 1/Ψ 2 and recall the definition of U 1 in (47). For any t ≥ U 1 , we have∥x t+1 -x * ∥ 2 ≤ ∥x t -x * ∥ 2 1 + 2αβ 3κ -1, where κ ≜ M 1 /µ is the condition number.Proof. See Appendix C.3.DefineEQUATIONMoreover, let ν ∈ (0, 1) be a parameter and defineEQUATIONFinally, letEQUATIONWhen w(t) = (t + 1) log(t+4) , we remark that J ′ = Õ(Υ 2 ) and thus J = Õ(κ + Υ 2 ). Moreover, similar to the derivation in [1] , it can be shown that U 2 = O(J ) = Õ(κ + Υ 2 ). Theorem 4. Let ν ∈ (0, 1) be a parameter satisfyingEQUATIONand recall the definition of U 2 in (50). For any t ≥ U 2 , we have∥x t+1 -x * ∥ ≤ 1 10β 2(1 -α 2 ) + 1 √ 2 θ t ∥x t -x * ∥,whereEQUATIONProof. See Appendix C.4.To simplify the notation, define the functionϕ(t) = 8ΨΥ E max log d(t + 1) δ w ′ (t) w(t) , log d(t + 1) δ w ′ (t) w(t) .Then we can rewrite (46) as ∥ Ēt ∥ ≤ ϕ(t) for all t ≥ 0. Similar to Lemma 7, we have the following result. Lemma 15. If t ∈ B, then we have η t ≥ αβ/(2M 1 + ∥ Ēt ∥) ≥ αβ/(2M 1 + ϕ(t)). Lemma 16. Assume that β ≤ 1/Ψ. For any t ≥ 0, we haveEQUATIONProof. We prove this lemma by induction. For t = 0, we distinguish two subcases. If 0 ∈ B, then by Lemma 7 we obtain that η 0 ≥ αβ 2M1+ϕ(0) . Otherwise, if 0 / ∈ B, we have η 0 = σ 0 . In both cases, we observe that (53) is satisfied for the base case t = 0. Now assume that (53) is satisfied for t = s. For t = s + 1, we again distinguish two subcases. If s + 1 ∈ B, then by Lemma 15 we obtain that η s+1 ≥ αβ 2M1+ϕ(s+1) , which implies that (53) is satisfied. Otherwise, if s + 1 / ∈ B, then we haveEQUATIONwhere we used the induction hypothesis in the last inequality. Furthermore, note thatϕ(s) ϕ(s + 1) ≤ w ′ (s)w(s + 1) w ′ (s + 1)w(s) ≤ Ψ ≤ 1 β ,and henceα M 1 + ϕ(s) ≥ αβ βM 1 + ϕ(s + 1) ≥ αβ M 1 + ϕ(s + 1).Therefore, (54) implies that η s+1 ≥ min αβ M1+ϕ(s+1) , σ0 β s+1and thus (53) also holds in this case. This completes the induction and we conclude that (53) holds for all t ≥ 0. Proof. This follows from Lemmas 13 and 14.Lemma 18. If t ∈ B and t ≥ J , thenEQUATIONand alsoEQUATION)Proof. Similar to the proof in Lemma 10, note that1 √ µη t ≤ L 2 ∥x t -x * ∥ 2αβ 2(1 -α 2 )βµ + ∥H t -Ht ∥ α √ 2βµ + ∥ Ēt ∥ α √ 2βµ .For the second term, note that Ht = t i=0 z i,t H i , where z i,t = w(i) -w(i -1) w(t) .Hence, by Jensen's inequality, we have∥H t -Ht ∥ ≤ t i=0 z i,t ∥H t -H i ∥ = J -1 i=0 z i,t ∥H t -H i ∥ + t i=J z i,t ∥H t -H i ∥.When 0 ≤ i ≤ J -1, we use Assumption 2 to bound ∥H t -H i ∥ ≤ M 1 and thusJ -1 i=0 z i,t ∥H t -H i ∥ ≤ M 1 J -1 i=0 w(i) -w(i -1) w(t) = M 1 w(J -1) w(t) .Moreover, for J ≤ i ≤ t, we use Assumption 3 to get∥H t -H i ∥ ≤ L 2 ∥x t -x i ∥ ≤ L 2 (∥x t -x * ∥ + ∥x i -x * ∥) ≤ 2L 2 ∥x i -x * ∥ ≤ 2νµ. Thus, t i=J z i,t ∥H t -H i ∥ ≤ 2νµ t i=I z i,t ≤ 2νµ.Combining the above inequalities, we arrive at∥H t -Ht ∥ ≤ M 1 w(J -1) w(t) + 2νµ.This leads to the first result in (55). To show (56), we note that ∥ Ēt ∥ ≤ ϕ(t) and ∥x t -x * ∥ ≤ νµ/L 2 for all t ≥ J .Lemma 19. Assume that β ≤ 1/Ψ 2 . For any t ≥ J , we haveEQUATIONProof. We shall use induction to prove Lemma 19. For t = J , note that by Corollary 3, we have η J ≥ αβ/(3M 1 ). Thus, this implies that1 √ µη J ≤ √ 2κ √ αβ ≤ √ 2κ α √ β ,where we used κ ≥ 1, α, β < 1 and I ≥ 2. This proves the base case where t = J . Now assume that (57) holds for t = s, where s ≥ J . For t = s + 1, we distinguish two cases. If s + 1 ∈ B, then by Lemma 18 we obtain that (57) is satisfied for t = s + 1. Otherwise, if s + 1 / ∈ B, then we have η s+1 = σ s+1 = η s /β. Hence, by using the induction hypothesis, we have1 √ µη s+1 = √ β √ µη s ≤ ν 2αβ 2(1 -α 2 )β + 3ν α √ 2β + √ 2βκw(J ) α √ βw(s) .Since β ≤ 1/Ψ 2 , we have w(s + 1)/w(s) ≤ Ψ ≤ 1/ √ β. Thus, we further have1 √ µη s+1 ≤ ν 2αβ 2(1 -α 2 )β + 3ν α √ 2β + √ 2κw(J ) α √ βw(s + 1).This shows that (57) also holds in this case.Recall that w(U 2 ) = w(J ) κ ν . Then by Lemma 19, for t ≥ U 2 we have1 √ µη t ≤ ν 2αβ 2(1 -α 2 )β + 3ν α √ 2β + √ 2κw(J ) α √ βw(U 2 ) ≤ ν 2αβ 2(1 -α 2 )β + 3ν α √ 2β + √ 2ν α √ β = ν 2αβ 2(1 -α 2 )β + 5ν α √ 2β .We will choose ν such thatEQUATIONTherefore, this further implies that ∥x t+1 -x * ∥ ≤ ∥x t -x * ∥/(2Ψ) for t ≥ U 2 .Lemma 20. If t ∈ B and t ≥ U 2 , thenEQUATIONProof. Recall that we have1 √ µη t ≤ L 2 ∥x t -x * ∥ 2αβ 2(1 -α 2 )βµ + ∥H t -Ht ∥ α √ 2βµ + ∥ Ēt ∥ α √ 2βµ .For the second term, we can write∥H t -Ht ∥ ≤ t i=0 z i,t ∥H t -H i ∥ = J -1 i=0 z i,t ∥H t -H i ∥+ U2 i=J z i,t ∥H t -H i ∥+ t i=U2+1 z i,t ∥H t -H i ∥.For the first part,J -1 i=0 z i,t ∥H t -H i ∥ ≤ M 1 w(J -1) w(t) ≤ M 1 w(J )w(t) . For the second part,U2 i=J z i,t ∥H t -H i ∥ ≤ 2νµ U2 i=J z i,t ≤ 2νµ w(U 2 ) w(t) = 2M 1 w(J ) w(t) ,where we used the fact that w(U 2 ) = w(J ) κ ν . For the third part,t i=U2+1 z i,t ∥H t -H i ∥ ≤ t i=U2+1 2L 2 w(i) -w(i -1) w(t) ∥x i -x * ∥ ≤ t-U2 j=1 2L 2 w(U 2 + j) -w(U 2 + j -1) w(t) ∥x U2 -x * ∥ (2Ψ) j ≤ t-U2 j=1 2L 2 w(U 2 + j) w(t) ∥x U2 -x * ∥ (2Ψ) j ≤ t-U2 j=1 2νµ w(U 2 ) w(t) w(U 2 + j) w(U 2 )(2Ψ) j ≤ 2νµ w(U 2 ) w(t) t j=1 w(U 2 + j) w(U 2 )(2Ψ) j ≤ 2νµ w(U 2 ) w(t) t j=1 1 2 j ≤ 2νµ w(U 2 ) w(t) = 2M 1 w(J ) w(t) .Therefore, we conclude that∥H t -Ht ∥ ≤ 3M 1 w(J ) w(t) + 2M 1 w(J ) w(t) = 5M 1 w(J ) w(t) .This completes the proof.Lemma 21. Recall the definition of θ t in (52). For any t ≥ 0, we haveEQUATIONProof. Note that by Proposition 1, we have∥x U2+t+1 -x * ∥ ≤ ∥x U2+t -x * ∥(1 + 2η U2+t µ) -1/2 ≤ ∥x U2+t -x * ∥ √ 2η U2+t µ . By Lemma 18, if U 2 + t ∈ B, then 1 √ µη U2+t ≤ L 2 ∥x U2+t -x * ∥ 2αβ 2(1 -α 2 )βµ + θ t .We will prove (60) by induction. First consider the base case t = 0. We note that θ 0 ≥ 5κw(J )α √ 2βw(U2) = 5ν α √ 2β . On the other hand, since ∥x U2 -x * ∥ ≤ νµ/L 2 , we obtain that 5L2∥x U 2 -x * ∥ α √ 2βµ ≤ 5ν α √2β ≤ θ 0 . Moreover, by Lemma 19, we have1 √ µη U2 ≤ ν 2αβ 2(1 -α 2 )β + 3ν α √ 2β + √ 2κw(J ) α √ βw(U 2 ) ≤ ν 2αβ 2(1 -α 2 )β + 5ν α √ 2β ≤ 1 10β √ 1 -α 2 + 1 θ 0 .This shows that (60) holds for t = 0. Now assume that (60) holds for t = s ≥ 0. For t = s + 1, by using the induction hypothesis, we obtain5L 2 ∥x U2+s+1 -x * ∥ α √ 2βµ ≤ 5L 2 ∥x U2+s -x * ∥ α √ 2βµ √ 2η U2+s µ ≤ 1 √ 2 1 10β √ 1 -α 2 + 1 θ 2 s . 1 10β √ 1-α 2 + 1 θ 2 s ≤ θ s /Ψ, which is equivalent to Ψ √ 2 1 10β √ 1-α 2 + 1 θ s ≤ 1.Furthermore, note that θ s is non-increasing andθ 0 ≤ ν α √ 2β + 5ν α √ 2β = 6ν α √2β . Thus, we only need to requireΨ √ 2 1 10β √ 1 -α 2 + 1 6ν α √ 2β ≤ 1 ⇔ 3 10αβ (1 -α 2 )β + 3 α √ β ν ≤ 1 Ψ ,which is satisfied due to (51). This proves the first inequality in (60).To prove the second inequality in (60) for t = s + 1, we distinguish two cases. If s + 1 ∈ B, then by Lemma 20, we have1 √ µη U2+s+1 ≤ L 2 ∥x U2+s+1 -x * ∥ 2αβ 2(1 -α 2 )βµ + ∥ ĒU2+s+1 ∥ α √ 2βµ + 5κw(J ) α √ 2βw(U 2 + s + 1) ≤ 1 10β √ 1 -α 2 θ s+1 + θ s+1 ≤ 1 10β √ 1 -α 2 + 1 θ s+1 .Otherwise, if s + 1 / ∈ B, then we have η U2+s+1 = η U2+s /β and hence1 √ µη U2+s+1 = √ β √ µη U2+s ≤ 1 10β √ 1 -α 2 + 1 βθ s . Since θ s /Ψ ≤ θ s+1 and √ β ≤ 1/Ψ, this implies that 1 √ µη T 2 +s+1 ≤ 1 10β √ 1-α 2 + 1 θ s+1 . This completes the induction.Proof of Theorem 4. By Proposition 1, we have∥x T2+t+1 -x * ∥ ≤ ∥x T2+t -x * ∥(1 + 2η T2+t µ) -1/2 ≤ ∥x T2+t -x * ∥ √ 2η T2+t µ .The rest follows from Lemma 21.Note that for t ≥ U 3 = Õ(Υ 2 + κ), we have ∥x t+1 -x * ∥ = Õ( Υ √ t )∥x t -x * ∥ by Theorem 2. By unrolling the inequality, this implies that∥x t+1 -x * ∥ ≤ ∥x U3 -x * ∥ Õ t s=U3 Υ √ s ≤ ∥x 0 -x * ∥ Õ t s=U3 Υ √ s .Further, for t ≥ 2U 3 , we can upper bound Υ √ s as follows: Using the Lambert W function 1 , the solution can be expressed as(i) Υ √ s ≤ Υ √ U3 ≤ 1 for any s ∈ [U 3 , t/2], (ii) Υ √ s ≤ Υ √ t/log t 2Υ 2 = W ( 2 Υ 2 log 1 ϵ ) ⇒ t = 2Υ 2 e W ( 2 Υ 2 log 1 ϵ ). Finally, by applying the bound e W (x) ≤ 2x+1 1+log(x+1) for any x ≥ 0, we conclude that t = O log(ϵ -1 ) log(Υ -2 log(ϵ -1 )) . Note that in the above derivation, we ignore the additional logarithmic factor log(t) in our superlinear convergence rate. However, a more careful analysis will show that it does not affect the final complexity bound. We also refer the reader to a similar derivation in [22, Appendix D.2] , where the authors provide the same complexity bound for a similar convergence rate of (1 + O( √ t)) -t .In Figure 3 , we test the effect of the extragradient step in our proposed SNPE method. We observe that in all cases, the variant without an extragradient step outperforms the original version, suggesting that the extragradient step may not be beneficial for minimization problems. Nevertheless, the SNPE method with the extragradient step, which is the one analyzed in our paper, still outperforms the stochastic Newton method in [1] .Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.Justification: There is no societal impact of the work performed.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).38th Conference on Neural Information Processing Systems (NeurIPS 2024).https://en.wikipedia.org/wiki/Lambert_W_function 


Decision: Overall, all reviewers were positive about this work and I support their decision.

Reviews: review: paper_summary: paper_summary: Newton method is well-known for its local quadratic convergence. However, the use of Hessian introduces additional computation challenges. One approach to tackling this issue is an inexact approximation of the Hessian. In this paper, the authors consider the finite sum minimization problem. They propose a stochastic approximation of the Newton method, where instead of the full Hessian they use its subsample approximation. The idea of the proposed method is based on another known approach, that uses weighted average over a subsample of Hessians. In this paper, the authors improve this approach using the Hybrid Proximal Extragradient (HPE) framework for strongly convex problems. As a result, the proposed algorithm achieves both linear and superlinear convergence areas in fewer iterations by improving the dependence on condition number of convergence rate.
strength_weakness: strength_weakness: Authors consider only a strongly convex setup. However they mention extending their approach to convex case, this seems to me as a major limitation of this work. Additionally, to check condition (3) authors employ line-search, which introduces an additional logarithmic factor in the convergence rate. However this can be the burden of the HPE framework, it seems not very significant, but still an issue. 

## Minor remarks
1. Line 168: remove "follows" from the end
2. Line 168: not Step 5, but Step 6
3. Sometimes you write "stepsize", sometimes - "step size". Please, be consistent.
questions: questions: 1. Why do you use $\sigma_{t+1} = \eta_t/\beta$? Please provide more details.
limitations: limitations: No limitations.
paper_summary: paper_summary: This paper uses the hybrid proximal extragradient framework to accelerate the convergence of Hessian average. The theoretical results significantly reduce the number of iterations to enter the linear phase, initial superlinear phase, and final superlinear phase when compared to the initial Hessian average method.
strength_weakness: strength_weakness: This paper does not provide empirical results for the proposed methods against AGD, which is necessary. 
The proposed methods still require exact gradient oracle and its iteration complexity is $O({\kappa}+\log(1/\epsilon))$, while the classical AGD method require only $O(\sqrt{\kappa}\log(1/\epsilon))$. It is very important to use numerical results to exhibit the benefits of using second-order information.
questions: questions: Refer to weakness part.
limitations: limitations: N/A
paper_summary: paper_summary: This work proposes a novel algorithm called the Stochastic Newton Proximal Extragradient method. Authors claim that their method reaches a superlinear convergence rate after $\mathcal{O}(\kappa)$ iterations, in contrast to the $\mathcal{O}(\kappa^2)$ iterations proved in previous work.
strength_weakness: strength_weakness: - The assumptions are very restrictive. For example, assumptions 3 and 5 are very strong, and few machine learning applications can satisfy such assumptions. 
- I am confused about how parts (c) and (d) of Theorem 1 prove superlinear convergence. From superlinear convergence, I expected to see $\rho^{2^t}$ where $\rho < 1$. How do the results show superlinear convergence?
- In Table 2, you mentioned the complexity of Damped Newton where the dependence on $\epsilon$ is $\log \log \epsilon$. However, for SNPE, the iteration complexity $\log \epsilon$ is similar to AGD. The iteration complexity of Damped Newton is called superlinear, while that of SNPE is called just linear. 
- In lines 342- 344, the dominating term in the complexity is the one that depends on $\epsilon$. The Damped Newton has better dependence than SNPE.
- Section 7 needs more numerical experiments. I expect to see a comparison between SNPE and Damped Newton. In that comparison, putting **real-time (and not iterations)** on the $x$-axis will be fair, as SNPE does a line search in every iteration. 
- In lines 334-337, the authors ignore the complexity of line-search in algorithm 1 (BLS). In the BLS, you compute $(1 + \eta \tilde{H})^{-1}g$ several times, which is very expensive, and authors ignore this while computing the complexity of SNPE.
questions: questions: Check weaknesses.
limitations: limitations: NA
paper_summary: paper_summary: The paper develops an accelerated scheme for strongly convex problems. As feedback it requires a deterministic first order oracles, but only an inexact Hessian estimator –– the main assumption being that the Hessian noise is mean zero and sub-exponential (implied by e.g. Hessian subsampling and $\Vert\nabla^2 f_i(x)\Vert$ bounded or Hessian sketching with a Gaussian sketch). The main idea is in combining the Hessian averaging scheme [1] with the hybrid proximal extragradient framework [18,19].
strength_weakness: strength_weakness: The main weaknesses seems to be that:

- It is unclear whether the scheme has relevance in practice, since the experiments seems to be run with a nonaccelerated version. For completeness I would suggest running the accelerated version, even if it performs worse (at least in the appendix). The theoretical results are strong and sufficient in themselves, so even though it would be better to find an instance that matches the theoretically guarantees, I don't see it as necessary.
- there are no guarantees for convex as also pointed out by the authors.

Neither are critical concerns – I mainly have a few remarks and questions regarding comparison with other methods and the experiments (see below).
questions: questions: Method comparison & theory:

- It would be informative to compare theoretically with existing (accelerated) results using exact Hessian. Does the proposed method recover the existing complexities?
- l. 369 can you really claim better condition number dependency when you choose $x_{k+1}=\hat{x}_k$, which is not covered by the theory? What convergence guarantees can be shown then, since the interpolation in Eq. 4 otherwise appears crucial for the acceleration?
- How does the structure of the scheme compare with [1] when $x_{k+1}=\hat{x}_k$? It seems that the main difference is the stepsize selection (error condition vs Armijo style backtracking line search). If this is true, the experimental comparison can seem a bit synthetic.
- If you have full gradients, is it not possible to use Pearlmutter's implicit Hessian-vector product to (inexactly) do the second order update while still maintaining a $\mathcal{O}(nd)$ like first-order methods? See e.g. https://www.cs.toronto.edu/~jmartens/docs/Deep_HessianFree.pdf
- Is it possible to relax the requirement of exact gradients? For first order methods it is still possible to achieve exponential convergence under e.g. relative noise (see e.g. Thm 5.1 in https://arxiv.org/pdf/2102.02921).

Experimentally:

- What are choices of hyperparameters $\alpha, \beta$ of the method? How are the baselines tuned in comparison?
- The proposed method can also be used with a deterministic Hessian. How does the stochastic version compare with an exact Hessian? This would show the influence of the first (slow) phase and provide an idealized baseline.

Minor:

- Is it possible to shave off the logarithmic factors hiding in the $\widetilde{\mathcal O}$-notation?
- Figure 1: is it a labeling mistake or is SNPE-UnifAvg consistently better than SNPE-WeightAvg?
- How come stochastic NPE suffers an addition $\Upsilon$ dependency in the superlinear rate (in comparison with stochastic Newton does not)?
limitations: limitations: See weaknesses.

Comments: title: Discussion
comment: I thank the authors for the answers. I don't have any further questions and increased my rating to 7.
Dear Reviewers,

Thank you for your effort in reviewing this paper. Please go over the authors rebuttal (and I suggest also other reviews and rebuttals to those) and see if it changes your mind or answers any of your concerns. In particular please engage in discussion with authors on points which have been clarified in the rebuttal or discuss why the rebuttal does not answer your concerns.

Thanks you!
AC
I thank the authors for their response and have no further questions.
Thank you for the details. Please add them in the updated version. 

I will raise my score to 5.
Thank you for reading our rebuttal and for the follow-up comments. 

**SNPE v.s. Newton in the regime $n = O(d)$.**

Thank you for the suggestion. We will include the additional plots for the regime where $n = O(d)$ in the appendix of our revision.

---

**The iteration complexity of $\log(\epsilon^{-1})/\log\log(\epsilon^{-1})$.**

Thank you for the question. Due to space constraints, we did not include the proof in the rebuttal. However, we will provide a proof sketch below.

Note that for $t \geq \mathcal{U}\_3 = \tilde{O}(\Upsilon^2 + \kappa)$, we have $\\|x_{t+1}-x^\*\\| = \tilde{O}(\frac{\Upsilon}{\sqrt{t}}) \\|x_{t}-x^\*\\|$ by Theorem 2. By unrolling the inequality, this implies that $$\\|x_{t+1}-x^\*\\| \leq \\|x_{\mathcal{U}\_3}-x^\*\\| \tilde{O}\left(\prod_{s = \mathcal{U}\_3}^t \frac{\Upsilon}{\sqrt{s}}\right) \leq \\|x_0 - x^\*\\| \tilde{O}\left(\prod_{s = \mathcal{U}\_3}^t \frac{\Upsilon}{\sqrt{s}}\right).$$ Further, for $t \geq 2\mathcal{U}\_3$, we can upper bound $\frac{\Upsilon}{\sqrt{s}}$ as follows: 
- $\frac{\Upsilon}{\sqrt{s}} \leq \frac{\Upsilon}{\sqrt{\mathcal{U}\_3}} \leq 1$ for any $s \in [\mathcal{U}\_3, t/2]$, 
- $\frac{\Upsilon}{\sqrt{s}} \leq \frac{\Upsilon}{\sqrt{t/2}}$ for any $s \in [t/2,t]$. 

Thus, we have $$\prod_{s = \mathcal{U}\_3}^t \frac{\Upsilon}{\sqrt{s}} \leq \prod_{s = t/2}^t \frac{\Upsilon}{\sqrt{s}} \leq \left(\frac{\Upsilon}{\sqrt{t/2}}\right)^{\frac{t}{2}}.$$ 

To derive a complexity bound, we upper bound the required number of iterations $t$ such that $ \left(\frac{\Upsilon}{\sqrt{t/2}}\right)^{\frac{t}{2}} = \epsilon$. Taking the logarithm of both sides and with some algebraic manipulation, we obtain $$\frac{t}{2\Upsilon^2} \log \frac{t}{2\Upsilon^2} = \frac{2}{\Upsilon^2} \log \frac{1}{\epsilon}.$$ Using the [Lambert W function](https://en.wikipedia.org/wiki/Lambert_W_function), the solution can be expressed as $ \log \frac{t}{2\Upsilon^2} = W(\frac{2}{\Upsilon^2} \log \frac{1}{\epsilon}) \Rightarrow t = 2\Upsilon^2 e^{ W(\frac{2}{\Upsilon^2} \log \frac{1}{\epsilon})}$. Finally, by applying the bound $e^{W(x)} \leq \frac{2x+1}{1 + \log (x+1)}$ for any $x \geq 0$, we conclude that $t = O\left(\frac{\log(\epsilon^{-1})}{\log(\log(\epsilon^{-1}))}\right)$. Note that in the above derivation, we ignore the additional logarithmic factor $\log (t)$ in our superlinear convergence rate. However, a more careful analysis will show that it does not affect the final complexity bound. We also refer the reviewer to a similar derivation in Appendix D.2 of Jiang et al. (2023), where the authors provide the same complexity bound for a similar convergence rate of $(1+ O(\sqrt{t}))^{-t}$. 

R. Jiang, Q. Jin, and A. Mokhtari. "Online learning guided curvature approximation: A quasi-Newton method with global non-asymptotic superlinear convergence." COLT 2023.
I thank the authors for their response and have no further questions.
Thank you for your explanation of superlinear convergence and the plots. 

So SNPE performs better than Newton in the regime $n >> d$. Adding plots for other regimes to the appendix will be beneficial (just a suggestion to check the performance of SNPE compared to Newton).

I will raise my score to 4. 

I am still not convinced about the iteration complexity in the table. You mentioned the iteration complexity of SNPE has dependence $\frac{\log(\varepsilon^{-1})}{\log \log(\varepsilon^{-1})}$. Is there a proof of this in the appendix that I can check? Else, can you add the proof here?