Newton Meets Marchenko-Pastur: Massively Parallel Second-Order Optimization with Hessian Sketching and Debiasing
Authors: Elad Romanov
Fangzhao Zhang
Mert Pilanci
Venue: ICLR 2025 Poster (2025)
Keywords: parallel and distributed optimization
Newton method
Marchenko-Pastur law
Hessian sketching
Primary area: optimization

Abstract: Motivated by recent advances in serverless cloud computing, in particular the ``function as a service'' (FaaS) model, 
we consider the problem of minimizing a convex function in a massively parallel fashion, where communication between workers is limited.
Focusing on the case of a twice-differentiable objective subject to an L2 penalty, we propose a scheme where the central node (server) effectively runs a Newton method, 
offloading its high per-iteration cost---stemming from the need to invert the Hessian---to the workers. 
In our solution, workers produce independently coarse but low-bias estimates of the inverse Hessian, using an adaptive sketching scheme. The server then averages the descent directions produced by the workers, yielding a good approximation for the exact Newton step. The main component of our adaptive sketching scheme is a low-complexity procedure for selecting the sketching dimension, an issue that was left largely unaddressed in the existing literature on Hessian sketching for distributed optimization. Our solution is based on ideas from asymptotic random matrix theory, specifically the Marchenko-Pastur law. For Gaussian sketching matrices, we derive non asymptotic guarantees for our algorithm which do not depend on the condition number of the Hessian nor a priori require the sketching dimension to be proportional to the dimension, as is often the case in asymptotic random matrix theory. Lastly, when the objective is self-concordant, we provide convergence guarantees for the approximate Newton's method with noisy Hessians, which may be of independent interest beyond the setting considered in this paper.

Published as a conference paper at ICLR 2025
NEWTON MEETS MARCHENKO-PASTUR. MASSIVELY PARALLEL SECOND-ORDER OPTIMIZATION WITH HESSIAN SKETCHING AND DEBIASING
Elad Romanov, Fangzhao Zhang, and Mert Pilanci
Stanford University
fercmanov, zfzhao, pilanci (alaanii@ttanford..ell
ABSTRACT
Motivated by recent advances in serveriess cloud computing, in particular the "function as a service" (FaaS) model, we consider the problem of minimizing convex function in massively parallel fashion, where communication between workers is limited. Focusing on the case a twice-differentiable objective subject to an L2 penalty, runs a Newton method, offloading its high per-iteration post-stemming from the need to invert the Hessian- to the workers. In solution, workers produce inde- pendently coarse but low-bias estimates of the inverse Hessian, using an adaptive sketching scheme. The server then averages the descent directions produced by the workers, yielding a good approximation for the exact Newton step. The main component of our adaptive sketching scheme is low-complexity procedure for selecting the sketching dimension, an issue that was left largely unaddressed in the existing literature on Hessian sketching for distributed optimization. Our so- lution is based on ideas from asymptotic random matrix theory, specifically the Marchenko-Pastiee law. For Gaussian sketching matrices, we derive non asymp- totic guarantees for our algorithm which do not depend the condition number the Hessian nor priori require the sketching dimension to be proportional the dimension, as is often the case in asymptotic random matrix theory. Lastly, when the objective is coff-ooncor,an,, we provide convergence guarantees for the approximate Newton's 's method with noisy Hessians, which may be of independent interest beyond the setting considered in this paper.
propose scheme where the central node (server) effectively
INTRODUCTION
Consider minimizing a convex, twice-diiferentiable function F Rdd - R, subject to an L2 regular- ization penalty:
min F(0) 8E34
11013
where Such problems frequently appear in machine learning and statistics. Common ex- amples include ridge regression, logistic regression with regularization, support vector machines (SVMs) and kernel machines, among others (Hastie et al., 2009).
Consider a scenario where a large number of workers (processors) are available to collaboratively solve (1), with the assumptions that each worker (i) has full access to the objective function; (ii) op- erates with limited individual computational resources; and (iii) cannot communicate directly with other workers, only with central server, with the latter responsible for orchestrating the work of the workers (a "star" network topology). This setting is motivated by recent developments in serverless cloud computing, particularly the "function as service" (FaaS) model (Jonas et al., 2017). Server- less computing and FaaS are well-suited for computing gradients and Hessians on large datasets, enabling parallelization of these computations across large number workers. Additionally, this approach provides resilience failed or straggler workers, ensuring robust and efficient processing (Bartan & Pilanci, 2019). In designing a solution of this setting, we are essentially guided by two key principles: (a) minimizing the number of communication worker deployment rounds, with

---

Published as a conference paper at ICLR 2025
each round utilizing large number of workers in parallel; and (b) reducing the computational load on individual workers by distributing smaller tasks across many workers concurrently.
This paper proposes and analyzes scheme for solving (1) in massively parallel manner. In our scheme, the server aims to solve (1) using Newton's ' method, a second-order iterative algorithm which is known converge extremely fast for sufficiently smooth and strongly convex objectives (quadratic convergence rate). Its fast convergence is attained by incorporating curvature (Hessian) information when searching for a descent direction, and specifically requires an inverse Hessian- gradient product at every iteration. This a priori results in massive per-iteration computational cost: exactly computing the descent direction (assuming the Hessian is accessible) costs practically (dd) flops; in high- sigh--dmenssonal problems (large d), this may be prohibitively expensive.
We wish to offset the high per-iteration cost of Newton's 's method, by deploying many workers in parallel to collaboratively compute the Newton direction. In the setting we consider, workers cannot easily communicate with one another, and therefore directly inverting the Hessian by is not feasible. In our scheme, each worker independently computes coarse but unbiased estimate of the exact Newton direction, using randomized sketching. These estimates are then aggregated (averaged) at the server, producing a good approximation of the true Newton step.
Newton's method:
The exact Newton method for minimizing (1) has the following form:
0, = 0,-1 ar Wi9t.
(t = I')..
(2)
where is a step size (typically chosen by line search) and the gradient and inverse Hessian are
9t FF((-11) + A0,-1,
(3)
We = (H, A H. V²F(0--))
(4)
For sufficiently smooth convex objectives F. Newton's method is known converge at a quadratic rate: T = Olog log(1 /e)) iterations are sufficient to approach the minimum within error e, as opposed T (1lgg(1 for gradient descent and its accelerated variants; see (Boyd & Van- denberghe, 2004, Chapter 9). A fast convergence rate is highly desirable in our setting, as the number of iterations corresponds directly to cummunication/worke. deployment rounds.
In our scheme, the server effectively runs the Newton method with an approximate inverse Hessian:
A. 0, a, W.9t:
(5)
where Win obtained by averaging q independent local estimates produced by the workers:!
9
W
(6)
Debiased inverse Hessian sketching: We wish to avoid directly inverting the full d x d Hessian, which practically (when done exactly) costs (d To this end, we have every worker compute, in parallel, a cheap but low-bias estimate of the inverse Hessian by sketching. Specifically, at every deployment round, every worker independently samples random sketching matrix S E m (where m. < d), and approximates the inverse Hessian by its sketched version
W = S (SHST + ÅT)-'S.
where > is a modified regularization parameter. Sketches of the type in (7) have been stud- ied before in the literature (also in the context of high-dimensional regression and optimization). Variations have appeared under various names: "dual sketching" (Zhang et al., 2013), "sketch-and- project" (Gower & Richtárik, 2015), "right-sketch" (Murray et al., 2023), and "feature sketching" (Patil & LeJeune, 2023). Note that the cost of forming Wg is O(m² + .), where M is the cost of multiplying SHST In this paper, we exclusively consider the case where S is a random dense matrix with i.i.d. entries.? In this case M d'm, which dominates the overall per-iteration cost. Naturally, we are particularly interested in settings where m can be chosen very small compared to d, so that cost, m))m)) substantially smaller than (dd))
In practice, the workers send the vector (+) St E R to the server, rather the whole matrix (i)(k) e dxdd BB "dense", we mean that S has M(md) many non-zero entries typically, though it's not necessary that all the entries, or majority of them, are non-zero.
2

---

Published as a conference paper at ICLR 2025
Remark 1.1. In this paper we do not assume any particular form for H. (i.e., that it is readily fac- torizable H = X x, X E RN a setting which is quite common in the sketching and distributed optimization literature; for example (Pilanci & Wainwright, 2015)). In fact, sketched measurements SHS may be obtainable even in settings when F is not available in analytic form, but only through a program (oracle) that calculates it for example by automatic differentiation (Paszke et al., 2017).
key question is how to choose m (and in accordance λ). Define the X- effective dimension of H,
(A == (r(H(( >-)
(8)
As we shall see in Section 2, results from asymptotic random matrix theory imply that, when S is an i.i.d. dense sketching matrix and provided that m > dæ (A), the choice
1 X==(1--dd(() da) A= -1- m
results in W which is low-bias estimator of the true inverse Hessian W. Note that exactly comput- ing the effective dimension is resource-intensive and essentially requires inverting the Hessian the very operation we wanted to avoid. A key component of our proposed scheme is novel low-complexity procedure, based on ideas from asymptotic random matrix theory, for adaptively selecting m = O(dH (A)) and X, by only observing sketched mcasurements of H.
Remark 1.2. When H is effectively low-rank, that is, it exhibits fast spectral decay, dh (A) may be much smaller than d. For example, when the spectrum has power decay X (H) x k and == O11,, dh (A) aa) when a < 1, ()( = Ollog d) when a = 1 and (A) = O(1)
when > 1. Hessian spectral with power decay common in regression problems involving kernels ( (Wainwright, 2019; Yang et al., 2017).
Outline and Contributions:
In Section 2 we review fundamental results from asymptotic random matrix theory (RMT), specif- ically the Marchenko-Pasee Pastur law and its formulation using deterministic equivalents, and demon- strate how they lead naturally to an asymptotically unbiased inverse Hessian estimator, from i.i.d. sketching matrices S. The connection to the Marchenko- Pastur law treats in a uniform framework-and and recovers almost immediately- earlier results on Hessian sketch debiasing for distributed Newton methods (Derezinski et al., 2020; Zhang & Pilanci, 2023).
In Section we propose a novel data-driven method for selecting m and À using only the sketched Hessian, avoiding the need to compute the effective dimension exactly. The approach is natural outgrowth of the RMT framework described previously, and is fully adaptive: no tuning param- eters needed. Section 3.1 further provides non-asymptotic error guarantees for our method, focusing on the case where is an i.i.d. Gaussian matrix. Notably, the error bounds obtained in this setting do not depend the Hessian's condition number, require m to be proportional to d in any particular way (in RMT, such bounds are sometimes called "dimension-free In Section 4 we derive convergence guarantees for the Newton method with inexact Hessians; to wit, when the exact inverse Hessian W in (2)-(4) is replaced by an naacurraee estimate W. By ensuring a sufficiently small 7). the algorithm achieves an arbitrarily fast linear convergence rate. Our results are proved for self-concordant functions F(-), and to our knowledge are novel and may be of independent interest. Combined with the results of the previous section, we obtain a non-asymptotic, end-to-end convergence guarantee for the entire parallel method (with Gaussian sketches), which remarkably does not depend on the condition number of the Hessians. Lastly, Section 5 is dedicated to experiments, on both synthetic and real-world data.
1.1 RELATED WORKS
Distributed second-ordes optimization methods have been widely studied for large-scale machine learning. Notable examples include DANE (Shamir et al., 2013), which addresses communication efficiency, AIDE (Reddi et al., 2016), designed for accelerated convergence, and DiSCO (Zhang & Lin, 2015), which focuses on distributed second-order methods for convex optimization problems.
d (A) is also called the "effective degrees of freedom" in ridge regression (Hastie et al., 2009), when the design matrix X Nxd H x X.
3

---

Published as a conference paper at ICLR 2025
The most closely related body of literature to our work centers around the averaging of inexact Newton steps, which has proven effective when combined with randomized linear algebra. GIANT (Wang et al., 2017) uses the average of local Newton steps the global step and Determinan- tal Averaging (Derezinkii et al., 2020) provides an unbiased averaging technique for distributed Newton methods. Surrogate sketching (Derezinski al., 2020) was introduced for distributed New- ton's method, revealing improvements via a simple shrinkage strategy. The work (Zhang & Pilanci, 2023) introduced an optimal debiasing method for distributed Newton's method when the effective dimension of the Hessian is known. We remark that the setting being addressed in distributed optimization is typically different than ours: in the former, the objective (which is often the empir- ical loss over data set, i.e. in federated learning (Kairouz et al., 2021)) is spread over different machines, with the goal of minimizing the overall transfer of data between machines. In constrast, in our setting all machines have access to the full data, and the goal is to leverage parallelism to reduce the overall runtime and number of deployment rounds. The literature on distributed opti- mization, and distributed second-order methods in particular, is vast and growing; additional recent references include (Safaryan et al., 2022; Agafonov et al., 2022; Qian et al., 2022: Elgabli et al., 2022: Bischoff et al., 2021), among others. Lastly, there exists a large related corpus of works in optimization involving Hessian subsampling or sketching, for cases where the Hessian admits good low-rank approximation, though not necessarily in a distributed setting. Examples include (Roosta-Khore) Khorasani & Mahoney, 2019; Gower et al., 2019; Frangella et al., 2022), among others. The basic idea of approximating the Hessian in second-ordes optimization gquasi-Newtor methods) is essentially classical, cf. the book (Nocedal & Wright, 1999). Additional recent references explor- ing related themes include (Gupta et al., 2018; Vyas et al., 2024; Jahani & Rusakov, 2022; Doikov et al., 2023; Na et al., 2022), among many others.
2 DEBIASING FROM L.L.D SKETCHES AND RANDOM MATRIX THEORY
This section concerns inverse sketched Hessian debiasing for sketches S E Rg which have i.i.d. entries with mean zero EjS,,] = O. We normalize the variance as E[S,,] = 1/m; this normalization is such that E[ST S] = I, equivalently, S: Rdd + R is an isometry in expectation.
Let H > 0 be positive definite matrix (e.g., the full Hessian at the current iteration). The spectral properties of the matrix SHST E R=xx which appears in the sketch (7), are of utmost impor- tance for what follows. Note that this matrix has the same non-zero eigenvalues* as the matrix H SH2/2; the latter is a sample covariance matrix (making the analogy explicit: number of samples, H is the population covariance, and 2; == vmll are independent "sam- ples",r1..... ri.... I'm being the rows of S): accordingly, S SHST is sometimes called the companion sam- ple covariance. Sample covariance matrices, and their spectral properties, have been central object of investigation in random matrix theory (RMT); ef. (Bai & Silverstein, 2010). A foundational re- sult is the Marchenko-Pastie Pastur law (Marčenko & Pastur, 1967), which describes the high-dimensional global behavior of the eigenvalues of HT 1/2 (limiting spectral distribution).
is the
The Marchenko-Pastur equation: For 2 E C+ let ((2) E C+ be the unique solution of
1 - 1 (H(L+s()))))))) m
1
(10)
s(2)
subject to the constraint s(x) E C+. It is known that (see (Bai & Silverstein, 2010) s(z) is the Stieltjes transform" of compactly supported probability measure, whose support lies in Rzo: this measure is the (companion) Marchenko-Passe PPsttur law associated with the population covariance H (and depends only on its cigenvalues) and the number of samples
Moreover: 1) The function 2 + s(2) extends continuously to RR, except possibly at 2 = 0; 2) for x E IR outside the support, the function x + s(x) is increasing; in particular, for negative æ < 0, it is positive and inereasing. In this paper, we shall restrict our attention to only negative real valued = in (10), understanding that it
*This is a standard consequence of the elementary determinant identity det + AB) = det(I + BA). 5 'erre C+ is the complex upper halfplane, (2) > O.
For finite measure on R. its Stieltjes transform is the function s, C+ + C+ given by s,(2) du(xx) Using the Stieltjes inversion formula, the measure μ can be recovered from s. in the following
sense: for every sp(-) bounded and continuous, J (t))d((tt) lim 4 (s (2 + in)
4

---

Published as a conference paper at ICLR 2025
holds also for 0 real. For more details, cf. (Silverstein & Choi, 1995; Bai & Silverstein, 2010).
Finally, note that the function 2 + s(2) is invertible, the inverse function given by
1 m
v (s)
tr(H(I +
(11)
The Marchenko-Pasee equation (10) may equivalently be written as 2 = (s(2)).
In RMT, it is typical to consider the proportional growth asymptotic regime: m, d 0o, m (d).
Definition 1. For (small) constant 0, assumption set [Asymp(T, H, holds if: (1) bounded
aspect ratio (proportional growth): - < m < 1/TT (2) bounded operator norm: || < L/T:: (3) H is invertible; (4) H has non-degenerate spectrum: at least rd eigenvalues are larger than T.(5)
VmS., ~ E, with their law &, satisfying: E[g] = 0, Bi5 = 1, and E[[{]*] < co for all l > 0.
Denote by 5(2) the Stieltjes transform of the empirical eigenvalue distribution of SHST:
s(z) tr(SHS - -11
(12)
((SHST) -
Note that 6(z) is random function, since S is random. The following is (quantitative) statement of the Marchenko- Pastur law (Knowles & Yin, 2017); for simplicity, we focus on the case = < O. Theorem 2.1. Set n, D > Assume [Asymp(T, E).. W.p. (D), simultaneously for all
-1/T S 2 S -T,
[5(2) - s(2) = 0 (m (m--//1 11/21/+
(13)
The constants in the O() notation depend on T, 71 D.
Proof. See Theorem 3.7 and Remark 3.10 in (Knowles Yin, 2017).
Deterministic equivalent for the sample covariance: The Marchenko-Pasee Pastur law (Theorem 2.1) describes the eigenvalue distribution of SHS equivalently H SH and relates it to that of H (the connection made via the Marchenko- equation (10)). One can show something stronger: the entire (random) resolvent S SH1/22 =/)-1 is in fact close (entrywise, not in operator norm) to the resolvent of H (which is deterministic), with appropriate shifting and shrink- age. A result of the following kind is sometimes referred in the RMT literature as a "deterministic equivalent", cf. (Couillet & Debbah, 2011; ; Bai & Silverstein, 2010).
Theorem 2.2. Set 7) D > O. Assume [Asymp(T, H, E).. For any u. E R4 |ull. ||l = 1. the following holds. W.p. O(m -D), simultaneously for all 11/7 S 2
u 5 SH -1)- = u (- =1) '1 u + 0 (m-2//)
(14)
Above, the constants in the O(-) notation depend on T, D and are uniform in u, u. In particular,
( SH - =1)- (- =1) = O(m 11/2hn)
(15)
Proof. See Eq. (3.3) and Theorem 3.7 in (Knowles & 2017).
Deterministic equivalent for the sketched Hessian: Recently, deterministic equivalent for the sketched matrix (7) was obtained by (LeJeune et al., 2022). We cite a quantitative version of their result (LeJeune et al., 2022, Theorem 4.1).
Theorem 2.3. Set 77, D > 0. Assume [Asymp(T, H, and furthermore )rank(H) (H) 2 r. For any u, V E Idd, [|u]]. ||m||| = 1, w.p. 1 O(m ((mm)) -D), simultaneously for all 117 S 2 S-T. -T,
a T (SHST - -21) -I1--Se = u T (H (u -1/2411)
(16)
In particular,
E S (SHST- -1] (H (s( 0(m 011-117
(17)
5

---

Published as a conference paper at ICLR 2025
Proof. The result is consequence of Theorem 2.2; for completeness, let us repeat the argument of (LeJeune et al., 2022), assuming H is full rank (see their paper for a generalization), Denote = u, i == u,,,, so that || |i S $1/vT=0(1) 0(1). We have w.p. Omm-),
=
u RSHSS -
ū 4 (SHST - -x1)-
SH a= 2u (H SH - -- =1)-1 16 = =ūJo - i -
=a] =H! (
+ 1/2277) = uT (H + ( 0 (m--1/2+)) 1/2+t).
AN ADAPTIVE SKETCHING AND DEBIASING PROCEDURE
Theorem 2.3 suggests a clear path for obtaining low-bias estimate of the inverse Hessian from the
sketch (7), namely: one should choose À > O such that s--)) = 1/2, when such a root exists.
Lemma 3.1. A solution À >> O satisfying s(- 1/A exists if and only if m > dz(A). When this is the case, À is given by (9).
Proof. By merit of being the Stieltjes transform of a finite measure supported on Rg0. s(-) maps bijectively (-o0, 0) to (0, s(0),, hence solution À > O exists if and only if s0) > 1/A. Taking the limit 2 - 0+ in (10), we find (after routine algebraic manipulation) that s(0) solves dæ (1/s(0)
m. Since the function f + dr (p) is decreasing for positive μr, we deduce that 1/s(0) < X if and
only if m dæ (A). When this is the case, Ā = (1/A) where Vi is the explicit inverse (11).
O
Remark 3.2. Earlier works on distributed Newton's method (Derezinski et al., 2020; Zhang & Pilanci, 2023), considered debiasing for a closely related but different Hessian sketch. Specifically, they assume the Hessian has the form H = X X where X E RNxd, such Hessians naturally appear in machine learning problems, for example in the (reguarlized) empirical loss minimization for a generalized linear model (GLM). (Zhang & Pilanci, 2023) considered (among others), the sketch X TT SX, and calculated the appropriate shrinkage factor for asymptotically debiasing the inverse: (7:X7 S SX + x1)-1 (H + Al) -11 where 7 := 1/ (1 m > (A), under an asymptotic setting as in Theorem 2.2 (though without making the connection to the Strchenkk-Pustur law explicit). In the scheme they considered, the effective dimension was assumed to be known. We remark that Theorem 2.2 holds verbatim if H1/2 is replaced by XT, and the ratio T < d/N < 1/T is bounded (Knowles & Yin, 2017). Rewritting (15) as (H + (s(2))- ~ 20-1==(((222))))) (1-25(2))-1XXT S ST + (s(z))=11) -1 setting - recovers their proposed shrinkage factor. Accordingly, the adaptive procedure we propose below applies to that sketching scheme as well, with the obvious modifications.
m d (A)), provided that
Note that we do not have access to the function s(.) directly (to calculate it numerically via (10), one needs to invert the full Hessian)- but we can estimate it from the sketched measurements SHST, using Theorem 2.1. Accordingly, we can estimate À by the solution I of the (random) equation
5(- 1/2, λ>0. 0.
(18)
if it exists. To have â R 2, need that m > (A), and we do not know this a priori therefore we also need a way to select large enough m prior to solving (18). Moreover, to ensure small error in Theorems 2.1.2.3, one needs to work with bounded away from O. As Î = A(1 d (A)), to ensure that this is Q(A) we need to have m 2 (1 + ((1))d (A). Specifically, we shall aim to choose m between 1.5d (A) < m < 4 (A), the constants chosen somewhat arbitrarily.
Testing for a good m. We devise a test that, given m, with high probability: (a) rejects if m < 1.5d (A); and (b) accepts if m 2 (A). Set 20 -5/112. The tests accepts if 6(20) > 1/X, and rejects otherwise.
If mm < 1 (A) then either there is no solution Ā > to s(o)((-)) = 1/A, or there is one, specifically = X(1 - ب dh (A)) < A/3 < -00 Since z + s(z) is positive and increasing for z < 0, necessarily s(z0) < s( s(-)) 1/X. Provided that is large enough, w.h.p. š(z0) < 1/> as well,
6

---

Published as a conference paper at ICLR 2025
hence the test rejects. If m 2d (A) then À = X(1 (A)) X/2 -20. hence s(zo) s(-x) 1/x. Provided that m is large enough, w.h.p. $(20)
m > 1/2 as well, hence the test accepts.
mg large enough; proceeds iteratively, doubling m after each iteration m := 2m. At every iteration, apply the above test, halting if it accepts and continuing if it rejects. W.h.p., (a) we not halt on any m 1. (A); (b) we halt on the first m satisfying m 2 2d (A), and therefore return satisfying m I (A).
Searching for by doubling. Start with
Computational complexity. W.h.p., the doubling procedure stops within iterations. Each iteration costs most O(d-mz), my 2'm the total cost is Odd min (A), and dominated by the cost of multiplying SHST for S E R at the last iteration.
Having found large enough m, we can solve the equation (18) by binary search, and form the sketch (7) with Ả in place of A. For the reader's convenience, the overall adaptive sketching procedure- which is performed (in parallel) by every worker- is summarized in Algorithms 1-2.
Algorithm
procedure Choose Aihloee-Seecchiin---iieennnnnnns X, mo)
2: 1
mo 3: while m < d do
4: Sample i.i.d. sketch S E Rmmxd and form SHST
5: Compute eigenvalues (1HSS ... 2 Am (SHST).
6:
if > 1/A then
S(2) defined in (12).
7:
return m
> Guarantee: w.h.p., (.ddr(() S m & (A).
8:
else
9:
m = 2m
end if'
11: end while
return
13: end procedure
Algorithm
1: procedure Estimate- teaeee-eeeeeeeeessann--eee- -eeee----siiinDDe (H, X, m)
2:
Sample sketch S Rmd and form SHS Compute eigenvalues (SHST)) .... Am ( (S
3:
4:
if 5m. 1/A then
5:
ERROR; set Ā 52/12.
6: else
7:
Find 52/12 < À < À such that 5( 5(--) = 1/A.
> Can be done by binary search.
8:
end if
9: return (Á, W) where V = (SHST + +
D (In practice, return W g where
g Idd is the current gradient.)
10: end procedure
is straightforward to establish the asymptotic validity of our adaptive sketching procedure. Theorem 3.3. Set 17, particular mo (dd)), and furthermore that Arank(H) (H) 2 Let vi be the output of Algorithm 1, and (), w) be the output of Algorithm (the latter given input m = m). Then w.p. 1 - O(m P): max(mo. 1.5d (X)) < m S max(mo, 4 (X)}-
let mo be an initial sketch size. Assume [Asymp(T, H, mos E)] (in
. X \ 0 1/2+9)
|R[ - W|| 0(m -1/24")
The constants in the O-) notation may depend on n, D and also X >
Proof. This is straighforward consequence of Theorems 2.1-2.3. We omit the details.

---

Published as a conference paper at ICLR 2025
3.1 NON-ASYMPTOTI GUARANTEES FOR GAUSSIAN SKETCHES
The guarantees provided by Theorem 3.3 are asymptotic, in that they do not depend explicitly on crucial problem parameters (such X, H [): in particular, implicitly assume that H is well- conditioned. When E Rmxd is an i.i.d. Gaussian matrix, Si ( 1/m), we are able to provide non-asymptotic eπтог guarantees for our adaptive sketching algorithm. Our bounds are dimension-frees and remarkably do not depend explicitly on the condition number of H. The theo- rems below are the main technical part of this paper; their proofs are deferred to the appendix.
Theorems 3.4 and 3.5 respectively provide non asymptotic guarantees for Algorithms and 2:
Theorem 3.4. Assume that S has i..dd Gaussian entries. Let å E (0, Suppose that mo C(1 + (1g(1/)) for large enough constant C. W.p. 1 6, Algorithm outputs vn such that max(1.5d (A); mo) S m S (A)
2
(19)
Theorem 3.5. Assume that S has i.idd Gaussian entries. Let 8 E (0,1) and & > 0. Suppose Algorithm is run with 2 (A) and that moreover m C (11+lg(1/8) for large enough constant C. W.p. 1 1- 8, the output satisfies i - \ & where x X (1- 1 m (A)).
=²
Our final non -on-asympootic result concerns the concentration of the aggregate inverse Hessian esti- be the outputs of independent runs of Algorithm 2, representing workers producing estimates of the Newton step in parallel. Denote their average W = > (((k), and the error matrix:
mate used by the server. To wit, fix m 2 1.5d and let W1(1) ..... W
: (H W X1)
(20)
Later on, we will obtain error bounds for the quasi Newton method in terms of ||E]].
Theorem 3.6. Assume that S has i.i.d. Gaussian entries. Suppose that m 2 1.5d (A), and let q 2 1, 5 E (0, 1) satisfy 1/8 & exp(O(d)), S exp(O(d)). Let exp(- I = (()) be small enough. If moreover
1
dlog d log(1/6)
2 C q C
(21)
for large enough C then wpp. 1- 1, 1181
Remark 3.7. Up the logarithmic factor (which we believe is artifact), the dimensional de- pendence of q in (21) is optimal. To see this, suppose that H = 0, = 1, so W = s(t) gik), qm, therefore ||--1 I|] 2 when qm < d. In fact, E is Wishart (Gaussian sample covariance) matrix corresponding to mq samples, therefore by standard estimates
In particular rank(W)
EI -Il
mg V d
hence the dependence on is correct as well.
my4
4 CONVERGENCE RESULTS FOR THE APPROXIMATE NEWTON METHOD
In this section we present guarantees for the Newton method with approximate Hessian estimates. These bounds are relevant beyond the parallel scheme considered in this paper, and accordingly are stated in general terms. Let G d + R be convex and tiice-differeetiable We consider generally the Newton-type iteration,
0: = VG(O,-1),
0,-1 - ai WeSt:
(22)
(23)
where W. approximates the inverse Hessian, H, V²G(O,-1), W, H.¹ and or step size. Definition 2. We say (22)-(23) implement 1paccurate Newton method if the following holds. Define the error matrix, which measures how well We approximates the true inverse Hessian H,
:z W
(24)
Then ||E < foror every iteration t.
We denote TR
R (no). With this, the domain of G Dom(G) (a: F(x) < oo).
8

---

Published as a conference paper at ICLR 2025
For desired precision
0, namely when desires 0 with G(Ø) G( < 0
argmin G(Ø), let T(=) be the smallest 2 1 (a priori, if it exists) such that G( G( S for all t 2 T(z). This notation suppresses the dependence We consider a variation where the step size is chosen by backtracking line search, see Algorithm 3.
G and the sequence of step sizes at-
Algorithm (Backtracking Line Search)
1: procedure LineSearch (0--1-9t: 19: Wi; G(-))
2: Parameters: a, E (0, 1) line search parameters.
3: =
4: while G(0,-1 -- gt) > G(O+-1) - agi (aW gr) do
5:
== ba:
6: end while
7: return
8: end procedure
This paper focuses on the case when concocann,,,,,,,, which is a standard setting in the literature on second order convex optimization, cf. (Boyd & Vandenberghe, 2004). We remark that previous works, for example (Wang et al., 2018; Dereziński & Mahoney, 2019), derived convergence guar- antees for the approximate Newton method under different assumptions, namely when the Hessian of is Lipschitz. Due to space constraint, we do not present their results here, though they are applicable for our setting well. For self-concordant objectives, are not aware of similar results previously being written down in the literature, hence our results may be of independent interest.
Definition 3. We say univariate convex function G + R self-concordant if is thrice- differentiable and (x)" S 2G"(x=3// for any x in its domain. A multivariate convex function G Rd + R is self-concordant if any I-dimensional restriction Gs,y(t) G(x + ty) self- concordant.
When f is self-concordant, the objective G(0) = f(0) 110) in (1) is self-concordant as well, as the sum of self-concordant function is itself self-concordant, ef. (Boyd & Vandenberghe, 2004).
Theorem 4.1. Suppose that (a) G is self-concordant; (b) the Newton method is n-accurate with 1 < 1/5; and (c) ax is chosen by backtracking line search. For large enough numerical C 0,
To C G((o) - G((")
Ti(e) = CCoo To(E) c log(n/E)
ab
log(1/n)
where a, b are the parameters used in backtracking line search (Algorithm 3). Then
To (T) To + +T2(e) if E <<n
if E 2
T(e) <
(25)
The proof of Theorem 4.1 builds on and extends an argument of (Pilanci & Wainwright, 2015).
Interestingly, Theorem 4.1 suggests that the approximate Newton method exhibits two stages: (1) While the error is above the noise floor 2 1), the error decays quadratically fast. (2) Upon reaching the noise floor, & S 77, the convergence rate becomes linear; however, the linear rate proportional to 17. that is the error decays
(O(n))*. When is very small this decay rate is very fast.
Finally, Theorem 4.1 may be readily combined with Theorem 3.6 to yield error bounds for the entire end-to-end parallel method: under the setting of Theorem 3.6, w.h.p. (and up to polylog(d) factors),
1
d
17
mg
With this, let us roughly bound the runtime of our method with parallel workers. We focus on total flops. ignoring that in practice, deployment rounds are typically more expensive than computation. Every round, each worker performs, in parallel, (dd2)) flops, where m O(d (A)) is assumed for simplicity be roughly constant throughout. We bound the number of deployment rounds by T(e) = ((lg((1/)) log(1/n)), assuming that E is much smaller than yy. Consider two regimes:
Massive parallelization: When q 2 d, the Hessian approximation error is 77 O(1/vm) (bias dominated), and the total runtime is (dd m log(1/))/ log m). Note that this quantity de- creases with meheene it is always preferable to pick to be the smallest possible (which is

---

Published as a conference paper at ICLR 2025
Q
(x)
mg
Success rate
Avg. dim.
8.8
10
1.0
20
2/3
59.7
10
1.0
160
1/2
190.4
10
1.0
640
Table 1: Success rate of Algorithm 1.
3500
dl
3 1MAAI11 WAI dedketnng
1000
-MA WWY " debiiiingg
sketch dimansion
3000
1.5dl
1.5
1Mar1- W
7 -M W1 10Hiwing!
2500
4dl
800
2000
1.0
bies
600
1500
400
2000
0.5
500
200
1.5
2000
4000
6000
8000
10000
50
300
500
200
600
Figure 1: The bias proxy of the bias-corrected inverse Hessian estimator is substantially lower than without bias correction. Rightmost plot: ensemble (R); left and middle plots: ensemble (L). Left- most plot: the sketching dimension found by Algorithm I on ensemble (L). Shaded area: 20%-80% confidence interval: We take T 10 Monte-Carlo trials.
what we do essentially). Note when dx(A) < d, this runtime substantially improves on Odd log log(1/s)), the runtime of the exact Newton method (with no parallelism). Importantly, this bound does not depend on the conditioning of the objective, unlike in first order methods.
Moderate parallelization: when d/m S q < d, the Hessian approximation error is n7 =
O( (d/mg) (variance dominated), and the total runtime is O log(1/e)/ log(mg/d)). At the extreme end of this range, when q d/m, let us crudely approximate log(mg/d) R mg/d 1, so that log(mg/d) d) Z mq/d - Z mg/d. Thus, the runtime is roughly O log(1/s)/g). For mod- erate number of workers and realistic precision 6, this improves substantially upon the runtime of exact Newton, (d log 1og(1/e)) We again stress that this bound does not depend on the condition number of the objective, unlike for first order methods.
(d (1/s) =) log(mg/d))
Ol log(l/))/mqg//////
NUMERICAL EXPERIMENTS
We demonstrate the utility and validity of our results via numerical experiments. Due to space contraints, we defer some details and additional experiments to the appendix.
Our first set of experiments aims to demonstrate our adaptive sketching procedure (Algorithms 1-2) independently of any specific optimization context. For the first experiment, we consider diagonal matrices H with cigenvalues (hi...., ha) exhibiting polynomial spectral decay h4 = k- where a > 0 is a parameter. The dimension = 10 is fixed and large, and X = 1. The sketching matrix is i...d. Gaussian (in the appendix, we present results for two more i.i.d. sketching matrices). We run Algorithm 1, repeating for T = 20 Monte-Carle trials, and report the sketching dimension found, see Table 1. Remarkably, our algorithm consistently succeeds, namely finds sketching dimension that is within 1. (A) < m < max 4 (A)} per Theorems 3.3 and 3.4.
Next, we demonstrate that Algorithm 2 produces low-bias estimates of W. Note that the average Frobenius error ||W- W /df is a reasonable proxy for the bias of W, see Theorems 2.2-2.3. In the next experiment, we consider two random ensembles of matrices H. which we call here (L) and (R), and exhibits respectively slow and fast spectral decay; we defer the details to the appendix, due to space constraints. Figure 1 plots the bias proxy for increasing, starting from m = 1r (A), for the Gaussian sketch. As m increases, the bias decreases (roughly at rate 1/vm). We also plot the bias proxy for a naive estimator of W, without bias correction (λ = A); elearly, our algorithm has substantially smaller bias throughout.
We choose H be diagonal w.l.o.g. since the Gaussian sketch is invariant to orthogonal transformations.
10

---

Published as a conference paper at ICLR 2025
Our final set of experiments concerns the performance of the proposed optimization strategy in its entirety, focusing on the benefits of bias reduction. We present experimental results for real-world optimization tasks for both ridge and logistic regression- (Chang & Lin, 2011). Due to space constraints, the details appear in the appendix. We compare our method with approximate parallel Newton method where Hessians are sketched (m chosen by Algorithm 1), but debiasing is done; that is, Ā A. It is evident from our results, summarized in Figure 2, that bias correction accelerates convergence, often substantially.
publicly available UCI datasets
oegmart regrederi
Innsechere liége sqgeedar
(déie regreecee
vebide mndeet
Detenrig
Destaing Decliiing
eiing
biblsiing
Menton Stege
Nowdon Stops
Newtan Siass
Nenton Stago
Juberis Boystk regrsvel
Tear dbyztt repris en!
ruma flysst repoose
Iversioodees Boysst reprosien)
babusiin (bun)
2m
Dealaring
10t
(Dun)
233
31-
20"
28
Figure 2: Improved convergence of our parallel sketched Newton method with bias correction. The title of sub-figure corresponds to the dataset used. Repeating for T = 10 Monte-Carlo trials, the curve corresponds to the median and the shaded part to 20%%-80% confidence interval.
ACKNOWLEDGMENTS
We are grateful to Daniel Lejeune for inspiring conversations about his work (LeJeune et al., 2022).
Mert Pilanci and Fangzhao Zhang were supported in part by National Science Foundation (NSF) under Grant DMS-2134248; in part by the NSF CAREER Award under Grant CCF-2236829; in part by the U.S. Army Research Office Early Career Award under Grant W911NF-21-1-0222. in part by the Office of Naval Research under Grant N00014-24-1-2164. In addition, Fangzhao Zhang was supported in part by a Stanford Graduate Fellowship.
The work of Elad Romanov Donoho), and in part by the personal faculty funds of Prof. David Donoho, for whom he wishes to extend his sincere appreciation.
supported in part by the NSF under Grant DMS-181 (PI:
REFERENCES
Artem Agafonov, Dmitry Kamzolov, Rachael Tappenden, Alexander Gasnikov, and Martin Takáč. Flees: federated learning second-order framework via compression and sketching. arXiv preprint arXiv:2206.020000 2022. 4
Zhidong Bai and Jack W Silverstein. Spectral analysis of large dimensional random matrices, volume 20. Springer, 2010. 19, 28, 29, 39
Burak Bartan and Mert Pilanci. Straggler resilient serverless computing based polar codes. In 2019 57th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 276-286620 783..88. IEEE, 2019.
Sebastian Bischoff, Stephan Günnemann, Martin Jaggi. and Sebastian U Stich. On second-order optimization methods for federated learning. arXiv preprint arXiv:2109.0225 02388 2021.
Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004. 2, 37 37
Chih- Chih-Chung Chang and Chih-Jen Lin. Libsvm: A library for support vector machines. ACM trans-
actions intelligent systems and technology (TIST), 2(3):1-27, 2011. 11, 17
11

---

Published as a conference paper at ICLR 2025
Romain Couillet and Merouane Debbah. Random matrix methods for wireless communications. Cambridge University Press, 2011.
Michal Derezinski, Burak Bartan, Mert Pilanci, and Michael W Mahoney. Debiasing distributed second order optimization with surrogate sketching and scaled regularization. Advances Neural Information Processing Systems, 33:6684-6699 2020. 3, 4, 6
Michal Dereziński and Michael W. Mahoney. Distributed estimation of the inverse hessian by de- terminantal averaging, 2019. 9, 17
Nikita Doikov, Martin Jaggi, et al. Second-order optimization with lazy hessians. In International Conference on Machine Learning, pp. $138-8161. PMLR, 2023. 4
Anis Elgabli, Chaouki Ben Issaid, Amrit Singh Bedi, Ketan Rajawat, Mehdi Bennis, and Vaneet Aggarwal. Fednew: A communicaton--fficcent and privacy-preserving newton-type method for federated learning. In International conference machine learning, pp. 5861-5877. PMLR, 2022. 4
Zachary Frangella, Pratik Rathore, Shipu Zhao, and Madeleine Udell. Sketchysgd: reliable stochas- tic optimization via randomized curvature estimates. arXiv preprint arXiv:2211.08597, 2022.
4
Robert M Gower and Peter Richtárik. Randomized iterative methods for linear systems. SIAM Journal on Matrix Analysis and Applications, 36(4):1660-169)) 2015. 2
Robert M. Gower, Dmitry Kovalev, Felix Lieder, and Peter Richtárik. Rsn: Randomized subspace newton, 2019. URL https:
10874.
Vineet Gupta, Tomer Koren, and Yoram Singer. Shampoo: Preconditioned stochastic tensor opti- mization. In International Conference on Machine Learning, pp. 1842-1850. PMLR, 2018.
Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Science & Business Media, 2009, 1,3
Majid Jahani and Sergey Rusakov. Doubly adaptive scaled algorithm for machine learning using 2nd order information. In International Conference Learning Representations, 2022. 4
Eric Jonas, Qifan Pu. Shivaram Venkataraman, Ion Stoica, and Benjamin Recht. Occupy the cloud: Distributed computing for the 99%. In Proceedings of the 2017 symposium cloud computing, pp. 441- 451, 2017.
Peter Kairouz, H Brendan McMahan, Brendan Avent, Aurélien Bellet, Mehdi Bennis, Arjun Nitin Bhagoji, Kallista Bonawitz, Zachary Charles, Graham Cormode, Rachel Cummings, al. Ad- vances and open problems in federated learning. Foundations and trends@ in machine learning, 14(1-2):1-21) 2021.
Antti Knowles and Jun Yin. Anisotropic local laws for random matrices. Probability Theory and Related Fields, 169:257-352, 2017. 5, 6, 20
Daniel LeJeune, Pratik Patil, Hamid Javadi, Richard G Baraniuk, and Ryan Tibshirani. Asymp- totics of the sketched pseudoinverse. arXiv preprint arXiv:2211.037755 2022. 5, 6.
Vladimir Marčenko and Leonid Pastur. Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR- OOSR-Sbornik, 1(4):457, 1967. 4
Riley Murray, James Demmel, Michael W Mahoney, N Benjamin Brichson, Maksim Melnichenko, Osman Asif Malik, Laura Grigori, Piotr Luszczek, Michal Dereziáski, Miles E Lopes, et al. Randomized numerical linear algebra: perspective the field with eye software. arXiv preprint 3rXiv:2.022111777 2023. 2
Sen Na. Michal Dereziáski, and Michael W. Mahoney. Hessian averaging in stochastic newton methods achieves superlinear convergence, 2022. URL https: http:///axiivvo/// www.:iiiirrrr.ooooooooooooooo org/ abs//bb/////e/// /204. 09266. 4
12

---

Published as a conference paper at ICLR 2025
Jorge Nocedal and Stephen J Wright. Numerical optimization. Springer, 1999. 4
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation in pytorch. 2017.
Pratik Patil and Daniel Leleune. Asymptotically free sketched ridge ensembles: Risks, cross-
validation, and tuning. arXiv preprint arXiv:2310.04357, 2023. 2
Mert Pilanci and Martin J. Wainwright. Newton sketch: A linear-time optimization algorithm with quaaar-quadratic convergence, 2015. 3, 3.
Xun Qian, Rustem Islamov, Mher Safaryan, and Peter Richtarik. Basis matters: Better communicationn-effccaann second order methods for federated learning. In International Confer- ence on Artificial Intelligence and Statistics, pp. 680-720. PMLR, 2022.
Sashank J. Reddi, Jakub Konečný, Peter Richtárik, Barnabás Póczós, and Alex Smola. Aide: Fast and communication efficient distributed optimization, 2016. 3
Farbod Roosta-Khorasani and Michael W Mahoney. Sub-sampled newton methods. Mathematical Programming, 1 174:293-326, 2019, 4
Mher Safaryan, Rustem Islamov, Xun Qian, and Peter Richtarik. Fednl:
Making newton-type methods applicable to federated learning. In International Conference on Machine Learning, pp. 18959-8999911 19010. PMLR, 2022. 4
Ohad Shamir, Nathan Srebro, and Tong Zhang. Communication efficient distributed optimization using an approximate newton-typs method, 2013.
Jack W Silverstein and Sang-Il Choi. Analysis of the limiting spectral distribution of large dimen- sional random matrices. Journal of Multivariate Analysis, ((()):955309, 1995. 5
Joel A Tropp et al. An introduction to matrix concentration inequalities. Foundations and Trends in Machine Learning, 8(1-2):1-230, 2015. 31
Roman Vershynin. High-dimensional probability: An introduction with applications in data science, volume 47. Cambridge university press, 2018. 27, 34, 39, 40
Nikhil Vyas, Depen Morwani, Rosie Zhao, Itai Shapira, David Brandfonbrener, Lucas Janson,
and Sham Kakade. Soap: Improving and stabilizing shampoo using adam. arXiv preprint
arXiv:2409.11321, 2024. 4
Martin Wainwright. High-dimensiona statistics: non-asymptotic viewpoint, volume 48. Cam- bridge university press, 2019. 3
Shusen Wang, Farbod Roosta- Khorasani, Peng Xu, and Michael W. Mahoney. Giant: Globally
improved approximate newton method for distributed optimization, 2017. 4, 17
Shusen Wang, Fred Roosta, Peng Xu, and Michael W Mahoney. Giant: Globally improved approx- imate newton method for distributed optimization. Advances in Neural Information Processing Systems, 31, 2018.
Yun Yang, Mert Pilanci, and Martin Wainwright. Randomized sketches for kernels: Fast and optimal nonparametric regression. 2017. 3
Fangzhao Zhang and Mert Pilanci. Optimal shrinkage for distributed second-order optimization. 2023. 3, 4, 6, 17
Lijun Zhang, Mehrdad Mahdavi, Rong Jin, Tianbao Yang, and Shenghuo Zhu. Recovering the optimal solution by dual random projection. In Conference Learning Theory, pp. 135-157. PMLR, 2013. 2
Yuchen Zhang and Xiao Lin. Disco: Distributed optimization for self-concordant empirical loss. In Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference Machine Learning, volume 37 Proceedings of Machine Learning Research, pp. 362-370, Lille, France, 07-09 Jul 2015. PMLR. 3, 17
13

---

Published as a conference paper at ICLR 2025
Appendix
Table of Contents
A Additional Experiments and Details
15
A.1 Extension of Table
15
A.2 Details for Figure
15
A.3 Additional Synthetic Optimization Tasks
15
Details for Experiments on UCT Data Sets (Figure 2)
A.5 On the Merit of Centralized Data Access Versus Distributed Newton Methods
B Non-Asymptotic Random Matrix Theory Results
19
C Proof of Theorem 3.4
21
D Proof of Theorem 3.5
25
E Proof Theorem B.1
F Proof of Lemma B.2
28
G Proof of Theorem B.3
29
H Proof of Theorem 3.6
31
Proof Lemma
32
J Proof of Lemma H.3
33
K Proof of Lemma H.4
34
L Proof Theorem 4.1
37
Auxilliary Technical Results
14

---

Published as a conference paper at ICLR 2025
A
ADDITIONAL EXPERIMENTS AND DETAILS
A.1 EXTENSION OF TABLE
We run experiments with random sketching matrices with i.i.d. entries of the following types: (1) Gaussian (G): St,5 ~ (0, 1/m);; Rademacher (R): Stg Uni((+1//m) (3) sparse Rademacher (SR): S 0 w.p. - p and St,; ~ U /pm) w.p. p for p = 1/10. The following is an extension of Table 1, to include also the Rademacher (R) and sparse Rademacher (SR) sketch types:
dæ (A)
Sketch
Success rate
Avg. dim.
G
1.0
20
2/3
59.7
10
G
1.0
160
1/2
190.4
10
G
1.0
640
8.8
10
R
1.0
20
2/3
59.7
1.0
160
1/2
190.4
10
R
1.0
640
8.8
10
SR
1.0
20
2/3
59.7
10
SR
1.0
160
1/2
190.4
SR
1.0
640
Table 2: Success of Algorithm 1. T = 20 Monte-Carlo trials.
Our method appears to be remarkably robust, per the finding in Table 2, to the extent that the table may seem somewhat pointless include. We nonetheless provide it (and code to run this experi- ment) as a sort of sanity check.
A.2 DETAILS FOR FIGURE
We provide additional details for the experiments in Figure 1. the middle and rightmost plots, the estimator W is obtained by averaging q 500 local sketched Hessians, both with and without bias correction.
Ensemble (L): We use L 10-3 3. The matrix H H where X Roxd has the singular value decomposition (SVD) X UDV E Rand E Rdddd have orthonormal columns and D R is diagonal. U, V are sampled uniformly (Haar) from the corresponding Stieffel manifolds (have a rotationally invariant distribution). For the leftmost plot = 10 and d varies. For the middle plot, we fix 10 The diagonal entries of D are random, and satisfy
Dr. = (0.9 + &x) k/2 & ~ NNO,
Ensemble (R): We use À = 10-5 Similarly to ensemble (L), H = -x X for X E Raxd with the singular vectors of X similarly uniformly random. We use n 10 d = 103 The singular values of X are non-random, with DE. (k/d)².
A.3 ADDITIONAL SYNTHETIC OPTIMIZATION TASKS
In this section we present additional experiments on synthetic optimization tasks, aiming to test the performance of our proposed end-to-end parallel optimization method. Specifically, consider ridge and logistic regression. We test both Gaussian and Rademacher sketches (see Section A.
Distribution of covariates (design). We generate designs E Raxd with n = 10 and = 500. Denote the SVD X T The singular vectors U paxd V Rdxd are uniformly random. The singular values decay exponentially, DR 0.99
Parameters. In all our experiments, À 10 3 For experiments with ridge regression, the number of workers is q = 10; for logistic regression, it is = 20.
Ridge regression. The responses are generated according y: x. T + NO 0.01), where
the ground-truth regressor 0* Rd is generated according to 0* (O, 7). The unregularized
15

---

Published as a conference paper at ICLR 2025
Gaussian Sketch (Ridge Regression)
Gaussian Sketch (Logistic Regression)
Rademacher Sketch (Ridge Regression)
Openatt 5x0
10
Optimakiy
33
Newrion
Newton Steos
Nevitn Shaps
Monton
Rademacher Sketch
Sparse Rademacher Sketch
no debiasing
debiasing (ours)
Figure 3: Convergence of the parallel Newton method with Hessian sketching, synthetic opti- mization tasks.
objective F(0) is the L2 loss:
F(0)
Σ
y1) = ||0 y|1
n
where 21...., In E IRd are the rows of X. Recall that aim to minimize the regularized objective G(()) = F(0) + 1012.
Logistic regression.
We generate labels according to = (sign [ 0" + (N 10t)| + 1)/2 E {0,1),
where 0* (0, The function F(0) is the log loss:
(1+e-176 + (1 1o (1 -1- 2
1
F(0)
log
Figure 3 plots the error G(O,) argmin, G(f) over the iterations t of the algorithm. Each plot is obtained by T 10 Monte-Carlo trials; the shaded area corresponds to a 20%-80% confidence interval. evident that bias correction improves the the convergence speed of the algorithm over the uncorrected (A A) alternative.
A.4 DETAILS FOR EXPERIMENTS ON UCI DATA SETS (FIGURE 2)
tridines
Cettoayy dee
38
4-1
397
30-3
20
3
Meatten Strpe
Kestin Steps
Newetnn Stage
Meattn Strpe
pyrin
3.5
Q-10
0.5
044
Nanton Stare
Nowton Sicte
Monion Dogg
Noweos Stape
Figure 4: Convergence rate for optimization tasks on additional UCI data sets. Shaded area corre- sponds to 90%-80% confidence interval.
16

---

Published as a conference paper at ICLR 2025
We provide additional details for the experiments reported in Figure 2, and present results for addi- tional data sets from the UCI data repository (Chang & Lin, 2011).
We rum experiments directly the available data, without additional preprocessing. The objec- tives, for ridge and logistic regression, are as described in Section A.3, but now the covariates and reponses/labels are given and not generated.
Parameters. We use À = 10-3 3 for all experiments. We vary the number of workers between 4 E (1, 10, 50). For every data set, we run T = 10 Monte-Carlo trials.
Results. The experiments are reported in Figure 4. The plots on the top correspond to the opti- mality gap G(O,) argmin, G(#), G' being the regularized objective, function of the iteration number f. The bottom plots are the inverse Hessian normalized Frobenius error |[, - W.1l / We see that as q grows larger, the algorithm tends converge faster; though the benefit of further increasing q tends to diminish when is large.
A.5 ON THE MERIT OF CENTRALIZED DATA ACCESS VERSUS DISTRIBUTED NEWTON METHODS
A.5.1 BASELINE REVIEW
Here, include comparisons with existing distributed Newton's methods, which operate under slightly different paradigm as they focus on data splitting among workers. These experiments are presented to illustrate that our debiasing approach, which operates under centralized data access within a multi-worker setting, achieves superior acceleration for Newton's methods.
We provide review all baseline methods we comparing against. Consider any loss function F(x) where x is the desired parameter. Given a dataset with n samples, we first split the dataset uniformly across q workers in total, which results in each worker holding n/q² data points. Then, each worker (i E [g]), computes the i-th component of the Hessian (V²F), F); only based on its local data, where the full Hessian is Σx= PV² E==(""") + XL. The difference of all baseline methods lies only in the process of finding approximate Newton directions in each iteration. In iteration t, GIANT (Wang et al., 2017). Determinantal Averaging (Derezinski & Mahoney, 2019) ("Determinant" be- low), Optimal Shrinkage (Zhang & Pilanci, 2023) ("Shrinkage" below), and DiSCO (Zhang & Lin, 2015) methods all take the next Newton step as (((+1) = x(t) ((()))((()
(+)) g( where r7(*) is the step size found by line search current iteration, (()) the global gradient and AT0) is an estimate of the Hessian at (!!)) Specifically,
. GIANT takes estimated Hessian inverse average over all local Hessian inverses, i.e., h = = Σ. ((V F); -
Determinantal Averaging exploits the expectation identity ((de)())(((AA)) E(det(A)A-)) which is valid under mild conditions for invertible A and approximates
the Hessian inverse by - Σ, {dde((("P))+AA)
Optimal Shrinkage method approximates Hessian inverse with A(0 ; Σ (V" where dx is the effective dimension of (V²F)t and
it requires n/qq > d=. In the experiments, follow the approach in (Zhang & Pilanci, 2023) and take d² used by worker as effective dimension of local (V²F),
DiSCO approximates Hessian inverse by using the first worker's statistics A0) = ((V
DANE considers slightly different update rule which involves solving
3 || 200
x (++1) argmin (.(x)
マ | ((VF)(z()) + xx()) (+))
x
2
Here, we assume is an integer. For real-world datasets, any remainder when dividing n by q is disre- garded for simplicity.
17

---

Published as a conference paper at ICLR 2025
then it takes the next iterate x((+1) Ei experiments.
(e+1), We take 1 and 3 =2 in our
We compare these methods across two scenarios, namely ridge regression and logistic regression, detailed in the following subsections.
A.5.2 RIDGE REGRESSION
Due to the varying sizes of samples in each dataset, we set number of workers to = 10 for "bodyfat", "splice", and "pyrim" datasets; q = for "segment". "eunite2001" and "iris" datasets; q 30 for "mg" dataset. For each method, we repeat the optimization process over 10 random seeds, with the data split (or sketch) redrawn for each iteration, and plot the mean and one standard deviation of the logarithmic optimality gap. Figure 5 shows the convergence progress for different datasets. Since for quadratic loss, DANE and GIANT coincide, we thus write "DANE/GIANT" to represent both methods. It can be observed that our adaptive de-biased Hessian inverse estimation accelerates convergence significantly and converges the fastest in almost all experiments. The most significant improvement over baseline methods happens on "eunite20011, "bodyfat", and "pyrim" datasets, where all baseline methods take at least ~ 20 iterations to reach certain accuracy and our method takes S 5 iterations. On "segment" and "mg" datasets, though most baseline methods are already good and converge within 10 steps, our de-biased Hessian inverse estimator still provides improvement in this case. Moreover, our method exhibits the smallest variance, which highlights its robustness and stability.
33"
BANDDIGGT
CANUUDANT
Ditermetant
Strincage
Strieaage
Shrinlgge
Shrriage
hanton Stupa
hewbon Stapa
Nowton Stopa
Nexton Stare
162
1.2
35
Ditammmnntt
2
2-
16
10"
Desaminant
Merkayg
Detamrinntt
labiiiing (Cun)
21
Declasing (Oun
Catliaing 1Ound
Dabkiing 10us)
23-
0
74
36
28
15 124 355
Monion Staps
Nowdon Stuga
Mowion Siepa
Moxtan Stote
Figure 5: Comparison of Distributed Newton Methods for Ridge Regression: Our method demon- strates the fastest convergence and the lowest variance across nearly all experiments, highlighting its effectiveness and robustness.
A.5.3 LOGISTIC REGRESSION
Here set q for "liver-disorders" and "ionosphere" datasets; 10 for "heart" dataset; q 20 for "diabetes" and "sonar" datasets; 50 for "germannnumer" numer" dataset; q = 250 for "fourclass" dataset; q 500 for "svmguide3" dataset. For each method, we repeat the optimization process over 10 random seeds, with the data split (or sketch) redrawn for each iteration, and plot the mean and one standard deviation of the logarithmic optimality gap. Figure 6 shows the convergence progress for various methods on different datasets. We observe that the convergence behavior for logistic regression across different methods is more challenging to achieve compared to ridge regression. This difference may be attributed to the non-quadratic nature of the logistic regression objective function. From the results, one can observe that our de-biasing method always converges faster than all baseline methods and the improvement is most significant on "german ggerman_numer", "ionosphere" and "sonar" datasets where our method halves the iteration number compared to the second best method. Moreover, DANE diverges on datasets such "german "german_numer" and "fourclass", Determinantal averaging encountered NaN values during optimization and we thus skip it for "ionosphere" and "sonar" datasets. This may be due to the numerical issues in the determinant caleulations. Although most of the methods suffer from high variability, our method still behaves stably with little variation. This again proves the robustness of our Hessian inverse de-biaser.
18

---

Published as a conference paper at ICLR 2025
dabvtes
Feet
feerdleerde'
IRANT
SMANT
DIPAT
Bettrn
16-
Drttiiigg (bung
20
22
To JTTIBAA
11:
1P
195
Optmata
Los opptenai
Desaaiing 10uno
26
20
20
35
23
25
25 20
Figure 6: Comparison with distributed Newton methods for logistic regression. Our method con- verges the fastest and exhibits the smallest variance for almost all experiments, which shows its effectiveness and robustness.
B NON-ASYMPTOTIC RANDOM MATRIX THEORY RESULTS
The proofs of Theorems 3,4-3.6, rely non-asymptotic results from random matrix theory that that we state in this section and prove later on. The results below concern the behavior of random matrices SHST S SH where H and S max is an i.i.d. Gaussian matrix Sig ~ (00, 1/m).
The following result concerns the smallest eigenvalue of the matrix SHST E Rmxm. Theorem B.1. Suppose S has i.i.d. Gaussian entries. There are constants C1, C2 C, C such that the following holds. Set any rj > 0. If m S (T) then w.p. 1 - Ce-om
Amin Main(SSS (SHST) c2ll m
d
(n)
1
The following two results concern the closedness of the empirical (companion) Stieltjes transform
$(z) == tr (SHS - -21)
(26)
m
to the deterministic function s(z), the solution of the Marchenko-Pane Pastur equation (10).
The following concentration inequality well known, cf. Bai & Silverstein (2010). Lemma B.2. Suppose that S _
mm [ril-.. ... [mm]] I where the rows Fi...., Tm I2d are independent
random variables. For any < and 2 0,
Pr([S(z) RS(2) 2 1) < 2e 2
We provide self-contained proof for completeness, see Section F.
Next, we show that EE(() approximately satisfies (10), with explicit non-asymptotic error bounds. Note that (10) can be be equivalently rewritten (assuming 2 # 0) as
s(z) - 1 m
(1-
dn
L
(27)
Theorem B.3. Suppose that S has i.i.d. Gaussian entries. For any 2 < 0, BS(2) satisfies:
BS 2)
1
du (I(E(((=))) + Cm (m; H)
(28)
z
m
where the error term em (m; H) satisfies
8 Cm (2; H)
(29)
a 2 1 m
19

---

Published as a conference paper at ICLR 2025
with am H2; B Pm H) defined as
= > [2)(
(30)
= m
d (1 (=))
(31)
The proof of Theorem B.3, given Section follows along a well-known computation technique in random matrix theory, and explicitly leverages the orthogonal invariance of the Gaussian dis- tribution. Note that the error terms at 3 depend on [Eš(z) itself. This feature, which somewhat complicates the proofs of Theorems 3.4-3.5, is generally not mere artifact of the calculation, see e.g. the bounds of (Knowles & Yin, 2017) under asymptotic conditions.
20

---

Published as a conference paper at ICLR 2025
C
PROOF OF THEOREM 3.4
Denote for brevity, per the description of Algorithm 1,
5
20 == 12
(32)
The proof of Theorem 3.4 consists of two parts. First, we show that if 1. (A), then with high probability (w.h.p.) Sm (=o) < 1/x, hence such m would not be returned. Then, we show that when m 2 2d then w.h.p. ( 1/A, which would guarantee that the algorithm must return such m (if it has not halted before). Consequently, the output of the algorithm must satisfy m S (A).
The next lemma shows that Algorithm does not return m which is much smaller than dæ (A). Lemma C.1. Suppose that has ii.dd Gaussian entries. There a (small) constant co such that if m cody (X) then
Pr($(=o) > 1/2) S Cm
Proof. The lemma is consequence of Theorem B.1, which bounds the smallest eigenvalue of SHST. Set cp min(er, 12 for cı, e2 from Theorem B.1. Assume that m < cod (A). Since in particular < ei crd (A), the following event holds w.p. Ce-em
1
S(20) S
Amin (SHST) 20
C2 m
dH (A) x 12
Moreover, since c2/c0 2 7/12,
7
C2 dH (A)
C2 Co
12
Hence, under this event, $(20) 1/A.
Next, we treat the regime of moderate m, namely when cod (A) S m < (A). The main tool here is the non-asymptotic estimate Theorem B.3. The reason we have to assume m Z dn (A) is that only then we can guarantee the error term is controlled.
The following simple lemma useful.
Lemma C.2. For any 0 and y N 1.
(H())) s (p) S 7dw(pp).
(33)
Proof. The inequality dur (1r(7)4) < d (1) is trivial. As for the other one,
d tr(H(H + - = ir( ur(H( - + = (x) 2 + = ddr(u),
where (+) follows because H ž 0. hence 1 pl I H + pll and therefore H ul) 1 -
H(H + -
0
Lemma C.3. Suppose that has i...d. Gaussian entries. Moreover, suppose that cnd (A) S m <
1.5dg (A). There is numerical constant M such that if moreover 2 M then
14 15 X
E$(20)
In particular, Pr(š(20) > 1/X) S Ce-cm
Proof. Note that the second part follows directly from the first part and Lemma B.2; we therefore focus bounding Eš(2g). Suppose it were the case that (2)) n) where 7 S is some constant. Apply Theorem B.3 with 2 = 2o. Let us bound the parameters a, 3 used for error control. First,
5
1
5
>
11
2
12 7
(34)
12
21

---

Published as a conference paper at ICLR 2025
Next, since E5(2))
1 Fol
12 we have 1/E5(20) 2
1.2 so
B = m
1 n (1/B5(za))
v ddr
))
(+) 121
d 5 m
(14)
12
(35)
m 12
5 cD
where (+) follows from Lemma C.2 and (**) from the assumption m 2 (A). particular, the enror 6m (=o H) in (29) can be bound as 6m(20, H) = Ogl x/m) Furthermore, as ES(zo) 2 77 by
assumption,
B
11 dn dn (A) > z". m
(36)
m
where we used Lemma C.2 and the assumption m < 1.5d (A). With these estimates in mind, consider (28),
1
1
12 (1 3 20em (20: H))
1 12 (1
2
1 Onl
(2()
(1 3 m (z0, H))
vm)) ))
20
X
3
m
As ES5m (20) 2 3, we deduce from the above
12
2 3
77 5
(1
+ O,
Vm ))
5
so that rearranging,
12 (1+ O ))
7 S 13 (1
13
By requiring that m 2 CG for large enough Cy,, we can ensure that the error term is as small a constant as we like, so that necessarily (for example) n S 14 15
Next, we show that if the algorithm has reached m 2 2dz(A), then it will halt in that iteration. We start with two helper lemmas.
Lemma C.4. Suppose that [SS 5) = T. Then for all = < 0,
Es(2) 2 1 (1 du (-2)
(37)
(Note that since IES(2) 1122, 2, we have d=((=) 2 d (1/ [B(()), so that (37) is perfectly consis- tent with (28).)
Proof. Recall that the matrices SHST E mx and / SH1/2 have the same non-zero eigen- values, and therefore
ES(=) == L EErrSHST - 1 -1 (d - m)
SH
m
m
m
Note that since < O, the matrix function A + tr(A =1) is convex on the cone of PSD matrices A 1 0. Thus, by Jensen's 's inequality,
22 SH2/2 =1)-1 2 m
tr Lr SH 21) m
tr(H 21)
m
so that
E$(2)
1 2 m
1 (d (d-m)
- tr((H -2 -1 ---11 - 1 - 1 (1- 1 2
(--d-(-2)
L
The next lemma formally establishes that the inverse function of s(2) is increasing on appropriate interval, As in (11), denote
(0, 20) I,
8) (1 (n(ss))
1
(38)
m
so that for all & < O, ((((2)) z. Note however that in general, Range(s) may only be a proper subset of (0, 20); nonetheless, (38) is well-defined.
22

---

Published as a conference paper at ICLR 2025
Lemma C.5. Let μs O, and suppose that 2 (p). Then v is increasing (0, 1/p).
Proof. Write
((()) =
tr(H(I (1+ + +H)
m
so that its derivative is
1
('(()
tr HH + sH) -)
22
m
1
tr(H2 + 1/5)
g2
m
1
11- dl(s)/)
2
82
1
m
t + 1/3)
= z²
m
(39)
where the inequality follows since that (1/8) < du (M), hence (V) 0.
HHH + 1/8) T. When < 1/p,, we have 1/8 l so
We now prove:
Lemma C.6. Suppose that has i.i.d. Gaussian entries, and moreover that m 2 (d) There is numerical M such that whenever m M,
24
ES(20) 2
23
In particular, Pr(š(20) I 1/A) S Cm
Proof. Suppose that (=)) S fー x for some (0,1); we shall show that r) has to be somewhat large so to not generate a contradiction. First, by Lemma C.4 (with = --),
Es(za) EES (-(- 1
1 dz (A))
X
where the last inequality follows since m 2 (d)( Of course, this lower bound does not suffice for our purposes: we want bound which is > (strictly). It does allow us, however, control the error term in Theorem B.3 at zo (specifically, to lower bound a). We can bound the parameters a, 3 in Theorem B.3 as
11
5
Q: = 2a ES
20
2
12
2
24
1
1
5
12 1
6 5
8 - dH (1 (1/B5(z00) <
-
-d
dH
X)
<
d
m
m
12
m
where, in upper bounding 9, we used Fš(zo) implies
- = and Lemma C.2.
1
With these, Theorem B.3
(5())) = 1
d (1/ES,. + 0
m
A/m
Multiplying by 2a and dividing by EE((2))
1 11-dx((/EB((a) + +0 A +0(-
[
Zy = ES(())
m
(E$(2))
+0(
mm
where in modifying the enror term we used Eš(z0) Z x.
Note that 2 2d (A) implies 2 ( by Lemma C.2. By Lemma C.5, this implies that 1// is increasing 20, Thus, if we further assume that 7) < 1/2, then V is increasing (0, - x)
23

---

Published as a conference paper at ICLR 2025
and (EE(z)) I | - x)) Using this with the previous display, and writing 20
-
1 - }
- V
+
O
12
1
= =-(1 - 1i}
dH((1 - -))a) + +0
m
mm
((1 -
7)
(1
1 1
dH d(((A)
+ +0(
1 - mmm
mm
3 (11 7)n
1
1 - 1 1-n n
+ +0
mm
-1 1 11 0
Vm
where (+) uses dæ<((1 - 1)1)
d- (A) (Lemma C.2) and (**) uses the assumption m
2d (A). Dividing by (-)),
5 0
12
Vm
For large enough m 2 M, the 0()) term is S 21 therefore 17 2 J.
That is, we deduce that
21'
necessarily (Z)(2)) 2 1- X X
241
We are ready to conclude the proof of Theorem 3.4.
Proof. (Of Theorem 3.4.) Let m, 2/m 0, D,..... be the value of at iteration j of the algorithm. Set
J 11.50115d (A)/mo)))
J [log (A)/mo)l.
chosen that: J- is the smallest j such that my 2 1.5d (A): 2) 2d) I m,++ < (A). Thus, to establish (19) it suffices to show that the algorithm halts at an iteration j satisfying J S j S J Accordingly, the error probability can be bounded
Pr((19) does not hold) <
Pr(Sm, (zo) 1/A) + Pr(Sm, (zo) 111).
Assume that mo 2 M is large enough. By Lemmas C.I and C.3, for all j J - 1, Pr(Sm, (2o) > 1/A) I Cm-om, By Lemma C.6, Pr($m,, (20) S 1/A) & Ce-emit, Thus,
Jー
Pr((19) does not hold) -
Σ
Ce -2- + Ce at ma <
Σ
Ce-c2 -e--e2mm I e como
1=0
5=0
We deduce that if ma 2 Log(1/5), then the error probability
S 8.
24

---

Published as a conference paper at ICLR 2025
D PROOF OF THEOREM 3.5
Recall the (deterministic) Marchenko- Pastur equation (10); rearranging yields
2 s(z)
dd((/s((2))
(40)
m
so that setting s(- 1/A yields Á X(1 d (A)). Further recall that do not have access to s(z), but instead to its random counterpart 5(z). We know that its expectation EE(2) approximately satisfies the Marcheko-Pasiee Pasttur equation with an error term (Theorem B.3), and that we have concentration around the expectation (Lemma B.2). The proof of Theorem 3.5 down to controlling and propagating these error terms.
Lemma D.1. Suppose that m 2 1.5d (A) and m 2 M for large enough M.
unique root z < O such that Es(=*) exists, and 2 < -cx small enough c. > 0.
Consider a sufficiently small neighbooiooo of 2", I* = [22 -1,22- 2" + 17) for small enough 7) > Then,
1. [Es(2) 2 1/ for all all E
2. [B((2) E8(2') 2o
for z
Proof. The existence of 2², well as the bound (2) 2 n can be shown by an argument similar to Lemma C.6; we omit the technical details. Specifically, one can show that for small
enough c, Es( (A) 2 (1 + c) Since Es(-) is continuous increasing, and (E)(-00) 0, root E$(2*) x exists and 2* -c.
Next, we prove Item above. We have
x'
1 x
|Eš (2) Es, (z') E- m
(S. 21) (1) 2/1)~ 211
2'1.
22
since
|=|, こ X, as 2, o' E In-
Item 1 follows immediately from Item 2.
We need to establish a result of the following kind: if & < 0 is such that Es(2) Ro x. then necessarily 2". This necessitates lower bound (rather than upper bound, such as given in Lemma D.1, Item 2) on the derivative of [=)()) around z².
Lemma D.2. Suppose that 2 1.5d (A) and m M for large enough M.
For a sufficiently small -eighhhoiooo of 2", =" E TT,
1
[E$(2) EE(('') x - O X
(41)
Proof. By Theorem B.3, having established that R$(2) 2 1 for all E 1*, uniformly for = T*,
Eš(2) (1(11 1 Hp d1 + O A/m
uniformly for all E I*. Dividing by Eš(2) and multiplying by
z = ((E((2)) + O( vm
(42)
where (s) - ( (1/s)) is the inverse function of s(z) (see e.g. (38)). Note that for 2, I", [E(() Es(2') I < 1z - 2 ). Thus, to conclude the proof of the lemma, suffices to give an upper bound on the derivative (or Lipschitz constant) of in a small O(1/A) neighborhood of FE(='))
25

---

Published as a conference paper at ICLR 2025
Indeed, if s (1 =) for small (1), then by Lemma C.2, I (A) - d1((/s) 1 - Since # (A) < (1.5)-1 - 0(1) - (1/((/A) = O(s). Thus, for $1, 92 in a sufficiently small 0 ( neighboroodd of t. m d
= ((1)))
and so (v1) ((22) = O |1 $2)). Using (42), we deduce that for 21, R2 in small enough
(()) neighborhood of 2*.
12 - 22| = (ES((2)) )) - (ES(22)) + 0
Vm
2 X2 ==) - F=5((2)) + 0
m
The desired estimate (41) follows dividing by X2.
We are ready to conclude the proof of Theorem 3.5. Let & > 0, 5 E (0, be, respectively, the precision on confidence parameters. Suppose that = (1) is small, such that [a* -2λe, 2" +2Xe] S I", therefore setting
=_ = A
(43)
we have =1 E I*. By Lemma D.2, Eq. (41), provided that m 2 1/6² is large enough,
1
1
Es(x1)
S X Cie + +0
S X
C2E
A/m
A
Bs(24) 2 x
+ CiE 0
2 X
+ CZE X
A/m m
By Lemma B.2, if m 2 1 log(1/0) 16) is large enough, then w.p. 6, |5(=)) - (l.
22
Es(2t)) )I E c/2 On this event,
1 1 バ' s(z2) 7 2
1 1
1
1 1
$(2)) 2 x 2 C2S
ス² > x
In particular, the output Å of Algorithm 2, being the solution so 1 (-2*)] S e. Finally, by Eq. (42),
s(- 2)= = 1/x, satisfies 2n < - $ 27,
A/mm )=
x
オ = (1/A) + O
= - + 0
= x + O(As)
Vm
Vm
whenever m 2 1/E² large enough. Since Ā 2 A(1 - (1.5) 1) = 1/2²/0x = 22(2), this implics
Å \1
O(E)
Ā
Thus, Theorem 3.5 is proved.
In the sequel, the following bound on Σ| - Àl will be useful.
Lemma D.3. Suppose that m 2 1.5d n (A) and m 2 M for large enough M. For any p >> 0.
E -I0 S mp/2m-1122
(44)
Proof. Let to > 0 be small enough constant, such that [2* - Ato, 2" + Atol C I*. By the above,
Pr(IA- - À 2 tod) Pr(-A & [a* Xt., 2* + Atul)= < Ce-cm Set (t) 2² = Xt, so that by
the above, (r; < 1/A), (r(s > 1/A) < Ce-cmt for all < to- We have
Pr(Ii- 2 tA) = Pr(-- & |- (t), (f)) Pr (( (2) (t)) < < (= (t)) 1/11) < 2Ce -eme?
Thus,
EJ Åpr = ElIĀ - ÀPr A - + BlJ - X 1 A < EIパĀ 7 1 CXee- -CN'
and
EIJÀ- X = // t0
1Pp( - - 2 xt)at S &DAP/ 1/ p-le-om² dt < CpXPm m/2
And so, BAÅ - ālp S -1 mm + amm = 0(1
26

---

Published as a conference paper at ICLR 2025
E PROOF OF THEOREM B.1
By definition, if H 3 (in PSD order) then SHST I SHSS, in particular Amiu (SHST) > Amin (SHST) Choose H = HuH + I)=,, noting that ||]|| < 1/P can be bounded irrespective of ||]]
We now lower bound Amin (SHS) using a standard net argument, cf. (Vershynin, 2018). For any
> 0, let We be e-net of 8mh1 of minimum size. One show that [ (1 (1+33 -) We have
Amin (SHST)
min x' SHS min T SHS - 21|SSS |e.. ||||==1
(45)
wENe
Note that for fixed E -5 y ST has independent, mean zero, sub-Gaussian entries with maxigiga lball = (1) and variance 1/m.
We first give a bogh-probability upper bound on ||SST ||. Let W1/4 a 1/4-net of minimum size, so that
||SHs
max T SHSx S max x TSSSS* +2 2 7 #EN1/a
||]|=|
hence ||SHST| S 12maxxENi, x² Sm HnS. By the Hanson Hanson-Wright inequality, Lemma 95 (see
also (Vershynin, 2018, Theorem 6.2.1)),
Pr(]ISHS 2 2m tr(H) + Wi/a1 sEN1/4
max Pr Pr(x TTTTTHHe SHSr 2 m tr(E) + t/22)
| | 1
m2t²
mt
< 2gm2ex)
-emin!
Note that
1
tr(H)= dd (1/M))
| $ 1/4
| S V (r(Ht vdu(1//).
Set / C1 max
Ir, | large enough. For such choice, w.p. Ce
2
1
||S v m
tr (T) + C1 max
M |1|,, |1
Vm
d (1/M)
V
d2in
mmxx()
(46)
m
1
Next, we give a high-probability lower bound minneN, x. SHST x in (45). Again by Hanson- Wright,
Pr(min min x T SHST 2 m tr(#t) - t) - 1 2 E) zeN,
2cxp -eminn
m222 mt
|
H
so that w.p. 1 Cemem cm for small E 0,
(1/ log(1 + 1/E) 1
max(1, (H(1/))(H) (47)
min x TTSS x x> NEN,
m
Inserting (46)-(47) into (45) implies that w.p. 1 Ce-om
Amin (SHST) (1 - (())) 1 dn/(p/) O(E + +og((++ max / (dn()/ }. (48) /
P
m
By choosing 0 sufficiently small numerical constant, can deduce the following: there are numerical constants c, C2, > 0 such that du (1/p)/m > 1/01, then w.p. - Ce-em Amin (H) 2 c2 dn(1/M) To recover the claimed result, use μ 1/n.
27

---

Published as a conference paper at ICLR 2025
F PROOF OF LEMMA B.2
The following argument standard (cf. (Bai & Silverstein, 2010)) and presented for completeness.
Let Tm Rdd the rows of S; that is, ST = [ri]... |m..../mm....m...... Let
mm S SH
H 11/2 1/2
(49)
This is sample covariance matrix, corresponding to samples with population covariance E[Em] H. Recall that Em E Iddd and HS R have the same non-zero eigenvalues. Denote
(m( (Em + =1)-
(50)
that
5(z) m tr(Pm(()) + m m=-((dd m)
(51)
Central to the proof is the following leave-one-out decomposition:
1 P Pm (2) /
(=) = =P (2)
(52)
1+ H1/2 1 (2)
(=) = (2.m H rrrt 1/2 + =1)=
k 1, m.
(53)
The above can be readily shown by the MMorison lemma (Lemma M.1).
We shall prove Lemma B.2 by standard martingale concentration argument. Let F& be the o-algebra generated by rı,..., Tke = 0, 1, m....,,,, m, and E = E[- Decompose into sum of martingale differences.
$(2) -[((2)] =
EEc Egk-1)(tr Pm (z)]
Rm
m
(Eck Egk-1)
(54)
k=1
where, using (52),
m r (P
1
Dm,((2)
(55)
1 m +
r H1/2Pmm (2)
Note that 0 3 (x) I. hence
1 T H (3) k
D (m)((2)] S
1 1 112 (2) (2)
-2
Thus, (54) is of bounded martingale differences |EEg Eg () (m,, < 2. To conclude, use the Azuma-Hoefiding inequality (Lemma M.4).
Remark: The above calculation (regarding the boundedness of (2)) is, essentially, the proof of the well-known low-rank perturbation bound for resolvents, Lemma M.3, which we shall use later.
28

---

Published as a conference paper at ICLR 2025
G PROOF OF THEOREM B.3
We implement well-known computation technique in random matrix theory (cf. Bai & Silverstein (2010)), while carefully keeping track the error terms.
We rely explicitly on properties of the Gaussian distribution. Denote the eigendecompposition H Ei=1 TeVeti and
Hy H TeVe St Sug-
for L = .....dd Note that ,1,...,8d N(0, are independent (since uh..... are orthogonal).
Recall Eq. (51), so that
1
1
S(2) =
(d - m.
+
tr( SH - =1)=
1
1 -
1
m (d m)
0 12 S11 =1)-11z
(56)
m
We now compute an expression for 1g (H Sm H -1) the f-th diagonal element of the resolvent ( S, H1/2 - =1)-1 written the population covariance eigenbasis V [val.../va)
Upon a coordinate permutation (where l becomes the first), we have
Tells2|12 -
SH. 1/2 V 21
V (H =
! se H ST SH 1/2
By the block matrix inverse formula (Lemma M.2),
ve T ST S/2/2 =11- = (V 12 1/2 =1)V -
(Tellsel2 T1s 12 (11 ST SH - =1)= H!! //22STT ST 88
(Tellsell2- - - (SH - - S se) sa)
(1 + (S -2]) -1
(57)
where in (x) we used the identity X X(XXXXX) x ((XTxxx X which holds for any matrix X and analytic function f(-).
We now wish to estimate the expectation of (57). Denote
D: sr(Sl -21)- E8(2)
Dz,1 D2.2 + Dt.3
where
Dr,1 = SISS - -21) 1 tr - -1,
m
1
1
D2.2 tr(SH ST -21)-1 tr(SHS =11
m 1
m
D3.3 Tr(SHS -=1) Eur(SHS =11)-1
We have
TrD;
do ( 1/2 - =1)- = - (1 (1+ 71())
(1 (SH -1)-1ss)(11 Ts(s
29

---

Published as a conference paper at ICLR 2025
so that using (56),
Bs(z)
- (d m)
1
1
+
1
Σ
Bug S SH -11)-1u 'e
m
(d m
(1 + T TT)) + (m)
m
2
t=1
1 1
+
[I
(I BS(:) em
2
2 m
dn (1/ (I/BB(=)) cm))),
(58)
m
where
m
TeDi (1 - (1) ( T+BB((2))
Cm
(2)
E
2m71
(1 +
We may now bound
len(2))
TF.
max
D
m|2 t=1
1 + (x)S(2)
15t5d
d (1r(1/BB(z))
max E E|
(59)
1= (ES(2)
Isesd
We now decompose D| - E De,1
De,zl Dt,). By Lemma M.8 (Item 1),
E
Dr.1 < E 11 El((Sm H = S llym
By the low rank resolvent perturbation lemma, Lemma M.3, almost surely
De.2l S 7
Finally, by Lemma B.2,
E De.al = B[$(z) ESm (=)l S |z|ym
Thus,
Jem (z)l= =0 (S(E(() d (1 (I/E3(z)) 12l
(60)
and so the proof is concluded.
30

---

Published as a conference paper at ICLR 2025
H PROOF OF THEOREM 3.6
Let the output of Algorithm 2, and denote
= (H / 21)
: := S (SHS" +S
(61)
We have & = q Ei= 800, where E(1) El lid E.
The proof proceeds in two parts. First show that the expectation [[E] is small; then, show that & concentrates around [[ = Elg] using matrix concentration inequalities.
The following lemma, proven in Section I, bounds B[E].
Lemma H.1. Suppose that S has i..dd. Gaussian entries. Assume that m 2 1.5d (A). We have
||[ O
To establish concentration, we use the matrix Bernstein inequality (Tropp et al., 2015, Theorem 6.6.1). We cite it here.
Lemma H.2. Let XI....,X. Xn Raxd be independent Hermitian random matrices. Assume that
E[X,]] = 0, , ll & L, and denote
N
E
E
(62)
Then for all t 0.
-12/2 12 + Lt/3
Pr Σ Xn
21)
S dexp
(63)
We shall apply Lemma H.2 with a truncated version of & El Eeess E[E1 e for suitably chosen L. To this end, we give the following high-probability bound on (| Lemma H.3. Suppose that S has i.i.d. Gaussian entries. For ō E (0, and q 2 ser
L6.4 = - + + log(g/0)
(64)
for large enough C' 0. Then,
1. W.p. 8/g, ||||| < Lö.q.
2. 1| 0 0 dd)
3. || [1 0 1 ed))
We prove Lemma H.3 in Section
The last component needed to apply matrix Bernstein is the variance proxy (62). Note that we always have 2² < I,2N, so that
N
2 < dexpp cN min ( = L (
Pr
X,
(65)
N
However, if the variance E[ < L.2, one may obtain substantially sharper bounds for small < L. This is the case in our setting. The following is proved in Section K.
Lemma H.4. Suppose that S has i...d. Gaussian entries. Assume that m 2 1.5d (A). Then
E| = Odd m).
With the above, ready to conclude the proof of Theorem 3.6.
31

---

Published as a conference paper at ICLR 2025
Proof. (Of Theorem 3.6.)
Consider the truneated matrices 2(1)
LA."" L La. with Ló.g as in Lemma H.3. By Item of
g(1)
Lemma H.3 (and union bound over 1,... D. with probability 5, we have
3
= 1
a(t)
Lo.0
+ E [EL [
(66)
9 t=1
Note that || [Elpens [[[E] + E TEETELL
{ |II [[[[]]] + E [| 0(1//m)
(67)
where the first inequality uses Jensen's inequality, and the second inequality follows from Lemma H.3 Item 2 and Lemma H.1. Similarly,
|| EL.
|||
S | +I|||++ + 1e O(d/m),
(68)
where we used Lemma H.4.
Now, by the matrix Bernstein inequality (Lemma H.2), the bound || E + ((1/mm)
(69)
holds with probability at least
2 8
dexp (-cq min! min
d/m Lp.9 1)
(70)
When q I <c)p(O(d)) and 1/8 S exp(O(d)) we have Ls.q = Od/mm), so for O(1) the smaller
term is 2² And so, q 2 d'oogd d log(1/5) ensures that the probability above is 2 1 - 6.
d/m
2²
1 PROOF OF LEMMA H.1
We use the machinery from Section G. As before, we denote the eigendecomposition of H,
H = V (iiag(T.....,Ta)))) TT),
V = [1] ... |va]
and
11e 11 ToUETT St Sug,
$ S8 l.d. N(O,, mm-
We shall compute the expectation of W = SSHST + +1- in the eigenbasis V, that is, [[VTWV]
Denote, for brevity.
T( = (SHST +X))1 Te((A) (SH, +\11-1,
so that by the Sherman-Morrison formula (Lemma M.1), writing SHST = SH + TEESST
т(ィ) = TTeA) 1+ ++TT
(71)
Let j # f. Note that ', (AA) do not change under a sign flip sy + -j. Since sy has a symmetric
distribution (s, d -5)), we deduce that
E Wo,) B[sTT[iiiii T())s,) 0,
(72)
that is, E[V T WV] is diagonal. Consequently, E EV] = E (H + 12 + 1/3 is diagonal well (since V is eigenbasis of H).
32

---

Published as a conference paper at ICLR 2025
Let us calculate the diagonal elements. By (71),
1/ Wuz = T()ss
8/ Re())8
1
(73)
1 Tes (T)(S)se
(s (T)ss) -1 +
So
A n
(8) (c))ss()
A- (0] (e))))))))))-
(s (TT(((Ass))+++
- -(A8 T (e){ss 11)-
(74)
Let be given by (9), so that by Theorem 3.5, is close to Ā with high probability. Let
A TT. As! T (()sss
1
42 As]
A trT
A3 = tr (T
1
1
T()),
m
m
4 - T())
1
-
tr
1
m
m
m
so that
- As -
We have
E|A
E| T. --1-- =
A
A
1ee11 31|
s 11/
vm
-
EIUTI, - 5 i A1
x 1
/m..
m
1
E S
mm
where we used Lemma D.3 and that Á = 2((2) w.p. 1. As for A2, by Lemma M.8 Item 1,
1 m
1
(BI
1/2 s mm
Thus, we conclude that | = (EV = 0 (mm
J PROOF OF LEMMA H.3
Clearly, ||E|| < 1 1 [|(H + //2 ( By the inequality Va+b s Va + vb (for < 1,
a, 2 0), we have [[
|| < [1 + 1/ 1/ 1/2 + X/ S 2| /2 -
where the second inequality follows from the fact that V is PSD, and therefore for any u, v, (u + w) W) + v) = I|u 01 < 2| 2|| || where ||u|| 2 w : T Wa is norm. We have
1/2 5 (SHST& + Xπ) - 1/2 = 1 SH + I
1/1
122 5 + x
33

---

Published as a conference paper at ICLR 2025
so that || s 2. Moreover, [[W]| |SST + - N I S||. Recalling
that 2 52/12 (by construction), X||||| < 12 |S S|l. Combining all the above, ||E|| 24 |SST
5 +
By (Vershynin, 2018, Theorem 4.6.1), for every 2 0, w.p. 1-2e" -e-
d 1 V
1
||S S|| < 1 C1 max 1²) where
+
m
m
for C1 large enough. Since m & d, note that max 2²} = 1² 2 2 + 2 setting ~ log(1/8).
Item follows by
As for Item 2,
B||sT 5|1 ST 511
Pr(IIS S/| 2 s)ds
1440
1
d
t
201
Pr([IST S|| 2 1 201 202 jt
m
d
m
m
0c
e dt
e
d
Deducing Item 2 of the lemma is now easy from the above estimate.
Finally, for Item
Eijist 51
1ST 111144111
2s Pr([S S|I? 2 $2)ds
1+10 #
(/00 00
(1+20) m
d
t
2C1 m J
+ 2C'
Pr(IIST S|| 2 + 201
+
2C,
dt
m
m
m
1 (
d
1
c d' O d
+ 20
e a
m dd m
m
m²
K PROOF OF LEMMA H.4
We continue from Section I, expressing & in the eigenbasis of H.
Recall that (H 1)) 2)2/2- - I. Then starting from (71), have for any { = I,...,d. d,
1 Ev, = (Tx + x)1/2 + 1/1/2 8, 8yTRe{Ass
-z=j
(75)
1 TesETT
Consequently, j 전 then note that Blu] Ehraa Eus] = 0, since 85, 8h have symmetric distribu-
tion.
Recall, we are interested in bounding E || =
1l We have
(V EV)
EV EV) EV),.).
t=1
which is zero when j * k; that is, E EV)? is diagonal. Let compute the diagonal clements, j = We have
B[(VT (ETEEV B(VV- &V)..) + Σ (VV) TVV).e] (76)
34

---

Published as a conference paper at ICLR 2025
The first term is
( (1101110101000000100000110111111
27
EI(V k = E
O(1As Tie 112),
(77)
where we used (74), Similar to Section I, one can show this is O(1/m). For the second term, using (75),
(S)T)
Σ (VTEV)a((VTEVV VV)ee===
Σ (Te
A) A)
T1 T (A)sa
(=1,#s
8=1,b
d Σ
( + +1)
8k (a(A)
Tr((A)sa.
(78)
(1 + TRSTT (())E))
1=1,
For small > 0, let be the event that sy Ti (A)sx 2 d/A:; A; by similar calculations as in Section I, for some c, C > 0, Pr((2)) & Ce om Then
X²
Σ (VTEV)...(T &e),elln (VT
S
8k RA(A)
Σ (=1,tpk
(Te A)sEs T
T (A)sx11
The
1=1,p
X²
T k R(((A)
Σ
(Te + A)8E88
Ra(()ss
Th +
(=1,EE
so taking the expectation,
d Σ
(VTeV),,e(V)) &V),,e (VVV)alVVVV
S + A m
X Etr Tu(()
d Σ (Te + Nassk
Tacai)
1=1,
1=1,
Tucis)
Etr
Σ (Te + Naselk
T + İ m
Em
Σ
(Te + A)Es T
1=1,6$
Using (71),
8 T Time(A)s
Els T (
Es (A)ss)s EE
(79)
1 TES T (d)ss
1 + T3 57 1: T
where we used Jensen's inequality with the concave function x + 14 x >
By calculation
similar Section I, E[s, ({)s(()) + O ( hence
[ T Tx. (A)sa) Es, Tx S 1 E(s T (A)s()
$
Te + A
Consequently,
2 d
E
Σ
EV)aeln (VV))all
x m
Σ
(re A)Bsy T (()())
(80)
(=1.
t=1,
m
35

---

Published as a conference paper at ICLR 2025
Finally, towards bounding the expectation under the complement event !2" note we write
Σ (V
(Tk + A)
Σ (Te TT
(1 TR (o (2)s8)
1=1,
111,
8/ To
$ (Tk
Σ
(Te + T (
(1 + TASTI (s)
1=1,
8 TT
S (Tx
Σ
(Te + Ti (81)
+ T 1=1,
where the first inequality, (s TT v (s (1) (5 T(()ss))) follows by Cauchy- Schwartz. We have
of (2))s
Tk + x Tk + (5 Re( X (5 Tv Th + (s) T (d)ss)
(Th (Tx + A)
1 Ths
< 1 + [1 T (A)s)) - 11.
Recall 2, the event that Asy TT d E 1 (VTEV), e Under Ω 0
d': we now need to bound
TT
c, and so
(Th- A) ('Ddoo lo 0(1). Plugging this into (81),
d
E
Σ
Σ
(Tx +
(82)
1=1,
E11
Using (79),
A
(Te +
= E
+ (s Ti, ( A)st) 1
(1 |2(s T u, (5)se) - -11 1
EE
< E[(1 + ||s ( 11 S Pr(ff)) + /||ss|||PP((S) S Se
as Pr(8*) S e-m Thus, E (VT (V)..((V (V S om And so, we finally
conclude
E
(V (V).ve [V)) = E Σ
+ E
Σ
EV) VV
EV
EV
L=1,6
8=1,8
t=
d
+ de cm (d/mm).
(83)
m
36

---

Published as a conference paper at ICLR 2025
L PROOF OF THEOREM 4.
Define the Newton decrement at a point Rd:
(()) = (1VC)))))(((())))(
(84)
Denote the approximate Newton decrement, where W(0) is an approximation of the true inverse Hessian:
N(Ø)
1/2
(85)
Note that W is n-accurate in the sense of Definition 2, then
VI - (N - Ñ(O)= - VI+ v
(86)
For self-concordant functions, it is known that the Newton decrement yields upper bound on the suboptimality gap G(Ø) G(
Lemma L.1. (0) < for some numerical constant No S 0.68, then
- -(0)) < (0).
Proof. See (Boyd & Vandenberghe, 2004, Section 9.6.3), specifically Eq. (95).
Clearly, by (86), the approximate Newton decrement (() provides an upper bound on the optimality gap provided that the inverse Hessian is 17-accurate, albeit a looser one.
The key tool in our analysis are the following estimates, due to (Pilanci & Wainwright, 2015).
Lemma L.2. Operate under the conditions of T'eorre 4. I. There exist numerical constants A, 0 with A < 1/16 such that the following holds. Assume that 0 S ll < 1/2.
If at an iteration t, one has N(ft-1) 2 A, then at the next (approximate) Newton step: G(Ot) - GG(")) I -abv, where a, b O are the parameters for backtracking line search (see Algorithm 3).
((0--)) & A then at the next iteration N(O.) & N(O--1).
Proof. See (Pilanci & Wainwright, 2015, Lemma 6.4).
Note that Lemma L.2 implies that within at most
To = (G(0)) - G(8*))/(abv)
(87)
iterations, guaranteed to reach 0, such that N(O.) A; furthermore, have achieved that, N(Ox) < holds for every subsequent iteration 1' 2 f.
The following is the main estimate for the remainder of the analysis: Lemma L.3. Operate under the conditions of Theorem 4.1, and assume that rj < 1/2. Then
N(O.) < (1 + v)) (N((0-1))
(88)
(1 (1 + N(8--1))
Proof. This result is (Pilanci & Wainwright, 2015, Lemma 6.6).
One can further simplify (88) assuming that N(0,-1) < 1/16. Coarsely lower bounding the denom- inator and upper bounding the numerator:
(
(4N2() if N(O 1 (4)N(0--1 N(O,-1) E1
N(O.) < (0.10) 27)N(0,
(89)
provided that N(O,-1) 1) < 1/16.
With the above estimates in hand, ready prove Theorem 4.1.
37

---

Published as a conference paper at ICLR 2025
Proof. (Of Theorem 4.1.) Recall that assume that 7 1/5,, and in particular 1 1/2. By Lemma L.2, within To iterations, the approximate Newton method reaches 0, such that N(O,.) < 1/16 for all t' 2 t. In particular, also G(Ø) - G( (1 + n)1/162 I 3/500. Let 0 be the desired precision; suppose that < 3/500. We consider two cases.
First, assume that 1) 3/500. Let (e) be the smallest integer t such that the dynamic
A: 442 Ao 1/16
(90)
satisfies A, < E. It is easy to verify that A, 1/4(1/4)2 (e.g. by induction), hence Ti Ollog log(1/e) /e)). We conclude that when & 2 7/. T(e) To + (e).
The second case is when E 7), therefore are in the regime of (89) where decay is linear rather than quadratic. Let (2) be the smallest integer t such that the dynamic
B. B1 Bo = "
(91)
satisfies B, S e. Clearly B, = (4n)) Bo, hence (2 0 log(n/s) 1/n) (here we used that 4nj < 1, by
assumption). We conclude that when 17, T(E) S To (n) + T2(e).
38

---

Published as a conference paper at ICLR 2025
M AUXILLIARY TECHNICAL RESULTS
Lemma M.1 (Sherman-Morison) Let A E Ruxn be invertible, and E IR.. The matrix uu is invertible if and only if T A o O. In that case,
A A- uv' T A-1
(A + uv ) A-1
(92)
1+ 0 A-1u
Lemma M.2 (Block matrix inversion). Let A, B, C, D be matrices of conforming dimensions. Pro- vided that A, D are both invertible,
A B (A
CD-
0 (D
-BD-1
C
D
0
(D-CA-1B)-1
CA-1
(93)
Lemma M.3. (Resolvent low-rank perturbation.) Let A, ž 0 be PSD matrices, and x > 0.
Then
tr -r B) rank(A - B) |||
2
For reference, see for example (Bai & Silverstein, 2010).
Let {0.2) =: Fo < S < be filtration of g-algebras over some given probability space Q. Recall that random process {X.,] 1 is martingale difference (adapted to this filtration) if 1) for every n, Xn, is Fx -=-mmasurabl;; and 2) Fn-1-almost surely, IE[ Fn-1l 0.
Lemma M.4 (Azuma-Houe Hoeff--oe Suppose that (x),Fr) is martingale difference, and (Br) JN constants. Suppose that for all S n N, almost surely, |Xal S Bn- Then for all
t $ 0.
1
Pr( Σ 5al20)2ee
2
N
B
1==1
The sub-Gaussian and sub-Exponential norms of a (scalar) random variable are, respectively
[ = inf frfo>0 Ee /0 S $2
| = inf {o fo>00 >0
E < 52]
(94)
We call X sub- Gaussiausin resp. sub- Exponential, if || l oc, resp. ||
For random vector define ||rlls sup / rrll, That is, it is the largest ev-morm of one-dimensional projection of r.
Lemma M.5. The following holds.
1. || Ilve is a norm on the subspace of random variables X such that ||X|| < co.
2. If X,Y Y are sub- Gaussian (possible statistically dependent), then XY is sub- Subb-Exonential, with |XYYl, S ||X ll |Y Ilwa-
3. If ( < B almost surely for B 2 O constant, ||XY I BB|
4. (Centralization.) For any X, |X [[X]|I, S ||
We refer to the book (Vershynin, 2018) for reference.
Lemma M.6 (Bernstein). Suppose that (X,,, Fn) is a martingale difference. Then for all t. 2 0,
Pr( X 2) < 2exp (-eminn max,=1..... 1x
t
t
EN || 112
39

---

Published as a conference paper at ICLR 2025
Next, cite the following version of the Hanson- Wright inequality, see (Vershynin, 2018). (Hanson-Wright). Suppose that has independent Sub-Gaussian entries, with E[r;] 1, max; ||ri K. For any matrix and O,
Lemma M.7 (Hanson- Wright),
B
Pr(I T Ar Ar ir(A)) 2 1) S exp (-emint K2IA K IA
(95)
where e > is a universal constant.
Lastly, the following lemma collects, for convenience, some properties satisfied by any isotropic random vector that satisfies concentration inequality of the form (95).
Lemma M.8. Suppose r is istropic E[r] = O, 1.rrr = I, and satisfies (95) for some K > 0.
1. For any matrix (independent rf B[[r" T Ar - tr(A) IA
2. For any vector v (independent r), -a1 rlls, S K||l|l- That is, is Sub-Gaussian with ||lle 2 K.
Proof. Start with the first item. Write B|r T Ar - tr(A)l = S Pr(p Ar tr(A))2 > s)ds
2 2.0 tPr(I Ar tr(A) t)dt. Using (95), crudely bound
tPr(I Ar - tx()))² t)dt & [0/ texp(- (-c A*
t²
dt + /0 texpl -c Ka
)t S SK
As for the second item, note that ||w" = aaal ||5 = ||/||2 Hence,
Pr([o Pr(loT/1 2 Pr(rT vo T 2 12) S exp( min K2||0||? K( 1)-
That is, 11 (uTr) T|] has a Gaussian tail with variance proxy K21|0|12.
40


Reviews:
[{"review_id": "d8c9oSEDeJ", "reviewer": "Reviewer_K6ze", "summary": "The paper presents a new optimization method for distributed centralized optimization, where all workers have access to the full data. The method is based on damped Newton method and inexact approximation of inverse Hessian with sketching by Marchenko-Pastur scheme. The authors present extensive theoretical support of their method with a focus on debiasing properties of the proposed approximation.", "questions": "1) Could you clarify the distributed setting paradigm? As I understand, all workers have access to the whole function and they have the same gradient. As this setting is quite unusual to me, especially from Federated Learning and privacy perspective.\n\n2) Does this method and approximation work for pure convex functions without regularisation? Also, may they work for non-convex functions in practice? \n\n3) How good is the proposed approximation compared to Krylov and Nystorm-based approximations?\n\n4) It might be interesting to apply this approximation to the new branch of methods based on cubic regularisation:\n\n [6] Scieur, Damien. \"Adaptive Quasi-Newton and Anderson acceleration framework with explicit global (accelerated) convergence rates.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2024 \n\n[7] Kamzolov, Dmitry, et al. \"Cubic Regularization is the Key! The First Accelerated Quasi-Newton Method with a Global Convergence Rate of $ O (k^{-2}) $ for Convex Functions.\" arXiv preprint arXiv:2302.04987 (2023)\n\n[8] Jiang, Ruichen, et al. \"Krylov Cubic Regularized Newton: A Subspace Second-Order Method with Dimension-Free Convergence Rate.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2024.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "I think that the results are new. I appreciate the proposed approximation and guarantees by itself without the need of a distributed setting. This approach is fresh for the optimization community. The paper is complicated but well-written.", "weaknesses": "1. One of the main weaknesses from my point of view is the experiments section. \na) Firstly, there is no comparison with other methods. So, it is unclear how fast is the proposed method. For example, it can be interesting how far is it from the exact Newton method. \nb) Secondly, the Ridge Regression seems to be a very simple task for Newton-based methods, as it is a quadratic function with a relatively small noise($\\sigma = 0.01$) and big regulariser $\\lambda = 0.001$.  So, it might be much more useful to show how good is the inverse approximation depending on the conditioning of the Hessian and size of $m$, similar to Fig.1. \n\n2. Small remark: I recommend changing the indexing in formulas 2, 3, 4. As $g_t$ is a gradient in point $\\theta_{t-1}$ and should be $g_{t-1}$, also Hessian $H_t$ is calculated in $\\theta_{t-1}$.\n\n3. Also, the authors may improve the related literature section by some distributed Newton and Quasi-Newton methods. Some possible references:\n\n [1] Safaryan, M. et.al. \"FedNL: Making Newton-type methods applicable to federated learning.\" (2021) arXiv preprint arXiv:2106.02969. \n\n[2] Agafonov, A. et.al. \"Flecs: A federated learning second-order framework via compression and sketching\". (2022) arXiv preprint arXiv:2206.02009. \n\n[3] Qian, X. et.al  \"Basis matters: better communication-efficient second order methods for federated learning.\" (2021) arXiv preprint arXiv:2111.01847. \n\n[4] Elgabli, A. et.al. \"FedNew: A communication-efficient and privacy-preserving Newton-type method for federated learning\". (2022) In International conference on machine learning (pp. 5861-5877). PMLR. \n\n[5] Bischoff, S. et.al. \"On second-order optimization methods for federated learning.\" (2021) arXiv preprint arXiv:2109.02388.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "cg4kSkuhaE", "reviewer": "Reviewer_wj6U", "summary": "The key idea of the paper is to relate random matrices and one of the key optimization methods, Newton's method. Convergence for strongly convex self condordant functions is proved.", "questions": "See 3 in Weaknesses", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1) A new (best of my knowledge) matrix sketching scheme is presented.\n2) The proofs look correct.", "weaknesses": "1) Convergence only for self condordant functions. The result would be more solid if the authors consider general strongly convex functions.\n\n2) There is no comparison with other methods in the experiments. There is a huge family of Newton-like methods where the Hessian is approximated. Starting from quasi-Newton technique to newer Hutchinson [1] or Fisher approximations [2]. The same Adam can be considered as a method with a diagonal approximation of the Hessian. I could also mention Shampoo or Soap or [5]\n\n[1] Jahani, M., Rusakov, S., Shi, Z., Richt´arik, P., Mahoney, M. W., and Tak´aˇc, M. (2021). Doubly adaptive scaled algorithm for machine learning using second-order information.\n\n[2] Amari, S. (2021). Information geometry.\n\n[3] Gupta, V., Koren, T., and Singer, Y. (2018). Shampoo: Preconditioned stochastic tensor optimization.\n\n[4] Vyas, N., Morwani, D., Zhao, R., Shapira, I., Brandfonbrener, D., Janson, L., and Kakade, S. (2024). Soap: Improving and stabilizing shampoo using adam\n\n[5] Doikov, N., Jaggi, M., et al. (2023). Second-order optimization with lazy hessians. \n\n3) There's an important point in the proof of Theorem 3.6: $\\delta$ is small probability, what do we need to do if $d \\cdot \\log\\delta^{-1} \\cdot e^{-2}$ is greater than $\\delta \\cdot e^d$? Maybe I misunderstood something.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "9g5eQvbWRH", "reviewer": "Reviewer_UJF8", "summary": "This paper proposes a novel sketching Newton method. This paper also provides the convergence analysis of the algorithm and gives the non asymptotic guarantees.", "questions": "No", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper proposes a novel sketching Newton method. This paper also provides the convergence analysis of the algorithm and gives the non asymptotic guarantees. The massively parallel fashion is consider in their method.\nThe theory in this paper is good.", "weaknesses": "This paper claims that the proposed method is motived by serverless cloud computing, in particular the ``function as a service'' (FaaS) model. However, Eq. (7) is an biased estimation of $W_t$. Thus,  even $m$ machines run in parallel. This will not help to obtaim an $m$-times speed-up.\nThe convergence rate of the proposed algorithm with $m$ machines running in parallel is the same to one only runing in  a single machine.\n\nThe experiments prove the biases of the proposed algorithm is lower than other agorithm and the algorithm can achieve lower iteration complexity.\nHowever, the coresponding computation cost is not mentioned. Thus, the comparisons in this paper seem not enough to prove the advantages of the proposed algorithm.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}, {"review_id": "QzYiihaRyW", "reviewer": "Reviewer_1hmD", "summary": "This paper propoed a distributed Newton method, where the server performs the Newton update,  and the workers use Hessian skeching to return a neatly unbiased inverse Hessian to the server. The solution is based on random matrix theory, i.e. the Marchenko-Pastur law.", "questions": "1. The paper says \"The connection to the Marchenko-Pastur law treats in a uniform  framework—and recovers almost immediately—earlier results on Hessian sketch debiasing for distributed Newton methods (Derezinskietal.,2020;Zhang&Pilanci,2023). Can the authors compare the works of  (Derezinskietal.,2020;Zhang&Pilanci,2023) and this work for me to better clarify the contributions of this paper. \n\n2. Why the authors claim their result on non-asymptotic error guarantees in Section 3.1 is dimension-free. I can not uderstand this because Theorem 3.6 requires $q = \\Omega (d / (m \\epsilon^2)) $. I can understand that if we have enough workers, i.e, $q = \\Omega (d / (m \\epsilon^2)) $ then we can set $m = O(1/\\epsilon^2)$, which is dimension-free. But in applications, I think $q$ is considered to be fixed and then to ensure the error is smaller than $\\epsilon$, one should take $m = O( d / (q \\epsilon^2 )$ if $q <d$. In this view, the gurantee seems to depend on dimension $d$.\n\n3. (continued my second question) What is the total complexity to minimize a self-concordant function? By Theorem4.1 if we have $\\eta < \\epsilon$, then the algorithm can recovers the $O(\\log \\log (1/\\epsilon))$ rate of Newton's method. By according to Theorem 3.6, $m$ should be $m = O(  1/ \\epsilon^2 \\lor d / (q \\epsilon^2 ))$, which may be very large. The advantage of Newton's method is that it can have superlinear convergence rate, but as $m = \\Omega(1/\\epsilon^2)$, the total complexity has polynomial dependency on $1/\\epsilon$. Then the total complexity results seems to be sublinear.  If so, what is the advatange of Newton's method compared to gradient methods? In other words, can the authors compare the total complexity of their Newton's method and distributed gradient methods (including the accleration) for me?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The solution in this paper is novel to me. I think the idea is very interesting.", "weaknesses": "See the questions part.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)