Distributed Least Squares in Small Space via Sketching and Bias Reduction

Authors: sachin  garg
kevin  tan
michał  derezi
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 05:37:42

Abstract: Matrix sketching is a powerful tool for reducing the size of large data matrices. Yet there are fundamental limitations to this size reduction when we want to recover an accurate estimator for a task such as least square regression. We show that these limitations can be circumvented in the distributed setting by designing sketching methods that minimize the bias of the estimator, rather than its error. In particular, we give a sparse sketching method running in optimal space and current matrix multiplication time, which recovers a nearly-unbiased least squares estimator using two passes over the data. This leads to new communication-efficient distributed averaging algorithms for least squares and related tasks, which directly improve on several prior approaches. Our key novelty is a new bias analysis for sketched least squares, giving a sharp characterization of its dependence on the sketch sparsity. The techniques include new higher moment restricted Bai-Silverstein inequalities, which are of independent interest to the non-asymptotic analysis of deterministic equivalents for random matrices that arise from sketching.

Main paper: 
 introduction: 
 Matrix sketching is a powerful collection of randomized techniques for compressing large data matrices, developed over a long line of works as part of Randomized Numerical Linear Algebra [RandNLA, e.g., 45, 26, 37, 39, 21] . Sketching can be used to reduce the large dimension n of a data matrix A ∈ R n×d by applying a random sketching matrix (operator) S ∈ R m×n to obtain the sketch Ã = SA ∈ R m×d where m ≪ n. For example, sketching can be used to approximate the solution to the least squares problem, x * = argmin x L(x) where L(x) = ∥Ax -b∥ 2 , by using a sketched estimator x = argmin x ∥ Ãx -b∥ 2 , where Ã = SA and b = Sb.Perhaps the simplest form of sketching is subsampling, where the sketching operator S selects a random sample of the rows of matrix A. However, the real advantage of sketching as a framework emerges as we consider more complex operators S, such as sub-Gaussian matrices [1] , randomized Hadamard transforms [3] , and sparse random matrices [12] . These approaches ensure higher quality and more robust compression of the data matrix, e.g., leading to provable ϵ-approximation guarantees for the estimate x in the least squares task, i.e., L(x) ≤ (1 + ϵ)L(x * ). Nevertheless, there are fundamental limitations to how far we can compress a data matrix using sketching while ensuring an ϵ-approximation. These limitations pose a challenge particularly in space-limited computing 38th Conference on Neural Information Processing Systems (NeurIPS 2024). environments, such as for streaming algorithms where we observe the matrix A, say, one row at a time, and we have limited space for storing the sketch [11] .One strategy for overcoming the fundamental limitations of sketching as a compression tool is to look beyond the single approximation guarantee provided by a sketching-based estimator x, and consider how its broader statistical properties can be leveraged in a given computing environment. To that end, many recent works have demonstrated both theoretically and empirically that sketching-based estimators often exhibit not only approximation robustness but also statistical robustness, for instance enjoying sharp confidence intervals, effectiveness of statistical inference tools such as bootstrap and cross-validation, as well as accuracy boosting techniques such as distributed averaging [e.g., 35, 25, 32, 33 ]. Yet, these results have had limited impact on the traditional computational complexity analysis in RandNLA and sketching literature, as many of them either impose additional assumptions, focus on sharpening the constant factors, or require using more expensive sketching techniques. In this work, we demonstrate that the statistical properties of sketching-based estimators can in fact have a substantial impact on the computational trade-offs that arise in RandNLA.Our key motivating example is the above mentioned least squares regression task. It is well understood that for an n × d least squares task, to recover an ϵ-approximate solution out of an m × d sketch, we need sketch size at least m = Ω(d/ϵ). This has been formalized in the streaming setting with a lower bound of Ω(ϵ -1 d 2 log(nd)) bits of space required, when all of the input numbers use O(log(nd)) bits of precision [11] . One setting where this can be circumvented is in the distributed computing model where the bits can be spread out across many machines, so that the per-machine space can be smaller.Here, one could for instance hope that we can maintain small O(d) × d sketches in q = O(1/ϵ) machines and then combine their estimates to recover an ϵ-approximate solution. A simple and attractive approach is to average the estimates xi produced by the individual machines, returning x = 1 q q i=1 xi , as this only requires each machine to communicate O(d log(nd)) bits of information about its sketch. This approach requires the sketching-based estimates xi to have sufficiently small bias for the averaging scheme to be effective. While this has been demonstrated empirically in many cases, existing theoretical results still require relatively expensive sketching methods to recover low-bias estimators, leading to an unfortunate trade-off in the distributed averaging scheme between the time and space complexity required.In this work, we address the time-space trade-off in distributed averaging of sketching-based estimators, by giving a sharp characterization of how their bias depends on the sparsity of the sketching matrix. Remarkably, we show that in the distributed streaming environment one can compress the data down to the minimum size of O(d 2 log(nd)) bits at no extra computational cost, while still being able to recover an ϵ-approximate solution for least squares and related problems. Importantly, our results require the sketching matrix to be slightly denser than is necessary for obtaining approximation guarantees on a single estimate, and thus, cannot be recovered by standard RandNLA sampling methods such as approximate leverage score sampling [27] .Before we state the full result in the distributed setting, we give our main technical contribution. which is the following efficient construction of a low-bias least squares estimator in a single pass using only O(d 2 log(nd)) bits of space, assuming all numbers use O(log(nd)) bits of precision. Below, γ > 0 denotes an arbitrarily small constant.Total runtime Parallel passes Folklore Gaussian sketch O(nd ω-1 ) 1 [43] Leverage Score Sampling nnz(A) + Õ(d ω / √ ϵ) 1 [5] Determinantal Point Process nnz(A) + Õ(d ω ) log 3 (n/ϵ) [9] Weighted Mb-SGD (sequential) nnz(A) + Õ(d 2 /ϵ) 1/ϵ This work (Thm. 1) Leverage Score Sparsification nnz(A) + Õ(d 2 /ϵ) 1Table 1 : Comparison of time complexities and parallel passes over the data required for different methods to obtain a (1 + ϵ)-approximation in O(d 2 log(nd)) bits of space for an n × d least squares problem (A, b), given a preconditioner P such that κ(AP) = O(1) (see Section 3 for our computational model). We include the fully sequential Weighted Mb-SGD as a reference.Theorem 1. Given streaming access to A ∈ R n×d and b ∈ R n , and direct access to a preconditioner matrix P ∈ R d×d such that κ(AP) ≤ α, within a single pass over (A, b), in O(γ -1 nnz(A) + ϵ -1 αd 2+γ polylog(d)) time and O(d 2 log(nd)) bits of space, we can construct a randomized estimator x for the least squares solution x * = argmin x ∥Ax -b∥ 2 such that:(Bias) AE[x] -b 2 ≤ (1 + ϵ)∥Ax * -b∥ 2 , (Variance) E ∥Ax -b∥ 2 ≤ 2 ∥Ax * -b∥ 2 .Remark 1. The above construction assumes access to a preconditioner matrix P with κ(AP) ≤ α (where κ denotes the condition number). Such matrix can be obtained efficiently with α = O(1) in a separate single pass, leading to a two-pass algorithm described later in Theorem 2.Our estimator x is constructed at the end of the data pass from a sketch ( Ã, b) where Ã = SA and b = Sb, by minimizing ∥ Ãx -b∥ 2 using preconditioned conjugate gradient. Here, S is a carefully constructed sparse sketching matrix which is inspired by the so-called leverage score sparsified (LESS) embeddings [18] . Leverage scores represent the relative importances of the rows of A which are commonly used for subsampling in least squares (see Definition 1), and their estimates can be easily obtained in a single pass by using the preconditioner matrix P.Our time complexity bound of Õ(nnz(A) + d 2 /ϵ) matches the time it would take (for a single machine) to subsample Õ(d/ϵ) rows of A according to the approximate leverage scores and produce an estimator x that achieves the ϵ-error bound ∥Ax -b∥ 2 ≤ (1 + ϵ)∥Ax * -b∥ 2 . However, this strategy requires either maintaining Õ(d 2 /ϵ) bits of space for the sketch, or computing x directly along the way, blowing up the runtime to Õ(d ω /ϵ). Since approximate leverage score sampling leads to significant least squares bias, averaging can only improve this to Õ(d ω / √ ϵ) (see Table 1 ). An alternate strategy would be to combine leverage score sampling with a preconditioned mini-batch stochastic gradient descent (Weighted Mb-SGD), with mini-batches chosen so that they fit in Õ(d 2 ) space. This achieves the same time and space complexity as our method, but due to the streaming access to A and the sequential nature of SGD, it requires O(1/ϵ) data passes. Instead, our algorithm essentially mixes an Õ(d/ϵ) size leverage score sample into an O(d) size sketch, merging Õ(1/ϵ) rows of A into a single row of the sketch (see Figure 1 ). This results in better data compression compared to direct leverage score sampling, with only Õ(d 2 ) bits of space, while retaining the same Õ(nnz(A) + d 2 /ϵ) runtime complexity as the above approaches and requiring only a single data pass. The resulting estimator x can no longer recover the ϵ-error bound, but remarkably, its expectation E[x] still does. To turn this into an improved estimator in a distributed model, we can simply average q = 1/ϵ such estimators, i.e., x = 1 q q i=1 xi , obtaining E∥Ax -b∥ 2 ≤ (1 + 2ϵ)∥Ax * -b∥ 2 . As shown in Table 1 , ours is the first result in this model to achieve Õ(d 2 ) space in a single pass and faster than current matrix multiplication time O(d ω ).Finally, by incorporating a preconditioning scheme, we illustrate how our construction can be used to design the first algorithm that solves least squares in current matrix multiplication time, constant parallel passes and O(d 2 log(nd)) bits of space. We note that the O(d ω ) cost comes only from the worst-case complexity of constructing the preconditioner P, which can often be accelerated in practice. The computational model used in Theorem 2 is described in detail in Section 3. Theorem 2. Given A ∈ R n×d and b ∈ R n in the parallel computing model, using two parallel passes with q machines, we can compute x such that with probability 0.9∥Ax -b∥ ≤ 1 + ϵ + O(1/q) ∥Ax * -b∥ in O(γ -1 nnz(A) + d ω + ϵ -1 d 2+γ polylog(d)) time, O(d 2 log(nd)) bits of space and O(d log(nd)) bits of communication. In particular, choosing q = 1/ϵ we recover an O(ϵ)-approximation. Remark 2. While the parallel computing model in Theorem 2 assumes that all machines have streaming access to the entire data matrix, this result can be easily extended to the setting where A has been randomly down-sampled or partitioned into separate size n chunks A 1 , ..., A q , and each machine constructs an estimate xi based on a sketch of its own chunk. Then, with the same computational guarantees as in Theorem 2, the averaged estimator x = 1 q q i=1 xi with probability 0.9 enjoys a guarantee of:∥Ax -b∥ ≤ 1 + O(ϵ + 1/q + B chunk ) ∥Ax * -b∥,where B chunk is the bias that would be incurred if we solved each chunk exactly and averaged those solutions. Using existing guarantees for uniform down-sampling of least squares [ Theorem 20, 43] , one can bound this bias asB chunk = Õ µ qn + ( µ n ) 2, where µ is the coherence of the data matrix and n is the chunk size. For sufficiently large chunks, this bias is negligible compared to the error ϵ. We prove the above high probability guarantee in Theorem 6 in Appendix E.Further applications. Our least squares analysis can be extended to other settings where prior works [e.g., 17, 18, 16, 22] have analyzed randomized estimators based on sparse sketching via techniques from asymptotic random matrix theory. The primary and most direct application involves correcting inversion bias in the so-called sketched inverse covariance estimate ( Ã⊤ Ã) -1 , which was the motivating task of [18] , with applications including distributed second-order optimization and statistical uncertainty quantification, where quantities like ( Ã⊤ Ã) -1 x are approximated. This result should be compared with ϵ = Õ (1 + d/s) √ d m obtained by [18] . Thus, we get a direct improvement for very sparse sketches, i.e., s = o(d). This can be immediately translated into an improved local convergence guarantee for Distributed Newton Sketch which is a second-order convex minimization algorithm used in settings where the Hessian matrix can be expressed as A ⊤ A for a tall matrix A, e.g., in generalized linear models like logistic regression. In this method, following Corollary 16 of [18] as well as related results [16, 44, 19] , we use sketching and averaging to estimate a Newton step:x t+1 = x t - 1 q q i=1 H-1 i g t ,where g t is the gradient at x t and H1 , ..., Hq are Hessian sketches constructed by independent machines. The following corollary, which is a direct improvement over Corollary 16 of [18] , shows that Distributed Newton Sketch on a generalized linear model task can achieve a fast local convergence rate of the formf (xt+1)-f (x * ) f (xt)-f (x * ) = o(1) with O(d ω ) time per iteration (see Appendix B.1 for details). Corollary 1 (Distributed Newton Sketch). Consider f (x) = 1 n n i=1 ℓ i (x ⊤ ϕ i ) + λ 2 ∥x∥ 2, where ℓ i are convex twice continuously differentiable functions, such that f has a Lipschitz Hessian, λ > 0, and ϕ ⊤ i is the ith row of an n × d data matrix Φ. Given ϵ > 0, there is a neighborhood U ϵ around the minimizer x * = argmin x f (x) such that, for any x t ∈ U ϵ , using two parallel passes with Õ(1/ϵ) machines, we can compute a Distributed Newton Sketch update x t+1 such thatf (x t+1 ) -f (x * ) ≤ ϵ • f (x t ) -f (x * ) , in O(γ -1 nnz(Φ) + d ω + ϵ -1 d 2+γ polylog(d)) time, O(d 2 log(nd)) bits of space and O(d log(nd)) bits of communication.Our Techniques. At the core of our analysis are techniques inspired by asymptotic random matrix theory (RMT) in the proportional limit [e.g., see 6] . Here, in order to establish the limiting spectral distribution (such as the Marchenko-Pastur law) of a random matrix Ã⊤ Ã whose dimensions diverge to infinity, one aims to show the convergence of the Stieltjes transform of its resolvent matrix ( Ã⊤ Ã -zI) -1 . Recently, [18] showed that these techniques can be adapted to sparse sketching matrices (via leverage score sparsification) in order to characterize the bias of the sketched inverse covariance ( Ã⊤ Ã) -1 , where Ã = SA.Our main contribution is two-fold. First, we show that a similar argument can also be applied to analyze the bias of the least squares estimator, x = ( Ã⊤ Ã) -1 Ã⊤ b. Unlike the inverse covariance, this estimator no longer takes the form of a resolvent matrix, but its bias is also associated with the inverse, which means that we can use a leave-one-out argument to characterize the effect of removing a single row of the sketch on the estimation bias. Our second main contribution is to improve the sharpness of the bounds relative to the sparsity of the sketching matrix by combining a careful application of Hölder's inequality with a higher moments analysis of the restricted Bai-Silverstein inequality for quadratic forms. Those improvements are not only applicable to the least squares analysis, but also to all existing RMT-style results for LESS embeddings, including the aforementioned inverse covariance estimation, as well as applications in stochastic optimization, resulting in the sketching cost of LESS embeddings dropping below matrix multiplication time. 
related work: 
 Randomized numerical linear algebra. RandNLA sketching techniques have been developed over a long line of works, starting from fast least squares approximations of [42] ; for an overview, see [45, 26, 37, 39, 21] among others. Since then, these methods have been used in designing fast algorithms not only for least squares but also many other fundamental problems in numerical linear algebra and optimization including low-rank approximation [13, 34] , l p regression [14] , solving linear systems [28, 24] and more. Using sparse random matrices for matrix sketching also has a long history, including data-oblivious sketching methods such as CountSketch [12] , OSNAP [40] , and more [38] . Leverage score sparsification (LESS) was introduced by [18] as a data-dependent sparse sketching method to enable RMT-style analysis for sketching (see below).Unbiased estimators for least squares. To put our results in a proper context, let us consider other approaches for producing near-unbiased estimators for least squares, see also Table 1 . First, a well known folklore result states that the least squares estimator computed from a dense Gaussian sketching matrix is unbiased. The bias of other sketching methods, including leverage score sampling and OSNAP, has been studied by [43] , showing that these methods need a √ ϵ-error guarantee to achieve an ϵ-bias which leads to little improvement unless ϵ is extremely small and the sketch size is sufficiently large. Another approach of constructing unbiased estimators for least squares, first proposed by [23] , is based on subsampling with a non-i.i.d. importance sampling distribution based on Determinantal Point Processes [DPPs, 31, 20] . However, despite significant efforts [30, 7, 5] , sampling from DPPs remains quite expensive: the fastest known algorithm requires running a Markov chain for polylog(n/ϵ) many steps, each of which requires a separate data pass and takes O(d ω ) time. Other approaches have also been considered which provide partial bias reduction for i.i.d. RandNLA subsampling schemes in various regimes that are are either much more expensive or not directly comparable to ours [2, 44] .Statistical and RMT analysis of sketching. Recently, there has been significant interest in statistical and random matrix theory (RMT) analysis of matrix sketching; see [21] for an overview. These approaches include both asymptotic analysis via limiting spectral distributions and deterministic equivalents [35, 25, 32, 33] , and non-asymptotic analysis under statistical assumptions [36, 41, 4] . A number of works have shown that the RMT-style techniques based on deterministic equivalents can be made rigorously non-asymptotic for certain sketching methods such as dense sub-Gaussian [17] , LESS matrices [16, 18] , and other sparse matrices [9] , which has been applied to low-rank approximation, fast subspace embeddings and stochastic optimization. Our new analysis can be viewed as a general strategy for directly improving the sparsity required by LESS embeddings (and thereby, the sketching time complexity) in many of these applications, specifically those that rely on analysis inspired by the calculus of deterministic equivalents via generalized Stieltjes transforms.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [NA]Justification: The results proved in the paper hold in a very general setting, without any limitations or requiring strong assumptions.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
preliminaries: 
 Notations. In all our results, we use lowercase letters to denote scalars, lowercase boldface for vectors, and uppercase boldface for matrices. The norm ∥ • ∥ denotes the spectral norm for matrices and the Euclidean norm for vectors, whereas ∥ • ∥ F denotes the Frobenius norm for matrices. We use ⪯ to denote the p.s.d. ordering of matrices.Computational model. We next clarify the computational model that is used in Theorem 2. We consider a central data server storing (A, b), and q machines. The jth machine has a handle Stream(j), which can be used to open a stream and to read the next row/label pair (a i , b i ) in the stream. After a full pass, the machine can re-open the handle and begin another pass over the data. The machines can operate their streams entirely asynchronously, and each has its own limited local storage space, e.g., in Theorem 2 we use O(d 2 log(nd)) bits of space per machine. At the end, they can communicate some information back to the server, e.g., in Theorem 2, they communicate their final estimate vectors xj , using O(d log(nd)) bits of communication. Then, the server computes the final estimate, in our case via averaging, x = 1 q q i=1 xi , which can be done either directly or via a map-reduce type architecture.We define the parallel passes required by such an algorithm as the maximum number of times the stream is opened by any single machine. We analogously define time/space/communication costs by taking a maximum over the costs required by any single machine (for communication, this refers only to the number of bits sent from the machine back to the server).Definitions and useful lemmas. In our framework, we construct a sparse sketching matrix S where sparsification is achieved using a probability distribution over rows of data matrix A, that is proportional to the leverage scores of A. The next definition [following, e.g., 9] provides the explicit definition of exact and approximate leverage scores for our setting. Definition 1 ((β 1 , β 2 )-approximate leverage scores). Fix a matrix A ∈ R n×d and consider matrix U ∈ R n×d with orthonormal columns spanning the column space of A. Then, the leverage scores l i , 1 ≤ i ≤ n are defined as the row norms squared of U, i.e., l i = ∥u i ∥ 2 , where u ⊤ i is the ith row of U. Furthermore, consider fixed β 1 , β 2 > 1. Then li are called (β 1 , β 2 )-approximate leverage scores for A if the following holds for all i l i β 1 ≤ li andn i=1 li ≤ β 2 • d.The approximate leverage scores can be computed by first constructing a preconditioner matrix P ∈ R d×d such that κ(AP) = O(1), which takes O(nnz(A) + d ω ) in a single pass, and then relying on the following norm approximation scheme. Lemma 1 (Based on Lemma 7.2 from [10] ). Given A ∈ R n×d and P ∈ R d×d , using a single pass over A in time O(γ -1 (nnz(A) + d 2 )) for small constant γ > 0, we can compute estimates l1 , ..., ln such that with probability ≥ 0.95:n -γ ∥e ⊤ i AP∥ 2 ≤ li ≤ O(log(n))∥e ⊤ i AP∥ 2 ∀i and i li ≤ O(1) • ∥AP∥ 2 F .In the next definition, we give the sparse sketching strategy used in our analysis. This approach is similar to the original leverage score sparsification proposed by [18] , except: 1) we adapted it so that it can be implemented effectively in a single pass, and 2) we use it in a much sparser regime (fewer non-zeros per row). Definition 2 ((s, β 1 , β 2 )-LESS embedding). Fix a matrix A ∈ R n×d and some s ≥ 0. Let the tuple ( l1 , • • • , ln ) denote (β 1 , β 2 )-approximate leverage scores for A. Let p i = min{1, sβ1 li d }. We define a (s, β 1 , β 2 )-approximate leverage score sparsifier ξ as follows.ξ = b 1 √ p 1 , • • • , b n √ p n where b i ∼ Bernoulli(p i ).Moreover, we define the (s, β 1 , β 2 )-leverage score sparsified (LESS) embedding of size m as matrix S ∈ R m×d with i.i.d. rows 1 √ m x i such that x i = diag(ξ i )y i where ξ i denotes a randomly generated (β 1 , β 2 )-approximate leverage score sparsifier and y i ∈ R n consist of random ±1 entries.A key property of a sketching matrix is the subspace embedding property, defined below. It was recently shown by [9] that LESS embeddings require only polylogarithmically many non-zeros per row of S to prove that S is a subspace embedding for the data matrix A with the optimal m = O(d) sketching dimension. For our main results, it is sufficient to use η = O(1) below. Lemma 2 (Subspace embedding for LESS, Theorem 1.3, [9] ). Fix η, δ > 0. Consider β 1 , β 2 > 1 and a full rank matrix A ∈ R n×d . Then for a (β 1 , β 2 )-leverage score sparsified embedding S ∈ R m×n with s ≥ O(log 4 (d/δ)/η 4 ) and m = O((d + log 1/δ)/η 2 ), with probability 1 -δ we haveEQUATIONQuestion: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
least squares bias analysis: 
 In this section we provide an outline of the bias analysis for the sketched least squares estimator constructed using a LESS embedding, leading to the proofs of our main results, Theorems 1 and 2. In particular, we prove the following main technical result (detailed proof in Appendix C). Theorem 4 (Bias of LESS-sketched least squares). Fix A ∈ R n×d and let S be an (s, β 1 , β 2 )-LESS embedding of size m for A. Let S satisfy (1) with η = 1 2 and probability 1 -δ where δ < 1 m 4 . Then there exists an event E with probability at least 1 -δ such thatL(E E [x]) -L(x * ) = O d m 2 1 + d s log 9 (n/δ) • L(x * ).Remark 3. Thus, the biasL(E E [x]) -L(x * ) of the LESS estimator using O(β 1 β 2 s) = Õ(s) non- zeros per row of S is of the order Õ( d 2 sm 2 + d m 2 ) • L(x * ). By comparison, the standard expected loss bound which holds for sketched least squares (including this estimator) isE[L(x)] -L(x * ) ≤ Õ( d m )L(x *), and the best known bound on the bias of most standard sketched estimators (e.g., leverage score sampling) is Õ( d 2 m 2 )L(x * ), given by [43] . So, our result recovers the standard bias bound for s = 1 and improves on it for s ≫ 1 by a factor of min{s, d}. At the end of the section, we discuss how to deal with the lower order term Õ( d m 2 ) to reduce the bias further.Proof sketch. Using a standard argument, we can replace the matrix A with the matrix U ∈ R n×d consisting of orthonormal columns spanning the column space of A, and assume that n = poly(d).Let S be an (s, β 1 , β 2 )-LESS embedding for U. Also, let b ∈ R n be a vector of responses/labels corresponding to n rows in U. Let x = argmin x ∥SUx -Sb∥ 2 . Furthermore for any x ∈ R d we can find the loss at x as L(x) = ∥Ux -b∥ 2 . Additionally, we use r to denote the residual b -Ux * . We also define Q = (γU ⊤ S ⊤ SU) -1 as the sketched inverse covariance matrix with scaling γ = m m-d representing the standard correction accounting for inversion bias. We condition on the high probability event E guaranteed in Lemma 2 and consider L(E E (x))-L(x * ). By Pythagorean theorem, we haveL(E E [x]) -L(x * ) = ∥U(E E [x]) -Ux * ∥ 2 .Note that by the normal equations we have x = (U ⊤ S ⊤ SU) -1 U ⊤ S ⊤ Sb = γQU ⊤ S ⊤ Sb, and alsoS ⊤ S = 1 m m i=1 x i x ⊤i . These two facts lead to writing the bias as follows:L(E E [x]) -L(x * ) = ∥γ • E E [QU ⊤ x i x ⊤ i r]∥ 2 .Using a leave-one-out technique , we replace Q withQ -i = (γU ⊤ S ⊤ -i S -i U) -1, where S -i denotes matrix S without the ith row, by noting thatQ = (γU ⊤ S ⊤ -i S -i U + γ m U ⊤ x i x ⊤ i U) -1and applying the Sherman-Morrison formula. This leads to the following relation:L(E E [x]) -L(x * ) ≤ 2 ∥E E [Q -i U ⊤ x i x ⊤ i r]∥ 2 ∥Z0r∥ 2 +2 E E γ γ i -1 Q -i U ⊤ x i x ⊤ i r 2 ∥Z2r∥ 2whereγ i = 1 + γ m x ⊤ i UQ -i U ⊤ x i .Due to the subspace embedding assumption and assuming m large enough, we have ∥Q∥ = O(1) and also ∥Q -i ∥ = O(1). The first term ∥Z 0 r∥ 2 is quite straightforward to bound since, if not for the conditioning on the high probability event E, we would haveE[Q -i U ⊤ x i x ⊤ i r] = E[Q -i U ⊤ r] = 0, which follows from U ⊤ (b -Ux * ) = 0. We get an upper bound on ∥Z 0 r∥ 2 as O d 2 log(d/δ) sm 2 + d m 2 • ∥r∥ 2 , which is sufficient for us.The central novelty of our analysis lies in bounding ∥Z 2 r∥ 2 for (s, β 1 , β 2 )-LESS embeddings, which is the dominant term. Our key observation is that, when examining a random variable of the form x ⊤ i v for some vector v, the dependence on the sparsity of row x i only arises when considering moments higher than 2 + 1 O(log(n)) , because otherwise we can simply rely on the fact that E[x i x ⊤ i ] = I. Thus, when decomposing ∥Z 2 r∥ 2 , we must carefully separate the contribution of near-second moments vs the contribution of higher moments to the overall bound.To obtain this separation, we start by applying Hölder's inequality on ∥Z 2 r∥ with p = O(log(n))and q = 1 + 1 O(log(n)) to get ∥Z 2 r∥ ≤ E E [ γ γ i -1 p ] 1/p • sup ∥v∥=1 E E [v ⊤ Q -i U ⊤ x i x ⊤ i r] q 1/q .Furthermore applying Cauchy-Schwarz inequality on the second term leads to∥Z 2 r∥ ≤ E E [ γ γ i -1 p ] 1/p • E E ∥Q -i U ⊤ x i ∥ 2q 1/2q • E E ∥x ⊤ i r∥ 2q 1/2q.Unlike [18] , we exploit the fact that∥x i ∥ 1/O(log(n)) = O(1)and get a constant upper bound on(E E ∥Q -i U ⊤ x i ∥ 2q). However, this results in a much more careful argument, requiring now an upper bound onEQUATIONwhereγ = 1 + γ m E E (x ⊤ i UQ -i U ⊤ x i ).In particular, for the second term, we haveEQUATIONTo bound the first of these two terms, we prove a new version of the Restricted Bai-Silverstein inequality (Lemma 3) for (s, β 1 , β 2 )-LESS embeddings. Unlike [18] , we provide a proof with any p and any (β 1 , β 2 ) values. Furthermore, utilizing the subspace embedding guarantee from Lemma 2, we prove a much more general result where the number of non-zeros in the approximate leverage score sparsifier ξ can be much smaller than d (proof in Appendix D). Lemma 3 (Restricted Bai-Silverstein for (s, β 1 , β 2 )-LESS embeddings). Let p ∈ N be fixed and U ∈ R n×d be such that U ⊤ U = I. Let x i = diag(ξ)y i where y i ∈ R n has independent ±1 entries and ξ is an (s, β 1 , β 2 )-approximate leverage score sparsifier for U. Then for any matrix with 0 ⪯ C ⪯ O(1) • I and any δ > 0 we have for an absolute constant c > 0.E [tr(C) -x ⊤ i UCU ⊤ x i ] p 1/p < c • √ dp 3 • 1 + dp log(d/δ) s .Using Lemma 3, we upper bound the first term squared in (3) asÕ d m 2 1 + d s. Moreover, also using Lemma 3, we get a matching upper bound on |γ -γ|. We then design a martingale concentration argument to prove a high probability upper bound on the last remaining term, |tr(Q -i )-E E tr(Q -i )|,which implies the desired moment bound (proof in Appendix B), concluding the proof of Theorem 4. Lemma 4. For given δ > 0 and matrix Q -i we have with probability 1 -δ:|tr(Q -i ) -Etr(Q -i )| ≤ c ′ γ • d √ m log 4.5 (m/δ).Completing the proof of Theorem 1. First, suppose that ϵ ≥ O(polylog(d)/d) so that the bias bound can be achieved from Theorem 4. Our implementation is mainly based on the online construction of approximate leverage scores, given the preconditioner P, using Lemma 1. Briefly, this construction proceeds by first sketching P using a d × O(1/γ) Gaussian matrix G to produce the matrix P = PG, and then, for each observed row a i of A, we compute li = ∥a ⊤ i P∥ 2 . Assuming without loss of generality that d = poly(n) and adjusting γ, the estimates satisfyβ 1 β 2 = O(αd γ ).Next, we sample the non-zero entries of S corresponding to the observed row a i , i.e., the i-th column of S. Note that for this we only need to know the single leverage score estimate li . Crucially for our analysis, the entries of this column need to be sampled i.i.d., which can be done in time proportional to the number of non-zeros in that column by first sampling a corresponding Binomial distribution to determine how many non-zeros we need, then picking a random subset of that size, and then sampling the random ±1 values. Altogether, the cost of constructing the sketch isO(γ -1 nnz(A)+β 1 β 2 sd 2 ) = O(γ -1 nnz(A)+αϵ -1 d 2+γ polylog(d)) by setting s = O(polylog(d)/ϵ).Finally, once we construct the sketch, at the end of the pass we can run conjugate gradient preconditioned with P on the sketched problem, which takes Õ(αd 2 ).We note that in the (somewhat artificial) regime where we require extremely small bias, i.e., ϵ = o(polylog(d)/d), the bound claimed in Theorem 1 can still be obtained, since in this case for small enough γ we haved 2+γ /ϵ = O(d ω / √ ϵ)with ω < 2.5, so we can rely on direct leverage score sampling (which corresponds to s = 1), and instead of maintaining the sketch, we compute the estimator x = ( Ã⊤ Ã) -1 Ã⊤ b directly along the way. This involves performing a separate d × d matrix multiplication after collecting each d leverage score samples, to gradually compute Ã⊤ Ã, and then inverting the matrix at the end. From Theorem 4, we see that it suffices to set sketch size m = Õ(d/ √ ϵ), which leads to the desired runtime.Completing the proof of Theorem 2. Here, we use a slightly modified variant of Lemma 2, given as Theorem 1.4 in [9] , which shows that using a single pass we can compute a sketch Ã in time O(nnz(A) + d ω ), which satisfies the subspace embedding property (1) with η = 1 2 . Then, we can perform the QR decomposition Ã = QR and set P = R -1 in additional time O(d ω ) to obtain the desired preconditioner. Next, we use Theorem 1 to construct q i.i.d. estimators xi in a second parallel pass, and finally, the estimators are aggregated to compute x = 1 q q i=1 xi , which satisfies E∥Ax -b∥ 2 ≤ 1 + ϵ + O(1/q) ∥Ax * -b∥ 2 . Applying Markov's inequality concludes the proof. 
experiments: 
 In this section, we illustrate empirically how our results point to a practical free lunch phenomenon in distributed averaging of sketching-based estimators. As mentioned in Section 1, our construction from Theorem 1 essentially works by taking a subsample of the data and then mixing groups of those rows together to produce an even smaller sketch (see Figure 1 ). According to our theory, while the small sketch does not recover the same ϵ-small error as the larger subsample, it does recover an ϵ-small bias. Moreover, this happens without incurring any additional computational cost, as the cost of the sketching is proportional to the cost of simply reading the subsampled rows. This suggests that we can use sparse sketching to compress a data subsample down to a small size while retaining the least squares performance in a distributed averaging environment.To verify this, we evaluate the effectiveness of distributed averaging of sketched least squares estimators on several benchmark datasets. Specifically, we visualize the relative error of the averaged sketch-and-solve estimatorL(x)-L(x * ) L(x * ), against the number of machines q used to generate the estimate x = 1 q q i=1 xi . Each estimate xi is constructed with the same sparsification strategy used by LESS, except that instead of sparsifying the sketch with leverage scores, we instead sparsify them with uniform probabilities (which is often sufficient in practice). Following [16] , we call the resulting method LESSUniform. Within each dataset, we perform four simulations, designed so that the total sketching cost stays the same for all four test cases, by simultaneously changing sketch size and sparsity. Concretely, we vary these so that the product (sketch size × nnz per row) stays the same in each case, so as to ensure that the total cost of sketching is fixed in each plot.In Figure 2 , on the X-axis we plot the number q of estimators being averaged, so that the bias of a single estimator appears on the right-hand side of the plot (large q), whereas the variance (error) appears on the left-hand side (q = 1). In each plot, the line with nnz per row = 1 (and large sketch size) corresponds to uniform subsampling, whereas the remaining ones are sketches produced by compressing that subsample. The plot shows that decreasing the sketch size (i.e., compressing the sample) does increase the error of a single estimator (as expected), however it also shows that the bias of these estimators remains essentially unchanged regardless of the sketch size (since all lines showing that sparse sketching can be used to compress the data while preserving near-unbiasedness without increasing the estimation cost (see Appendix F for results on the Boston dataset). meet as q → ∞), confirming that suitable sparse sketches that "compress" rows of data can preserve near-unbiasedness without increasing the cost. This phenomenon may not occur for all sketching methods. Figure 4 within Appendix F showcases a few more interesting results. We first further demonstrate that suitable sketches that "compress" rows of data A into a single row of Ã can preserve near-unbiasedness without increasing the cost. In particular, this desirable phenomenon that LESSUniform enjoys also extends to LESS proper, as well as the Gaussian and Subgaussian (Rademacher) sketches. In fact, we observe that LESS enjoys similar desirable performance as the Gaussian and Subgaussian sketches and virtually no least squares bias, while retaining the computational speedups of sparse sketching, suggesting that it attains the best of both worlds.However, when we decrease the number of subsamples within leverage score subsampling, the bias introduced by subsampling increases as expected. This happens as the number of subsamples is reduced without increasing the amount of "compression" as one would with LESS or LESSUniform. 1We also show that the subsampled randomized Hadamard transform (SRHT) can exhibit some amount of least squares bias as the sketch size decreases. The numerical results shown for the bias introduced by leverage score subsampling and the SRHT complement the lower bounds established in [18] . 
conclusions: 
 We gave a new sparse sketching method that, using two passes over the data, produces a nearlyunbiased least squares estimator, which can be used to improve upon the space-time trade-offs of solving least squares in parallel or distributed environments via simple averaging. In particular, our algorithm is the first to require only O(d 2 log(nd)) bits of space and current matrix multiplication time O(d ω ) while obtaining an ϵ = o(1) approximation in few passes. Our techniques are of broader interest to sketching-based optimization algorithms, including Distributed Newton Sketch. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: All experiments were run on an i9-13900k with 128GB of RAM and a NVIDIA RTX 3090 GPU.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: This work aligns with the NeurIPS code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: Our theoretical work provides a Matrix sketching algorithm, recovering a nearly unbiased least squares estimator. Solving linear systems approximately is an important component of a large number of algorithms from varied domains. Our work can thus provide a way to recover approximation solutions with a smaller bias. This does not lead to any negative societal impact.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This work does not release any data or model, and carries no risk of misuse.• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We test our theoretical results on a few publicly available datasets that have been cited in the paper as well.Guidelines:• The answer NA means that the paper does not use existing assets.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: We start by providing several classical results, used in our analysis. The following formula provides a way to compute the inverse of matrix A after a rank-1 update, given the inverse before the update.A ∈ R d×d and vector u, v ∈ R d , A + uv ⊤ is invertible if and only if 1 + v ⊤ A -1 u ̸ = 0. If this holds then, (A + uv ⊤ ) -1 = A -1 - A -1 uv ⊤ A -1 1 + v ⊤ A -1 u .In particular,(A + uv ⊤ ) -1 u = A -1 u 1 + v ⊤ A -1 u .The following inequality provides a crucial tool for writing expectation of the product of two random variables as the product of higher individual moments.Lemma 6 (Hölder's inequality). For real-valued random variables X and Y ,E[|XY |] ≤ (E[|X| p ]) 1/p • (E[|Y | q ]) 1/qwhere p, q > 0 are Hölder's conjugates, i.e. 1 p + 1 q = 1.The following technical lemmas provide concentration results for the sum of random quantities. We collect these results here and then refer to them while using in our analysis.Lemma 7 (Matrix Chernoff Inequality).For i = 1, 2, • • • , n consider a sequence Z i of d × d positive semi-definite random matrices such that E[ 1 n i Z i ] = I d and ∥Z i ∥ ≤ R.Then for any ϵ > 0, we havePr λ max 1 n n i=1 Z i ≥ (1 + ϵ) ≤ d • exp - nϵ 2 (2 + ϵ)R .where λ max denotes taking the maximum eigenvalue.Lemma 8 (Azuma's inequality).If {Y 0 , Y 1 , Y 2 , • • • } is a martingale with |Y j -Y j-1 | ≤ c j then for any m, λ > 0 we have Pr (|Y m -Y 0 | ≥ λ) ≤ 2 • exp - λ 2 2 m j=1 c 2 j .Lemma 9 (Rosenthal's inequality ([29], Theorem 2.5 and Corollary 2.6)). Let 1 ≤ p < ∞ and X 1 , X 2 , • • • , X n are nonnegative, independent random variables with finite p th moments then,E i X i p 1/p ≤ 2p log(p) • max    i E[X i ], i E[X p i ] 1/p    .Furthermore, for mean-zero independent and symmetric random variables we haveE i X i p 1/p ≤ 2p log(p) • max    i E[X 2 i ] 1/2 , i E[X p i ] 1/p    .Lemma 10 (Bai-Silverstein's Inequality Lemma B.26 from [6] ). Let B be a d × d be a fixed matrix and x be a random vector of independent entries. Let E[x i ] = 0 and E[x 2 i ] = 1,and E|x j | l ≤ ν l . Then for any p ≥ 1,E|x ⊤ Bx -tr(B)| p ≤ (2p) p • (ν 4 tr(BB ⊤ )) p/2 + ν 2p tr(BB ⊤ ) p/2 .We use the notion of an (ϵ, δ) unbiased estimator of A as defined in [18] . Definition 3 ((ϵ, δ)-unbiased estimator). For ϵ, δ > 0, a random positive definite matrix B ∈ R d×d is called an (ϵ, δ) unbiased estimator of A if there exists an event E with Pr(E) ≥ 1 -δ such that,1 1 + ϵ A ⪯ E E [B] ⪯ (1 + ϵ)A and, B ⪯ O(1) • A,when conditioned on the event E.In this section, we give a formal statement and proof for Theorem 3, which is then used in the proof of Theorem 4. We replace A with U such that U consists of d orthonormal columns spanning the column space of A. Here S m ∈ R m×n denotes a LESS sketching matrix with independent rows1 √ m x ⊤ , x ⊤ = y ⊤ • diag(ξ)where y consists of ±1 Rademacher entries and ξ is an (s, β 1 , β 2 )approximate leverage score sparsifier. Note that E[xx ⊤ ] = I n . We assume that the sketching matrix S m consists of m ≥ 10d i.i.d rows and 3 divides m. Also, we assume that the S m satisfies the subspace embedding condition for U (Theorem 2) with η= 1 2 . Let Q = (γU ⊤ S ⊤ m S m U) -1 where γ = m m-d . Theorem 5 (Small inversion bias for (s, β 1 , β 2 )-LESS embeddings). Let δ > 0 satisfy δ < 1 m 4 and m ≥ O(d). Let S m ∈ R m×n be an (s, β 1 , β 2 )-LESS embedding for data matrix U ∈ R n×d such that U ⊤ U = I.Then there exists an event E with Pr(E) ≥ 1 -δ such that,1 1 + ϵ • I ⪯ E E [Q] ⪯ (1 + ϵ) • I and 1 2 I ⪯ Q ⪯ 2I when conditioned on E where ϵ = O √ d log 4.5 (n/δ) m 1 + d s .Proof. Let S -i denote S m without the i th row, and S -ij denote S m with the i th and j th rows removed. LetQ -i = (γU ⊤ S ⊤ -i S -i U) -1 and Q -ij = (γU ⊤ S ⊤ -ij S -ij U) -1 .We proceed with the same proof strategy as adopted in [18] . We define the events E j as follows:E j = 3 m U ⊤   tj i=t(j-1)+1 x i x ⊤ i   U ⪰ 1 2 I, j = 1, 2, 3 , E = ∧ 3 j=1 E j .Note that event E j means that the sketching matrix with just (1/3) rd rows (scaled to maintain unbiasedness of the sketch) from S m satisfies a lower spectral approximation of U ⊤ U = I. Also we notice that events E 1 , E 2 , E 3 are independent, and for any pair (i, j) there exists at least one event E k , k ∈ {1, 2, 3} such that E k is independent of both x i and x j . Furthermore conditioned on E k we haveQ -i ⪯ 6 • I d and Q -ij ⪯ 6 • I d .Note that as guaranteed in Theorem 2, we have Pr(E k ) ≥ 1 -δ ′ for all k and therefore Pr(E) ≥ 1 -δ with δ ′ = δ/3. Let E E denote the expectation conditioned on the event E.I -E E [Q] = -E E [Q] + γE E [QU ⊤ S ⊤ m S m U] = -E E [Q] + γE E [QU ⊤ x i x ⊤ i U] = -E E [Q] + γE E Q -i U ⊤ x i x ⊤ i U 1 + γ m x ⊤ i UQ -i U ⊤ x i = E E [Q -i U ⊤ (x i x ⊤ i -I)U] Z0 + E E [Q -i -Q] Z1 + E E γ γ i -1 Q -i U ⊤ x i x ⊤ i U Z2whereγ i = 1 + γ m x ⊤ i UQ -i U ⊤ x i .The second equality follows by noting that S ⊤ m S m = 1 m m i=1 x i x ⊤ i and using linearity of expectation. The third equality holds due to the application of Sherman Morrison's (Lemma 5) formula onQ = (γU ⊤ S ⊤ -i S -i U + γ m U ⊤ x i x ⊤ i U) -1. We start by upper bounding Z 0 and use the following from [18] .Lemma 11 (Upper bound on ∥Z 0 ∥). For any k > 0 we have,∥E E [Q -i U ⊤ (x i x ⊤ i -I)U]∥ ≤ 12 kδ ′ + ∞ k Pr(x ⊤ i UU ⊤ x i ≥ x)dx .Using Chebyshev's inequality we have,Pr(x ⊤ i UU ⊤ x i ≥ x) ≤ Var(x ⊤ i UU ⊤ xi) x 2. Using Restricted Bai-Silverstein inequality for (s, β 1 , β 2 )-approximate LESS embeddings i.e., Lemma 3 (proved in Lemma 16 ) with p = 2, we have Var(x ⊤ i UU ⊤ x i ) ≤ cd • 1 + d log(d/δ)s for some absolute constant c. Let k = m 2 and δ ′ < 1 m 4 we get,EQUATION)We use the bound on the term ∥Z 1 ∥ directly from [18] , provided below as a Lemma.Lemma 12 (Upper bound on ∥Z 1 ∥, [18] ).EQUATIONIt remains to upper bound ∥Z 2 ∥.∥Z 2 ∥ = ∥E E γ γ i -1 Q -i U ⊤ x i x ⊤ i U ∥ ≤ sup ∥v∥=1, ∥z∥=1 E E [| γ γ i -1| • |v ⊤ Q -i U ⊤ x i x ⊤ i Uz|].Applying Hölder's inequality with p = O(log(n)) and qEQUATIONwhere we used Cauchy-Schwarz inequality on E E |v ⊤ Q -i U ⊤ x i x ⊤ i Uz| q . Now note that x i = diag(ξ) • y i , we have ∥x i ∥ = poly(n) and therefore ∥x i ∥ ∆ = O(1). We now show the terms involving exponents depending on q are O(1) as highlighted in the inequality (6) .sup ∥v∥=1 E E |v ⊤ Q -i U ⊤ x i | 2q 1/2q ≤ sup ∥v∥=1 E E |v ⊤ Q -i U ⊤ x i | 2 • v ⊤ Q -i U ⊤ x i | 2∆ ] 1/2q ≤ sup ∥v∥=1 E E |v ⊤ Q -i U ⊤ x i | 2 • v ⊤ Q -i U ⊤ ∥ 2∆ • ∥x i ∥ 2∆ ] 1/2q ≤ O(1) • E E ∥Q -i U ⊤ x i ∥ 2 • ∥Q -i U ⊤ ∥ 2∆ 1/2q ≤ O(1) • 2 • E E ′ ∥Q -i U ⊤ x i ∥ 2 • ∥Q -i U ⊤ ∥ 2∆ 1/2q .Here E ′ is an event independent of x i . Without loss of generality, we can assume that E ′ = E 1 ∧ E 2 . We first condition on Q -i and take expectation over x i . Also note that Q -i and x i are independent and event E ′ is independent of x i , and furthermoreE E ′ [x ⊤ i x i ] = E[x ⊤ i x i ] = 1. We get, sup ∥v∥=1 E E |v ⊤ Q -i U ⊤ x i | 2q 1/2q ≤ O(1) • E E ′ [∥Q -i U ⊤ ∥ 2q ] 1/2q .Now we use that conditioned on E ′ , ∥Q -i ∥ ≤ 6 and ∥U ⊤ ∥ = 1, we get,sup ∥v∥=1 E E |v ⊤ Q -i U ⊤ x i | 2q 1/2q= O(1).i ∥ ∆ = O(1) and E E ′ [x ⊤ i x i ] = E[x ⊤ i x i ] = 1, sup ∥z∥=1 E E |x ⊤ i Uz| 2q 1/2q = O(1).E γ γi -1 p 1/p. Without loss of generality, we assume that p is even. We have,E E γ γ i -1 p ≤ 2 • E E ′ γ γ i -1p where E ′ is event independent of x i . The above can be upper bounded as,EQUATIONwhereγ = 1 + γ m E E ′ (x ⊤ i UQ -i U ⊤ x i ). As E ′ is independent of x i and Q -i is independent of x i , we get E E ′ (x ⊤ i UQ -i U ⊤ x i ) = E E ′ tr(Q -i ). Therefore, γ = 1 + γ m E E ′ tr(Q -i ). We now aim to upper bound (E E ′ [(γ i -γ) p ])1/p as,(E E ′ [(γ i -γ) p ]) 1/p ≤ γ m • E E ′ (tr(Q -i ) -x ⊤ i UQ -i U ⊤ x i ) p 1/p + (E E ′ [tr(Q -i ) -E E ′ tr(Q -i )] p ) 1/p .Using our new Restricted Bai-Silverstein inequality from Lemma 3 (restated as Lemma 16 and proven in Appendix D), we haveE E ′ (tr(Q -i ) -x ⊤ i UQ -i U ⊤ x i ) p 1/p < c • p 3 √ d • 1 + dp log(d/δ) s .We now consider(E E ′ [tr(Q -i ) -E E ′ tr(Q -i )] p )1/p . In Lemma 4 (restated below as Lemma 13), we 5 (m/δ) with probability at least 1-δ. Conditioned on this high-probability event we have,show that |tr(Q -i )-E E ′ tr(Q -i )| ≤ c ′ γ √ m •d log 4.(E E ′ [tr(Q -i ) -E E ′ tr(Q -i )] p ) 1/p ≤ c ′ γ √ m • d log 4.5 (m/δ)for an absolute constant c ′ > 0. Therefore we get with probability at least 1 -δ,(E E ′ [(γ i -γ) p ]) 1/p ≤ γ m c • p 3 √ d • 1 + dp log(d/δ) s + c ′ γ √ m • d log 4.5 (m/δ) . As m > d, we get (E E ′ [(γ i -γ) p ]) 1/p = O √ d log 4.5 (n/δ) m • 1 + d s. Also using the analysis in [18] for upper bounding |γ -γ|, we get a matching upper bound on |γ -γ| as follows:|γ -γ| = O √ d log 4.5 (n/δ) m • 1 + d s .Substituting these bounds in (7) and then in (6) we get,EQUATIONCombining the upper bounds for Z 0 , Z 1 and Z 2 using relations (4,5,8) we conclude our proof.We now provide the proof of Lemma 4, which we restate in the following Lemma.Lemma 13. For given δ > 0 and matrix Q -i we have with probability 1 -δ:|tr(Q -i ) -E E ′ tr(Q -i )| ≤ c ′ γ √ m • d log 4.5 (m/δ)for an absolute constant c ′ > 0.Proof. Writing tr(Q -i ) -E E ′ tr(Q -i ) as a finite sum, we have,tr(Q -i ) -E E ′ tr(Q -i ) = m j=1 E E ′ ,j tr(Q -i ) -E E ′ ,j-1 tr(Q -i ).DenotingX j = E E ′ ,j tr(Q -i ) with X 0 = E E ′ tr(Q -i ),we have the following formulationtr(Q -i ) -E E ′ tr(Q -i ) = m j=1 X j -X j-1 with E E ′ ,j-1 [X j ] = X j-1 .The random sequence X j forms a martingale and tr(Q -i ) - E E ′ tr(Q -i ) = X m -X 0 .We find an upper bound on |X j -X j-1 |. To achieve that we note,X j -X j-1 = E E ′ ,j tr(Q -i ) -E E ′ ,j-1 tr(Q -i ) = -(E E ′ ,j -E E ′ ,j-1 )(tr(Q -ij -Q -i ) -tr(Q -ij )).Therefore withψ j = (E E ′ ,j -E E ′ ,j-1 )tr(Q -ij -Q -i ) and χ j = -(E E ′ ,j -E E ′ ,j-1 )tr(Q -ij ),we have,|X j -X j-1 | ≤ |ψ j + χ j | ≤ |ψ j | + |χ j |.From [18] , we have |χ j | ≤ 1 m . We now prove an upper bound on ψ j .0 ≤ tr(Q -ij ) -tr(Q -i ) = tr γ m Q -ij U ⊤ x j x ⊤ j UQ -ij U ⊤ U 1 + γ m x ⊤ j UQ -ij U ⊤ x j = γ m x ⊤ j (UQ -ij U ⊤ ) 2 x j 1 + γ m x ⊤ j UQ -ij U ⊤ x j ≤ γ m x ⊤ j UQ 2 -ij U ⊤ x j .Now look at the term x ⊤ j UQ 2 -ij U ⊤ x j . For any a > 0 and any k > 0, by Markov's inequality we have,Pr x ⊤ j UQ 2 -ij U ⊤ x j ≥ a ≤ E[|x ⊤ j UQ 2 -ij U ⊤ x j | k ] a k ≤ 2 k-1 • E[|x ⊤ j UQ 2 -ij U ⊤ x j -tr(Q 2 -ij )| k a k + 2 k-1 • E[(tr(Q 2 -ij )) k ] a k.Let E 1 be an event independent of both x i and x j and have probability at least 1 -δ ′ . Therefore we haveEQUATIONWe now upper bound both terms on the right-hand side separately. Considering the termE E1 [|x ⊤ j UQ 2 -ij U ⊤ x j -tr(Q 2 -ij )| k and using Lemma 3 we get,EQUATION)Now considering the second term in (9), i.e., E E1 [(trEQUATIONSubstituting ( 10) and ( 11) in (9) ,Pr x ⊤ j UQ 2 -ij U ⊤ x j ≥ a ≤ 2 k • c k • k 3k • d 2 k log(d/δ)) s + d k/2 a k + 2 k • 6 k • d k a k Pr x ⊤ j UQ 2 -ij U ⊤ x j ≥ a ≤ 2 k • c k • k 3k • d k • (k log(d/δ)) k/2 a kfor some potentially different constant c. Consider k = log(m/δ) log( 2)and a = 4 • c • k 3 • d • k log(d/δ)and we have with probability at least 1 -δ/m,x ⊤ j UQ 2 -ij U ⊤ x j ≤ 4 • c • k 3 • d • log(kd/δ). This implies that for an absolute constant c ′ we have,|tr(Q -ij ) -tr(Q -i )| ≤ c ′ • γ m • d • log 3.5 (m/δ).Therefore we now have an upper bound for|ψ j | |ψ j | = |E E ′ ,j -E E ′ ,j-1 (tr(Q -ij -Q -i ))| ≤ 2c ′ • γ m • d • log 3.5 (m/δ).This means for all j, we have, 3.5 (m/δ). We get with probability at least 1 -δ and for potentially different absolute constant c ′ > 0:|X j -X j-1 | ≤ 4c ′ • γ m • d • log 3.5 (m/δ) with probability at least 1 -δ. Consider c j = 4c ′ • γ m • d • log 3.5 (m/δ). Then m j=1 c 2 j = γ 2 m • 16c ′2 • d 2 log 7 (m/δ). Applying Azuma's inequality (Lemma 8) with λ = γ √ m • 4c ′ • d • log|X m -X 0 | ≤ c ′ γ √ m • d log 4.5 (m/δ).This concludes our proof.We consider the following Distributed Newton Sketch method:x t+1 = x t -p t , p t = 1 q q i=1 H-1 i g t ,where g t = ∇f (x t ) and Hi = m m-d A ⊤ t S ⊤ i S i A t is a sketched estimate of the Hessian H t = ∇ 2 f (x t ), when expressed as the matrix product H t = A ⊤ t A t for the appropriately chosen n × d matrix A t (as can be done for any generalized linear model). This update can be viewed as an approximate Newton step with the Hessian inverse estimate Ĥ-1t = 1 q q i=1H-1 i . As long as each H-1t is an ( √ ϵ, δ/2q)-unbiased estimator of H -1 t , then using Lemma 34 from [18] we get that the averaged Hessian inverse with q = Õ(1/ϵ) satisfies:1 1 + ε H -1 t ⪯ Ĥ-1 t ⪯ (1 + ε)H -1 t , ε = √ ϵ + Õ(1/ √ q) = O( √ ϵ).Using standard approximate Newton analysis (e.g., see Lemmas 1 and 3, and the related discussion in [15] ), this implies that when x t is in a sufficiently small neighborhood around x * (determined solely by the strong convexity and Lipschitz constants of the Hessian of f ), we have:f (x t+1 ) -f (x * ) f (x t ) -f (x * ) = O(ε 2 ) = O(ϵ).Adjusting the constants appropriately, relying on Theorem 3 for constructing Hi , and following the complexity analysis from Theorem 2, we recover the claim.In this section, we aim to prove Theorem 4. Let x * = argmin x ∥Ux -b∥ 2 where U ∈ R n×d is the data matrix containing n data points and b ∈ R n is a vector containing labels corresponding to n data points. We adopt the same notations as used in the proof of Theorem 5. Let x = argmin x ∥S m Ux -S m b∥ 2 . Furthermore for any x ∈ R d we can find the loss at x as L(x) = ∥Ux -b∥ 2 . Additionally, we use r to denote the residual b -Ux * . We aim to provide an upper bound on the bias introduced due to this sketch and solve paradigm, i.e. L(E(x)) -L(x * ). Similar to Theorem 4 we condition on the high probability event E and consider L(E E (x)) -L(x * ). By Pythagorean theorem, we haveL(E E [x]) -L(x * ) = ∥U(E E [x]) -Ux * ∥ 2 . Also, L(E E [x]) -L(x * ) = ∥U(E E [x]) -Ux * ∥ 2 = ∥E E [x] -x * ∥ 2 = E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m b] -U ⊤ b 2 = E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m ](b -UU ⊤ b) 2 . Note that b -UU ⊤ b = b -Ux * = r. We get, L(E E [x]) -L(x * ) = E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m ]r 2 . Consider Q = (γU ⊤ S ⊤ m S m U) -1 , Q -i = (γU ⊤ S ⊤ -i S -i U) -1 and γ = m m-d , E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m r] = E E [γ(γU ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m r] = γE E [QU ⊤ S ⊤ m S m r] = γE E [QU ⊤ x i x ⊤i r] where we used linearity of expectation in the last line combined withS ⊤ m S m = 1 m m i=1 x i x ⊤ i . Using Sherman-Morrison formula (Lemma 5) we have QU ⊤ x i = Q-iU ⊤ xi 1+ γ m x ⊤ i UQ-iU ⊤ xi . Denote γ i = 1 + γ m x ⊤ i UQ -i U ⊤x i and substitute we get,E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m r] = E E γ γ i Q -i U ⊤ x i x ⊤ i r E E [(U ⊤ S ⊤ m S m U) -1 U ⊤ S ⊤ m S m r] = E E [Q -i U ⊤ x i x ⊤ i r] + E E γ γ i -1 Q -i U ⊤ x i x ⊤ i r .So we get the following decomposition:EQUATION)Note that a similar decomposition was considered in the proof of Theorem 5 (see Appendix B) with slightly different Z 0 and Z 2 . We first bound ∥Z 0 r∥ 2 in the following argument. Without loss of generality, we assume that events E 1 and E 2 are independent of x i andEQUATIONNow note that in the first term Q -i is independent of x i and x i is also independent of the event E ′ . Using this with the fact that E[x i x ⊤ i ] = I we get E E ′ [Q -i U ⊤ x i x ⊤ i r] = E E ′ [Q -i U ⊤ r] = 0, since U ⊤ r = 0. Therefore, ∥Z 0 r∥ 2 ≤ 2 ∥E E ′ [Q -i U ⊤ x i x ⊤ i r • 1 ¬E3 ]∥ 2 ≤ 72∥E E ′ [U ⊤ x i x ⊤ i r • 1 ¬E3 ]∥ 2 .The last inequality holds because conditioned on E ′ , we know that ∥Q -i ∥ ≤ 6. Using Cauchy-Schwarz inequality we have,∥Z 0 r∥ 2 ≤ 72 E E ′ [x ⊤ i UU ⊤ x i • 1 ¬E3 ] • E E ′ [r ⊤ x i x ⊤ i r] 2 .Note thatE E ′ [r ⊤ x i x ⊤ i r] = ∥r∥ 2 . Also, we can bound E E ′ [x ⊤ i UU ⊤ x i • 1 ¬E3] as following:E E ′ [x ⊤ i UU ⊤ x i • 1 ¬E3 ] = ∞ 0 Pr(x ⊤ i UU ⊤ x i • 1 ¬E3 > x)dx = y 0 Pr(x ⊤ i UU ⊤ x i • 1 ¬E3 > x)dx + ∞ y Pr(x ⊤ i UU ⊤ x i • 1 ¬E3 > x)dx ≤ yδ ′ + ∞ y Pr(x ⊤ i UU ⊤ x i > x)dx. Using Chebyshev's inequality, Pr(x ⊤ i UU ⊤ x i ≥ x) ≤ Var(x ⊤ i UU ⊤ xi) x 2. By Restricted Bai-Silverstein, Lemma 3 with p = 2, we have Var(x ⊤ i UU ⊤ x i ) ≤ c • d + d 2 log(d/δ) sfor some absolute constant c.Let y = m 2 and δ ′ < 1 m 4 we get, E[x ⊤ i UU ⊤ x i • 1 ¬E3 ] = O 1 m 2 + d m 2 + d 2 log(d/δ) sm 2 = O d m 2 + d 2 log(d/δ) sm 2 .This finishes upper bounding ∥Z 0 r∥ 2 as:EQUATION)Now we proceed with ∥Z 2 r∥ 2 ,∥Z 2 r∥ = E E γ γ i -1 Q -i U ⊤ x i x ⊤ i r .Applying Hölder's inequality with p = O(log(n)) and q = 1 + ∆ where ∆ = 1 O(log(n)) , we get,∥Z 2 r∥ ≤ E E γ γ i -1 p 1/p • sup ∥v∥=1 E E [v ⊤ Q -i U ⊤ x i x ⊤ i r] q 1/q ≤ E E γ γ i -1 p 1/p • E E ∥Q -i U ⊤ x i ∥ 2q 1/2q • E E ∥x ⊤ i r∥ 2q 1/2q ≤ E E γ γ i -1 p 1/p • 2 • E E ′ [∥Q -i U ⊤ x i ∥ 2 • ∥Q -i U ⊤ x i ∥ 2∆ ] 1/2q • 2 • E E ′ [∥x ⊤ i r∥ 2 • ∥x ⊤ i r∥ 2∆ ] 1/2q.Since ∥x i ∥ = poly(n), we have ∥x i ∥ 2∆ = O(1) and therefore we have∥Q -i U ⊤ x i ∥ 2∆ = O(1) ∥Q -i U ⊤ ∥ 2∆ . Also E E ′ ∥Q -i U ⊤ x i ∥ 2 = E E ′ ∥Q -i ∥ 2 . Using that conditioned on E ′ we have ∥Q -i ∥ ≤ 6, we get E E ∥Q -i U ⊤ x i ∥ 2q 1/2q = O(1). Similarly using ∥x i ∥ 2∆ = O(1) we get E E ∥x ⊤ i r∥ 2q 1/2q= O(1) ∥r∥. This gives us:∥Z 2 r∥ ≤ O(1) • E E γ γ i -1 p 1/p• ∥r∥ .Now using (8) from the proof of Theorem 5, we get,E E γ γ i -1 p 1/p = O √ d log 4.5 (n/δ) m • 1 + d s .Finally the bound for ∥Z 2 r∥ 2 follows as:EQUATION)Combining ( 13) and ( 14) we conclude our proof.In this section, we prove Lemma 3. We need the following two auxiliary lemmas to derive the main theorem of this section. The first lemma uses Matrix-Chernoff concentration inequality to upper bound the spectral norm of U ⊤ ξ U ξ where ξ is a (s, β 1 , β 2 )-approximate leverage score sparsifier. Lemma 14 (Spectral norm bound with leverage score sparsifier). Let U ∈ R n×d has orthonormal columns. Let ξ be a (s, β 1 , β 2 )-approximate leverage score sparsifier for U, and denote U ξ = diag(ξ)U. Then for any δ > 0 we have,Pr U ⊤ ξ U ξ ≥ 1 + 3d log(d/δ) s ≤ δ if s < d, Pr U ⊤ ξ U ξ ≥ (1 + 3 log(d/δ)) ≤ δ if s ≥ d.Proof. Writing U ⊤ ξ U ξ as a sum of matrices we have,U ⊤ ξ U ξ = n i=1 ξ 2 i u i u ⊤ i = n i=1 b i p i u i u ⊤ i = n i=1 Z iwhere p i = min{1, sβ1 li d } and Z i = bi pi u i u ⊤ i . Note that Z ′ i s are independent random variables andE[Z i ] = u i u ⊤ i . Also n i=1 EZ i = U ⊤ U = I d . If p i = 1 then Z i = u i u ⊤ i and therefore ∥Z i ∥ = ∥u i ∥ 2 ≤ 1. If p i < 1, we have ∥Z i ∥ ≤ 1 pi ∥u i ∥ 2 = d sβ1 li • l i . As l i ≤ β 1 li , we get ∥Z i ∥ ≤ d s . Therefore ∥Z i ∥ ≤ max{1, d s } for all i. Denote R = max{1, d s }.We use Matrix Chernoff (Lemma 7) to upper bound the largest eigenvalue of U ⊤ ξ U ξ . For any ϵ > 0, we have,Pr λ max n i=1 Z i ≥ (1 + ϵ) ≤ d • exp - ϵ 2 (2 + ϵ)R .With R = max{1, d s } and depending on the case whether s ≤ d or s > d we getPr λ max n i=1 Z i ≥ 1 + 3d log(d/δ) s ≤ δ if s ≤ d, Pr λ max n i=1 Z i ≥ (1 + 3 log(d/δ)) ≤ δ if s > d. Let A ξ denote the event U ⊤ ξ U ξ ≤ 1 + 3d log(d/δ) s, holding with probability at least 1 -δ, for small δ > 0. In the next result we upper bound the higher moments of the trace of U ξ CU ⊤ ξ for any matrix C ⪯ O(1) • I. We first prove the upper bound in the case when the high probability event A ξ does not occur.Lemma 15 (Trace moment bound over small probability event). Let k ∈ N be fixed. Let U ∈ R n×d have orthonormal columns. Let ξ be a (s, β 1 , β 2 )-approximate leverage score sparsifier for U. Let U ξ = diag(ξ)U. Also let Pr(A ξ ) ≥ 1 -1 (12d) 4k and event E ′ be independent of the sparsifier ξ.Then we have,E E ′ (tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ )) k |¬A ξ ≤ (4k) 4k for any fixed matrix C such that 0 ⪯ C ⪯ 6I.Proof.E E ′ (tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ )) k |¬A ξ = ∞ 0 Pr((tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ )) k • 1 ¬A ξ ≥ x|E ′ )dx = ∞ 0 Pr((tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ )) k • 1 ¬A ξ ≥ x)dx.The last equality holds because E ′ is independent of ξ. Consider some fixed y > 0.EQUATION)The last inequality holds because by Chebyshev's inequality Pr((tr(U ξ CU ⊤ ξ )) 2k ≥ x) ≤ E[(tr(U ξ CU ⊤ ξ )) 4k ] x 2. Also note that,(tr(U ξ CU ⊤ ξ )) 4k = n i=1 ξ 2 i u ⊤ i Cu i 4k = n i=1 b i p i u ⊤ i Cu i 4k .For 1 ≤ i ≤ n, let R i be random variables denoting bi pi u ⊤ i Cu i . Then note that R i are independent random variables with E[R i ] = u ⊤ i Cu i . Also R i are non-negative random variables with finite (4k) th moment. Using Rosenthal's inequality (Lemma 9) we get,E   n i=1 b i p i u ⊤ i Cu i 4k   ≤ 2 4k • (4k) 4k •   n i=1 E[R 4k i ] + n i=1 E[R i ] 4k   . Now n i=1 E[R i ] = tr(UCU ⊤ ) and E[R 4k i ] can be found as follows: if p i = 1 then R i = u ⊤ i Cu i and therefore R 4k i = (u ⊤ i Cu i ) 4k ≤ ∥u i ∥ 2(4k-1) u ⊤ i C 4k u i ≤ u ⊤ i C 4k u i , if p i < 1, we have, E[R 4k i ] = p i • 1 p 4k i (u ⊤ i Cu i ) 4k = d 4k-1 s 4k-1 • 1 (β 1 li ) 4k-1 (u ⊤ i Cu i ) 4k ≤ d 4k-1 s 4k-1 • 1 (β 1 li ) 4k-1 • ∥u i ∥ 2(4k-1) u ⊤ i C 4k u i = d 4k-1 s 4k-1 • 1 (β 1 li ) 4k-1 • l 4k-1 i u ⊤ i C 4k u i . Now using l i ≤ β 1 li we get, E[R 4k i ] ≤ d 4k-1 s 4k-1 • u ⊤ i C 4k u i .Therefore,E   n i=1 b i p i u ⊤ i Cu i 4k   ≤ 2 4k • (4k) 4k • max 1, d 4k-1 s 4k-1 • tr(C 4k ) + (tr(C)) 4k .Using the above inequality along with using C ⪯ 6I we get,E   s t=1 1 sp it u ⊤ it Cu it 4k   ≤ 2 4k • (4k) 4k • 6 4k • d 4k = (12) 4k • (4k) 4k • d 4k .Substituting the above bound in (15) , it follows that:E E ′ (tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ )) k |¬A ξ ≤ y • δ + (12) 4k • (4k) 4k • d 4k • 1 y .For y > (12d) 4k and δ < 1 y , we get the desired result.We are now ready to prove the main result of this section. The following result, which is central to our analysis, upper bounds the high moments of a deviation of a quadratic form from its mean. Lemma 16 (Restricted Bai-Silverstein for (s, β 1 , β 2 )-LESS embedding). Let p ∈ N be fixed and U ∈ R n×d have orthonormal columns. Let x i = diag(ξ)y i where y i ∈ R n has independent ±1 entries and ξ is a (s, β 1 , β 2 )-approximate leverage score sparsifier for U. Let U ξ = diag(ξ)U.Then for any matrix 0 ⪯ C ⪯ 6I and any δ > 0 we have,E [tr(C) -x ⊤ i UCU ⊤ x i ] p 1/p < c • p 3 √ d • 1 + dp log(d/δ) sfor an absolute constant c > 0.Proof. Let x i = diag(ξ)y i where y i is vector of Rademacher ±1 entries. Denote U ξ = diag(ξ)U.EQUATIONFirst, consider T 1 , substitute ξ 2 i = bi pi , and assume exponent p to be even.E tr(C) -tr(U ξ CU ⊤ ξ ) p = E tr(UCU ⊤ ) -tr(U ξ CU ⊤ ξ ) p = E n i=1 (ξ 2 i -1)u ⊤ i Cu i p = E n i=1 ( b i p i -1)u ⊤ i Cu i p . For 1 ≤ i ≤ n consider random variables R i where R i = bi pi u ⊤ i Cu i . Furthermore let Y i = R i -u ⊤ i Cu i . Then n i=1 Y i = n i=1 ( bi pi -1)u ⊤ i Cu i and E[Y i ] = 0.Y i are independent mean zero random variables with finite p th moments and therefore we can use Rosenthal's inequality (for symmetric random variables, Lemma 9) to get,EQUATIONwhere A(p) is a constant depending on p. We bound T 1 1 and T 2 1 separately, starting with T 1 1 as,T 1 1 = n i=1 E[Y i ] p . Recall that Y i = R i -u ⊤ i Cu i . R i is always non-negative because C is a positive semi-definite matrix. Therefore we have, E(Y i ) p ≤ E(R i ) p . We find the p th moment of R i , E(R i ) p = (u ⊤ i Cu i ) p if p i = 1. If p i < 1 then, E(R i ) p = p p-1 i (u ⊤ i Cu i ) p = d p-1 s p-1 • 1 (β 1 li ) p-1 (u ⊤ i Cu i ) p ≤ d p-1 s p-1 • 1 (β 1 li ) p-1 • ∥u i ∥ 2(p-1) u ⊤ i C p u i ≤ d p-1 s p-1 u ⊤ i C p u iwhere in the last inequality we use∥u i ∥ 2(p-1) = l p-1 i ≤ (β 1 li ) p-1 .Summing over i from 1 to n we get an upper bound for T 1 1 as the following:EQUATIONNow we upper boundT 2 1 . Note that E(Y i ) 2 = E[R i ] 2 = (u ⊤ i Cu i ) 2 if p i = 1 and if p i < 1 we have E[R i ] 2 = 1 pi • (u ⊤ i Cu i ) 2 ≤ d s • 1 β1 li ∥u i ∥ 2 u ⊤ i C 2 u i ≤ d s u ⊤ i C 2 u i .Summing over i from 1 to n we get an upper bound for T 2 1 as,EQUATIONSubstituting ( 18) and ( 19) in (17) we have,EQUATIONNow we aim to upper bound the term T 2 in (16) i.e., E tr(U ξ CU ⊤ ξ ) -y ⊤ i U ξ CU ⊤ ξ y i p. First, we condition over ξ and take expectation over y i . This requires using standard Bai-Silverstein inequality (Lemma 10), we get,E tr(U ξ CU ⊤ ξ ) -y ⊤ i U ξ CU ⊤ ξ y i p ≤ B(p) • E ν 4 tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ ) p/2 + ν 2p tr(U ξ CU ⊤ ξ U ξ CU ⊤ ξ ) p/2where B(p) is a constant depending on p. Since y i consists of ±1 entries we have ν 4 , ν 2p ≤ 1. Also using tr(AB) ≤ tr(A)tr(B) and considering the high probability event A ξ capturing U ⊤ ξ U ξ ≤ 1 + max 3d log(d/δ) s , 3 log(d/δ) . We get the following,EQUATIONThe first term in the last inequality follows from the Matrix-Chernoff (Lemma 7) and the second term follows from Lemma 15 by considering k = p/2 (assuming that δ is small enough so that Lemma 15 is satisfied). We now upper-bound E tr(U ξ C 2 U ⊤ ξ ) p/2 using Rosenthal's inequality for uncentered (non-symmetric) random variables,tr(U ξ C 2 U ⊤ ξ ) = n i=1 ξ 2 i u ⊤ i C 2 u i = n i=1 b i p i u ⊤ i C 2 u i . For 1 ≤ i ≤ n consider independent random variables R ′ i = bi pi u ⊤ i C 2 u i . We have, n i=1 R ′ i = tr(U ξ C 2 U ⊤ ξ ).Here R ′ i are positive random variables with finite (p/2) th moment. Using Rosenthal's inequality we get,EQUATIONIt is straightforward to upper bound E(R ′ i ) p/2 as,E(R ′ i ) p/2 = (u ⊤ i C 2 u i ) p/2 if p i = 1, ≤ d p 2 -1 s p 2 -1 u ⊤ i C p u i if p i < 1.Summing over i from 1 to n we get upper bound for T 1 2 as,EQUATIONNow we consider T 2 2 . It is simply given as,EQUATIONCombining ( 23) and ( 24) and substituting in (22) we get,EQUATION)Substituting (25) in (21) and let w = max(1, d/s),EQUATIONCombining the bounds for T 1 (20) and T 2 (26) substituting in (16) , and noting that tr(UC k U ⊤ ) = tr(C k ) for any k, we get,E (tr(C) -x ⊤ i UCU ⊤ x i ) p ≤ 2 p-1 A(p) • w p-1 • tr(C p ) + w • tr(C 2 ) p/2 +2 p B(p)C(p/2)• (1 + 3w log(d/δ)) p/2 • w p 2 -1 tr(C p ) + tr(C 2 ) p/2 +2 p B(p)(2p) 2p .Now we specify the various constants depending on p. We haveA(p) ≤ (2p) p , B(p) ≤ (2p) p , C(p/2) ≤ p p/2 . Also we use tr(C k ) ≤ 6 k • d since C ⪯ 6I.This implies for an absolute constant c we have,E [tr(C) -x ⊤ i UCU ⊤ x i ] p 1/p ≤ c • p 3 √ d • 1 + dp log(d/δ) swhere δ > 0 is now arbitrary.In this section, we extend our main result to settings where the dataset is partitioned across multiple machines. We illustrate this in the distributed setting considered in [43] . Let A ∈ R N ×d and b ∈ R N with rank(A) = d. In practice, we randomly shuffle the rows of A and b, and partition the data uniformly across q machines, with every machine getting chunks of size n denoted by (A i , b i ). In this distributed setup, A i can be considered a uniformly subsampled sketch of A with sketch size n. However, similar to the analysis in [43] , in the following analysis we assume that each machine constructs its sketch (A i , b i ) by uniformly sampling rows from A and b with replacement and independently from other machines. This setting is not ideally a partition but we believe similar results can also be shown for the setting where we partition the dataset after reshuffling the rows. Let µ = N max n i=1 ℓ i (A) denote the matrix coherence of A. Let S ∈ R n×N be a uniform sampling matrix with size n, where we will assume that S also scales the samples by N/n, which implies that E[S ⊤ S] = I. Then the following holds: Lemma 17 (Subspace embedding using uniform subsampling; based on Theorem 12 in [43] ). Let S ∈ R n×N be a uniform sampling matrix of size n where n ≥ O(µ log(d/δ)). Then with probability 1 -δ we have,1 1 + η • A ⊤ A ⪯ A ⊤ S ⊤ SA ⪯ (1 + η) • A ⊤ A for η = O( µ log(d/δ)/n).Let x * = A † b denote the minimizer of least squares loss ∥Ax -b∥2 and x * i = A † i b i denote the solution to the sketched least squares problem at i th machine. Here A i = S i A and b i = S i b, where S i denotes the uniform subsampling matrix for i th machine. We can use the model averaging result from [43] to prove the following result. Lemma 18 (Model averaging result; adapted from Theorem 20 in [43] ). Let there be q machines and each machine constructs its sketch (A i , b i ) independently by uniformly sampling n rows from A and b, where n = O(µ log(dq/δ)). Then with probability 0.99,EQUATION)where x * = 1 q q i=1 x * i and c > 0 is an absolute constant 2 .Since computing x * i at any machine could be potentially expensive due to large n, we instead find xi at i th machine using the sketching procedure from Theorem 2. Note that this requires applying our debiasing techniques to the least squares problem min ∥A i xb i ∥ 2 at every machine. To that end, let x = 1 q q i=1 xi and consider E P ∥Ax -b∥ 2 , where E P means conditional expectation given the sketches (A i , b i ) at every machine. The following theorem states that if the bias from partitioning is much smaller than the sketching error in each machine, then our approach can be successfully used to reduce the bias of the sketching estimate so that averaging will work as desired. Theorem 6. Let there be q machines and each machine constructs its sketch (A i , b i ) independently by uniformly sampling n rows from A and b, where n ≥ O(µ log(dq/δ)). Also, let i th machine construct xi via locally sketching A i and b i such that xi satisfies Theorem 1 (with A i and b i ) at i th machine. Then, the estimator x = 1 q q i=1 xi averaged across the machines satisfiesE P ∥Ax -b∥ 2 ≤ 1 + c ′ ϵ + 1 q + µ log(d/δ) qn + µ 2 log 2 (d/δ) n 2 ∥Ax * -b∥ 2 .with probability at least 0.9 and an absolute constant c ′ > 0.Proof. By the Pythagorean theorem we have:EQUATION)We proceed to upper bound E P ∥A(x -x * )∥ 2 as,E P ∥A(x -x * )∥ 2 ≤ 2 E P ∥A(x -x * )∥ 2 + ∥A(x * -x * )∥ 2 .We use Lemma 18 to upper bound the second term and get,E P ∥A(x -x * )∥ 2 ≤ 2 E P ∥A(x -x * )∥ 2 + c µ log(d/δ 1 ) qn + µ 2 log 2 (d/δ 1 ) n 2 ∥Ax * -b∥ 2 .(29) On the other hand,EQUATION)We upper bound both terms in (30) separately. Proceeding with the first term we have,EQUATIONEQUATIONEQUATIONEQUATION)The equality (31) holds as all x′ i s are independent. The inequality (32) holds due to the subspace embedding property (Lemma 17). The inequality (33) holds due to the variance bound in Theorem 1. The last inequality holds because x * i minimizes the least squares loss ∥A i xb i ∥ 2 at i th machine. We now upper bound the second term in (30) asEQUATION≤ ϵ(1 + η) q • q i=1 ∥A i x * i -b i ∥ (36) ≤ ϵ(1 + η) q • q i=1 ∥A i x * -b i ∥The inequality (35) holds due to the subspace embedding property (Lemma 17) and inequality (36) holds due to the bias bound in Theorem 1. Therefore we get,EQUATION)Combining (28, 29, 34, 37) and assuming η < 1 we get,EQUATIONAs E 1 q q i=1 ∥A i x * -b i ∥ 2 = ∥Ax * -b∥ 2 , we just use Markov's inequality to upperbound the second term in (38) by 160 q ∥Ax * -b∥ 2 with probability 0.95. For the last term, we haveq i=1 ∥A i x * -b i ∥ 2 = q • 1 q q i=1 ∥A i x * -b i ∥ 2 + (q 2 -q) • 1 q 2 -q i,j,i̸ =j ∥A i x * -b i ∥ • ∥A j x * -b j ∥ Again 1 q q i=1 ∥A i x * -b i ∥ 2 ≤ 20∥Ax * -b∥ 2 withprobability 0.95. The last term is an average of q 2 -q random variables and we can provide uniform bound for the expectation of all of them using the independence among different machines. We haveE[∥A i x * -b i ∥ • ∥A j x * -b j ∥] ≤ E[∥A i x * -b i ∥ 2 ] • E[∥A j x * -b j ∥ 2 ] = ∥Ax * -b∥ 2 .Again using Markov's inequality we can upper bound the average of these q 2 -q random variables with 20 times their expectation with probability 0.95, finishing the proof of the theorem.Here, we provide a small set of numerical experiments to empirically examine the relative error of distributed averaging estimates from individual machines to return an estimator x = 1 q q i=1 xi . First, we show that when the sketching-based estimates xi do have a non-negligible bias, so that the distributed averaging estimator x remains inconsistent in the number of machines -even with an unlimited number of machines, as long as the space on each machine is limited, the averaged estimator's performance will be limited by the bias of the individual estimates. Second, we show that one can use sketching to compress a data subsample at no extra computational cost, without increasing its bias, which we refer to as the free lunch in distributed averaging via sketching.We examine three benchmark regression datasets, Abalone (4177 rows, 8 features), Boston (506 rows, 13 features), and YearPredictionMSD (truncated to the first 2500 rows, with 90 features), from the libsvm repository from [8] . All the experiments were run on an i9-13900k processor with 128GB of RAM and an NVIDIA RTX 3090 GPU. We repeat the experiment 100 times, with the shaded region representing the standard error. We visualize the relative error of the averaged sketchand-solve estimator L(x)-L(x * )L(x * ), against the number of machines q used to generate the estimate x = 1 q q i=1 xi . Each estimate xi was constructed with the same sparsification strategy used by LESS, except that instead of sparsifying the sketch with leverage scores, we instead sparsify them with uniform probabilities. Following [16] , we call the resulting method LESSUniform. Within each dataset, we perform four simulations, each with different sketch sizes and different numbers of nonzero entries per row. We vary these so that the product (sketch size × nnz per row) stays the same, so as to ensure that the total cost of sketching is fixed in each plot.As expected, decreasing the sketch size while increasing the number of nonzeros per row (effectively increasing the amount of "compression" occurring here by sparse sketching) increases the error in all three datasets. However, remarkably, increasing the amount of "compression" does not seem to increase the amount of bias. We can therefore conclude that sparse sketches preserve nearunbiasedness, while enabling us to reduce the sketch size from subsampling without incurring any additional computational cost. The increase in error can be mitigated in a distributed setting by increasing the number of estimates/machines.To further illustrate the phenomenon that suitable sketched least squares estimators enjoy small bias, we provide a further set of numerical experiments below in Figure 4 for the Boston and Abalone Figure 3 : Comparison of the relative error of the distributed averaging estimator of sketch-andsolve least squares estimates where the sketches are constructed with sparse sketching matrices with uniform probabilities (LESSUniform) on libsvm dataset Boston (see Figure 2 for results on YearPredictionMSD and Abalone). For each dataset, the computational cost of sketching is the same in all four parameter settings. Remarkably, sketching to a smaller size appears to preserve near-unbiasedness without incurring any additional computational cost. datasets. Recall that in Figures 2 and 3 , we show that increasing the amount of "compression" of the sample, by reducing the sketch size when using sparse sketching matrices with uniform probabilities, does not increase the bias of the sketched estimator. Figure 4 further demonstrates that Gaussian and Subgaussian sketches, which are computationally much more expensive than our proposed sparse sketches, exhibit virtually no bias in sketched least squares (as expected). A similar conclusion can be made for the proper LESS method, which uses leverage score estimates, and is therefore more expensive than the LESSUniform method we use in other experiments, but still potentially much cheaper than Subgaussian sketches.However, the unbiasedness does not hold for all sparse sketches, as we demonstrate below in Figure 4 . At a high level, we show that sketches constructed using leverage score subsampling and the SRHT can exhibit a non-negligible level of least squares bias.Gaussian and Subgaussian sketches, as we see in the figure, enjoy near-unbiasedness while not introducing additional bias as the sketch size decreases. This does not hold for leverage score subsampling. When we decrease the number of subsamples within leverage score subsampling, the bias introduced by subsampling increases as suggested by the lower bound in Theorem 10 of [18] . Intuitively, this happens as the number of subsamples is reduced without increasing the amount of "compression" as one would with LESS or LESSUniform.We also show that the subsampled randomized Hadamard transform (SRHT) also introduces increased bias as the sketch size decreases. This complements the lower bound established for the failure of SRHT and other data-oblivious sparse embeddings like CountSketch to satisfy the restricted Bai-Silverstein property that is a structural condition for near-unbiasedness sketches [18] .For completeness, we also demonstrate the desirable performance of LESS proper. Figure 4 demonstrates that the desirable phenomenon that LESSUniform enjoys also extends to LESS proper. In fact, we observe that LESS enjoys similar desirable performance as the Gaussian and Subgaussian sketches and does not increase (in fact, it has minimal) least squares bias, while enjoying the computational speedups that sparse sketches also enjoy, suggesting a best-of-both-worlds property. 4 : Distributed averaging experiment for the Boston and Abalone datasets using five different sketching techniques: Leverage Score Subsampling, Subsampled Randomized Hadamard Transform (SRHT), Gaussian sketches, Subgaussian sketches, and LESS. Both Gaussian and Subgaussian sketches exhibit no observable least squares bias, whereas Leverage Score Subsampling and SRHT exhibit a small but measurable level of bias. LESS enjoys the similar high performance as the Gaussian and Subgaussian sketches and does not increase (in fact, it has minimal) least squares bias, but as a sparse sketch also allows for the significant improvements in runtime that sparse sketches enjoy, demonstrating that it enjoys the best-of-both-worlds. Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.One can think of LESS and LESSUniform as generalizations of leverage score subsampling and uniform subsampling respectively, where we mix Õ(1/ϵ) subsamples into a single row of the sketch. This suggests that sparse sketching, as an extension of subsampling, yields desirable bias reduction without significantly increasing computational costs.The result Theorem 20[43] also requires that ∥S ⊤ i Si∥ 2 ≤ N n for all i, which holds trivially if Si is a sampling matrix without replacement. However, in our setting, we are implicitly assuming that with high probability ∥S ⊤ i Si∥ 2 = O( N n ) and absorb the extra resulting constant in c. 


Decision: This paper proves an interesting headline theoretical result on faster small space distribution regression, requires non-trivial analysis, and also provides nice initial experiments on the method. There are some limitations — for example, that the distributed regression setup the authors consider is fairly niche. Nevertheless, we recommend accepting this paper to NeurIPS.

Reviews: review: paper_summary: paper_summary: This paper presents new techniques for distributed least squares regression using matrix sketching. The key contributions are:

1. A sparse sketching method that produces a nearly-unbiased least squares estimator in two passes over the data, using optimal space and current matrix multiplication time.
2. Improved communication-efficient distributed averaging algorithms for least squares and related tasks.
3. A novel bias analysis for sketched least squares, characterizing its dependence on sketch sparsity. This includes new higher-moment restricted Bai-Silverstein inequalities.

The theoretical results are backed by experiments on real datasets showing the practical benefits of the approach.
strength_weakness: strength_weakness: 1. Experiments are conducted on only a few datasets. And it does not explore a wide range of problem sizes or distributed computing scenarios.
2. Primarily focused on least squares regression, with limited discussion of extensions to other problems.
3. Some of the theoretical results rely on assumptions (e.g., about leverage score approximation) that may not always hold in practice.
questions: questions: 1. Are there natural extensions of this work to other loss functions beyond least squares?
2. Do you expect the techniques developed here to be applicable to other sketching problems beyond least squares? If so, which ones?

Minor Comments:

1. The abstract could more clearly state the key theoretical results/bounds achieved
2. Some additional discussion of practical implications and potential applications would be valuable
3. A few typos noted (e.g. line 330)
limitations: limitations: yes
paper_summary: paper_summary: The paper studies the least squares regression task and improves the space and communication amount in a distributed setting to be independent of $\epsilon$. The key to achieve this is by sketching the data in blocks that have an $\epsilon$ dependence which can be reduced to only d-dependencies by aggregating them. For the covariance A^TA this becomes only a d x d matrix and for the other required term A^Tb it can be handled by communicating and aggregating just the solutions which have a dimension of d. Although the idea seems very simple described this way, the analysis seems highly non-trivial and has very interesting aspects and novel techniques (or at least new to me). It analyzes in this setting not only the standard least squares error but also a bias term, which allows for smaller sketch dimensions, and only 1/sqrt(eps) dependencies.
strength_weakness: strength_weakness: * motivation of the model, seems to be a niche
* very low improvements but there seems no big gap that can be leveraged
questions: questions: * I am slightly confused by the definition of the bias. Is $||AE(\tilde x)-b||$ standard in some literature? I would say this is the variance of the expected sketched estimator. Why would I be interested in this expected estimator, instead of the actual outcome after sketching?
* I would rather define the bias as $||\tilde x -x^*||$ which I think does not allow for any improvements over standard sketching results, right?
limitations: limitations: f
paper_summary: paper_summary: Sketched least squares involve estimating the term $(X^TX)^{-1}$ which has a high bias when the sketch matrix $S$ is not sub-Gaussian. This paper gives a sparse sketching method using a LESS embedding which runs in optimal space and current matrix multiplication time, where $S$ is sparse, and constructed based on the leverage scores of the data matrix $A$ (Definition 2 $(s,\beta_1,\beta_2)$-LESS embedding). The paper also improves the sharpness of the probability bounds which is applicable to similar problems in RMT using LESS embeddings. For $s=1$, nothing is different, but when $s > 1$, the bias bound is reduced.
strength_weakness: strength_weakness: 1. It would be nice to have experiments with other sketching methods to (empirically) justify some statements on the bias, variance (although not stated) and computational time, e.g. estimators mentioned in Table 1. For example, are there sketching estimators with higher bias, but less variance? 

There are some minor typos, e.g. line 212 reference missing, lines 282 to Equation (2) at bottom of the page is missing a bracket for the numerous expectations (in contrast to lines 542 onwards in the appendix).

In Appendix A, notation for concentration inequalities should be looked at and made consistent, e.g. Lemma 6 / H{\"o}lder's inequality should have a $\frac{1}{q}$, Lemma 7 ($\lambda$ max isn't defined), Lemma 8 / Azuma's inequality ($\lambda$, $m$ should be consistent), Lemma 10 is missing a bracket for $\mathbb E[x_i^2]$. 

The presentation of the proof for Theorem 5 was slightly jarring (due to Lemma 11, Lemma 12 appearing in the proof), but there also doesn't seem to be a good way to include them (since referring to the two lemmas requires the upper bounds, and flipping back a page is also inconvenient). Maybe a solution is to indent the Lemmas, or box them up?

The dot before line 626 (after 72) should be removed.
questions: questions: 1. It would be nice to have experiments with other sketching methods to (empirically) justify some statements on the bias, variance (although not stated) and computational time, e.g. estimators mentioned in Table 1. For example, are there sketching estimators with higher bias, but less variance? 

2. Despite being clear to read, I had to go back and forth a bit to find out what the novelty is. I appreciate the clarity and thoroughness of explaining the bounds on the bias and variance, runtime, but I would like it if Definition 1 & 2 came much earlier (or at least maybe an informal Definition 2 after line 85?)
limitations: limitations: N/A
paper_summary: paper_summary: Sketching is a technique from randomized numerical linear algebra which compresses an input matrix $A$ by multiplication with a random matrix $S$. Recent works [18] have shown how to characterize the bias in the sketched inverse covariance $(\tilde A^\top\tilde A)^{-1}$, where $\tilde A = SA$ for an input matrix $A$ and a sketching matrix $S$. This work extends this machinery to construct approximately unbiased estimators $\tilde x$ for the least squares regression problem $\min_x \|Ax-b\|$. As the authors show in Theorem 2, this gives an efficient distributed algorithm by averaging the unbiased estimator over $q = O(1/\epsilon)$ servers. The techniques introduced by the authors also improve prior results on bias-free estimation of the inverse covariance.
strength_weakness: strength_weakness: The result and techniques seem to be a bit niche and incremental. The notion of distributed computation is different from ones that I see in other works in the sketching literature (e.g. https://arxiv.org/abs/1408.5823, https://arxiv.org/abs/1504.06729) where the input is partitioned across multiple servers. Instead, this work considers a setting where there is just one stream which contains the entire input, and the central server has access to $q$ machines that can each access this stream (see the Computational model section), which seems nonstandard. The number of servers needed is also rather large (e.g. $q = O(1/\epsilon)$), which also seems restrictive. I encourage the authors to include references that prove results in this setting if there are any others. The main technical novelty (line 275) also seems to be an improvement to a bound in the analysis of LESS embeddings from prior work, while the overall sketching framework is largely unchanged.
questions: questions: - The authors may consider discussion the related work https://arxiv.org/abs/2203.09755 on distributed least squares and https://arxiv.org/abs/2007.01327 on using debiasing approaches

Other comments

- Lemma 2 should have a probability statement (e.g., with probability at least $1-\delta$)
limitations: limitations: The authors claim that the work does not have limitations in the checklist. Some of the points I raised in the weaknesses section could be discussed in the work, such as restrictions on the computational model and the number of servers.

Comments: comment: Thank you for your response on the questions! After reading the other reviews, I would like to keep my current score of 5.
Thank you for the clarifications. I will keep my score.
Thank you for the rebuttal. I appreciate the additional comments which show how the limitations on the number of servers can be addressed with a smooth trade-off, as well as allowing the data to be partitioned if the chunks are partitioned uniformly (although perhaps it still doesn't work if the data is partitioned adversarially, which is still a limitation in my mind). Thank you also for the explanation of connections with the prior work. I have increased my score, in recognition that there is a growing body of work on RMT-style arguments in the sketching literature, and thus this work may have a broad impact in many future works.
Thank you for the rebuttal, and for providing the experiments. I was thinking more of biased estimators with lower variance though. My score remains unchanged.