Newton Meets Marchenko-Pastur: Massively Parallel Second-Order Optimization with Hessian Sketching and Debiasing

## Abstract

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 are essentially dimension-free. 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.

### 1 Introduction

Consider minimizing a convex, twice-differentiable function $F:\mathbb{R}^{d}\to\mathbb{R}$,
subject to an L2 regularization penalty:

$$ $\min_{\theta\in\mathbb{R}^{d}}F(\theta)+\frac{\lambda}{2}\|\theta\|^{2}\,,$ (1) $$

where $\lambda>0$.
Such problems frequently appear in machine learning and statistics. Common examples include ridge regression, logistic regression with regularization, support vector machines (SVMs) and kernel machines, among others \citephastie2009elements.

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) operates with limited individual computational resources; and (iii) cannot communicate directly with other workers,
only with a 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 a service" (FaaS) model \citepjonas2017occupy. Serverless computing and FaaS are well-suited for computing gradients and Hessians on large datasets, enabling parallelization of these computations across a large number of workers. Additionally, this approach provides resilience to failed or straggler workers, ensuring robust and efficient processing \citepbartan2019straggler.
In designing a solution of this setting, we are essentially guided by two key principles:
(a) minimizing the number of communication or worker deployment rounds, with each round utilizing a 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 a
scheme for solving
(1) in a massively parallel manner.
In our scheme, the server aims to solve (1) using Newton’s method, a second-order iterative algorithm which is known to 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 $O(d^{3})$ flops; in high-dimensional problems (large $d$), this may be prohibitively expensive.

We wish to offset the high per-iteration cost of Newton’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 a 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:

$$ $\displaystyle\theta_{t}=\theta_{t-1}-\alpha_{t}W_{t}g_{t},\qquad(t=1,\ldots,T)\,,$ (2) $$

where $\alpha_{t}$ is a step size (typically chosen by line search) and the gradient and inverse Hessian are

$$ $\displaystyle g_{t}$ $\displaystyle:=\nabla F(\theta_{t-1})+\lambda\theta_{t-1},$ (3) $\displaystyle W_{t}$ $\displaystyle:=(H_{t}+\lambda I)^{-1},\;H_{t}:=\nabla^{2}F(\theta_{t-1})\;.$ (4) $$

For sufficiently smooth convex objectives $F$, Newton’s method is known to converge at a quadratic rate: $T=O(\log\log(1/\varepsilon))$ iterations are sufficient to approach the minimum within error $\varepsilon$, as opposed to $T=O(\log(1/\varepsilon))$ for gradient descent and its accelerated variants; see [boyd2004convex, Chapter 9]. A fast convergence rate is highly desirable in our setting, as the number of iterations corresponds directly to communication/worker deployment rounds.

In our scheme, the server effectively runs the Newton method with an approximate inverse Hessian:

$$ $\theta_{t}=\theta_{t-1}-\alpha_{t}\bar{W}_{t}g_{t},$ (5) $$

where $\bar{W}$ is obtained by averaging $q$ independent local estimates produced by the workers:(^1^11In practice, the workers send the vector $\hat{W}^{(k)}g_{t}\in\mathbb{R}^{d}$ to the server, rather the whole matrix $\hat{W}^{(k)}\in\mathbb{R}^{d\times d}$.)

$$ $\bar{W}:=\frac{1}{q}\sum_{k=1}^{q}\hat{W}_{t}^{(k)}.$ (6) $$

##### Debiased inverse Hessian sketching:

We wish to avoid directly inverting the full $d\times d$ Hessian, which practically (when done exactly) costs $O(d^{3})$. 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 a random sketching matrix $S\in\mathbb{R}^{m\times d}$ (where $m<d$), and approximates the inverse Hessian by its sketched version

$$ $\hat{W}=S^{\top}(SHS^{\top}+\tilde{\lambda}I)^{-1}S,$ (7) $$

where $\tilde{\lambda}>0$ is a *modified* regularization parameter. Sketches of the type in (7) have been studied before in the literature (also in the context of high-dimensional regression and optimization). Variations have appeared under various names: “dual sketching” \citepzhang2013recovering, “sketch-and-project” \citepgower2015randomized, “right-sketch” \citepmurray2023randomized, and “feature sketching” \citeppatil2023asymptotically.
Note that the cost of forming $\hat{W}g$ is $O(m^{3}+M)$, where $M$ is the cost of multiplying $SHS^{\top}$.
In this paper, we exclusively consider the case where $S$ is a random dense matrix with i.i.d. entries.(^2^22By “dense”, we mean that $S$ has $\Omega(md)$ many non-zero entries typically, though it’s not necessary that all the entries, or even a majority of them, are non-zero.) In this case $M=d^{2}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, $O(d^{2}m)$, is substantially smaller than $O(d^{3})$.

###### Remark 1 .

In this paper we do not assume any particular form for $H$ (i.e., that it is readily factorizable $H=X^{\top}X$, $X\in\mathbb{R}^{N\times d}$, a setting which is quite common in the sketching and distributed optimization literature; for example \citeppilanci2015newton). In fact, sketched measurements $SHS^{\top}$ 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 \citeppaszke2017automatic.

A key question is how to choose $m$ (and in accordance $\tilde{\lambda}$). Define the $\lambda$- effective dimension of $H$,(^3^33$\mathsf{d}_{H}(\lambda)$ is also called the “effective degrees of freedom” in ridge regression \citephastie2009elements, when the design matrix $X\in\mathbb{R}^{N\times d}$ is $H=X^{\top}X$.)

$$ $\mathsf{d}_{H}(\lambda):=\mathop{\mathrm{tr}}(H(H+\lambda I)^{-1}).$ (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>\mathsf{d}_{H}(\lambda)$, the choice

$$ $\tilde{\lambda}:=\lambda\left(1-\frac{1}{m}\mathsf{d}_{H}(\lambda)\right)$ (9) $$

results in $\hat{W}$ which is a *low-bias* estimator of the true inverse Hessian $W$. Note that exactly computing 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 a novel low-complexity procedure, based on ideas from asymptotic random matrix theory, for adaptively selecting $m=O(\mathsf{d}_{H}(\lambda))$ and $\tilde{\lambda}$, by only observing sketched measurements of $H$.

###### Remark 2 .

When $H$ is effectively low-rank, that is exhibits fast spectral decay, $\mathsf{d}_{H}(\lambda)$ may be much smaller than $d$. For example, when the spectrum has power decay $\lambda_{k}(H)\propto k^{-\alpha}$ and $\lambda=O(1)$, $\mathsf{d}_{H}(\lambda)=O(d^{1-\alpha})$ when $\alpha<1$, $\mathsf{d}_{H}(\lambda)=O(\log d)$ when $\alpha=1$ and $\mathsf{d}_{H}(\lambda)=O(1)$ when $\alpha>1$. Hessian spectral with power decay are common in regression problems involving kernels \citepwainwright2019high, yang2017randomized.

###### Remark 1 .

###### Remark 2 .

##### Outline and Contributions:

- •
In Section 2 we review fundamental results from asymptotic random matrix theory (RMT), specifically the Marchenko-Pastur law and its formulation using deterministic equivalents, and demonstrate 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 recovers almost immediately—earlier results on Hessian sketch debiasing for distributed Newton methods
\citepderezinski2020debiasing,zhang2023optimal.
- •
In Section 3 we propose a *novel data-driven method* for selecting $m$ and $\tilde{\lambda}$ using only the sketched Hessian, avoiding the need to compute the effective dimension exactly. The approach is a natural outgrowth of the RMT framework described previously, and is fully adaptive: no tuning parameters are needed.
Section 3.1 further provides non-asymptotic error guarantees for our method, focusing on the case where $S$ is an i.i.d. Gaussian matrix. Notably, the error bounds obtained in this setting are essentially dimension-free and independent of the Hessian’s condition number.
- •
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 $\eta$-accurate estimate $\bar{W}$. By ensuring a sufficiently small $\eta$, the algorithm achieves an arbitrarily fast linear convergence rate. Our results are proved for self-concordant functions $F(\cdot)$, 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-order optimization methods have been widely studied for large-scale machine learning. Notable examples include DANE \citepshamir2013communication, which addresses communication efficiency, AIDE \citepreddi2016aide, designed for accelerated convergence, and DiSCO \citepzhang2015disco, which focuses on distributed second-order methods for convex optimization problems. 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 \citepwang2017giant uses the average of local Newton steps as the global step and Determinantal Averaging \citepderezinski2020debiasing provides an unbiased averaging technique for distributed Newton methods. Surrogate sketching \citepderezinski2020debiasing was introduced for distributed Newton’s method, revealing improvements via a simple shrinkage strategy.
The work \citepzhang2023optimal introduced an optimal debiasing method for distributed Newton’s method when the effective dimension of the Hessian is known.
Lastly, there exists a large related corpus of works in optimization involving Hessian subsampling or sketching, for cases where the Hessian admits a good low-rank approximation, though not necessarily in a distributed setting. Examples include \citeproosta2019sub,
gower2019rsnrandomizedsubspacenewton,frangella2022sketchysgd, among others.

### 2 Debiasing from I.I.D Sketches and Random Matrix Theory

This section concerns inverse sketched Hessian debiasing for sketches $S\in\mathbb{R}^{m\times d}$ which have i.i.d. entries with mean zero $\mathbb{E}[S_{ij}]=0$. We normalize the variance as $\mathbb{E}[S_{ij}]^{2}=1/m$; this normalization is such that $\mathbb{E}[S^{\top}S]=I$, equivalently, $S:\mathbb{R}^{d}\to\mathbb{R}^{m}$ is an isometry in expectation.

Let $H\succeq 0$ be a positive definite matrix (e.g., the full Hessian at the current iteration). The spectral properties of the matrix $SHS^{\top}\in\mathbb{R}^{m\times m}$, which appears in the sketch (7), are of utmost importance for what follows. Note that this matrix has the same non-zero eigenvalues(^4^44This is a standard consequence of the elementary determinant identity $\det(I+AB)=\det(I+BA)$.) as the matrix $H^{1/2}S^{\top}SH^{1/2}$: the latter is a *sample covariance* matrix (making the analogy explicit: $m$ is the number of samples, $H$ is the population covariance, and $x_{i}:=\sqrt{m}H^{1/2}r_{i}$ are independent “samples”, $r_{1},\ldots,r_{m}$ being the rows of $S$); accordingly, $SHS^{\top}$ is sometimes called the companion sample covariance.
Sample covariance matrices, and their spectral properties, have been a central object of investigation in random matrix theory (RMT); cf. \citepbai2010spectral. A foundational result is the Marchenko-Pastur law \citepmarvcenko1967distribution,
which describes the high-dimensional global behavior of the eigenvalues of $H^{1/2}S^{\top}SH^{1/2}$ (limiting spectral distribution).

##### The Marchenko-Pastur equation:

For(^5^55Here $\mathbb{C}^{+}$ is the complex upper halfplane, $\Im(z)>0$.) $z\in\mathbb{C}^{+}$ let ${\mathsf{s}}(z)\in\mathbb{C}^{+}$ be the unique solution of

$$ $\frac{1}{{\mathsf{s}}(z)}=-z+\frac{1}{m}\mathop{\mathrm{tr}}\left(H(I+{\mathsf {s}}(z)H)^{-1}\right)$ (10) $$

subject to the constraint ${\mathsf{s}}(z)\in\mathbb{C}^{+}$. It is known that (see \citepbai2010spectral) ${\mathsf{s}}(z)$ is the Stieltjes transform(^6^66For a finite measure $\mu$ on $\mathbb{R}$, its Stieltjes transform is the function ${\mathsf{s}}_{\mu}:\mathbb{C}^{+}\to\mathbb{C}^{+}$ given by ${\mathsf{s}}_{\mu}(z)=\int\frac{1}{x-z}d\mu(x)$. Using the Stieltjes inversion formula, the measure $\mu$ can be recovered from ${\mathsf{s}}_{\mu}$ in the following sense: for every $\varphi(\cdot)$ bounded and continuous, $\int\varphi(t)d\mu(t)=\lim_{\eta\to 0+}\frac{1}{\pi}\int\varphi(x)\Im{\mathsf{ s}}_{\mu}(x+\mathsf{i}\eta)dx$.) of a *compactly supported probability measure*, whose support lies in $\mathbb{R}_{\geq 0}$; this measure is the (companion) Marchenko-Pastur law associated with the population covariance $H$ (and depends only on its eigenvalues) and the number of samples $m$. Moreover: 1) The function $z\mapsto{\mathsf{s}}(z)$ extends continuously to $\mathbb{R}$, except possibly at $z=0$; 2) for $x\in\mathbb{R}$ outside the support, the function $x\mapsto{\mathsf{s}}(x)$ is increasing; in particular, for negative $x<0$, it is positive and increasing. In this paper, we shall restrict our attention to only negative real valued $z$ in (10), understanding that it holds also for $z<0$ real. For more details, cf. \citepsilverstein1995analysis,bai2010spectral.
Finally, note that the function $z\mapsto{\mathsf{s}}(z)$ is invertible, the inverse function given by

$$ $\Psi({\mathsf{s}})=\frac{1}{m}\mathop{\mathrm{tr}}(H(I+{\mathsf{s}}H)^{-1})- \frac{1}{{\mathsf{s}}}.$ (11) $$

The Marchenko-Pastur equation (10) may equivalently be written as $z=\Psi({\mathsf{s}}(z))$.

In RMT, it is typical to consider the *proportional growth* asymptotic regime: $m,d\to\infty$, $m=\Theta(d)$.

###### Definition 1 .

For (small) constant $\tau>0$, assumption set [Asymp($\tau,H,\xi$)]  holds if: (1) bounded aspect ratio (proportional growth): $\tau<\frac{m}{d}<1/\tau$; (2) bounded operator norm: $\|H\|<1/\tau$; (3) $H$ is invertible; (4) $H$ has non-degenerate spectrum: at least $\tau d$ eigenvalues are larger than $\tau$.(5) $\sqrt{m}S_{ij}\sim\xi$, with their law $\xi$, satisfying: $\mathbb{E}[\xi]=0$, $\mathbb{E}[\xi^{2}]=1$, and $\mathbb{E}[|\xi|^{k}]<\infty$ for all $k>0$.

Denote by ${\hat{\mathsf{s}}}(z)$ the Stieltjes transform of the empirical eigenvalue distribution of $SHS^{\top}$:

$$ ${\hat{\mathsf{s}}}(z)=\frac{1}{m}\mathop{\mathrm{tr}}(SHS^{\top}-zI)^{-1}= \frac{1}{m}\sum_{i=1}^{m}\frac{1}{\lambda_{i}(SHS^{\top})-z}.$ (12) $$

Note that ${\hat{\mathsf{s}}}(z)$ is a random function, since $S$ is random. The following is a (quantitative) statement of the Marchenko-Pastur law \citepknowles2017anisotropic; for simplicity, we focus on the case $z<0$.

###### Theorem 1 .

Set $\eta,D>0$. Assume [Asymp($\tau,H,\xi$)]. W.p. $1-O(m^{-D})$, simultaneously for all $-1/\tau\leq z\leq-\tau$,

$$ $|{\hat{\mathsf{s}}}(z)-{\mathsf{s}}(z)|=O\left(m^{-1/2+\eta}\right).$ (13) $$

The constants in the $O(\cdot)$ notation depend on $\tau,\eta,D$.

###### Proof.

See Theorem 3.7 and Remark 3.10 in \citepknowles2017anisotropic.
∎

###### Definition 1 .

###### Theorem 1 .

###### Proof.

##### Deterministic equivalent for the sample covariance:

The Marchenko-Pastur law (Theorem 1) describes the eigenvalue distribution of $SHS^{\top}$, equivalently $H^{1/2}S^{\top}SH^{1/2}$, and relates it to that of $H$ (the connection made via the Marchenko-Pastur equation (10)). One can show something stronger: the entire (random) resolvent $(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}$ is in fact close (entrywise, *not* in operator norm) to the resolvent of $H$ (which is deterministic), with appropriate shifting and shrinkage.
A result of the following kind is sometimes referred to in the RMT literature as a “deterministic equivalent”, cf. \citepcouillet2011random,bai2010spectral.

###### Theorem 2 .

Set $\eta,D>0$. Assume [Asymp($\tau,H,\xi$)]. For any $u,v\in\mathbb{R}^{d}$, $\|u\|,\|v\|=1$, the following holds. W.p. $1-O(m^{-D})$, simultaneously for all $-1/\tau^{-1}\leq z\leq-\tau$,

$$ $u^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}v=u^{\top}(-z{\mathsf{s}}(z)H-zI)^{-1 }v+O(m^{-1/2+\eta}).$ (14) $$

Above, the constants in the $O(\cdot)$ notation depend on $\tau,\eta,D$ and are uniform in $u,v$. In particular,

$$ $\left\|\mathbb{E}\left[(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}\right]-(-z{\mathsf{s} }(z)H-zI)^{-1}\right\|=O(m^{-1/2+\eta})$ (15) $$

###### Proof.

See Eq. (3.3) and Theorem 3.7 in \citepknowles2017anisotropic.
∎

###### Theorem 2 .

###### Proof.

##### Deterministic equivalent for the sketched Hessian:

Recently, a deterministic equivalent for the sketched matrix (7) was obtained by \citeplejeune2022asymptotics. We cite a quantitative version of their result \citep[Theorem 4.1]lejeune2022asymptotics.

###### Theorem 3 .

Set $\eta,D>0$. Assume [Asymp($\tau,H,\xi$)], and furthermore $\lambda_{\mathop{\mathrm{rank}}(H)}(H)\geq\tau$. For any $u,v\in\mathbb{R}^{d}$, $\|u\|,\|v\|=1$, w.p. $1-O(m^{-D})$, simultaneously for all $-1/\tau\leq z\leq-\tau$,

$$ $u^{\top}S^{\top}(SHS^{\top}-zI)^{-1}Sv=u^{\top}(H-({\mathsf{s}}(z))^{-1}I)^{-1 }v+O(m^{-1/2+\eta}).$ (16) $$

In particular,

$$ $\left\|\mathbb{E}\left[S^{\top}(SHS^{\top}-zI)^{-1}S\right]-(H-({\mathsf{s}}(z ))^{-1}I)^{-1}\right\|=O(m^{-1/2+\eta}).$ (17) $$

###### Proof.

The result is a consequence of Theorem 2; for completeness, let us repeat the argument of \citeplejeune2022asymptotics, assuming $H$ is full rank (see their paper for a generalization). Denote $\tilde{u}:=H^{-1/2}u$, $\tilde{v}:=H^{-1/2}v$, so that $\|\tilde{u}\|,\|\tilde{v}\|\leq 1/\sqrt{\tau}=O(1)$. We have w.p. $1-O(m^{-D})$,

$$ $\displaystyle u^{\top}S^{\top}(SHS^{\top}-zI)^{-1}Sv=\tilde{u}^{\top}H^{1/2}S^ {\top}(SHS^{\top}-zI)^{-1}SH^{1/2}\tilde{v}$ $\displaystyle\quad=\tilde{u}^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}H^{1/2}S^{ \top}SH^{1/2}\tilde{v}=\tilde{u}^{\top}\tilde{v}+z\tilde{u}^{\top}(H^{1/2}S^{ \top}SH^{1/2}-zI)^{-1}\tilde{v}$ $\displaystyle\quad=\tilde{u}^{\top}\tilde{v}-\tilde{u}^{\top}({\mathsf{s}}(z)H +I)^{-1}\tilde{v}+O(m^{-1/2+\eta})=u^{\top}(H+({\mathsf{s}}(z))^{-1}I)^{-1}v+O (m^{-1/2+\eta}).$ $$

∎

###### Theorem 3 .

###### Proof.

### 3 An Adaptive Sketching and Debiasing Procedure

Theorem 3 suggests a clear path for obtaining a low-bias estimate of the inverse Hessian from the sketch (7), namely: one should choose $\tilde{\lambda}>0$ such that ${\mathsf{s}}(-\tilde{\lambda})=1/\lambda$, *when such a root exists*.

###### Lemma 1 .

A solution $\tilde{\lambda}>0$ satisfying ${\mathsf{s}}(-\tilde{\lambda})=1/\lambda$ exists if and only if $m>\mathsf{d}_{H}(\lambda)$. When this is the case, $\tilde{\lambda}$ is given by (9).

###### Proof.

By merit of being the Stieltjes transform of a finite measure supported on $\mathbb{R}_{\geq 0}$, ${\mathsf{s}}(\cdot)$ maps bijectively $(-\infty,0)$ to $(0,{\mathsf{s}}(0))$, hence a solution $\tilde{\lambda}>0$ exists if and only if ${\mathsf{s}}(0)>1/\lambda$. Taking the limit $z\to 0+$ in (10), we find (after routine algebraic manipulation) that ${\mathsf{s}}(0)$ solves $\mathsf{d}_{H}(1/{\mathsf{s}}(0))=m$. Since the function $\mu\mapsto\mathsf{d}_{H}(\mu)$ is decreasing for positive $\mu$, we deduce that $1/{\mathsf{s}}(0)<\lambda$ if and only if $m>\mathsf{d}_{H}(\lambda)$. When this is the case, $\tilde{\lambda}=\Psi(1/\lambda)$ where $\Psi$ is the explicit inverse (11).
∎

###### Remark 3 .

Earlier works on distributed Newton’s method \citepderezinski2020debiasing,zhang2023optimal, considered debiasing for a closely related but different Hessian sketch. Specifically, they assume the Hessian has the form $H=X^{\top}X$ where $X\in\mathbb{R}^{N\times d}$; such Hessians naturally appear in machine learning problems,
for example in the (reguarlized) empirical loss minimization for a generalized linear model (GLM).
\citepzhang2023optimal considered (among others), the sketch $X^{\top}S^{\top}SX$, and calculated the appropriate shrinkage factor for asymptotically debiasing the inverse: $(\gamma X^{\top}S^{\top}SX+\lambda I)^{-1}\simeq(H+\lambda I)^{-1}$ where $\gamma:=1/\left(1-\frac{1}{m}\mathsf{d}_{H}(\lambda)\right)$, provided that $m>\mathsf{d}_{H}(\lambda)$, under an asymptotic setting as in Theorem 2 (though without making the connection to the Marchenko-Pastur law explicit). In the scheme they considered, the effective dimension was assumed to be known. We remark that Theorem 2 holds verbatim if $H^{1/2}$ is replaced by $X^{\top}$, and the ratio $\tau<d/N<1/\tau$ is bounded \citepknowles2017anisotropic. Rewritting (15) as $(H+({\mathsf{s}}(z))^{-1}I)^{-1}\simeq\left((-z{\mathsf{s}}(z))^{-1}X^{\top}S^
{\top}SX+({\mathsf{s}}(z))^{-1}I\right)^{-1}$, setting $z=-\tilde{\lambda}$ recovers their proposed shrinkage factor. Accordingly, the adaptive procedure we propose below applies to that sketching scheme as well, with the obvious modifications.

Note that we do not have access to the function ${\mathsf{s}}(\cdot)$ directly (to calculate it numerically via (10), one needs to invert the full Hessian)—but we *can* estimate it from the sketched measurements $SHS^{\top}$, using Theorem 1. Accordingly, we can estimate $\tilde{\lambda}$ by the solution $\hat{\lambda}$ of the (random) equation

$$ ${\hat{\mathsf{s}}}(-\hat{\lambda})=1/\lambda,\quad\hat{\lambda}>0,$ (18) $$

if it exists. To have $\hat{\lambda}\approx\tilde{\lambda}$, we need that $m>\mathsf{d}_{H}(\lambda)$, 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 1-3, one needs to work with $z$ bounded away from 0 0. As $\tilde{\lambda}=\lambda(1-\frac{1}{m}\mathsf{d}_{H}(\lambda))$, to ensure that this is $\Omega(\lambda)$ we need to have $m\geq(1+\Omega(1))\mathsf{d}_{H}(\lambda)$. Specifically, we shall aim to choose $m$ between $1.5\mathsf{d}_{H}(\lambda)\leq m\leq 4\mathsf{d}_{H}(\lambda)$, the constants chosen somewhat arbitrarily.

###### Lemma 1 .

###### Proof.

###### Remark 3 .

##### Testing for a good m 𝑚 m italic_m .

We devise a test that, given $m$, with high probability: (a) rejects if $m<1.5\mathsf{d}_{H}(\lambda)$; and (b) accepts if $m\geq 2\mathsf{d}_{H}(\lambda)$.
Set $z_{0}=-5\lambda/12$. The tests accepts if ${\hat{\mathsf{s}}}(z_{0})>1/\lambda$, and rejects otherwise.
If $m<1.5\mathsf{d}_{H}(\lambda)$ then either there is no solution $\tilde{\lambda}>0$ to ${\mathsf{s}}(-\tilde{\lambda})=1/\lambda$, or there is one, specifically $\tilde{\lambda}=\lambda(1-\frac{1}{m}\mathsf{d}_{H}(\lambda))\leq\lambda/3<-z_
{0}$. Since $z\mapsto{\mathsf{s}}(z)$ is positive and increasing for $z<0$, necessarily ${\mathsf{s}}(z_{0})<{\mathsf{s}}(-\tilde{\lambda})=1/\lambda$. Provided that $m$ is large enough, w.h.p. ${\hat{\mathsf{s}}}(z_{0})<1/\lambda$ as well, hence the test rejects.
If $m\geq 2\mathsf{d}_{H}(\lambda)$ then $\tilde{\lambda}=\lambda(1-\frac{1}{m}\mathsf{d}_{H}(\lambda))\geq\lambda/2>-z_
{0}$, hence ${\mathsf{s}}(z_{0})>{\mathsf{s}}(-\tilde{\lambda})=1/\lambda$. Provided that $m$ is large enough, w.h.p. ${\hat{\mathsf{s}}}(z_{0})>1/\lambda$ as well, hence the test accepts.

##### Searching for m 𝑚 m italic_m by doubling.

Start with $m=m_{0}$ 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 do not halt on any $m<1.5\mathsf{d}_{H}(\lambda)$; (b) we halt on the first $m$ satisfying $m\geq 2\mathsf{d}_{H}(\lambda)$, and therefore return $m$ satisfying $m\leq 4\mathsf{d}_{H}(\lambda)$.

##### Computational complexity.

W.h.p., the doubling procedure stops within $\log_{2}(4\mathsf{d}_{H}(\lambda)/m_{0}))$ iterations. Each iteration $t$ costs at most $O(d^{2}m_{t})$, $m_{t}=2^{t}m_{0}$, so the total cost is $O(d^{2}\min\{\mathsf{d}_{H}(\lambda),m_{0}\})$ and dominated by the cost of multiplying $SHS^{\top}$ for $S\in\mathbb{R}^{m_{t}\times d}$ at the last iteration.

Having found large enough $m$, we can solve the equation (18) by binary search, and form the sketch (7) with $\hat{\lambda}$ in place of $\tilde{\lambda}$.
For the reader’s convenience, the overall adaptive sketching procedure—which
is performed (in parallel) by every worker—is summarized
in Algorithms 1-2.

Figure: Algorithm 1

Figure: Algorithm 2

It is straightforward to establish the asymptotic validity of our adaptive sketching procedure.

###### Theorem 4 .

Set $\eta,D>0$, let $m_{0}$ be an initial sketch size. Assume [Asymp($\tau,H,m_{0},\xi$)] (in particular $m_{0}=\Omega(d)$), and furthermore that $\lambda_{\mathop{\mathrm{rank}}(H)}(H)\geq\tau$.
Let $\hat{m}$ be the output of Algorithm 1, and $(\hat{\lambda},\hat{W})$ be the output of Algorithm 2 (the latter given input $m=\hat{m}$).
Then w.p. $1-O(m_{0}^{-D})$:

- •
$\max\{m_{0},1.5\mathsf{d}_{H}(\lambda)\}\leq\hat{m}\leq\max\{m_{0},4\mathsf{d}
_{H}(\lambda)\}$.
- •
$|\hat{\lambda}-\tilde{\lambda}|=O(m_{0}^{-1/2+\eta})$.
- •
$\|\mathbb{E}[\hat{W}]-W\|=O(m_{0}^{-1/2+\delta})$.

The constants in the $O(\cdot)$ notation may depend on $\eta,D>0$ and also $\lambda>0$.

###### Proof.

This is a straighforward consequence of Theorems 1-3. We omit the details. ∎

###### Theorem 4 .

###### Proof.

#### 3.1 Non-Asymptotic Guarantees for Gaussian Sketches

The guarantees provided by Theorem 4 are asymptotic, in that they do not depend explicitly on crucial problem parameters (such as $\lambda,\|H\|$); in particular, we implicitly assume that $H$ is well-conditioned. When $S\in\mathbb{R}^{m\times d}$ is an i.i.d. Gaussian matrix, $S_{ij}\sim\mathcal{N}(0,1/m)$, we are able to provide non-asymptotic error guarantees for our adaptive sketching algorithm. Our bounds are dimension-free, and remarkably do not depend explicitly on the condition number of $H$. The theorems below are the main technical part of this paper; their proofs are deferred to the appendix.

Theorems 5 and 6 respectively provide non-asymptotic guarantees for Algorithms 1 and 2:

###### Theorem 5 .

Assume that $S$ has i.i.d. Gaussian entries.
Let $\delta\in(0,1)$. Suppose that $m_{0}\geq C(1+\log(1/\delta))$ for large enough constant $C$.
W.p. $1-\delta$, Algorithm 1 outputs $\hat{m}$ such that

$$ $\max\{1.5\mathsf{d}_{H}(\lambda),m_{0}\}\leq\hat{m}\leq\max\{4\mathsf{d}_{H}( \lambda),m_{0}\}\,.$ (19) $$

###### Theorem 6 .

Assume that $S$ has i.i.d. Gaussian entries.
Let $\delta\in(0,1)$ and $\varepsilon>0$. Suppose
Algorithm 2 is run with $m\geq 1.5\mathsf{d}_{H}(\lambda)$ and that moreover
$m\geq C\frac{1+\log(1/\delta)}{\varepsilon^{2}}$ for large enough constant $C$.
W.p. $1-\delta$, the output $\hat{\lambda}$ satisfies $|\hat{\lambda}-\tilde{\lambda}|\leq\varepsilon\tilde{\lambda}$, where $\tilde{\lambda}=\lambda\left(1-\frac{1}{m}\mathsf{d}_{H}(\lambda)\right)$.

Our final non-asymptotic result concerns the concentration of the aggregate inverse Hessian estimate used by the server. To wit, fix $m\geq 1.5\mathsf{d}_{H}(\lambda)$, and let $\hat{W}^{(1)},\ldots,\hat{W}^{(q)}$ be the outputs of $q$ independent runs of Algorithm 2, representing $q$ workers producing estimates of the Newton step in parallel. Denote their average $\bar{W}=\frac{1}{q}\sum_{k=1}^{q}\hat{W}^{(k)}$, and the error matrix:

$$ $\displaystyle\bar{\mathcal{E}}:=(H+\lambda I)^{1/2}\bar{W}(H+\lambda I)^{1/2}-I.$ (20) $$

Later on, we will obtain error bounds for the quasi Newton method in terms of $\|\bar{\mathcal{E}}\|$.

###### Theorem 7 .

Assume that $S$ has i.i.d. Gaussian entries.
Suppose that $m\geq 1.5\mathsf{d}_{H}(\lambda)$, and let $q\geq 1$, $\delta\in(0,1)$ satisfy $1/\delta\leq\exp(O(d))$, $q\leq\exp(O(d))$. Let $\exp(-O(d))\leq\varepsilon\leq O(1)$ be small enough. If moreover

$$ $\displaystyle m\geq C\frac{1}{\varepsilon^{2}}\,,\quad q\geq C\frac{d\log d}{m }\frac{\log(1/\delta)}{\varepsilon^{2}}$ (21) $$

for large enough $C>0$, then w.p. $1-\delta$, $\|\bar{\mathcal{E}}\|\leq\varepsilon$.

###### Remark 4 .

Up to the logarithmic factor (which we believe is an artifact), the dimensional dependence of $q$ in (21) is optimal. To see this, suppose that $H=0,\lambda=1$, so $\hat{W}^{(k)}={S^{(k)}}^{\top}S^{(k)}$. In particular $\mathop{\mathrm{rank}}(\bar{W})=qm$, therefore $\|\bar{\mathcal{E}}-I\|\geq 1$ when $qm<d$. In fact, $\bar{\mathcal{E}}$ is a Wishart (Gaussian sample covariance) matrix corresponding to $mq$ samples, therefore by standard estimates $\mathbb{E}\|\bar{\mathcal{E}}-I\|\asymp\frac{d}{mq}\vee\sqrt{\frac{d}{mq}}$, hence the dependence on $\varepsilon$ is correct as well.

###### Theorem 5 .

###### Theorem 6 .

###### Theorem 7 .

###### Remark 4 .

### 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:\mathbb{R}^{d}\to\bar{\mathbb{R}}$ be convex and twice-differentiable.(^7^77We denote $\bar{\mathbb{R}}=\mathbb{R}\cup\{\infty\}$. With this, the domain of $G$ is $\mathrm{Dom}(G):=\{x\,:\,F(x)<\infty\}$.) We consider generally the Newton-type iteration,

$$ $\displaystyle\theta_{t}$ $\displaystyle=\theta_{t-1}-\alpha_{t}\bar{W}_{t}g_{t},$ (22) $\displaystyle g_{t}$ $\displaystyle:=\nabla G(\theta_{t-1}),$ (23) $$

where $\bar{W}_{t}$ approximates the inverse Hessian,
$\mathsf{H}_{t}:=\nabla^{2}G(\theta_{t-1})$, $\bar{W}_{t}\approx\mathsf{H}^{-1}_{t}$ and $\alpha_{t}$ is a step size.

###### Definition 2 .

We say (22)-(23) implement an $\eta$-accurate Newton method if the following holds. Define the error matrix, which measures how well $\bar{W}_{t}$ approximates the true inverse Hessian $\mathsf{H}_{t}^{-1}$:

$$ $\bar{\mathcal{E}}_{t}:=\mathsf{H}^{1/2}_{t}\bar{W}_{t}\mathsf{H}^{1/2}_{t}-I\,.$ (24) $$

Then $\|\bar{\mathcal{E}}_{t}\|\leq\eta$ for every iteration $t$.

For a desired precision $\varepsilon>0$, namely when one desires $\hat{\theta}$ with $G(\hat{\theta})-G(\theta^{\star})\leq\varepsilon$, $\theta^{\star}:=\operatorname*{\arg\!\min}G(\theta)$,
let $\mathsf{T}(\varepsilon)$ be the smallest $t\geq 1$ (a priori, if it exists) such that $G(\theta_{t})-G(\theta^{\star})\leq\varepsilon$ for all $t\geq\mathsf{T}(\varepsilon)$. This notation suppresses the dependence on $G$ and the sequence of step sizes $\alpha_{t}$. We consider a variation where the step size is chosen by backtracking line search, see Algorithm 3.

Figure: Algorithm 3 (Backtracking Line Search)

This paper focuses on the case when $G$ is self-concordant, which is a standard setting in the literature on second order convex optimization, cf. \citepboyd2004convex. We remark that previous works, for example \citepwang2018giant,dereziński2019distributed, derived convergence guarantees for the approximate Newton method in a different setting, namely when the Hessian of $G$ is Lipschitz. Due to space constraint, we do not present their results here, though they are applicable in our setting as well. For self-concordant objectives, we 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 a univariate convex function $G:\mathbb{R}\to\bar{\mathbb{R}}$ is self-concordant if it is thrice-differentiable and $|G^{\prime\prime\prime}(x)|\leq 2G^{\prime\prime}(x)^{3/2}$ for any $x$ in its domain. A multivariate convex function $G:\mathbb{R}^{d}\to\bar{\mathbb{R}}$ is self-concordant if any 1-dimensional restriction $G_{x,y}(t)=G(x+ty)$ is self-concordant.

When $f$ is self-concordant, the objective $G(\theta)=f(\theta)+\frac{\lambda}{2}\|\theta\|^{2}$ in (1) is self-concordant as well, as the sum of self-concordant function is itself self-concordant, cf. \citepboyd2004convex.

###### Theorem 8 .

Suppose that (a) $G$ is self-concordant; (b) the Newton method is $\eta$-accurate with $\eta<1/5$; and (c) $\alpha_{t}$ is chosen by backtracking line search. For large enough numerical $C>0$, let

$$ $T_{0}=C\frac{G(\theta_{0})-G(\theta^{\star})}{ab},\quad T_{1}(\varepsilon)=C \log\log(1/\varepsilon),\quad T_{2}(\varepsilon)=C\frac{\log(\eta/\varepsilon) }{\log(1/\varepsilon)}\,,$ $$

where $a,b$ are the parameters used in backtracking line search (Algorithm 3).
Then

$$ $\mathsf{T}(\varepsilon)\leq\begin{cases}T_{0}+T_{1}(\varepsilon)\quad&\textrm{ if}\quad\varepsilon\geq\eta\\ T_{0}+T_{1}(\eta)+T_{2}(\varepsilon)\quad&\textrm{if}\quad\varepsilon<\eta\end {cases}\,.$ (25) $$

The proof of Theorem 25 builds on and extends an argument of \citeppilanci2015newton. Interestingly, Theorem 25 suggests that the approximate Newton method exhibits two stages: (1) While the error is above the noise floor $\varepsilon\geq\eta$, the error decays quadratically fast. (2) Upon reaching the noise floor, $\varepsilon\leq\eta$, the convergence rate becomes linear; however, the linear rate is proportional to $\eta$, that is the error decays as $\sim(O(\eta))^{t}$. When $\eta$ is very small this decay rate is very fast.

Finally, Theorem 25 may be readily combined with Theorem 7 to yield error bounds for the entire end-to-end parallel method: under the setting of Theorem 7, w.h.p. (and up to $\mathrm{polylog}(d)$ factors),

$$ $\eta=\tilde{O}\left(\frac{1}{\sqrt{m}}+\sqrt{\frac{d}{mq}}\right)\,.$ $$

###### Definition 2 .

###### Definition 3 .

###### Theorem 8 .

### 5 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(^8^88We choose $H$ to be diagonal w.l.o.g. since the Gaussian sketch is invariant to orthogonal transformations.)
$H$ with eigenvalues $(h_{1},\ldots,h_{d})$ exhibiting polynomial spectral decay $h_{k}=k^{-\alpha}$, where $\alpha>0$ is a parameter. The dimension $d=10^{4}$ is fixed and large, and $\lambda=1$.
The sketching matrix is i.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-Carlo trials, and report the sketching dimension found, see Table 1. Remarkably, our algorithm consistently succeeds, namely finds a sketching dimension that is within $1.5\mathsf{d}_{H}(\lambda)\leq m\leq\max\{m_{0},4\mathsf{d}_{H}(\lambda)\}$ - per Theorems 4 and 5.

**Table 1: Success rate of Algorithm 1.**
| $\alpha$ | $\mathsf{d}_{H}(\lambda)$ | $m_{0}$ | Success rate | Avg. dim. |
| --- | --- | --- | --- | --- |
| $1$ | $8.8$ | $10$ | $1.0$ | $20$ |
| $2/3$ | $59.7$ | $10$ | $1.0$ | $160$ |
| $1/2$ | $190.4$ | $10$ | $1.0$ | $640$ |

Next, we demonstrate that Algorithm 2 produces low-bias estimates of $W$.
Note that the average Frobenius error $\|\hat{W}-W\|_{F}^{2}/d^{2}$ is a reasonable proxy for the bias of $\hat{W}$, see Theorems 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 $m$ increasing, starting from $m=1.5\mathsf{d}_{H}(\lambda)$, for the Gaussian sketch. As $m$ increases, the bias decreases (roughly at rate $1/\sqrt{m}$). We also plot the bias proxy for a naive estimator of $W$, without bias correction ($\tilde{\lambda}=\lambda$); clearly, our algorithm has substantially smaller bias throughout.

Figure: 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). Leftmost plot: the sketching dimension found by Algorithm 1 on ensemble (L). Shaded area: 20%-80% confidence interval; We take $T=10$ Monte-Carlo trials.
Refer to caption: extracted/5895313/figures/fig1_revise2-01.png

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—using publicly available UCI datasets \citepchang2011libsvm. Due to space constraints, the details appear in the appendix.
We compare our method with an approximate parallel Newton method where Hessians are sketched ($m$ chosen by Algorithm 1), but no debiasing is done; that is, $\tilde{\lambda}=\lambda$. It is evident from our results, summarized in Figure 2, that bias correction accelerates convergence, often substantially.

Figure: 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 a $20\%$-$80\%$ confidence interval.
Refer to caption: x1.png

### Acknowledgments

We are grateful to Daniel Lejeune for inspiring conversations about his work \citeplejeune2022asymptotics.

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-0242; 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 was supported in part by the NSF under Grant DMS-1811614 (PI: Donoho), and in part by the personal faculty funds of Prof. David Donoho, for whom he wishes to extend his sincere appreciation.

## Appendix

### Appendix A Additional Experiments and Details

#### A.1 Extension of Table 1

We run experiments with random sketching matrices with i.i.d. entries of the following types: (1) Gaussian (G): $S_{i,j}\sim\mathcal{N}(0,1/m)$; Rademacher (R): $S_{i,j}\sim\mathrm{Unif}(\pm 1/\sqrt{m})$; (3) sparse Rademacher (SR): $S_{i,j}=0$ w.p. $1-p$ and $S_{i,j}\sim\mathrm{Unif}(\pm 1/\sqrt{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:

**Table 2: Success of Algorithm 1. $T=20$ Monte-Carlo trials.**
| $\alpha$ | $\mathsf{d}_{H}(\lambda)$ | $m_{0}$ | Sketch | Success rate | Avg. dim. |
| --- | --- | --- | --- | --- | --- |
| $1$ | $8.8$ | $10$ | G | $1.0$ | $20$ |
| $2/3$ | $59.7$ | $10$ | G | $1.0$ | $160$ |
| $1/2$ | $190.4$ | $10$ | G | $1.0$ | $640$ |
| $1$ | $8.8$ | $10$ | R | $1.0$ | $20$ |
| $2/3$ | $59.7$ | $10$ | R | $1.0$ | $160$ |
| $1/2$ | $190.4$ | $10$ | R | $1.0$ | $640$ |
| $1$ | $8.8$ | $10$ | SR | $1.0$ | $20$ |
| $2/3$ | $59.7$ | $10$ | SR | $1.0$ | $160$ |
| $1/2$ | $190.4$ | $10$ | SR | $1.0$ | $640$ |

Our method appears to be remarkably robust, per the finding in Table 2, to the extent that the table may seem somewhat pointless to include. We nonetheless provide it (and code to run this experiment) as a sort of sanity check.

#### A.2 Details for Figure 1

We provide additional details for the experiments in Figure 1. In the middle and rightmost plots, the estimator $\bar{W}$ is obtained by averaging $q=500$ local sketched Hessians, both with and without bias correction.

Ensemble (L): We use $L=10^{-3}$. The matrix $H$ is $H=X^{\top}X$, where $X\in\mathbb{R}^{n\times d}$ has the singular value decomposition (SVD) $X=UDV^{\top}$; $U\in\mathbb{R}^{n\times d},V\in\mathbb{R}^{d\times d}$ have orthonormal columns and $D\in\mathbb{R}^{d\times d}$ is diagonal. $U,V$ are sampled uniformly (Haar) from the corresponding Stieffel manifolds (have a rotationally invariant distribution).
For the leftmost plot $n=10^{4}$, and $d$ varies. For the middle plot, we fix $d=10^{3}$. The diagonal entries of $D$ are random, and satisfy

$$ $D_{k,k}=(0.9+\varepsilon_{k})^{k/2},\quad\varepsilon_{k}\sim\mathcal{N}(0,10^{ -4})\,.$ $$

Ensemble (R): We use $\lambda=10^{-5}$. Similarly to ensemble (L), $H=X^{\top}X$ for $X\in\mathbb{R}^{n\times d}$, with the singular vectors of $X$ of similarly uniformly random. We use $n=10^{4}$, $d=10^{3}$. The singular values of $X$ are
non-random, with $D_{k,k}=(k/d)^{2}$.

#### A.3 Additional Synthetic Optimization Tasks

Figure: Figure 3: Convergence of the parallel Newton method with Hessian sketching, on synthetic optimization tasks.
Refer to caption: extracted/5895313/figures/fig3_revise-01.png

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, we consider ridge and logistic regression. We test both Gaussian and Rademacher sketches (see Section A.1).

##### Distribution of covariates (design).

We generate designs $X\in\mathbb{R}^{n\times d}$ with $n=10^{4}$ and $d=500$. Denote the SVD $X=UDV^{\top}$. The singular vectors $U\in\mathbb{R}^{n\times d},V\in\mathbb{R}^{d\times d}$ are uniformly random. The singular values decay exponentially, $D_{k,k}=0.99^{k/2}$.

##### Parameters.

In all our experiments, $\lambda=10^{-3}$. For experiments with ridge regression, the number of workers is $q=10$; for logistic regression, it is $q=20$.

##### Ridge regression.

The responses are generated according to $y_{i}=x_{i}^{\top}\theta^{\star}+\mathcal{N}(0,0.01)$, where the ground-truth regressor $\theta^{\star}\in\mathbb{R}^{d}$ is generated according to $\theta^{\star}\sim\mathcal{N}(0,I)$.
The unregularized objective $F(\theta)$ is the L2 loss:

$$ $F(\theta)=\frac{1}{n}\sum_{i=1}^{n}(x_{i}^{\top}\theta-y_{i})^{2}=\|X\theta-y \|^{2},$ $$

where $x_{1},\ldots,x_{n}\in\mathbb{R}^{d}$ are the rows of $X$.
Recall that we aim to minimize the regularized objective $G(\theta):=F(\theta)+\frac{\lambda}{2}\|\theta\|^{2}$.

##### Logistic regression.

We generate labels according to

$$ $y_{i}=(\mathop{\mathrm{sign}}\left[x_{i}^{\top}\theta^{\star}+\mathcal{N}(0,10 ^{4})\right]+1)/2\in\{0,1\}\,,$ $$

where $\theta^{\star}\sim\mathcal{N}(0,I)$. The function $F(\theta)$ is the log loss:

$$ $F(\theta)=-\frac{1}{n}\sum_{i=1}^{n}\left[y_{i}\log\left(\frac{1}{1+e^{-x_{i}^ {\top}\theta}}\right)+(1-y_{i})\log\left(1-\frac{1}{1+e^{-x_{i}^{\top}\theta}} \right)\right]\,.$ $$

##### 

Figure 3 plots the error $G(\theta_{t})-\operatorname*{\arg\!\min}_{\theta}G(\theta)$ 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. It is evident that bias correction improves on the the convergence speed of the algorithm over the uncorrected ($\tilde{\lambda}=\lambda$) alternative.

#### A.4 Details for Experiments on UCI Data Sets (Figure 2 )

Figure: Figure 4: Convergence rate for optimization tasks on additional UCI data sets. Shaded area corresponds to 20%-80% confidence interval.
Refer to caption: x9.png

We provide additional details for the experiments reported in Figure 2, and present results for additional data sets from the UCI data repository \citepchang2011libsvm.

We run experiments directly on the available data, without additional preprocessing. The objectives, 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 $\lambda=10^{-3}$ for all experiments. We vary the number of workers between $q\in\{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 optimality gap $G(\theta_{t})-\operatorname*{\arg\!\min}_{\theta}G(\theta)$, $G$ being the regularized objective, as a function of the iteration number $t$. The bottom plots are the inverse Hessian normalized Frobenius error $\|W_{t}-\bar{W}_{t}\|_{F}^{2}/d^{2}$. We see that as $q$ grows larger, the algorithm tends to converge faster; though the benefit of further increasing $q$ tends to diminish when $q$ is large.

### Appendix B Non-Asymptotic Random Matrix Theory Results

The proofs of Theorems 5-7, rely on 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 the random matrices $SHS^{\top},H^{1/2}S^{\top}SH^{1/2}$ where $H\succeq 0$ and $S\in\mathbb{R}^{m\times d}$ is an i.i.d. Gaussian matrix $S_{ij}\sim\mathcal{N}(0,1/m)$.

The following result concerns the smallest eigenvalue of the matrix $SHS^{\top}\in\mathbb{R}^{m\times m}$.

###### Theorem 9 .

Suppose $S$ has i.i.d. Gaussian entries.
There are constants $c_{1},c_{2},c,C>0$ such that the following holds.
Set any $\eta>0$.
If
$m\leq c_{1}\mathsf{d}_{H}(\eta)$ then w.p. $1-Ce^{-cm}$,

$$ $\lambda_{\min}(SHS^{\top})\geq c_{2}\eta\frac{1}{m}\mathsf{d}_{H}(\eta)\,.$ $$

The following two results concern the closedness of the empirical (companion) Stieltjes transform

$$ $\displaystyle{\hat{\mathsf{s}}}(z):=\frac{1}{m}\mathop{\mathrm{tr}}\left(SHS^{ \top}-zI\right)^{-1},$ (26) $$

to the deterministic function ${\mathsf{s}}(z)$, the solution of the Marchenko-Pastur equation (10).

The following concentration inequality is well known, cf. [bai2010spectral].

###### Lemma 2 .

Suppose that $S=\frac{1}{\sqrt{m}}[r_{1}|\cdots|r_{m}]^{\top}$ where
the rows $r_{1},\ldots,r_{m}\in\mathbb{R}^{d}$ are independent random variables. For any $z<0$ and $t\geq 0$,

$$ $\Pr(|{\hat{\mathsf{s}}}(z)-\mathbb{E}{\hat{\mathsf{s}}}(z)|\geq t)\leq 2e^{- cmz^{2}t^{2}}\,.$ $$

We provide a self-contained proof for completeness, see Section F.

Next, we show that $\mathbb{E}{\hat{\mathsf{s}}}(z)$ approximately satisfies (10), with explicit non-asymptotic error bounds. Note that (10) can be be equivalently rewritten (assuming $z\neq 0$) as

$$ ${\mathsf{s}}(z)=-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{H}(1/{\mathsf{s}}(z ))\right)\,.$ (27) $$

###### Theorem 10 .

Suppose that $S$ has i.i.d. Gaussian entries.
For any $z<0$, $\mathbb{E}{\hat{\mathsf{s}}}(z)$ satisfies:

$$ $\displaystyle\mathbb{E}{\hat{\mathsf{s}}}(z)$ $\displaystyle=-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{H}({1/\mathbb{E}{\hat {\mathsf{s}}}(z)})\right)+e_{m}(z;H)\,,$ (28) $$

where the error term $e_{m}(z;H)$ satisfies

$$ $\displaystyle e_{m}(z;H)$ $\displaystyle\lesssim\frac{\beta}{\alpha}\frac{1}{|z|\sqrt{m}}\,,$ (29) $$

with $\alpha=\alpha_{m}(z;H),\beta=\beta_{m}(z;H)$ defined as

$$ $\displaystyle\alpha$ $\displaystyle:=|z|\mathbb{E}{\hat{\mathsf{s}}}(z)\,,$ (30) $\displaystyle\beta$ $\displaystyle:=\frac{1}{m}\mathsf{d}_{H}({1/\mathbb{E}{\hat{\mathsf{s}}}(z)})\,.$ (31) $$

The proof of Theorem 10, given in Section G, follows along a well-known computation technique in random matrix theory, and explicitly leverages the orthogonal invariance of the Gaussian distribution. Note that the error terms $\alpha,\beta$ depend on $\mathbb{E}{\hat{\mathsf{s}}}(z)$ itself. This feature, which somewhat complicates the proofs of Theorems 5-6, is generally not a mere artifact of the calculation, see e.g. the bounds of \citepknowles2017anisotropic under asymptotic conditions.

###### Theorem 9 .

###### Lemma 2 .

###### Theorem 10 .

### Appendix C Proof of Theorem 5

Denote for brevity, per the description of Algorithm 1,

$$ $z_{0}:=-\frac{5}{12}\lambda\,.$ (32) $$

The proof of Theorem 5 consists of two parts. First, we show that if $m<1.5\mathsf{d}_{H}(\lambda)$, then with high probability (w.h.p.) ${\hat{\mathsf{s}}}_{m}(z_{0})<1/\lambda$, hence such $m$ would not be returned. Then, we show that when $m\geq 2\mathsf{d}_{H}(\lambda)$, then w.h.p. ${\hat{\mathsf{s}}}_{m}(z_{0})\geq 1/\lambda$, 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\leq 4\mathsf{d}_{H}(\lambda)$.

The next lemma shows that Algorithm 1 does not return $m$ which is *much* smaller than $\mathsf{d}_{H}(\lambda)$.

###### Lemma 3 .

Suppose that $S$ has i.i.d. Gaussian entries. There is a (small) constant $c_{0}$ such that if $m<c_{0}\mathsf{d}_{H}(\lambda)$ then

$$ $\displaystyle\Pr({\hat{\mathsf{s}}}(z_{0})>1/\lambda)\leq Ce^{-cm}\,.$ $$

###### Proof.

The lemma is a consequence of Theorem 9, which bounds the smallest eigenvalue of $SHS^{\top}$.
Set $c_{0}=\min\{c_{1},12c_{2}/7\}$ for $c_{1},c_{2}$ from Theorem 9. Assume that $m<c_{0}\mathsf{d}_{H}(\lambda)$. Since in particular $m<c_{1}\mathsf{d}_{H}(\lambda)$, the following event holds w.p. $1-Ce^{-cm}$:

$$ ${\hat{\mathsf{s}}}(z_{0})\leq\frac{1}{\lambda_{\min}(SHS^{\top})-z_{0}}\leq \frac{1}{c_{2}\frac{1}{m}\mathsf{d}_{H}(\lambda)\cdot\lambda+\frac{5}{12} \lambda}\,.$ $$

Moreover, since $c_{2}/c_{0}\geq 7/12$,

$$ $c_{2}\frac{1}{m}\mathsf{d}_{H}(\lambda)>c_{2}\frac{1}{c_{0}}\geq\frac{7}{12}\,.$ $$

Hence, under this event, ${\hat{\mathsf{s}}}(z_{0})<1/\lambda$.
∎

Next, we treat the regime of moderate $m$, namely when $c_{0}\mathsf{d}_{H}(\lambda)\leq m<1.5\mathsf{d}_{H}(\lambda)$. The main tool here is the non-asymptotic estimate of Theorem 10. The reason we have to assume $m\gtrsim\mathsf{d}_{H}(\lambda)$ is that only then we can guarantee the error term is controlled.

The following simple lemma is useful.

###### Lemma 4 .

For any $\mu>0$ and $\gamma\geq 1$,

$$ $\mathsf{d}_{H}(\gamma\mu)\leq\mathsf{d}_{H}(\mu)\leq\gamma\mathsf{d}_{H}( \gamma\mu)\,.$ (33) $$

###### Proof.

The inequality $\mathsf{d}_{H}(\gamma\mu)\leq\mathsf{d}_{H}(\mu)$ is trivial. As for the other one,

$$ $\displaystyle\mathsf{d}_{H}(\gamma\mu)=\mathop{\mathrm{tr}}(H(H+\gamma\mu I)^{ -1})=\frac{1}{\gamma}\mathop{\mathrm{tr}}(H(\frac{1}{\gamma}H+\mu I)^{-1}) \overset{(\star)}{\geq}\frac{1}{\gamma}\mathop{\mathrm{tr}}(H(H+\mu I)^{-1})= \frac{1}{\gamma}\mathsf{d}_{H}(\mu)\,,$ $$

where $(\star)$ follows because $H\succeq 0$, hence $\frac{1}{\gamma}H+\mu I\preceq H+\mu I$ and therefore $H(\frac{1}{\gamma}H+\mu I)^{-1}\succeq H(H+\mu I)^{-1}$.
∎

###### Lemma 5 .

Suppose that $S$ has i.i.d. Gaussian entries. Moreover, suppose that $c_{0}\mathsf{d}_{H}(\lambda)\leq m<1.5\mathsf{d}_{H}(\lambda)$.
There is a numerical constant $M$ such that if moreover $m\geq M$ then

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\leq\frac{14}{15}\cdot\frac{1}{\lambda}\,.$ $$

In particular, $\Pr({\hat{\mathsf{s}}}(z_{0})>1/\lambda)\leq Ce^{-cm}$.

###### Proof.

Note that the second part follows directly from the first part and Lemma 2; we therefore focus on bounding $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})$.
Suppose it were the case that $\mathbb{E}{\hat{\mathsf{s}}}_{m}(z_{0})\geq\eta\frac{1}{\lambda}$ where $\eta\leq 1$ is some constant. Apply Theorem 10 with $z=z_{0}$.
Let us bound the parameters $\alpha,\beta$ used for error control. First,

$$ $\displaystyle\alpha:=|z_{0}|\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\frac{5}{12 }\lambda\cdot\eta\frac{1}{\lambda}\geq\frac{5}{12}\eta\,.$ (34) $$

Next, since $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\leq\frac{1}{|z_{0}|}=\frac{12}{5}\frac{1}{\lambda}$ we have $1/\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\frac{5}{12}\lambda$,
so

$$ $\displaystyle\beta:=\frac{1}{m}\mathsf{d}_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z _{0}))\leq\frac{1}{m}\mathsf{d}_{H}(\frac{5}{12}\lambda)\overset{(\star)}{\leq }\frac{12}{5}\frac{1}{m}\mathsf{d}_{H}(\lambda)\overset{(\star\star)}{\leq} \frac{12}{5}\frac{1}{c_{0}}\,,$ (35) $$

where $(\star)$ follows from Lemma 4 and $(\star\star)$ from the assumption $m\geq c_{0}\mathsf{d}_{H}(\lambda)$. In particular, the error $e_{m}(z_{0},H)$ in (29) can be bound as $e_{m}(z_{0},H)=O_{\eta}(\frac{1}{\lambda\sqrt{m}})$. Furthermore, as $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\eta\frac{1}{\lambda}$ by assumption,

$$ $\beta\geq\frac{1}{m}\mathsf{d}_{H}(\frac{1}{\eta}\lambda)\geq\eta\frac{1}{m} \mathsf{d}_{H}(\lambda)>\frac{2}{3}\eta\,,$ (36) $$

where we used Lemma 4 and the assumption $m<1.5\mathsf{d}_{H}(\lambda)$. With these estimates in mind, consider (28),

$$ $\displaystyle\mathbb{E}{\mathsf{s}}(z_{0})=-\frac{1}{z_{0}}(1-\beta-z_{0} \varepsilon_{m}(z_{0},H))=\frac{1}{\lambda}\cdot\frac{12}{5}(1-\beta-z_{0}e_{m }(z_{0},H))\leq\frac{1}{\lambda}\cdot\frac{12}{5}(1-\frac{2}{3}\eta+O_{\eta}( \frac{1}{\sqrt{m}}))\,.$ $$

As $\mathbb{E}{\hat{\mathsf{s}}}_{m}(z_{0})\geq\eta\frac{1}{\lambda}$, we deduce from the above

$$ $\displaystyle\eta\leq\frac{12}{5}(1-\frac{2}{3}\eta+O_{\eta}(\frac{1}{\sqrt{m} }))\,,$ $$

so that rearranging,

$$ $\displaystyle\eta\leq\frac{12}{13}(1+O_{\eta}(\frac{1}{\sqrt{m}}))\,.$ $$

By requiring that $m\geq C_{\eta}$ for large enough $C_{\eta}$, we can ensure that the error term is as small a constant as we like, so that necessarily (for example) $\eta\leq\frac{14}{15}$.

∎

Next, we show that if the algorithm has reached $m\geq 2\mathsf{d}_{H}(\lambda)$, then it will halt in that iteration. We start with two helper lemmas.

###### Lemma 6 .

Suppose that $\mathbb{E}[S^{\top}S]=I$. Then for all $z<0$,

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z)\geq-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{ H}(-z)\right)\,.$ (37) $$

(Note that since $\mathbb{E}{\hat{\mathsf{s}}}(z)\leq-1/z$, we have $\mathsf{d}_{H}(-z)\geq\mathsf{d}_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z))$, so that (37) is perfectly consistent with (28).)

###### Proof.

Recall that the matrices $SHS^{\top}\in\mathbb{R}^{m\times m}$ and $H^{1/2}S^{\top}SH^{1/2}$ have the same non-zero eigenvalues, and therefore

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z):=\frac{1}{m}\mathbb{E}\mathop{\mathrm{tr}}(SHS ^{\top}-zI)^{-1}=\frac{1}{m}\mathbb{E}\mathop{\mathrm{tr}}(H^{1/2}S^{\top}SH^{ 1/2}-zI)^{-1}+\frac{1}{m}(d-m)\frac{1}{z}\,.$ $$

Note that since $z<0$, the matrix function $A\mapsto\mathop{\mathrm{tr}}(A-zI)^{-1}$ is convex on the cone of PSD matrices $A\succeq 0$. Thus, by Jensen’s inequality,

$$ $\frac{1}{m}\mathbb{E}\mathop{\mathrm{tr}}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}\geq \frac{1}{m}\mathop{\mathrm{tr}}(\mathbb{E}[H^{1/2}S^{\top}SH^{1/2}]-zI)^{-1}= \frac{1}{m}\mathop{\mathrm{tr}}(H-zI)^{-1}\,,$ $$

so that

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z)\geq\frac{1}{m}\mathop{\mathrm{tr}}(H-zI)^{-1}+ \frac{1}{m}(d-m)\frac{1}{z}=\frac{1}{m}\mathop{\mathrm{tr}}[(H-zI)^{-1}+\frac{ 1}{z}I]-\frac{1}{z}=-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{H}(-z)\right)\,.$ $$

∎

The next lemma formally establishes that the inverse function of ${\mathsf{s}}(z)$ is increasing on an appropriate interval. As in (11), denote $\Psi:(0,\infty)\to\mathbb{R}$,

$$ $\Psi(s)=-\frac{1}{s}(1-\frac{1}{m}\mathsf{d}_{H}(1/s))\,,$ (38) $$

so that for all $z<0$, $\Psi({\mathsf{s}}(z))=z$. Note however that in general, $\mathop{\mathrm{Range}}({\mathsf{s}})$ may only be a proper subset of $(0,\infty)$; nonetheless, (38) is well-defined.

###### Lemma 7 .

Let $\mu>0$, and suppose that $m\geq\mathsf{d}_{H}(\mu)$. Then $\Psi$ is increasing on $(0,1/\mu)$.

###### Proof.

Write

$$ $\Psi(s)=-\frac{1}{s}+\frac{1}{m}\mathop{\mathrm{tr}}(H(I+sH)^{-1})\,,$ $$

so that its derivative is

$$ $\displaystyle\Psi^{\prime}(s)$ $\displaystyle=\frac{1}{s^{2}}-\frac{1}{m}\mathop{\mathrm{tr}}(H^{2}(I+sH)^{-2})$ $\displaystyle=\frac{1}{s^{2}}\left(1-\frac{1}{m}\mathop{\mathrm{tr}}(H^{2}(H+1 /s)^{-2})\right)$ $\displaystyle\geq\frac{1}{s^{2}}\left(1-\frac{1}{m}\mathop{\mathrm{tr}}(H(H+1/ s)^{-1})\right)=\frac{1}{s^{2}}\left(1-\frac{1}{m}\mathsf{d}_{H}(1/s)\right)\,,$ (39) $$

where the inequality follows since $0\preceq H(H+1/s)^{-1}\preceq I$. When $s<1/\mu$, we have $1/s>\mu$ so that $\mathsf{d}_{H}(1/s)<\mathsf{d}_{H}(\mu)$, hence $\Psi^{\prime}(s)>0$. ∎

We now prove:

###### Lemma 8 .

Suppose that $S$ has i.i.d. Gaussian entries, and moreover that $m\geq 2\mathsf{d}_{H}(\lambda)$.
There is a numerical $M$ such that whenever $m\geq M$,

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\frac{24}{23}\cdot\frac{1}{\lambda}\,.$ $$

In particular, $\Pr({\hat{\mathsf{s}}}(z_{0})\leq 1/\lambda)\leq Ce^{-cm}$.

###### Proof.

Suppose that $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\leq\frac{1}{(1-\eta)}\frac{1}{\lambda}$ for some $\eta\in(0,1)$; we shall show that $\eta$ has to be somewhat large so to not generate a contradiction.
First, by Lemma 6 (with $z=-\lambda$),

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\mathbb{E}{\hat{\mathsf{s}}}(-\lambda)= \frac{1}{\lambda}(1-\frac{1}{m}\mathsf{d}_{H}(\lambda))=\frac{1}{2}\cdot\frac{ 1}{\lambda}\,,$ $$

where the last inequality follows since $m\geq 2\mathsf{d}_{H}(\lambda)$. Of course, this lower bound does not suffice for our purposes: we want a bound which is $>\frac{1}{\lambda}$ (strictly). It does allow us, however, to control the error term in Theorem 10 at $z=z_{0}$ (specifically, to lower bound $\alpha$). We can bound the parameters $\alpha,\beta$ in Theorem 10 as

$$ $\displaystyle\alpha$ $\displaystyle:=|z_{0}|\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\frac{5}{12} \lambda\cdot\frac{1}{2}\frac{1}{\lambda}\geq\frac{5}{24}\,,$ $\displaystyle\beta$ $\displaystyle:=\frac{1}{m}\mathsf{d}_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z_{0}) )\leq\frac{1}{m}\mathsf{d}_{H}(-z_{0})=\frac{1}{m}\mathsf{d}_{H}(\frac{5}{12} \lambda)\leq\frac{12}{5}\frac{1}{m}\mathsf{d}_{H}(\lambda)\leq\frac{6}{5}\,,$ $$

where, in upper bounding $\beta$, we used $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\leq\frac{1}{-z_{0}}$ and Lemma 4.
With these, Theorem 10 implies

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})=-\frac{1}{z_{0}}(1-\frac{1}{m}\mathsf{d}_{ H}(1/\mathbb{E}{\hat{\mathsf{s}}}_{m}(z_{0})))+O(\frac{1}{\lambda\sqrt{m}})\,.$ $$

Multiplying by $z_{0}$ and dividing by $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})$,

$$ $\displaystyle z_{0}$ $\displaystyle=-\frac{1}{\mathbb{E}{\hat{\mathsf{s}}}(z)}\left(1-\frac{1}{m} \mathsf{d}_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z))\right)+O(\frac{\lambda}{ \sqrt{m}})$ $\displaystyle=\Psi(\mathbb{E}{\hat{\mathsf{s}}}(z))+O(\frac{\lambda}{\sqrt{m}} )\,,$ $$

where in modifying the error term we used $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\gtrsim\frac{1}{\lambda}$.

Note that $m\geq 2\mathsf{d}_{H}(\lambda)$ implies $m\geq\mathsf{d}_{H}(\frac{1}{2}\lambda)$ by Lemma 4. By Lemma 7, this implies that $\Psi$ is increasing on $(0,2\frac{1}{\lambda})$. Thus, if we further assume that $\eta<1/2$,
then $\Psi$ is increasing on $(0,\frac{1}{1-\eta}\frac{1}{\lambda})$, and so $\Psi(\mathbb{E}{\hat{\mathsf{s}}}(z))\leq\Psi(\frac{1}{1-\eta}\frac{1}{\lambda})$. Using this with the previous display, and writing $z_{0}=-\frac{5}{12}\lambda$,

$$ $\displaystyle-\frac{5}{12}\lambda$ $\displaystyle\leq\Psi(\frac{1}{1-\eta}\frac{1}{\lambda})+O(\frac{\lambda}{ \sqrt{m}})$ $\displaystyle=-(1-\eta){\lambda}\left(1-\frac{1}{m}\mathsf{d}_{H}((1-\eta) \lambda)\right)+O(\frac{\lambda}{\sqrt{m}})$ $\displaystyle\overset{(\star)}{\leq}-(1-\eta){\lambda}\left(1-\frac{1}{1-\eta} \frac{1}{m}\mathsf{d}_{H}(\lambda)\right)+O(\frac{\lambda}{\sqrt{m}})$ $\displaystyle\overset{(\star\star)}{\leq}-(1-\eta){\lambda}\left(1-\frac{1}{1- \eta}\frac{1}{2}\right)+O(\frac{\lambda}{\sqrt{m}})$ $\displaystyle=-\lambda(\frac{1}{2}-\eta)+O(\frac{\lambda}{\sqrt{m}})\,,$ $$

where $(\star)$ uses $\mathsf{d}_{H}((1-\eta)\lambda)\leq\frac{1}{1-\eta}\mathsf{d}_{H}(\lambda)$ (Lemma 4) and $(\star\star)$ uses the assumption $m\geq 2\mathsf{d}_{H}(\lambda)$. Dividing by $(-\lambda)$,

$$ $\displaystyle\frac{5}{12}\geq\frac{1}{2}-\eta+O(\frac{1}{\sqrt{m}})\,.$ $$

For large enough $m\geq M$, the $O(\cdot)$ term is $\leq\frac{1}{24}$, therefore $\eta\geq\frac{1}{24}$. That is, we deduce that necessarily $\mathbb{E}{\hat{\mathsf{s}}}(z_{0})\geq\frac{1}{1-\frac{1}{24}}\frac{1}{
\lambda}=\frac{24}{23}\frac{1}{\lambda}$.
∎

We are ready to conclude the proof of Theorem 5.

###### Proof.

(Of Theorem 5.)
Let $m_{j}=2^{j}m_{0}$, $j=0,1,\ldots$, be the value of $m$ at iteration $j$ of the algorithm. Set

$$ $J^{-}=\lceil\log_{2}(1.5\mathsf{d}_{H}(\lambda)/m_{0})\rceil,\qquad J^{+}= \lceil\log_{2}(2\mathsf{d}_{H}(\lambda)/m_{0})\rceil\,,$ $$

chosen so that: 1)
$J^{-}$ is the smallest $j$ such that $m_{j}\geq 1.5\mathsf{d}_{H}(\lambda)$; 2) $2\mathsf{d}_{H}(\lambda)\leq m_{J^{+}}<4\mathsf{d}_{H}(\lambda)$. Thus, to establish (19) it suffices to show that the algorithm halts at an iteration $j$ satisfying $J^{-}\leq j\leq J^{+}$.
Accordingly, the error probability can be bounded as

$$ $\displaystyle\Pr(\textrm{\eqref{eq:thm:Alg1-Gauss} does not hold})\leq\sum_{j= 0}^{J^{-}-1}\Pr({\hat{\mathsf{s}}}_{m_{j}}(z_{0})>1/\lambda)+\Pr({\hat{\mathsf {s}}}_{m_{J^{+}}}(z_{0})\leq 1/\lambda)\,.$ $$

Assume that $m_{0}\geq M$ is large enough. By Lemmas 3 and 5, for all $j\leq J^{-}-1$, $\Pr({\hat{\mathsf{s}}}_{m_{j}}(z_{0})>1/\lambda)\leq Ce^{-cm_{j}}$. By Lemma 8, $\Pr({\hat{\mathsf{s}}}_{m_{J^{+}}}(z_{0})\leq 1/\lambda)\leq Ce^{-cm_{J^{+}}}$. Thus,

$$ $\displaystyle\Pr(\textrm{\eqref{eq:thm:Alg1-Gauss} does not hold})\leq\sum_{j= 0}^{J^{-}-1}Ce^{-c2^{j}m_{0}}+Ce^{-c2^{J^{+}}m_{0}}\leq\sum_{j=0}^{\infty}Ce^{ -c2^{j}m_{0}}\lesssim e^{-cm_{0}}\,.$ $$

We deduce that if $m_{0}\gtrsim\log(1/\delta)$, then the error probability is $\leq\delta$.
∎

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

###### Lemma 7 .

###### Proof.

###### Lemma 8 .

###### Proof.

###### Proof.

### Appendix D Proof of Theorem 6

###### .

Recall the (deterministic) Marchenko-Pastur equation (10); rearranging yields

$$ $z=-\frac{1}{{\mathsf{s}}(z)}\left(1-\frac{1}{m}\mathsf{d}_{H}(1/{\mathsf{s}}(z ))\right)\,,$ (40) $$

so that setting ${\mathsf{s}}(-\tilde{\lambda})=1/\lambda$ yields $\tilde{\lambda}=\lambda(1-\frac{1}{m}\mathsf{d}_{H}(\lambda))$. Further recall that we do not have access to ${\mathsf{s}}(z)$, but instead to its random counterpart ${\hat{\mathsf{s}}}(z)$. We know that its expectation $\mathbb{E}{\hat{\mathsf{s}}}(z)$ approximately satisfies the Marcheko-Pastur equation with an error term (Theorem 10), and that we have concentration around the expectation (Lemma 2). The proof of Theorem 6 down to controlling and propagating these error terms.

###### Lemma 9 .

Suppose that $m\geq 1.5\mathsf{d}_{H}(\lambda)$ and $m\geq M$ for large enough $M$.

A unique root $z^{\star}<0$ such that $\mathbb{E}{\hat{\mathsf{s}}}(z^{\star})=\frac{1}{\lambda}$ exists, and $z^{\star}<-c\lambda$ for small enough $c>0$.

Consider a sufficiently small $O(\lambda)$-neighborhood of $z^{\star}$, $I^{\star}=[z^{\star}-\eta\lambda,z^{\star}+\eta\lambda]$ for small enough $\eta>0$. Then,

- 1.
$\mathbb{E}{\hat{\mathsf{s}}}(z)\gtrsim\frac{1}{\lambda}$ for all $z\in I^{\star}$.
- 2.
$|\mathbb{E}{\hat{\mathsf{s}}}(z)-\mathbb{E}{\hat{\mathsf{s}}}(z^{\prime})|
\lesssim|z-z^{\prime}|\frac{1}{\lambda^{2}}$ for $z,z^{\prime}\in I^{\star}$.

###### Proof.

The existence of $z^{\star}$, as well as the bound $\mathbb{E}{\hat{\mathsf{s}}}_{m}(z)\geq\eta\frac{1}{\lambda}$, can be shown by an argument similar to Lemma 8; we omit the technical details. Specifically, one can show that for small enough $c$,
$\mathbb{E}{\hat{\mathsf{s}}}(-c\lambda)\geq(1+c)\frac{1}{\lambda}$. Since $\mathbb{E}{\mathsf{s}}(\cdot)$ is continuous increasing, and $\mathbb{E}{\hat{\mathsf{s}}}(-\infty)=0$, a root $\mathbb{E}{\hat{\mathsf{s}}}(z^{\star})=\frac{1}{\lambda}$ exists and $z^{\star}<-c\lambda$.

Next, we prove Item 2 above. We have

$$ $\displaystyle|\mathbb{E}{\hat{\mathsf{s}}}_{m}(z)-\mathbb{E}{\hat{\mathsf{s}}} _{m}(z^{\prime})|=\mathbb{E}\frac{1}{m}\mathop{\mathrm{tr}}[(S_{m}HS_{m}^{\top }-zI)^{-1}(S_{m}HS_{m}^{\top}-z^{\prime}I)^{-1}]|z-z^{\prime}|\leq\frac{|z-z^{ \prime}|}{|zz^{\prime}|}\lesssim\frac{1}{\lambda^{2}}|z-z^{\prime}|\,,$ $$

since $|z|,|z^{\prime}|\gtrsim\lambda$, as $z,z^{\prime}\in I_{\eta}$.

Item 1 follows immediately from Item 2.
∎

We need to establish a result of the following kind: if $z<0$ is such that $\mathbb{E}{\hat{\mathsf{s}}}(z)\approx\frac{1}{\lambda}$, then necessarily $z\approx z^{\star}$. This necessitates a lower bound (rather than an upper bound, such as given in Lemma 9, Item 2) on the derivative of $\mathbb{E}{\hat{\mathsf{s}}}(z)$ around $z^{\star}$.

###### Lemma 10 .

Suppose that $m\geq 1.5\mathsf{d}_{H}(\lambda)$ and $m\geq M$ for large enough $M$.

For a sufficiently small $O(\lambda)$-neighborhood of $z^{\star}$, $z^{\star}\in I^{\star}$,

$$ $|\mathbb{E}{\hat{\mathsf{s}}}(z)-\mathbb{E}{\hat{\mathsf{s}}}(z^{\prime})| \gtrsim\frac{1}{\lambda^{2}}|z-z^{\prime}|-O(\frac{1}{\lambda\sqrt{m}})\,.$ (41) $$

###### Proof.

By Theorem 10, having established that $\mathbb{E}{\hat{\mathsf{s}}}(z)\gtrsim\frac{1}{\lambda}$ for all $z\in I^{\star}$, uniformly for all $z\in I^{\star}$,

$$ $\mathbb{E}{\hat{\mathsf{s}}}(z)=-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{H}( 1/\mathbb{E}{\hat{\mathsf{s}}}(z))\right)+O(\frac{1}{\lambda\sqrt{m}})\,.$ $$

uniformly for all $z\in I^{\star}$. Dividing by $\mathbb{E}{\hat{\mathsf{s}}}(z)$ and multiplying by $z$,

$$ $z=\Psi(\mathbb{E}{\hat{\mathsf{s}}}(z))+O(\frac{\lambda}{\sqrt{m}})\,,$ (42) $$

where $\Psi(s)=-\frac{1}{s}(1-\frac{1}{m}\mathsf{d}_{H}(1/s))$ is the inverse function of ${\mathsf{s}}(z)$ (see e.g. (38)).
Note that for $z,z^{\prime}\in I^{\star}$, $|\mathbb{E}{\hat{\mathsf{s}}}(z)-\mathbb{E}{\hat{\mathsf{s}}}(z^{\prime})|
\lesssim|z-z^{\prime}|\frac{1}{\lambda^{2}}=O(\frac{1}{\lambda})$. Thus, to conclude the proof of the lemma, it suffices to give an upper bound on the derivative (or Lipschitz constant) of $\Psi$ in a small $O(1/\lambda)$ neighborhood of $\mathbb{E}{\hat{\mathsf{s}}}(z^{\star})=\frac{1}{\lambda}$.

Indeed, if $s=\frac{1}{\lambda}(1+\varepsilon)$ for small $\varepsilon=O(1)$, then by Lemma 4, $\frac{1}{1+|\varepsilon|}\mathsf{d}_{H}(\lambda)\leq\mathsf{d}_{H}(1/s)\leq
\frac{1}{1-|\varepsilon|}\mathsf{d}_{H}(\lambda)$. Since $\frac{1}{m}\mathsf{d}_{H}(\lambda)\leq(1.5)^{-1}=O(1)$, $|\frac{1}{m}\mathsf{d}_{H}(1/s)-\frac{1}{m}\mathsf{d}_{H}(1/\lambda)|=O(\varepsilon)$. Thus, for $s_{1},s_{2}$ in a sufficiently small $O(\frac{1}{\lambda})$-neighborhood of $\frac{1}{\lambda}$, $|\frac{1}{m}\mathsf{d}_{H}(1/s_{1})-\frac{1}{m}\mathsf{d}_{H}(1/s_{2})|=O(
\lambda|s_{1}-s_{2}|)$, and so $|\Psi(s_{1})-\Psi(s_{2})|=O(\lambda^{2}|s_{1}-s_{2}|)$.
Using (42), we deduce that for $z_{1},z_{2}$ in a small enough $O(\lambda)$ neighborhood of $z^{\star}$,

$$ $\displaystyle|z_{1}-z_{2}|$ $\displaystyle=\left|\Psi(\mathbb{E}{\hat{\mathsf{s}}}(z_{1}))-\Psi(\mathbb{E}{ \hat{\mathsf{s}}}(z_{2}))+O(\frac{\lambda}{\sqrt{m}})\right|$ $\displaystyle\lesssim\lambda^{2}|\mathbb{E}{\hat{\mathsf{s}}}(z_{1})-\mathbb{E }{\hat{\mathsf{s}}}(z_{2})|+O(\frac{\lambda}{\sqrt{m}})\,.$ $$

The desired estimate (41) follows dividing by $\lambda^{2}$.
∎

We are ready to conclude the proof of Theorem 6.
Let $\varepsilon>0$, $\delta\in(0,1)$ be, respectively, the precision on confidence parameters.
Suppose that $\varepsilon=O(1)$ is small, such that $[z^{\star}-2\lambda\varepsilon,z^{\star}+2\lambda\varepsilon]\subseteq I^{\star}$, therefore setting

$$ $\displaystyle z^{\star}_{\pm}=z^{\star}\pm\lambda\varepsilon$ (43) $$

we have $z^{\star}_{\pm}\in I^{\star}$. By
Lemma 10, Eq. (41), provided that $m\gtrsim 1/\varepsilon^{2}$ is large enough,

$$ $\displaystyle\mathbb{E}{\hat{\mathsf{s}}}(z_{-}^{\star})$ $\displaystyle\leq\frac{1}{\lambda}-C_{1}\varepsilon\frac{1}{\lambda}+O(\frac{1 }{\lambda\sqrt{m}})\leq\frac{1}{\lambda}-C_{2}\varepsilon\frac{1}{\lambda}\,,$ $\displaystyle\mathbb{E}{\hat{\mathsf{s}}}(z_{+}^{\star})$ $\displaystyle\geq\frac{1}{\lambda}+C_{1}\varepsilon\frac{1}{\lambda}-O(\frac{1 }{\lambda\sqrt{m}})\geq\frac{1}{\lambda}+C_{2}\varepsilon\frac{1}{\lambda}\,.$ $$

By Lemma 2, if $m\gtrsim\frac{1}{\varepsilon^{2}}\log(1/\delta)$ is large enough, then w.p. $1-\delta$, $|{\hat{\mathsf{s}}}(z_{-}^{\star})-\mathbb{E}{\hat{\mathsf{s}}}(z_{-}^{\star})
|,|{\hat{\mathsf{s}}}(z_{+}^{\star})-\mathbb{E}{\hat{\mathsf{s}}}(z_{+}^{\star
})|\leq C_{2}\frac{1}{\lambda}\varepsilon/2$. On this event,

$$ $\displaystyle{\hat{\mathsf{s}}}(z_{-}^{\star})\leq\frac{1}{\lambda}-\frac{1}{2 }C_{2}\varepsilon\frac{1}{\lambda}<\frac{1}{\lambda}\,,\qquad{\hat{\mathsf{s}} }(z_{+}^{\star})\geq\frac{1}{\lambda}+\frac{1}{2}C_{2}\varepsilon\frac{1}{ \lambda}>\frac{1}{\lambda}\,.$ $$

In particular, the output $\hat{\lambda}$ of Algorithm 2, being the solution ${\hat{\mathsf{s}}}(-\hat{\lambda})=1/\lambda$, satisfies $z^{\star}_{-}\leq-\hat{\lambda}\leq z^{\star}_{+}$, so $|\hat{\lambda}-(-z^{\star})|\leq\varepsilon$. Finally, by Eq. (42),

$$ $z^{\star}=\Psi(1/\lambda)+O(\frac{\lambda}{\sqrt{m}})=-\tilde{\lambda}+O(\frac {\lambda}{\sqrt{m}})=-\tilde{\lambda}+O(\lambda\varepsilon)$ $$

whenever $m\gtrsim 1/\varepsilon^{2}$ large enough.
Since $\tilde{\lambda}\geq\lambda(1-(1.5)^{-1})=\frac{1}{3}\lambda=\Omega(\lambda)$, this implies

$$ $\frac{|\hat{\lambda}-\tilde{\lambda}|}{\tilde{\lambda}}=O(\varepsilon)\,.$ $$

Thus, Theorem 6 is proved.

∎

In the sequel, the following bound on $\mathbb{E}|\hat{\lambda}-\tilde{\lambda}|$ will be useful.

###### Lemma 11 .

Suppose that $m\geq 1.5\mathsf{d}_{H}(\lambda)$ and $m\geq M$ for large enough $M$. For any $p>0$,

$$ $\displaystyle\mathbb{E}|\hat{\lambda}-\tilde{\lambda}|^{p}\leq C_{p}\lambda^{p }m^{-p/2}\,.$ (44) $$

###### Proof.

Let $t_{0}>0$ be a small enough constant, such that $[z^{\star}-\lambda t_{0},z^{\star}+\lambda t_{0}]\subseteq I^{\star}$. By the above, $\Pr(|\hat{\lambda}-\tilde{\lambda}|\geq t_{0}\lambda)=\Pr(-\hat{\lambda}\notin
[z^{\star}-\lambda t_{0},z^{\star}+\lambda t_{0}])\leq Ce^{-cm}$. Set $z^{\star}_{\pm}(t)=z^{\star}\pm\lambda t$, so that by the above, $\Pr({\hat{\mathsf{s}}}_{m}(z^{\star}_{+}(t))<1/\lambda),\Pr({\hat{\mathsf{s}}}
_{m}(z^{\star}_{-}(t))>1/\lambda)\leq Ce^{-cmt^{2}}$ for all $t\leq t_{0}$. We have

$$ $\Pr(|\hat{\lambda}-\tilde{\lambda}|\geq t\lambda)=\Pr(-\hat{\lambda}\notin[z^{ \star}_{-}(t),z^{\star}_{+}(t)])=\Pr(\{{\hat{\mathsf{s}}}_{m}(z^{\star}_{+}(t) )<1/\lambda\}\cup\{{\hat{\mathsf{s}}}_{m}(z^{\star}_{-}(t))>1/\lambda\})\leq 2 Ce^{-cmt^{2}}\,.$ $$

Thus,

$$ $\displaystyle\mathbb{E}|\hat{\lambda}-\tilde{\lambda}|^{p}=\mathbb{E}[|\hat{ \lambda}-\tilde{\lambda}|^{p}\mathds{1}_{|\hat{\lambda}-\tilde{\lambda}|\leq \lambda t_{0}}]+\mathbb{E}[|\hat{\lambda}-\tilde{\lambda}|^{p}\mathds{1}_{| \hat{\lambda}-\tilde{\lambda}|>\lambda t_{0}}]\leq\mathbb{E}[|\hat{\lambda}- \tilde{\lambda}|\mathds{1}_{|\hat{\lambda}-\tilde{\lambda}|\leq\lambda t_{0}}] +C\lambda^{p}e^{-cm}\,,$ $$

and

$$ $\displaystyle\mathbb{E}[|\hat{\lambda}-\tilde{\lambda}|^{p}\mathds{1}_{|\hat{ \lambda}-\tilde{\lambda}|\leq\lambda t_{0}}]=\int_{0}^{t_{0}}\lambda^{p}pt^{p- 1}\Pr(|\hat{\lambda}-\tilde{\lambda}|\geq\lambda t)dt\lesssim p\lambda^{p}\int _{0}^{t_{0}}t^{p-1}e^{-cmt^{2}}dt\leq C_{p}\lambda^{p}m^{-p/2}\,.$ $$

And so, $\mathbb{E}|\hat{\lambda}-\tilde{\lambda}|^{p}\lesssim\lambda\frac{1}{\sqrt{m}}
+\lambda e^{-cm}=O(\lambda\frac{1}{\sqrt{m}})$.
∎

###### .

###### Lemma 9 .

###### Proof.

###### Lemma 10 .

###### Proof.

###### Lemma 11 .

###### Proof.

### Appendix E Proof of Theorem 9

###### .

By definition, if $\tilde{H}\preceq H$ (in PSD order) then $S\tilde{H}S^{\top}\preceq SHS^{\top}$, in particular $\lambda_{\min}(SHS^{\top})\geq\lambda_{\min}(S\tilde{H}S^{\top})$. Choose $\tilde{H}=H(\mu H+I)^{-1}$, noting that $\|\tilde{H}\|\leq 1/\mu$ can be bounded irrespective of $\|H\|$.

We now lower bound $\lambda_{\min}(S\tilde{H}S)$ using a standard net argument, cf. \citepvershynin2018high.
For any $\varepsilon>0$, let $\mathcal{N}_{\varepsilon}$ be an $\varepsilon$-net of $\mathcal{S}^{m-1}$ of minimum size. One can show that $|\mathcal{N}_{\varepsilon}|\leq\left(1+\frac{2}{\varepsilon}\right)^{m}$.
We have

$$ $\displaystyle\lambda_{\min}(S\tilde{H}S^{\top})=\min_{\|x\|=1}x^{\top}S\tilde{ H}S^{\top}x\geq\min_{x\in\mathcal{N}_{\varepsilon}}x^{\top}S\tilde{H}S^{\top}x -2\|S\tilde{H}S^{\top}\|\varepsilon\,.$ (45) $$

Note that for fixed $x\in\mathcal{S}^{m-1}$, $y=S^{\top}x$ has independent, mean zero, sub-Gaussian entries with $\max_{1\leq i\leq d}\|y_{i}\|_{\psi_{2}}=O(1)$ and variance $1/m$.

We first give a high-probability upper bound on $\|S\tilde{H}S^{\top}\|$. Let $\mathcal{N}_{1/4}$ a $1/4$-net of minimum size, so that

$$ $\displaystyle\|S\tilde{H}S^{\top}\|=\max_{\|x\|=1}x^{\top}S\tilde{H}Sx\leq\max _{x\in\mathcal{N}_{1/4}}x^{\top}S\tilde{H}Sx+2\cdot\|S\tilde{H}S^{\top}\|\cdot \frac{1}{4}\,,$ $$

hence $\|S\tilde{H}S^{\top}\|\leq 2\max_{x\in\mathcal{N}_{1/4}}x^{\top}S_{m}\tilde{H}
S_{m}x$.
By the Hanson-Wright inequality, Lemma 95 (see also [vershynin2018high, Theorem 6.2.1]),

$$ $\displaystyle\Pr(\|S\tilde{H}S^{\top}\|\geq 2m^{-1}\mathop{\mathrm{tr}}(\tilde {H})+t)$ $\displaystyle\leq|\mathcal{N}_{1/4}|\max_{x\in\mathcal{N}_{1/4}}\Pr\left(x^{ \top}S\tilde{H}Sx\geq m^{-1}\mathop{\mathrm{tr}}(\tilde{H})+t/2\right)$ $\displaystyle\leq 9^{m}2\exp\left(-c\min\{\frac{m^{2}t^{2}}{\|\tilde{H}\|_{F}^ {2}},\frac{mt}{\|\tilde{H}\|}\}\right)\,.$ $$

Note that

$$ $\displaystyle\mathop{\mathrm{tr}}(\tilde{H})=\frac{1}{\mu}{\mathsf{d}_{H}(1/ \mu)},\qquad\|\tilde{H}\|\leq 1/\mu,\qquad\|\tilde{H}\|_{F}\leq\sqrt{\|\tilde{ H}\|\mathop{\mathrm{tr}}(\tilde{H})}\leq\frac{1}{\mu}\sqrt{\mathsf{d}_{H}(1/ \mu)}\,.$ $$

Set $t=C_{1}\max\{\frac{1}{\sqrt{m}}\|\tilde{H}\|_{F},\|\tilde{H}\|\}$ large enough. For such choice, w.p. $1-Ce^{-cm}$,

$$ $\displaystyle\|S\tilde{H}S\|$ $\displaystyle\leq\frac{2}{m}\mathop{\mathrm{tr}}(\tilde{H})+C_{1}\max\{\frac{1 }{\sqrt{m}}\|\tilde{H}\|_{F},\|\tilde{H}\|\}$ $\displaystyle\lesssim\frac{1}{\mu}\frac{{\mathsf{d}_{H}(1/\mu)}}{m}+\frac{1}{ \mu}\max\{1,\sqrt{\frac{d_{1/\mu}}{m}}\}$ (46) $$

Next, we give a high-probability lower bound on $\min_{x\in\mathcal{N}_{\varepsilon}}x^{\top}S\tilde{H}S^{\top}x$ in (45). Again by Hanson-Wright,

$$ $\displaystyle\Pr(\min_{x\in\mathcal{N}_{\varepsilon}}x^{\top}S\tilde{H}S^{\top }x\leq m^{-1}\mathop{\mathrm{tr}}(\tilde{H})-t)\leq\left(1+\frac{2}{ \varepsilon}\right)^{m}2\exp\left(-c\min\{\frac{m^{2}t^{2}}{\|\tilde{H}\|_{F}^ {2}},\frac{mt}{\|\tilde{H}\|}\}\right)\,,$ $$

so that w.p. $1-Ce^{-cm}$, for small $\varepsilon>0$,

$$ $\displaystyle\min_{x\in\mathcal{N}_{\varepsilon}}x^{\top}S\tilde{H}S^{\top}x \geq\frac{1}{\mu}\frac{{\mathsf{d}_{H}(1/\mu)}(H)}{m}-C_{2}\log(1+1/ \varepsilon)\frac{1}{\mu}\max\{1,\sqrt{\frac{{\mathsf{d}_{H}(1/\mu)}(H)}{m}}\}$ (47) $$

Inserting (46)-(47) into (45) implies that w.p. $1-Ce^{-cm}$,

$$ $\lambda_{\min}(SHS^{\top})\geq(1-O(\varepsilon))\frac{1}{\mu}\frac{\mathsf{d}_ {H}(1/\mu)}{m}-O(\varepsilon+\log(1+1/\varepsilon))\frac{1}{\mu}\max\{1,\sqrt{ \frac{{\mathsf{d}_{H}(1/\mu)}(H)}{m}}\}.$ (48) $$

By choosing $\varepsilon>0$ a sufficiently small numerical constant, we can deduce the following: there are numerical constants $c,C,c_{1},c_{2}>0$ such that if $\mathsf{d}_{H}(1/\mu)/m>1/c_{1}$, then
w.p. $1-Ce^{-cm}$, $\lambda_{\min}(\tilde{H})\geq c_{2}\frac{1}{\mu}\frac{\mathsf{d}_{H}(1/\mu)}{m}$. To recover the claimed result, use $\mu=1/\eta$.

∎

###### .

### Appendix F Proof of Lemma 2

The following argument is standard (cf. \citepbai2010spectral) and presented for completeness.

###### .

Let $r_{1},\ldots,r_{m}\in\mathbb{R}^{d}$ be the rows of $S$; that is, $S^{\top}=[r_{1}|\ldots|r_{m}]$. Let

$$ $\hat{\Sigma}_{m}=H^{1/2}S^{\top}SH^{1/2}=\frac{1}{m}\sum_{i=1}^{n}H^{1/2}r_{i} r_{i}^{\top}H^{1/2}\,.$ (49) $$

This is a sample covariance matrix, corresponding to $m$ samples with population covariance $\mathbb{E}[\hat{\Sigma}_{m}]=H$. Recall that $\hat{\Sigma}_{m}\in\mathbb{R}^{d\times d}$ and $S_{m}HS_{m}^{\top}\in\mathbb{R}^{m\times m}$ have the same non-zero eigenvalues. Denote

$$ $\displaystyle P_{m}(z)=(\hat{\Sigma}_{m}+zI)^{-1}\,,$ (50) $$

so that

$$ $\displaystyle{\hat{\mathsf{s}}}(z)=m^{-1}\mathop{\mathrm{tr}}(P_{m}(z))+m^{-1} (d-m)\frac{1}{z}\,.$ (51) $$

Central to the proof is the following leave-one-out decomposition:

$$ $\displaystyle P_{m}(z)$ $\displaystyle=P_{m}^{-k}(z)-\frac{\frac{1}{m}P_{m}^{-k}(z)H^{1/2}r_{k}r_{k}^{ \top}H^{1/2}P_{m}^{-k}(z)}{1+\frac{1}{m}r_{k}^{\top}H^{1/2}P_{m}^{-k}(z)H^{1/2 }r_{k}},$ (52) $\displaystyle P_{m}^{-k}(z)$ $\displaystyle=(\hat{\Sigma}_{m}-\frac{1}{m}H^{1/2}r_{k}r_{k}^{\top}H^{1/2}+zI) ^{-1},\qquad k=1,\ldots,m\,.$ (53) $$

The above can be readily shown by the Sherman-Morrison lemma (Lemma 92).

We shall now prove Lemma 2 by a standard martingale concentration argument. Let $\mathcal{F}_{k}$ be the $\sigma$-algebra generated by $r_{1},\ldots,r_{k}$, $k=0,1,\ldots,m$, and $\mathbb{E}_{\leq k}[\cdot]:=\mathbb{E}[\cdot|\mathcal{F}_{k}]$. Decompose into a sum of martingale differences,

$$ $\displaystyle{\hat{\mathsf{s}}}(z)-\mathbb{E}[{\hat{\mathsf{s}}}(z)]$ $\displaystyle=m^{-1}\sum_{k=1}^{m}(\mathbb{E}_{\leq k}-\mathbb{E}_{\leq k-1})[ \mathop{\mathrm{tr}}P_{m}(z)]$ $\displaystyle=m^{-1}\sum_{k=1}^{m}(\mathbb{E}_{\leq k}-\mathbb{E}_{\leq k-1})[ D_{m,k}(z)]\,,$ (54) $$

where, using (52),

$$ $\displaystyle D_{m,k}(z)=\frac{\frac{1}{m}r_{k}^{\top}H^{1/2}(P_{m}^{-k}(z))^{ 2}H^{1/2}r_{k}}{1+\frac{1}{m}r_{k}^{\top}H^{1/2}P_{m}^{-k}(z)H^{1/2}r_{k}}\,.$ (55) $$

Note that $0\preceq P_{m}^{-k}(z)\preceq\frac{1}{-z}I$, hence

$$ $\displaystyle|D_{m,k}(z)|\leq\frac{1}{-z}\cdot\frac{\frac{1}{m}r_{k}^{\top}H^{ 1/2}(P_{m}^{-k}(z))H^{1/2}r_{k}}{1+\frac{1}{m}r_{k}^{\top}H^{1/2}P_{m}^{-k}(z) H^{1/2}r_{k}}\leq\frac{1}{-z}\,.$ $$

Thus, (54) is a sum of bounded martingale differences $|(\mathbb{E}_{\leq k}-\mathbb{E}_{\leq k-1})[D_{m,k}(z)]|\leq\frac{2}{-z}$. To conclude, use the Azuma-Hoeffding inequality (Lemma 22).

Remark: The above calculation (regarding the boundedness of $D_{m,k}(z)$) is, essentially, the proof of the well-known low-rank perturbation bound for resolvents, Lemma 21, which we shall use later.

∎

###### .

### Appendix G Proof of Theorem 10

###### .

We implement a well-known computation technique in random matrix theory (cf. [bai2010spectral]), while carefully keeping track of the error terms.

We rely explicitly on properties of the Gaussian distribution. Denote the eigendecomposition $H=\sum_{\ell=1}^{d}\tau_{\ell}v_{\ell}v_{\ell}^{\top}$ and

$$ $\displaystyle H_{\setminus\ell}=H-\tau_{\ell}v_{\ell}v_{\ell}^{\top}\,,\qquad s _{\ell}=Sv_{\ell}\,.$ $$

for $\ell=1,\ldots,d$. Note that $s_{1},\ldots,s_{d}\sim\mathcal{N}(0,m^{-1}I_{m})$ are independent (since $v_{1},\ldots,v_{d}$ are orthogonal).

Recall Eq. (51), so that

$$ $\displaystyle{\hat{\mathsf{s}}}(z)$ $\displaystyle=-\frac{1}{m}(d-m)\frac{1}{z}+\frac{1}{m}\mathop{\mathrm{tr}}(H^{ 1/2}S^{\top}SH^{1/2}-zI)^{-1}$ $\displaystyle=-\frac{1}{m}(d-m)\frac{1}{z}+\frac{1}{m}\sum_{\ell=1}^{d}v_{\ell }^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}v_{\ell}\,.$ (56) $$

We now compute an expression for $v_{\ell}^{\top}(H^{1/2}S_{m}^{\top}S_{m}H^{1/2}-zI)^{-1}v_{\ell}$, the $\ell$-th diagonal element of the resolvent $(H^{1/2}S_{m}^{\top}S_{m}H^{1/2}-zI)^{-1}$ written in the population covariance eigenbasis $V=[v_{1}|\dots|v_{d}]^{\top}$.

Upon a coordinate permutation (where $\ell$ becomes the first), we have

$$ $\displaystyle V^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)V=\begin{bmatrix}\tau_{\ell} \|s_{\ell}\|^{2}-z&\sqrt{\tau_{\ell}}s_{\ell}SH^{1/2}_{\setminus\ell}\\ \sqrt{\tau_{\ell}}H_{\setminus\ell}^{1/2}S^{\top}s_{\ell}&H_{\setminus\ell}^{1 /2}S^{\top}SH_{\setminus\ell}^{1/2}-zI\end{bmatrix}\,.$ $$

By the block matrix inverse formula (Lemma 93),

$$ $\displaystyle v_{\ell}^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}v_{\ell}$ $\displaystyle=(V^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)V)^{-1}_{\ell,\ell}$ $\displaystyle=\left(\tau_{\ell}\|s_{\ell}\|^{2}-z-{\tau_{\ell}}s_{\ell}SH^{1/2 }_{\setminus\ell}(H_{\setminus\ell}^{1/2}S^{\top}SH_{\setminus\ell}^{1/2}-zI)^ {-1}H_{\setminus\ell}^{1/2}S^{\top}s_{\ell}\right)^{-1}$ $\displaystyle\overset{(\star)}{=}\left(\tau_{\ell}\|s_{\ell}\|^{2}-z-{\tau_{ \ell}}s_{\ell}(SH_{\setminus\ell}S^{\top}-zI)^{-1}SH_{\setminus\ell}S^{\top}s_ {\ell}\right)^{-1}$ $\displaystyle=-\frac{1}{z}\left(1+\tau_{\ell}s_{\ell}(SH_{\setminus\ell}S^{ \top}-zI)^{-1}s_{\ell}\right)^{-1}\,,$ (57) $$

where in $(\star)$ we used the identity $X^{\top}f(XX^{\top})X=f(X^{\top}X)X^{\top}X$ which holds for any matrix $X$ and analytic function $f(\cdot)$.

We now wish to estimate the expectation of (57). Denote

$$ $\displaystyle D_{\ell}$ $\displaystyle=s_{\ell}(SH_{\setminus\ell}S^{\top}-zI)^{-1}s_{\ell}-\mathbb{E}{ \hat{\mathsf{s}}}(z)$ $\displaystyle=D_{\ell,1}+D_{\ell,2}+D_{\ell,3}$ $$

where

$$ $\displaystyle D_{\ell,1}$ $\displaystyle=s_{\ell}(SH_{\setminus\ell}S^{\top}-zI)^{-1}s_{\ell}-\frac{1}{m} \mathop{\mathrm{tr}}(SH_{\setminus\ell}S^{\top}-zI)^{-1}\,,$ $\displaystyle D_{\ell,2}$ $\displaystyle=\frac{1}{m}\mathop{\mathrm{tr}}(SH_{\setminus\ell}S^{\top}-zI)^{ -1}-\frac{1}{m}\mathop{\mathrm{tr}}(SHS^{\top}-zI)^{-1}\,,$ $\displaystyle D_{\ell,3}$ $\displaystyle=\frac{1}{m}\mathop{\mathrm{tr}}(SHS^{\top}-zI)^{-1}-\frac{1}{m} \mathbb{E}\mathop{\mathrm{tr}}(SHS^{\top}-zI)^{-1}\,.$ $$

We have

$$ $\displaystyle v_{\ell}^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}v_{\ell}=-\frac{ 1}{z}\left(1+\tau_{\ell}\mathbb{E}{\hat{\mathsf{s}}}(z)\right)^{-1}+\frac{1}{z }\frac{\tau_{\ell}D_{\ell}}{(1+\tau_{\ell}s_{\ell}(SH_{\setminus\ell}S^{\top}- zI)^{-1}s_{\ell})(1+\tau_{\ell}\mathbb{E}{\mathsf{s}}(z))}\,,$ $$

so that using (56),

$$ $\displaystyle\mathbb{E}{\mathsf{s}}(z)$ $\displaystyle=\frac{1}{m}(d-m)\frac{1}{z}+\frac{1}{m}\sum_{\ell=1}^{d}\mathbb{ E}v_{\ell}^{\top}(H^{1/2}S^{\top}SH^{1/2}-zI)^{-1}v_{\ell}$ $\displaystyle=\frac{1}{m}(d-m)\frac{1}{z}-\frac{1}{m}\sum_{\ell=1}^{d}\frac{1} {z}\left(1+\tau_{\ell}\mathbb{E}{\hat{\mathsf{s}}}(z)\right)^{-1}+e_{m}(z)$ $\displaystyle=-\frac{1}{z}+\frac{1}{z}\frac{1}{m}\mathop{\mathrm{tr}}\left[I-( I+\mathbb{E}{\hat{\mathsf{s}}}(z)H)^{-1}\right]+e_{m}(z)$ $\displaystyle=-\frac{1}{z}\left(1-\frac{1}{m}\mathsf{d}_{H}(1/\mathbb{E}{\hat{ \mathsf{s}}}(z))\right)+e_{m}(z)\,,$ (58) $$

where

$$ $\displaystyle e_{m}(z)$ $\displaystyle=\frac{1}{zm}\sum_{\ell=1}^{m}\mathbb{E}\left[\frac{\tau_{\ell}D_ {\ell}}{(1+\tau_{\ell}s_{\ell}(SH_{\setminus\ell}S^{\top}-zI)^{-1}s_{\ell})(1+ \tau_{\ell}\mathbb{E}{\mathsf{s}}(z))}\right]\,.$ $$

We may now bound

$$ $\displaystyle|e_{m}(z)|$ $\displaystyle\leq\frac{1}{m|z|}\sum_{\ell=1}^{d}\frac{\tau_{\ell}}{1+\tau_{ \ell}\mathbb{E}{\hat{\mathsf{s}}}(z)}\max_{1\leq\ell\leq d}\mathbb{E}|D_{\ell}|$ $\displaystyle=\frac{1}{|z|\mathbb{E}{\hat{\mathsf{s}}}(z)}\frac{1}{m}\mathsf{d }_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z))\max_{1\leq\ell\leq d}\mathbb{E}|D_{ \ell}|\,.$ (59) $$

We now decompose $\mathbb{E}|D_{\ell}|\leq\mathbb{E}|D_{\ell,1}|+\mathbb{E}|D_{\ell,2}|+\mathbb{
E}|D_{\ell,3}|$. By Lemma 26 (Item 1),

$$ $\displaystyle\mathbb{E}|D_{\ell,1}|\leq\sqrt{\mathbb{E}|D_{1,\ell}|^{2}} \lesssim\frac{1}{m}\mathbb{E}\|(S_{m}H_{\setminus\ell}S_{m}-zI)^{-1}\|_{F}\leq \frac{1}{|z|\sqrt{m}}\,.$ $$

By the low rank resolvent perturbation lemma, Lemma 21, almost surely

$$ $\displaystyle|D_{\ell,2}|\lesssim\frac{1}{|z|m}\,.$ $$

Finally, by Lemma 2,

$$ $\displaystyle\mathbb{E}|D_{\ell,3}|=\mathbb{E}|{\hat{\mathsf{s}}}(z)-\mathbb{E }{\hat{\mathsf{s}}}_{m}(z)|\lesssim\frac{1}{|z|\sqrt{m}}\,.$ $$

Thus,

$$ $|e_{m}(z)|=O\left(\frac{1}{|z|\mathbb{E}{\hat{\mathsf{s}}}(z)}\frac{1}{m} \mathsf{d}_{H}(1/\mathbb{E}{\hat{\mathsf{s}}}(z))\frac{1}{|z|\sqrt{m}}\right)\,,$ (60) $$

and so the proof is concluded.

∎

###### .

### Appendix H Proof of Theorem 7

Let $\hat{\lambda}$ the output of Algorithm 2, and denote

$$ $\displaystyle\mathcal{E}=(H+\lambda I)^{1/2}\hat{W}(H+\lambda I)^{1/2}-I, \qquad\hat{W}:=S^{\top}(SHS^{\top}+\hat{\lambda}I)^{-1}S\,.$ (61) $$

We have $\bar{\mathcal{E}}=\frac{1}{q}\sum_{\ell=1}^{q}\mathcal{E}^{(\ell)}$, where $\mathcal{E}^{(1)},\ldots,\mathcal{E}^{(q)}\overset{\mathrm{i.i.d.}}{\sim}
\mathcal{E}$.

The proof proceeds in two parts. First we show that the expectation $\mathbb{E}[\mathcal{E}]$ is small; then, we show that $\bar{\mathcal{E}}$ concentrates around $\mathbb{E}[\bar{\mathcal{E}}]=\mathbb{E}[\mathcal{E}]$ using matrix concentration inequalities.

The following lemma, proven in Section I, bounds $\mathbb{E}[\mathcal{E}]$.

###### Lemma 12 .

Suppose that $S$ has i.i.d. Gaussian entries. Assume that $m\geq 1.5\mathsf{d}_{H}(\lambda)$. We have

$$ $\|\mathbb{E}[\mathcal{E}]\|=O(\frac{1}{\sqrt{m}})\,.$ $$

To establish concentration, we use the matrix Bernstein inequality [tropp2015introduction, Theorem 6.6.1]. We cite it here.

###### Lemma 13 .

Let $X_{1},\ldots,X_{n}\in\mathbb{R}^{d\times d}$ be independent Hermitian random matrices.
Assume that $\mathbb{E}[X_{n}]=0$, $\|X_{n}\|\leq L$, and denote

$$ $\displaystyle\nu^{2}=\left\|\sum_{n=1}^{N}\mathbb{E}[X_{n}^{2}]\right\|\,.$ (62) $$

Then for all $t\geq 0$,

$$ $\displaystyle\Pr\left(\left\|\sum_{n=1}^{N}X_{n}\right\|\geq t\right)\leq d \exp\left(\frac{-t^{2}/2}{\nu^{2}+Lt/3}\right)\,.$ (63) $$

We shall apply Lemma 13 with a truncated version of $\mathcal{E}$, $\tilde{\mathcal{E}}_{L}=\mathcal{E}\mathds{1}_{\|\mathcal{E}\|\leq L}-\mathbb{
E}[\mathcal{E}\mathds{1}_{\|\mathcal{E}\|\leq L}]$ for a suitably chosen $L$.
To this end, we give the following high-probability bound on $\|\Phi\|$.

###### Lemma 14 .

Suppose that $S$ has i.i.d. Gaussian entries. For $\delta\in(0,1)$ and $q\geq 1$, set

$$ $\displaystyle L_{\delta,q}=C\frac{d}{m}+C\frac{\log(q/\delta)}{m}$ (64) $$

for large enough $C>0$. Then,

- 1.
W.p. $1-\delta/q$, $\|\mathcal{E}\|\leq L_{\delta,q}$.
- 2.
$\mathbb{E}[\|\mathcal{E}\|\mathds{1}_{\|\mathcal{E}\|\geq L_{\delta,q}}]=O(
\frac{1}{m}e^{-d})$.
- 3.
$\mathbb{E}[\|\mathcal{E}\|^{2}\mathds{1}_{\|\mathcal{E}\|\geq L_{\delta,q}}]=O
(\frac{d}{m^{2}}e^{-d})$.

We prove Lemma 14 in Section J.

The last component needed to apply matrix Bernstein is the variance proxy (62). Note that we always have $\nu^{2}\leq L^{2}N$, so that

$$ $\displaystyle\Pr\left(\left\|\frac{1}{N}\sum_{n=1}^{N}X_{n}\right\|\geq \varepsilon\right)\leq d\exp(-cN\min\left\{\left(\frac{\varepsilon}{L}\right)^ {2},\frac{\varepsilon}{L}\right\})\,.$ (65) $$

However, if the variance $\mathbb{E}[X^{2}]\ll L^{2}$, one may obtain substantially sharper bounds for small $\varepsilon\ll L$.
This is the case in our setting. The following is proved in Section K.

###### Lemma 15 .

Suppose that $S$ has i.i.d. Gaussian entries. Assume that $m\geq 1.5\mathsf{d}_{H}(\lambda)$. Then

$$ $\mathbb{E}[\mathcal{E}^{2}]=O(d/m)\,.$ $$

With the above, we are ready to conclude the proof of Theorem 7.

###### Proof.

(Of Theorem 7.)

Consider the truncated matrices $\tilde{\mathcal{E}}^{(1)}_{L_{\delta,q}},\ldots,\tilde{\mathcal{E}}^{(q)}_{L_{
\delta,q}}$ with $L_{\delta,q}$ as in Lemma 14. By Item 1 of Lemma 14 (and union bound over $\ell=1,\ldots,q$), with probability $1-\delta$, we have

$$ $\displaystyle\bar{\mathcal{E}}=\frac{1}{q}\sum_{\ell=1}^{q}\tilde{\mathcal{E}} ^{(\ell)}_{L_{\delta,q}}+\mathbb{E}\left[{\mathcal{E}}\mathds{1}_{\|\mathcal{E }\|\leq L_{\delta,q}}\right]\,.$ (66) $$

Note that

$$ $\displaystyle\left\|\mathbb{E}\left[{\mathcal{E}}\mathds{1}_{\|\mathcal{E}\| \leq L_{\delta,q}}\right]\right\|=\left\|\mathbb{E}[\mathcal{E}]+\mathbb{E} \left[{\mathcal{E}}\mathds{1}_{\|\mathcal{E}\|>L_{\delta,q}}\right]\right\| \leq\|\mathbb{E}[\mathcal{E}]\|+\mathbb{E}\left[\|\mathcal{E}\|\mathds{1}_{\| \mathcal{E}\|>L_{\delta,q}}\right]=O(1/\sqrt{m})\,,$ (67) $$

where the first inequality uses Jensen’s inequality, and the second inequality follows from Lemma 14 Item 2 and Lemma 12. Similarly,

$$ $\displaystyle\left\|\mathbb{E}\left[\tilde{\mathcal{E}}_{L_{\delta,q}}^{2} \right]\right\|\lesssim\left\|\mathbb{E}\left[{\mathcal{E}}^{2}\right]\right\| +\mathbb{E}\left[\|\mathcal{E}\|^{2}\mathds{1}_{\|\mathcal{E}\|>L_{\delta,q}} \right]=O(d/m)\,,$ (68) $$

where we used Lemma 15.

Now, by the matrix Bernstein inequality (Lemma 13), the bound

$$ $\displaystyle\|\bar{\mathcal{E}}\|\leq\varepsilon+O(1/\sqrt{m})$ (69) $$

holds with probability at least

$$ $\displaystyle 1-d\exp\left(-cq\min\{\frac{\varepsilon^{2}}{d/m},\frac{ \varepsilon}{L_{p,q}}\}\right)\,.$ (70) $$

When $q\leq\exp(O(d))$ and $1/\delta\leq\exp(O(d))$ we have $L_{\delta,q}=O(d/m)$, so for $\varepsilon=O(1)$ the smaller term is $\frac{\varepsilon^{2}}{d/m}$. And so, $q\gtrsim\frac{d\log d}{m}\frac{\log(1/\delta)}{\varepsilon^{2}}$ ensures that the probability above is $\geq 1-\delta$.

∎

###### Lemma 12 .

###### Lemma 13 .

###### Lemma 14 .

###### Lemma 15 .

###### Proof.

### Appendix I Proof of Lemma 12

###### .

We use the machinery from Section G. As before, we denote the eigendecomposition of $H$,

$$ $\displaystyle H=V\mathop{\mathrm{diag}}(\tau_{1},\ldots,\tau_{d})V^{\top}, \qquad V=[v_{1}|\cdots|v_{d}]$ $$

and

$$ $\displaystyle H_{\setminus\ell}=H-\tau_{\ell}v_{\ell}v_{\ell}^{\top}\,,\qquad s _{\ell}=Sv_{\ell},\qquad s_{1},\ldots,s_{d}\overset{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,m^{-1}I_{m})\,.$ $$

We shall compute the expectation of $\hat{W}=S^{\top}(SHS^{\top}+I)^{-1}S$ in the eigenbasis $V$, that is, $\mathbb{E}[V^{\top}\hat{W}V]$.

Denote, for brevity,

$$ $\displaystyle T(\hat{\lambda})=(SHS^{\top}+\hat{\lambda}I)^{-1},\qquad T_{ \setminus\ell}(\hat{\lambda})=(SH_{\setminus\ell}S^{\top}+\hat{\lambda}I)^{-1}\,,$ $$

so that by the Sherman-Morrison formula (Lemma 92), writing $SHS^{\top}=SH_{\setminus\ell}S^{\top}+\tau_{\ell}s_{\ell}s_{\ell}^{\top}$,

$$ $\displaystyle T(\hat{\lambda})=T_{\setminus\ell}(\hat{\lambda})-\frac{\tau_{ \ell}T_{\setminus\ell}(\hat{\lambda})s_{\ell}s_{\ell}T_{\setminus\ell}(\hat{ \lambda})}{1+\tau_{\ell}s_{\ell}T_{\setminus\ell}(\hat{\lambda})s_{\ell}}\,.$ (71) $$

Let $j\neq\ell$. Note that $\hat{\lambda},T_{\setminus\ell}(\hat{\lambda})$ do not change under a sign flip $s_{j}\mapsto-s_{j}$. Since $s_{j}$ has a symmetric distribution ($s_{j}\overset{d}{=}-s_{j}$), we deduce that

$$ $\displaystyle\mathbb{E}[v_{\ell}^{\top}\hat{W}v_{j}]=\mathbb{E}[s_{\ell}^{\top }T(\hat{\lambda})s_{j}]=0\,,$ (72) $$

that is, $\mathbb{E}[V^{\top}\hat{W}V]$ is diagonal. Consequently, $\mathbb{E}[V^{\top}\mathcal{E}V]=\mathbb{E}[V^{\top}(H+\lambda I)^{1/2}\hat{W}
(H+\lambda I)^{1/2}V]$ is diagonal as well (since $V$ is an eigenbasis of $H$).

Let us calculate the diagonal elements. By (71),

$$ $\displaystyle v_{\ell}^{\top}\hat{W}v_{\ell}=s_{\ell}^{\top}T(\hat{\lambda})s_ {\ell}=\frac{s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell}}{1+\tau_{ \ell}s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell}}=\frac{1}{(s_{ \ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell})^{-1}+\tau_{\ell}}\,,$ (73) $$

so

$$ $\displaystyle v_{\ell}^{\top}\Phi v_{\ell}$ $\displaystyle=\frac{\lambda+\tau_{\ell}}{(s_{\ell}^{\top}T_{\setminus\ell}( \hat{\lambda})s_{\ell})^{-1}+\tau_{\ell}}-1$ $\displaystyle=\frac{\lambda-(s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_ {\ell})^{-1}}{(s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell})^{-1}+ \tau_{\ell}}$ $\displaystyle=O(|\lambda s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{ \ell}-1|)\,.$ (74) $$

Let $\tilde{\lambda}$ be given by (9), so that by Theorem 6, $\hat{\lambda}$ is close to $\tilde{\lambda}$ with high probability. Let

$$ $\displaystyle\Delta_{1}$ $\displaystyle=\lambda s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell}- \lambda s_{\ell}^{\top}T_{\setminus\ell}(\tilde{\lambda})s_{\ell}\,,$ $\displaystyle\Delta_{2}$ $\displaystyle=\lambda s_{\ell}^{\top}T_{\setminus\ell}(\tilde{\lambda})s_{\ell }-\lambda\frac{1}{m}\mathop{\mathrm{tr}}T_{\setminus\ell}(\tilde{\lambda})\,,$ $\displaystyle\Delta_{3}$ $\displaystyle=\lambda\frac{1}{m}\mathop{\mathrm{tr}}T_{\setminus\ell}(\tilde{ \lambda})-\lambda\frac{1}{m}\mathop{\mathrm{tr}}T(\tilde{\lambda})\,,$ $\displaystyle\Delta_{4}$ $\displaystyle=\lambda\frac{1}{m}\mathop{\mathrm{tr}}T(\tilde{\lambda})-\lambda \frac{1}{m}\mathop{\mathrm{tr}}T(\hat{\lambda})=\lambda\frac{1}{m}\mathop{ \mathrm{tr}}T(\tilde{\lambda})-1\,,$ $$

so that

$$ $\lambda s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{\ell}-1=\Delta_{1}+ \Delta_{2}+\Delta_{3}+\Delta_{4}\,.$ $$

We have

$$ $\displaystyle\mathbb{E}|\Delta_{1}|$ $\displaystyle=\mathbb{E}[\lambda|s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda })T_{\setminus\ell}(\tilde{\lambda})s_{\ell}||\hat{\lambda}-\tilde{\lambda}|]= \mathbb{E}\left[\frac{\lambda}{\hat{\lambda}\tilde{\lambda}}\|s_{\ell}\|^{2}| \hat{\lambda}-\tilde{\lambda}|\right]\lesssim\frac{1}{\lambda}(\mathbb{E}\|s_{ \ell}\|^{4}\|)^{1/2}(\mathbb{E}|\hat{\lambda}-\tilde{\lambda}|^{2})^{1/2} \lesssim\frac{1}{\sqrt{m}}\,,$ $\displaystyle\mathbb{E}|\Delta_{3}|$ $\displaystyle=\lambda\frac{1}{m}\mathbb{E}[\mathop{\mathrm{tr}}(T_{\setminus \ell}(\hat{\lambda})T_{\setminus\ell}(\tilde{\lambda}))|\hat{\lambda}-\tilde{ \lambda}|]\leq\mathbb{E}\left[\frac{\lambda}{\hat{\lambda}\tilde{\lambda}}| \hat{\lambda}-\tilde{\lambda}|\right]\lesssim\frac{1}{\sqrt{m}}\,,$ $\displaystyle\mathbb{E}|\Delta_{4}|$ $\displaystyle\lesssim\frac{1}{\sqrt{m}}$ $$

where we used Lemma 11 and that $\hat{\lambda}=\Omega(\lambda)$ w.p. $1$. As for $\Delta_{2}$, by Lemma 26 Item 1,

$$ $\displaystyle\mathbb{E}|\Delta_{2}|\lesssim\lambda\frac{1}{m}(\mathbb{E}\|T_{ \setminus\ell}(\tilde{\lambda})\|^{2}_{F})^{1/2}\lesssim\frac{1}{\sqrt{m}}\,.$ $$

Thus, we conclude that $\|\mathbb{E}\mathcal{E}\|=\|\mathop{\mathrm{diag}}(\mathbb{E}V^{\top}\mathcal{
E}V)\|=O(\frac{1}{\sqrt{m}})$.

∎

###### .

### Appendix J Proof of Lemma 14

###### .

Clearly, $\|\mathcal{E}\|\leq 1+\|(H+\lambda I)^{1/2}\hat{W}(H+\lambda I)^{1/2}\|$. By the inequality $\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$ (for $a,b\geq 0$), we have $\|(H+\lambda I)^{1/2}(H^{1/2}+\lambda^{1/2}I)^{-1}\|\leq 1$, so

$$ $\displaystyle\|\Phi\|\leq\|(H^{1/2}+\lambda^{1/2})\hat{W}(H^{1/2}+\lambda^{1/2 })\|\leq 2\|H^{1/2}\hat{W}H^{1/2}\|+2\lambda\|\hat{W}\|\,,$ $$

where the second inequality follows from the fact that $\hat{W}$ is PSD, and therefore for any $u,v$, $(u+v)\hat{W}(u+v)=\|u+v\|_{\hat{W}}^{2}\leq 2\|u\|_{\hat{W}}^{2}+2\|v\|_{\hat{
W}}^{2}$ where $\|u\|_{\hat{W}}^{2}:=u^{\top}\hat{W}u$ is a norm.
We have

$$ $\displaystyle H^{1/2}\hat{W}H^{1/2}=H^{1/2}S^{\top}(SHS^{\top}+\hat{\lambda}I) ^{-1}SH^{1/2}$ $\displaystyle=(H^{1/2}S^{\top}SH^{1/2}+\hat{\lambda}I)^{-1}H^{1/2}S^{\top}SH^{ 1/2}$ $\displaystyle=I-\hat{\lambda}(H^{1/2}S^{\top}SH^{1/2}+\hat{\lambda})^{-1}\,,$ $$

so that $\|H^{1/2}\hat{W}H^{1/2}\|\leq 2$.
Moreover, $\|\hat{W}\|=\|S^{\top}(SHS^{\top}+\hat{\lambda}I)^{-1}S\|\leq\frac{1}{\hat{
\lambda}}\|S^{\top}S\|$.
Recalling that $\hat{\lambda}\geq 5\lambda/12$ (by construction),
$\lambda\|\hat{W}\|\leq\frac{12}{5}\|S^{\top}S\|$. Combining all the above, $\|\mathcal{E}\|\leq 5+\frac{24}{5}\|S^{\top}S\|$.

By [vershynin2018high, Theorem 4.6.1], for every $t\geq 0$, w.p. $1-2e^{-t}$,

$$ $\displaystyle\|S^{\top}S\|\leq 1+C_{1}\max\{\mu,\mu^{2}\}\quad\textrm{where} \quad\mu=\sqrt{\frac{d}{m}}+\sqrt{\frac{t}{{m}}}\,.$ $$

for $C_{1}$ large enough.
Since $m\leq d$, note that $\max\{\mu,\mu^{2}\}=\mu^{2}\leq 2\frac{d}{m}+2\frac{t}{m}$. Item 1 follows by setting $t\sim\log(q/\delta)$.

As for Item 2,

$$ $\displaystyle\mathbb{E}[\|S^{\top}S\|\mathds{1}_{\|S^{\top}S\|>1+4C_{1}\frac{d }{m}}]$ $\displaystyle=\int_{1+4C_{1}\frac{d}{m}}^{\infty}\Pr(\|S^{\top}S\|\geq s)ds$ $\displaystyle=2C_{1}\frac{1}{m}\int_{d}^{\infty}\Pr(\|S^{\top}S\|\geq 1+2C_{1} \frac{d}{m}+2C_{1}\frac{t}{m})dt$ $\displaystyle\lesssim\frac{1}{m}\int_{d}^{\infty}e^{-t}dt=\frac{1}{m}e^{-d}\,.$ $$

Deducing Item 2 of the lemma is now easy from the above estimate.

Finally, for Item 3,

$$ $\displaystyle\mathbb{E}[\|S^{\top}S\|^{2}\mathds{1}_{\|S^{\top}S\|>1+4C_{1} \frac{d}{m}}]=\int_{1+4C_{1}\frac{d}{m}}^{\infty}2s\Pr(\|S^{\top}S\|^{2}\geq s ^{2})ds$ $\displaystyle\qquad=2C_{1}\frac{1}{m}\int_{d}^{\infty}\left(1+2C_{1}\frac{d}{m }+2C_{1}\frac{t}{m}\right)\Pr(\|S^{\top}S\|\geq 1+2C_{1}\frac{d}{m}+2C_{1} \frac{t}{m})dt$ $\displaystyle\qquad\lesssim\frac{1}{m}\int_{d}^{\infty}\left(1+2C_{1}\frac{d}{ m}+2C_{1}\frac{t}{m}\right)e^{-t}dt=O(\frac{d}{m^{2}}e^{-d})\,.$ $$

∎

###### .

### Appendix K Proof of Lemma 15

###### .

We continue from Section I, expressing $\mathcal{E}$ in the eigenbasis of $H$.

Recall that $\mathcal{E}=(H+\lambda)^{1/2}S^{\top}T(\hat{\lambda})S(H+\lambda)^{1/2}-I$.
Then starting from (71),
we have for any $\ell,j=1,\ldots,d$,

$$ $\displaystyle v_{\ell}^{\top}\mathcal{E}v_{j}=(\tau_{\ell}+\lambda)^{1/2}(\tau _{j}+\lambda)^{1/2}\frac{s_{j}T_{\setminus\ell}(\hat{\lambda})s_{\ell}}{1+\tau _{\ell}s_{\ell}T_{\setminus\ell}(\hat{\lambda})s_{\ell}}-\mathds{1}_{\ell=j}\,.$ (75) $$

Consequently, if $j\neq k$ then note that $\mathbb{E}[v_{j}^{\top}\mathcal{E}v_{\ell}v_{k}^{\top}\mathcal{E}v_{\ell}]=0$, since $s_{j},s_{k}$ have a symmetric distribution.

Recall, we are interested in bounding $\|\mathbb{E}[\mathcal{E}^{2}]\|=\|\mathbb{E}[(V^{\top}\mathcal{E}V)^{2}]\|$. We have

$$ $\displaystyle\mathbb{E}[(V^{\top}\mathcal{E}V)^{2}_{j,k}]=\sum_{\ell=1}^{d} \mathbb{E}[(V^{\top}\mathcal{E}V)_{j,\ell}(V^{\top}\mathcal{E}V)_{k,\ell}]\,,$ $$

which is zero when $j\neq k$; that is, $\mathbb{E}[(V^{\top}\mathcal{E}V)^{2}]$ is diagonal. Let us compute the diagonal elements, $j=k$. We have

$$ $\displaystyle\mathbb{E}[(V^{\top}\mathcal{E}V)^{2}_{k,k}]$ $\displaystyle=\mathbb{E}[(V^{\top}\mathcal{E}V)_{k,k}(V^{\top}\mathcal{E}V)_{k ,k}]+\sum_{\ell\neq k}\mathbb{E}[(V^{\top}\mathcal{E}V)_{k,\ell}(V^{\top} \mathcal{E}V)_{k,\ell}]\,.$ (76) $$

The first term is

$$ $\displaystyle\mathbb{E}[(V^{\top}\Phi V)_{k,k}^{2}]$ $\displaystyle=\mathbb{E}\left[\left(\frac{\lambda-(s_{\ell}^{\top}T_{\setminus \ell}(\hat{\lambda})s_{\ell})^{-1}}{(s_{\ell}^{\top}T_{\setminus\ell}(\hat{ \lambda})s_{\ell})^{-1}+\tau_{\ell}}\right)^{2}\right]$ $\displaystyle=O(|\lambda s_{\ell}^{\top}T_{\setminus\ell}(\hat{\lambda})s_{ \ell}-1|^{2})\,,$ (77) $$

where we used (74). Similar to Section I, one can show this is $O(1/m)$. For the second term, using (75),

$$ $\displaystyle\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{ \top}\mathcal{E}V)_{k,\ell}=\sum_{\ell=1,\ell\neq k}^{d}(\tau_{k}+\lambda)( \tau_{\ell}+\lambda)\frac{(s_{k}T_{\setminus k}(\hat{\lambda})s_{\ell})^{2}}{( 1+\tau_{k}s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k})^{2}}$ $\displaystyle\qquad=(\tau_{k}+\lambda)\frac{1}{(1+\tau_{k}s_{k}T_{\setminus k} (\hat{\lambda})s_{k})^{2}}s_{k}^{\top}T_{\setminus k}(\hat{\lambda})\left[\sum _{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)s_{\ell}s_{\ell}^{\top}\right]T_{ \setminus k}(\hat{\lambda})s_{k}\,.$ (78) $$

For small $c^{\prime}>0$, let $\Omega$ be the event that $s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k}\geq c^{\prime}/\lambda$; by similar calculations as in Section I, for some $c,C>0$, $\Pr(\Omega^{c})\leq Ce^{-cm}$. Then

$$ $\displaystyle\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{ \top}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega}$ $\displaystyle\lesssim\frac{\lambda^{2}}{\tau_{k}+\lambda}s_{k}^{\top}T_{ \setminus k}(\hat{\lambda})\left[\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+ \lambda)s_{\ell}s_{\ell}^{\top}\right]T_{\setminus k}(\hat{\lambda})s_{k} \mathds{1}_{\Omega}$ $\displaystyle\leq\frac{\lambda^{2}}{\tau_{k}+\lambda}s_{k}^{\top}T_{\setminus k }(\hat{\lambda})\left[\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)s_{\ell }s_{\ell}^{\top}\right]T_{\setminus k}(\hat{\lambda})s_{k}\,,$ $$

so taking the expectation,

$$ $\displaystyle\mathbb{E}\left[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V )_{k,\ell}(V^{\top}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega}\right]$ $\displaystyle\lesssim\frac{\lambda^{2}}{\tau_{k}+\lambda}\frac{1}{m}\mathbb{E} \mathop{\mathrm{tr}}\left(T_{\setminus k}(\hat{\lambda})\left[\sum_{\ell=1, \ell\neq k}^{d}(\tau_{\ell}+\lambda)s_{\ell}s_{\ell}^{\top}\right]T_{\setminus k }(\hat{\lambda})\right)$ $\displaystyle\lesssim\frac{\lambda}{\tau_{k}+\lambda}\frac{1}{m}\mathbb{E} \mathop{\mathrm{tr}}\left(\left[\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+ \lambda)s_{\ell}s_{\ell}^{\top}\right]T_{\setminus k}(\hat{\lambda})\right)$ $\displaystyle\leq\frac{1}{m}\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+\lambda) \mathbb{E}s_{\ell}^{\top}T_{\setminus k}(\hat{\lambda})s_{\ell}\,.$ $$

Using (71),

$$ $\displaystyle\mathbb{E}s_{\ell}^{\top}T_{\setminus k}(\hat{\lambda})s_{\ell}= \mathbb{E}\left[\frac{s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda})s_{ \ell}}{1+\tau_{\ell}s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda})s_{\ell} }\right]\leq\frac{\mathbb{E}[s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda} )s_{\ell}]}{1+\tau_{\ell}\mathbb{E}[s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{ \lambda})s_{\ell}]}\,,$ (79) $$

where we used Jensen’s inequality with the concave function $x\mapsto\frac{x}{1+\tau_{\ell}x}$, $x>0$. By a calculation similar to Section I, $\mathbb{E}[s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda})s_{\ell}]=\frac{1
}{\lambda}+O(\frac{1}{\lambda\sqrt{m}})$, hence

$$ $\displaystyle\mathbb{E}s_{\ell}^{\top}T_{\setminus k}(\hat{\lambda})s_{\ell} \leq\frac{\mathbb{E}[s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda})s_{\ell }]}{1+\tau_{\ell}\mathbb{E}[s_{\ell}^{\top}T_{\setminus k,\ell}(\hat{\lambda}) s_{\ell}]}\lesssim\frac{1}{\tau_{\ell}+\lambda}\,.$ $$

Consequently,

$$ $\displaystyle\mathbb{E}\left[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V )_{k,\ell}(V^{\top}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega}\right]$ $\displaystyle\lesssim\frac{1}{m}\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+ \lambda)\mathbb{E}s_{\ell}^{\top}T_{\setminus k}(\hat{\lambda})s_{\ell} \lesssim\frac{d}{m}\,.$ (80) $$

Finally, towards bounding the expectation under the complement event ${\Omega^{c}}$ note we can write

$$ $\displaystyle\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{ \top}\mathcal{E}V)_{k,\ell}$ $\displaystyle=(\tau_{k}+\lambda)\frac{1}{(1+\tau_{k}s_{k}T_{\setminus k}(\hat{ \lambda})s_{k})^{2}}\sum_{\ell=1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)(s_{\ell} ^{\top}T_{\setminus k}(\hat{\lambda})s_{k})^{2}$ $\displaystyle\leq(\tau_{k}+\lambda)\frac{s_{k}^{\top}T_{\setminus k}(\hat{ \lambda})s_{k}}{(1+\tau_{k}s_{k}T_{\setminus k}(\hat{\lambda})s_{k})^{2}}\sum_ {\ell=1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)s_{\ell}^{\top}T_{\setminus k}( \hat{\lambda})s_{\ell}$ $\displaystyle\leq(\tau_{k}+\lambda)\frac{s_{k}^{\top}T_{\setminus k}(\hat{ \lambda})s_{k}}{1+\tau_{k}s_{k}T_{\setminus k}(\hat{\lambda})s_{k}}\sum_{\ell= 1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)s_{\ell}^{\top}T_{\setminus k}(\hat{ \lambda})s_{\ell}\,,$ (81) $$

where the first inequality, $(s_{\ell}^{\top}T_{\setminus k}(\hat{\lambda})s_{k})^{2}\leq(s_{\ell}^{\top}T_
{\setminus k}(\hat{\lambda})s_{\ell})(s_{k}^{\top}T_{\setminus k}(\hat{\lambda
})s_{k})$, follows by Cauchy-Schwartz. We have

$$ $\displaystyle(\tau_{k}+\lambda)\frac{s_{k}^{\top}T_{\setminus k}(\hat{\lambda} )s_{k}}{1+\tau_{k}s_{k}T_{\setminus k}(\hat{\lambda})s_{k}}$ $\displaystyle=\frac{\tau_{k}+\lambda}{\tau_{k}+(s_{k}^{\top}T_{\setminus k}( \hat{\lambda})s_{k})^{-1}}$ $\displaystyle=1+\frac{\lambda-(s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k} )^{-1}}{\tau_{k}+(s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k})^{-1}}$ $\displaystyle\leq 1+|\lambda(s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k})- 1|\,.$ $$

Recall $\Omega$, the event that $\lambda s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k}\geq c^{\prime}$; we now need to bound $\mathbb{E}[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{\top
}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega^{c}}]$. Under $\Omega^{c}$, $0\leq\lambda s_{k}^{\top}T_{\setminus k}(\hat{\lambda})<c^{\prime}$, and so $(\tau_{k}+\lambda)\frac{s_{k}^{\top}T_{\setminus k}(\hat{\lambda})s_{k}}{1+
\tau_{k}s_{k}T_{\setminus k}(\hat{\lambda})s_{k}}\mathds{1}_{\Omega^{c}}=O(1)$. Plugging this into (81),

$$ $\displaystyle\mathbb{E}[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k, \ell}(V^{\top}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega^{c}}]\lesssim\sum_{\ell =1,\ell\neq k}^{d}(\tau_{\ell}+\lambda)\mathbb{E}[s_{\ell}^{\top}T_{\setminus k }(\hat{\lambda})s_{\ell}\mathds{1}_{\Omega^{c}}]\,.$ (82) $$

Using (79),

$$ $\displaystyle(\tau_{\ell}+\lambda)\mathbb{E}[s_{\ell}^{\top}T_{\setminus k}( \hat{\lambda})s_{\ell}\mathds{1}_{\Omega^{c}}]$ $\displaystyle=\mathbb{E}\left[\frac{\tau_{\ell}+\lambda}{\tau_{\ell}+(s_{\ell} ^{\top}T_{\setminus k,\ell}(\hat{\lambda})s_{\ell})^{-1}}\mathds{1}_{\Omega^{c }}\right]$ $\displaystyle\leq\mathbb{E}\left[\left(1+|\lambda(s_{\ell}^{\top}T_{\setminus k ,\ell}(\hat{\lambda})s_{\ell})-1|\right)\mathds{1}_{\Omega^{c}}\right]$ $\displaystyle\lesssim\mathbb{E}[(1+\|s_{\ell}\|^{2})\mathds{1}_{\Omega^{c}}] \lesssim\Pr(\Omega^{c})+\sqrt{\mathbb{E}\|s_{\ell}\|^{4}\Pr(\Omega^{c})} \lesssim e^{-cm}\,,$ $$

as $\Pr(\mathcal{E}^{c})\lesssim e^{-cm}$. Thus, $\mathbb{E}[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{\top
}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega^{c}}]\leq de^{-cm}$. And so, we finally conclude

$$ $\displaystyle\mathbb{E}[\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\Phi V)_{k,\ell}( V^{\top}\Phi V)_{k,\ell}]$ $\displaystyle=\mathbb{E}\sum_{\ell=1,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k, \ell}(V^{\top}\mathcal{E}V)_{k,\ell}\mathds{1}_{\Omega}+\mathbb{E}\sum_{\ell=1 ,\ell\neq k}^{d}(V^{\top}\mathcal{E}V)_{k,\ell}(V^{\top}\mathcal{E}V)_{k,\ell} \mathds{1}_{\Omega^{c}}$ $\displaystyle\lesssim\frac{d}{m}+de^{-cm}=O(d/m)\,.$ (83) $$

∎

###### .

### Appendix L Proof of Theorem 25

Define the Newton decrement at a point $\theta\in\mathbb{R}^{d}$:

$$ $\mathsf{N}(\theta)=\left((\nabla G(\theta))^{\top}(\nabla^{2}G(\theta))^{-1}( \nabla G(\theta))\right)^{1/2}.$ (84) $$

Denote the approximate Newton decrement, where $\bar{W}(\theta)$ is an approximation of the true inverse Hessian:

$$ $\tilde{\mathsf{N}}(\theta)=\left((\nabla G(\theta))^{\top}\bar{W}(\theta)( \nabla G(\theta))\right)^{1/2}.$ (85) $$

Note that if $\bar{W}$ is $\eta$-accurate in the sense of Definition 2, then

$$ $\sqrt{1-\eta}\mathsf{N}(\theta)\leq\tilde{\mathsf{N}}(\theta)\leq\sqrt{1+\eta} \mathsf{N}(\theta)\,.$ (86) $$

For self-concordant functions, it is known that the Newton decrement yields an upper bound on the suboptimality gap $G(\theta)-G(\theta^{\star})$:

###### Lemma 16 .

If $\mathsf{N}(\theta)<\mathsf{N}_{0}$ for some numerical constant $\mathsf{N}_{0}\leq 0.68$, then

$$ $G(\theta)-G(\theta^{\star})\leq\mathsf{N}^{2}(\theta)\,.$ $$

###### Proof.

See \citep[Section 9.6.3]boyd2004convex, specifically Eq. (95).
∎

Clearly, by (86), the approximate Newton decrement $\tilde{\mathsf{N}}(\theta)$ provides an upper bound on the optimality gap provided that the inverse Hessian is $\eta$-accurate, albeit a looser one.

The key tool in our analysis are the following estimates, due to \citeppilanci2015newton.

###### Lemma 17 .

Operate under the conditions of Theorem 25. There exist numerical constants $\Lambda,\nu>0$ with $\Lambda<1/16$ such that the following holds. Assume that $0\leq\eta<1/2$.

- •
If at an iteration $t$, one has $\mathsf{N}(\theta_{t-1})\geq\Lambda$, then at the next (approximate) Newton step: $G(\theta_{t})-G(\theta^{\star})\leq-ab\nu$, where $a,b>0$ are the parameters for backtracking line search (see Algorithm 3).
- •
If $\mathsf{N}(\theta_{t-1})\leq\Lambda$ then at the next iteration $\mathsf{N}(\theta_{t})\leq\mathsf{N}(\theta_{t-1})$.

###### Proof.

See \citep[Lemma 6.4]pilanci2015newton.
∎

Note that Lemma 17 implies that within at most

$$ $T_{0}=(G(\theta_{0})-G(\theta^{\star}))/(ab\nu)$ (87) $$

iterations, we are guaranteed to reach $\theta_{t}$ such that $\mathsf{N}(\theta_{t})\leq\Lambda$; furthermore, once we have achieved that, $\mathsf{N}(\theta_{t^{\prime}})\leq\Lambda$ holds for every subsequent iteration $t^{\prime}\geq t$.

The following is the main estimate for the remainder of the analysis:

###### Lemma 18 .

Operate under the conditions of Theorem 25, and assume that $\eta<1/2$. Then

$$ $\mathsf{N}(\theta_{t})\leq\frac{(1+\eta)\mathsf{N}^{2}(\theta_{t-1})+\eta \mathsf{N}(\theta_{t-1})}{\left(1-(1+\eta)\mathsf{N}(\theta_{t-1})\right)^{2}}\,.$ (88) $$

###### Proof.

This result is \citep[Lemma 6.6]pilanci2015newton.
∎

One can further simplify (88) assuming that $\mathsf{N}(\theta_{t-1})<1/16$. Coarsely lower bounding the denominator and upper bounding the numerator:

$$ $\mathsf{N}(\theta_{t})\leq 2\mathsf{N}^{2}(\theta_{t-1})+2\eta\mathsf{N}( \theta_{t-1})\leq\begin{cases}4\mathsf{N}^{2}(\theta_{t-1})\quad&\textrm{if} \quad\mathsf{N}(\theta_{t-1})\geq\eta,\\ 4\eta\mathsf{N}(\theta_{t-1})\quad&\textrm{if}\quad\mathsf{N}(\theta_{t-1}) \leq\eta\end{cases}\,,$ (89) $$

provided that $\mathsf{N}(\theta_{t-1})<1/16$.

With the above estimates in hand, we are ready to prove Theorem 25.

###### Proof.

(Of Theorem 25.)
Recall that we assume that $\eta<1/5$, and in particular $\eta<1/2$. By Lemma 17,
within $T_{0}$ iterations, the approximate Newton method reaches $\theta_{t}$ such that $\mathsf{N}(\theta_{t^{\prime}})<1/16$ for all $t^{\prime}\geq t$. In particular, also $G(\theta)-G(\theta^{\star})\leq(1+\eta)1/16^{2}\leq 3/500$.
Let $\varepsilon>0$ be the desired precision; suppose that $\varepsilon<3/500$.
We consider two cases.

First, assume that $\eta\leq\varepsilon<3/500$. Let $T_{1}(\varepsilon)$ be the smallest integer $t$ such that the dynamic

$$ $A_{t}=4A_{t-1}^{2},\quad A_{0}=1/16$ (90) $$

satisfies $A_{t}\leq\varepsilon$. It is easy to verify that $A_{t}=1/4(1/4)^{2^{t}}$ (e.g. by induction), hence $T_{1}(\varepsilon)=O(\log\log(1/\varepsilon))$. We conclude that when $\varepsilon\geq\eta$, $\mathsf{T}(\varepsilon)\leq T_{0}+T_{1}(\varepsilon)$.

The second case is when $\varepsilon<\eta$, therefore we are in the regime of (89) where decay is linear rather than quadratic.
Let $T_{2}(\varepsilon)$ be the smallest integer $t$ such that the dynamic

$$ $B_{t}=4\eta B_{t-1},\quad B_{0}=\eta$ (91) $$

satisfies $B_{t}\leq\varepsilon$. Clearly $B_{t}=(4\eta)^{t}B_{0}$, hence $T_{2}(\varepsilon)=O(\frac{\log(\eta/\varepsilon)}{\log(1/\eta)})$ (here we used that $4\eta<1$, by assumption). We conclude that when $\varepsilon<\eta$, $\mathsf{T}(\varepsilon)\leq T_{0}+T_{1}(\eta)+T_{2}(\varepsilon)$.

∎

###### Lemma 16 .

###### Proof.

###### Lemma 17 .

###### Proof.

###### Lemma 18 .

###### Proof.

###### Proof.

### Appendix M Auxilliary Technical Results

###### Lemma 19 (Sherman-Morrison) .

Let $A\in\mathbb{R}^{n\times n}$ be invertible, and $u,v\in\mathbb{R}^{n}$.

The matrix $A+uv^{\top}$ is invertible if and only if $1+v^{\top}A^{-1}u\neq 0$. In that case,

$$ $(A+uv^{\top})^{-1}=A^{-1}-\frac{A^{-1}uv^{\top}A^{-1}}{1+v^{\top}A^{-1}u}\,.$ (92) $$

###### Lemma 20 (Block matrix inversion) .

Let $A,B,C,D$ be matrices of conforming dimensions. Provided that $A,D$ are both invertible,

$$ $\begin{bmatrix}A&B\\ C&D\end{bmatrix}^{-1}=\begin{bmatrix}(A-BD^{-1}C)^{-1}&0\\ 0&(D-CA^{-1}B)^{-1}\end{bmatrix}\begin{bmatrix}I&-BD^{-1}\\ -CA^{-1}&I\end{bmatrix}\,.$ (93) $$

###### Lemma 21 .

(Resolvent low-rank perturbation.)
Let $A,B,C\succeq 0$ be PSD matrices, and $z>0$. Then

$$ $\left|\mathop{\mathrm{tr}}C(zI+A)^{-1}-\mathop{\mathrm{tr}}C(zI+B)^{-1}\right| \leq\mathop{\mathrm{rank}}(A-B)\frac{\|C\|}{z}\,.$ $$

For reference, see for example \citepbai2010spectral.

###### Lemma 19 (Sherman-Morrison) .

###### Lemma 20 (Block matrix inversion) .

###### Lemma 21 .

##### 

Let $\{\emptyset,\Omega\}=:\mathcal{F}_{0}\leq\mathcal{F}_{1}\leq\mathcal{F}_{2}\leq\ldots$ be a filtration of $\sigma$-algebras over some given probability space $\Omega$. Recall that a random process $\{X_{n}\}_{n=1}^{\infty}$ is a *martingale difference* (adapted to this filtration) if 1) for every $n$, $X_{n}$ is $\mathcal{F}_{n}$-measurable; and 2) $\mathcal{F}_{n-1}$-almost surely, $\mathbb{E}[X_{n}|\mathcal{F}_{n-1}]=0$.

###### Lemma 22 (Azuma-Hoeffding) .

Suppose that $(X_{n},\mathcal{F}_{n})_{n=1}^{N}$ is a martingale difference, and $\{B_{n}\}_{n=1}^{N}$ constants. Suppose that for all $1\leq n\leq N$, almost surely, $|X_{n}|\leq B_{n}$. Then for all $t\geq 0$,

$$ $\Pr(|\sum_{n=1}^{N}X_{n}|\geq t)\leq 2\exp\left(\frac{-\frac{1}{2}t^{2}}{\sum_ {n=1}^{N}B_{n}^{2}}\right)\,.$ $$

The sub-Gaussian and sub-Exponential norms of a (scalar) random variable are, respectively

$$ $\displaystyle\|X\|_{\psi_{2}}=\inf\left\{\sigma>0\;:\;\mathbb{E}e^{X^{2}/ \sigma^{2}}\leq 2\right\}\,,\qquad\|X\|_{\psi_{1}}=\inf\left\{\sigma>0\;:\; \mathbb{E}e^{|X|/\sigma}\leq 2\right\}\,.$ (94) $$

We call $X$ sub-Gaussian, resp. sub-Exponential, if $\|X\|_{\psi_{2}}<\infty$, resp. $\|X\|_{\psi_{1}}<\infty$.

For a random vector $r$, we define $\|r\|_{\psi_{i}}=\sup_{\|v\|=1}\|v^{\top}r\|_{\psi_{i}}$. That is, it is the largest $\psi_{i}$-norm of a one-dimensional projection of $r$.

###### Lemma 23 .

The following holds.

- 1.
$\|\cdot\|_{\psi_{i}}$ is a norm on the subspace of random variables $X$ such that $\|X\|_{\psi_{i}}<\infty$.
- 2.
If $X,Y$ are sub-Gaussian (possible statistically dependent), then $XY$ is sub-Exponential, with $\|XY\|_{\psi_{1}}\leq\|X\|_{\psi_{2}}\|Y\|_{\psi_{2}}$.
- 3.
If $|Y|\leq B$ almost surely for $B\geq 0$ constant, $\|XY\|_{\psi_{i}}\leq B\|X\|_{\psi_{i}}$.
- 4.
(Centralization.) For any $X$, $\|X-\mathbb{E}[X]\|_{\psi_{i}}\leq\|X\|_{\psi_{i}}$.

We refer to the book \citepvershynin2018high for reference.

###### Lemma 24 (Bernstein) .

Suppose that $(X_{n},\mathcal{F}_{n})_{n=1}^{N}$ is a martingale difference. Then for all $t\geq 0$,

$$ $\Pr(|\sum_{n=1}^{N}X_{n}|\geq t)\leq 2\exp\left(-c\min\{\frac{t}{\max_{n=1, \ldots,N}\|X_{n}\|_{\psi_{1}}},\frac{t^{2}}{\sum_{n=1}^{N}\|X_{n}\|_{\psi_{1}} ^{2}}\}\right)\,.$ $$

###### Lemma 22 (Azuma-Hoeffding) .

###### Lemma 23 .

###### Lemma 24 (Bernstein) .

##### 

Next, we cite the following version of the Hanson-Wright inequality, see \citepvershynin2018high.

###### Lemma 25 (Hanson-Wright) .

Suppose that $r$ has independent sub-Gaussian entries, with $\mathbb{E}[r_{i}]=0,\mathbb{E}[r_{i}^{2}]=1$, $\max_{i}\|r_{i}\|_{\psi_{2}}\leq K$.
For any matrix $A$ and $t\geq 0$,

$$ $\displaystyle\Pr(|r^{\top}Ar-\mathop{\mathrm{tr}}(A)|\geq t)\leq\exp\left(-c \min\{\frac{t}{K^{2}\|A\|},\frac{t^{2}}{K^{4}\|A\|_{F}^{2}}\}\right)\,,$ (95) $$

where $c>0$ is a universal constant.

Lastly,
the following lemma collects, for convenience, some properties satisfied by any isotropic random vector $r$ that satisfies a concentration inequality of the form (95).

###### Lemma 26 .

Suppose $r$ is istropic $\mathbb{E}[r]=0,\mathbb{E}[rr^{\top}]=I$, and satisfies (95) for some $K>0$.

- 1.
For any matrix $A$ (independent of $r$), $\mathbb{E}[|r^{\top}Ar-\mathop{\mathrm{tr}}(A)|^{2}]\lesssim K^{4}\|A\|_{F}^{2}$.
- 2.
For any vector $v$ (independent of $r$), $\|v^{\top}r\|_{\psi_{2}}\lesssim K\|v\|$. That is, $r$ is sub-Gaussian, with $\|r\|_{\psi_{2}}\lesssim K$.

###### Proof.

Start with the first item.
Write $\mathbb{E}[|r^{\top}Ar-\mathop{\mathrm{tr}}(A)|^{2}]=\int_{0}^{\infty}\Pr(|r^{
\top}Ar-\mathop{\mathrm{tr}}(A)|^{2}>s)ds=2\int_{0}^{\infty}t\Pr(|r^{\top}Ar-
\mathop{\mathrm{tr}}(A)|^{2}>t)dt$. Using (95), we crudely bound

$$ $\displaystyle\int_{0}^{\infty}t\Pr(|r^{\top}Ar-\mathop{\mathrm{tr}}(A)|^{2}>t)dt$ $\displaystyle\leq\int_{0}^{\infty}t\exp(-c\frac{t^{2}}{K^{4}\|A\|_{F}^{2}})dt+ \int_{0}^{\infty}t\exp(-c\frac{t}{K^{2}\|A\|})dt\lesssim K^{4}\|A\|_{F}^{2}\,.$ $$

As for the second item, note that $\|vv^{T}\|=\|vv^{T}\|_{F}=\|v\|^{2}$. Hence,

$$ $\Pr(|v^{\top}r|\geq t)=\Pr(r^{\top}vv^{\top}r\geq t^{2})\lesssim\exp(-\min\{ \frac{t^{2}}{K^{2}\|v\|^{2}},\frac{t^{4}}{K^{4}\|v\|^{4}}\})\,.$ $$

That is, $\|v^{\top}r\|$ has a Gaussian tail with variance proxy $K^{2}\|v\|^{2}$.
∎

###### Lemma 25 (Hanson-Wright) .

###### Lemma 26 .

###### Proof.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
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.

2) Does this method and approximation work for pure convex functions without regularisation? Also, may they work for non-convex functions in practice? 

3) How good is the proposed approximation compared to Krylov and Nystorm-based approximations?

4) It might be interesting to apply this approximation to the new branch of methods based on cubic regularisation:

 [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 

[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)

[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.
soundness: 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.
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.
weaknesses: 1. One of the main weaknesses from my point of view is the experiments section. 
a) 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. 
b) 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. 

2. 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}$.

3. Also, the authors may improve the related literature section by some distributed Newton and Quasi-Newton methods. Some possible references:

 [1] Safaryan, M. et.al. "FedNL: Making Newton-type methods applicable to federated learning." (2021) arXiv preprint arXiv:2106.02969. 

[2] Agafonov, A. et.al. "Flecs: A federated learning second-order framework via compression and sketching". (2022) arXiv preprint arXiv:2206.02009. 

[3] Qian, X. et.al  "Basis matters: better communication-efficient second order methods for federated learning." (2021) arXiv preprint arXiv:2111.01847. 

[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. 

[5] Bischoff, S. et.al. "On second-order optimization methods for federated learning." (2021) arXiv preprint arXiv:2109.02388.

---

confidence: 4
contribution: 2
initial_rating: 6
presentation: 2
questions: See 3 in Weaknesses
soundness: 3
strengths: 1) A new (best of my knowledge) matrix sketching scheme is presented.
2) The proofs look correct.
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.
weaknesses: 1) Convergence only for self condordant functions. The result would be more solid if the authors consider general strongly convex functions.

2) 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]

[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.

[2] Amari, S. (2021). Information geometry.

[3] Gupta, V., Koren, T., and Singer, Y. (2018). Shampoo: Preconditioned stochastic tensor optimization.

[4] Vyas, N., Morwani, D., Zhao, R., Shapira, I., Brandfonbrener, D., Janson, L., and Kakade, S. (2024). Soap: Improving and stabilizing shampoo using adam

[5] Doikov, N., Jaggi, M., et al. (2023). Second-order optimization with lazy hessians. 

3) 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.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: No
soundness: 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.
The theory in this paper is good.
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.
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.
The convergence rate of the proposed algorithm with $m$ machines running in parallel is the same to one only runing in  a single machine.

The experiments prove the biases of the proposed algorithm is lower than other agorithm and the algorithm can achieve lower iteration complexity.
However, the coresponding computation cost is not mentioned. Thus, the comparisons in this paper seem not enough to prove the advantages of the proposed algorithm.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 2
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. 

2. 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$.

3. (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?
soundness: 3
strengths: The solution in this paper is novel to me. I think the idea is very interesting.
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.
weaknesses: See the questions part.

Meta-review: 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.

All 4 reviewers (scores 8, 8, 6, 6) agreed that the paper should be accepted. I have read the reviews, the rebuttals, and the discussion. Moreover, I have had a quick read of the main body of the paper. I agree that the paper develops some very interesting new results, with potential for further extensions. 

- The method and the theoretical results are new and interesting. 
- The paper is well written.
- Related work is (after rebuttal) in a good shape.

I recommend the paper for acceptance.

Decision: Accept (Poster)