Stochastic Variance-Reduced Newton: Accelerating Finite-Sum Minimization with Large Batches
Authors: Michal Derezinski
Venue: Accepted by TMLR ()

Abstract: Stochastic variance reduction has proven effective at accelerating first-order algorithms for solving convex finite-sum optimization tasks such as empirical risk minimization. Incorporating second-order information has proven helpful in further improving the performance of these first-order methods. Yet, comparatively little is known about the benefits of using variance reduction to accelerate popular stochastic second-order methods such as Subsampled Newton. To address this, we propose Stochastic Variance-Reduced Newton (SVRN), a finite-sum minimization algorithm that provably accelerates existing stochastic Newton methods from $O(\alpha\log(1/\epsilon))$ to $O\big(\frac{\log(1/\epsilon)}{\log(n)}\big)$ passes over the data, i.e., by a factor of $O(\alpha\log(n))$, where $n$ is the number of sum components and $\alpha$ is the approximation factor in the Hessian estimate. Surprisingly, this acceleration gets more significant the larger the data size $n$, which is a unique property of SVRN. Our algorithm retains the key advantages of Newton-type methods, such as easily parallelizable large-batch operations and a simple unit step size. We use SVRN to accelerate Subsampled Newton and Iterative Hessian Sketch algorithms, and show that it compares favorably to popular first-order methods with variance~reduction.

Published in Transactions on Machine Learning Research (01/2025)
Stochastic Variance-Reduced Newton: Accelerating Finite-Sum Minimization with Large Batches
Michat Dereziński
dereain@umich..omm edu
Deportment of Electrical Engineering and Computer Science
University of Michigan
Reviewed on OpenReview: https: openreview. net/ forun? id= daQCRHKRdC
Abstract
Stochastic variance reduction has proven effective at accelerating first-order algorithms for solving convex finite-sum optimization tasks such as empirical risk minimization. Incorpo- rating second-order information has proven helpful in further improving the performance of these first-order methods. Yet. comparatively little is known about the benefits of using vari- ance reduction to accelerate popular stochastic second-order methods such as Subsampled Newton. To address this, we propose Stochastic Variance-Reduued Newton (SVRN), a finite- sum minimization algorithm that provably accelerates existing stochastic Newton methods from (a log(L/))) to O n is the number of sum components and a is the approximation factor in the Hessian esti- mate. Surprisingly, this acceleration gets more significant the larger the data size n, which unique property of SVRN. Our algorithm retains the key advantages of Newton-type methods, such easily parallelizable large-batch operations and a simple unit step size. We SVRN to accolerate Subsampled Newton and Iterative Hessian Sketch algorithms, and show that it compares favorably to popular first-order methods with variance reduction.
log(1/c) passes over the data, i.e,, by a factor of ((a log(n)), where
log(n)
1 Introduction
Consider a convex finilesum minimization task:
find x" = argmin f(x) for f(x) = Σ (x))..
(1)
XERd
This optimization lask naturally arises machine learning through empirical risk minimization, where x is the model parameter vector and each function (x)) corresponds to the loss incurred by the model on the i-th element in training data set (e.g., square loss for regression, logistic loss for classification). Many other optimization tasks, such as solving semi-definite programs and portfolio optimization, this general form. Our goal is to find an 6-approximate solution, i.e., X such that f(x) f(x*) <.
be cast in
Classical iterative optimization methods (such as gradient descent and Newton's 's method) use first rres//scood- order information function f to construct a sequence Xg. XI.... that converges to x" However, this does not leverage the finite-sum structure of the problem. Thus, extensive literature has been dedicated to efficiently solving finite-sum minimization tasks using stochastic optimization methods, which use firso information of randomly sampled component functions v/,, that can often be computed much faster than the entire function f. Among first-order methods, variance-reduction techniques such as SAG (Roux et al., 2012), SDCA (Shalev-Shwartz & Zhang, 2013), SVRG (Johnson & Zhang, 2013), SAGA (Defazio et al., 2014), Katyusha (Allen-Zhu, 2017) and others (Frostig et al., 2015; Konecný el al., 2015; Allen-Zhu & Yuan, 2016), have proven particularly effective. One of the most popular variants of this approach is Stochastic Rrrancc-RRedued Gradient (SVRG), which achieves variance reduction by combining frequent stochastic gradient queries with occasional full batch gradient queries, to optimize the overall cost of finding an e- approximate solution, where the cost measured by the total number of queries the components (xl.(x).

---

Published in Transactions on Machine Learning Research (01/2025)
Many stochastic second-order methods have also been proposed for solving finite-sum minimization, including Subsampled Newton (Frdogdu & Montanari, 2015; Roosta-Khorason Khorasani & Mahoney, 2019; Bollapragada et al., 2018; Berahas et al., 2020; Dereziński & Mahoney, 2019) Newton Sketch (Pilanci & Wainwright, 2016; 2017; Dereziński et al., 2021), and others (Kovalev et al., 2019; Moritz et al., 2016; Tripuraneni et al., 2018; Mokhtari et al.; 2018; Gupta al., 2021). These approaches are generally less sensitive to hyperparameters such as the step size, and they typically query larger random batches of component gradients/ Boddiens/HHesiaas Hessiias at a time, as compared first-order methods. The larger queries make these methods less sequential, allowing for more effective vectorization and parallelization.
A number of works have explored whether second-order information can be used to accelerate stochastic variance-reduced methods, resulting in several algorithms such as Preconditioned SVRG (Gonen et al., 2016), SVRG2 (Gower et al., 2018) and others (Gower et al., 2016; Liu et al., 2019). However, these are still primarily stochastic first-order methods, highly sequential and with a problem-deeendeent step size. Comparatively little work has been done on using variance roduction accelerate stochastic Newton-type methods for convex finite-sum minimization (see discussion in Section 2.3). To that end, we ask:
Can variance reduction accelerate local convergence of Stochastic Newton in conver finite-sum minimazaiion?
We show that the answer this question is positive. The method that we use demonstrate this, which we call Stochastic Variance-Reduced Newton (SVRN), retains the positive characteristics of second-order methods, including easily parallelizable large-batch gradient queries, well as minimal hyperparameter tuning (e.g., accepting unit step size). We prove that, when the number of components 1/; is sufficiently large, SVRN achieves a better parallel complexity than SVRG, and a better sequential complexity than the corresponding Stochastic Newton method (see Table 1).
2
Main result
In this section, we present our main result, which is the parallel and sequential complexity analysis of the local convergence of SVRN. The algorithm itself is discussed detail in Section 3.
We now present the assumptions needed for our main result, starting with p-strong convexity of and X-smoothness of each the These are standard for establishing linear convergence rate of SVRG. Our result also requires Hessian regularity assumptions (Definition 1), which are standard for Newton's method and only affect the size of the local convergence neighborhood.
We assume that f(x) = n E. ve(z) has continuous first and second derivotives, as well as a bounded condition number A/P, where a and X are defined as follows:
Assumption 1
1. Punction f is u-strongly convex, i.e.,
f(x) 2 f(x) + f(('') - x') - 1|x - x'112
2. Each of the components the is X-smooth, i.e.,
(x) S (x() + V (x-x) + || -
To highlight the parallelizability of SVRN due to large mini- hini-batches, as well as the effect of variance reduction on its performance, we will consider two standard complexity measures:
1. Parallel complexity: Number of batch gradient queries, i.e., times the algorithm computes a gradient at an iterate, over the full batch or a mini-batch. This corresponds to the PRAM model.
2. Sequential complexity: Number of queries to component gradients (x)()), normalized by the total number of components n. Thus, one can think of this as the number of "data passes", where one "data pass" may possibly query component gradients at different itcrates. This is a natural measure of complexity for stochastic first-order methods, which coincides with the parallel complexity when only query full-batch gradients.
2

---

Published in Transactions on Machine Learning Research (01/2025)
The motivation for this dual-complexity perspective is centered on the benefits of using large-batch gradient queries. These benefits come from the fact that on most computing architectures it takes far less time to compute a single gradient estimate on a batch of component functions at one location, than it takes to compute component gradients at different locations, one at a time. We highlight this by the notion of parallel complexity, which measures the number of batch gradient queries. Note that the benefits of large mini-batch sizes are not limited to parallel computing- Even standard single- and multi-core architectures benefit substantially from the vectorization of gradient computations, which is only effective when using large batches (as we observe in our experiments, see Section 5).
Despite these benefits, parallelization and vectorization still come with some overhead, which is why it is natural to also optimize over the sequential complexity of the algorithms. Altogether, to express this in our model, we optimize both over parallel and sequential complexity, while prioritizing the parallel one. Note that for any algorithm using only full gradients (such as the standard ver- sions of gradient descent or Subsampled Newton), the two notions of complexity are exactly equivalent. For example, in gradient descent (GD), both parallel and sequential complexity is O(s (og(/e)))) By introducing stochastic gradients and variance reduction, as in SVRG, we can improve upon the sequential complexity of GD, while preserving (but not improving) its parallel complexity. Specifically, when > K, then SVRG with optimal muiii-batch size takes (loo((//e) sequential time to find an e-approximate solution, however it still needs ( ((xloo)(/))) parallel time (Table 1).
We can avoid the dependence of parallel complexity on the condition number k by using second-order information. In particular, suppose that in each iteration we compute the full gradient V f(x) and are given access to (typically stochastic) Hessian estimate H such that for some 1 S a & K.
(Hessian a-approximation)
v2s(x) x H (VV
In fact, our condition can be stated in even more general way, by asking that & VV HH for some lixed 0 < 31 < 82 such that = 82/31. In this case, (1/vB132)H is an a-approximation in the sense of (2), which means that we can apply our results and algorithms given any such Hessian estimates.
A standard Stochastic Newton (SN) update, given below in (3), can achieve parallel and sequential complexity of (a log(1/e)) locally in the neighborhood of the optimum. It is thus natural to ask whether we can use stochastic gradients and variance reduction to accelerate the local sequential complexity of this method, while preserving its parallel complexity. Our main result shows not only that this is possible, but also, remarkably, that this acceleration gets more significant the larger the data size n. See also Theorem 3 for algorithmic details and convergence analysis.
Theorem (informal Theorem 3) Suppose that Assumplion 1 holds and: (a) I has a Lipschilz Hessian, or (b) f is self-concordant. Moreover, let n > k > a. There is an algorithm (SVRN) and an open neighborhood U such that, given any x E U with a corresponding Hessian a a-approzimation as in (2), the cost of returning & such that f - f(x*) ) S e (f(x) - f(x')) as follows:
Sequential time = o log(1/c) data passes. Tog(n)
Parallel time = log(//e) batch queries
and
Remark 1 SVRN improves on the sequential complexity of Stochastic Newton by O(a log(n)) while retain- ing the same parallel complezity. Moreover, if S 2, then the algorithm accepts unit step size, and still achieves ((lgg()) acceleration. Note that this acceleration improves with the problem size n, which is unique property of SVRN.
Remark 2 To find an initialization point for SVRN, one can simply run a few iterations of a Subsampled Newton method with line search. In Section 4, we propose a globally convergent algorithm based on this approach -SVRN-HA; see Algorithm 1), and in Section 5 we show empirically that it substantially accelerates Subsampled Newton.
3

---

Published in Transactions on Machine Learning Research (01/2025)
Second-order
Parallel (batch queries)
Sequential (data passes)
Gradient Descent
x
O (O))og(1/e)
( (oo(((1/)))
Accelerated GD
X
0 x log(1/e)
0 log(1/e))
Stochastic Newton
( log(1/e))
( log(1/e))
Mini-batch SVRG
X
( (1)()
(log(1//e)
Mini-batch Katyusha
0 (V log(1/e))
(log(1//e)
SVRN (this work)
(og/og()))
0
log(1/c) (og(
Table 1: Comparison of local convergence behavior for SVRN and related stochastic methods in the big data regime, i.e., > K, along with full-batch Gradient Descent (GD), and Accelerated GD. Time complexities are obtained by first optimizing parallel time (batch queries), and then optimizing sequential time (data passes). For the second-order methods, we assume a Hessian a-approximation (2) where & K.
2.1 Discussion
In this section, we compare the local convergence complexity of SVRN standard stochastic first-order and socond-order algorithms. In this comparison, we focus what call the big data regime, i.c., K, which is of primary interest in the literature on Subsampled Newton methods. Then, in Section 2.2, illustrate how SVRN can be further improved via sketching and importance sampling, when solving problems with additional structure, such as least squares.
Comparison to SVRG and Katyusha. As we can see in Table 1, first-order algorithms, including variance-reducce methods such as SVRG (Johnson & Zhang, 2013), and its accelerated variants like Katyusha (Allen-Zhu, 2017), suffer from a dependence on the condition number & in their parallel complexity, Namely, they require either log(//)) or log(/)) batch gradient queries, compared to O(a log(1/c)) for SVRN and Stochastic Newton, where is the Hessian approximation factor, which is often much smaller than Vk. This is because, unlike SVRN, these methods do not scale well to large mini-batches, making them less parallelizable
Another difference between SVRN and SVRG or Katyusha is that, when the Hessian approximation is sufficiently accurate (a S 2), then SVRN accepts unit step size, which leads to optimal convergence rate without any tuning. On the other hand, the optimal step size for SVRG depends on the strong convexity and smoothness constants and X, and thus, requires tuning.
Comparison to Stochastic Newton. We next compare SVRN with Stochastic Newton methods such as Subsampled Newton and Newton Sketch. Here, the most standard proto-algorithn considered in the literature is the following update:
8041= = -, v)(.)).
(3)
As mentioned earlier, this update uses only full gradients, so both its parallel and sequential complexity is (aa log(1/c)) (see Stochastic Newton in Table 1). On the other hand, SVRN provides a direct accoleration of the sequential complexity without sacrificing any parallel complexity.
The Hessian a-approximation H V2 f(t.) can be produced in a number of ways, but perhaps the most relevant for this discussion is Hessian subsampling, a.k.a. Subsampled Newton (c.g., & Mahoney, 2019). In this setting, given & component Hessians sampled uniformly at random, can with high probability construct an estimate H with approximation factor a = 1+ 1100 (K (olog())/ + Vnlog(d)/K (see Appendix D.2), so that a & k for any k > log(d), and a < 2 for k = Ok log(d)). However, if we wanted to recover SVRN's sequential complexity of O log(1//))
Roosta-Khorasse
log(n) purely by improving the Hessian approximation in Subsampled Newton, the required Hessian sample size i would become least as large as n, meaning that would ealentially have to use the exact Hessian (i.e., Newton's method), which is highly undesirable.
4

---

Published in Transactions on Machine Learning Research (01/2025)
We note that some of the literature on Subsampled Newton proposes to subsample both the Hessian and the gradient (e.g., Bollapragada et al., 2018), which would be akin to Xx+1 = X, 7/s 1 Ei= Vilr, (X.). However, as is noted in the literature, to maintain linear convergence of such a method, one has to keep increasing the gradient sample size at an exponential rate, which means that, for finite-sum minimization, we quickly revert back to the full gradient. (see experiments in Section 5.2).
2.2 Accelerating SVRN with sketching and importance sampling
When the minimization task possesses additional structure, then we can combine SVRN with Hessian and gradient estimation techniques other than uniform subsampling. For example, one such family of tech- niques, called randomized sketching (Drineas & Mahoney, 2016; Woodruff, 2014; Murray et al., 2023; Dereziński & Mahoney, 2024), is applicable when the Hessian can be represented by a decomposition fx) = (x)) Ax(x) + C, where Ax(x) is a tall d matrix and C a fixed d x d matrix. This setting applies for many empirical risk minimization lasks, including linear and logistic regression, among others.
Sketching can be used to construct an estimate of the Hessian by applying a randomized linear transformation to A (x), represented by a k x n random matrix S, where k = ( is much smaller than n. Using standard sketching techniques, such as Subsampled Randomized Hadamard Transforms (SRHT, Ailon & Chazelle, 2009), Sparse Johnsoon-Lindenstrauus Transforms (SJLT, Clarkson & Woodruff, 2017; Nelson & Nguyên, 2013; Chenakkod et al., 2024) and Leverage Score Sparsified embeddings (LESS, Dereziński et al., 2021), we can construct Hessian estimate that satisfies the requirements of Theorem 1 at the cost of Ö(nd + dº), which corresponds to nearly-constant number of data passes and d x d matrix multiplies. In particular, this eliminates the dependence of Hessian estimation on the condition number.
Another way of making SVRN more efficient is to use importance sampling in the stochastic gradient mini- batches. Importance sampling can be introduced to any finitc-sum minimization task (1) by specifying an r-dimensional probability vector p = (pI,...,P.)) such that E. P4 = 1, and sampling the component gradient 1 (x) so that the index 1, is drawn according to p. With the right choice of importance sampling, we can substantially roduce the smoothness parameter X, and thereby, the condition number и of the problem (sce Appendix A.4).
The above techniques can be used to accelerate SVRN, for instance, in the important task of solving least squares regression. Here, given an n x d matrix with rows a and an n-dimensional vector y, the objective being minimized is the following quadratic function:
f(x)
1 _ 2n
||Ax Ax
y/1²
1 E (ajx 22 n
(4)
One of the popular methods for solving the least squares task, known as the Iterative Hessian Sketch (IHS, Pilanci & Wainwright, 2016), is exactly the Stochastic Newton update (3), where the Hessian estimate H is constructed via sketching, In this context, SVRN can be viewed as an accelerated version of IHS. To fully leverage the structure of the least squares problem, we use popular importance sampling technique called leverage score sampling (Drineas et al., 2006; 2012), where the importance probabilities are (approximately) proportional to Pa & a,' (A Aa, Through an adaptation of our main result, we show that a version of SVRN for least squares, with sketched Hessian and leverage score sampled gradients, improves on the complexity for reaching O(ndlo(11/)) (Rokhlin & Tygert, 2008; Avron et al., 2010; Meng et al., 2014) to O(dd log(1/e) ) Sec ing(n/d))-
solution to preconditioned least squares task from
Appendix A.4 for proof and further discussion.
Theorem 2 (Fast least squares solver) Given A E Raxd and y E R after O(d n + d log d) pre- processing cost to find the sketched Hessian estimate H and an approzimate leverage score distribution, SVRN finds & so that
S(x) < (1 (1+))((x))) in o ndd log(1/c)
lime.
log(n/a]
5

---

Published in Transactions on Machine Learning Research (01/2025)
Crucially, the SVRN-based least squares solver only requires preconditioner H that is constant factor approximation of the Hessian, i.e., 0(1). Interestingly, our approach of transforming the problem via leverage score sampling appears to be connected to a weighted and preconditioned SGD algorithm of (Yang et al., 2017) for solving a more general class of ép-rogression problems. We expect that Theorem 2 can be similarly exlended beyond least squares regression.
2.3 Further related work
As mentioned in Section 1, a number of works have aimed to accelerate first-order variance reduction meth- ods by preconditioning them with second-order information. For example, (Gonen et al., 2016) proposed Preconditioned SVRG for ridge regression. The effect of this proconditioning, as in related works (Liu et al., 2019); reduced condition number of the problem. This is different from Theorem 1, which uses preconditioning to make variance reduction effective for large mini-batee and with unit step size.
Some works have shown that, instead of preconditioning, one can use momentum to accelerate variance reduction, and also to improve its convergence rate when using mini-batches. These methods include Catalyst (Lin ct al., 2015) and Katyusha (Allen-Zhu, 2017). However, unlike SVRN, these approaches are still limited to fairly small mini-batches, as demonstrated in Table
A number of works have proposed applying techniques inspired by variance reduction to stochastic Newton- type methods in scttings which are largely incomparable to ours. First, (Rodomanov & Kropotov, 2016; Ko- valev et al., 2019) consider algorithms where the Hessian and gradient information is incrementally updated with either individual samples or mini-batches. However, the approximate Hessian information required by these methods is quite different than the one used in SVRN: they require the Hessian estimate to be initialized with all n component Hessians, possibly computed at different locations (compared to, e.g., a subsampled estimate). For example, in the case of least squares, this means computing the exact Hessian, which costs O(nd)) time and renders the task trivial (compare this to our Theorem 2, where the Hessian estimate required by SVRN can be approximated efficiently). Sctting this aside, we still compare the local convergence rate of SVRN with the Stochastic Newton method of (Kovalev et. al., 2019, Theorem 1) using the same mini-batch size. Assuming n > к and using the setup from Theorem 1, their method achieves (llo((//e) data passes, whereas SVRN obtains the accelerated complexity of O (ogi(/e)) data passes.
In the non-convex sctting, variance reduction was used by (Zhou et al., 2019; Zhang et al., 2022) to accelerate Subsampled Newton with cubic regularization. They use variance reduction both for the gradient and the Hessian estimates. Also, (Wang et al., 2017) incorporate variance reduction into a stochastic quasi-Newton method. However, due to the non-convex setting, these results are incomparable to ours, as are focusing on strongly convex optimization.
3 Local convergence analysis of SVRN
In this section, we present the convergence analysis for SVRN, leading to the proof of Theorem
1.
Notation. For x d positive semidefinite matrices A and B, we define ||v|a = v Av, and we say that xe B, when (1 B)B x A 3 (1 + e)B, where 3 denotes the positive semidefinile ordering (we define analogous notation a R b for non-negative scalars a, b). We use and C to denote positive absolute constants, and let ~~ [n] denote a uniformly random sample from nI.
We will present the analysis in slightly more general setting of expected risk minimization, i.e., where fxx = Eg [x)))( Here, D is a distribution over convex functions 1 Rdd - RR Clearly, this setting subsumees (1), since we can let D be a uniformly random sample v.. Thanks to this extension, our results can apply to importance sampling of component functions, as in Theorem 2.
Definition 1
We define the local convergence neighborhood Uy(e), parometerized by e E (0, 1), as: 1. Tf f has an L-Lipschitz Hessian, then Uy(c) = {x ||x - x" [[ <cp //////111
2. T f is self-concordant, then use Uy(t)= (x |x==x"||==(xx < </4).

---

Published in Transactions on Machine Learning Research (01/2025)
Complexity of SVRN
Parallel complexity
Sequential complexity
(lalog()/))
butl/s] Togin
Mini-batch size m
Figure 1: Illustration of the local convergence complexity analysis for SVRN, function of the mini-batch size m, with the number of inner iterations set to fmax n/m. As we decrease the mini-batch size from n (standard Stochastic Newton; SN) downto m log(n//)) (optimal SVRN), the sequential complexity (number of passes over the data) improves by O(a log(n)), while the parallel complexity (number of batch gradient queries) remains optimal.
Our local convergence analysis is captured by the following theorem, which provides the rate of convergence after one outer iteration of SVRN (stated below), for range of mini-batch sizes m.
Theorem 3 (Convergence rate of SVRN) Suppose that Assumption I holds, a 2 1, and either: (a) f has a Lipschitz continuous Hessian, or (b) f is self-comcordant. There is an absolute constant > 0 such that if &g E (1/c/c)), and we are given the gradient g. = ff(e) as well as a Hession a-approzimation, Va 2 ((x,) HH vo V2 ((x)), then, letting Xy X, and:
i.e., H such that
m
Xx+1 = X
nH
(xi( (x,) +8
in Um D,
after iterations with mini-batch size m 2 2 log(t/6) and step size = min min(v2/a, 1}, the iterate X,+1 = Xx (t.e,, one outer iteration of SVRN) with probability satisfies:
f(x++1) -fo-)) f(xx) f(xx) f(8,) f(x*) <
(1-22))
2a = + ca log(t/5)
The dependence on the condition number k in our result comes from sub-sampled gradient es- timation. This is consistent with Subsampled Newton works such as Roosta-Khorasani & Mahoney (2019): to recover their fast condition number-free convergence guarantees they require the gradient somple size to be sufficiently larger than K. We showed how this can be avoided for certain losses (least squares; see Theorem and Lemma 6) by relying on importance sompling.
Remark 3
The proof of Theorem 3, which given in Appendix A, relies on a new high-probability bound for the error of the variance-reduced gradient estimates in the large mini-batch regime, measured using the vector norm defined by the inverse Hessian at the optimum (Lemma 3). Unlike results from prior work, which hold in expectation, this bound crucially relies the iterate being in the local neighborhood. Also, unlike standard SVRG analysis, we achieve our convergence guarantee for the last ilerate of SVRN's inner loop (as opposed a random or averaged iterate), which is again enabled by exploiting local second-order information.
Discussion. For simplicity, let us fix the number of inner itcrations tmax = n n/m, m, so that single outer iteration of SVRN always takes two passes over the data. Then, we can define the linear convergence rate after one outer iteration as function of mini-batch size m:
Pm. (1-
1 n/m
Ō(k/m).
2x
7

---

Published in Transactions on Machine Learning Research (01/2025)
Input: iterate Xo, gradient batch size m, Hessian sample size k; and local iterations maaxi Initialize step size 1/-1 = 0 and Hessian estimate H-1 0: for 8 = 1, do
Compute the subsampled Hessian: H, = - (&.)( for 11...., v/ D; Compute the Hessian average: H, " -,,-1 + s+1 H
s+1 Compute the full gradient: g. fxx,): if 7/s-1 < then
Compute the descent direction v. by solving: H.V. = -
else
Initialize X = Xsi
for == O,..., mmax - do
Compute gr (x,) and B. (x.); for =xxx= V $1...., vm D:
m
Compute aariance-reduced gradient g = (x,) (x,) + gai
Compute the descent direction V, by solving: H.V; = -
Update X = Xx + Vt
end
Compute the descent direction: V. Xtmax X,i
end
Compute 7s for iterate X, and direction v, using the Armijo condition; Update 8041 = X, Vsi
end
Algorithm 1: SVRN with Hessian Averaging (SVRN-HA)
Let us assume the big data regime, i.e. n >> K. If we only use full- batl-batch gradients (m= = n), then the first term in the rate dominates, and we have Pm Stochastic Newton (3). As we decrease m (and change mmxx accordingly), the first term in Pm decreases, whereas the second term increases. As result, the overall rate rapidly improves, reaching its optimal value of Pm Ó(s/n) for m s login/k3
2n which is similar to what we would get using standard
Complexity analysis. The complexity analysis given in Theorem 1 follows directly from the above dis- cussion, since the sequential complexity (number of data passes needed to improve by factor e) is given by
(( Tog(
log(l/e) log(1/e)_ whereas the parallel complexity (number of batch gradient queries) is Of Tog(1/pm)) Figure 1, we illustrate how these quantities change as function of m. In particular, we observe that the batch gradients essentially stay flat at (( log(1/c)) as we decrease m, until reaching Ton(na/a) On the other
hand, the data pass complexity decreases lincarly with m, until it reaches the optimal value of O (og(n/x) ). which, for sufficiently large n, recovers Theorem 1.
4 Globally convergent algorithm
We next present practical stochastic second-order method (see Algorithm 1, called SVRN-HA) which uses SVRN accelerate its local convergence phase.
The key in implementing SVRN is that the algorithm is guaranteed to converge with unit step size only once reach local neighborhood of the optimum, and if we have sufficiently accurate Hessian estimate. For this reason, we introduce an initial phase of the algorithm, in which standard Stochastic Newton method is ran, using the Armijo line search to select the step size. Once the method reaches the local convergence neighborhood, as long as the Hessian estimates are accurate enough, the line search is guaranteed to return a unit step size. At this point, the algorithm switches to SVRN and achieves acceleration. Finally, to ensure that we reach a sufficiently accurate Hessian estimate, our Stochastic Newton method should gradually increase the accuracy of the Hessian estimates.
8

---

Published in Transactions on Machine Learning Research (01/2025)
Based on these insights, propose an algorithm called Stochastic Variance-Reduccee Newton with Hessian Averaging (SVRN-HA). In the initial phase, this algorithm is variant of Subsampled Newton, based on a method proposed by (Na et al., 2022), where, at each iteration, we construct subsampled Hessian estimate based on a fixed sample size k. To increase the accuracy over time, all past Hessian estimates are averaged logether, and the result. is used to precondition the full gradient. At each iteration, we check whether the last line search returned a unit step size. If yes, then we start running SVRN with local iterations Imax [ogg((/d)] and gradient batch size m = Ln/ logg(n/d) (n/d),, where n is the number of data points and d is the dimension. This is motivated by our theory (see discussion below Theorem 3), using d as a proxy for the condition number τ. This choice has proven effective in all of the evaluated datasets, which indicates that the theoretical dependence of the mini-batch size on the condition number is likely quite pessimistic.
In the following result, we establish global convergence of SVRN-HA, by showing that the global phase of this method will not only reach any local neighborhood, but also that the Hessian estimate will get progressively more accurate, eventually reaching the desired approximation accuracy.
Theorem 4 Let f be as in Theorem 3. For any neighborhood U around the optimum, Algorithu 1 will almost surely reach a point where: (a) X, belongs to the neighborhood U, and (b) the Hessian estimate H, satisfies the condition in Theorem 3, At this point, the line search will return Be = 1.
5 Experiments
We next demonstrate numerically that SVRN can be effectively used to accelerate stochastic Newton methods in practice. We also show how variance reduction can be incorporated into a globally convergent. Subsampled Newton mothod in a way that is robust to hyperparameters and preserves its scalability thanks to large-batch operations." 1
5.1 Logistic regression experiment
In this section, we present numerical experiments for solving regularized logistic loss minimization task. For x d data matrix A with rows a,", n-dimensional target vector y (with #1 entries y4) and a regularization parameter 7, our task is to minimize:
f(x)
log(1 +e aannaaa x) 2 ||x| 2
(5)
As a dataset, we used the Extended MNIST dataset of handwritten digits (EMNIST, Cohen et al., 2017) with n = 500k datapoints, transformed using a random features map (with dimension = 1000). Experimental details, as well as further results on the CIFAR-10 dataset and several synthetic data matrices, are presented in Appendix B.
In Figure 2, we compared SVRN-HA to three baselines which are most directly comparable: (1) the classical Newton's method; (2) SVRG with the step size and number of inner iterations tuned for best wall-clock time; and (3) Subsampled Newton with Hessian Averaging (SN-HA), ie., the method we use in the global phase of Algorithm 1 (without the SVRN phase). All of the convergence plots are averaged over 10 runs. For both SVRN-HA and SN-HA we use Hessian sample size = 4d.
From Figure 2(a), we conclude that as soon as SVRN-HA exits the initial phase of the optimization, it accelerates dramatically, to the point where it nearly matches the rate of classical Newton. This accelera- tion corresponds to the improvement in sequential complexity from O(a log(1/e)) for Stochastic Newton to 0 log((/e) for SVRN. In all our experiments, the transition to the SVRN phase in SVRN-HA occurred very quickly, generally within 1-2 iterations, which indicates that the method casily reaches local convergence. Finally, we observed that the convergence of SVRG is initially quite fast, but over time, it stabilizes at a slower rate, indicating that the Hessian information plays significant role in the performance of SVRN-HA.
log(n)
The code available at https://github: com/ svrneuton/svTn.

---

Published in Transactions on Machine Learning Research (01/2025)
EMSIST 1d:0000
EMNIST ( 1d:1000
Neaton
167
SUTC
22% 11A
10-1
Stitt644
22"4
(a) Convergence
(b) Runtime (seconda)
Convergence and iuntime comparison of SVRN-HA on the EMNIST dataset against three baselines: classical Newton, SVRG (after parameter tuning), and Subsampled Newton with Hessian Averaging (SN- HA), i.., the global phase of Algorithm 1, ran without switching to SVRN. Further results on the CIFAR-10 dataset are in Appendix
Figure 2:
In Figure 2(b), we plot the wall clock time of the algorithms. Here, SVRN-HA also performs better than all of the baselines, despite some additional per-iteration overhead. We expect that this be further optimized. Finally, we note that Newton's method is drastically slower than all other methods due to the high cost of solving a large linear system, and the per-iteration time of SVRG is substantially slowed by its sequential nature.
5.2 Further investigations on least squares task
We next study the setting of least squares regression (4) to analyze the trade-offs convergence for different implementations of SVRN, we vary the gradient and Hessian estimation schemes. We evaluated the algorithms synthetic data matrices, defined in Appendix B.
Communication cost of gradient resampling. Our theoretical analysis requires that for each small step of SVRN, a fresh sample of components vh is used to compute the gradient estimates. However, in Lemma 3 showed that, after variance reduction, the gradient estimates accurate with high probability, which suggests that we might be able to reuse previously sampled components. While this technically does not improve the number of required gradient queries, it can substantially reduce the communication cost for some practical implementations.
As an example, let us consider the setting where each full/mini-bateh gradient computation requires reading the corresponding data chunk from the server onto the computing core. Then, the theoretical version of SVRN requires reading the entire dataset roughly times (2n data points) in the course of one stage (outer iteration): data points for computing the full gradient, and then mmax all of the mini-batch steps together. On the other hand, if we were to reuse the same mini-batch for all of the inner iterations in one stage, then we require reading the dataset only + o(1) times.
(n/m)
data points for
In Figure 3(a), we investigate how much the convergence rate of SVRN-HA is affected by the frequency of component resampling for the gradient estimates. Recall that in all our experiments, we use gradient sample size of nn/ logg(n/d)] We consider the following variants:
Sampling once: extreme policy of sampling one set of components and reusing them for all gradient estimates;
Sampling per stage: an intermediate policy of resampling the components after every full gradient computation.
Sampling per step: the policy which used in our theory, i.c., resampling the gradients every step of the inner loop of the algorithm.
10

---

Published in Transactions on Machine Learning Research (01/2025)
Synthetic Least Squares 500k, d=1000)
Synthatic Least Squares [n 500k, d-1000)
100
SVM HA (sampling once)
Si-HA
10
50:35 MA
3VM HA Gammlligg stage)
Sy10 mA
10
Gompling
10
10
10"4-
10
10
10"h
10.45
10-10
Renakion
Beration
(a) Comparison of three variants SVRN-HA (alongside SN-HA), depending on how frequently we resample data points used to compute the gradient estimate. We consider three variants of SVRN-HA: (1) sampling once for the entire optimization, (2) sampling once for each full gradient stage (per stage), (3) sampling in each small step (per step).
(b) Comparison of SVRN-HA (alongside SN-H HA) against Subsampled Newton with Gradient Subsampling (SNGS- HA), which is implemented exactly like SVRN HA except without the variance-reducing correction.
Figure 3: How different types of gradient estimation affect the convergence propertics of SVRN.
From Figure 3(a) we conclude that, while all three variants of SVRN-HA converge and are competitive with SN-HA, the extreme policy of sampling once leads to substantial degradation in convergence rate, whereas sampling per stage and sampling per step perform very similarly. Thus, overall recommendation is to resample the components at every stage of SVRN-HA, but reuse the sample for the small steps of the algorithm (this is what we used for the EMNIST and CIFAR-10 experiments).
Effect of variance reduction. We next investigate the effect of variance reduction the convergence rate of SVRN. While gradient subsampling has been proposed by many works in the literature on Subsampled Newton (e.g., see Roosta-Khone Khorasani & Mahoney, 2019), these works have shown that the gradient sample size must be gradually increased to retain fast local convergence (which means that after a few iterations, we must use the full gradient). On the other hand, SVRN, instoad of incrcasing the gradient sample size, we use variance reduction with a fixed sample size, which allows us to retain the acceleratod convergence indefinitely.
To illustrate this point, in Figure 3(b) we plot how the convergence behavior of our algorithm changes if we take variance reduction out of it. The resulting method is called Subsampled Newton with Gradient Subsampling (SNGS-HA). For this experiment, we resample the gradient estimate at every small step (for both SNGS-HA HA and SVRN-HA). For the sake of direct comparison, all of the other parameters are retained from SVRN-HA. In particular, one itcration of SNGS-HA corresponds to [nog)(n/d) steps using resampled gradients, and Hessian averaging occurs once every such iteration. As expected, we observe that, while initially converging at fast rate, eventually SNGS-HA reaches point where the subsampled gradient estimates are not sufficiently accurate, resulting in a sudden dramatic drop-off in the convergence rate, to the point where the method virtually stops converging altogether. On the other hand, SVRN-HA continues to converge at the same fast rate throughout the optimization procedure without any reduction in performance. This indicates that variance reduction does improve the accuracy of gradient estimates, especially when our goal is to converge lo a high-precision solution.
Effeet of Hessian accuracy. In our experiments, for both SVRN and SN, we used Hessian averaging (Na et al., 2022) to construct the Hessian estimates. This approach is desirable practice, since gradually increases the accuracy of the Hcssian estimate as we progress the optimization. As a result, it is more robust to the Hessian sample size and we are guaranteed to reach sufficient accuracy for SVRN work well. In the following experiment, we take Hessian averaging out of the algorithms provide better sense of how
11

---

Published in Transactions on Machine Learning Research (01/2025)
Synthatic Least Squares ( 5==000,, d=1000)
High Cohererce Least Squares n- 500k, d=1003)
109
(h 1681
160
SN--A
SVANN It=10k]
S6/0N MA
sPN (b=20k)
50V HA WH1
10
SVAR [h=20kd
10:2
SN (h- 30kl
SURN (h-30k]
10"4
10t4
10
10"1
10-20
10-10
Reration
Reration
(a) Convergence comparison SVRN and SN using fixed Hessian estimates (i.e., without Hessian averaging). Here, denotes the number Hessian samples used to generate the estimate.
(b) Convergence comparison of SVRN-HA and SN-HA, with and without preconditioning using Randomized Hadamard Transform (RHT), for a high-coherence least squares dataset.
Figure 4: How Hessian sample size and data coherence affect the convergence properties of SVRN.
the performance of SVRN and SN depends on the accuracy of the provided Hessian estimate. For simplicity, we focus here on least squares, where the Hessian the same everywhere, so we can simply construct an initial Hessian estimate and then it throughout the optimization. However, our insights apply more broadly to local convergence for general convex objectives. In Figure 4(a), we plot the performance of the algorithms we vary the accuracy of the subsampled Hessian estimates, where h denotes the number of samples used to construct the estimate. In all the results, we keep the gradient sample size and local steps in SVRN fixed as before.
Remarkably, the performance of SVRN is affected by the Hessian accuracy very differently than SN. We observe that SVRN requires certain level of Hessian accuracy provide any acceleration over SN. As soon this level of Hessian accuracy is reached (by increasing the Hessian sample size h), the peformance of SVRN quickly jumps to the fast convergence the accuracy no longer appears to affect the convergence rate of SVRN. This is in contrast to SN, whose convergence slowly improves as we increase the Hessian sample size. This intriguing phenomenon is actually fully predicted by our theory for SVRN (together with prior convergence analysis for SN). Our convergence result (Theorem 3) requires a sufficiently accurate Hessian inverse estimate for SVRN to work with a unit step size (which is what is used SVRN-HA), but the actual rate of convergence independent of the Hessian accuracy (only the required number of small steps affected). We conclude that SVRN is more desirable than SN when we have small budget for Hessian samples.
observed in the other experiments. Further increasing
Effect of high coherence. We next analyze the performance of SVRN and SN on a slightly modified least squares task. For this experiment, we modify the data matrix A, by multiplying the ith row by 1/VTT for each i, where 9i is an independent random variable distributed according to the Gamma distrution with shape and scale 1/2. This high leverage score. Recall that the leverage score of the ith row of A defined as l = (AT A)-laa see Appendix A.4. This can be viewed as affecting the component-wise smoothness of the objective, which hinders subsampling- based estimators of the Hessian and the gradient.
standard transformation designed to produce a matrix with many rows having
In Figure 4(b), we illustrate how the performance of SVRN-HA and SN-HA degrades for the high-coherence least squares lask, and we also show how this can be addressed by relying on the ideas developed in Section 2.2 (and further discussed in Appendix A.4). First, notice that not only is the convergence rate of both SVRN- HA and SN-HA worse on the high-coherence dataset than on the previous least squares examples (c.g., compare with Figure 3(b)), but also, the acceleration coming from variance reduction is drastically reduced to the point of being negligible. The former effect is primarily caused by the fact that uniform Hessian
12

---

Published in Transactions on Machine Learning Research (01/2025)
subsampling is much less effective at producing accurate approximations for high-coherence matrices, and this affects both algorithms similarly (we note that one could construct an even more highly coherent matrix, for which these methods would essentially stop converging altogether). The latter effect is the consequence of the fact that gradient subsampling is also adversely affected by high coherence, so it becomes nearly impossible to produce gradient estimates with uniform sampling that would lead to an accelerated rate, even with variance reduction. This corresponds to the regime of k 2 n our theory.
Fortunately, for least squares regression, this phenomenon one of two strategies: (1) use importance sampling proportional to the leverage
be addressed easily. As outlined in Ap-
pendix A.4, we scores of A for both the Hessian and gradient estimates; or (2) precondition the problem using the Random- ized Hadamard Transform (RHT) to uniformize all the leverage scores, and then use uniform subsampling, Both of these methods require roughly Ond n) preprocessing cost and eliminate dependence on the con- dition number for both SVRN-HA and SN-HA. The latter strategy somewhat more straightforward since it does not require modifying the optimization algorithms, and we apply it here for our high coherence least squares lask: we let SVRN-HA-RHT and SN-HA-RHT denote the optimization algorithms ran after applying the RHT preconditioning to the problem. Note that this not only improves the convergence rate of both methods but also brings back the accelerated rate enjoyed by SVRN-HA in the previous experiments. In fact, our least squares results (Theorem 2 and Lemma 6) can be directly applied to SVRN-HA-RHI so this method and its accelerated convergence rate of ((d/n) is provably unaffected by any high-coherence matrices.
6 Conclusions
We propose and analyze Stochastic Variance-Reduced Newton (SVRN), provably effective stralegy of incorporating variance reduction into popular stochastic Newton methods for solving finite-sum minimization tasks. We show that SVRN improves the local convergence complexity of Subsampled Newton (per data pass) from (( log(//)) to O
log(x) ), while retaining all the benefits of second-order optimization, such as simple unit step size and easily scalable batch operations.
References
Nir Ailon and Bernard Chazelle. The fast Johnson Lindenstrauss transform and approximate nearest neigh- bors. SIAM Journal on computing, (9):)3302 -22, 2009.
Zeyuan Allen-Zhu. Katyusha: The first direct accoleaation of stochastic gradient methods. The Journal of Machine Learning Research, 18(1):8194 8244, 2017.
Zeyuan Allen-Zhu and Yang Yuan. Improved svrg for non-strongly-conee In International conference on machine learning, pp. 1089. PMT.R, 2016.
sum-of-non-comeee objectives.
Haim Avron, Petar Maymounkov, and Sivan Toledo. Blendenpik: Supercharging lapack's least-squares solver. SIAM Journal on Scientífic Computing, 3217 2010.
Albert S Berahas, Raghu Bollapragada, and Jorge Nocodal. An investigation of newton-sketch and subsam- pled newton methods. Optimization Methods and Software, 35(4):661-66) 680, 2020.
Raghu Bollapragada, Richard H Byrd, and Jorge Nocedal. Exact and inexact subsampled newton methods for optimization. IMA Journal of Numerical Analysis, 39(2), 2018.
Stephen Boyd and Lieven Vandenberghe. Conver optimization. Cambridge university press, 2004.
Shabarish Chenakkod, Michal Dereziński, Xiaoyu Dong, and Mark Rudelson. Optimal embedding dimension for sparse subspace embeddings. In 56th Annual ACM Symposium on Theory of Computing, 2024.
Kenneth L. Clarkson and David Woodruff. Low-rank approximation and regression input sparsity time. J. ACM, 63(6):54:1-54:40 January 2017.
13

---

Published in Transactions on Machine Learning Research (01/2025)
Gregory Cohen, Saced Afshar, Jonathan Tapson, and Andre Van Schaik. Emnist: Extending mnist to handwritten letters. In 2017 international joint conference on neural networks (IJCNN), pp. 2921-2926. LEEE, 2017.
Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradient method with support for non-strongly convex composite objectives. Advances in neural information processing systems, 27, 2014.
Michal Dereziński and Michael W Mahoney. Distributed estimation of the inverse hessian by determinantal averaging. In Advances in Newral Information Processing Systems 32, pp. 11401-11411. Curran Associales, Tnc., 2019.
Michal Dereziński and Michael W Mahoney. Determinantal point processes in randomized numerical linear algebra. Notices of the American Mathematical Society, 68(1):34-45, 2021.
Michal Dereziáski and Michael W Mahoney. Recent and upcoming developments in randomized numerical linear algebra for machine learning. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 2470-6479, 2024.
Michal Dereziússki, Jonathan Lacotte, Mert Pilanci, and Michael W Mahoney. Newton-LESS: Sparsification without trade-offs for the sketched newton update. Advances in Neural Information Processing Systems, 34, 2021.
Michal Dereziński, Zhenyu Liao, Edgar Dobríban, and Michael W Mahoney. Sparse sketches with small inversion bias. In Proceedings of the 34th Conference Learning Theory, 2021.
Petros Drineas and Michael W. Mahoney, RandNLA: Randomized numerical linear algebra. Communications of the ACM, 59:80-90, 2016.
Petros Drincas, Michael W Mahoney, and S Muthukrishhan. Sampling algorithms for &g regression and applications. In Proceedings of the seventeenth anmual ACM-SIAM symposium on Discrele algorillm,, pp. 1136, 2006.
Petros Drineas, Malik Magilon-Tsmail, Michael W. Mahoney, and David P. Woodruff. Fast approximation of matrix coherence and statistical leverage. J. Mach. Learn. Res., 13(1):3475-3506, December 2012.
Murat A Erdogdu and Andrea Montanari. Convergence rates of sub-sampled Newton methods. Advances in Neural Information Processing Systems, 28:3052-3060 30600 2015.
Gerald Folland. Advanced calculus. Upper Saddle River, NJ Prentice Hall, 2002.
Roy Frostig, Rong Ge, Sham M Kakade, and Aaron Sidford. Competing with the empirical risk minimizer in single pass. In Conference on learning theory, pp. 728-763. PMLR, 2015.
Alon Gonen, Francesco Orabona, and Shai Shalov-Shwarzz. Solving ridge regression using sketchod precon- ditioned svrg. In International conference on machine learning, pp. 1397 1405. PMLR, 2016.
Robert Gower, Donald Goldfarb, and Peter Richtárik. Stochastic block bfgs: Squeczing more curvature out data. In International Conference on Machine Learning, pp. 1869-1878. PMLR, 2016.
Robert Gower, Nicolas Le Roux, and Francis Bach. Tracking the gradients using the hessian: A new look at variance reducing stochastic methods. In International Conference on Artificial Intelligence and Statistics, pp. 707-715. PMLR, 2018.
Vipul Gupta, Avishek Ghosh, Michal Dereziáski, Rajiv Khanna, Kannan Ramchandran, and Michael W Mahoney, Localnewton: Reducing communication rounds for distributed learning. In Uncerlainty in artificial intelligence, pp. 642-642. PMLR, 2021.
Robert Jerrard. Multivariable calculus. notes/Contents. html, 2018. Course Notes.
toronto. htt.:atttttttttttttoooeeeeeeeeee htt.:///oooees//a/33/////////
14

---

Published in Transactions on Machine Learning Research (01/2025)
Ric Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Aduances in neural informalion processing systems, 26, 2013.
Jakub Konecný, Jie Liu, Peter Richtárik, and Martin Takác. Mini-batch semi-stochastic gradient descent. the proximal setting. IEEE Journal of Selected Topics in Signal Processing, 10(2):242 255, 2015.
Dmitry Kovalev, Konstantin Mishchenko, and Peter Richtárik. Stochastic newton and cubic newton methods with simple local lincar-quadratic rates. orXiv preprint arXi::1912.01599 2019.
Hongzhou Lin, Julien Mairal, and Zaid Harchaoui. A universal catalyst for first-order optimization. Advances in neural information processing systems, 28, 2015.
Yanli Liu, Fei Feng, and Wotao Yin. Acceloration of svig and katyusha x by inexact preconditioning. International Conference on Machine Learning, pp. 4003 44012. PMLR, 2019.
Xiangrui Meng, Michael A Saunders, and Michael W Mahoney. LSRN: A parallel iterative solver for strongly over-or underdetermined systems. SIAM Journal on Scientific Computing, 36(2):C95-C11) 2014.
Stanislav Minsker. On some extensions of bernstein's inequality for self-adjoint operators. Statistics & Probability Letters, 127:111-119, 2017.
Aryan Mokhtari, Mark Fisen, and Alejandro Ribeiro. Tqn: An incremental quasi-newton method with local superlinear convergence rate. SIAM Journal on Optimization, 28(2):1670 1698, 2018.
Philipp Moritz, Robert Nishihara, and Michael Jordan. A linearly-convergent stochastic I-bfgs algorithm. In Arlificial Intelligence and Stalistics, pp. 249 258. PMLR, 2016.
Riley Murray, James Demmel, Michael W Mahoney, Benjamin Erichson, Maksim Melnichenko, Osman Asif Malik, Laura Grigori, Piotr Luszezek, Michal Dereziński, Miles E Lopes, al. Randomized numerical linear algebra: A perspective on the field with an eye to software. arXiv preprint arXiv:2902.11 11474, 2023.
Sen Na, Michal Dereziński, and Michael W Mahoney. Hessian averaging in stochastic newton methods achieves superlinear convergence. Mathematical Programming, 2022.
Jelani Nelson and Huy L Nguyēn. Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings. In 2013 feee 54th annual symposium on foundations of compuler science, pp. 117 126. IEEE, 2013.
Mert Pilanci and Martin Wainwright, Iterative hessian sketch: Fast and accurate solution approximation for constrained least-squares. The Journal of Machine Learning Research, 17(1):1442-1879, 2016.
Mert Pilanci and Martin Wainwright. Newton sketch: A near linear-time optimization algorithm with linear-quadratic convergence. SIAM Journal on Optimization, 27(1):205 245, 2017.
Anton Rodomanov and Dmitry Kropotov. A superlinearly-couvereett proximal newton-type method for the optimization of finite sums. In International Conference Machine Learning, pp. 2597 2605. PMLR, 2016.
Vladimir Rokhlin and Mark Tygert. A fast randomized algorithm for overdetermined linear least-squares regression. Proceedings of the National Academy of Sciences, 105(36):13212-1332) 2008
Farbod Goosta-Khorasan and Michael W Mahoney. Sub-sampled newton methods. Mathemalical Program- ming, 174(1):293 326, 2019.
Nicolas Roux, Mark Schmidt, and Francis Bach. A stochastic gradient method with an exponential conver- gence rate for finite training sets. Advances in neural information processing systems, 25, 2012.
Tamas Sarlos. Improved approximation algorithms for large matrices via random projections. Proceedings of the 47th Annual IEEE Symposium on Foundations of Computer Science, FOCS '06, pp. 143-152. Washington, DC, USA, 2006. IEEE Computer Society.
15

---

Published in Transactions on Machine Learning Research (01/2025)
Shai Shalev-Shwaie and Tong Zhang. Stochastic dual coordinate ascent methods for regularized loss mini- mization. Journal of Machine Learning Research, 12(,), 2013.
Nilesh Tripuraneni, Mitchell Stern, Chi Jin, Jeffrey Regier, and Michael Jordan. Stochastic cubic regular- ization for fast nonconvex optimization. Advances in neural information processing systems, 31, 2018.
Joel A. Tropp. User-friendly tail bounds for sums of random matrices. Foundations of Computational Mathemalics, 12(4))389 434, August 2012.
Xiao Wang, Shiqian Ma, Donald Goldfarb, and Wei Liu. Stochastic quasi-newton methods for nonconvex stochastic optimization. SIAM Journal on Optimization, (7 2017.
David P Woodruff. Sketching as a tool for numerical lincar algebra. Foundations and Trendsfy in Theoretical Computer Science, 10(1-2):1-177, 2014.
Jiyan Yang, Yin-Lam Chow, Christopher Ré, and Michael W Mahoncy. Weighted sgd for Ip regression with randomized preconditioning. The Journal of Machine Learning Research, 18(1):7811 7853, 2017.
Junyu Zhang, Lin Xiao, and Shuzhong Zhang, Adaptive stochastic variance reduction for subsampled newton method with cubic regularization. INFORMS Journal on Optimization, 4(1):45-64 64, 2022.
Dongruo Zhou, Pan Xu, and Quanquan Gu. Stochastic variance-reduced cubic regularization methods. J. Mach. Learn. Res., 20(134):1 47, 2019.
A Proofs for local convergence analysis of SVRN
In this section, we provide the proofs of our main technical results, i.., local convergence analysis for SVRN. First, we prove the result for the general case (Theorem 3), then (Theorem 2).
prove the result. for least squares
A.1 Preliminaries
First, let us recall the formal definitions of the standard Hessian regularity assumptions used in Theorem 3. For all our results, it is sufficient that the function f satisfics either one of these assumptions.
Assumption 2 Function f Rd + R has Lipschitz continuous Hessian with constant L, i.e., [V s2(((() s T. ||x x'Il for all x, x' E Rdd.
Assumption Function f Rdd. > is self-concordant, i.e., for all x, x' e Rd, the function s(t) = (x++xx) satisfies: (v) < 2(0"
In the proof, lary 4.1 in Minsker, 2017).
the following version of Bernstein's concentration incquality for random vectors (Corol-
Lemma 1 Let Vi,...; Vm E Red be independent random vectors such that E[vs] 0 and ||v/|| S R almost surely. Denote o² ;= >m [|v:] Then, for all f2 2 o² + tR/3, we have
m
< 28 exp (- 0² o2+tR/3)
f2/2
Pr
V;
We also use the following lemma convert from convergence in the norm, ||x x x"I to convergence excess loss,
in the neighborhood around the optimum x". The proof of this lemma, given in Appendix A.3, uses Quadratic Taylor's Theorem.
Lemma 2 If satisfies Assumption 1 and either Assumption 2 or 3, then for any E [0, 1] and x E Uy(a),
we have:
V2 F2f(x) Re, V2f(x") and f(x) - f(x") ~e, ||x x" x2
16

---

Published in Transactions on Machine Learning Research (01/2025)
A.2 Proof of Theorem 3
To simplify the notation, we will drop the subscript 8, 80 that X = X, and g g.. Also, let us define g(x) 1 Ee ii (x). We use g(x) = Vx(x,, gr = g(x)), H, = fxx), H = fxx(xx), &. = g(x)), and g. g g(x) + as shorthands. Also, we will use A = x x" We start by splitting up the error bound into two terms: the first one is an error term that would arise if we were using the exact gradient ge instead of the gradient estimate ge; and the second term addresses the error coming from the noise the gradient estimate. Initially, we use the error || analyze the convergence rate in one step of the procedure, where recall that ||v||M / Mv. We then convert that to get convergence in function value via Lemma 2.
We first address the assumption that H is an a-approximation of H,, as defined by the condition (2). Note that, via Lemma 2, for any X₁ E (f(er) we have that H, Re, H, which for a sufficiently small el implies that in the sense of (2). This, in turn implies that is a I.lo-approximation
H. is a L.laaprooximation of H. because H N vaH, I VI.laH (the other direction is analogous). For the sake of simplicity, we will replace 1.1a with a and say that H an a-approximation of H (this can be casily accounted for by adjusting the constants at the end).
Now, suppose that after inner iterations, get x, (e() satisfying ||el| S | [H. Our decomposition of the error into two terms procceds as follows, where we use P.- = H-GHH ge:
| H |1(x; r/p:) x H
= || - 1 g: + 1/ H g: - 1/H-'B.ll
A "
g: |||
+ + H (g: -g:)]
(6)
To bound the second term in (6), we first observe that H- 3 VoH-- -1 which in turn yields -22
Val. I. Thus, we can write H 1 /1 S va and we get:
H gr g:) H H 1/ H
H
1/2 gr -g))
H
11 | 1/2 (g: ||
Va |gg: 8:|
We now break l|g: -grlle 1- down into two parts, introducing g(x*) and separating g(x) from g::
g:
grl H-1 ||g: (g: [ 3) H
S g: - (8. g(xx)) + - (g(x) (Ixx))
We bound the above two terms using the following lemma, which gives a new high-probability error bound for the variance reduced gradient estimates in the large mini-batch regime which, unlike results from prior work that hold in expectation, crucially relies on the iterate being in the local neighborhood (1) around the optimum x".
Lemma 3 There is an absolute constant C > such that for any x (L), letting H V² f(x)), the gradient estimate g(x) m Em (x), using m k xlog(1/)) samples, with probability 1 5 satisfies:
I (x) g( V (x) H- C'loo(1/0) K |x -
Proof We will apply Bernstein's concentration inequality for random vectors (Lemma 1) to v, Vx((x) (x)x(x)) [x()). First, observe that E vxx(x) Vx(x) and E vxx((")) vx(x*) = 0, so in particular,
E[v.] = 0.
In the next step, we will use the fact that for any X-smooth function 9. we have [| S 2X (g(x) minx 9(x'),, which follows because:
min x
gxx) g(xx Vx()))
1 X |1(2))1
A 1
S g((x)
[xx))]
2 x2
1
-x() 2X [x)()))
17

---

Published in Transactions on Machine Learning Research (01/2025)
We will use this fact once on f, and also second time, the function g(x) xix(x) xx(xx*) - (x
x²) "Vv((""" which is X-smooth because the is X-smooth, observing that Vg(x) = VVd((x) - Vole(x*) and that minx g(x') g(x*) = O. Thus, we have
||v:|| < ( +
1 4A (vi - -x((x') - (x-xxx) Vxxxx))V(((x')) + +4A (J(x)
- 2x - x2I² + 2A |x 4x-xxx1==>>|1xx--- 1x
where in the last step we used again that v4 and are X-smooth. To bound the expectation E [vss]?, we use the intermediate inequality from the above derivation, obtaining:
v. 1|2 = B V b/ (x) Vel(x - 2 V (x) [oI(x)] VVxx)
= B V by (x) Vxl(x - [x)]1
< [[| V (x) (x)(x)
s E 2A (v+(x) <x(x") (x x²) (x)(()))
2A (f(x) -(x*) (x (xxxx)V[(x'))
2A (f(x) - -(xx))
We now use the assumption that x E (1), which implies via Lemma that f(x)- - -xxx'xxx)+>>> 22-
||x xx-x"1li I Thus, we can use Lemma with R = 2A|[x x*1l and o² 2mAl/x x" ll as well p-stroog convexity of f, obtaining that, for some absolute constant C, with probability 1 5, we have:
o xx) - [(x") - fx(x)
0² log(1/0)
R log (1/0)
m²
m²
C
2A||x -
x* log(1/0)
+ -
-x" 1 2x"11"1o''
(116)
/
m
m2
k: /llog(1/0)
+ x²2 log? (1/6)
< 4C
|x - 1xx-x116
m²
k ||x x 1|
< 8C log(1/8)
where in the last step we used that m log(1/6))
Letting ===vo (2C/log(/)))//m. Lemma 3 implies that with probability 1 1 -8/F? ==22C'lo((f
Ilge (g. g(x) g) S eg (1A.l| + [|AolHH) S 2egl Aol|H.
Finally, return to the first term in (6), i.c., 14 n 'g I|.- To control this term introduce the following lemma which is potentially of independent interest to the local convergence analysis of Newion-type methods.
Lemma 4 Suppose that f satisfies Assumption and either one of the Assumptions or 3, and take any x x U((()) (see Definition 1) for 4 S 1/ca for a sufficiently large absolute constant > Let H f(x) and consider a pd malria H that salisfies 1 H H H. Then, for == min{ 22/a, 1}, we have: Va
xx 1x-18 1x-7 - - x S (1 1.9a ||x x" I|
18

---

Published in Transactions on Machine Learning Research (01/2025)
Proof Let Au x x" and A1 := (Vfxx) x'. Using that Vf(x*) 0, have:
A1 = Au 7 VV(x)
Ao nH (V)(xx) vVf(x)))
= A n //" v2 v2
-Do))ad0
= [I n AH
where we defined H : So' 42 v2 f(x* + 0 d0. follows that can bound the norm of A using norm defined by the matrix H:
Ailla | 7 A)
= =I - 1 / H (2 12
$ | 11-1HH122 1/2 " 1/2 1
Observe that for any 0 E [0, 1], the vector x" + +0o belongs to Uy(e), which via Lemma 2 implies that
+ Re H v0 E [0.
where H v²ʃ(x*. In particular, this means that H Rig H, which, combined with the a-approximation property of H, allows us to write the following:
H VaH 3 va(1 H and H- = va
H- H
Putting these inequalities together, we obtain that:
a I nH HH' 1/2 =nva(1+
Va
Now, using the fact that 17 = m 2/a, 1}, we conclude that:
I - H 1/ & max V1 + + - 1, 1 - El
Va
V2(1 4))
1 - 1
{V2(1+ +e) 11,1
max
V2 ʃ
a
(1
1
1
< max
1
<
1.8
1.8x
where used that a < 1/c for sufficiently large constant such that max 22(1++++)))))) 1,, -1,1 < Ig. Now, analyse convergence in the norm induced by H, instead of H, by relying again the fact that H mqg H. obtaining:
v2
v1 Ailla I 1 (1- _1) |
| H
-1-q - a 1.8a
q (1- (1-18a) El
) |loll S
1 _ l| H:
where the last step requires e; 1/cor for sufficiently large absolute constant
Using Lemma 4 to bound the first term in (6), we obtain that:
| n, |1 < (1 1.9a
19

---

Published in Transactions on Machine Learning Research (01/2025)
Putting everything together, we obtain the following bound for the error of the update that uses the stochastic variance-reduced gradient estimate:
|
|4,
7/PrlH
||
11 || + 17 1 (g: g
| H + 21/ / aEe || |H
1.9a
|A,|| 3eg ||
1.9a
Note that, long as 3eg 2 L (which can be ensured by our assumption on m), this implies that || lla 2a
||Ao|| H and so X (e(e)). Thuss our analysis be applied recursively at each inner ileration. To expand the error recursion, observe that if we apply union bound over the high-probability events in Lemma 3 for each inner iteration using failure probability & 5/t², then they hold for all with probability at least - - 8/1² 2 d We obtain:
|A (1 1.9a
)1 A 1|H 3eg Aoo1|Ao||H
1-1
1.9a) ) I|oollH
Σ (1-1.90))
3eg Aoll
1 _111
+ 9acg) ||AollH.
1.9a
Applying Lemma 2, we convert this to convergence in function value:
1
f(x;) - f(x))
A
11
- 6
2
$ 1
1
(1
1
+ +90rg
A 2
q 2
1.9a
+
(1-
1
+ gg
(f(xo) - fxx"))
El
1.9a
(1 2a
)
C'a
2 log(t/8)
(f(x)) f(xx)),
where C" is an absolute constant, and we again used that S 1/ea for a sufficiently large e, thus concluding the proof.
A.3 Proof of Lemma 2
First, we show that the Hessian at x E (f(e)) is an q-approximation the Hessian at the optimum x". This is broken down into two cases, depending on which the two Assumptions and are satisfied.
Case 1: Assumption 2 (Lipschitz Hessian). Using the shorthand convexity (Assumption 1) implies that 42 3²f(x) 2 μIl, we have:
V fxx(x" and the fact that strong
H f(x - HJH 1/2 = - [ P
L
H|l s x - x
xx x"Il S El
13/22
which implics that f2f(x) Re H.
Case 2: Assumption 3 (Self-concorlancee. The fact that v²f(x) Rey H follows from the following property of self-concordant functions (Boyd & Vandenberghe, 2004, Chapter 9.5), which holds when ||x x" x"Il < 1:
(1 xx - x" 11 H V2 f(x) = (1 - x|x - x" (H) -2) H.
where we again let H fx*))
20

---

Published in Transactions on Machine Learning Research (01/2025)
We next use version of Quadratic Taylor's Theorem, as given below. See Theorem 3 in Chapter 2.6 of (Jerrard, 2018) and Chapter 2.7 in (Folland, 2002).
Lemma 5 Suppose that f : Rd - lK has continuous first and second derivatives. Then, for any a and v, there exists E (0,1) such that:
f(a +a+v) v) [ + VV((a) v 1
+ VVTV2/(a Ov)v.
Applying Talyor's thoorem with a = x* and v x", there is a z = + - such that:
- 2 - xx"xx"p x
f(x) f(x")
where we use that fx(") = Since we assumed that x E Uf()),, and naturally also x* E EUy(a). this means that Us(ez), given that Uy(ez) is convex. Thus, using that 2 (z) Ra H, we have ||x x (()))))) Re
||x x" I
A.4 Proof of Theorem 2
In this section, we discuss how the convergence analysis of SVRN can be adapted to using leverage score sampling when solving a least squares lask (proving Theorem 2).
apy byy and I is an index from [I,...,n) sampled according to some importance sampling distribution More specifically, consider a least squares task, where the components are given by (>x = 1 (a x - ni)2. Then, the overall minimization task becomes:
Consider an expected risk minimization problem f(x) = [x)(, where & = _
1 n
Bxx+ = LE 19-- w (x)] = Σ
(a x y1)
(7)
rppT
2n
1=1
Moreover, we have f(x) [(x*) 2/ 11 x" x1)11= = 2/1l1xx -
x* 1 where H (2](x- LATA. Also,
Vd((x) = (a. y1)a
v'vx(x)= = a
Naturally, since the Hessian is the same everywhere for this task, the local convergence neighborhood Uy is simply the entire Euclidean space Rd Let us first recall our definition the condition number for this task. Assumption 1 states that each vi is X-smooth, ie., (x)]l, = ||a/| < À and f is p-strongly convex, i.c., Amiu (H) o miinnn 2 μt and the condition number of the problem is defined as = A/μ Z max((n||a.|)))) o (A). Can we reduce the condition number of this problem by importance sampling?
Consider the following naive strategy which can be applied directly with convergence result. Here, we let the importance sampling probabilitics be & ||, so that the smoothness of the new reweighted problem will be Ā = _ E ||a;|| -2 In other words, it will be the average smoothhess of the original problem, instead of the worst-case smoothness. Such importance sampling strategy can theoretically be applied to a general finite-sum minimization task with some potential gain, however we may need different sampling probabilities at each step. For least squares, the resulting condition number is ă X/μ = (E, [|a.l3// o (A). This is still worse than what claimed for least squares, but it is still potentially much better than
Next, we will show that by slightly adapting our convergence analysis, we can use leverage score sampling to further improve the convergence of SVRN for the least squares task. Recall that the ith leverage score of A is defined as 4, = ||a:l/AT -)- = || l and the laverage scores satisfy Σ ls = d. This result will require showing a specialized version of Lemma 3, which bounds the error in the variance-reduced subsampled gradient. In this case we show bound in expectation, instead of with high probability.
Lemma 6 Suppose that f defines a least squares task (7) and the sampling probabilities satisfy p. 2 || l2 T AC)=:(CCd). Then g(x)= = m Em opI, (x),()), where II,... Im ~ p, satisfies:
E [g) g((x)) - f(x) - sc d 1|x x" 1I. m
21

---

Published in Transactions on Machine Learning Research (01/2025)
Proof We define V,
_1_ (Vii, (x) - (x)((x))) - Vx(x).
Note that E[v,] 0, so we have:
upr;
m
1
(x) -g(x*) - Vf(x)|
BlVi
m
H
m
1
<
E
V4l,, (x) Vx,,("))
m
P
11
1
|la,, 2
-E
(a), (x -x'))²
m
n n²pi,
2
1
(aj, (x x²))
d
< m
Cd E
C x
x l|
np,,
m
where we used that ||a:|| = n ||aill/A ATAATAAT AJ
S Cndpi.
Since the above bound obtained in expectation, to insert into our high probability analysis, we apply Markov's inequality. Namely, it holds with probability 1 - 5 that:
Cd
[ - ((x))
S
6m ||x
x|x-x x
Compared to Lemma 3, the dependence on the condition number ss is completely climinated in this result. Letting m = n/ log(n/d) and the number of local iterations of SVRN to be = Ollo(n/dd), we can apply the union bound argument from the proof Theorem 3 by letting 8 = 1/(Ct), so that with probability 1/C, one stage of leverage score sampled SVRN satisfies:
0 p
f(x,+1) f(x)) S p. (f(x))-f(e
for
p p=
n
Alternatively, our main convergence analysis can be adapted (for least squares) to convergence in expectation, obtaining that B[f(x++1) - fx==) < ý (x,) - f(x")] for p p= (oggn///d))
The time complexity stated in Theorem 2 comes from the fact that constructing a preconditioning matrix H that is an a-approximation of H with a = O(1,, together with approximating the leverage scores, takes (nd log n + d8 log d) (Drineas et al., 2012), whereas one stage of SVRN takes O(nd+ dd log(n/d)). Here, the preconditioning matrix can be formed by applying a k sketching transformation S to the data matrix A, and then computing the Hessian estimate AT & H. For example, if use the Subsampled Randomized Hadamard Transform (SRHT, Ailon & Chazelle, 2009), then it suffices to use k = O(d log d). Finally, the initial iterate Xn can be constructed using the same sketching transformation via the so-called Aketch-and-sole technique (Sarlos, 2006):
X argmin |AAx -
With à: = Od d), this initial iterate will satisfy fxx) < ((1)-f(x*), so the mumber of iterations of SVRN needed to obtain f(x,) (1 c))(x*) is only 0 lug (lg(Lc))
log (x//))
We note that another way to implement SVRN with approximate leverage score sampling is first precon- dition the entire least squares problem with Randomized Hadamard Transform (i.c., SRHT without the subsampling):
= HDA
and
y = HDy,
(8)
where H is a Hadamard matrix scaled by 1/vū and D a diagonal matrix with random sign entries. This is popular technique in Randomized Numerical Linear Algebra (Woodruff, 2014; Drineas & Mahoney, 2016; Dereziński & Mahoney, 2021). The cost of this transformation is O(nd log n)), thanks to fast Fourier transform techniques, and the resulting least squares task is equivalent to the original onc, because ||Ax - - |A y for all x. Moreover, with high probability, all of the leverage scores of Á are nearly uniform, so, after this preconditioning, we can simply implement SVRN with uniform gradient subsampling and still enjoy the leverage score sampling when A is a dense matrix, but is less effective when we want to exploit data sparsity.
convergence rate from Thcorem 2. This strategy is as effficient as direct
22

---

Published in Transactions on Machine Learning Research (01/2025)
Synthotic Logistic Regression n- m00k d-1000
Synthetic Logistic Regression n=-00k, d-1000
Newton
Newton
100
107
SVhN-1A
SVNN--
10
10
30"1
10
30-3
10
20*4
30
Beration
Iterator
(a) Synthetic LR, A
(b) Synthetic LR, A 10
FI d=500
CIFAR-1D (d=1000)
Newton
30a
Nowton SNHHA SVAN 4A
34-W4
SVAN
20-4
10
20-1
10
20
10
10
Berral Fon
(c) CIPAR-10, d = 500
(d) CIFAR-10, d 1000
Figure 5: Convergence comparison of SVRN-HA against SN-HA and Newton for a synthetic logistic regression task as we vary the condition number of the data matrix, and for the CIFAR-10 dataset.
B Further experimental details
In this section we provide additional details regarding our experimental setup in Section 5, as well as some further results on logistic regression with several datasets.
As dataset, in Section 5, we used the Extended MNIST dataset of handwritten digits (EMNIST, Cohen et al., 2017) with 500k datapoints. Here, we also include results on the CIFAR-10 image dataset with = 50k datapoints. Both datasets are preprocessed in the same way: Each image is transformed by random features map that approximates a Gaussian kernel having width 0.002, and we partitioned the classes into two labels 1 and -1. We considered two feature dimensions: d = 500 and d = 1000, and we used the regularization parameter 7 10-8, To measure the error in the convergence plots, we use ||x; - x* / x² 12, where H = fxx(xx).
We next present further results, studying the convergence properties of SVRN on synthetic datasots with varying properties, for the logistic regression task as in (5). To construct our synthetic data matrices, we first generate an d Gaussian matrix G, and let = UDV be the reduced SVD of that matrix (we used 500k and d = 1000). Then, we replace diagonal matrix D with matrix D that has singular values spread linearly from to KA- We then let A = UDV be our data matrix. To generate the vector for logistic regression, we first draw a random vector x N(0, 1/d Ia), and then we let y = signn(x))
For the least squares tasks Section 5.2 we generated the same synthotic matrices A, but with the vector y generated as follows: y = Ax+ & where & the Gaussian noise & ~ A/ 0.1 I,). Here, we observed little difference in convergence behavior when varying K (we show the results for h^ 103).
23

---

Published in Transactions on Machine Learning Research (01/2025)
Logistic regression with varying condition number. To supplement the EMNIST logistic regression experiments in Section 5, we present convergence of SVRN-HA on the CIFAR-10 dataset, as well as for the syuthetic logistic regression task while varying the squared condition number x² of the data matrix. Note that, while A is not the same as the condition number of the finite-sum minimization problem, it is correlated with it, by affecting the convexity and smoothness of J. From Figure 5, we observe that SVRN-HA outperforms SN-HA for both values of the data condition number. However, the convergence of both algorithms gets noticeably slower after increasing A while it does not have as much of an effect on the Newton's method. Given that the increased condition number affects both methods similarly, we expect that the degradation in performance is primarily due to worse Hessian approximations, rather than increased variance in the gradient estimates. This may be because we are primarily affecting the global convexity of f. as opposed to the smoothness of individual components 4/. See our higl-coherence least squares experiments for a discussion of how the smoothness of component functions affects the performances of SVRN and SN very differently.
C Related work on Subsampled Newton
In this section, we discuss several important prior works on Subsampled Newton methods to put our results in context. Specifically, we aim to illustrate how the Hessian approximation condition used in Theorems 1 and 3, i.e., J/2 (x) H 3 vaV²f((x), relates to the Hessian estimates used in this line of works when showing fast local convergence rates. Also, in Appendix D.2, we show that to recover our condition with 2 via uniform Hessian subsampling, one needs (x log(d)) samples. Throughout this section, use notation from the respective references.
First, we consider the Hessian averaging method studied by (Na el al., 2022). It is important to distinguish between the condition they impose on the stochastic Hessian oracle H, and the Hessian approximation guarantee that they obtain for the actual estimate H, that they use (and that we use in SVRN-HA). The stochastic oracle is only required to have a sub-exponential tail (see their Assumption 2.1, and also Example 2.3 illustrating this for Subsampled Newton). However, the actual estimate H, is a result of averaging many samples from that oracle. Their local convergence analysis is only deployed once enough oracle samples are averaged so that H, achieves the approximation guarantee given in their Lemma 3.5, i.e. (1 (1- HH H, < (1 + H This approximation guarantee is strictly stronger than ours, but it becomes equivalent once 2.
Next, we consider (Roosta-Khorasaa & Mahoney, 2019) which analyzes a broad class of Subsampled Newton methods. In this paper, the most relevant results are Lemma 2 (Hessian approximation guarantee) and Theorem 5 (local convergence result). The guarantee reduces to our condition with a S 2, with the only difference being that. their Hessian approximation is restricted to the "cone of feasible directions", defined in (3). This restriction is only present in a constrained optimization setting (we focus on unconstrained optimization). The lemma also shows that the required Hessian sample size is again larger than the condition number of the problem (their condition number K1 matches our for unconstrained optimization).
Next, we look at (Bollapragada et al., 2018) which presents convergence analysis of Subsampled Newton under slightly different assumptions. Here, the key statements for local convergence are Lemma 2.4 in the journal version (Lemma 2.3 in arxiv), and equation (2.17). The lemma gives an approximation guarantee for the subsampled Hessian. This guarantee is in some sense weaker than our condition, because it only requires the Hessian approximation to be good in one direction, i.e., Ws - w² in their notation. However, examining the convergence bound in (2.17), for the local convergence analysis to hold, the Hessian sample size must still satisfy [St] 2 /²/p² where o² is effectively the upper bound on the component Hessians (potentially large as our X-smoothness) and ji is the lower bound on the component Hessians. The latter is effectively a component-wise strong convexity constant, which can be much smaller (i.e., worse) than our global strong convexity p. In summary, their Hessian approximation condition for local convergence analysis, while slightly different, also requires the Hessian sample size to be larger than a condition number of the problem. Their condition number can be much larger than our condition number x, or even infinite (for problems as simple least squares), and is less standard in the literature.
24

---

Published in Transactions on Machine Learning Research (01/2025)
Finally, examine (Erdogdu & Montanari, 2015), where the authors consider Subsampled Newton with possibly low-rank approximation of the Hessian. The most relevant result. in that work is Lemma. 3.1. Here, the standard version of Subsampled Newton is recovered when we let Q' = 55 =1 Then, the standard Hessian approximation condition appears implicitly through the fact that č has to be less than for the bound to be non-vacuous. To see the condition more clearly, we point to Equation (B.1) the appendix (Equation A.1 the arxiv version), which requires that |["'|l |Hs, - H|| < 1. For Q' Hg, this is essentially equivalent to our condition with a S 2. Also, from the bound in Lemma 3.1, we see that once again the condition requires Hessian sample size to be larger than condition number of the problem (which is for them K [I'|l, and after some effort, this can be seen comparable to our condition number). In the main algorithm of the paper, NewSamp, the authors aim to reduce the required sample size by using a different Q' computed from a low-rank approximation of Hs. This roughly corresponds to constructing a Hessian a-approximation with < a <
Omitted proofs
Here, we include the proofs of the auxiliary results stated the paper. First, we discuss in detail the global convergence analysis of SVRN-HA (Theorem 4). Then, we illustrate how the Hessian approximation required in Theorem 3 can be obtained via subsampling.
D.1 Global convergence of SVRN-HA
Here, we show how the proof of Theorem 4, i.e., global convergence of SVRN-HA (Algorithm 1), follows from global convergence analysis of Hessian averaging (Na al., 2022). They show Lemma 3.5 that we were to run the global phase of SVRN-HA exclusively, then for any e, E (0,1) there is T := T(e, such that with probability 1- 8 for all T we have & (c),), H, = v2f(x,), and n7 = 1. This means that, for any e, the probability that the above event does not happen with any T < co is less than any 8 > O, so it must be 0. This implies that SVRN-HA will eventually switch to the local phase (i.e., to SVRN). Note that it is possible that the switch will occur before the local neighborhood and Hessian approximation conditions are met. But if this causes SVRN to produce a poor descent direction, will be caught by the line search (resulting in 7/ < 1) and the method will simply revert back to the global phase. Eventually, the global phase will ensure that both conditions met, and we can rely on Theorem for the local convergence analysis.
D.2 Hessian approximation via subsampling
Here, we illustrate how the Hessian a-approximation condition (2), used in Theorems 1 and 3, can be obtained via uniformly subsampling O d) component Hessians. This result follows from Bernstein's concentration inequality for random matrices, given below (Tropp, 2012).
Lemma 7 (Matrix Bernstein's inequality) Let Zi,..., Zy be independent random symmetric d x d ma- trices such that >, E[Z,] Z. Suppose that:
1 Σ k
(iz E|Z)) 5² and 1Z - -E|Za] S R.
Then, for any e 0
k
kk/2 72 + eR/3
Pr
Σ Z. -
2 2e < 2d exp (-
We are now ready to show the approximation guarantee for subsampled Hessian estimate.
Lemma 8 Suppose Assumption 1, and let D be the sampling distribution for component functions s, as in Theorem 3. Let tha,.., th ~ D be i...d. samples from this distribution. There is absolute consland such
25

---

Published in Transactions on Machine Learning Research (01/2025)
that for any Rd, with probability 1 8, the matric
H = (1+ 1/2 V (x))
-mmax with 7 max { 122 (22/1og(2//0)/k
is an a-approzimation of 4 f(x) as in (2) with
a= 1 + 0 (log(d/8)/k vklog(d/o)/k)
Proof Let H., = fxf(x) + 71 for some 7 2 0. We will use Lemma 7 with Z: = H,, ( 1/2 First, note that [[)) f() so that Z Z H, 22 1/2 = Next, using that | 1/(n + 7), we compute the boundedness parameter R:
21 H (x +
IZ. E|Z||
R.
#+ ?
Now, we similarly bound the variance parameter σ²:
[[ - E|Z,) < | < |E||ZIZ]]|| <
|B| < + =: 0²
Thus, Lemma 7 implies that z
3A μ+ク log(2d/0)/2 then with probability 5 the Hessian estimate
H (1+ ) (1 (x) + +I satisfies:
H 1/2 (1+ 1) 1/2
HH,
1/2 1/2 = k i=1
Σ Z.
2
sc
sc
We can rewrite this as:
(1-6)I H 1/2 H2 1/21111 ( 1/2
HH 1/2 3 (1+
1
which is equivalent to:
- (1- 1/2 1
H =1 3)
-1/2
H,,
1
Moreover, note that the regularized Hessian H. satisfies:
fxx v2 Y 1I = (1+ 2 42 f(x).
-1=
/
Putting this together, and assuming that 5 1/2, for k
32 aty
log log(2d/))/e 0) we get:
1 v((())
ú H vav2 (x), with a= (1+ 2e) (1 2).
Va
Thus, if we set 7 = max 12A (2d/0)/kk mb, then there are two cases:
If < 12k log then we have = 12 log(2d/6) 2 3A n+ク log(2d/8)/4 for e = 1/2, which implies
that H an a-approximation of fx) with a Ox log(d/ (d//)///k))) k)/k).
2. If l 2 12xd)g), then let = vn/2u, so that we have = 12 log(2d/6) = 3 lo (2d/6)/ee
which implies that H is an a-approximation of (² with a = (1 2c)s = O(1 + //xlo(d//)//
26


Reviews:
[{"review_id": "IWOPpt0JiY", "reviewer": "Reviewer_8WgW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* The exposition is relatively clear; the paper is sufficiently well written.\n* The proposed method is conceptually simple to understand and a simple extension combining SVRG and Subsampled Newton.\n* The assumptions are relatively standard for second order methods.\n* The theoretical results demonstrate improved local convergence compared to vanilla subsampled newton.\n* Numerical results provide convincing evidence that the proposed SVRN improves local convergence of Subsampled Newton, while also notably reducing wall clock time relative to full-batch Newton.\n* Some limitations of the method are discussed; ablations on leverage scores in linear least squares highlight the sensitivity of the proposed approach to subsampling for building sufficiently accurate Hessian estimates, where the number of sampled components needed to produce sufficiently accurate estimates will depend on the component-wise smoothness of the objective.\n\nWeaknesses:\n* The motivation for the proposed approach and comparison to related work could be better articulated. Why would inexactly preconditioned SVRG (e.g., Liu et al., 2019) not yield a similar effect as the proposed SVRN (improving local convergence with large mini-batches). Some positives; the comparison to full-batch second order methods is clear (SVRN using stochastic estimates). Moreover, the comparison to subsampled second order methods is clear (they don't have variance reduction). The comparison to certain accelerated variance-reduction methods is also somewhat clear; some of these approaches rely on momentum (such as Katyusha).\n* While the proposed local convergence theory has removed the dependence on the condition number relative to accelerated variance-reduced methods, the minimum mini-batch size required for local accelerated convergence still depends on the condition number of the objective.\n* Novelty is somewhat limited, initial phase is just subsampled newton, latter phase is SVRG using the hessian estimates from subsampled newton; however, as novelty is not a TMLR criterion, this will not affect by judgement of this work.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work proposes SVRN (Stochastic Variance-Reduced Newton) for finite-sum minimization of smooth strongly convex functions. \n\nAssumptions:\n* Specifically, it is assumed that the global objective is strongly convex and the component objectives are each L-smooth and convex.\n* It is further assumed that the objective is twice continuously differentiable and is either self-concordant, or has a Lipschitz Hessian.\n\nMethod:\n* The proposed approach is to first commence with a variant of Subsampled Newton (Na et al., 2022), where the full batch gradient must be computed in each iteration, and the Hessian estimates are produced by random sampling of individual components, and Hessian Averaging across time to improve the estimates (i.e.; directly leveraging Subsampled Newton with Hessian Averaging; SN-HA).\n* Once SN-HA has converged to a sufficiently small neighborhood of the minimizer (as determined by the SN-HA Armijo line-search), the proposed SVRN method kicks in.\n  * The proposed SVRN method, which only operates in the local neighboorhood of the minimizer, is essentially SVRG, scaled by the inverse Hessian estimate, which is computed using random sampling and Hessian Averaging.\n\nResults:\n* It is proven that SN-HA proposed in previous work will indeed converge to a sufficient local neighborhood of the minimizer, such that SVRN will be able to kick in. \n* It is also proven that, given sufficiently accurate Hessian estimates, SVRN will achieve accelerated local convergence to the global minimizer. Specifically, the parallel time (number of parallel batch queries) for SVRN is shown to be similar to Stochastic Newton, while Sequential time (number of passes over the dataset) is reduced by a factor of ~log(n).\n* Numerical results are provided on least squares and logistic regression tasks, comparing both the iteration-wise and wall-clock convergence of a) full-batch Newton, b) Subsampled Newton with Hessian Averaging, c) SVRG, and d) SVRN (where SN-HA is used to enter the local neighborhood). Ablations are conducting in the least squares setting examining the effect of leverage scores on the proposed approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "* Please add comparisons to accelerated variance-reduced methods in the numerical results, both momentum-based, such as Katyusha, and with inexact preconditioning.\n* It is not clear to me why SN-HA and SVRN-HA have non-identical convergence early on in training (Figure 2).\n* Please consider annotating the points in the plots at which the SVRN update kicks in, and the SN-HA updates halt.\n* Please elaborate on how the problem conditioning still impacts the derived theory; ideally second order methods should be robust to problem conditioning, which is reflected in your rates, but somehow still appears to pop out from the bounds on the required hyper-parameters.\n* While a main motivation of this work is explaining how SVRN improves local convergence with large batch sizes, it is not clear to me in the exposition or the numerical results where the improved dependence on batch size comes in relative to momentum accelerated variance-reduction methods.", "broader_impact_concerns": "N.A.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "iBuwhkLi28", "reviewer": "Reviewer_X5u3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. The introduction of SVRN is a significant contribution, as it combines the benefits of variance reduction with second-order information, leading to faster convergence rates for finite-sum minimization tasks.\n2. The paper provides a thorough theoretical analysis of SVRN, proving that it accelerates the convergence rate of existing stochastic Newton methods. \n3. The proposed initialization method using a few iterations of a Subsampled Newton method with line search is practical and effective. This initialization significantly accelerates the convergence of Subsampled Newton.\n\nWeaknesses:\n1. While the theoretical analysis requires fresh samples of components for each small step, the practical implementation suggests that resampling per stage can be more efficient. The authors should provide a more detailed discussion on the trade-offs between theoretical guarantees and practical efficiency in terms of communication cost.\n2. The performance of SVRN is sensitive to the Hessian sample size, especially in the initial stages. The authors should provide more insights into how to choose an appropriate Hessian sample size for different types of problems and datasets.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper \"Stochastic Variance-Reduced Newton: Accelerating Finite-Sum Minimization with Large Batches\" introduces Stochastic Variance-Reduced Newton (SVRN), a novel algorithm designed to accelerate finite-sum minimization tasks, particularly in the context of large datasets. The main contributions and new knowledge presented by the submission can be summarized as follows:\n1. The authors propose SVRN, a finite-sum minimization algorithm that leverages variance reduction techniques to accelerate existing stochastic Newton methods. This is achieved by reducing the number of passes over the data from $O(\\alpha log(1/\\epsilon))$ to $O(\\frac{log(1/\\epsilon)}{log(n)})$, where $n$ is the number of sum components and $\\alpha$ is the approximation factor in the Hessian estimate. This improvement is significant, especially for large datasets.\n2. The paper provides a rigorous theoretical analysis of SVRN, demonstrating that it achieves a faster convergence rate compared to existing stochastic Newton methods. Specifically, the sequential complexity is improved by a factor of $O(\\alpha log(n))$, while maintaining the same parallel complexity.\n3. The authors propose a practical initialization method for SVRN by running a few iterations of a Subsampled Newton method with line search. This approach is shown to substantially accelerate Subsampled Newton.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. About the communication cost: It would be better to provide a detailed discussion on the trade-offs between theoretical guarantees and practical efficiency, particularly regarding the frequency of gradient resampling, and explain the impact of resampling per stage versus per step.\n2. About Hessian sample size: I think conducting a sensitivity analysis of the Hessian sample size, especially in the initial stages of the algorithm is important. And it would be better to provide guidelines on how to choose an appropriate Hessian sample size for different types of problems and datasets.", "broader_impact_concerns": "This paper primarily focuses on theoretical research and does not present any specific ethical concerns that need to be highlighted.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "FNVxJD58D6", "reviewer": "Reviewer_DCx8", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper seems strong, and like a good fit for TMLR.\n\nIt proposes a simple algorithm with a clean motivation -- merging variance reduction with stochastic newton methods. My personal expertise is not specifically in netwon methods, so I cannot personally attest to the significance of this result relative to the prior work. That is, I am not specifically aware if there is or is not any other variance-reduced stochastic newton method. That said, taking the author's word for it, I feel that this core pitch is compelling.\n\nThe theorems are well written, and I have no reason to doubt their accuracy (though I did not verify the proofs in detail). The presentation is all around very solid. The experiments are interesting and carry a nice strong message. It's overall a good paper. I do have some nitpicks and points where I lack clarity, but overall I'm happy to recommend publishing this paper at TMLR.\n\nNotably, both of my nitpicks can be resolved with cleaner citations to the prior work. In general, that's really my biggest issue with the paper. As someone who understands all of the words on the paper but isn't an expert in finite-sum minimization or stochastic newton methods, I don't know exactly where to look in the prior work to either A) fully appreciate the novelty of this result or B) verify that the assumptions made in their paper are standard or C) verify that their assertions about the prior work are correct. It's easily solved by adding more clear and direct citations.\n\n---\n### Nitpicks\n\nI have two key technical points where I lack some clarity.\n\nThe first is about the metrics used in the paper. On the bottom of page 2, the authors discuss the two complexity metrics used in the paper: parallel complexity and sequential complexity. Parallel complexity is completely fine by me and makes a lot of sense (it's the number of \"rounds of adaptivity\" or \"batches of gradient computations\"). Sequential complexity feels a bit odd in the way it's phrased. They describe it as \"passes over the data\", but in a kinda weird way imo:\n> One pass corresponds to $n$ queries of component gradients.\n\nWhich strikes me as kinda odd because it seems that n queries to the first gradient $\\nabla \\phi_1(x)$ would count as \"1 pass over the data\". That is, the sequential complexity seems to just be $\\frac{\\text{Total number of gradient computations}}{n}$, which doesn't sound a lot like \"passes over the data\" to me. This is why I phrased my \"Summary of Contributions\" in terms of the total number of gradient computations instead of the sequential complexity. A few sentences clearing up why it's okay to think of $n$ queries to $\\nabla \\phi_1(x)$ as a single pass over the data would be nice. Or perhaps just a citation to a paper or two where that particular definition of sequential complexity is used.\n\n\nThe second point where I lack clarity is in the $\\alpha$-approximation to the Hessian. On page 3 in equation (2), the authors say that their analysis is built on having access to an approximation to the Hessian $\\tilde H$ such that\n$$\n\\frac1{\\sqrt\\alpha} \\nabla^2 f(x) \\lesssim \\tilde H \\lesssim \\sqrt{\\alpha} \\nabla^2 f(x).\n$$\nThis notion of an approximation slightly surprised me (namely, why is there a square root on the $alpha$ -- spectral approximations tend to just have a constant there without the square root). But the authors then show in Table 1 on page 3 that the prior work on Stochastic Newton methods does use this notion of $\\alpha$-approximation (see the third row of the table, titled \"Stochastic Newton\"). I wanted to verify that the prior work uses this exact notion of $\\alpha$-approximation, but I had a hard time finding this exact expression. I checked in [Roostra Mahoney], [Berhas et al], and [Derezinkski et al]. I would love to see either A) a clear citation to a paper which uses this exact notion of $\\alpha$-approximation or B) a discussion of why this notion of $\\alpha$-approximation is right.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes a new algorithm, called Stochastic Variance-Reduced Newton (SVRN), for finite-sum minimization problems. It assumes we are minimizing the average of $n$ functions, and that this overall loss function is strongly convex and smooth (with condition number $\\kappa$). Further, we assume that we can compute each of the $n$ functions, their gradients, and their Hessians.\n\nThe core claim of the paper is that, if we can access an $\\alpha$-approximation to the Hessian (i.e. can access a matrix $\\tilde H$ such that $\\frac1{\\sqrt\\alpha} \\ \\nabla^2 f(\\vec x) \\lesssim \\tilde H \\lesssim \\sqrt{\\alpha} \\ \\nabla^2 f(\\vec x)$), then SVRN converges using just $O(\\alpha \\log(1/\\varepsilon))$ batches of gradient computations / rounds of adaptivity and $O(\\frac{n}{\\log(n)} \\log(1/\\varepsilon))$ total gradient computations. These bounds are local convergence bounds, assuming that we start iterating from a point $x$ not too far from the global comptimum $x^*$.\n\nThis rate matches the best known bound on rounds of adaptivity used in prior work, and outperforms the best known bound on total gradient computations by a factor of $\\log(n)$.\n\nThe core premise in the algorithms is to integrate variance reduction techniques with a stochastic newton method.\n\nA combination of both theoretical analysis and empirical evidence is provided.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The bullet points with a $\\star$ below are important to me. Most of these are fairly minor (add citations; clarify the meaning of a column in a table, etc). The bullet points without a $\\star$ are recommendations and can be dropped.\n\n- $\\star$ [Page 2, Parallel Complexity at the bottom] add a citation to the standard PRAM model\n- $\\star$ [Page 2, Sequential Complexity at the bottom] add a citation to where this is used (see \"nitpicks\" above)\n- $\\star$ [Page 3, Table 1] clarify the speedup column. Is this a speed-up in sequential complexity? The speedup of SVRN wrt mini-batch SVRG and mini-batch Katyusha is entirely in the parallel complexity. The speedup of SVRN wrt stochastic newton is only in sequential complexity. This is kinda confusing to me overall. I don't really get what speedup you're talking about.\n- $\\star$ [Page 3, Table 1] would be nice to add some citations to which particular papers your getting these particular rates from (see \"nitpicks\" above wrt $\\alpha$).\n- $\\star$ [Page 3, equation 2] Cite somewhere or otherwise argue that this notion of an $\\alpha$-approximation is standard.\n- [Page 3, Theorem 1] If $n > \\frac1\\varepsilon$, then the Sequential time becomes $<1$ passes over the data. Is there some lower order term that dominates the sequential complexity when $n > \\frac1\\varepsilon (i.e. moderate accuracy regime)? Or do you really perform less than one pass over the entire dataset?\n- $\\star$ [Page 4, below equation 3] Cite the $O(\\alpha \\log(1/\\varepsilon))$ rate.\n- $\\star$ [Page 4, \"In this setting\" sentence] Doesn't sampling $k = O(\\kappa \\log(d))$ component Hessians get you back to the rates of SVRG? What's the upside here? The following sentence mentions a result of [Erdogdu Montanari] which uses far fewer samples. Can we get a value for $\\alpha$ here which is much better than $\\kappa$? There any theoretical rates you can place down to show this?\n- $\\star$ [Page 9, Figure 2] Add units for the wall clock time (seconds? minutes? hours?). Add confidence intervals over the multiple runs of the algorithm.\n- $\\star$ [Page 9, \"All of the convergence plots\" sentence] How many times did you run each experiment? That number matters.\n- [Sections 5.1 and 5.2] Could you write down ballpark values of $\\kappa$ and $\\alpha$ for this problem? It's okay if it ain't possible, but it'd be nice to see what values we're really dealing with in relation to the preceding theory.\n- [Page 11, Figure 4] Would be nice to also see the wall clock times for these experiments. Using a smaller Hessian or no preconditioning should make the iterations faster. At least that's my intuition. So, I'd really like to see how those all balance out. Can easily be in the appendix though.\n- [Page 11, \"Effect of Hessian Accuracy\"] I'd really like to see the value of $\\alpha$ in a plot or paragraph here. At the minimum, I'd like to see ballpark values for $\\alpha$ attached to the trends in figure 4a. Preferably, I'd love to see a plot with something like $\\alpha$ on the x-axis and \"iterations needed for convergence\" on the y-axis. Or perhaps a plot with iterations on the x-axis and $\\alpha$ on the y-axis, with a different series for each value of $h$ used in Figure 4(a). In essence, I don't really feel like I've seen the real dependence on the accuracy of the Hessian. Instead, I've seen that different sample complexities can matter. These claims are not totally equivalent.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Authors already responded to the reviews thoroughly. There are a few minor issues that remain to be addressed.

- please include a discussion on the effect of condition number on the minimum mini-batch size (which was also raised by other reviewers).
- please clarify the numerical results in Figure 2 that the SVRN update kicks in after only a couple of iterations, which also explains why the SN-HA and SVRN updates are not identical in the early phases of the plot.
- please emphasize the motivation of the work in the introduction.
- please discuss differences between theoretical guarantees and practical efficiency (e.g. in terms of communication cost).

Once these are addressed, the paper will be ready for publication at TMLR.
Decision: Accept (poster)