Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton
Authors: Chengmei Niu
Zhenyu Liao
Zenan Ling
Michael W. Mahoney
Venue: ICML 2025 oral (2025)
Keywords: Inversion bias
random matrix theory
randomized numerical linear algebra
random sampling
sub-sampled Newton
Primary area: theory
TL;DR: Use RMT to characterize inversion bias for random sampling, propose debiasing, and apply to establish problem-independent convergence rate for SSN.

Abstract: A substantial body of work in machine learning (ML) and randomized numerical linear algebra (RandNLA) has exploited various sorts of random sketching methodologies, including random sampling and random projection, with much of the analysis using Johnson--Lindenstrauss and subspace embedding techniques.  Recent studies have identified the issue of *inversion bias* -- the phenomenon that inverses of random sketches are *not* unbiased, despite the unbiasedness of the sketches themselves. This bias presents challenges for the use of random sketches in various ML pipelines, such as fast stochastic optimization, scalable statistical estimators, and distributed optimization. In the context of random projection, the inversion bias can be easily corrected for dense Gaussian projections (which are, however, too expensive for many applications). Recent work has shown how the inversion bias can be corrected for sparse sub-gaussian projections. In this paper, we show how the inversion bias can be corrected for random sampling methods, both uniform and non-uniform leverage-based, as well as for structured random projections, including those based on the Hadamard transform. Using these results, we establish problem-independent local convergence rates for sub-sampled Newton methods.

Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton
Chengmei Niu Zhenyu Liao|
Zenan Ling'
Michael W. Mahoney
Abstract
As pivotal technique in RandNLA, random sketching which encompasses both random projection and random sampling is becoming increasingly critical in many mod- ern large-scale machine learning (ML) applications.
A substantial body of work in machine learning (ML) and randomized numerical linear algebra (RandNLA) has exploited various sorts of random sketching methodologies, including random sam- pling and random projection, with much of the analysis using Johnson-Lindenstrau and sub- space embedding techniques. Recent studies have identified the issue of inversion bias the phe- nomenon that inverses of random sketches are not unbiased, despite the unbiasedness of the sketches themselves. This bias presents challenges for the use of random sketches in various ML pipelines. such as fast stochastic optimization, scalable sta- tistical estimators, and distributed optimization. In the context of random projection, the inversion bias can be easily corrected for dense Gaussian projections (which are, however, too expensive for many applications). Recent work has shown how the inversion bias can be corrected for sparse sub-gaussian projections. In this paper, we show how the inversion bias can be corrected for ran- dom sampling methods, both uniform and non- uniform leverage-based, as well as for structured random projections, including those based on the Hadamard transform. Using these results, tablish mroblem-indepencee local convergence rates for sub-sampled Newton methods.
More precisely, for a tall matrix A Raxd with n random projection proposes to obtain sketch À E Rmxx of size < n of A by randomly and linearly combin- ing the rows of A; while random sampling, on the other hand, carefully selects small subset (of size say) of the rows and rescales them to obtain E TRmxd For both approaches, follows from (JL) type analysis (Johnson Lindenstrauss, 1984) that the random sketch A can be used "proxy" of in many downstream ML tasks, e.g., A A A are close in some sense with high probability, even with < n. This leads to a significant boost in the running time, communica- tion time, and/or memory cost, for many numerical meth- ods (Drineas et al., 2006b; 2011; 2012; Avron et al., 2017: Dereziński Mahoney, 2021; Lacotte & Pilanci, 2022). See also Mahoney (2011); Halko et al. (2011); Woodruff (2014); Drineas & Mahoney (2018); Martinsson & Tropp (2020); Dereziáski & Mahoney (2021); Murray et al. (2023); Dereziński & Mahoney (2024) and reference therein for an overview of RandNLA and the applications in modern ML.
Despite this promising achieved with random sketching, in many MI. pipelines ranging from linear/ridge regression scalable statistical estimation and fast stochastic optimization, the object of di- rect interest is the sketched matrix inverse (A - for some (perhaps all-zeros or diagonal) positive (p.s.d.) C (instead of Á Á itself). Since the matrix inverse is a nonlinear operator, the fact that À T Á is an unbiased or nearly unbiased estimator of A A (i.e., EIÀ A = A A, which is major guide forming the sketch Á) does not, in general, imply the unbiasedness of its inverse, i.e., " 인 (A A)-1. This phenomenon of inversion
cemplexity-accuracee trade-off
1. Introduction
Randomized numerical linear algebra (RandNLA) signifi- cantly reduces computation, communication, and/or storage overheads by using randomness as an algorithmic resource.
'School of Electronic Information & Communications, Huazhong University of Science & Technology. Wuhan, Hubei, China. *ICSI, LBNL, and Department of Statistics, University of California, Berkeley, USA. Correspondence to: Zhenyu Liao <zhenyu_liae liao hust.edu.en>.
[[ A) bias has been long known in the literature: in the case of Gaussian random projection with A = SA for S e R having i.i.d. (N 1/m) entries, (ÀT is known to follow the inverse Wishart distribution (Haff, 1979) with ({ A) ever, much less known beyond the Gaussian setting, Build-
Proceedings of the 42nd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 20255 by the author(s).
m-d-i (A T A) for m d 1. How-

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
ing upon recent progress in non-asymptotic random matrix theory (RMT), it has recently been shown by Dereziński et al. (2021b) that a similar inversion bias holds, and that it can be corrected, with E[(A T the sense of partial order of p.s.d. matrices, for sub-gausse -aassaassiinn and the so-called LEverage Score Sparsified (LESS) pro- jections. The latter can be of greater practical interest than, e.g. dense Gaussian or sub- sub-gaussian projections, as it is sig- nificantly sparser and can thus be evaluated more efficiently.
1.2. Related Work
Inversion bias. Given a matrix A R the ma- trix inverse (A (ATA) AA)= is fundamental in ML, numeri- cal computation, and statistics. Examples include lin- ear functions methods (Boyd & Vandenberghe, 2004), quadratic forms af (ATA) A 'a, (with a. the 2thh row of A) in computing matrix leverage scores (Drineas et al., 2012), and trace forms, tr L(A A) AA for some given L, of interest in uncer- tainty quantification (Kalantzis al., 2013) and experimen- tal designs (Pukelsheim, 2006). In the case of tall matrices with n > d, random sketching applies to efficiently reduce the computational overhead of (A A) A)- by using a sketch Á = SA of A for random matrix S E RRnx with m & n. It has been shown recently (Dereziński et al., 2020; 202 la:b) that these sketched inverses are biased for unstructured sub- gaussian S, with [E A } &6 (A T A-- In this paper, we consider the (often practically more interesting) case of structured random matrix S. including random sampling matrices (Definition 2.1) and randomized Hadamard trans- forms (Definition 3.6).
m.. (A A) in
m--d
AA) that are crucial to Newton's 's
This precise characterization of inversion bias has direct implications for RandNLA and ML. As a telling example, Dereziński et al. (2021a) applied LESS sketches to New- ton Sketch (Pilanci & Wainwright, 2017) and showed that this NeSS approach enjoys almost the same local convergence rates as Newton Sketch with dense Gaussian projections. This leads to a significantly better "complexity- convergence" trade-off than the vanilla Gaussian projections in stochastic second-order optimization methods.
1.1. Our Contributions
In this paper, we consider the inversion bias of random sam- pling, including uniform and non-uniform sampling, as well as structured random projections such as the Subsampled Transform (SRHT) (Ailon & Chazelle, 2006). The analysis framework in Dereziński et al. (2021a;b) does not apply (to get non-vacuous results; see, e.g., Proposition 2.8 and Remark 2.9 below). Instead, we exploit novel and non-trivial connections between non- asymptotic RMT and RandNLA to show that this inversion bias can be precisely characterized and numerically cor- rected. We also show how this inversion bias result can be used to improve the local convergence rates of the popu- lar sub-sampled Newton (SSN) method (Yao et al., 2018; Roosta- Roosta-Khorasani Khorasani & Mahoney, 2019; Xu al., 2020).
Different from our approach that explicitly modifies the sketch to correct the bias. another line of work proposes to use shrinkage-basse bassed correction techniques (Zhang & Pilanci, 2023; Romanov al., 2024).
Randomized (Walsh-)
Random sampling. Random sampling is at the core of RandNLA (Drineas et al., 2006a; Mahoney, 2011; Ma et al., 2015; Drineas & Mahoney, 2016; Dereziński & Ma- honey, 2021; 2024), and it plays central role in fast ma- trix multiplication (Drineas et al., 2006a), approximate re- gression (Drineas et al., 2006b), and low-rank approxima- tion (Cohen et al., 2017), to name a few. It of particular interest in scenarios where the dataset massive and can- not be stored and/or computed on a single machine, e.g., the census data (Wang et al., 2018) and online network data (Deng et al., 2024). See Definition 2.1 below for a formal definition of random sampling and discussions there- after for commonly-used sampling schemes including (exact and approximate) leverage score sampling (Mahoney, 2011; Cohen et al., 2017), shrinkage leverage sampling (Ma et al., 2015), as well as optimal subsampling (Wang al., 2018; Wang & Ma, 2021; Yu et al., 2022; Ma et al., 2022). In this paper, instead of providing classical JL and subspace embedding-tyre results on random sampling, we precisely characterize (and correct) the inversion bias for a variety of commonly random sampling schemes.
Our main contributions can be summarized as follows.
1. We provide precise characterization of the inversion bias for general random sampling (in Theorem 3.1) and the corresponding de-biased approach (in Propo- sition 3.2). The proposed analysis and debiasing tech- nique hold for exact and approximate leverage-based sampling (Corollary 3.4), as well as the structured SRHT (Corollary 3.7) as special cases.
2. With this precise inversion bias result, we further es- tablish (in Theorem 4.3), the first problem-independent local convergence rates for sub-sampled Newton that approximately matches the dense Gaussian Newton Sketch scheme. Numerical results are provided in Sec- tion to support these findings.
Sub-sampled Newton. Sub-sampled Newton (SSN) meth- ods propose to approximate the Hessian in Newton's 's method using a small subset of samples, and they have been ex- tensively studied within the fields of ML, RandNLA, and optimization (Xu et al., 2016; Bollapragada et al., 2019;
2

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Roosta-Khorasani & Mahoney, 2019; Xu et al., 2020; Ye al., 2021). Although these fast optimization methods are easy to implement, their convergence rates are challenging to analyze. Existing results often depend on the Hessian condition number or the Lipschitz constant and fall short of, e.g., the problem-independent convergence rates achieved by sub-gaussian Newton Sketch (Lacotte & Pilanci, 2019; Dereziński et al., 2021a). In this paper, we establish the first problem-indes independent local convergence rates for SSN that closely align with Newton Sketch. This addresses the convergence guarantee gap identified in Iterative Hessian Sketch (Pilanci & Wainwright, 2016) for random sampling.
2. Preliminaries on Random Sampling
In this section, we introduce a few definitions that will be used in the remainder of this paper.
Definition 2.1 (Random sampling). For a matrix A E RRaxd with n 2 d, a sketch A E <R of A can be constructed by randomly sampling with replacement t from the n rows of A with an importance sampling distribution, {m²)== Li=1 Ti = 1, and then rescaling by 1/ //V, This can be expressed as A = SA, with random sampling matrix S E R having only one nonzero entry per row.
includes commonly-used random sampling schemes such as uniform sampling (with _; = 1/n1), row- norm- based sampling (with मं = || / Il), exact or approximate leverage and ridge leverage score sam- pling (Mahoney, 2011; EI Alaoui & Mahoney, 2015) defined below, as well as a mix between them, e.g., the shrinkage leverage sampling (Ma et al., 2015).
Definition 2.1
Random matrix theory (RMT). RMT studies the (lim- iting) eigenspectra of large-dimensional random matri- ces (Anderson et al., 2010) and finds its applications in signal processing and communication (Couillet & Debbah, 2011), statistical finance (Plerou et al., 2002), optimiza- tion (Paquette et al., 2021; 2023), and more recently in large-scale ML (Pennington & Worah, 2017; Fan & Wang, 2020; Mci & Montamari, 2022; Couillet & Liao, 2022). A re- cent line of work (Liao et al., 2020; 2021; Liao & Mahoney, 2021) has highlighted non-trivial connections between RMT and RandNLA that this paper further develops.
Definition 2.2 (Leverage score sampling, Mahoney (2011)). For a matrix A E Raxd of rank d with n 2 and p.s.d. matrix C E Rdxd the 7th leverage score 6f of A given C, is defined as e = a. (AT A + C)-aa, E {I....,n] n]. The exact leverage score sampling refers to the random sampling approach in Definition 2.1 with n. E ///ett. for def = E;= E the effective dimension of A given C.
1.3. Notations and Organization of the Paper
The leverage score sampling has been extensively studied in RandNLA and ML. By taking C O₂ in Definition 2.2. we obtain the standard leverage scores (Mahoney, 2011; Drineas et al., 2012); and by taking C ALa, we ob- tain the X-ridge leverage scores (El Alaoui & Mahoney, 2015). Given A E RRadd its leverage scores can be ap- proximately computed in ((nz((A) (A) log n + d'og d)² + dl log n) time, for anz(A) the number of non-zero en- tries in A, see Drineas et al. (2012); Clarkson & Woodruff (2017); Cohen al. (2017).
We denote scalars by lowercase letters, vectors by bold lowercase, and matrices by bold uppercase. For a matrix A E Raxdd we denote AT, a!" E RR, and |A|| the trans- pose, 2th row, and spectral norm of A, respectively. We denote A B if B A is p.s.d., and use La for the iden- tity matrix of size d. For vector E Re and matrix with the conven- tion ||x|| = | For a random variable x, we denote B[c] the expectation of x and [<[x]] the expectation of x, conditional on the event Ç. We use 6(.), and (-) notations as in standard computer science literature.
B E Bdddd we denote ||x|| =
We also introduce an "approximation factor" to measure the extent to which one importance sampling distribution approximates another importance sampling distribution (the exact leverage score distribution, in our case).
The remainder of this paper is organized follows. Sec- tion 2 presents preliminaries random sampling and coarse-grainet characterization of its inversion bias, by di- rectly (and naively) adapting the proof approach from Derez- iński et al. (2021b) (which turns out to be vacuous in our setting). Section 3 delivers a fine- fine-grained analysis of this inversion bias and proposes an efficient non-vacuous debi- asing approach. Section 4 demonstrates how these technical results apply establish problem-indepentee independent local conver- gence rates for SSN. Section 5 provides numerical results that support our theoretical findings. Section 6 provides a conclusion, summarizing our findings and discussing fu- ture perspectives. Additional material can be found in the appendices.
Definition 2.3 (Importance sampling approximation factor). For a matrix A E gnxd with 2 d, p.s.d. matrix C e Rdxd and a random sampling matrix S E TRmx as in Def- inition 2. with importance sampling distribution {mr) the (min and max) importance sampling approximation factors of the random sampling scheme S is defined as the pair (Pmins Pmax),, with Pmin = minigisn 4 / (T,dee) and Pmax = maxigisn ( (T:deet)
For us, the importance sampling approximation factors (Pmin Pmax) in Definition 2.3 provide qualitative charac- terization on how the random sampling scheme under study
3

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
differs from the exact leverage score sampling in Defini- tion 2.2. This extends the classical notion of sampling approximation factor in Drineas et al. (2006a) to include both the maximum and the minimum. While prior work primarily focuses on the max factor, here we focus on the inversion bias, where the min factor also plays natural and significant role; see below in Section 3 and also Ma et al. (2015), who first noted the importance of the min factor in statistical style analysis. By the (generalized) median inequality, we have Pmin S 1 S Pmax, with equality for exact leverage score sampling.
tion of "unbiasedness" specifically refers to inversion bias. Furthermore, while subspace embeddings automatically en- sure e-unbiasedness up to the same level of error e, they are generally not guaranteed to remain unbiased for a smaller e (Dereziński et al., 2021b).
With these definitions and notations at hand, we are ready to assess the statistical properties of random sampling. A first quantity of interest to the design of random sampling is m, the number of trials needed to construct (c,8)- subspace embedding, for some given importance sampling distribution {T A slightly more general result is given as follows.'
It follows from Definition 2.1 that E[ĀTÁ] = ATA, so that the randomly sampled matrix ÄT Á is an unbiased estimator of the true A A. This, together with controls on the higher-order moments of A T'A A, allows one conclude that A A fluctuates, with high probability and within small "distance," around the true A A of interest. This can be made precise using the relative error approximation (for scalars and matrices) defined as follows.
Lemma 2.7 (Subspace embedding for random sampling). Given A E Raxd of rank d with n. 2 d and p.s.d. C E Rdxd, let S be a random sampling matrix with number of trials m and importance sampling distribution { as in Definition 2.1, and let dest ttrAAAA) be the effective dimension of A given C with Ac = A C)-1/2 + Then, there exists C > O independent of n, deft such that for m 2 Cuuasdell log(den/ /8) failure probability 8 e (0, 12222 e > 0, and Pmax in Definition 2.3, ATS SAC is an (c, 8) -approximation of ATACc.
Definition 2.4 (Relative error approximation), For a non- negative scalar ă 2 0, we say ă is an e-approximation of (another scalar) x. denoted ä me x; if
(1+c)-²x 2 S (1 c)z.
(1)
The proof of Lemma 2.7 uses standard matrix concentration techniques and is given in Appendix B. Note that Lemma 2.7 includes existing results of both leverage (C 0,) and ridge leverage score (C- = Ala)) sampling as special cases, see Chowdhury et al. (2018, Theorm 3).
For & being random, we say 3 is an (c, 6) approxxaaatio of x if (1) holds with probability at least 1 8. Similarly, for a p.s.d. matrix X, we say X is an e-approximation (or an (e, approxxaaaiiinnnnnn when being random) if
X X + (1+e) x X < (1++)X.
(2)
With Lemma 2.7, we are now ready to evaluate the inver- sion bias of random sampling. Since the matrix inverse is nonlinear, for A T SA with B[A S SA] = ATA A, one should, a priori, not expect that (A " (A) -1 is an unbiased or nearly unbiased estimator of (ATA) the following result, we show (by adapting, in an almost straightforward fashion, the scalar debiasing proof approach of Dereziński et al. (2021a:b) that this inversion bias can be corrected, but only to some extent, using the same scalar factor as for uub-gaussian or LESS projections. The proof of Proposition 2.8 is given in Appendix C for completeness.
Remark 2.5 (Subspace embedding). For A sketch of A, the property that A A A R. A A holds with probability 1- 8 is known as (e, -suuspaee embedding RandNLA. This concept was introduced by Drineas et al. (2006b); see also Mahoney (2011) for a history. It was subsequently used in data-oblivious form by Sarlós (2006); Drineas al. (2011), and then popularized in data-oblivious form (and mis-attributed to Sarlós (2006)) by Woodruff (2014). It plays a central role in the statistical characterization of ran- dom sketching techniques.
In
Proposition 2.8 (Coarse-grained debiasing of random sam- pling). Given A E Raxd of rank d with n 2 d and p.s.d. C E Rdddd let S E R.xx be a random sampling matrix with importance sampling distribution { as in Defini- tion 2.1 and max importance sampling approximation factor Pmax in Definition 2.3. Then, there exists C' > 0 indepen- dent of n. deft such that if m 2 C Crmnse (log(det) /8) + / /dot/// with ō < m mn mn-deg
The focus of this paper is to go beyond the subspace embedding-tye results in Definition 2.4 and Remark 2.5, and to assess the inversion bias of the form E[(A T A + -)-11 (versus the true inverse (A T C) To this end, we need the following measure of unbiased estimators. Definition 2.6 (Unbiased estimator). We say a random p.s.d. matrix X is an (e, d) -unbiased estimator of X if, conditioned on an event Ç that happens with probability at least 5,
m A A SA + -1 is an
(e, 8) -unubiide estimator of (A A CC
(1+c) = [X] = (1++)X, and X 3 01)x. (3)
Remark 2.9 (On Proposition 2.8). While the debiasing factor
Herr we present the subspace embedding result in Lemma 2.7 the (regularized) matrix Ac = A(A A + C) This is of direct use in analyzing sub-sampled Newton methods in Section
Note that the error parameter é in subspace embeddings quantifies spectral approximation error, whereas the no-
4

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
is the same as that proposed for random (e.g., sub- gaussian LESS) projections (Dereziński al., 2021a;b), the resulting inversion bias is significantly larger. In par- ticular. we have, in the case of Proposition 2.8 and for m = (pmaxelet log des)). that (m-dam
Heuristic derivation of T'Thherem 3.1. For a more transpar- ent understanding of the self-consistent equation in (4)
in m-day
of Theorem 3.1, we provide here a heuristic deriva- tion. The detailed proof of Theorem 3.1 is deferred to Appendix D. Denote x! = e, A.A/viT Tr,- Q =
A S SA + C)-1 has an inversion bias of order O (Vy / log dott), This is vacuous bound. It follows from Lemma 2.7 that for the same +
m
(A SA +
Q- = (2j C)-1 = L m x,
[ E _ + +C)= and for which we have Σ a,a? AA Then, we fol-
m
+ +0
choice of m. (A (ATS C) - is, without the debiasing factor, an (O(1), approxxmmatonoonn of (A A C) and thus has an inversion bias of order ((1). See Lemma B. in Appendix B for a proof of this fact.
Σ m [[xx,x low the deterministic equivalent framework (see, c.g., Couil- let & Liao (2022, Chapter 2) for an introduction) and show that [[Q - H-1 = 0, for H A DAA + C, for D E Ruxx given in (4). First, note that || - 1|
From Remark 2.9, it appears that the inversion bias result in Proposition 2.8 is disappointing: random sampling, in con- trast with sub-gaussian or LESS random projections, while being numerically attractive and easy to implement, does not lead to a small inversion bias, at least with the m-dwit debiasing factor and the proof approach in Dereziúski et al. (2021 la:b) under Proposition 2.8. One may thus wonder:
E[QA S SAJH I| it then fol- lows from Averman-Morrisee formula (Lemma A.3) that
[[Q A DAAH A
m
xX,x H-1
S SAJI
E
1 X
Q xxx//x m
Qa H
E
Is it possible to get sharper control on the inver- sion bias of random sampling, either by introduc- ing different debiasing scheme and/or by using more refined proof than Proposition 2.8?
+ a Q-, /m;;
Using the rank-one perturbation lemma of matrix inverse. see, e.g., Silverstein & Bai (1995, Lemma 2.6), we obtain
Qa,a H-
Below, we show that such improvement indeed possible.
EQAT STATSSSAHH
E
1 1+ H
a, a//mm,
3. Fine-grained Analysis of Inversion Bias for Random Sampling
E
QQ A DAH
for D diag diag(ma:/mmm + H 1 This leads to the self-consistent equation in (4) of Theorem 3.1.
We have seen in Proposition 2.8 and Remark 2.9 that the scalar debiasing and the proof approach in Dereziński al. (2021a;b) do not lead, in the case of random sampling, non-vacuous small inversion bias. In the following result, we provide fine-grained analysis of the inversion bias of random sampling (finer than that in Proposition 2.8), and we show that the inverse (A T SA + C)-1 for random sampling S biased in more involved fashion than random projections studied in Dereziński et al. (2021a:b)
Theorem 3.1 says that the (conditional) expectation Es SA S "I instead of being close to (A A + C)-1, is in fact close to (A DA DADAA- + pending on A and the random sampling scheme per m and {r in an implicit fashion. While seemingly un- interpretable and unusable at first sight, Theorem 3.1 can be tuned to design a de-biased random sampling approach. This is given in the following result.
C)-1, with D de-
Theorem 3.1 (Inversion bias for random sampling: fine- grained analysis). Given A E Raad of rank d with n 2 d and p.s.d. C E daxd let S E Rm be random sampling matrix with importance sampling distribution { ) in Definition 2.1 and (Pusin Piaax) as in Definition 2.3. Then, for diagonal matrix D = diag Dx the solution to?
Proposition 3.2 (Fine- trne--gained debiasing for general random sampling). Under the settings and notations of Theorem 3.1, for C the it leverage score of A as in Definition 2.2 and standard random sampling matrix S as in Definition 2.1, define the de- de-biased sampling matrix S E pmxx as
m
DE = m a. (ATDA + C)-la,/
(4)
$_diag { (mm/mm - (ammmmmmmmmm ((/TT,))
m
S. (5)
there exists C > 0 independent of n, deit so that for m M Cumux dest (lo 1 < m -3, (A (ATSTSA+ +
s=1
+
Then, there exists constant C > 0 independent of n, deit such that for m 2 C Cp (oo(ddot/) + 1/22/), S
c)-1 is an (e, -Junbiased estimator of (A (ATDa DA -1
II can be checked that 20t In 3 D m4 deff In
m 33 (A T STAA + C) is an (e, 6) --unbiase estimator
with Paiin Pmax in Definition 2.3. See Lemma D.3 in Appendix D.
of (A A + C)
5

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Heuristic derivation of Proposition 3.2. To make the intu- ition behind Proposition 3.2 more accessible, we present here heuristic derivation of (5). We refer the reader to Appendix E for the detailed proof of Proposition 3.2. Let STS = E==1 Fe,i, e,,e_// /mi, for some de- terministic Fi to be specified, Q = (ATSTSA SA + + C) -1 and similarly 0)- = ( Em x
ef in (5). In Corollary E.1 of Appendix E, consider approximate leverage scores (which are much faster com- pute (Drineas et al., 2012; Clarkson & Woodruff, 2017; Cohen et al., 2017)). We show that for given sampling scheme { replacing exact leverage scores with their approximate counterparts in the de-biased sampling matrix $ in (5) increases the inversion bias, but only very slightly.
Q- = ( [n 2izs XIX] + +C- 1 as in the heuris- tic derivation of Theorem 3.1 above. We thus have 1
Note from Proposition 3.2 that the proposed fine-grained de- de-biasing matrix $ depends on the importance sampling distribution only via & / /T.. (See Appendix E. for the RMT intuition on how the exact leverage scores arise from the derivation.) As such, for any random sampling method with T. R / /dett close to those of exact leverage score sampling in Definition 2.2, we have S ~ m S. This coincides with the scalar debiasing scheme in the coarse-grainet result of Proposition 2.8, but it has a much smaller inversion bias. This special case is discussed in the following result, proven in Appendix E.3.2.
m =.==1 [[x,] = En= F. Our goal is to determine Q (and FR)) such that [E H ~ 0, for H- = (A A + C)- To this end, observe that
[[Q] - H TI--BE-O EQQAA AHHATAHH' EIQA SOATSTSAJHE SAJH = 0.
deft
By Sherman-Morrison formula (Lemma A.3), we obtain
Fii, TT, cT, /T,//,, AH 1 1 Fi,, 0 AQ A ATe TTc.//ma /mx,.
E[QATST SRQATSTSAHH-LS
Corollary 3.4 (Inversion bias using scalar debiasing un- der approximate leverage). Under the settings and nota- tions of Theorem 3.1, for random sampling scheme with sampling distribution ni E [ /[dett E. Pddet) with Pmin E [1/2, 1] as in Definition 2.3,² there exists C C> >>> 0, 2 loga (ogg(det//)) 5 < m m3 such that m
Q- T e,e, AH
E
1=1 1+ Paee AQ- A ATe. es/ es/mm
where we see the exact leverage score e, AHAH-IAT ATe a. (A A + Ca a, e as in Definition 2.2 naturally appears in the denominator from the derivation. Invoking the rank-one perturbation lemma once more, we have that
A SA + C)-1 is an for m 2 (Pt" (m-das A S
(c.n)--n estimator of (ATA A + C) with inver- sion bias c = max max(O(d of
3//22 ((cptt)))
and cp =
ST SAJH-''B[OO Σ
F. AH
max (pmi - 1, - pmax)
I Fall Fall//mm
Remark 3.5 (Inversion bias for exact versus approximate leverage score sampling). It follows from Corollary 3.4 that for exact and/or approximate leverage score sampling with Pmaxx 2 Pmnin/ Pmmin/(2pmin - 1) 2 and Pmin > (1 + O(drr > 1/2 (so that e, 1 - pmmae). the inversion bias induced by the scalar debiasing m-d- establishes the following phase transition behavior:
=E[QJATAHE AH
where we take the debiasing factor Fu mT./ (m/; 1
such that (1 + F mmr) = 1. This leads to the form of the debiasing matrix S as in (5) of Proposition 3.2. D
Comparing the fine-grained results in Proposition 3.2 to the coarse-grained results in Proposition 2.8, we see that the large inversion bias in Proposition 2.8 is indeed conse-
1. if the random sampling scheme is sufficiently close to exact leverage sampling, in that Pmax E 1), - O(d=r" (or equiva- lently the importance sampling probabilities satisfy (44{ then the inversion n; E [(1= =
quence of the proof approach adapted from Derezińiki et al. (2021a;b), that is inadequate for random sampling and for structured random projections such as the SRHT.
|Pain// /2pmi/
Remark 3.3 (Š as a random sampling scheme). Note that S E Rmxx in Proposition 3.2 is nothing but another ran- dom sampling matrix: it features exactly one nonzero entry per row that is equal to (mT,, - (9)-1/2, as opposed to (mTT,)-1 -1/2 for the standard random sampling S in Defi- nition 2.1. This non-standard re- weigghhting (that uses the leverage scores of A) ensures that (A S SA + C)-1 is a nearly unbiased estimate of (A A + C -1, per Proposi- tion 3.2. From a computational perspective, the re- re-weighted random sampling Š in (5) can be computationally demand- ing due to the need for exact computation of leverage scores
bias under scalar debiasing is the same as that (of the fine- iine-grained matrix debiasing) in Proposition 3.2; but
2. if the random sampling scheme significantly devi- ates from exact leverage sampling with fmax > O (2) (or equivalently | - / ///dttt > ((d=v/3)) / /dea), then the inversion bias under scalar
"Note that by Definition 2.3 we have Pmiin A (F(((r(aa))) S Piaxx for all i, which, together with Pmnin 2 1/2 yields that |T: ( /deetl S E /den S e /deett.
6

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
debiasing becomes larger than that in Proposition 3.2, increases with Pmaxx and saturates at Puna (1))
for uniform random sampling matrix S E Rmx* T. = 1/n as in Definition 2.1, H,, TR the Walsh-Hadamard matrix of size n. and diagonal Dn E gnxa having i.i.d. Rademacher random variables on its diagonal.
This phase transition behavior is visualized in Figure 1, See also Figure in Appendix G for the numerical comparison of inversion bias using scalar debiasing between exact and approximate leverage score sampling.
The SRHT in Definition 3.6 enjoys the following proper- ties: the randomized of A is known to have approximately uniform leverages scores, that is (HH,, D,,,,,,, D,,A/ & /n, see Drineas et al. (2011) and Tropp (2011, Theorems 3.1 and 3.2), as well as Lemma E,4 in Appendix E.3 in our setting; and since H. Hn/n = I, and D = In. one has A D, n A have the same effective dimension. These lead to the fol- lowing fine-grained debiasing result for SRHT with scalar debiasing, proven in Appendix E.3.4.
transform H,,
Inversion bias F
o )
A = ATA, so that A,,,,,/ V and
9(1)
XPaun -I
1-
Pmax
Corollary 3.7 (Fine-grained debiasing for SRHT us- ing scalar debiasing). Under the setting and notations of Theorem 3.1, for ASRHT E Rmxn the SRHT of A as in Definition 3.6, then there exists C >> 0. V 2 0, n.exp(- m 2 Cpmax 11 An-da m-ddn off
Figure I: The phase transition behavior of inversion bias as function of Pmax discussed in Remark 3.5 with scalar debiasing.
<5 5 < m-3 such that for m
As side remark, it is known from Derezinski et al. (2021b, Theorem 10) that for approximate leverage sampling, and any scalar y >> 0, m. > 0, (7 S SA + C)-1 with C == O₂ is not an (c, -uuubiieed estimator of (A T A + C)-1, with any e S cdett/m /m and e > an absolute constant. Thus,
A SSHT ASRHT + O)-² is an (0(d + Op /lg(n/ ( (0)))///0 1012), 8) -uuuiased
estimator of (A T A +
4. Application to De-biased Sub-sampled Newton with Improved Convergence
1. in the case of approximate leverage score sampling with Pinax 3/2 and Pmnin = 1/2, for any m 2 Cpmaxelott log detf, it follows from the proof of Corol- lary 3.4 that the inversion bias is upper bounded by (dett/m) and this coincides with the lower bound in Dereziński et al. (2021b, Theorem 10); and
In this section, we show that the precise characterizations of random sampling inversion bias and the debiasing tech- niques in Section apply to establish local convergence rates of SSN methods.
Consider the following optimization problem:
2. in the case of exact leverage score sampling with Punxx = -min = 1, the inversion bias can be made smaller than dett / mm under scalar debiasing.*
B arg min F(3) = arg min f(3) - +D()), (7) BEC
BEC
for some smooth function F; Rd + R that can be de- composed into f and b and C S Fd convex set. This decomposition naturally arises in MI. when, e.g., the loss function F is the sum of the empirical risk f over set of training samples and some regularization penalty +.
As an important consequence, Proposition 3.2 also applies to effectively de-bias another commonly-used data-oblivious sketching scheme, the SRHT (Ailon & Chazelle, 2006).
Definition 3.6 (Sub-sampled randomized Walsh- Hadamaad transform, SRHT, Ailon & Chazelle (2006)). For given matrix A E Raxd of rank d with 2 d, assume without loss of generality that = 2n for some integer p. Then, the SRHT of A is given by
Newton's 's method solves (7) by performing iterative updates B = 3, - -HHH, - ((,)), with μ4, the step size, VF(3,) E Rdd the gradient, and (3,) E Raddd the Hes- sian of F at Bg that can be decomposed as
ÄSRHT SH, RTX"
(6)
H2, = A(B)) A(B + C(3,),
(8)
'Notably, using exact (instead of approximate) leverage score sampling in the same setting of Dereziński et al. (2021b, Theo- rem 10), the inversion bias (conditioned on any event ( that ensures invertibility) can be made zero by taking tributed as Binomial(m, 1/d). This aligns with our conclusion of possibly smaller inversion bias than dea/ / /m. See Corollary E.2 Appendix for a proof of this fact.
with A(3,) Baxdd and some p.s.d. matrix C(B,)
V2 E gaxd that takes a simple form, c.g., C(3,)
2 Ex(11/b for b dis-
2 in the case of L2 regularization ((3) = \||BI|k
'Recall from Definition 2.3 that here for SRHT, we have n, /nn and thus Pmax = maxigign ni, (H, A/
7

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Despite having a locally super-linear convergence rate, New- ton's method suffers from heavy computational burden in forming the Hessian matrix H,))),, particularly when the training samples n is large, e.g., n >> d. In this case, the major computational bottleneck of Newton's method lies in the construction of ((3 (33,) for the inverse Hessian matrix, as this takes (Ind) time. Many random- ized second-orde methods have been proposed replace the exact Hessian inverse by some computationally efficient estimate. Here, we consider SSN methods, that randomly sample the Hessian (Yao et al., 2018; Roosta- KhostaaKoooasani & Mahoney, 2019; Xu et al., 2020) as follows.
Under Assumption 4.2, we evaluate the local convergence rate of the following de-biased SSN iterations:
3.41-3.-4 (A(3,) STS,A(()) + C(3.)) -1 gr, (10)
with de- S. diag {v 1//mmmm - (B,)) ( ==1 S as in Proposition 3.2, with (,,(()) the ttth leverage score of (3,) given C(B,). This leads to the following result. Theorem 4.3 (Local convergence of de-biased SSN). Let Assumption 4.2 hold. For p.r A(B") ((8)) = (33) and p.s.d. C( = (f*( there exists a neighborhood U of such that the de-biased SSN iteration in (10) start- ing from Go E EU satisfies, for U = =B3: ||3 - -3'1 || < (Pmaxdett miin/m)) L3, step size br
Definition 4.1 -Sub-sampled Newton, SSN), To solve the optimization problem in (7), the SSN method performs the following iteration:
Pranx
m/ dorr t Praax
m 2 d and 2 logder (log(denT//)) that
(A(3) S,A(3)) + C(B)) -1
gr: (9)
1 - B8 H
1/T
Pmaxdell
(1 1++0 t), (1I)
for t- 0, with g: (F,)) Rd the gradient of F at B,, Mr the step size at time 4, and random sampling matrix S, E Rmxx as in Definition 2.1.
1B3
B
1|
m
holds for e. = OOd and conditioned on an event Ç that happens with probability at least 1 1-6 8. Here, a m the smallest singular value of H = A(3 A + C(3"), Pmax is the max importance sampling approximation fac- tor in Definition 2.3 for EC deft = maxigisT dea (3,) with (B,) and dea (3,) the leverage scores and effective dimension of (B,) given C(3,), respectively.
By randomly sampling the (computationally intense com- ponent of the) Hessian matrix, each SSN step in (9) takes only Ommd) time. The (local global) convergence rates of SSN have been extensively studied in the literature of optimization, RandNLA, and ML; see, e.g., Xu et al. (2016); Bollapragada et al. (2019); Roossa-KKorasani & Mahoney (2019); Ye et al. (2021); Lacotte et al. (2021) and Section 1.2 above. In particular, Lacotte et al. (2021) proposed an adap- tive SSN-type algorithm that achieves a quadratic conver- gence rate by dynamically adjusting the sketch size, improv- ing upon the well-established qunear-quadrati rates of tradi- tional SSN methods. Incorporating their adaptive algorithm into our proposed de-biased approach could potentially ac- celerate convergence further; however, this exploration lies beyond the scope of the present work.
maxigigT (3,) and
The proof of Theorem 4.3 can be found in Appendix F.2. The proof relies on a precise characterization of the second inverse moment of the randomly sampled Hessian matrix, extending beyond the first inverse moment result in Propo- sition 3.2. Due to space limitation, this technical result is presented in Proposition F.1 of Appendix F. For the sake of practical implementation, we also present, in Corollaries F.5 and F.6 of Appendix F, local convergence rates of SSN it- erations using scalar debiasing m/ (m - der) under exact and approximate leverage score sampling as well as SRHT.
Due to the absence of precise characterizations of the sub- sampled Hessian inverse, as in Theorem 3.1 and Proposi- tion 3.2 (that allow to, e.g., prove the near-unbiasednes of SSN iteration E| + 2 B, H (B,)))) it is tech- conver- gence rates for SSN. In the following, we fill this gap by showing how our inversion bias results in Section apply to establish problem-independen local convergence rates for de-biased SSN. We position ourselves under the following standard assumption on the objective function F.
5. Numerical Experiments of De-biased SSN
nically challenging to obtain
In this section, we provide empirical evidence show- ing the improved convergence (and a better "ccoplexitt- convergence" trade-off as consequence) of different de- biased SSN methods proposed in Section 4, with respect to several first- and second-order baseline optimization meth- ods. We solve the following logistic regression problem
Assumption 4.2 (Lipschitz Hessian). Fff: I2d + RR in (7) have Lipschitz continuous Hessian with Lipschitz con- stant L. that is, for any B, B E Rdd Fa) F(F(')),, 2 f (() 2 ((3')) < L B'1).
a Elog log (1 + exp(-yAa, 3)) + +
A 11312 (12)
min BER 12
of the form (7), for regularization À > 0, where a e Rd is the thh row of data matrix A E Raxd sampled from both
'Of course, this does not need to done explicitly.
8

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
MNIST (LeCun al., 1998) and CIFAR-IO (Krizhevsky, 2009) datasets, and y [+1]² plementation details are provided in Appendix G. Note that the time reported in Figures 2 and 3 include both the input data pre-processing time (e.g,, the computation of exact or approximate leverage scores, and Haldmadd trans- form) and the computational overhead associated with the sketching process.
GD
LESS SRITT
ARLey
SLev
th the response vector. Im-
SGD
ALev
Relather erme
Relathe
10
10-3
9
10-0 0
Wall-clock time (S)
Wall-clock time (S)
Figure 2 assesses the impact of the sketch size m on both relative error ()) Br - 3 | - (3 and running time of de-biased SSN employing approximate ridge leverage score sampling (SSN-ARLev), in comparison to the Newton- LESS method (Dereziński et al., 2021a) based on random projection. The results in Figure 2 demonstrate that the proposed de-biased SSN consistently outperforms Newton- LESS across all tested sketch size m, exhibiting superior trude-off. Notably, while the run- ning time for Newton-LESs LESS increases significantly with that of SSN-ARLev remains approximately constant.
(a) MNIST data
(h) CIFAR-IOO data
Figure 3: Convergence-compllatt complexity trade-off between vari- ous optimization methods on logistic regression for MNIST and CIFAR-10 10 data, with sketch size m = 300 for MNIST and m = 400 for CIFAR-10 data. Results are obtained by averaging over 10 independent runs (except for GD that is deterministic).
comvegence-commleene
we see that SSN methods, when prop- erly de-biased, exhibit a significantly better convergence- complexity trade-off than both first-order methods and the Newton-LESS approach. The SRHT sampling outperforms Newton-LESS but still lags slightly behind all other SSN variants in speed. Among these, the SLev scheme edges out ALev yet remains slower than ARLev. Across both datasets tested, de-biased SSN with ARLev consistently delivers the optimal convergence-compler trade-off among all methods evaluated.
From Figure
Newton LESS
SSN ARLey
Wall time (s)
10--
Wall-cook time (s)
Rebaive croe
10
Relutve ernr
10 10-8
10-9
1.6
10
10
0.00
50U
1,000
500
1,000
1,500
Sketch size
Sketch sioe
(a) MNIST data
CIFAR-I0
Figure 2: Relative errors (in solid lines) and wall-clock time (in dashed lines) as a function of the sketch size m, for Newton-LESS and the proposed de-biased SSN-ARLev methods, applied to logistic regression on both MNIST and CIFAR-10 data. Relative errors are obtained after a fixed number of iterations (T 5 for MNIST data and T = 7 for CIFAR-10 data). Results are obtained by averaging over 30 independent runs.
6. Conclusions and Perspectives
In the work, we investigate the inversion bias inherent in various random sampling schemes, including uniform and non-uniform leverage-based sampling, wel struc- tured random projections (c.g., the Hadamard transform- based SRHT). Leveraging recent advances in RMT and RandNLA, we provide precise characterization of this inversion bias and propose corresponding de-biasing tech- niques. Notably, for approximate leverage sampling and SRHT, this de-biasing reduces to a simple scalar rescaling. We further show that our results enable an improved SSN method, achieving local convergence rates comparable to those of Newton Sketch with dense Gaussian projections. Our theoretical insights are complemented by numerical results on MNIST and CIFAR-10 10 datasets, underscoring the practical relevance of the proposed approach.
Figure 3 compares the relative errors as a function of wall- clock time across various optimization methods:
1. First-order baselines: Gradient Descent (GD) and Stochastic GD (SGD).
2. De-biased SSN methods using different sampling schemes: Approximate X-Ridge Leverage Score (AR- Lev), Approximate Leverage Score (ALev), Shrinkage Leverage Score (SLev) sampling (Ma et al., 2015), and SRHT (see Definition 3.6 above).
It would be of future interest to see whether the proposed debiasing technique, when wisely combined with adap- tive sampling schemes, can achieve or even improve the quadratic convergence in Lacotte et al. (2021). In addition, it is also worthwhile to extend our debiasing framework to dependent sampling methods (Cortinovis & Kressner, 2024), such as volume sampling, which may further improve SSN's efficiency and convergence.
3. Newton Sketch with LESS-uniform sketch (LESS) (Dereziński et al., 2021a), which achieves significantly shorter running times compared to the original Newton Sketch with dense Gaussian sketches.

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Acknowledgements
chine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 989-998. PMLR, 2018.
Z. Liao would like to acknowledge the National Natural Sci- 101) and the Guangdong Provincial Key Laboratory of Mathematical Foundations for Artificial Intelligence (2023B 12120 1212010001) for providing partial support. Z. Ling is supported by the National Natural Science Foundation of China (via NSFC- 62406119 19) and the Natural Science Foundation of Hubei Province (2024AFB074). M. W. Mahoney would like to ac- knowledge the DARPA, DOE, NSF, and ONR for providing partial support of this work.
ence Foundation of China (via funn
Clarkson, K. L. and Woodruff, D. P. Low-rank approxima- tion and regression in imput sparsity time. Journal of the ACM (JACM),
Cohen, M. B., Musco, C., and Musco, C. Input sparsity time low-rank approximation via ridge leverage score sampling. In Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 1758- 1777. SIAM, 2017.
Cortinovis, A. and Kressner D. Adaptive randomized piv- oting for column subset selection, DEIM, and low-rank approximation. arXiv preprint urXiv:2412.13399 2024.
Impact Statement
This paper presents work whose goal is to advance the field of MI. RandNIA, and their connection to RMT. The work is primarily theoretical, and we do not see any immediate negative societal impact. If anything, providing methods that come with smaller theory-practice gap will make it easier to identify and mitigate unintended negative impact.
Couillet, R. and Debbah, M. Random Matrix Methods for Wireless Communications. Cambridge University Press, 2011.
Couillet, R. and Liao, Z. Random Matrix Methods for Machine Learning. Cambridge University Press, 2022.
References
Deng, J., Huang, D., Ding, Y., Zhu, Y., Jing, B., and Zhang, Subsampling spectral clustering for stochastic block models in large-scale networks. Computational Statistics & Data Analysis, 189:107833 2024.
Ailon, N. and Chazelle, B. Approximate nearest neighbors and the fast Johnson-Lindenstraus transform. In Pro- ceedings of the Thirty-Eighth Annual ACM Symposium on Theory of Computing, pp. 557-563, 2006.
Anderson, G. W., Guionnet, A., and Zeitouni, O. An Intro- duction Random Matrices, volume 118 of Cambridge Studies in Advanced Mathemutics. Cambridge University Press, 2010.
Dereziński, M. and Mahoney, M. W. Determinantal point processes in randomized numerical linear algebra. No- tices of the AMS, 38(1):34 45, 2021.
Dereziński, M. and Mahoney, M. W. Recent and upcoming developments in randomized numerical linear algebra for machine learning. In Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD '24, pp. 6470-6479, New York, NY, USA, 2024. Association for Computing Machinery.
Avron, H., Clarkson, K. L., and Woodruff, D. Faster ker- nel ridge regression using sketching and preconditioning SIAM Journal Matrix Analysis and Applications, 38 (4):1116-1138, 2017.
Bollapragada, R., Byrd, R. H., and Nocedal, J. Exact and in- exact subsampled Newton methods for optimization. IMA Journal of Numerical Analysis, 39(2):545-577 2019.
Dereziński, M., Liang, F. T., Liao, Z., and Mahoney, M. W. Precise expressions for random projections: Low-rank approximation and randomized Newton. In Advances in Neural Information Processing Systems, volume 33, pp. 18272-1828 Curran Associates, Inc., 2020.
Bonnans, J.-F., Gilbert, J. C., Lemaréchal, C., and Sagas-
tizábal, C. A. Numerical Optimization: Theoretical and
Practical Aspects. Springer Science & Business Media, 2006.
Dereziński, M., Lacotte, J., Pilanci, M., and Mahoney, M. W. Newton-less: Sparsification without trade-offs for the sketched Newton update. In Advances in Neural Infor- mation Processing Systems, volume 34, pp. 2835-2847. Curran Associates, Inc., 2021a.
Boyd, S. P. and Vandenberghe, I.. Convex Optimization. Cambridge University Press, 2004.
Burkholder, D. L. Distribution function inequalities for martingales. The Annals of Probability, 1(1):19 42, 1973.
Dereziński, M., Liao, Z., Dobriban, E., and Mahoney, M. Sparse sketches with small inversion bias. In Proceedings of Thirty Fourth Conference on Learning Theory, volume 134, pp. 1467-1510. PMLR, 2021b.
Chowdhury. A., Yang, J., and Drineas, P. An iterative, sketching-baat based framework for ridge regression. In Pro- ceedings of the 35th International Conference on Ma-
10

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Drineas, P. and Mahoney, M. W. RandNLA: Randomized numerical linear algebra. Communications of the ACM, 59:80-90, 2016.
Lacotte, J. and Pilanci, M. Faster least squares optimization. arXiv preprint /9t 2019.
Lacotte, J. and Pilanci, M. Adaptive and oblivious random- ized subspace methods for high-dimensional optimiza- tion: Sharp analysis and lower bounds. IEEE Transac- tions on Information Theory, ((15):321-33000 2022.
Drineas, P. and Mahoney, M. W. Lectures random- ized numerical linear algebra. In Mahoney, M. W., Duchi, J. C., and Gilbert, A. C. (eds.), The Mathematics of Data, IAS/Park City Mathematics Series, pp. 1-48. AMS/IAS/SIAM, 2018
Lacotte, J., Wang, Y., and Pilanci, M. Adaptive Newton sketch: Linear-time optimization with quadratic conver- gence and effective Hessian dimensionality. In Proceed- ings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learn- ing Research, pp.
Drineas, P., Kannan, R., and Mahoney, M. W. Fast Monte Carlo algorithms for matrices I: Approximating matrix multiplication. SIAM Journal on Computing, 36(1):132- 157, 2006a.
PMLR, 2021.
Drineas, P., Mahoney, M. W., and Muthukrishnan, S. Sam- pling algorithms for l2 regression and applications. In Proceedings of the Seventeenth Annual ACM-SIAM Sym- posium on Discrete Algorithm. pp. 1127-1136, 2006b.
LeCun, Y., Bottou, L.,, Bengio, Y., and Haffner, Gradient- based learning applied to document recognition. Proceed- ings of the IEFE, 86(1 2-227722224 1998.
Drineas, P., Mahoney, M. W., Muthukrishnan, S., and Sar- 16s, T. Faster least squares approximation. Numerische Mathematik, 117(2):219-2) 249, 2011.
Liao, Z. and Mahoney, M. W. Hessian eigenspectra of more realistic nonlinear models. In Advances in Neural Information Processing Systems, volume 34, pp. 20104- 201 17. Curran Associates, Inc., 2021.
Drineas, P., Magdon-Ismail, M., Mahoney, M. W., and Woodruff, D. P. Fast approximation of matrix coher- ence and statistical leverage. The Journal of Machine Learning Research, 13(1):3475-355) 2012.
Liao, Z., Couillet, R., and Mahoney, M. W. A random matrix analysis of random Fourier features: Beyond the Gaussian kernel, precise phase transition, and the corresponding double descent. In Advances in Neural Information Pro- cessing Systems, volume 33, pp. 13939- 13950. Curran Associates, Inc., 2020.
El Alaoui, A. and Mahoney, M. W. Fast randomized kernel ridge regression with statistical guarantees. In Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
Liao, Z., Couillet, R., and Mahoney, M. W. Sparse quan- tized spectral clustering. In International Conference on Learning Representations, 2021.
Fan, Z. and Wang, Z. Spectra of the Conjugate Kernel and Neural Tangent Kernel for linear-width neural networks. In Advances in Neural Information Processing Systems, volume 33, pp. 7710-7721. Curran Associates, Inc., 2020.
Ma, P., Mahoney, M. W., and Yu, B. A statistical perspective on algorithmic leveraging. Journal of Machine Learning Research, 16:861-911, 2015.
Haff, L. R. An identity for the Wishart distribution with applications. Journal of Multivariate Analysis, 9(4):531- 544, 1979.
Ma, P., Chen, Y., Zhang, X., Xing. X., Ma, J., and Mahoney, M. W. Asymptotic analysis of sampling estimators for randomized numerical linear algebra algorithms. Journal of Machine Learning Research, 23(177):1-45 2022.
Halko, N., Martinsson, P.-G., and Tropp, J. A. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217-25 2011.
Mahoney, M. W. Randomized algorithms for matrices and data. Foundations and Trends@ in Machine Learning, 3 (2):123-224, 2011.
Johnson, W. B. and Lindenstrauss, J. Extensions of Lip- schitz mappings into a Hilbert space. Contemporary Mathematics, pp. 189-206, 1984.
Martinsson, P.-G. and Tropp, J. A. Randomized numer- ical linear algebra: Foundations and algorithms. Acta Numerica, 29:403-572, 2020.
Kalantzis, V., Bekas, C. Curioni, A., and Gallopoulos, F. Accelerating data uncertainty quantification by solving linear systems with multiple right-hand sides. Numericul Algorithmss, 62(4))637-653, 2013.
Mci, S. and Montanari, A. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75(4):667-766, 2022.
Krizhevsky, A. Learning Multiple Layers of Features from Tiny Images. PhD thesis, 2009.

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Murray, R., Demmel, J., Mahoney, M. W., Erichson, N. B., Melnichenko, M., Malik, O. A., Grigori, L.., Luszczek, P., Derezióski, M., Lopes, M. E., Liang, T., Luo, H., and Dongarra, J. Randomized Numerical Linear Algebra perspective on the field with an eye to software. Technical Report Preprint: arXiv:2302. 1 1474v2, 2023.
Tropp, J. A. Improved analysis of the subsampled random- ized Hadamard transform. Advances in Adaptive Data Analysis, 3(01n02):
2011.
Tropp. J. A. An introduction to matrix concentration inequal- ities. Foundations and Trends@ in Machine Learning, 8 (1-2):1-230, 2015.
Paquette, C., Lee, K., Pedregosa, F., and Paquette, E. SGD in the large: Average-case analysis, asymptotics, and step- size criticality. In Proceedings of Thirty Fourth Confer- ence on Learning Theory, pp.
Vershynin, R. High- Soah-Diiensioall Probability: An Introduc- tion with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematies. Cambridge University Press, 2018
PMLR, 2021. Paquette, C., van Merriënboer, B., Paquette, E., and Pe- dregosa, F. Halting time is predictable for large models: A universality property and average-case analysis. Foun- dations of Computational Mathematics, 23(2):597-673, 2023.
Wang, H. and Ma, Y. Optimal subsampling for quantile regression in big data, Biometrika, 108(1):99-112, 2021.
Wang, H., Zhu, R., and Ma, P. Optimal subsampling for large sample logistic regression. Journal of the American Statistical Association, 113(522):829 844, 2018.
Pennington, J. and Worah, P. Nonlinear random matrix theory for deep learning. In Advances in Neural Infor- mation Processing Systems, volume 30 of NIPS' 17, pp- 2637-2646. Curran Associates, Inc., 2017.
Woodruff, D. P. Sketching as tool for numerical linear alge- bra. Foundations and Trends in Theoretical Computer Science, 10(1 10(1-2):1-157 2014.
Pilanci, M. and Wainwright, M. J. Iterative Hessian sketch: Fast and accurate solution approximation for constrained least-squares. Journal of Machine Learning Research, 17 (-3):1-38, 2016.
Xu, P., Yang. J., Roosta-KKorasaaa Khorasani, F., Ré, C., and Mahoney, M. W. Sub-sampled Newton methods with non-uniform sampling. In Advances in Neural Information Processing Systems, volume 29 of N6,S pp. 3000-3008 Curran Associates, Inc., 2016.
Pilanci, M. and Wainwright, M. J. Newton sketch: near linear-time optimization algorithm with linear-quadratic convergence. SIAM Journal on Optimization, 27(1):205- 245, 2017.
Xu, P., Roosta- Khorasani, F., and Mahoney, M. W. Newton- type methods for non-convex optimization under inexact Hessian information. Mathematical Programming, 184 (1 (1-2):35-70,, 2020.
Plerou, V., Gopikrishnan, P., Rosenow, B., Amaral, L, A. N., Guhr, T., and Stanley, H. E. Random matrix approach to cross correlations in financial data. Physical Review E, 65(6):066126 2002.
Yao, Z., Xu, P., Khosst-Khoaasnn,, F., and Mahoney, M. W. Inexact non-convex Newton-type methods. arXiv preprint arXiv: /802.66955 2018.
Pukelsheim, F. Optimal Design of Experiments. SIAM, 2006.
Ye, H., Luo, L., and Zhang. Z. Approximate Newton meth- ods. Journal of Machine Learning Research, 22(66): 441 2021.
Romanov, E., Zhang, F., and Pilanci, M. Newton meets Marchenko- Pastur: Massively parallel second-order op- timization with Hessian sketching and debiasing. arXiv preprint arXiv22410.01374, 2024.
Yu, J., Wang, H., Ai, M., and Zhang, H. Optimal distributed subsampling for maximum quasi-likelihood estimators with massive data. Journal of the American Statistical Association, 117(537):265-276 2022.
Koosta-Khorasani F. and Mahoney, M. W. Sub-sampled Newton methods. Mathematical Programming. 174(1-2):
293-326, 2019.
Zhan, X. Singular values of differences of positive semidef- inite matrices. SIAM Journal on Matrix Analysis and Applications, 22(3):8 119-8899,,,1 2001.
Sarlós, T. Improved approximation algorithms for large ma- trices via random projections. In Proceedings of the 47th
Annual IEEE Symposium on Foundations of Computer Science, pp. 143-152, 2006.
Zhang, F. and Pilanci, M. Optimal shrinkage for distributed second-order optimization. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pp. 41523 41549. PMLR, 23-29 Jul 2023.
Silverstein, J. W. and Bai, Z. On the empirical distribution of eigenvalues of a class of large dimensional random matrices. Journal of Multivariate Analysis, 54(2):175-
192, 1995.
12

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Supplementary Material of Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton
The technical appendices of this paper organized as follows.
We recall few technical lemmas that will used in our proofs Appendix
The proof of Lemma 2.7 on the subspace embedding property for random sampling is given in Appendix B.
The proof of Proposition 2.8 on coarse-grained debiasing of random sampling is given in Appendix C.
The proof of Theorem 3.1 on fine- fineggriieed analysis of inversion bias for random sampling given Appendix
The proof of Proposition 3.2 on fine-grained debiasing for random sampling is given in Appendix E.
The proof concerning results the application to de-bbased SSN in Section 4 is given in Appendix F. Implementation details for the numerical results in Section 5 are given Appendix
A. Useful Lemmas
In this section, we introduce few technical lemmas to be used in subsequent sections.
Lemma A.1 (Singular value bounds of symmetric p.s.d. matrices, Zhan (2001, Theorem 2.1)). For real symmetric p.s.d.
matrices J, K E R=x* having (ordered) singular values 0 (J] 2 (J) 2 2 (J( and 01 (K) 2 o2(K) 2 a, (KKK,, the singular values of the difference
K are bounded as
(
0) K)
o,(J [O)-K) s oj
i1.....,
In particular, we obtain
J - K S max [|]]|, K|]]
Lemma A.2 (Intrinsic matrix Bernstein, Tropp (2015, Theorem 7.3.1)). Let X I..... X,, be n independent real symmetric random matrices such that E.X,] O, maxigign |X,|| S Pl, and Σ X I for some symmetric p. p.s.d. matrix B. Then, for any e Z ||B||2/2
12 + pI/3, we have
4r B B
2/22
Pr
X
exp
|
pie/3
Lemma A.3 SSherman-Morrsee formula). For invertible matrix A Raxx and two vectors u. v E R", A uv is invertible and only if + v - 0 and
(A + - A-1 A - Uv T -
1+ A 'u
Besides, it also follows that
A-1m u (A uv u 1+v A
Lemma A.4 (Burkholder inequality, Burkholder (1973)). For fir Pi a real martingale difference sequence with respect to the increasing a o-field F,, we have
T/2
E
Ui
< LT E
for any integer and some constant Lor > independent of
13

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
We will use Lemma A.4 establish concentration results for resolvent/inverse matrices. In that context, the increasing a-field Fi is defined with respect the independent trials of random sampling (i.e., the rows of the random sampling matrix S E 1 in Definition 2.1, see Lemma D.I in Appendix D.3 and the discussion thereafter.
Theorem A.5 (Hanson-Wrinn Wright inequality, Vershynin (2018)). Assume (II,..... ,nn IR" be random vector with independent, mean-cero, sab-gaassian variables. For matrix M E x and any t 2 0, then we have
t K
t )
Pr (xx Mx E[x Mx 2 t) S exp
-emin
M
K
where
1| denotes the Frobenius norm, K maxigisa 0: Blexp(a? //////]]/
2), and is universal constant.
B. Proof of Lemma 2.7
In this section, we present the proof of Lemma 2.7. Recall that our objective is, for A E Raxd of rank d with n 2 d, p.s.d. C Rdxd, some random sampling matrix with number trials m and importance sampling distribution (miii= as in Definition 2.1, deff = (r() the effective dimension of A given C with Ac A(A A + C) CC-1// as Definition 2.2, and max importance sampling approximation factor Pmax as in Definition 2.3, to show that for Cem log (dett / /8)/e, ATS SAc is an (c, aappproxxinnnnnnnn of ACAc. That
(1 - 3 AS SAC 3 (1+e)A
(13)
holds with probability at least 1- 6.
First note that, for Ac A(A A A C)- Rexd have
ac;, aci,
ATAc
m = =D F.,
A SAC ATAC = Σ
MT,
m
where a.; E Rd denotes the ith row of Ac, Given Ac, follows from Definition 2.1 that the indices i, are i.i.d. drawn with replacement from the index set {1......)) .n]. We would like apply the Intrinsic matrix Bernstein, Lemma A.2, to bound the sum of random matrices J TSAC ATAc Σ. F.. To this end, note first that given Ac we have
E F = Oa.
(14)
then remains bound (i) the spectral norm |[ and (ii) the sum Em A.I and the triangle inequality, for I..... m, we get
TE [F the sense of p.s.d. matrices. By Lemma
accac,
1
|accll
Pmax P1,
F. < max ISi5D
max m
mm;;
Isisn
mi
where we recall Acll < 1 and Pmax = maxisisn |lac (T:der)) 2 1 is the max importance sampling approxima- tion factor in Definition 2.3.
Further note that, per its definition F, ec;, ac;,act
AGAC
we have
AT ATAC
|lacs, 1| aci, aci, F2 F. m2m?
ATAC
A ATAC
(ATAc)²
F
F +
(15)
m
m
m
m²
Taking expectations both sides of (15) and applying (14), we get
(A]Ac)2
|ac
1| ac
Pmax dotf m2
F
E
m²
m
m²m?
A ATACC
Thus,
Pmax deff
EF
A
A
=P,
m
14

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
for which we have
|P|| Amaxdatf
A Acll and tr(P) Pmmaedalt so that
tr [P]
dat AAA
Recall that A < I, applying Lemma A.2, we obtain
1 - Pr
ACAc I ATS SAC I (1+
1
m2
< Pr (I SAC ACACl > c) = Pr E F, >
22/2
6 <4|AAAA11 -eeft exp
Prnzy thaff
A A
+ Pmax
e/3
m
2
= =4|AAAC!| defl exp
Penaxx deft
A + +2c/3)
n7
Clearly, Pr( Pr. I ATSTSAC 1 (1+ $1
holds if
(- 4| A|A|| ddtt exp (- Pmaxddatt (2ll + 22/3)
"Pemeedatt
5..
(16)
that can be rewritten as
Pmax detf (2/ (2 + 26/3)
4 A4ATTAC | deet
m 2
log
(17)
2
5
Since e < and ||Al < 1, (17) holds if m 2
Spdar
log(
AAAAAcl datt
3c
Lastly, let check the condition 2 ||| /2 + 3//3 in Lemma A.2 is satisfied. Solving (16) for we get
2
4| A A
- dottexp
2
-mmaadall (211 A 22/3)
m
2 doff
log( (|| [A - /o)t/50cc - @Pmaxdeft deff || AT Acg(4))AAAcll1-66 /dd///8) 2
2
3m
3m
2 ZuE 6/ ATAC 10 2
2 +
2 + 611 A A
where v
Panax
log(
A|ATTAC 'dott dalf ). As 8 and A||ATAcll- > e, we have log( A||AT.Acll' d'itt
> 1. As
3m
such, get / and 6l| A Acclv |P]] so that 2 |P] /2 + P1/3 holds. This concludes the proof of Lemma 2.7.
Given the subspace emhedding result in Lemma 2.7, can be shown that the inversion (A T SA + C) also satisfies a similar subspace embedding property and close to the (population) inversion (A + C)-1, This is given in the following result.
Lemma B.1. Under the settings and notations of Lemma 2.7, we have that ST SA + C) -1 is an (c, approximmati of (AA A + C
Proof of Lemma B.I. Denote H ATA + (so that H H S 1), it then follows from Lemma 2.7 that
15

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
A is an (a:o)-approximation of A Ac that is,
1
A E ASSSAAC 3 A)
1+c
1 (ACAc + H H CH 8) A SAC H 31+ - HH- H
1
1
H A
S
SA C3 e)+
1
1
H (A C) 3 (1+
1+
where we recall the definition Ac = A A + C) 1/2, This concludes the proof Lemma B.I.
C. Proof of Proposition 2.8
The proof of Proposition 2.8 follows roughly the same line that of Dereziñski et al. (2021a, Theorem 6) and Dereziński et al. (2021b, Theorem 11). It is provided here for completeness.
Here, we focus on the differences from the proof in Dereziński et al. (2021a;b) by considering now under study is a random sampling matrix as Definition 2.1.
To with, we introduce the following condition that is erucial to the proof in Dereziński et (2021a;b) and of our Proposition 2.8 here.
Condition C.1 (Concentration property of random vector s). Given V Rxxd the n-dimensional random vector satisfies Varls VBV sj S a tr( (VBBvTT for all p.s.d. matrices B Rdxd and some > 0.
The proof of Proposition 2.8 then comes in the following two steps:
1. construct a high probability event Ç (via subspace embedding-type results Lemma 2.7), on which the inverse S SA + C)-1 exists (in particular for = 0g); and
| H] via "leave-one-out" type analysis, for 2. conditioned on that event Ç, bound the spectral norm - [
A + C and Q (>ATs) TTATSSSA SA+ C)-1 with 7 m-d
Let us start with the first step, for random sampling in Definition 2.1, denote x =eT./ VT. E Rd the ith row the sketch A, so that E E[x,x,] = AT A. Denote
- A 2 SA
X,X, +C = m
(18)
del
and
Σ E x + iike
Q-
that is independent of X
Without loss of generality, assume that = m/3 is an integer, and define the following events:
G:
Σ
x,X,
ATA 111,233 <= N G3-
(19)
1
Recall that 7= m-dei 1, the events Çj imply
X,x, = A TA j= 1,233
8= 1)+1
16

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Here, the event S presents that the (weighted) average of the rank-one matrices x,X, over the corresponding j-th third of indices 1.....7 forms sketch (of size m/3) of A A that is "lower" " / of A A sense of Definition 2.4.
A in the
In the case of random sampling in Definition 2.1, the events (1, C2, and {3 are independent. As such, for each E {I...., m). there exists an index j = (ss) E (1,2,3) such that
1. Cj independent of x3; and
2. conditioned on S one has Q z Q- = 6H which, for m > 2den, further leads to / "2,"HH v 12I
Following the decomposition of [ in Dereziński et al. (2021a;b) with R. 1 xxx Q we write
(2-
Ia - E Ec Bc[Q-.(e -
A ATA)] Ex - T A+B,
Q- X
(20)
(75 7/2
Zo
21
Recall that Q in (18) is p.s.d. symmetric, to establish the result in Proposition 2.8, suffices to bound the spectral norm
Ia 1 (QJH Z. H $ |H H
HH +11
(21)
for Zo. Z Ze defined in (20).
Without loss of generality, assume that the events <i and C2 are independent of x.. and set Ç' = ( as well as 63- = -- To bound (21), we first recall the following results from the proof of Dereziński et al. (2021a:b).
Lemma C.2. Under the settings and notations of Proposition 2.8, we have,
1. for p.s.d. random matrix M (or non-negative random variable) living in the probability space of S.
BI[TT= ;= MM
1
B
M 3 2 [M].
(22)
Pr(()
1-8
where 1,, is the indicator of the event Çz.
2. ||H 0
Next, bound the terms invoking Zo and 2, in (21), particularly by emphasizing the difference between our proof of Proposition 2.8 here from that of Dereziński et al. (2021a;b). By 83 < and
2o 1 83
E
(x,x
A ATA
1.h.
it follows that
H Z, + 12E
[1 + H 1-1x,x) x.)) 1-cal-
We now need to bound the (conditional) expectation Eç.[(1 + x H x.) 1-(s]- To do this, we resort to Condition C.I with an appropriate choice of a. Note that
E| H
H x.] tr(HH-AA A) = deff.
per Definition 2.2, it follows that
Var(x] H x.] Varle, //////TT AH A e,, VT,
TTe,
(a. H (1--aa,) )2
= [x H H xx,) 21 [E[x] H xx,))
deff
1=1
Smaxx La (pmaax
17

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN for max importance sampling approximation factor Pmax in Definition 2.3, which implies that or in Condition C.I satisfies
a (Pmaxdett)
Further, by Chebyshev's inequality, one has, for 2 2detf, that Pr(xg H 2 x 53) S cdoot/ xx". This further leads to, for 83 S 1/m².
Zo $ S 12 / Pr(x, H 'xx, 1-co 2 x (')dx + 0
1
m
1
2 24m 8g + 12
Pr(x. H H 2 x (')dx + 0
2m2
m3
adeff m2
a'V doeff
0
m
We then move on to bound the last term concerning 2, in (21). Applying Cauchy-Schwarz inequality twice, we get
H +1 = V Bc 1(5
sup ||||=
E
(u
sup ||||==1
Ex
(u H $x,)
(23)
Ti
T2
T3
We start by bounding the term T3 using Condition C.I with B
AH u/ H 1/2 AT Noting that tr(BB
H-1/2 A AA--1/2 1/2u S 1, we get
Eg (u T H 1x, xx
s 2 (u H x,) = 2T [ H uu "H ixx,) )
Var [xx]] H uu H xx- + [2]
H 3 uu H xx.))
(a. H Lu T H aa) 2(rr(
2
T4
2p deff
H a H uu 2(tr
< 2pmaxede tr BB + 2(tr
- 2(pdet + 1) 0(a -),
for a O(pmaxdett) in Condition C.I. This results in 7; = = I). Similarly, we bound T2 by taking B =
AQ
AQ_.H!/² 1/22
H T in Condition C.I. It follows tr(B) S T 2H H2/2)/u S 62 and tr(B²) S 64, so
that T₃- O(Va+I).
It thus remains to bound the first term Ti in (23). Noting 7 =E 17 = 1 z tr(B 12 -A A), we write
2₂2
Ec - ? < 2(7 - 7) +
m² E. (
E 1 A A)
272
E ((tr(AQQ A) xx (..x.)9
m
Analogously as above, letting B = AQ- T so that tr
< 36deff in Condition C.I, we have, for m 2 2d that
Ey m²
272
(IT(AQ_.A)
x Q Q-
27
2n²
Ee
Σ
a. aTQ .a )
m
T
i==1
2x2 m Be [(T[AQ AT)))
+
2x2
Be
(a.
Q Q-
422
E
[tt(AQQ.AA) 1)
m
T
m
72-? Pinax m²
de
272
Be (I((AQ.AT))
m²
72-2 Pinaxe Pott
72,2Punxx
288pmaxiles
0 ( (addetf m²
<
m² 2
m
18

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
with again O ( in Condition C.I.
This, following the line of arguments in Dereziński et al. (2021a), further leads to
ad
2 [(% 7) $ 0
m²
which, together with 1 -1 71 O(v (0(Vadd//m) yields that Tj (vvdde mm)
This allows to conclude that
(avdeft = m
(
H ||| T3 Ts
Putting everything together, we conclude that
[ H H& -0
m
with (efaxxeet) in Condition C.I. This concludes the proof of Proposition 2.8.
D. Proof and Discussions of Theorem 3.1
In this section, we start by presenting in Appendix D.1 the intuition for the self-consistent equation (4) Theorem 3.1. The detailed proof of Theorem 3.1 is given in Appendix D.2. In Appendix D.3, we provide discussions and auxiliary results on Theorem 3.
D.1. RMT Intuition on the Self-consistent Equation in Theorem 3.1
Here, we present a heuristic derivation of the ceff-consistent equation in (4) Theorem 3.1, Let recall some notations from Theorem 3.1. Let x - /Ti A as in Appendix C. For the ease of further use, denote
_ x + +C) ))²
Q (A
and Q- =0x==Ejta 1 x + CC- for which we get
Σ 8=1
E m
Σ a.a = ATA.
Then, follow the deterministic equivalent framework (see Couillet & Liao (2022, Chapter 2) for introduction) and
show that [ H 1 = 0, for H A DA C, where D E uu is the diagonal matrix in Theorem 3.1. Note that
Q H
E Q D E TIQATS SAJH 20
and using Sherman- Sherman-Morrison formula in Lemma A.3, we further ascertain that
X.
a-aa, H a. Q-
E[QA STATTSAARE
'1=
E
E
Q
--.x//m
Q
a;a
H-1
E
a
H
1=1
Using the rank-one perturbation lemma in Silverstein & Bai (1995, Lemma 2.6), we obtain
SS
Qa
H
E
a,a
1+a H
aa/m/
a
a mm,
E Q A DAH-1
so that D diag
This leads to the self-consistent equation in (4) of Theorem 3.
19

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
D.2. Detailed Proof of Theorem 3.1
As outlined in Appendix C, the proof of Theorem 3.1 also comes in the following two steps:
1. construct high probability event in (19); and
2. conditioned on that event ç, bound the spectral norm ||_ |Eç[Q using "leave-one-out" analysis.
Furthermore, for each {l,...,m), m., there also exists an index j - j(s)) {1,2,3) such that, conditioned Cj. we have
Q Q- 1 6H
To complete the proof of Theorem 3.1, we first rewrite
[La H [Q] H (Eç Q Ia) H
Taking Yss = + 1 x Q Xx, Q-, 8
1, ,,...,m.. m, and D.
1+ x, H-1x," we then get
R2 Q Ia = (Rc(Q) H H Egg Q
SASTAA H
= 1Q (A DA - S SA)]
1
= =E Q- A
DA E -
DA X,x,x, +B
Xx
7s
21
Z2
7g
which yields
[L - - E
-2 [ Zz)
||
H - H
(24)
Then, we bound the first term H 2 H 1/ Together with the fact that the event Ç is independent of x,, and using the Sherman-Moriison formula in Lemma A.3 and (22), we get
Ee Qc - Q] = Eel m77
T Q-al Q-al=
Q-xx,x
m
2 R
2
1Q_xx,x T Q-a)
xx,
2
c IQ
a,a [aa]Q=.]] Q-1]
2 E. |Q_
2 E. [Q,
AQ AQ_
m
m
which, by incorporating the fact that H - H and the event indicates H /2 Q 1/2 = 61,, leads
H Z
H 3
H
. A
DAH b11
Q Q H
H H
H 1 A DAH :|
QQ-
H
I m
E H .HH.Q_HHH.
72
(25)
m
Further, we move on to bound the second term in (24). Recalling the assumption regarding <1, (2, C2. 5' and 63 from Appendix C, we have
E 10- (A DA D.x.x. =
E.. (E.r Q (A DA - - DA D.x.x.) 1-col)
83
Ec 1Q- DA - D.x.x.) -1-cll
83
20

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Then, follows that
H || 2 - x 1 1 $1
< 2 H 3 (A DA x,x,x]) 1 H
s 12 H
H (A DA - D,x,x
H
$ 12E 12E A DA - D,x,x,] ) 11 } |||
< <2
[(1 + H D,x,x!1 H $ 11
g [(1 D,x, H x,) 1 [- 125y 12E (D, H 'xxx. 1-sa]-
follows from Lemma D.3 and m 2 2pmaxdelf that
1
m + Pemaxdelff
H
H
m
2Pmax deff
S 2,
(26)
m
m
so that
2
ma H
Vare (D., xx) s E2. [DD,( H x,)
π
mny + a, H 1 a,
j=1
D2
Day (a) 2ay)²
H 12 H ag) 2
Jj
< 4 Apmax
Σ
a T H 3T''a
aj S 4pmax
==1
and
Be 1Dx H xx]]
Mm
H
a, 5 1=1
H
H [1,a
H 1 ay S 2deff.
j=1 m + a H aa aj
Using Chebyshev's inequality, we get, for x > 2detf.
4pmax deef
Pr(D.x, H 2 (0)
x22
This, together with 83 d < m3
and > Pmaxx dest further leads to
1
H Z
0
+ 12
Pr(DDx,HH H
''', x, 1-co 2 x (')
m3
< 0
1
24m 2dz + 12 j
Pr(D,x] T x. 2 x d'
m3
< 0
1
24
48 18
&
_d do s0(=)
24 24p La
3
m
2 x2
m
m2
=0
2 ds (V
3
21

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Using (22) and (26), bound the last term in (24) follows:
H
ZgH
Zz
sup
E
[1D
7
H X H
vv]
10/=11/|v]
21 [10 7 " H Q x, H [v]
sup
|| 1-
Q__QQ- N m m
H H
sup ||||-|||v~=
21E
D3 7
Q a a, H aj
(u X, Hin H i H Hd)
sup
R
|||==|||v|==
m
7 a Q aj a H aj
sup
E
11 H Q ,,a Q H.u
||
j
Q a a H
H
+ sup Ex¹ |v||=1
H à aaa H iv
m
Q aj a
< sup Egy ||u|=1
Σ H Q-
1515
mmj
a ajj a H Laj
+ sup E |v||=|
max 15j5m
H a H
M
j=1
Q , a H laj
< sup ||u|=1
Eg
A'H AQ ii
J
a
aj a aj
+ sup
Eyy
max 15J5n
A AH 'V
|v||=1
< Ee'
max
a
a H
aj
Q . AH
15/šn
m
+ Exy
max
a
aj
a H aj
A AH I|1
15J5n
mmy
E (a) Q- a, H 1a,1),)
(0) max 38E,
a Q s aj
a, H
aj
(M)
max 38 15j5n
15J5R
mTy
(mTT))²
Ey (a] Q a-
aj
E. [a] a212
E [a] Q_ a,l ]
a. H
Si'ayl
max 38
max 38 13j5n
(27)
18j5n
(mT,)² )2
mR;
M
M2
where in (a) we use the fact that A A = H. ~ H and the event Ç indicates H 1/2 - 6I, along with (26), and in (b) we use Cauchy-Schwall inequality. Then, by Lemma D.I, we bound the term Mj in (27) as
Pmax (a. H aa
(a. H H 11-1a,)2 (MRT) )2
(
(V
ds
Mi < max 2736L; 15j5
/ m (m
max 645L 15
0
m²
m
22

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Now, it remains to bound the second term M2 in (27). Using again (26), we get
2 a -aaa
- 2 H
a
H
H E (Q H
اayl
<a H (
H
< 2a H
-
a, H (B [--2]
-
< 2a, H aa'all
H (( - 12 |2 Q
Q
+ Ex Q -
x 2a, H''aj aj
H (Bp - B H
H
B
Q- Q
H
H
(28)
G1
G2
Noting 83 < 1 m3 and
1
EC Q
(Rg- QQ -]
1 83
we have, for the term 11 defined in (28), that
63 H
G1 1 1 83
1Q-1 1Q
[Q-1]
B.
H
24 m
20 H Q
H + H [Q -5
H
S
248s
For the term G2 in (28), it follows from (25) that G2 O 1
We thus have
M2 S
max
a
JH''a
+
Z
H
H
22H
2
| 2gH 21)
15J5n
m3
m
76 dlf
1
+
+
|
Z
H
||
||
H
H
m
m3
m
dff
1
(v
Pinas coff mi
+ +0 (V
Pinux det m3
+
+
+M2
m
m3
which, for m > 76 deft, yields that
M2 = o
(V-m
Pinax de m3
d2 (m 76pmnxxdett
762
= -VV
desf m3
max
Putting everything together, we conclude that H 2/ 0 (V Pinx of Theorem 3.1.
d /
This concludes the proof
D.3. Discussions on Theorem 3.1 and Auxiliary Results
Lemma D.I. Under the settings and notations Theorem 3.1 and let Var [1] denote the variance conditioned the event ç with x, independent of Ç. Then, we have
5184L2pmax (a) H H2 288Lz(a H a,)
Var [a Q 2-
Proof of Lemma D.I. Use Q to denote the matrix (ATST S where S at is the matrix S without the tt and Tt rows; let that, for each pair s, I, one of (1, (a is independent of both x, and x;- Without loss of generality, assume that is Si. Also, let [[]] be the expectation conditioned on and the o-field Fi generated by the basic events X Xi- To apply the Burkholder inequality Lemma A.4, we rewrite a. Q- E [a] Q Q_sa,) as the martingale
23

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
difference sequence
a Q- - _ [a a-] = Ee, m mlaj aalaa_Q_00 QQQQ_aay - Be, Qaa])
Σ B E i- Q- - aj Q-s/a/]++
+ 1=1
Σ
[ Ee (--a-aaaa Q-sayy
m
Ot + Pi),
where di = (Be' Be' (-)) at 3 a- and ip: = - -Bgg' - (e(,-))aa Q-11ay]. Analogous to (22), we have the following result, the proof of which is straightforward.
Lemma D.2. Under the settings and notations of Theorem 3.1, we have, for a p.s.d. random matrix M (or a non-negative random variable) living in the probability space of S.
EIII j=1 M Pr((")
1 E[1, M] 3 2B,, [M)
Eg M
3 1 -
(29)
where 1e is the indicator of the event Cj-
Conditioned on Ç', using (29) and the fact that
1 1aj0 Q x/ at
a Q .a a, a-
m 1 m xi Q-
we obtain the following bound on the second moment of Φr as
(
[ Laj Q x Q Q-slaj
a Q ,/a,a,
staj
E [o]] - 2Ec, [of] S 2Bc,
< 21
n
m
1 m x
Q
mmT; 1 mes. a Q sa
< 2E Σ m²r:
(a] a Q sta)) )
a) aJQ Q st HQ -/aaa H''a.
<
Bxi
T Q a Q a
?
2 deff m²
Bei
Σ
a T Q aaa Q-aaaaa Q HQ aa
1=1
2eft
Eç, [a] Q AAQ a, Q HH-HH
m²
2Pmax
Be, [a]] Q
m²
This, combined with the fact that H 12 i H/ /2 1 61,, results in
2592p Pmaxx deif
E [0]]
(a H HH'HHaa 2,)2
m
Subsequently, we proceed to bound lpx... Considering that &a is independent of x,, it follows that E,,,(a, -//] [a]] a/ay) Then, noting that, for &z Pr(-62) m
1
82 E E laj Q Q-s1aj/--44
Ey [aa] Qus1a,)
Ex Eq 12, --a
1- 82
1 82
24

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
and
82
R (21111aa 1 Q-s/aj]=
E (1
ja
Q-slaj]
Rx 11.--11 1- (a] Q-
-t -(],
1 82
82
so that
82
|0:1 = (( 1=) ---aaalll
(B 1 -/aj |
1- 62
12a H-
S Rs [ Q s/ / S 12 H aj
m
Consequently, let y1 - + pi) such that
1 V1 = a
Q- - E. Q- it follows from the Burkholder
inequality in Lemma A.4 (for T = 2) that
Varg |a] eargla, Q- = E
Σ
<
L2
Ey ((b1 + (0) - 2L
E [ + 2L
B [+f1
1=1
1=1
1=1
1=1
5184 5184T.2fmaxdeet
H 2a)²
28 H a
m
m
This concludes the proof Lemma D.I.
Below the proof the result Footnote 2.
Lemma D.3 (On the self-consise consisenstt D). For a given matrix E Ruxd S be random sampling matrix with number of trials m and importance sampling distribution )n in Desinition 2.1, and let C IRdxd be p.s.d. matrix and delf Ei=1 e with leverage score E in Definition 2.2. We have
In 3
Ins
m 2pmax deff
m + Pefi
in the sense of p.s.d. matrices, where we recall Pmax = maxigign e. (T;deet and Puin minigign / (T:dee))
Proof of Lemma D.3. For each Dis i = I,...,n, n. it follows from its definition in Theorem 3.1 that
mT;
mm;
D&
mJ; + a. (A DA + Ca
HMIT + D min a (ATA AATA+ CC-
where Dmin minj Kisn Dj < 1. Without loss of generality, let Dmin Dun and use the fact that m 2 Cpmaxx 2 CCC / = Ca, (A A C) an/ Thi we obtain lower bound on Dan as
mr
Min
Dun
mon a (A DA +
-Ia
mmn D a A an
-la
MRR
C
D mon
C+D
C+
so that Dan
A> 1/2, for some C <>} 2. Then, it follows that, for i 1, 1,...;
mT;;
m
m
D N
>
mm: + A l (ATA A + Cla m+A-la (A
C)-la,/m
m + A Pmaxxelll
This together with C' > 2 results in
m
DR 2 + 2pmax
On the other hand, we have
D4 m a (A DA + C)-1a//T, Consequently, we have that m m+20masdet In DD
17 + a (A AATA
+ Pmii desf
anttnanndat In this concludes the proof of Lemma D.3.
25

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
E. Proof of Proposition 3.2
In this section, we start by presenting the RMT intuition for Proposition 3.2, and particular how the exact leverage scores e comes into play the analysis and debiasing, Appendix E.I. The detailed proof of Proposition 3.2 is then given in Appendix E.2; and finally we provide discussions and auxiliary results Proposition 3.2 in Appendix E.3.
E. RMT Intuition on Proposition 3.2
Here, we present the heuristic derivation of Proposition 3.2. First, let us recall some notations from Proposition 3.2: let xg e VT. A as Appendix for the ease of further use, let
C; C mm
S
Pasts
s=1
for some deterministic Fj to be specified, and
Q = (AATSSSAACC)
Σ Fi,, x,X, +
and similarly Q- - ( E140 XIX +C) for which we have
X P.aaa.
Our objective here is find Q (and Fj²) such that |E H H 0, H (ATA A C)
This is
E E -
E [[Q]A AH EIQATS = 0.
Using the Sherman-Morrison formula in Lemma A.3, we further ascertain that
Q- A Ae, TT,TT,,HH"" 1 P. c, AQ A ATc, /mmT.
n
F.A. Te,e, T AH
EQQA SIQAI " =E
E
1 Fac AQ-,AT AT://m;;
A e AH-I
E 1+ e AH e, /mT;
i=1
where we observe that the exact leverage score e, AH AHH'ATO naturally appears in the denominator from this leave-one-out derivation.
AATe; e; a (AT A A+C) -aa EC given C as in Definition 2.2
This, together with the rank-one perturbation lemma in Silverstein & Bai (1995, Lemma 2.6), gives that
eRe.e + F.ce EC //m
EIÂA BRATSTSAH'SS =JQQA
AH EIQ]ATAH".
where we take the debiasing factor Fi; mTTs///// (mT+ ( such that Pa/(1 Puel ///..//mm))))m) 1. This thus leads
$ = diag { /mm(mm-e//mm))) m
S, E (1...,n))
E.2. Detailed Proof of Proposition 3.2
The proof of Proposition 3.2 follows the same line of arguments as in that of Theorem 3.1. In the following, presented the proof for completeness. We specifically highlight the difference in the proof of Proposition 3.2 from that of Theorem 3.1, where Q is replaced by Q.
The proof of Proposition 3.2 also comes in the following two steps:
26

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
1. construct high-probability event ç, in (19); and
2. conditional the event derive bound for the spectral norm [ [] analysis.
using again the "leave-one-out" type
For the second term above, we have, for m > 2pmax dest that
mi,,
Fis. < Moii, 2 -miii
2: 8 I,...,m.
(30)
This yields that, for given (1,...,m,, m.,, there exists index j jj(s) (1,2,3) such that Q 3 Q- 6H holds on the event Çz-
Denote %/2 = 1 1 m F,
x, Q- X,,, we then obtain
m A
Eg QQHH- Ia (Eg) Q
H
H
Q
F, X
H-1 H
m
m A A
Eg Q
F
EEIQAA A - -Bg
Q
F
m
9
Q A
xIx +B₃ (1-
Figio
Q Q_
9/3
22
Za
which leads to
- H H
| (Eç H H- & |H 1|
2
H 211
(31)
We first bound the first term |H / "2,H-112 in (31). Adapting the bound for 12 and H 1/ 1 in the proof of Theorem 3.1 in Appendix D.2, and using (30) along with the fact that 1/2 Q2, 61 when conditioned on
< we get
Fist 7am
Ex[Q_ - Q <
Ec
Q __ 3
(32)
and
H 2
AT 211 | | Q- H. A AH $1 =0
1
m
Next, we bound the second term
|/ H 1/2 in (31). Note that E [x,
H x,] = delf and
H a,
)2
(a) H aay)
Vare [x_ H
xx,
< Eg [x
xx,xx |
T
S
Pmax &a..
Then, using Chebyshev's inequality, we have, for 2d that
Pmax &
Pr(x, H x, 2 x 1 5') <
x
Analogous to the bound on LE< [Ex[Q- (A DA D.x,x.11 )] and H 2 1/2 in Appendix D.2, we get
DAA
E (A T A x.X, =
Ex (-. (A TATA xxx!) 1-sgl:
1 8g
27

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
and
| 22H 211 120g c [x, H H 1-(3] 12dg 12 Pr(x, H Hx, <a 2 (')da
(
(v Pinne ell m
= 0
Now, we move on bound the last term
HH/2 2 in (31). Considering (30), it follows that
Fi,
a a C
a Q 2--aa, as,
c
a
aa C
F
7s
m
CC + a Q-
mR.
CC
miT
i Q Q- pC
a
a a,
ai.
< 2
(33)
mA;;
Following the bound on || 11/2|| in Appendix D.2, and recalling (30), we further have
H 2g
Fi, %
u H x
X H iv
sup
|||||=1|||||
E (a]] Q a a, HH-H-ay aa) )2
< max 74 15/5n
(MTT)²
Egg [a]Q Q aj
150 [a Q_
E [a] @yaaalll a. 1SSHap ag
max 74
max 74
(34)
15
(mT,) )2
1855
M2
M2
Subsequently, like the bounds for M1 and M2 in Appendix D.2, and using (30), we obtain the following bounds for M1 and M2:
Mi =0 0 (( Pinaz
m
and
74a, 'a a
M2 < max 18i5n
H (E. E=) | HH QJ
H (Eş[QQ --)
miy
= =0
d3 (V Piiax
3
-0(Viaex Pf / and thus complete the proof of Combining the above, we conclude that [ H2/2E, Q2 Proposition 3.2.
=0 -00
E.3. Discussions on Proposition 3.2 and Auxiliary Results
Here, we present some auxiliary results in addition to the fine-grained debiasing results in Proposition 3.2. Precisely, we show in Appendix E. that by substituting exact leverage scores with some (computationally efficient) approximations leads to a controlled inversion bias. We then provide the proof of Corollary 3.4 in Appendix E.3.2. A "counter-example" to Dereziński et al. (2021b, Theorem 10) that uses sealar debiasing to achieve zero inversion bias in the case of exact leverage score sampling is given in Appendix E.3.3. Finally, we establish the proof of Corollary 3.7 in Appendix E.3.4.
28

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
E.3.1. DISCUSSIONS DEBIASING $ USING APPROXIMATE LEVERAGE SCORES
Corollary E.1 (Fine-grained debiasing $ using approximate leverage scores), Under the setting and notations of Theorem 3.1, if one approximate leverage scores é. in the debiasing matrix, instead of the exact leverage scores eg in (5) in Proposition 3.2, that is
$ diag V mm/mm-m////
/. (../mT.))) S, (1-w), E s (L+W)) w [0101
Then, there exists C > O independent n, deft such that for m 2 Cfmaxdett (log /6) + max fw/e, with 8 S m (A SA C)-1 is an (e, (enbiaud estimator of C)
Proof of Corollary E.I. The proof of Corollary E.1 largely follows the proof of Proposition 3.2, and is presented here for completeness. In the following, we emphasize the difference from the proof of Proposition 3.2, particularly regarding the matrix S, which is constructed using the approximate leverage scores denoted by E as
-m/(mm-9./mm $ diag { -0an---g//1 (.//m;) s=1
m
S, is nl.
This implies, in the proof of Proposition 3.2 that
MA
Fus mT. CC
(35)
Note that in this setting, we have, for w E [0, that pep e. I wi? and &p s (1 w)ef s 20 As such, for del CC///T, with C 4, we have
mA, 2
Fun S mn,,
<
Thus, adapting the proof of Proposition 3.2, we get the same bounds on | H 1/21| and |H 1/2 as in Appendix E.2, that
H H | = O H Z, H 3 -
1
(V
(Vinna m
eff
Next, we establish a bound for the term
| 12 1/2
in the setting of Corollary E.1. Using (35), we obtain
F
a Q- cC
a Q- a,
a
Q-
TC
<
<
C + a a
(1 C
mT,
Q ai
e is
CC
a; Q
Q_ (C
2 2wpm
< 2
+ 2
< 2
mi;
MT,
This, together with (34) yields that
2 =
WPmaax
2 d3 m3
ell
m
Putting these bounds together, we conclude the proof of Corollary E.I.
E.3.2. PROOF OF COROLLARY 3.4
The proof Corollary 3.4 largely mirrors that Proposition 3.2, and is included here for thoroughness. For m > Cdeat with C > we have
mm; (1 C- mmn,
Fi. m dest mR;
2.
Tia
29

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Then, recalling ny E [ Pmox)) Pmax)) ], and (33) can be rewritten
( ((iit) such that |. -///deet S spli/den with cp mmax ( 1), (1
Fi
a aa di
a Q aa,
deft
a
cC deft As,
1
+2
1s
m - Ti, + +a, Qa
(1 C -
MTT,
mTs,
a ,aa -6e
L./delf The
2 a, 5āi
/e/dea
2detf
22
miii
mmii,
mi
a ,aa,,
2epPmax
< 2
m
This, combined with the bound in (34), gives that
2 =0
CoPmuxdett
P da m33
+
m
Recall the proof of Proposition 3.2 that
H #1= (1)
(
(iinneaat m3
H 2 =1
m
Putting these together, we conclude the proof of Corollary 3.4.
In particular, taking Piax 3/2 and Pmin 1/2, have cp and thus inversion bias of order (dett /m).
E.3.3. A "COUNTEREXAMPLE" TO DEREZIŇSKI ET AL. (2021B, THEOREM 10)
Corollary E.2 ("Counterexample" to the lower bound in Dereziński et al. (2021b, Theorem 10) using exact leverage sampling). For any n 2 2d 2 4, there exists E Read an exact leverage score sampling matrix S E Rmxn in Definition 2.2 and a high probability event (that ensures the invertibility of A " SA), such that when conditioned on Ç. SA + C)-1 an unbiased estimator (A A C) for = m B [1/b, b distributed as Binomial(m, 1/d). and C 0g.
(ATA
Proof of Corollary E.2. We begin by recalling the setup of the matrix A in Dereziíski al. (2021b, Theorem 10). Without loss of generality, assume 2d (otherwise, pad A with zeros). The matrix A consists of 2d scaled standard of for i 2 2, and the first two rows are /
basis vectors, where consecutive rows are defined as a (2(111)+1 a 22(1-1)42
a
L a / eTT This is
v
0
0
000
0
0
0
0
0
In this case we have A A = La and 02, so that the leverage scores satisfy & 4.
c 3. and CC
2(1-1)+1
cC "2(i-1)+2 2 for i 2 2. The leverage score sampling distribution {mi) in Definition 2.2, used to construct the sampling matrix S of size m 2 d, is thus given by
for i = 1,
4d
3
for i 2,
1d
otherwise.
27
30

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
For any 7 > O, the matrix 7A S 'A is diagonal, with diagonal entries given by
[ATST i = yyb, for ,,
m
where by's 's are all identically (but not independently) distributed as Binomial (m, 1/d). Then, conditioning on any high probability event Ç that ensures the invertibility A T
TSSSAA SA, it follows that
E TS SA)
Ex
for i 1,...,d,
yd
bi
so that
Ex (7A | = (ATA)
with = E [1/b], where is any one of the aforementioned (b.) This completes the proof of Corollary E.2.
E.3.4. PROOF OF COROLLARY 3.7
The proof of Corollary 3.7 comes the following two steps:
1. construct two independent high probability events: Sw which the randomized of A as in Definition 3.6 has approximately uniform leverages scores; and Ç as in (19) to subspace embedding, but for x H,,,,, VIOTs E IRd the ith row of ASRIT: and
transform H,, A
2. conditioned on the event SN apply the bounds in Appendix E.3.2 to obtain inversion bias for is the effective dimension of A given
dett ASHHTAARRE + +C) where dert = Ei= l
We start by recalling some notations from (the proof of) Proposition 3.2. Denote Ac = AH-1/2 with H - A A + C and let e. (H, = |ee A,ClA be the leverage score of H, A/ / /n. Define the following event tw:
Sw : e (HDD doff
S w.
AA vi n
n
Note that this is equivalent to, for uniform sampling with ni 1/n that
e (H,,,, A/VA) down;
W visisn
mm;
m
Next, we verify that the event fu holds with a controlled probability. Precisely, some C > 0, it follows from Lemma E.4 below and |Acl| = dea that
deff Pr (H,DAA//V))
S 5. 2 2w,1sisn
(36)
n
with w 2 max{C den log(n/ /0)/n, C'lo((n/8)))) which is equivalent to
(HH de
Cod log(n/8) Cllo(n/8))
max
for 1 S S n.
MT:;
Further recall from Lacotte et al. (2021, Lemma 1) that sampling size m 2 Cpmax (deff + looool(// (og()/(Ed)) log (og(def/ e)),
Pmax = maxigien ne (H,, D, /ndddet and E (0, 1/2] suffices to ensure that H 1/2 1//AASRHT A H 1/2 is an
(c, apprrxxaatiiiooo of A T c. This leads Pr(() < 8/2. Combining the above results, we obtain Pr(( (w)) S 6. This. together with
(V
Pinax m
13 off
C CHH doff
H 2gH
+ +0
max Isisn
Mi:
H 2 H -
0
1
|
=
(v
Pinnx ds off m3
m
in Appendix E.3.2, concludes the proof of Corollary 3.7.
31

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Remark E.3 (SRHT inversion bias with different sampling sizes). Recalling (36), we observe that
max ne, (H,D;AA/M)/der 2 1s
Pmax
holds with the probability at least 1 1- 6/2. Consequently, choosing m 2 C(dett + log(1/(c6)) log(den/ (6)/// suffices to ensure that H-1/2 A AltArr-ASHTT 1/2 an (e, appppooxxinnnooonn (or subspace embedding) of AcAc. On the other hand, note from Corollary 3.7 that
1. by taking
0, a number of m (dder) samples with n exp( < 6 m-3 we have that m- dots A ASHHT + C)-1 is an (O(1), -)-uniiseed estimator of (A (ATA) A C)-l: This, in particular, agrees with the above subspace embedding condition up logarithmic factor.
2. taking 1 in Corollary 3.7 allows one to (further) reduce the inversion bias (to level that is significantly smaller than ((1) by increasing the sample size m (ddatt")
Lemma E.4 (Row norms). Let Hn E TRux be the Walsh- Hadamadd matrix of size 2 4 as in Definition 3.6 and D diag(v) E Ruxx with v E IR" a Rademacher random vector. Then, we have, for a matrix X E BRaxd with n 2 d and 2 max [|xx XX | (ong(22))/(cr2)) (c)(()))),)))))))))))))))))))))) ||2 log(2n/0) 8) (cn)) that
c H D
X
Pr
2 2tisisn Sb.
where ||p denotes the Frobenius norm, and > 0 is universal constant.
Proof of Lemma E.4. Recall from Definition 3.6 that both H,, and D,, are orthogonal matrices such that H D H /n = In-
Fix a row index i {I......0 n. and consider
D2 = ||0 EX|
π
where E di H /n) a diagonal matrix formed from the 8th row of H /n. Observe that E I we further have
X
E 1
EXX
Euj = (E Eve EX]) = tr E
11
Note that
X|
EXX E X X
n
XX 1
EXX = tr((EEEXX E E) tr E E²)
n²
and
K = max inf > 0: E. siiaa
E[exp(ef (uf/s²)] 2 (log(2)) 1.
where is the ith variable of v. max(log(2)]
Apply Hanson- Wright inequality in Theorem A.5 with t 2 ( log(2n/5)) (0/55)/((n))) (cn)), for each i = 1..... we have
Hanson-Wright
(eex/ol/llg221 (on//oo((2n)))) (2n)(nn))/((nn)) [1] X
e
H
X
8
Pr
Taking a union bound over these n events, we conclude the proof of Lemma E.4.
32

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
F. Proof of the Results in Section 4
In this section, we present proof the results Section Precisely, to establish convergence guarantees for the bias- corrected SSN iteration in (9), we first introduce in Proposition F.1 fine-graing non-asymptotic bound on second moment of the normalized sub-sampled inverse matrix in Proposition F.1 and present its proof in Appendix F.I. Stemming from this second inverse moment analysis, also provide non-asymptotic bounds the second inverse moment for the exact and/or approximate leverage sampling and SRHT in Corollary F.2 and Corollary F.3, respectively, under scalar debiasing. Then, in Appendix F.2, we provide the detailed proof of Theorem 4.3. Finally, further discussions and additional results related to Theorem 4.3 are provided in Appendix F.3.
Proposition F.1 (Fine-graiie cene-grained analysis second inverse moment). For a given matrix E Raxd let be a random sampling matrix with number of trials m and importance sampling distribution { [ as Definition 2. I, and let C E Raxd be a p.s.d. matrix and CA (ATA+ A C)-1/²C(A A C) [g Define doff tr (A tr(AAcc) with m
Ac A+C) 1/2 E gaxd and the fine- era-igr de- ee-aiased sampling matrix $ diag { mm/mm - m./mT,)) ==1
.
as in Proposition 3.2. Then, for diagonal matrix F = diag ( with
F. ac; [A SAT + CA) -lac
mm.
and ac; T E Rd the 7th row of Ac, there exists universal constant C > independent of doff such that for m 2 Camaxden (log((d (dea/)) /0) + +/22/3) with 8 S m and max factor Pmax max e (Tidda) in Definition 2.3, (A SSTTAAA++ + CA) is (e, 6) -ubbaased estimator of Ia Ac FAc-
Note that Proposition F.1 can be seen as second-order extension of (the first-order inversion bias in) Proposition 3.2. As we shall see below in Appendix F.2, this result instrumental establishing the convergence rate for SSN in Theorem 4.3.
F.. Proof of Proposition F.1
Following the methodology of the proofs of Theorem 3.1 and Proposition 3.2, the proof of Proposition F.1 also comes in the following two steps:
1. construct an high probability event Ç as (19); and
2. conditional on the event Ç. derive a bound for the spectral norm [| SAc+ using again "leave-one-ou" type analysis.
+ CA) (La AFF)))
First, let recall some notations from the proofs of Theorem 3.1 and Proposition 3.2. For the ease of further use, denote
Š diag (Fr 1 S, is E [l....,n. n). Pss mm/mm (g,/TTs);
= A AH 1/2 1/2 and x = =e,, Ac/ /Tt such that E[8,8]1 ATAc Further let
& + CA c₂)" -1
m F 1
m
+CA +CA )
Q (ATS TATSTAAC
SAC + -N Σ m
and Q-,
Σ
148
To prove Proposition F.1, first rewrite
Ex10* - (I. AFFC) = E. I.]
+ [a) AFFAC
(37)
T.
33

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Then, along with A AC + CA Ia and let 7.- -1+ m
F Q- ts. 1,- ..... mi, get
F L. Q- 8
E<[Q(Q- [a]] 2 QQQQIAAAAAA
SAc))
E₂ ATAd - EE
F L.,,,
F
-EcIQ'ATAc| - E
+ E
Sade
m²s
m
L==1
F Q &
F2
Q
X Q² X
=EGIOZATA - E3
+ IEg
As
m
= Q Q2.8.8
E&
Q
1
Pigio
8x8
+Br
F2 Q 28.8 _x,xx Lade
(38)
%s
m
Ts
where the second and third equalities follow from the Sherman- Morrison formula. Further define F F ac; [ //acc//mt m/m; Using (37) and (38), we get
=diag F F with
F2 Q- x Q
X
E<[ (La + ACFAc) T1 + T2 + +3+ T3+ 133
AFFAC
m
-T++T.
T3 1Q A|
Ecl(Q Q)A[FAc
T.
T.
Ec( (AFF Ac A[FAC)]
To
X Q &, Q X
+ Ex
EA( Ac))
F2
x5 Q2 x. X
+ E2
x2
Ts
At this point, note from Proposition 3.2 that || = 0 (V Pi &ss mm3 that Ç' implies that
Using Q -1/2QHH// and noting the fact
Q 6H
(39)
we have Fie S ac. (E<[Q"]acc acy/mm;; S 36pmaxdet/mm < 1, which, together with Proposition 3.2 yields that ||aa|| =
p Paanx & / mm3
Now, we proceed to bound T2. We first write
T22 = |E A - 12 A + [ - R 1E[
[[ Q" ATAAJI + | (ATAc
To
Txa
Recalling A A = H A AH 1/2 = I. we get, for the first term T21. that
Tai <
[2
22.11 < [ Q --10I + [ - 2-.1)-
(a)
Fi, QQ 2 Q_
Fi Q 2. Q² ||
Bc
E.
m²s
m/s
(b)
2
Fi.io Q &2,x,
LE
p2 P
2 Q 2 Q-)
(40)
mns
m273
34

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
where used (again) Sherman-Morrison formula twice in the step (a) and (6), respectively. For the first term in (40), using the fact that Fi,i,/9 Fi,i, 2, together with (22) and (39), we get
Fuc m²
F m
B5
-_ Q-
sup
B
QTT x, Q-
& x vv
sup
m ||||| |=||||||||||=
E
u
Q x
u u Q Q-
|
||||=||||||=
u Q ac ac, Q Q. aa,aac Q-sv
sup
m
||u|=-1||v|==1
1=1
2
u Q ACAc uuv Q
sup
||
o
Analogously as above, we have, for the second term in (40) that
F
Q
8
8
Q
x
Q
E
IQ
x
Q-
8 Q-
m
m
Q- ac, _Q ac, _Q_
E
m
Jj
288p
E. [
=
(fmaxdefft = m2
m²
We thus conclude that T21 = O(1/m).
Following the methodology used to bound 2 and |||//2 -1/2|| in Appendix F.2, together with (39), we
get similarly
T22 = -0
(
/inax Pimaxxx off m
(V
/ m3
ds
[T3l = O
Utilizing the fact that Fa S 36pmax /m. together with (32) and (39), it follows that
Prmaxtell m
|T/| S Q A So
Since (39) and T21 = O(1mm), we can bound [ as
ac, (E<
12 E&1Q9 Jacc
T s || Q- Ac
1
F Fll 6 max 15j5n
mn,
6pmax duff
Pmux duef m.
E. Q = 0
m
We then move on to bound |||]| We start by rewriting
x Q2,
& E2
T
Q_
ATF
AC
<
Q_
ATF Ac
Tr
x_ (Be Q²
Q² )8,
x_ (B Q² ,
Ee 12 &
+
ES
Q
+
Eş Q-
m
m
T72
T78
35

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Note that 22 Q- (ATF' Ac
Be Re )
0. By adapting the techniques of bounding | 1/2 in
m
Appendix E.2, for dz m further get
x, B2BB 12,,2,,,,,
T71 s Bc
Be Q-
ACF Ac
I-cs
m
0 (mmaxheett m+
x. E Q2 &. X
12E
1 53
m
Furthermore, applying the Chebyshev's inequality again, and considering E |X_ Qq[ _, x /m] S 36pmmax d m,, alongside
<s/m]
x [,
(x, ExTE 1. x2)
Pinax de m²
Varg
S Bey
m2
it follows that, for x 2 72pmax ber// Pr(x, 102 11, x/x
/ m 2 x 10 S 6*
Pmaxx Lee (m²xx² Subsequently. by
0g <. m we deduce that
x, Be. (2. xx.
x 1 --
00
1xx x
X
Eg''
=
Pr
1 1-53 2 2 5" dr
1
0
m
TT Bc 18,8 &
I 2m 83 +
Pr
2 x 5') dx
2m²
m
<
m
2
+
6*
Pimax doff
m²
2m x.2
1
du
m
+
6 Punxx doft
2m
0
(Veine
m3
This leads to T71 = =0 =0(V Pinax des /m3
Next, we show a bound on T72. Using again (22) and (39), we get
T72 S
sup
xx_ (E (
) x (B.r [Q2]]
11
||||| =
|||[|||
m
x v
< 2 sup
Er'
Q )x)
||||=11|||||=
m
sup
E
x (E&/ Q
²²,)x, )},,
(uTQ Q Q v v
m [|u|
=1.llv|
=1
< 1 sup
|ac] (B. Q2 Q Q2,)ac,l
(u a Q u u + ac acc vv
||u|
j=1
ac (B Q s1
)ayyl
Σ u ac, Q_,u
< m7 ||u|
sup Ee
max 1515a
=1
j=1
ac T (Ber Q2 )ac,
n v accaacy v
sup Ey + m |v||=1
max 15J50
Mj
37
|ac T (Ber [22.1 - ,
max
E.
1515 S m
J
E (ac) (Bxr (2
< max 37 IS
(mTT))²
36

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Observe that adapting the proof of Lemma D.I, can readily ascertain the following:
(ac)) ac,) Pmax doff
R [ac]] _ (Rer [Q_.] - )ac)) Vare ac, acTTQ__acc [c, =
m
so that T72 0 Pinnx Recalling (39) again, derive
T73 1 1 x (1Eg .1 1 x)..,
Jacac (2. ac (EE 1Q2 E.,(Q
acj
I 12
m
jml
12Pmaxeefff
<
ac T (Eq[Q² E. 1.ac;
m
Along with
83
53
E E₂( 13
Es 1q
1 83
83
we obtain, for 63 I m 33
144ac acj
|ac, (E,1Q2,) Ec.(Q² Dac, 28 ac, (1 Q .]1l +
[Q=.
11<c,ll) s
m
This results in T73 = O (pmax pn/mm²)
Now, it remains to bound
| Noting (39) and Fusee
S 2, follows that
F22 i,
Q²
8,8
T
sup
|||||
|u]
7/3
m
Fi
F.
-1 11
x Q 8, Q X,
sup
Ec
|u]
7x
7.3
m
Pigts
u X Q 2²,x,00 x v
sup
E&
||]|=]1|||]|=1
7ss
m
Figis
u T 8 Q 8, Q 8.81
3 sup
E5
|||||=1_|v||=1
108 deff
Fijis
sup
E5
1
Q x
m
||]|=
3/8
Further, applying the techniques of bounding
H |22H2 1/21| in Appendix E.2 again, we get [T&
o o Putting the above together, we conclude that
(V &
|E. |B (
(La ATFAa)
m3
This thus concludes the proof of Proposition F.I.
Based on the second inverse moment result in Proposition F.1, we follow the line of arguments in Corollary 3.4 and Corol- lary 3.7,to derive the following fine- grained second inverse moment results for exact and/or approximate leverage score sampling and SRHT using scalar debiasing. The proofs of these results follow from direct combination of the proofs of Corollary 3.4, Corollary 3.7 and Proposition F.I,, and are omitted here.
37

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Corollary F.2 (Second inverse moment using scalar debiasing under approximate leverage). Under the settings and notations of Proposition F.I, for any random sampling scheme with sampling distribution mi [C/ /d/ttt (dett E / [(etttmmnn) with some Pmin E [1/2, 1] as in Definition 2.3 and diagonal matrix F diagg FiiggFFije Fi
with
aci Bell m-
A ASS + CA)-
F.
Then, there exists universal constant C O independent of n, delt such that for m Plaaddel log(den/0) with
5 I m mm-dar A SAc + CA) is an (c, -unbiisd estimator of Ia + Ac FAc with inversion bias c= max {o Pinax dit / m3) (CppPmaxxdet /m) for Cp = max 1, Pmax]
Corollary F.3 (Second moment for SRHT using scalar debiasing ). Under the settings and notations of Proposition for ASRHT Rmxn the SRHT of as Definition 3.6, and diagonal matrix F Giag(FF with
H 1/2 A AAST H 1/2 + CA) - ADDDH H, da dart
FR
TUTTTT
then there exists universal constant C 0 independent of n. deft and nexp(- < 8 < m-3 such that for m Cpmmadett log(detr/)) with max factor Pmax = maxicien né, (H, D.A////////eee m--eet H H-112 A SRUT ASRITH-1/22 + CA)-² is (e, (e.8ba-udd estimator of La + 1 A H, H,,,H, D,, D, Ac with inversion bias e O (VPimx mm3 / + deet log(n/5)/m)
These results are extension the results Corollaries 3.4 and 3.7, and can be used to derive SSN local convergence rates under exact/approximate leverage and SRHT similar to Theorem 4.3, See Corollaries F.5 and F.6 in Appendix F.3.
F.2. Proof of Theorem 4.3
Here, we prove Theorem 4.3 by adapting the proof approaches of our Propositions 3.2 and F.1 and Dereziński al. (2021a, Theorem 10). The proof of Theorem 4.3 comes in the following two steps:
we first establish Lemma F.4 that connects the de-biased SSN iterations 3,41 to the true 3r4 and
perform detailed convergence analysis of de-biased SSN.
To begin, we define a high probability event Ç. for each ( = O,..., independently for each iteration, and we denote Ç = = S. Under the setting of Theorem 4.3, follows from Lemma 2.7 that Pr((t) 2 8/T and Pr(() 2 5. These events will be constantly exploited in the remainder of the proof of Theorem 4.3.
, T 1, as defined in (19). These events are established
Proceeding to the core of our proof, present auxiliary lemma, which uses Propositions 3.2 and F.1 to connect the de-biased iteration 3 to the true Newton iteration 3 This is pivotal to prove Theorem 4.3. Note that this result is universally applicable any de-biased B, and is independent of the smoothness assumption on F in, e.g., Assumption 4.2. and C(3,) (2 a p.s.d. matrix, and let B be the de-biased SSN iteration as (10) with de-biased 8, as Proposition 3.2. Denote D = G, 3 and å, = =B.- 3.. If the exact Newton step B = B, H, gg (with g: the gradient of F at B,) is a descent direction, i.e., || < || then, there exists universal constant C > 0 independent of n, deut such that for m 2 Cemax (log( T/ 8 + 1/62/3) with e > 0, 8 < m3 we have
Lemma F.4. Let H, = V2 (B,) + C(A,) with (²) =
B 14 21 < | PHH + A. + 36pmax dell
| - -B112= ((8,))"
Here, Pmax is the max importance sampling approximation factor in Definition 2.3 for &. = max SIST (3,) and dell = maxigi5T dea (3,( with (3,) and dea respectively.
dea(3,) (3,) the leverage scores and effective dimension of A(B,) given C(B,).
38

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Proof of Lemma F.4. In both this proof and the subsequent proof of Theorem 4.3, simplify the notation, we abbreviate Ac (3,) Ac- We first recall some notations from (the proof of) Proposition F.1. Let Pt Bo+1 3 H, g
and Q. = (ATST AAAc + Ca)=1 with the de-biased sampling matrix S, = diag { Vv mm (3.))))) mi.))
m
a=1 S.
as in Proposition 3.2 and CA = H 1/2 H,
1/2. Define the diagonal matrix F = diag(F;)] with Fa = ac; T IEc, (A AS SAc + CA) -2 a/cc///n mn; as in Proposition F.l and ac; E IRd the ith row Ac.
Building on the results from Proposition 3.2 and Proposition F.1, for c O (V Pmmaxee /m² we obtain
c [|| ++
2 ||
li. E,, ||| | B
|.
24, H 2, 3+1+1111 31+1
= 2/4 A T
H
R Q Ia)] H h. g+ + /
H
c (Q: - Ia) H H ge
++1
s 2| A IH. |P:
|l. [c. [q] Iæll + H E3
- 2(Q: I4) -H Pt
x 2 2|4 ||p: [E [q:] Lall + |||p F |Eg. Q - Iall + p. P. Q
Ia ATF,ACH H.P.
+ P: H AFF,,
s 2c| A++1 H. PIllH + 2cIlp: + 1lPP: |B [Q - I. AFF, IIp:l FFll
c|| [p.||a. cllpalli.
max Isisn
ac:
EGg [ac
mn;
S E|| |p: |H. +elpeeli 36pmaxdeff
| (p.l(>=1(0.)>
which combined with ||Pl|H. 5 <11 11A + +11 S 2/| Dr IlH, and ||
Il,eell S 211|Brlii,, leads to
36pmaxdeft
Ex, ||l A e+ [H.]] - ||44 I.. S 2e 22e||A + 4el1 +
m
36pmax (
= 2 +
m
This concludes the proof of Lemma F.4.
Proof of Theorem 4.3. First note that by Assumption 4.2, we have the following result
(V Piinn mi
H2 Bella,,' S v/ A H. H. Re H, for U O
off
(41)
Letting á, 3, 9 and A, = Br 3r, we then show (41) holds under Assumption 4.2. Consider B, E U for all t, we get the following condition:
(Puaxcler
3/2
||
Ir,
m
where mmin denotes the smallest eigenvalue of H. We further derive
L 1---1
L |l S
V
Pinas m3
H)H $1
H.
H|
113.
3
3/2 min
U,
Umin
Tmin
Jmii
which leads to that H. u H. We now advance to verify the second part of (41). Recalling that the standard analysis of Newton's method (Boyd & Vandenberghe, 2004), we have
|ell=|HAA H.S. - (''2F(C
+ (f'VFI8 V'F(88 + rr))dT) L. S الÂ,ا L'
H gell gill=
[VF V²F(8" + +
< <11ă:11 1/ (1 -1- (--T)l(A,,
~ 1
39

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
which together with v < 1/2 and the fact that He R H implies that |H
1-2 achieves Cnin
u) | 2
+
2vI+
H. Be BellH,1
H,A: - gell S
S
||
| H
Omin
40
46 min
2
2. + v)
|Á S S014 v AI
2
Thus, (41) holds.
Next, we proceed to the main part of the proof. We first rewrite || ll.
|| | lir = 4 H (- B, + B, 8*) AT H, MA 18
= (1 - A T g + 141HHH4 g:)
= (1 - m + H. (1- ()) (H, gg) A gg))
[ H. m88 A.) p M (H,, ga)
= (1 - m4) | H. + +c( + (1 - -)) (H, -().
Invoking (41) and triangle inequality, we get
| F. S (1 - 11) 1 + topel A+t1 IH. I + Uuue11 - rell
With the fact that if x² + b then x² < a² + 2b, it follows that:
A l.. S r² A lt. + 2(1 (4) 11, + 204u((1 1ll = (2(1 - 1=)² 201 + n2/00² 2
< 2 -,)) + vn )| A 112.
then follows from Lemma F.4 and the inequality 42 (², = H, that
36pmax deff
Ex, 1.lll | 12 S v H.
|Polle..
(42)
m
For the second term in (42), we rewrite
36
36 36pmax deft
|
(g) T Å. -g. H. (H, A. g())
36pmax dell
H - AT H.A. gg) - -g. H (H, g:))
36pmax
36pmax
1 A 1|
( + H. T (H, - 8t)).
m
m
Using again (41), we further get
36pmax
36pmax dotf m
12 +
(H,,,, - 8:)
(- 11 A. ()) - -g:)
m
36 dear
Me(m + 22)/
m
Then, putting the above together, we have
36pmax 11 36pmax Me Molse + 2)) I1
Eg. (14 [1,11 5 (2(1- -2)2 20p4 +u+
m
m
= ((211 (2(1--p)) m2) (01)= + 36p 12 + 22 11+ 36pmaxdeett
Me (u4 + 21)) ||
=
m
m
40

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
which along with 2 = O (v Pinaxx oft/ / and
1-Pmaechat/mm < 1 results in
E (1] [i.11 (2(11-1)) 1x)² 36pmax ui 20) | <
+ (3t/mmaeee
3T (T
m
m
Applying v < and the fact that He Rig H indicates lIvIl F livl we get the following bound:
E S (1 + EG A Iim 1 (1 + uu) 37 PTpaxx
+ 2) || I1
$ (1 + v)² 37 Praxcchet + |
-( (Mmmaehhett
| ii:
(43)
m
where the constant "(1 v)3 37 is absorbed into m.
Now, it remains to check that 3, E U for all when conditioned on the event Ç = T Se. Assuming that this holds
for i 0, then it remains prove that, conditioned the event Ç, || S | || holds for each t almost surely.
Following the fact that Proposition 3.2 implies that, conditioned on Çt, we obtain
Qu La S E,
where Lemma 2.7 guarantees that is small. Adapting the techniques used in the proof of Lemma F.4 and by the analysis of the exact Newton step A++1, we arrive
| 1+ IH,
S | H. + 2|| H. PPl H. |1 Lall + P. H²
H. (Q:
1H) P.
S | li, +111, 12 Qu -all $ (e Praaxlefff
+2) ||
m
From Lemma 2.7, it is deduced that e Pmazde + is sufficiently small. Using He so H again, we deduce | | < || Consequently, we infer that every iterate lies within U, and the result (43) holds for O, 1, ....................... the above together, conclude that, conditioned on the event & that holds with probability at least 1 5, the union bound (11) is achieved, thereby completing the proof.
. T 1. Putting
E33. Additional Results on SSN and Further Discussions
Based on the above analysis, for any sampling method whose sampling probabilities close the exact approximate leverage scores, we then show that using scalar debiasing m/(m - doft) yields a slightly weaker local convergence result for SSN than that in Theorem 4.3, yet achieves enhanced computational efficiency. The result naturally follows by recalling Corollary 3.4 and Corollary F.2, and adapting the proof of Theorem 4.3 with a neighborhood U = {(3: ||3
3/2 (3'1 < "Mve Mii// [.} in place of U (3: ||3 BH (PmaxdettOmim/m /L), and the corresponding proof is
omitted for brevity.
Corollary F.5 (Local convergence of SSN using scalar debiasing under approximate leverage ). Under the settings and notations of Theorem 4.3, for any random sampling scheme with sampling distribution Ti E as in Corollary 3.4, there exists a neighborhood of 3* such
[60 (3) (3,)), & (B, (3))
that the de-biased SSN iteration Be41 B, Me m ( TT S.A(B)) c(B.)) g: starting from B E U -
satisfies, for U = 1B: | - 13 < Ulor 0 / step size 14 = 1 that
m//dat + Pmaa' and m 2 Ceaaxdea log (/)).
Pinxx
1/T
113 -3 11
Pinaxeefff (1 (1+E+Ep),
Ex
B. - -13"
m
holds for e e0 - m , eep = max (riit - 1, - pmine } and conditioned on an event Ç that happens with probability at least 1 8. Here, uhev (pmaxedest 332 + +pPmaxdeett //m, amin is the smallest singular value of H = AT (18))A((*)++ + C(3*), delt = maxisesT dent (B,) with (B,) and den the leverage scores and effective dimension of A(B,) given C(3,). respectively, Pmax max STST Pmax (B,) and Pmin = min SiST Pmii (3,.) with (B,) and Puii (3,) the (max and min) importance sampling approximation factors (of A(B,) given C(B,) in Definition 2.3.

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
Similarly, for SRHT, applying the scalar debiasing m/(m den) SSN also leads to slightly weaker local convergence than that in Theorem 4.3 but offers greater computational efficiency. This follows readily by recalling Corollary 3.7 and Corollary F.3, and adapting the proof of Theorem 4.3 with a neighborhood U = (3: || - 3 || < uSRITT 3/2
min/ instead of U- = {B: |/3 f || < (fmax / 3//2 / [.), and we omit the detailed proof for brevity.
Corollary F.6 (Local convergence of SSN for SRHT using scalar debiasing). Under the settings and notations Theorem 4.3, for ASRHT(B, (3,) = (3)) /n E R the SRIT of A as in Definition 3.6, there exists neighborhood U
-1
3" such that the de-biased SSN iteration B,+1 = B. - Me mh ma (, (8,)
m
À (3,)) (.A.RHT(B.) + C(33))
g:
starting from Bo E U satisfies, for U {B: ||3 - ,3 || < uSRHTO min/ },, step size pe 1
3/2
Pimxx
m 2
m/ dete + Piaace
Peaa log(denT/0) and T' exp( < 5 < I'm 3, that
11 33" -311111 1111
1/T
Pmaxxeeff (1 + E)2
Bc
(44)
113 3 |||
m
holds for O (Vrnaxdea// m + Pmnxx V dell log(nI/6) and conditioned on an event Ç that happens with probability
at least 1 - ô. Here, uSRNT = (Pmax (Pmaxdet/mm)31 3//2 + v deff (og(n//)) /m,, m is the smallest singular value of H = A TATTAAA (())) + C(B*), delf maxisisT dell (3,) with e (B,) and deul (3,) the leverage scores and effective dimension A /n given C(3,), respectively, Pmuxx maxisesT Praax ( (3,) with Praxx (B,) the max importance sampling approximation factor (of H,, (,, vT given C(B,) in Definition 2.3.
Remark F.7 (Comparison between Corollary F.6 and Lacotte et al. (2021, Theorem 2)). Using Corollary 3.7 and Corollary F.3 with V= = 0, it follows, by adapting the proof of Theorem 4.3, that for a number m > Cdee of samples, the linear convergence linear
rate in Corollary F.6 holds but with a right-hand side term of (1) in (44). That is similar to the convergence rate obtained in e.g., Lacotte et al. (2021, Theorem 2) for self-concordant f, 4, and F in (7). Notably, using larger number of samples m as in Corollary F.6 (than O )) the (linear) convergence rate in Corollary F.6 becomes dependent on Paaaelett/ /m, whereas the linear or quadratic rates stated in Lacotte et al. (2021, Theorem 2) cannot be characterized using doff or
G. Additional Numerical Experiments and Implementation Details
Here, we provide in Appendix G.1 and Appendix G.2 implementation details of the numerical experiments on SSN in Section 5, and then in Appendix G.3 additional numerical results on the inversion bias under approximate versus exact leverage score sampling.
G.1. Sketching Matrices and Step Size
For data matrix A Raxd, the approximate ridge leverage scores and leverage scores are computed using the method provided in Drineas al. (2012) and Cohen et al. (2017). The L.ESS-unifom sketching matrix is constructed in Derezifikii et al. (2021a, Section E.1). Further implementation details, including those for the SRHT, are available in the repository provided by Derezinski et al. (2021a) at https: : // github. com/ lessket ching/ hwwp//wwtttoeeeooeeeooeeee wwowwtonskttcch ch. Additionally, the Shrinkage Leverage Score sampling probability is formulated by combining uniform sampling probability and approximate leverage score sampling probability in equal proportions.
In our experiments, the first-order methods, Gradient Descent and Stochastic Gradient Descent, are used with fixed step size. As done in Dereziński al. (2021a), second-order methods, specifically Sub-sampled Newton and Newton Sketch with
Sketches, employ step sizes that dynamically adjusted using a line search algorithm based on the Armijo condition (Bonnans et al., 2006, Chapter 3).
For a fair comparison of the "convergence-complee Figures and 3 include the time for input data pre-pree pre-processing, e.g., the computation of exact or approximate leverage scores, and
trade-off across different optimization methods, the time reported in
transform, as well as the computational overhead associated with the sketching process. And we fix in Figures 2 and 3 the ridge regularization parameter to À = 10-2 2 for both MNIST and CIFAR-10 data in Section 5.
42

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
G.2. Datasets
For MNIST data matrix, we have 2 and d 27, and for CIFAR-IO data, we have 2 and d 2 We torchvision.tranese transforms from PyTorch to pre-process each image. We divide the ten classes of MNIST and CIFAR-10 datasets into two groups, assigning them labels of -1 and +1, respectively.
G.3. Inversion Bias for Exact versus Approximate Leverage Score Sampling under Scalar Debiasing
In this section, we present empirical experiments compare the inversion bias of exact approximate leverage score sampling under scalar debiasing m/(m deft), shown in Corollary 3.4 and discussed in Remark 3.5, for both MNIST and CIFAR-1 10 data.
Figure 4 depicts the inversion bias (with and without the scalar debiasing m/(m - dett) measured in spectral norm: (E
m
A ATST TTSSSA SA + A]a) (ATA + ALa) X and |E( T SA + ALa) (A A +
detr
/2 H A ALa) as function of the sketch size m. using the following random sampling schemes: ALA)
1. RLev: Exact X-ridge leverage score sampling with scalar debiasing m/(m - den.)
2. NO-RLev: Standard exact X-ridge leverage score sampling without scalar debiasing m/(m doft).
3. ARLev:
Approximate X-ridge leverage score sampling (identical to that used Section 5) with scalar debiasing m/(m dett), where the approximate ridge leverage scores & lle: AAAAAATTTAA S XIa) are computed using sparse Johnson Lindenstrauss transform (SJLT) (Clarkson & Woodruff, 2017) S1 E T x of size mi.
4. NO-ARLev: Standard approximate X-ridge leverage score sampling without scalar debiasing m/(m deft).
5. DARLev: A more efficient double-sketches variant of approximate X-ridge leverage score sampling, together with scalar debiasing m/ (m def)), ), where the approximate ridge leverage scores E. z ||e, A((A S A + XLa) -1/2g_] ) | are constructed using two SJLT matrices S E IR X and S2 E mgd with m2 < mi, see (Drineas et al., 2012; Cohen et al., 2017).
6. NO-DARLev: As DARLev but without the scalar debiasing m/ dott).
We observe from Figure that the inversion bias of all random sampling methods consistently decreases as m increases. Comparing solid to dashed lines, significant margin. We also find that RLev yields lower inversion bias compared its approximate counterparts ARLev and DARLev, under the scalar debiasing. These results corroborate the conclusions in Corollary 3.4 and Remark 3.5.
that the proposed scalar debiasing effectively reduces the inversion bias by
43

---

Fundamental Bias in Inverting Random Sampling Matrices with Application to SSN
RLey NO-RLev
109
ARL
NO-eR
Inversion bias
Inversion bias
NO-DARLev
10-1
10
500
1,000
1,500
2,000
2,000
3,000
500
1,000
1,500
2,000
2,500
3,000
3,500
4,000
Sketch size m
Sketch size m
(a) MNIST data
CIFAR- -IFAR-10cc data
Figure 4: Inversion bias as a function of the sketch size m, for various sampling methods on both MNIST and CIFAR-10 data, with ridge parameter À = 10-1 for MNIST data and À 10-6 for CIFAR-10 data. The results are averaged over 500 independent runs.
44


Reviews:
[{"review_id": "Gyktdfbgvd", "reviewer": "Reviewer_2Kx2", "summary": "This paper analyzes and corrects the inversion bias arising from common random sketching methods, including uniform/weighted sampling (based on leverage scores) and structured fast Johnson-Lindenstrauss transform like the sub-sampled randomized Hadamard transform (SRHT). Via a fine-grained characterization of inversion bias, it is demonstrated that for certain structured methods like (approximate) leverage-based sampling and SRHT, the bias can be effectively corrected through careful diagonal-wise rescaling. Leveraging these theoretical insights, the paper further establishes the first problem-independent local convergence rates for sub-sampled Newton methods, closely matching those obtained with dense Gaussian sketches. Empirical evaluations support these theoretical findings.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her constructive and thoughtful comments. Please find our point-by-point responses below.\n\n**Other Comments Or Suggestions**:\n\n\t1. Are the constants C in the theorems, propositions, and lemmas in Section 2 and 3 the same, or at least on the same order? It would be helpful to clarify this in the paper.\n\n**Reply**:\nWhile the constants $C$ in the theorems, propositions, and lemmas in Sections 2 and 3 are not identical, they are all *absolute constants* independent of the dimensions $n$, $d_{\\mathrm{eff}}$, and $m$.\n\n\n\t2. The preliminary in Section 2 seems to be a bit long, especially given the 8-page limit of the main text. Compared to the preliminaries, some proof sketches and intuitions in the appendix seem to be more insightful and may deserve a place in the main text.\n\n**Reply**: We thank the reviewer for this constructive comment. \nWe agree that the current section of preliminaries are a bit too long and that the proof sketches and intuitions (currently in the appendix) should be moved to the main text to better highlight the connection between RMT and RandNLA.\nDue to the 8-page limit for the submission, we were unable to include them in the current version. \n\nIf the paper gets accepted, we will revise the content accordingly (with the one additional page)  in the camera-ready version by incorporating the proof sketches and intuitions from Appendices D.1 and E.1 into the main text, placing them directly after Theorem 3.1 and Proposition 3.2, respectively. \n\n\n**Questions For Authors**:\n\n\tIs the log factors in m in most of the theorems, propositions, and lemmas in Sections 2 and 3 a consequence due to independent sampling of the m rows (uniformly or according to leverage scores)? If so, would it be possible to remove the log factors via dependent sampling, e.g., volume sampling, or some faster alternatives like adaptive leverage score sampling (see, e.g., [1])? It would be insightful to provide some discussions on this, even just as future directions.\n[1] Cortinovis, Alice, and Daniel Kressner. \"Adaptive randomized pivoting for column subset selection, DEIM, and low-rank approximation.\" arXiv preprint arXiv:2412.13992 (2024).\n\n**Reply**:\nWe thank the reviewer for the constructive comment and the reference. \nWe agree with the reviewer's perspective. \nAs shown in Section B of the supplementary materials, the subspace embedding property in Lemma 2.7 relies on the independence of the random row sampling procedure. \nThis implies that the logarithmic factors in $m$ present in most of the theorems, propositions, and lemmas in Sections 2 and 3 arise from the independent sampling of the $m$ rows (including, but not limited to, uniform sampling or leverage score-based sampling).\n\nAs such, it is very possible that dependent sampling methods, such as volume sampling or more efficient techniques like adaptive leverage score sampling (see, e.g., [1]), may provide tighter bounds and potentially eliminate the logarithmic factors.\nFor now, we are not able to include these efficient methods in the proposed analysis framework, since our current proof approach strongly relies on the independence of random sampling.\n\nWe also agree with the reviewer that these dependent sampling methods could yield better computational efficiency and/or convergence rate for SSN.\nWe leave them for future work. \nWe will include a discussion of these methods in the conclusion section (Section 6) of the revised version, addressing the potential benefits and challenges of using volume sampling, adaptive leverage score sampling, or other dependent sampling approaches.", "claims_and_evidence": "Yes, the main theoretical claims regarding the characterization and correction of inversion bias, as well as the convergence rates for sub-sampled Newton, are clearly presented. While I couldn't verify all the proofs in the appendix, with a quick scan, the proofs seem reasonable, and the intuitions provided at the beginning are helpful.  The empirical results are also convincing and support the theoretical claims.", "methods_and_evaluation_criteria": "Yes, the application of the proposed de-biasing method to sub-sampled Newton is well-aligned, providing convincing evidence for the broader applicability of the theoretical results.", "theoretical_claims": "While I couldn't verify all the proofs in the appendix, the main theoretical results in the paper seem reasonable. The proofs are well-structured, and the intuitions provided at the beginning are helpful.", "experimental_designs_or_analyses": "This is a theoretical paper with a few synthetic numerical experiments. The experiments are simple but sufficient to support the theoretical claims.", "relation_to_broader_scientific_literature": "Yes, this work provides a fine-grained characterization of inversion bias in common random sketching methods leveraging RMT techniques. It extends the existing results in Derezinski et al. (2021a;b) on approximately unbiased sparse embedding to more common random sketching methods. Applying the resulted de-biasing method to sub-sampled Newton leads to considerably improvement compared to Lacotte & Pilanci, (2019), Derezi ́nski et al., (2021a)., etc., both theoretically in terms of the convergence guarantee and empirically.", "essential_references_not_discussed": "To my knowledge, the paper discussed the essential references in the field.", "other_strengths_and_weaknesses": "Overall, I think this is a valuable theoretical study on inversion bias in random sketching methods. The paper is well-structured and easy to follow. The proposed de-biasing method is theoretically sound and practically useful. The application to sub-sampled Newton provides a good example of the broader applicability of the proposed de-biasing method. Some minor weaknesses are discussed below in the \"Other Comments\" and \"Questions For Authors\" sections.", "other_comments_or_suggestions": "Here are some minor points/questions:\n1. Are the constants $C$ in the theorems, propositions, and lemmas in Section 2 and 3 the same, or at least on the same order? It would be helpful to clarify this in the paper.\n2. The preliminary in Section 2 seems to be a bit long, especially given the 8-page limit of the main text. Compared to the preliminaries, some proof sketches and intuitions in the appendix seem to be more insightful and may deserve a place in the main text.", "questions_for_authors": "Is the log factors in $m$ in most of the theorems, propositions, and lemmas in Sections 2 and 3 a consequence due to independent sampling of the $m$ rows (uniformly or according to leverage scores)? If so, would it be possible to remove the log factors via dependent sampling, e.g., volume sampling, or some faster alternatives like adaptive leverage score sampling (see, e.g., [1])? It would be insightful to provide some discussions on this, even just as future directions.\n\n[1] Cortinovis, Alice, and Daniel Kressner. \"Adaptive randomized pivoting for column subset selection, DEIM, and low-rank approximation.\" arXiv preprint arXiv:2412.13992 (2024).", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "ZgMe0KVbUN", "reviewer": "Reviewer_GJpF", "summary": "The paper studies debiasing scheme for sub-sampled random matrix sketching. It improves prior work by providing a novel random sampling method which has better convergence result. Experiments are presented to corroborate theoretic results.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her constructive and thoughtful comments. Please find our point-by-point responses below.\n\n**Essential References Not Discussed**:\n\t\n\tNewton Meets Marchenko-Pastur: Massively Parallel Second-Order Optimization with Hessian Sketching and Debiasing (ICLR 2025, https://arxiv.org/abs/2410.01374)\n\n\tOptimal Shrinkage for Distributed Second-Order Optimization (ICML 2023, https://arxiv.org/abs/2402.01956)\n\n**Reply**:  We thank the reviewer for pointing to the two helpful references. \nThey will be included and discussed in the revised version of the paper.\n\n\n**Questions For Authors**:\n\n\t1. In Definition 2.3,$ \\rho_{\\min} \\equiv \\min_{1\\leq i \\leq n}\\ell_i^{C}/(\\pi_i d_{\\mathrm{eff}})$, isn't $\\pi_i d_{\\mathrm{eff}}=\\ell_i^{C}$ and thus cancels the numerator? Also, how does this sampling factor related to  $S$, I think it's just how much we weigh each row of $A$ when we sample, but the Definition sounds like $\\rho_{\\min} $ is also related to  $S$?\n\n**Reply**:  We thank the reviewer for this question. \nIn Definition 2.3, we define $\\rho_{\\min} \\equiv \\min_{1 \\leq i \\leq n} \\frac{\\ell_i^{\\mathbf{C}}}{\\pi_i d_{\\mathrm{eff}}}$. \nWhile for **exact (ridge) leverage score** (when the random sampling probability $\\pi_i$ is chosen according to $\\ell_i^{\\mathbf{C}}/d_{\\mathrm{eff}}$), one has $\\pi_i d_{\\mathrm{eff}} = \\ell_i^{\\mathbf{C}}$ and thus cancels the numerator, this is **not** true for generic random sampling methods.\n\nThe parameter $\\rho_{\\min}$ is introduced to capture the interaction between the choice of sampling probabilities $\\pi_i$ and the data structure of $\\mathbf{A}$ (e.g., its leverage scores) and how it affects the associated random sampling inversion bias.\n\n\n\t2. In Figure 2, the authors present numerical results showing the proposed method is much more efficient compare to Newton-LESS, I feel it would be more helpful if the authors can discuss more about the baseline and what's the main computation overhead there.\n\t3. Also in the experiments, the authors mention they exclude the time for computation of exact or approximate leverage score, does is mean the main workload of computing  $\\check S$ is removed? Is there any justification for such exclusion, i.e., is it a fair thing to do?\n\n**Reply**:\nWe thank  the reviewer for the comments.\nIt is worth noting that various second-order baselines in Figures 2 and 3, including Newton-LESS, ARLev, ALev, SLev, and SRHT, **all** require preprocessing of the data. \n\nSpecifically, as performed in [Derezinski2021a], Newton-LESS first performs the Hadamard transform preprocessing on the data, and then generates LESS-uniform sketching matrix $\\mathbf{S}$ to get the sketch $\\tilde{\\mathbf{A}} = \\mathbf{S} \\mathbf{A}$. \nThe computational costs of these two steps are $O(nd \\log n)$ and $O(msd)$, where $s$ represents the sparsity level of (each row in) $\\mathbf{S}$.  To ensure a fair comparison, we treat SRHT as a uniform sampling method that undergoes the same random Hadamard transform preprocessing as Newton-LESS. \n\nThe three leverage-based sampling methods, ARLev, ALev, and SLev, all involve the computation of approximate leverage scores and the sampling process, and has complexities of $O(\\text{nnz}(\\mathbf{A})\\log n + d^3 \\log d)$ and $O(m + md)$, respectively. \n\nThus, the preprocessing time is in fact **comparable** across all baseline methods, and by\nexcluding the time for the Hadamard transform and leverage score computation, we ensure a **fair comparison** when evaluating these methods. \nConsidering the computational overhead associated with the sketching process itself (rather than the preprocessing steps) allows for a clearer illustration of the \"convergence-complexity\" trade-off across different stochastic optimization methods.\n\nWe have reproduced Figures 1 and 2 by adding (for all methods) the preprocessing time, and observed a similar trends.\nIn particular, SSN-ARLev achieves a significantly lower overall runtime compared to Newton-LESS and establishes better complexity-convergence trade-off than other baselines.\n\nThese additional numerical results and discussions will be included in the revision.\n\n[Derezinski2021a] Derezinski M, Lacotte J, Pilanci M, et al. Newton-LESS: Sparsification without trade-offs for the sketched Newton update. Advances in Neural Information Processing Systems, 2021, 34: 2835-2847.", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "Yes", "theoretical_claims": "Not in full details", "experimental_designs_or_analyses": "Not in full details", "relation_to_broader_scientific_literature": "The key contribution is mainly a new random sampling scheme achieving better debiasing result, and connection between the new scheme and prior method (the scalar debiasing scheme, and phase transition, etc). It's a more refined version of prior debiasing literature.", "essential_references_not_discussed": "Newton Meets Marchenko-Pastur: Massively Parallel Second-Order Optimization with Hessian Sketching and Debiasing (ICLR 2025, https://arxiv.org/abs/2410.01374)\n\nOptimal Shrinkage for Distributed Second-Order Optimization (ICML 2023, https://arxiv.org/abs/2402.01956)", "other_strengths_and_weaknesses": "The paper is well-written, theoretically sound, and the experiments are clearly presented. I personally appreciate that the authors provide a lot of comments and discussions / remarks linking to prior findings, pointing out what's the similarity and what's the difference, which offers a comprehensive perspective for readers.", "other_comments_or_suggestions": "n/a", "questions_for_authors": "1. In Definition 2.3, $\\rho_\\min=\\min_{1\\leq i\\leq n} \\ell_i^C/(\\pi_id_{\\text{eff}})$, isn't $\\pi_id_{\\text{eff}}=\\ell_i^C$ and thus cancels the numerator? Also, how does this sampling factor related to $S$, I think it's just how much we weigh each row of $A$ when we sample, but the Definition sounds like $\\rho_\\min$ is also related to $S$? \n\n2. In Figure 2, the authors present numerical results showing the proposed method is much more efficient compare to Newton-LESS, I feel it would be more helpful if the authors can discuss more about the baseline and what's the main computation overhead there. \n\n3. Also in the experiments, the authors mention they exclude the time for computation of exact or approximate leverage score, does is mean the main workload of computing $\\check S$ is removed? Is there any justification for such exclusion, i.e., is it a fair thing to do?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bvxI3O2Xz8", "reviewer": "Reviewer_LRsm", "summary": "The paper studies the row-sampling of matrices in the context of approximating the inverse. Specifically, given an $n \\times d$ matrix $A$ and a row-sampling matrix $S$ with $m$ rows, sketching is commonly used to approximate $A^T A$ by $A^T S^T S A$. Prior work by Derezinski et al. showed that the inverse of the (rescaled) sketched matrix, $(\\frac{m}{m-d} A^T S^T S A)^{-1}$ is a ($1+\\epsilon$)-approximation to $(A^T A)^{-1}$ with only $m = O(d + \\sqrt{d}/\\epsilon)$. In contrast, a direct approximation of $A^T A$ typically requires $m = d \\log d / \\epsilon^2$ rows.\n\nThis paper generalizes the result of Derezinski et al. by studying how well $(A^T S^T SA + C)^{-1}$ approximates $(A^T A + C)^{-1}$, where $C$ is a PSD matrix. When $C$ is the identity matrix, this encompasses ridge leverage scores. This paper shows that $(A^T S^T SA + C)^{-1}$ is an ($1+\\epsilon$)-approximation to (A^T D A + C)^{-1} for some diagonal matrix D depending on ridge leverage scores and sampling probabilities, when the number of sampled rows $m$ satisfies $m \\gtrsim d \\log(m d) + d/\\epsilon^{2/3}$. When $C = 0$ and the sampling probabilities are exactly leverage scores, $D = (m-d)/m \\cdot I$, recovering the result of Derezinski et al. The paper then obtains results for some specific constructions of S, and applies the new results to the update step of Newton’s method, where it is used to approximate the inverse of the Hessian, and evaluates its performance through numerical experiments.\n\nOverall, this submission makes a solid contribution to the study of the inverse of row-sampled matrices, which also benefits the acceleration of optimization methods. I recommend acceptance, provided that my questions are adequately answered.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her constructive and thoughtful comments. Please find our point-by-point responses below.\n\n**Theoretical Claims:**\n\n\tI did not have time to check all the proofs. The following are a few questions:\n\n\t1. Page 4, left column, Line 199: what does it mean to say “subspace embeddings … ensure unbiasedness up to the same level of epsilon but are not guaranteed to remain unbiased for a smaller epsilon”?\n\n**Reply:** Sorry for this discussion that may be misleading.\nHere, we meant that sketching matrices satisfying the subspace embedding property, with an error $\\epsilon$ (say of order $\\epsilon = O\\left(\\sqrt{\\rho_{\\max} d_{\\mathrm{eff}} \\log(d_{\\mathrm{eff}}/\\delta)/m}\\right)$ ), can guarantee the **same** level $\\epsilon$ of the inversion bias, see Lemma B.1 in the supplementary material for a proof. \nIt is, however, unclear whether this (upper bound of) inversion bias is tight and/or can be made smaller with refined analysis.\nThis is in fact the very motivation of the proposed  de-biased sampling matrix $\\check{S}$, which, under the same sampling size, achieves a **smaller** inversion bias of order $\\epsilon = O\\left(\\sqrt{\\rho_{\\max}^3 d_{\\mathrm{eff}}^3 / m^3}\\right)$, see Proposition 3.2, offering an improvement over the standard subspace embedding-type analysis.\n\n\t2. Page 21, inequality (a) (in the middle of the page), can you explain how you get the constant factor of 38?\n\n**Reply:** We thank the reviewer for the constructive comments.\nRegarding the constant factor of 38 in inequality (a) on Page 21, we will provide a detailed explanation. \nBy combining $H^{1/2} Q_{-s} H^{1/2} \\preceq 6 I_d$ and the bound $|| H^{\\frac12} \\widetilde{H}^{-1} H^{\\frac{1}{2}}\\ || \\leq 2$, along with $||v || = 1$ and $||u || = 1$, we can derive the following two inequalities: \n$$\n\\sum^n_{j=1} u^\\top \\widetilde{H}^{\\frac12} Q_{-s} a_{j} a_{j}^\\top Q_{-s} \\widetilde{H}^{\\frac12} u=u^\\top \\widetilde{H}^{\\frac12} Q_{-s} A^\\top A Q_{-s} \\widetilde{H}^{\\frac12} u \\leq 36\n$$\nand \n$$\n\\sum^n_{j=1} v^\\top \\widetilde{H}^{-\\frac12} a_{j} a_{j}^\\top \\widetilde{H}^{-\\frac{1}{2}} v = v^\\top \\widetilde{H}^{-\\frac12} A^\\top A \\widetilde{H}^{-\\frac12} v  \\leq 2.\n$$\n\nAdding these two results gives the constant factor of 38. We will provide these proof details in Section D.2 of the supplementary material in the revised version for clarification. \n\n**Experimental Designs Or Analyses:**\n\n\tI do not see issues in general but I do not think one can claim “improved convergence” because the convergence is not improved. The main thing seems to be improving the runtime at the same convergence. The numerical experiments also show this: the convergence rate of SSN-ARLev looks about the same (at least the rate looks similar) as Newton-LESS in Figure 1 on both datasets.\n\n**Reply:** \nWe greatly appreciate the reviewer's helpful feedback and agree that the claim of \"improved convergence\" is indeed a bit misleading. \nAs discussed in Section 4 of the submission, the proposed de-biased SSN improves the local convergence rate over **standard** sub-sampled Newton method, to achieve problem-independent convergence rates (that, to the best of our knowledge, has **never** been established for SSN). \nIt is worth noting that the theoretical convergence rate of de-biasd SSN can approximately match, but **not** exceed, that of Newton-LESS (that is a closely-related, but different approach replying on random projection), but will **significantly reduced** computational complexity.\n\nIn particular, the proposed SSN-ARLev (that replies on random sampling) offers significant computational efficiency advantages over Newton-LESS, and thus outperforms Newton-LESS in terms of the complexity-convergence trade-off, see confirmed by Figures 1 and 2.\n\nWhile the numerical experiments in Figure 1 suggest that SSN-ARLev may offer **even improved** convergence in certain practical scenarios over the Newton-LESS approach, we are unable yet provide a theoretical proof of such superiority in convergence rate.\n\nAdditionally, we believe that combining the proposed de-biased approach with more efficient sub-sampled newton algorithms, such as the adaptive SSN-type algorithm proposed by [Lacotte2021], could further  accelerate convergence. \nHowever, exploring this is beyond the scope of the present work.\n\n\n[Lacotte2021] Lacotte J, Wang Y, Pilanci M. Adaptive newton sketch: Linear-time optimization with quadratic convergence and effective hessian dimensionality. International Conference on Machine Learning. PMLR, 2021: 5926-5936.\n\n**Other Comments Or Suggestions:**\n\n**Reply:**  We thank the reviewer again for the thorough and thoughtful review.  All minor issues raised will be addressed in the revised version of the paper.", "claims_and_evidence": "The claims are theoretical and proofs are provided.", "methods_and_evaluation_criteria": "The numerical experiments look reasonable.", "theoretical_claims": "I did not have time to check all the proofs. The following are a few questions:\n- Page 4, left column, Line 199: what does it mean to say “subspace embeddings … ensure unbiasedness up to the same level of epsilon but are not guaranteed to remain unbiased for a smaller epsilon”?\n- Page 21, inequality (a) (in the middle of the page), can you explain how you get the constant factor of 38?", "experimental_designs_or_analyses": "I do not see issues in general but I do not think one can claim “improved convergence” because the convergence is not improved. The main thing seems to be improving the runtime at the same convergence. The numerical experiments also show this: the convergence rate of SSN-ARLeve looks about the same (at least the rate looks similar) as Newton-LESS in Figure 1 on both datasets.", "relation_to_broader_scientific_literature": "Research on the inverse of row-sampling matrices appears to be less developed than that on row-sampling matrices themselves. This submission is a valuable contribution to the literature on the inverse of row-sampling matrices and also makes a nice addition to the literature on using sketching to accelerate numerical methods.", "essential_references_not_discussed": "NIL", "other_strengths_and_weaknesses": "NIL", "other_comments_or_suggestions": "- Page 4, left column, line 204: add a comma after $m$\n- Page 4, left column, line 207: as follow -> as follows\n- Page 5, left column, line 264: “n, eff” -> “n and d_eff”\n- Page 7, right column, line 338: It is better to say “Suppose that F, f: …”. \n- Page 16, Eq.(20): A right square bracket is missing for \\tilde{Z}_0.\n- Page 16, Lines 862-871: to be consistent with (22), it is better to write the event in the subscript, that is, $\\mathbf{1}_{\\neg\\zeta_3}$. \n- Page 21, the long list of derivation: Once u and v are separated, there is no need to carry both u and v for supremum; it suffices to have one supremum over u and the other one over v.", "questions_for_authors": "Please see the sections above.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7piFr8e42M", "reviewer": "Reviewer_vs1h", "summary": "This paper looks at how inversion bias (where the inverses of random sketches of a matrix are not unbiased) can be corrected for several random sampling methods. The paper gives an outline of how this can be done, and gives bounds for an $(\\epsilon,\\delta)$ estimator in Theorem 3.1 The paper makes a commentary on Derezinski's bounds (if I understand this correctly, the paper's derivation uses an alternative proof which beats the bounds given by Derezinski et al (2021a;b)). The paper then gives practical applications of these bounds to de-biased sub-sampled Newton with a better convergence rate (or convergence-complexity) trade-off.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for his/her constructive and thoughtful comments. Please find our point-by-point responses below.\n\n**Methods And Evaluation Criteria**:\n\n\tThe datasets are standard benchmark datasets (MNIST and CIFAR-10).\n\tI do have a question about lines 400-402: \"Figures 2 and 3 do not include the time for input data pre-processing......\" What is the relative magnitude of the time for input data pre-processing? That is to say, if the relative magnitudes of pre-processing is greater than the reduction in wall clock time for Newton-LESS / SSN-ARLev, then I don't think wall-clock time would be a meaningful metric.\n\n**Reply**:\nWe thank the reviewer for the constructive question (which is in fact shared by Reviewer GJpF).\nIt is worth noting that various second-order baselines in Figures 2 and 3, including Newton-LESS, ARLev, ALev, SLev, and SRHT, **all** require preprocessing of the data. \n\nSpecifically, Newton-LESS requires a Hadamard transform preprocessing with a complexity of $O(nd \\log n)$, see [Derezinski2021a]. \nTo ensure a fair comparison, we treat SRHT as a uniform sampling method that undergoes the same random Hadamard transform preprocessing as Newton-LESS. \nThe leverage-based sampling methods, ARLev, ALev, and SLev, require computing approximate leverage scores before constructing the sketches $\\check{\\mathbf{S}}\\mathbf{A}$ via random sampling, with complexity $O(\\text{nnz}(\\mathbf{A}) \\log n + d^3 \\log d)$.\n\nAs such, the preprocessing time are in fact **comparable** across these methods, and it makes sense to exclude the time for the Hadamard transform and leverage score computation ensures a fair comparison. \nConsidering the computational overhead associated with the sketching process itself (rather than the preprocessing steps) allows for a clearer illustration of the \"convergence-complexity\" trade-off across different stochastic optimization methods.\n\nWe have reproduced Figures 1 and 2 by adding (for all methods) the preprocessing time, and observed **a similar trends**.\nIn particular, SSN-ARLev achieves a significantly lower overall runtime compared to Newton-LESS and establishes better complexity-convergence trade-off than other baselines.\nThese additional numerical results and discussions will be included in the revision.\n\nAdditionally, for certain random Newton-type solvers, preprocessing is done **only once** at the beginning of the iteration and then reused in subsequent iterations. \nFor instance, in [Li2024], distributed least squares regression based on federated Newton iterations only performs data preprocessing (such as the Hadamard transform and computing leverage scores) once, and then reuse the obtained results alongside the optimization procedure. \nIn such cases, the preprocessing time becomes negligible compared to the iteration computation and communication costs. \nThus, we believe wall-clock time remains a meaningful metric.\n\n\n[Derezinski2021a] Derezinski M, Lacotte J, Pilanci M, et al. Newton-LESS: Sparsification without trade-offs for the sketched Newton update. Advances in Neural Information Processing Systems, 2021, 34: 2835-2847.\n\n[Li2024] Li J, Liu Y, Wang W. Fedns: A fast sketching newton-type algorithm for federated learning. Proceedings of the AAAI Conference on Artificial Intelligence. 2024, 38(12): 13509-13517.\n\n\n**Experimental Designs Or Analyses**:\n\n\tI would be interested to see the boxplots of the relative errors over the 30 independent runs in Fig 2, and 10 independent runs in Fig 3 for both Newton-LESS and SSN-ARLev, as the relative error for Newton-LESS and SSN-ARLev seem close to each other as the sketch size m increases. This is not to say that the experimental design / analyses is flawed, but I'm curious if the relative errors are skewed in one direction or equally symmetric about the solid lines in both Figures.\n\n**Reply**:\nWe thank the reviewer for the insightful comment. \nOur theoretical results say that the convergence rate of SSN-ARLev can approximate, but **not** exceed, the convergence rate of Newton-LESS, see Theorem 4.3. \nMoreover, as $m$ increases, the convergence rate of the de-biased sub-sampled Newton method approaches that of Newton-LESS. \nTherefore, it is reasonable to see in the experimental results that the relative errors for both Newton-LESS and SSN-ARLev becoming increasingly similar as $m$ increases. \n\nOur contribution in fact relies in showing that SSN-type methods can indeed achieve (local) **problem-independent convergence rates** (which, to the best of our knowledge, has **never** been established for SSN-type methods) that\n\n(1) can be significantly faster than *standard* SSN; and \n\n(2) matches that of Newton-LESS, but with significantly reduced computation complexity, and thus better complexity-convergence trade-off, as confirmed by Figures 1 and 2.\n\nAs suggested by the reviewer, we will include boxplots of the relative errors for both methods also converge as $m$ increases.", "claims_and_evidence": "Each claim is supported by proofs in the Appendix. While I did not have the opportunity to check each line thoroughly, each part of the (proofs) in the appendix gives the outline of the proof, and clear explanation for each subsequent step. This makes it straightforward for readers with more detailed background knowledge to follow and check.", "methods_and_evaluation_criteria": "The datasets are standard benchmark datasets (MNIST and CIFAR-10).\n\nI do have a question about lines 400-402: \"Figures 2 and 3 do not include the time for input data pre-processing......\" \nWhat is the relative magnitude of the time for input data pre-processing? That is to say, if the relative magnitudes of pre-processing is greater than the reduction in wall clock time for Newton-LESS / SSN-ARLev, then I don't think wall-clock time would be a meaningful metric.", "theoretical_claims": "The proofs look reasonable, however I did not carefully check the correctness of them.", "experimental_designs_or_analyses": "I would be interested to see the boxplots of the relative errors over the 30 independent runs in Fig 2, and 10 independent runs in Fig 3 for both Newton-LESS and SSN-ARLev, as the relative error for Newton-LESS and SSN-ARLev seem close to each other as the sketch size $m$ increases. This is not to say that the experimental design / analyses is flawed, but I'm curious if the relative errors are skewed in one direction or equally symmetric about the solid lines in both Figures.", "relation_to_broader_scientific_literature": "The contributions extend the work of randomized numerical algebra, where the paper looks at reducing the inversion bias of random sketches for random sampling methods. Based on the theorems and discussion in the main paper, the paper builds upon and extends (and critiques?) the results in the papers Derezinski 2021a;b using leverage based sampling.", "essential_references_not_discussed": "No", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "I will admit I am not very familiar with this paper, so I am basing my review solely on the clarity of the writing, some papers I referred to from the references in this paper, and questions from the experimental side.  But what I like about this paper is that the proofs are \"reproducible\", in the sense that the very detailed appendix explains what is happening in each line of the proof, as well as a general outline, hence anyone with sufficient background should be able to check what is happening at each step - unfortunately I am not the best person to do this.\n\nHence, please take this review with a grain of salt.", "questions_for_authors": "Please see above for my questions", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper considers the problem of inversion bias for random sampling matrices, which play a central role in randomized numerical linear algebra. Though previous works have shown how to handle specific samplings instances (e.g., Gaussian projection), this paper addresses the problem for general instances of random sampling. This is a very nice and well-written paper that provides a clear overview of the previous techniques and the reasons they do not immediately extend to the more general setting, along with their novel approach to debiasing based on approximate leverage score sampling. The reviewers unanimously agree that this paper should be accepted, and I also believe it provides important contributions that should be included in the conference.
Decision: Accept (oral)