Computing Approximate $\ell_p$ Sensitivities
Authors: Swati Padmanabhan
David Woodruff
Qiuyi Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: $\ell_p$ sensitivities
Lewis weights
leverage scores
approximation algorithms
total sensitivity
TL;DR: TL;DR: We give algorithms for approximating $\ell_p$ sensitivities and related summary statistics, with both theoretical guarantees and numerical experiments.

Abstract: Recent works in dimensionality reduction for regression tasks have introduced the notion of sensitivity, an estimate of the importance of a specific datapoint in a dataset, offering provable guarantees on the quality of the approximation after removing low-sensitivity datapoints via subsampling. However, fast algorithms for approximating sensitivities, which we show is equivalent to approximate regression, are known for only the $\ell_2$ setting, in which they are popularly termed leverage scores. In this work, we provide the first efficient algorithms for approximating $\ell_p$ sensitivities and other summary statistics of a given matrix. In particular, for a given $n \times d$ matrix, we compute $\alpha$-approximation to its $\ell_1$ sensitivities at the cost of $n/\alpha$ sensitivity computations. For estimating the total $\ell_p$ sensitivity (i.e. the sum of $\ell_p$ sensitivities), we provide an algorithm based on importance sampling of $\ell_p$ Lewis weights, which computes a constant factor approximation at the cost of roughly $\sqrt{d}$ sensitivity computations, with no polynomial dependence on $n$. Furthermore, we estimate the maximum $\ell_1$ sensitivity up to a $\sqrt{d}$ factor in $O(d)$ sensitivity computations. We also generalize these results to $\ell_p$ norms.  Lastly, we experimentally show that for a wide class of structured matrices in real-world datasets, the total sensitivity can be quickly approximated and is significantly smaller than the theoretical prediction, demonstrating that real-world datasets have on average low intrinsic effective dimensionality.

Computing Approximate lp Sensitivities
Swati Padmanabhan Massachusetts Institute of Technology pswtêmit edu
David Woodruff
Carnegie Mellon University dwoodruf@ee cmu edu
Qiuyi (Richard) Zhang Google Research qiuyi20google. com
Abstract
Recent works in dimensionality reduction for regression tasks have introduced the notion sensitivity, an estimate of the importance of a specific datapoint dataset, offering provable guarantees the quality of the approximation after removing low-sensitivity datapoints via subsampling. However, fast algorithms for approximating sensitivities, which show is equivalent to approximate regression, are known for only the l2 setting, which they are popularly termed leverage scores. In this work, approximating Ep sensitivities and related summary statistics of a given matrix. In particular, for a given n d matrix, compute o-approximation sensitivities at the cost of (m/a::) sensitivity computations. For estimating the total Ep sensitivity (i.e. the sum of l sensitivities), we provide an algorithm based on importance sampling of E, Lewis weights, which computes a constant factor approximation to the total sensitivity at the cost of roughly O(vd) sensitivity computations. Furthermore, factor, using (dd) sensitivity computations. We generalizeall these results to Ep norms for p > 1. Lastly, we experimentally show that for wide class of matrices in real-world datasets, the total sensitivity can be quickly approximated and is significantly smaller than the theoretical prediction, demonstrating that real-world datasets have low intrinsic effective dimensionality.
provide the first efficient algorithms for
its
estimate the maximum l sensitivity, up vd
Introduction
Many modern large-scale machine learning datasets comprise tall matrices A R with d features and n. >> training examples, often making it computationally expensive to run them through even ba- sic algorithmic primitives. Therefore, many applications spanning dimension reduction, privacy, and fairness aim to estimate the importance of a given datapoint as first step. Such importance estimates afford a principled approach focus our attention on subset of only the most important examples.
In view of this benefit, several sampling approaches have been extensively developed in the machine learning literature. One of the simplest such methods prevalent in practice is uniform sampling. Prominent examples of this technique include variants of stochastic gradient descent methods [1-4] such [5]14] developed for the "finite-sum minimization" problem
minimizeeex Σn= (0).
(1.1)
where, in the context empirical risk minimization, each fs : X' + R>0 in (1.1) measures the loss incurred by the i-th data point from the training set, and in generalized linear models, each fa repre- sents a link function applied to a linear predictor evaluated at the i-th data point. The aforementioned algorithms randomly sample function fi and make progress via gradient estimation techniques.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

However, uniform sampling can lead to significant information loss when the number of important training examples relatively small. Hence, importance sampling methods that assign higher sampling probabilities to more important examples have emerged both in theory [15, 16] and practice [17-21]. One such technique is based on sensitivity scores, defined [15] for each coordinate E [x] as:
f((x) Z==1 (x(
(1.2)
X
Sampling each fi independently with probability Ps o, and scaling the sampled row with 1/ps preserves the objective in (1.1) in expectation for every x E X. Further, one gets a
(1 + appproxiiini to this objective by sampling 3(Gd=-)) functions [22, Theorem 1.1], where G -=> o; is called the total sensitivity, and d is an associated VC dimension.
Numerous advantages of sensitivity sampling over é Lewis weight sampling [23] were highlighted in the recent work [24], prominent among which is that sampling matrix built using the E, sensitivities of A is significantly smaller than one built using its & Lewis weights in many important cases, e.g.. when the total sensitivity G is small for p > (which is seen in many structured matrices, such as combinations of Vandermonde, sparse, and low-rank matrices as studied in [25] as well as those in many naturally occuring datasets (cf. Section 4)). Additionally, sensitivity sampling has found use in constructing c-approximators for shape-fitting functions, also called strong coresets [15, 16, 26-28), clustering [16, 29, 30], logistic regression [31, 32], and least squares regression [33, 34].
Despite their benefits, the computational burden of estimating these sensitivities has been known to be significant, requiring solving an expensive maximization problem for each datapoint, which is a limiting factor in their use [35]. A key exception is for [((x) = (a] x)², when the sensitivity score of S is exactly the leverage score of a,, which can be computed quickly in (log(n)) sensitivity calculations, as opposed to the naive (n) sensitivity calculations (cf. Section 2).
1.1 Our contributions
In this work, we initiate systematic study of algorithms for approximately computing the &p sensitivities of a matrix, generally up to a constant factor. We first define ép sensitivities:
Definition 1.1. [15] Given a full-rank matrix A E Bgxd with n 2 d and scalar p E (0, 00), let a, denote the i'th row of matrix A. Then the vector of Ep sensitivities of A is defined as (ef. Section 2 for the notation)
a x|"
(a.) o
der
max XER"
|| p. for all i e [n],
(1.3)
and the total Ep sensitivity is defined as GA) = Eie(n) (
These f, sensitivities yield sampling probabilities for l, regression, a canonical optimization problem that captures least squares regression, maximum flow, and linear programming, in addition to appearing in applications like low rank matrix approximation [36], sparse recovery [37], graph-based semi-supervised learning [38-40], and data clustering [41]. However, algorithms for approximating lp sensitivities were known for only p- = 2 (when they are called "leverage scores"). As is typical in large-scale machine learning, we assume >> d. We now state our contributions, which are three-fold.
(1) Approximation algorithms. We design algorithms for three prototypical computational tasks associated with Lp sensitivities of matrix A E IRaxd Specifically, for an approximation parameter
1, we simultaneously estimate all l, sensitivities up to additive error of O aa (A)) via O(n/a) /a) individual sensitivity calculations, reducing our runtime by sacrificing some accuracy. We limit most of our results in the main text to l sensitivities, but extend these to all 2 1 in Appendix C. To state our results for the l, setting, use mnz (A) to denote the number of non-zero entries of the matrix A, w for the matrix multiplication constant, and we introduce the following piece of notation.
Definition 1.2 (Notation for ép Results). We introduce the notation LP(m, d, p) to denote the cost of approximating one l, sensitivity of an m x d matrix up to an accuracy of a given constant factor.
From hereon, for notational conciseness, omit stating Ax 0 in the problem constraint.

---

Theorem 1.3 (Estimating all sensitivities; informal version of Theorem 3.1 and Theorem C.2). Given a matrix A Rnxd and S < n, there exists an algorithm (Algorithm 1) which returns e
a vector o E 22 such that, with high probability, for each E [nl, we have
(a(a)) S o.
S (( (A)).
d More generally, for p 2 1, we can compute (via Algorithm 6) an estimate ō E 1220 satisfying o (a.) I o S 50(a= of , (a.) + a G,,())) with high probability for each i [n], at a cost of O (Anz(A) + (/
Our algorithm's runtime is Ô (nnz((A)
Two closely related properties arising widely in algorithms research are the total sensitivity (A) and the maximum sensitivity (A)]l As stated earlier, one important feature of the total sensitivity is that it determines the total sample size needed for function approximation, thus making its efficient estimation crucial for fast & regression. Additionally, it was shown in [24] that the sample complexity of sensitivity sampling is much lower than that of Lewis weight sampling in many practical applications such as when the total sensitivity is small; therefore, an approximation of the total sensitivity can be used as a fast test for whether or not to proceed with sensitivity sampling (involving the costly task of calculating all sensitivities).
While the total sensitivity helps characterize the whole dataset (as explained above), the maximum sensitivity captures the importance of the most important datapoint and is used in, e.g., experiment design and in reweighting matrices for low coherence [42]. Additionally, it captures the maximum extent to which a datapoint can influence the objective function and is therefore used in differential privacy [43]. While one could naively estimate the total and maximum sensitivities by computing all sensitivities using Theorem 1.3, we give faster algorithms that, strikingly, have no polynomial dependence on n. the assumed regime of >> d, this is a significant runtime improvement.
Theorem 1.4 (Estimating the total sensitivity; informal version of Theorem 3.5). Given a matrix A E Raxd, scalar p 2 1, and small constant < 1, there exists an algorithm (Algorithm 2), which returns a positive scalar 3 such that, with a probability of at least 0.99, we have
6)(A) 3 (1+ + ((0)))
Our algorithm's runtime is Ō (nnz(A) + 3 1d1p/- (1) d.p)).
In Algorithm 4 (presented in Appendix B.2), we give an alternate method to estimate 6 (A) based on recursive leverage score sampling, without l, Lewis weight calculations, which we believe to be of independent interest. All of our algorithms utilize approximation properties of Ep Lewis weights and subspace embeddings and recent developments in fast computation of these quantities.
Theorem 1.5 (Estimating the maximum sensitivity; informal version of Theorem 3.7 and
Theorem C.3). Given matrix A E Rnxd, there exists an algoritim (Algorithm 3) which returns
a scalar δ > O such that
|o x - (A)]]o
Our algorithm's 's runtime is Ô ((nzz()) + d=t1). Furthermore, via Algorithm 7, this result holds for all p E [1,2] and for p > 2, this guarantee costs Ö (nnz(A) + dp/2 LPP(O))/2)) d)p).
(2) Hardness results. In the other direction, while is known that e, sensitivity calculation reduces to lp regression, we show the converse. Specifically, we establish hardness of computing E, sensitivities in terms of the cost of solving multiple corresponding lp regression problems, by showing that
(1 + approximation to E, sensiiiiiisss solves Ep multi-regression up to an accuracy factor of 1 + E.
Theorem 1.6 (Leave-One-Out Ep Multiregression Reduces to lp Sensitivities; informal ) Suppose that we are given an sensitivity approximation algorithm A, which for any matrix A' E R xd' and accuracy parameter e' E (0, 1), computes (1 = '')a,((A') in time T(r d', nnz( A'),e"), Given matrix A Raxd with 2 d, let OPT, == minyeRs .IAA,,y A and y arg minyegd-1 (AA..y+ A for all the i [d]. Then, there exists an algorithm that takes A e n and computes (1 + =POTT for all in time T T(n + d, d, nnz(A),
2 Pleease see Lemma D.2 for the formal version with small additive error.
3

---

For p * 2, approximating all E sensitivities to a high accuracy would therefore require improving 4p multi-regression algorithms with (2(d)) instances, for which the current best high-accuracy algorithms take poly(d) nnz( A) time [44, 45]. More concretely, our hardness results imply that in general one cannot compute all sensitivities as quickly as leverage scores unless there is major breakthrough in multi-response regression. In order to show this, we introduce reduction algorithm that efficiently re- gresses columns of A against linear combinations of other columns of A, by simply adding a row to A capturing the desired linear combination and then computing sensitivities, which will reveal the cost of the corresponding regression problem. We can solve ((d) of these problems by augmenting the matrix to increase the rows and columns by at most constant factor. We defer these details to Appendix D.
(3) Numerical experiments. We consolidate our theoretical contributions with a demonstration of the empirical advantage of our approximation algorithms over the naive approach of estimating these sensitivities using the UCI Machine Learning Dataset Repository [46] in Section 4. We found that many datasets have extremely small total sensitivities; therefore, fast sensitivity approximation algorithms utilize this small intrinsic dimensionality of real-world data far better than Lewis weights sampling, with sampling bounds often a factor of 2-5x better. We also found these sensitivities to be easy to approximate, with the accuracy-runtin runttmmme tradeoff much better than our theory suggests, with up to 40x faster in runtime. Lastly, we show that our algorithm for estimating the total sensitivity produces accurate estimates, up to small constant factors, with a runtime speedup of at least 4x.
1.2 Related work
Introduced in the pioneering work of [15], sensitivity sampling falls under the broad framework of 'iimportance sampling". It has found use in constructing e-approximators of numerical integrands of several function classes in numerical analysis, statistics (e.g., in the context of VC dimension), and computer science (e.g.. in coresets for clustering). The results from [15] were refined and extended by [16], with more advances in the context of shape-fitting problems in [26-28], clustering [16, 29, 30], logistic regression [31, 32], least squares regression [33, 34, 47], principal component analysis [48], reinforcement learning [49, 50], and pruning of deep neural networks [51-53]. Sampling algorithms for lp subspace embeddings have also been extensively studied in functional analysis [54-59] as well as in theoretical computer science [23, 60, 44, 45].
Another notable line of work in importance sampling uses Lewis weights to determine sampling probabilities. As we explain in Section 2, Lewis weights are closely related to sensitivities. Many modern applications of Lewis weights in theoretical computer science we introduced in [23], which gave input sparsity time algorithms for approximating Lewis weights and used them to obtain fast algorithms for solving lp regression. They have subsequently been used in row sampling algorithms for data pre pre-processing [33, 61, 62, 60, 23], computing dimension free strong coresets for m-median and subspace approximation [63], and fast tensor factorization in the streaming model [64]. Lewis weights are also used for & regression in: [65] for stochastic gradient descent pre-conditioning, [66] for quantile regression, and [67] to provide algorithms for linear algebraic problems in the sliding window model. Lewis weights have also become widely used in convex geometry [68], randomized numerical linear algebra [23, 42, 69, 64, 70], and machine learning [71-73].
2 Notation and preliminaries
We use boldface uppercase and lowercase letters for matrices and vectors respectively. We denote the ith row vector of matrix A by a, Given a matrix M, we use Tr(M) for its trace, rank(M) for its rank, M for its Moore-Penrose pseudoinverse, and MM] for the number of its rows. When æ is an integer, use [xr] for the set of integers 1, 2,. ..... x. For two positive scalars æ and y and scalar 10, we use 2 Ra y to denote (1 a)y S x < (1 + a)y; we sometimes use x A y to mean (1-c)y < x < (1 + c)y for some appropriate universal constant We use Ô to hide dependence on (1))) We acknowledge that there is some notation overloading between p (when used to denote the scalar in, for example, fp norms) and p. (when used to denote some probabilities) but the distinction is clear from context. We defer the proofs of facts stated in this section to Appendix A.
Notation for sensitivities. We start with a slightly general definition of sensitivities: a (a.) to denote the E, sensitivity of a, with respect to B; when computing the sensitivity with respect to the

---

same matrix that the row is drawn from, we omit the matrix superseript:
a x P
max a x Ax 1
o B (a:) max
max
Bx p
and pa((a)=
(2.1)
When referring to the vector of pp sensitivities of all the rows of matrix, omit the subscript in the argument: so, o (A) is the vector of Ep sensitivities of rows of the matrix A, each computed with respect to the matrix B (as in Equation (2.1)), and analogously, o (A) op(A ) is the vector of Ep sensitivities of the rows of matrix A, each computed with respect to itself. Similarly, the total sensitivity is the sum of bp sensitivities of the rows of A with respect to matrix B, for which we use the following notation:
G
B (A)
Σ of and 6 (A) iE[IA
Σ oa(a.). [E[AI]
(2.2)
where, analogous to the rest of the notation regarding sensitivities, we omit the superscript when the sensitivities are being computed with respect to the input matrix itself.
While of(a (a.) could be infinite, by appending a to B. the generalized sensitivity becomes once again contained in [0, II. Specifically, simple rearrangement of the definition gives us the identity o BBa. (a,) 1/(1 + 1/o0 ( (a,)). We now define leverage scores, which are a special type of sensitivities and also amenable to efficient computation [60], which makes them an ideal "building block" in algorithms for approximating sensitivities.
Leverage scores are Eg sensitivities. The leverage score of the 2th row a of A is (A) def a (A A) Ta, The zth leverage score is also the z diagonal entry of the orthogonal projection matrix ((A T A) A T, Since the eigenvalues of an orthogonal projection matrix are zero or one,
implies Ti A) = 1! A (A A) 1, S 1 for all i [ni] and Σx= T; S d (see [74]).
turns out that T; (A) og (2(a:) (see (75]). This may be verified by applying a change of basis to the
variable y Ax in the definition of sensitivity and working out that the maximizer the vector par- allel to y = A A'a, This also gives an explicit formula for the total l sensitivity: (A) = rank ()).
The fastest algorithm for constant factor approximation to leverage scores is by [60], stated next.
Fact 2.1 ((60, Lemma 7).. Given matrix A E Bx we can compute constant-factor approxima- tions to all its leverage scores in time O(nzz( A) log(n) + d log(d)).
In this paper, the above is the result we Indeed, since the error accumulates multiplicatively across different steps of our algorithms, a high-accuracy algorithm for leverage score computation does not serve any benefit over this constant-accuracy one. Leverage scores approximate & sensitivities up to a distortion factor of Vn, as seen from the below known fact we prove, for completeness, in Appendix A.
compute leverage scores (to constant accuracy).
Fact 2.2 (Crude sensitivity approximation via leverage scores). Given matrix A Raxd let o, (a:) denote the ith En sensifivity of A with respect to A, and let (a(A) denote its ith leverage
score with respect to A. Then have [T((A) , (a.) S (A).
In general, Ep sensitivities of rows of matrix suffer distortion factor of at most 1/12211// 1/ml from corresponding leverage scores, which follows from Hölder's inequality. The lp general- izations of leverages scores are called l, Lewis weights [54, 23], which satisfy the recurrence w; = Ti( 1W112-1/P A). However, unlike when p 2, &p Lewis weights are not equal to sensitivities, and from an approximation perspective, they provide only a one-sided bound on the sensitivities (a.) < d m=ax(0,p/2-1) A (a:) [75, 76). While we have algorithms for efficiently computing leverage scores [60] and Lewis weights [23, 77], extensions for utilizing these ideas for fast and accurate sensitivity approximation algorithms for all have been limited, which is the gap this paper aims to fill. A key regime of interest for many downstream algorithms is a small constant factor approximation to true sensitivities, such as for subsampling [78].
Subspace embeddings. We use sampling- baasdd casstantapprooximatio lp subspace embeddings that we denote by Sp. When the type of embedding is clear from context, omit the subscript.
5

---

Definition 2.3 (lp Subspace Embedding). Given A E Rnxd (rypically n >> d), we call S E R (typically d S < ni) an e-approximate Ep subspace embedding of A if |A -= p Vx E R
When the diagonal matrix S defined as S = L/vPi with probability pt = (A)) (and 0 otherwise), then with high probability SA is an E2 subspace embedding for A [33, 79]; further, S has at most O(ds- non-zero entries. For general p. we use Lewis weights to construct lp subspace embeddings fast [23, Theorem 7.1J. This fast construction has been made possible by the reduction of Lewis weight computation to few (polynomial in p) leverage score computations for p 4 [23] as well as p 2 4 [77]. This efficient construction is the reason for our choice of this subspace embedding. We note that there has been extensive amount of literature designing { subspace embeddings for example using Cauchy [80] and exponential [81] random variables.
For S 2, it is known that the expected sample complexity of rows is Σ = O(d) [34), and for p > 2, this is known to be most (dr/ [35]. While these bounds are optimal in the worst case, one important application of faster sensitivity calculations is to compute the minimal subspace embedding dimension for average total sensitivity 6, (A), which has been used to perform more sample efficient subspace embeddings [67, 27, 28]. as sampling Ö(G, (A)d/e) rows suffice to provide an ép subspace embedding.
applications, particularly those seen practice. This gap captured by the
Accurate sensitivity calculations. While leverage scores give crude approximation, constant- factor approximations of l, sensitivities seem daunting since we must compute worst-case ratios over the input space. However, we can compute the cost of one sensitivity specific row by reducing it to lp regression since by scale invariance, we have
mil a x=1
Bx[]
(a.)
This reduction allows us use recent developments in fast approximate algorithms for Ep regression [44, 25, 45], which utilize subspace embeddings to first reduce the number of rows of the matrix to (poly(dd)) For the specific case of p 1, which is the focus of the main body of our paper, we may reduce to approximate & regression on a O(d) x d matrix, which is a form of linear programming.
Fact 2.4. Given an x matrix, the cost of computing is of its & sensitivities is Ó(nnz + ks d
3 Approximating functions of sensitivities
We first present a constant-prote -robbbbiip algorithm approximating the l sensitivities in Algorithm 1. Our algorithm is a natural one: Since computing the é sensitivities of all the rows simultaneously is computationally expensive, we instead hash a-sized blocks of random rows into smaller (constant-sizeet buckets and compute the sensitivities of the smaller, generated, computing the sensitivities with respect to SA, the & subspace embedding of A. Running this algorithm multiple times gives our high-probability guarantee via the standard median trick.
-nized matrix P so
Theorem 3.1. Given full-rank matrix A E Raxd and an approximation factor 1 a & let o (as) be the l sensitivity of the ith row of A. Then there exists an algorithm that, in time δ (nnz(A A) + n. d returns ō E Ro such that with high probability, for each i E [n],
o, (a,) S S S O(a,(())) + (A)).
(3.1)
Proof Sketch of Theorem 3.1: full proof in Appendix B.I. We achieve our guarantee via Algorithm 1. To see this, first note that ||SAx||| ~ @[]] (since as per Line 1, SA is an é subspace embed- ding of A). Combining this with the fact that every row in mapped, via the designed randomness, to some row in the matrix P helps establish the desired approximation guarantee. The runtime follows from using Fact 2.4 to compute E sensitivities with respect to SA of PP = O(n/a) rows.
As seen in Appendix C.I, our techniques described above also generalize to the p 2 1 case. Specifically, in Theorem C.2, we show that reducing E, sensitivity calculations by an a factor gives an approximation guarantee of the form (apP-oo, + o x (A())
Remark 3.2. The sensitivities returned by our algorithm are upproximate, with relative and additive error, but are useful in certain settings as they preserve E, regression approximation guarantees while increasing total sample complexity by only a small poly(d) factor compared to true sensitivities while still keeping it much smaller than that obtained via Lewis weights. To see this with a simple toy

---

Algorithm Approximating & :--ensitiviitie: Row-wise Approximation
Input: Matrix A E Raxd approximation factor a E (1,n)
Output: Vector o R5o that satisfies, for each E [n]]. with probability 0.9, that
(a S ã, s (a.) (A)))
Compute for an & subspace embedding SA O(d) d (cf. Definition 2.3)
2: Partition A into blocks B....., Bn/a each comprising a randomly selected rows
3: for the f block B₃, with E [ do
4: Sample a-dimensional independent Rademacher vectors
F100
5: For each jE()) compute the row vectors r B Bg IRd
6: end for
7: Let P I. 1o xd be the matrix of all vectors from Line 5. Compute of (P) using Fact 2.4
8: for each i [r]] do
9: Denote by J the of row indices in that a, is mapped Line 5
10: Set ã; maxje (p,))
end for
12: Return σ
example, consider the case >> 2. Note that by Theorem 1.5 of [24], the sample complexity using approximate sensitivities is 20222-26) 2-22/p 22-2/P(AA) (A)), using the aP factor blowup via Theorem C.2, that using true sensitivities is Ō
2- 2-/ (A)), and that using Ep Lewis weights is (dd/2), Suppose also that the total sensitivity (A) d. Further assume that n d and n tūp dr; then the dimension- dimension-deeendenee of sample complexities given by our approximate sensitivities, true sensitivities, and Lewis weights are, respectively, O(d"), and (dp/2)) (ignoring small factors in the expo- nent that do not affect the asymptotic unalysis). Thus, our approximate sensitivities provide a tradeoff in compulational cost and sample complexity between using true sensitivities and Lewis weights.
3.1 Estimating the sum of lp sensitivities
In this section, we present a one-shot algorithm that provides a constant-pprroximation the total lp sensitivity for all 2 1. This algorithm is based importance sampling using Ep Lewis weights. For our desired guarantees, we crucially need the approximation bounds of lp Lewis weights, as provided by Fact 3.3, efficient computation of E,, Lewis weights as provided by [23] for < 4 and by [77] for p 2 4, and the fact that random 6, sampling matrix (cf. Fact 3.4) yields SA, constant-approximatinn subspace embedding to A [24, Theorem 1.8).
Algorithm 2 Approximating the Sum of f Sensitiviiiss for p p2
Input: Matrix A E padd approximation factor 7 (0.01, 1). and scalar p 2 Output: A positive scalar that satisfies, with probability 0.99, that
G I $ (1 (
1: Compute all 4p Lewis weights of A, denoting the 1th weight by wp(a:)
2: Define the sampling vector v E IR by 0; (a().) for all E [n]]
O (dll-p//2) rows with replacement, picking the 7th row with probability of 4: Construct a random f, sampling matrix S.A with probabilities (wi)a_ (cf. Fact 3.4)
3: Sample
a SpA
(m:,)
5: For each sampled row iz. j E [m], compute rj
6: Return 1 Σ
Fact 3.3 (Lewis Weights Bounds Sensitivities, [42, 82]). Given matrix A E Raxd and p 2 1, the lp sensitivity o (ai) and l, Lewis weight w,(a:) of the row a satisfy, for all [n].
0 (a.) (a.) < d mmax(00p/221)) wy(a))
7

---

Fact 3.4 (Sampling via Lewis Weights [23, 35]). Given A E Ruxd and O. Consider random diagonal ep sampling matrix E R=x= with sampling probabilities {p,} proportional to the lp Lewis weight of A, i.e., for each E [n], the 7th diagonal entry is independently set to be
1 1/P
1/p
with probability
S
0
otherwise
Then, with high probability, S with O (ile-?mmax(ii d_max(L_//2) (log d) log(d/e)) rows is by subspace
embedding for A (cf. Definition 2.3).
Theorem 3.5. Given a matrix A E RRadd and an approximation factor E (0, 1), there exists an algorithm, which returns a positive scalar 8 such that, with a probability 0.99, we have
(x(A) S 8 < (1 ((1)),
Our algorithm's 's runtime is Ô (nz(A) + (0 d..)).
Proof. Without loss of generality, we may assume A be full rank. Then, its Lewis weights satisfy E (P((a.) = d. Per Line of Algorithm 2, our sampling distribution is chosen to U; w.(n.)) for all i E [r]. We sample from A rows with replacement, with row picked with probability of U/. From the definition of v; in Line and r in Line 5 and the fact that S,A.A. A a constant factor 4 subspace embedding of A, our algorithm's 's output satisfies the following unbiasedness condition:
m
n
S,A.A.
E
Σ
12
Vig
G. (A) (A)
M11
m
Vig
jE(m)
By independence, we also have the following bound on the variance of m EjE(m]
1
S
n S,A.A (a,)
Σ
a
Var
T2
m
m
U;
Wp(a:)
JE(mJ
1=1
i=1
with the final step holding by the choice of n in Line 2. In the case that p 2 2, we have
d
aa.) 2
d
2dp/2
Var rj
o 112-1
(
(A),
(a,))
5E[m]
REIm]
where the first inequality uses Fact 3.3. Therefore, applying Chebyshev's 's inequality on 1 m Eje(m" Fj
(as defined in Line 5) with m = 0
", gives us a ?-multiipicaative accuracy in approximating
/p/2
G (A) (with the desired constant probability). For p 2, we additionally have [24, Theorem 1.7] the lower bound 6, (A) d, which when plugged into the value of m gives the claimed sample complexity, A similar argument may be applied for the case < 2: specifically, we have that
S,A
d
2d
Var
Tj
Σ TE|M]
S,A (a.)
G,(A), m
m
m
jE(m]
where the second step is by wp (p(a:) 2 o (a.)
(a.) from Fact 3.3. For p < 2, we also have 6, dp/2 from [24, Theorem 1.7]. Applying Chebyshev's inequality with this fact gives a sample complexity of m = (d-- p22 This completes the correctness guarantee.
Runtime. We first compute all Ep Lewis weights of A E RRaxd up to constant multiplicative accuracy, the cost of which is O (0 (1-1i-/a log(log(n))) leverage score computations for p < 4 [23] and (og(p)) leverage score computations for p 2 [77]. Next, we compute m 1p/2)) sensitivities with respect to S.A.A From [23, 35], Sp has, with high probability, (dd/2 log d) rows when > and (dd log d) rows when p S 2. Summing the cost of computing these sensitivities and that of computing the Lewis weights gives the claimed runtime.
Remark 3.6. We present in Appendix B.2 an alternate algorithm, Algorithm 4, for estimating the total Ep sensitivity for p = 1. Algorithm uses recursive computations of leverage scores, in contrast to Algorithm 2 which uses Lewis weights in one-shot manner, and may be of independent interest.
8

---

3.2 Approximating the maximum 4 sensitivity
In this section, we present an algorithm that approximates l|o (A)ll maximum of the & sensitivities of the rows of a matrix A. As alluded to in Section 1, a first approach to this problem would be to simply estimate all sensitivities and compute their maximum. To do better than this, one idea, inspired by the random hashing approach of Algorithm 1, is as follows.
|| o (a.). the
If the matrix has large number of high-sensitivity rows, then, intuitively, the appropriately scaled maximum sensitivity of a uniformly sampled subset of these rows should approximate ||or [r] (A) lloe- Specifically, assume that the matrix has least a rows of sensitivity at least [|gr(A) loc/ Va: uniformly sample n// a rows and return a the (appropriately scaled) maximum of 0 (a.). for the sampled rows i. Then, ||o (A) lee I oa S =Ovallo a (o) with a constant probability. Here the upper bound guarantee is without any condition on the number of rows with large & sensitivities.
In the other case, if the matrix does not have too many high-sensitivity rows, we could estimate the maximum sensitivity by hashing rows into small buckets via Rademacher combinations of uniformly sampled blocks of a rows each (just like in Algorithm 1). Then, d. the scaled maximum of the & sen- sitivities of these rows satisfies |o (A)llo S os O(vallo ( >oo) Here, it is the lower bound that comes for free (i.e., without any condition on the number of rows with large & sensitivities).
Despite the above strategies working for each case, there is way combine these approaches without knowledge of whether the imput matrix has a enough "iigh-sensttivtt"" rows or not. We therefore avoid this approach and instead present Algorithm 3, where make use of A, an lo subspace embedding of A. Our algorithm hinges on the recent development [75, Theorem 1.3] on efficient construction of such an embedding with simply subset of (dd) rows of A.
Algorithm 3 Approximating the Maximum of & -SSensitiviites
Input: Matrix A E 1
Output: Scalar 3 e Rzo that satisfies |o (Ax S
1: Compute, for A, an &or subspace embedding S E (d) such that S A a subset of the rows of A [75, Theorem 1.3]
2: Compute, for A, an é subspace embedding S A Rodxdd
3: Return v o S. S.
AS A)
Theorem 3.7 (Approximating the Maximum of & Sensitivities). Given matrix A E R there exists an algorithmm, which in time ((nnz(A) + d=+1), outputs a positive sealar that satisfies
(A)l(( s O(Vd)o
Proof Sketch of Theorem 3.7; full proof in Appendix B.3. We achieve our guarantee via Algorithm 3. Define x" and a,- as: x" - arg maxxeR", [E[|A]]
x Since S,A is l subspace embedding of A, we have, for any x E Rd that Ax|||| = ((I If S A contains the row a,, then |x|,Ax]]==(1 Ax|l1) implies || ( A) = (()] (A) A)llo In the other case, suppose
|a][x]
[Ax]]
ar- not included in A. Then we observe that
c x
Ax O(1) max XER
SSSee Ax 20
| S. S.A (Sx
max XER ee E
||
max XER
||
|Ax|]
(3.2)
where the the second step is by choosing specific vector in the numerator and the third step uses |xx|| = -x|]Ax]11. We further have,
SS Ax|ex
Ax* 1ao
|Ax*
|axx"| v Ax+
1 = Vāl 1 (A))lo,,, (3.3)
max XER
|
|AxlI
vdl|Ax
where the first step is by choosing = x", the second step is by the guarantee of foo subspace embedding, and the final step is by definition or (a..). Combining Inequality (3.2) and Inequal-
ity (3.3) gives the claimed lower bound on -|o5 S (S A) The runtime follows from the cost of
computing S A from [75] and that of computing (dd of l sensitivities with respect A.A. D

---

addler FRlt
Figure Average absolute log ratios e, sensitivity approximations for wine, theoretical additive bound (dashed).
4 Numerical experiments
We demonstrate our fast sensitivity approximations on multiple real-world datasets in the UCI Ma- chine Learning Dataset Repository [46], such as the wine and fires datasets, for different and vary- ing approximation parameters a. focus on the wine dataset here, with full details in Appendix
Experimental Setup. For each dataset and each p, first apply the Lewis weight subspace embedding to compute smaller matrix SA. Then, we run our sensitivity computation Algorithm and measure the average and maximum absolute log ratio [o mme) to capture the relative error of the sensitivity approximation. Note that this is stricter than our guarantees, which give only an O(a) additive error; therefore we have no real upper bound on the maximum absolute log ratio due to small sensitivities. We plot an upper bound of log(a")) which is the worst case additive error approximation although it does provide relative error guarantees with respect the average sensitivity. We compare our fast algorithm for computing the total sensitivity with the naive brute-force method by approximating all sensitivities and then averaging.
Analysis. In practice, we found most real-world datasets have uus--to-aprrxximte sensitivities with total lp sensitivity about 2-5 times lower than the theoretical upper bound. Figure shows that much smaller than those
the average absolute log ratios for approximating the l, sensitivities suggested by the theoretical upper hound, even for large a. Specifically, when a = 40, we incur only 16x accuracy deterioration in exchange for 40x faster algorithm. This is significantly better than the worst-case accuracy guarantee which for p = 3 would be o =03.
More importantly, we find that empirical total sensitivities upper bounds suggest, often by a factor of at least 5. especially for p > 2. This implies that real-world structure can be utilized for improved data compression and justifies the importance of sensitivity estimation for real-world datasets. Furthermore, our novel algorithm approximates the total sensitivity up overestimate within factor of 1.3, in less than quarter of the time of the brute-force algorithm. Our observation generalizes to other real-world datasets (see Appendix E).
much lower than the theoretical
Total Sensitivity Upper Bound
Frute-Force/True
Approximation
Brute-Force Runtime (s)
Approximate Runtime (s)
14
5.2
6.4
667
105
1.5
14
11.6
14.2
673
131
2.5
27.1
13.8
14.9
693
209
3
52.4
7.2
8.9
686
192
Table 1: Runtime comparison for computing total sensitivities for the wine dataset, which has matrix shape (177, 14).
Acknowledgements
We are very grateful to: Sagi Perel, Arun Jambulapati, and Kevin Tian for helpful discussions about related work; Taisuke Yasuda for his generous help discussing results in lp sensitivity sampling: and our NeurlPS reviewers for their time, effort, and many constructive suggestions. Swati gratefully acknowledges funding for her student researcher position from Google Research and research
10

---

assistantship at UW from the NSF award CCF-1749609 (via Yin Tat Lee). Part of this work was done while D. Woodruff was at Google Research. D. Woodruff also acknowledges support from a Simons Investigator Award.

---

References
[1] Herbert Robbins and Sutton Monro. A stochastic approximation method. The annals of mathematical statistics, 1951. 1
[2] Léon Bottou and Yann Cun. Large scale online learning. Advances in neural information processing systems, 16, 2003.
[3] Tong Zhang. Solving large scale linear prediction problems using stochastic gradient descent algorithms. In Proceedings of the twenty-first international conference Machine learning. page 116, 2004.
[4] Léon Bottou. Stochastic gradient descent tricks. In Neural networks: Tricks of the trade, pages 421 121-436.30 Springer, 2012.
[5] Chih-Chung Chang and Chih-Jen Lin. LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2:27:1-27:27, 2011. Software available at http://wwet csie. ntu. edu. sw//""jli//ibbvm
[6] Nicolas Roux, Mark Schmidt, and Francis Bach. A stochastic gradient method with an exponen- tial convergence rate for finite training sets. Advances in neural information processing systems, 25, 2012.
[7] Shai loss minimization. Journal of Machine Learning Research, 14(2), 2013.
and Tong Zhang. Stochastic dual coordinate ascent methods for regularized
[8] Rie Johnson and Tong Zhang. Accelerating stochastic gradient descent using predictive variance reduction. Advances in neural information processing systems, 26, 2013.
[9] Mehrdad Mahdavi, Lijun Zhang, and Rong Jin. Mixed optimization for smooth functions. Advances in neural information processing systems, 26, 2013.
[10] Aaron Defazio, Francis Bach, and Simon Lacoste-Julien. Saga: A fast incremental gradi- method with support for non-strongly convex composite objectives. Advances in neural information processing systems, 27, 2014.
[11] Julien Mairal. Incremental scale machine learning. SIAM Journal on Optimization, 25(2), 2015.
optimization with application to large-
[12] Zeyuan Allen-Zhu and Yang Yuan. Improved svrg for non-strongly-conve or sum-of-non- convex objectives. In International conference on machine learning, 2016.
[13] Elad Hazan and Haipeng Luo. Variance-reduucd and projection-free stochastic optimization. In International Conference on Machine Learning. 2016.
[14] Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Mathematical Programming, 162(1), 2017.
[15] Michael Langberg and Leonard J Schulman. Universal e-approximators for integrals. In Proceedings of the twenty-firs annual ACM-SIAM symposium Discrete Algorithmss, 2010.
4
[16] Dan Feldman and Michael Langberg. A unified framework for approximating and clustering data. In Proceedings of the forty-third annual ACM symposium on Theory of computing, 2011.
2, 4
[17] Deanna Needell, Rachel Ward, and Nati Srebro. Stochastic gradient descent, weighted sampling, and the randomized kaczmarz algorithm. Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014. URL https: httpppp://pprcceediigggooo neurips. cc/ cc/paper_aa. paper_files// files/ .dff.
[18] Peilin Zhao and Tong Zhang. Stochastic optimization with importance sampling for regularized loss minimization. Proceedings of the 32nd International Conference on Machine Learning, pages 1-9, 2015.
12

---

[19] Angelos Katharopoulos and François Fleuret. Biased importance sampling for deep neural network training, arXiv preprint arXiv:1706. NXX::17006000443, 2017.
[20] Tyler Johnson and Carlos Guestrin. Training deep models faster with robust, approximate importance sampling. In Advances in Neural Information Processing Systems 31, 2018.
[21] Sören Mindermann, Jan M Brauner, Muhammed T Razzak, Mrinank Sharma, Andreas Kirsch, Winnie Xu, Benedikt Höltgen, Aidan N Gomez, Adrien Morisot, Sebastian Farquhar, et al. Pri- oritized training on points that are learnable, worth learning, and not yet learnt. In International Conference on Machine Learning, 2022.
[22] Vladimir Braverman, Dan Feldman, Harry Lang, Adiel Statman, and Samson Zhou. New frameworks for offline and streaming coreset constructions. arXiv preprint arXiv: 1ol:::6610008
2016.
[23] Michael B Cohen and Richard Peng. Lp row sampling by lewis weights. In Proceedings of the forty-seventh annual ACM symposium on Theory of computing, 2015. 2, 4, 5, 6, 7, 8, 28, 29
[24] David Woodruff and Taisuke Yasuda. Sharper bounds for l sensitivity sampling. In- ternational Conference on Machine Learning, pages 37238 -37722. PMLR, 2023. 2, 3,
8
[25] Raphael Meyer, Cameron Musco, Christopher Musco, David P Woodruff, and Samson Zhou. Fast regression for structured inputs. In International Conference on Learning Representations (ICLR), 2022. 2, 6
[26] Kasturi Varadarajan and Xin Xiao. On the sensitivity of shape fitting problems. In 32nd International Conference on Foundations of Software Technology and Theoretical Computer Science, 2012. 2,
[27] Vladimir Braverman, Dan Feldman, Harry Lang, Adiel Statman, and Samson Zhou. Efficient coreset constructions via sensitivity sampling. Proceedings of Machine Learning Reseurch, 157 (2021), 2021.
[28] Murad Tukan, Xuan Wu, Samson Zhou, Vladimir Braverman, and Dan Feldman. New coresets for projective clustering and applications. In International Conference on Artificial Intelligence and Statistics, 2022. 2, 4, 6
[29] Olivier Bachem, Mario Lucic, and Andreas Krause. Coresets for nonparametric estimation-the case of dp-means. In International Conference on Machine Learning, 2015. 2, 4
[30] Mario Lucic, Olivier Bachem, and Andreas Krause. Strong coresets for hard and soft bregman clustering with applications to exponential family mixtures. In Artificial intelligence and statistics, 2016. 2,
[31] Jonathan Huggins, Trevor Campbell, and Tamara Broderick. Coresets for scalable bayesian logistic regression. Advances in Neural Information Processing Systems, 29, 2016. 2, 4
[32] Alexander Munteanu, Chris Schwiegelshohn, Christian Sohler, and David Woodruff. On coresets for logistic regression. Advances in Neural Information Processing Systems, 31, 2018.
[33] Petros Drineas, Michael W Mahoney, and Shan Muthukrishnan. Sampling algorithms for lz regression and applications. In Proceedings of the seventeenth annual ACM-SIAM symposium on Discrete algorithm, 2006. 2, 4,6
[34] Michael W Mahoney al. Randomized algorithms for matrices and data. Foundations and Trends@ in Machine Learning. 3(2), 2011. 2, 4,
[35] David P Woodruff and Taisuke Yasuda. Online lewis weight sampling. In Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), 2023. 2, 6, 8, 26, 28, 29
[36] Flavio Chierichetti, Sreenivas Gollapudi, Ravi Kumar, Silvio Lattanzi, Rina Panigrahy, and David P Woodruff. Algorithms for lp low-rank approximation. In International Conference on Machine Learning, 2017. 2
13

---

[37] Emmanuel Candes and Terence Tao. Decoding by linear programming. IEEE transactions on information theory, 51())), 2005.
[38] Morteza Alamgir and Ulrike Luxburg. Phase transition in the family of p-resistances. Advances in neural information processing systems, 24, 2011. 2
[39] Jeff Calder. Consistency lipschitz learning with infinite unlabeled data and finite labeled data. SIAM Journal on Mathematics of Data Science, 1(4), 2019.
[40] Mauricio Flores, Jeff Calder, and Gilad Lerman. Analysis and algorithms for fy-based semi- supervised learning on graphs. Applied and Computational Harmonic Analysis, 60, 2022.
2
[41] Abderrahim Elmoataz, Matthieu Toutain, and Daniel Tenbrinck. On the p-laplacian and oc- laplacian on graphs with applications in image and data processing. SIAM Journal on Imaging Sciences, 8(4), 2015. 2
[42] Kenneth Clarkson, Ruosong Wang. and David Woodruff. Dimensionality reduction for tukey regression. In International Conference on Machine Learning, 2019. 3, 7,
[43] Cynthia Dwork. Aaron Roth, et al. The algorithmic foundations of differential privacy. Founda- tions and Trends@ in Theoretical Computer Science, ((3-4):211- 1-407, 2014.
[44] Deeksha Adil, Richard Peng, and Sushant Sachdeva. Fast, provably convergent irls algorithm for p-norm linear regression. Advances in Neural Information Processing Systems, 32, 2019.
6, 30
[45] Arun Jambulapati, Yang P Liu, and Aaron Sidford. Improved iteration complexities for over- constrained p-norm regression. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, 2022. 4, 6, 30
[46] Arthur Asuncion and David Newman. Uei machine learning repository, 2007. 4, 10
[47] Anant Raj. Cameron Musco, and Lester Mackey, Importance sampling via local sensitivity. In International Conference on Artificial Intelligence and Statistics, 2020. 4
[48] Alaa Maalouf, Adiel Statman, and Dan Feldman. Tight sensitivity bounds for smaller coresets. In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 2051-2066 2020. 4
[49] Ruosong Wang, Russ R Salakhutdinov, and Lin Yang. Reinforcement learning with general value function approximation: Provably efficient approach via bounded eluder dimension. Advances in Neural Information Processing Systems, 33:6123-6135 2020.
[50] Dingwen Kong, Ruslan Salakhutdinov, Ruosong Wang, and Lin Yang. Online sub-sampling for reinforcement learning with general function approximation. arXiv preprint arXiv:2106.07203, 2021.4
[51] Lucas Liebenwein, Cenk Baykal, Harry Lang, Dan Feldman, and Daniela Rus. Provable filter pruning for efficient neural networks. In International Conference on Learning Representations, 2019.
[52] Murad Tukan, Loay Mualem, and Alaa Maalouf. Pruning neural networks via coresets and convex geometry: Towards no assumptions. Advances Neural Information Processing Systems, 2022.
[53] Ben Mussay, Dan Feldman, Samson Zhou, Vladimir Braverman, and Margarita Osadchy. Data- independent structured pruning of neural networks via coresets. IEEE Transactions on Neural Networks and Learning Systems, 33(12):7829-78) 2021. 4
[54] Lewis. Finite dimensional subspaces of ly. Studia Mathematica, 1978. 4,5 5
[55] Gideon Schechtman. More on embedding subspaces of Ep in fp. Compositio Mathematica,
1987.
14

---

[56] Jean Bourgain, Joram Lindenstrauss, and Vitali Milman. Approximation of zonoids by zono- topes. 1989.
[57] Michel Talagrand. Proceedings of the American Mathematical Society, 1990.
[58] Michel Ledoux and Michel Talagrand. Probability in Banach Spaces: isoperimetry and processes, volume 23. Springer Science & Business Media, 1991.
[59] Gideon Schechtman and Artem Zvavitch. Embedding subspaces of Ep into f.. 0 p < Mathematische Nachrichten, 2001. 4
[60] Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron Sidford. Uniform sampling for matrix approximation. In Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, 2015. 4, 5, 19
[61] Petros Drineas, Malik Magdon-Ismail, Michael W Mahoney, and David P Woodruff. Fast approximation of matrix coherence and statistical leverage. The Journal of Machine Learning Research, 13(1), 2012.
[62] Mu Li, Gary L Miller, and Richard Peng. Iterative row sampling. In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, 2013. 4
[63] Christian Sohler and David P Woodruff. Strong coresets for k-median and subspace approxima- tion: Goodbye dimension. In 2018 IEEE 59th Annual Symposium Foundations of Computer Science (FOCS), 2018. 4
[64] Rachit Chhaya, Jayesh Choudhari, Anirban Dasgupta, and Supratim Shit. Streaming coresets for symmetric tensor factorization. In International Conference on Machine Learning, 2020.
[65] David Durfee, Kevin A Lai, and Saurabh Sawlani. & regression using lewis weights precondi- tioning and stochastic gradient descent. Conference On Learning Theory, 2018.
[66] Yi Li, Ruosong Wang. Lin Yang, and Hanrui Zhang. Nearly linear row sampling algorithm for quantile regression. In Proceedings of the 37th International Conference on Machine Learning, 2020. 4
[67] Vladimir Braverman, Petros Drineas, Cameron Musco, Christopher Musco, Jalaj Upadhyay, David P Woodruff, and Samson Zhou. Near optimal linear algebra in the online and sliding window models. In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), 2020. 4, 6
[68] Aditi Laddha, Yin Tat Lee, and Santosh Vempala. Strong self-concordance and sampling. In Proceedings of the 52nd annual ACM SIGACT symposium on theory of computing, 2020. 4
[69] Yi Li, Ruosong Wang. Lin Yang, and Hanrui Zhang. Nearly linear row sampling algorithm for quantile regression. In International Conference on Machine Learning, 2020. 4
[70] Arvind V Mahankali and David P Woodruff. Optimal l column subset selection and fast ptas for low rank approximation. In Proceedings of the 202/ ACM-SIAM Symposium on Discrete Algorithms (SODA), 2021. 4
[71] Tung Mai, Cameron Musco, and Anup Rao. Coresets for classfiiationn-smpliiec and strength- ened. Advances in Neural Information Processing Systems, 34, 2021. 4
[72] Xue Chen and Michal Derezinski. Query complexity of least absolute deviation regression via robust uniform convergence. In Conference on Learning Theory, 2021.
[73] Aditya Parulekar, Advait Parulekar, and Eric Price. LI regression with lewis weights sub- sampling. Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, 2021. 4
[74] Frederie G Foster. On the stochastic matrices associated with certain queuing processes. The Annals of Mathematical Statistics, 24(3), 1953, 5
15

---

[75] David P Woodruff and Taisuke Yasuda. High- ligh-dimensioall geometric streaming in polynomial space. In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 732-743, 2022. 5, 28
[76] David P Woodruff Cameron Musco, Christopher Musco, David P Woodruff, and Taisuke Yasuda. Active sampling for linear regression beyond the 12 norm. CoRR, 2021. 5
[77] Maryam Fazel, Yin Tat Lee, Swati Padmanabhan, and Aaron Sidford. Computing lewis weights to high precision. In Proceedings of the 2022 Annual ACM-SIAM Symposium Discrete Algorithms (SODA), 2022. 5, 6, 8
[78] David P Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends@ in Theoretical Computer Science, 10(1-2), 2014. 5
[79] Mark Rudelson and Roman Vershynin. Sampling from large matrices: approach through geometric functional analysis. Journal of the ACM (JACM), 2007. 6
[80] Christian Sohler and David P Woodruff. Subspace embeddings for the E 1-norrm with applica- tions. In Proceedings of the forty-third annual ACM symposium on Theory of computing, pages 755-764, 2011.
[81] David Woodruff and Qin Zhang. Subspace embeddings and & f_p-regression using exponential random variables. In Conference on Learning Theory, pages 546-567. PMLR, 2013.
[82] Cameron Musco, Christopher Musco, David P Woodruff, and Taisuke Yasuda. Active linear regression for fp norms and beyond. In 2022 JEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), 2022. 7
[83] Michael B Cohen, Yin Tat Lee, and Zhao Song. Solving linear programs in the current matrix multiplication time. Journal of the ACM (JACM), 68(1), 2021. 17
[84] Roman Vershynin. science. Cambridge university press, 2020. 17
probability: An introduction with applications in data
16

---

A Omitted proofs: general technical results
Fact A.1 (Crude sensitivity approximation via leverage scores), Given a matrix A Ruxd let o (a.) denote the eth E1 sensitivity of A with respect to A, and let (A) denote its ith leverage
score with respect to A. Then we have / [((A) S o (a,) S (A)
Proof. The proof of this claim follows from simple application of standard norm inequalities, and we present one here for completeness. For any E II"", we have ||u/|2 I ||l]i s Vn/jul2 Therefore, for any x E IRd, we have
x a.
x a.l AA||
x a. || |
(A.1)
Ax
Let x, be the vector that realizes the ith leverage score, and x, be the vector that realizes the ith & sensitivity (i.e. they are the maximizers sensitivity). Then, we may conclude from Inequality (A.1) that
X a
2 Ax-
x a: Ax,
1x a:l
x a Ax
X
a Ax,
<
<
S
1
Ax,
11
which gives the claimed result.
Fact A.2. Given an n x d matrix, the cost of computing k of ( sensitivities is (A)nz((A)))) +k)
Proof. Given an X d matrix A, we first compute SA, its l subspace embedding of size O(d) d.
and compute the 2th sensitivity of A follows.
a x Ax
max a. x|
o (a,)
O(||SA]] 1) RO(()) o (a,),
max xER-
xER-
where the second step is by the definition of l subspace embedding (cf. Definition 2.3). To compute o (SA(a,) (a:), need to compute mina,x= x=1 |SAA|| which, by introducing
variable for each of the rows of S Ax, can be transformed into linear program with O(d) variables and d constraints. To see the claim on runtime. the cost of computing the subspace embedding is Inz(A). Having computed this once, we can then solve & linear programs, each cost d (which is the current fastest LP solver [83]). Putting together these two costs gives the claimed runtime.
Lemma A.3. Given positive numbers 01, 02,..., Oaro let r max (in] Ra mins as and A (true) := Liglm) 0g- Suppose we sample, uniformly at random with replacement, a set S E (m] of these numbers, and let (en)) == IST Eia,es Og. Then, i 2 10 r(1+>) log (t) for a large enough absolute constant C, we can ensure with at least a probability of 5, that [Atea) Altre)) S A(rue))
Proof. By construction, the expected value of sample a; in S is (tue), Denoting by P the uniform distribution over d1 .2..... am and let (PP) be the standard deviation of this distribution P. Then, we can apply Bernstein's 's inequality [84, Theorem 2.8.1] to see that the absolute error [A(aat) A (wno)] satisfies the following guarantee for all 1 >>0 O:
Pr Pr(ln(ea)
A ((mee) 2 21)
Pr
Σ (a:
A ((rue)
2 t
S
kaa ES
1s
2.
m
81.0 (P) +
Si
mm
<
m
1
(P) + m t--.x 0z
e
(A.2)
17

---

Since each of the a;s is positive, we note that
a(P)² m [m]]
1
Σ
a4
A(true)
0
1
max ai 16[ml
Σ iE(m]
0:
1
max dj A
m
m
m
m
RE ml
(m(m
(A.3)
Combining Inequality (A.2) and Inequality (A.3) for (tme) implies that
2 (tme)
Pr
{1A((s)) - ((mm))) 2 7
Se
m(1+ 3 7) ma 0z
(A.4)
holds for the choice of |S] 10 [3 (1 ?)
maxieeml as log (?) yields the claimed emror guarantee. A(me)
B Omitted proofs: 4 sensitivities
B.1 Estimating all 4 sensitivities
Theorem 3.1. Given a full-rank matrix A E Waxd and an approximation factor 1 a < n. let o (as) be the & sensitivity of the 8th row of A. Then there exists algorithm that, in time Ô (nnz(A) A) + a d returns ō <Roo such that with high probability, for each i [n].
(a()) S o. 100, (a.) 6
(3.1)
Proof. Note that in Line of Algorithm 1, the rows of A are partitioned into randomly created n/a blocks. Suppose a,, the 7t row of A, falls into the block B₁ Recall, that in Line 5, the rows from B₁ are mapped to those in P with row indices in the set J. Then, compute the 7th L sensitivity estimate as o; = maxjeJ 0 (py). We observe that for all x E IR.,
||AA||| O(] Ax|l))
(B.1)
We use Equation (B.1) below to establish the claimed bounds. Let x* arg maxxeR
|a] x] Ax
be the
vector that realizes the 2th & sensitivity of A. Further, let the jth row of P be r Bg. Then, with probability of at least 1/2, we have
(2)
Bax
(6)
Bex*
(1) Ir
Bax 2 O(1)
a x'
(p,) = max XER
R
|SAx*11 = 011)
= ( (a.),
SAX||1
Ax*
A
(B.2)
where the first step is by definition of of (p,). the second is by evaluating the function being maximized at specific choice of x, the third step is by applying ||SAx|| O(IIAx]I1) from Equation (B. and the final step is by definition of x* and o (a:). For the fourth step, we use the ar + Ejfs rk, (B,x") ),l 2 la x"| with a probability
fact that |r B,x" since the vector rf has coordinates that are +1 or -1 with equal probability. By a union bound over the |J] independent rows that block B, mapped to, we establish the claimed lower bound in Inequality (B.2) with probability of at least 0.9. To show an upper bound on a (p;), we observe that
least 1/2
( Bix
| O(1) max XeRa
B AA
Ja T x
SA (p,)
max XER
max XER SAx
I O(1)
Σ
max XERd
SAAI|
|Ax
j:ay EB
(B.3)
where the second step is by Hölder inequality and the fact that the entries of r are all bounded between -1 and +1, and the third step by Equation (B.1). The final term of the above chain of inequalities may be rewritten as @(1) a EB, (A), by the definition of B Because B, is group of a rows selected uniformly at random out of n rows and contains the row a, we have:
E
Σ jaaeBB
(aa))
Σ o (a,). #
n
18

---

Therefore, Markov inequality gives us that with probability of at least 0.9, we have
(p,) S o (ay) + 6 (A)O
(B.4)
The median trick then establishes the bounds Inequality (B.2) and Inequality (B.4) with high prob- ability. The claimed runtime is obtained by the cost of constructing SA plus ((n/a) computations of l sensitivities with respect to SA, for which we may use Fact 2.4.
B.2 Estimating the sum of l sensitivities
In this section, we present a randomized algorithm that approximates the total sensitivity 6 (A), up to a m-mmliipplccativ approximation for some Y (0, 1) that is not too small. This algorithm is significantly less general than Algorithm (which holds for all 2 1) but it uses extremely simple facts about leverage scores and what think is we believe could be of potential independent interest.
recursion technique in this context, which
Theorem B.1. Given matrix A E Rexd and an approximation factor y E (0,1) 1) such that 7 Ω (+),, there exists an algorithm, which in time Ô (nnz(A) द. max (vd, })) returns positive scalar & such that, with probability of 0.99, we have
(A) S 2 (1 + (1))))
Algorithm 4 Approximating the Sum of & -e-Sensitivities
Input: Matrix E RRad and approximation factor 7 E (0, 1)
Output: A positive scalar that satisfies, with probability 0.99, that
(A) s 3 (1 (()))
1: Compute SA Ro(dxxd /pp-opppox f subspace embedding of A (cf. Definition 2.3)
2: Set p= O01 (Iog( + d))). Compute S' A, a p-approximate & subspace embedding of A
3: Construct a matrix A comprising only those rows of A with leverage scores at least 1/n,
4: Let be the output of Algorithm with inputs A, SA, S' A, n. and p
5: Return 3 = (1+?) (s (AA-AA) A
B.2.1 Proof sketch.
We achieve our guarantee via Algorithm 4 and Algorithm 5, which are based on the following key principles. The first Fact 2.2: the hth sensitivity and 7th leverage score of A E Raxd satisfy
(.(A) (a,) S (T). The second idea, derived from Bernstein's 's inequality (and formalized in Lemma A.3), is: the sum of set of positive numbers with values bounded between a and b can be approximated to a ?-multiplicative factor using (((baa)> uniformly sampled (and appropriately scaled) samples. The third principle is that ar (a,) R o (a:) when SA a constant-aprroximation l subspace embedding of A.
Based on these ideas, we first divide the rows A into B = Olog n) submatrices based on which of the B buckets [1/2, 1], [1/4, 1/2], [1/8, 1/4], etc. their leverage scores (computed with respect to A) fall into. Since each of these submatrices comprises rows with leverage scores in the range [a, 2a] for some a, per Fact 2.2, we have / & or (a,) I V2a. Therefore, by Lemma A.3, the total sensitivity of this submatrix may be approximated by the appropriately scaled total sensitivity of O(vn) of its uniformly sampled rows. Our design choice to split rows based on leverage scores is based on the fact that leverage scores (unlike sensitivities) are efficiently computable (using [60]). Computing the total sensitivity of an Öivns-sized matrix thus costs ((nnz(A) + /n- d²) under this scheme.
Applying the idea sketched above recursively further reduces the cost. Suppose, after dividing the input matrix into B submatrices M, and subsampling M, from these matrices, we try to recurse on each of these B submatrices M, A key requirement for this idea to work is that the sum of the l sensitivities of the rows of M, in isolation be close enough to the true sum of L sensitivities of
19

---

Algorithm 5 Subroutine for Approximating the Sum of & --eesittivities
Input: Matrices M
Rrxd SA Ro(d)xd S'A RO(dxx, scalars and p (0, 1) Output: A positive scalar & that satisfies, with probability 0.99, that
( 8 (1 [1++0000 ((p)) (M).
1: Set B (log D Oloog(log(1 + d))). 5 g and Ō D'
D
7 vd)
72
max
2:
if the number of rows of M is most b then
3:
Let G A (M) be the sum of l sensitivities of the rows M with respect to S'A
4: Using Fact 2.4 on M, compute & such that 6 (M) < & < (1 - ( ( (M)
5: Return 8
6:
else
M SA
7:
Construct the matrix C =
Set O p)p *log ()).
8:
Compute TC), the vector of leverage scores of C
9:
Divide प,t,, 1] into geometrically decreasing sub-intervals [1/2, 1]. [1/4, 1/2], ..... etc. for each bucket [B] created in the previous step do
11:
Construct matrix M, with those rows of M whose leverage scores fall in the zth bucket Construct matrix M, composed of r uniformly sampled (with replacement) rows of M; Compute 3, the output of Algorithm 5 with inputs M. SA, S'A, n,, and p.
12:
13:
end for
15:
Return 8 := (1 p)) [
M
M
16: end if
those rows with respect to the original matrix A. general, this not true. However, we can meet our requirement using the l subspace embedding SA of the (original) matrix A. Specifically, if we vertically concatenate a matrix M, with SA and call this matrix C, then for each row M, E M, we have of (J) ~ of (M. Further, |C| I O(vn), and so in the step subsampling the rows of M,, we choose O(n²/4) samples. Recursing this exponentially decreases the row dimension of the input matrix at each level, with (logn) n matrices at the tthh level recursion.
B.2.2 Helper lemmas for proving Theorem B.1
For our formal proof, need Definition B.2, Fact B.3, and Definition B.4.
Definition B.2 (Notation for recursion tree of Algorithm 5). Recall that B = ((og(()) as in Line I of Algorithm 5. Then we have the following notation corresponding Algorithm
We define rooted B-ary tree T corresponding to the execution of Algorithum 5.
We use Tj denote the subtree starting from the root node with all the nodes up to and including those the jth level. We denote the ith node the juh level by Ty,,y Thus, every node is uniquely specified by its level in the tree and its index within that level.
The root node To.1) at level 1 corresponds to the first call to Algorithm 5 from Algorithm 4.
For the node To. we denote by ( the input matrix to the corresponding recursive call; note that all the other inputs are global and remain the same throughout the recursive call. Analogously, use ((3,8) to denote our estimate of G A (M(..)) the sum of Eı sensitivities of rows of the matrix M(A) in the node Tú.ay, defined with respect to the sketched matrix S'A.
Fact B.3. event A happens with probability TI and event B happens with probability T2 conditioned on event A, then the probability of both events A and B happening is at least Ti T2-
Definition B.4. We say that a node N with matrix M satisfies the (p, 8) approximaatin guarantee if the output of Algorithm on M, satisfies the following guarantee for the true sum @ (M) of sensitivities of the rows in M:
$ 2 G (M) with probability at least 1 -
20

---

(B.5) where the first step is because in Line of Algorithm we discard those rows ag of A with leverage scores smaller than 1/n and there are n rows in A, so every row in the input matrix M to Algorithm 5 satisfies Tx (x( 2 Aoo: the second step is by Fact 2.2 applied ag; the third step is by the definition of the & sensitivity of ar with respect to A; the fourth step is by observing that constant-factor (with for some approximation parameter p E (0,1) and an enror probability parameter 8 E (0.
Lemma B.5 (Partitioning of M). In Algorithm 5, even though we ignore the range [0, noσ when creating the submatrices M1, MM...., Mg from M, every row in M is part of exactly one of the submatrices M2..M...... Mg. In other words, the matrix M is partitioned into the matrices M1.M2. M2........ Mg, with no row missed out.
Proof. Our proof strategy is to show that for all row indices & E [[]]] each of the leverage scores (C( satisfy (C) 2 7-20- Therefore, creating buckets of rows with the range of leverage score starting at मा (instead of 0) does not miss out any rows.
To this end, we note that for every row in the input M to Algorithm 5, we have:
(A) A
x ag Ax
ak
T S
A (ax) S o
max XER"
max xER"
30 S 3 $3VTH(C)
|x
||Cx|]1 = |SAx|| + ||Mx|] I 2 Ax|]1 |Ax|| 1 =
Ax||1 since SA is
the constant assumed 1/2) & subspace embedding of A and because M is a submatrix of A; the fifth step by the definition of l sensitivity of aı with respect to C; the final step is by Fact 2.2 applied to within C. Observing the first and last terms of this inequality chain, we have (C) Z IT for all E [[C]]. This finishes the proof of the claim.
Lemma B.6. For the recursion tree T corresponding to Algorithm 5, let the size b of the input + /dlog(d)). Then, the depth of the recursion free is less than D = 1 + (og(log(21)) (2n 2dlo)()))
matrices at the nodes be b = 50D* log()00B) 22D* log(100B)
7
72
Proof. For now, let 3 be arbitrary but less than B² S + d', where d' is the row dimension of the subspace embedding SA. Then, the size of the input matrix to the nodes of the recursion tree as we increase in depth evolves
ovd BV 3vd + BV + BV + Bvd + n....
Thus, this sequence evolves as
It+1 = Bv d' + It, with xo = d'
To see the limit of this sequence, we set =3vdd + z. and solve for this quadratic equation as
3² /Bt + 4.2d 3(3 Vd))
(B.6)
2
We set the base case to be 6* := 73(3 + vd'), constant times larger than this asymptotic limit. We now compute the depth t" at which the recursion attains the value b* is attained. To do so, we consider sequence {ve} } which overestimates the sequence {2,} } of matrix sizes for big range of t. We define {ye} } as
30 + d', and Be+1 = B
and it can be checked that for l such that yı 2 d', we have Di+1 Ft+1- Moreover, for all / > 1, -1+1/244--412-- (226 = 32-2" -t (2yo) 1/2", The limit of this sequence can easily be seen to 32. We now split our analysis in two cases based on the limits of these two sequences.
Case 32 32 d' Since the limit of the sequence {ye) 8 it is easy to see that yr is always greater d' (since y: is monotonically decreasing yo 2 (39). Thus, at will always be less than y, in this case. We will now show that ye (and thus x²)) will be less than b* in lg log log(2yo) steps. Solving the inequality for & as follows: V1 < (2yo) 4, which happens if 2 to. Thus, t* < 2 log log 2yo-
< 482 < b*, it suffices to ensure that
21

---

Case 2: B <d d'. In this case, the sequence yr will eventually go below d' and thus the inequality Yr 2 zt might eventually break down. Let = t' denote the first time step at which y, S 5d' (which will be finite in this case). Then, by definition, Yt 2 Et for all 1. < L'. We will now upper bound the value of t' as follows: Y/ S (2 < d' (2 21/2" < 4d'; Thus, it suffices to ensure that 1/20 I 4, which happens if t 2 to = 2log (2yo). Hence, t' S to. Moreover, the value of be S 5d', implies that - S WV-1 S (5d')
23th We will now show that Ju+1 S b using the evolution of 2++1 Byd + I₁ as follows: First, IV S 3Vdd + 3 / S Bvd' + 5d' - 6d"; Then, 21+11 I 3Vd + /d6d' < 7 < b". Thus, S 1 2og log 2yo-
Therefore, in both of these cases, we have that after at most D = 1 + log log(2n + 2d') depth, the base case of b" = 7B(3 + vd') is achieved at the leaf node. Recall that the al- gorithm sets D = ln d') + 1 and p Y/D for d' dlog(d) which implies that CD (D( log((003*) I 2D Iog (100)) Since Algorithm 5 sets 12D (000))
B = C log(8-1)
7"
b = 50D log(100B)
+ logg(d)): which is larger than b*, we may conclude that the choice of b implies that the recursion tree T of Algorithm has depth less than
73
D = 14 log(log(2n + lolo)())).
Lemma B.7. Suppose the rool node T(1.1) of Algorithm 5 satisfies (p. 3) -ppproxiimaaionn guarantee as defined in Definition B.4, and denote by $ the output of Algorithm 4. Then, with a probability of at least 1 8, we have
I+y 5+o p+p GA
As a result of this lemma, in the final proof correctness of Algorithm 4 (appearing later as proof of Theorem B. it suffices to show that & satisfies a (p, -paaproxxiiiiooonn guarantee for some appropriate choices of these parameters.
Proof of Lemma B.7. Suppose the root node satisfies Definition B.4 with some (7, 8) approximmaio guarantee. This means that there exists an ž such that the output of Algorithm 5 on the input matrix A satisfies
6 A (A) with a probability at least 1-3. 8.
(B.7)
Since S'A is, by design, a papprooimmation & subspace embedding for A, it means for all x E IRd we have
(1 - []])]|Ax|]] < ||a s (1 p)]
(B.8)
Therefore, for every row a, E Å, have by the definition of É sensitivity and Inequality (B.8) that the & sensitivity of a; computed with respect to S' A is a p-approximation of the & sensitivity of a; with respect to A:
0 S'A (a,) = max
x
a x A
(a:).
max X
XER
S Ax|
Therefore, by linearity, this p-approximation factor transfers over in connecting the total sensitivity of A with respect to S'A to the total sensitivity of A with respect to A:
(A)
Σ S' (a:) Ep Σ -f
(B.9)
AE[[A]]
Therefore, chaining Equation (B.7) and Equation (B.9) yields the following approximation guarantee:
Sotp G (A) with a probability at least 1-8. 8.
(B.10)
There is no change in the error probability above from Equation (B.7) because Equation (B.9) holds deterministically. Next, since in Line 3 of Algorithm 4, we dropped rows of that have leverage scores less than n. for each of the dropped rows a,, we have
A S V (A) S 5
(B.11)
22

---

Therefore, the quantity we return, 3, satisfies the following approximation guarantee:
1+7 8 = 3 7 (|A]
A 2 (1 p) + Σ (1-
ia,
(B.12)
where the first step is by definition of 8 in Algorithm 4; the second step is by Equation (B.10); the third step is by applying Inequality (B.11) each of the rows not present in A. In the other direction,
1 - - n (IA) 1+7
N/4 n (1+ (p p))
A = (1 + (0 + ())
(B.13)
where the first step is by definition of 3 in Algorithm 4; the second step is by the fact that A S n; the final step is by combining the assumption 2 + and the facts that all sensitivities are non-negative and that the total é sensitivity is least one. Combining Inequality (B.12) and Inequality (B.13) finishes the proof.
Lemma B.8 (Single-Level Computation Recursion Tree). Consider a node N in the recursion tree corresponding to Algorithm S, and let M be the input matrix for N. Further assume that N is not a leaf node. Denote by N1,N22..... .BB (where B = (log(n)) is the number of recursive calls at each node in Algorithm 5) the children of node N, with each node N; with the input matrix M, (as constructed in Algorithm 5). Suppose each of these child nodes N; satisfies appprxxxxinnnaoo guarantee (cf. Definition B.4). Then N satisfies (p+ + p. (8 + 6) apppprooxinttioon guarantee.
Proof. We first reiterate the steps of Algorithm 5 that are crucial to this proof. As proved in Lemma B.5, Algorithm partitions the rows of the matrix M in node N into submatrices M1, M2,.... Mg based on their leverage scores. Next, for each matrix M,, we uniformly sample with replacement O( C[v|C((1 + 1og(6-1) rows and term this set of rows matrix M, which we then pass to Algorithm 5 recursively in node N₁ (recall that we assume that the node N is not a leaf node). In other words, each child node N; has as its input matrix M,.
Since we assume that all the child nodes N2..N....... NB satisfy a (p. 8) apprroxiimtio guarantee as in Definition B.4, it means that, for all i E [B], the output & of Algorithm on node M, satisfies the following guarantee:
S up G (M,) with a probability at least 1 3.
(B.14)
Next, we recall that each M, is composed of only those rows of M whose leverage scores (computed with respect to C) lie in the ith bucket [2 2,2-++1] Combining this with Fact 2.2, this implies that the l sensitivity of the jth row M, of the matrix M,, computed with respect to C, satisfies
2 i C
C (M,)
|C| S of ((MM) S (MM.) S V2-141
(B.15)
Additionally, we claim that
1
(M.())
S0 (M, 33 + o)
(B.16)
3(1 + P)
To see this, we observe that
[x M.L [|x]|]
x Mil
x M,
(M)
C (M,Li) = max XER
max XER
3 Ax|11
XER 3(1 + p) S'AA||
max
3(1 p)
(B.
where the first step is by the definition of l sensitivity of [.] (i.c., the 7th row of matrix M,) with respect to C: the second step is by the definition of C. the fact that SA is 1/2-factor L
subspace embedding of A, and because M is subset of rows of A; the third step is because S' A is a p-approximate & subspace embedding for A, and the final step is by the definition of a (M, (]).
In the other direction, we have
M.LiIl
MIIl
(1 p) M. S 33(+p)o S' (M.)),
of ((,
XER4
max
Cx||11
max XER 0.5
Ax|11
max XER-
S.5|S'AA
(B.18)
23

---

where the second step is by the definition of C as a vertical concatenation of SA and M and further using that SA a constant factor E subspace embedding of A and dropping the non-negative term |[Mx]| 11- Combining Inequality (B.15) and Inequality (B. (B.16), we see that for all rows j E [M we have
2 C
S Sa S'A
(M. S 3(1 + p)V2-
(B.19)
3(1
Since Inequality (B.19) shows upper and lower bounds on each sensitivity of A (M) of the matrix M. we may use Lemma A.3 to approximate the sum of these sensivities. In particular, by Lemma A.3,
we may construct a matrix M, composed of 201//01(1 log(ō-") rows of M, sampled uniformly at random with replacement, and are guaranteed
Mi
(1
G A (M,) =, A with probability at least - 8.
(B.20)
M,
Note that the output 3 of Algorithm 5 when applied to node N is defined as
M
x :1
8,-
(B.21)
M
Therefore, by union bound over the failure probability (ef. Fact B.3), we can combine Equation (B.14) and sum Equation (B.20) over all i B child nodes to conclude that 3 from Equation (B.21) satisfies (p p. (6 + 8 8)B) approximation guarantec of Definition B.4.
B.2.3 Proof of Theorem B.1
Proof of Theorem B. We first sketch our strategy to prove Theorem B.1's correctness guarantee. Essentially, our goal is establish Definition B.4 for the root node TT.1). That this goal suffices to prove the theorem is justified in Lemma B.7. We now show that the node T(1,)) satisfies (p, approxxiinnn guarantee for some P and 8. We use the method of induction.
Correctness guarantee of Theorem B.1. We know that the leaf nodes all satisfy the (p, 0)- approximation guarantee in Definition B.4 since at this level, the base case is triggered, where compute p-approximation estimate of the total f sensitivity with respect to the subspace embedding S' A. Having established this approximation guarantee at the leaf nodes, we can use Lemma B.8 to inductively propagate the property of (p. approxipppoinitiii guarantee up the recursion tree T. From Lemma B.6, we have that the recursion tree T of Algorithm is of depth at most D Olog log(n d)), and the number of child nodes at any given node is at most B (llg(n))) we factor these into the approximation guarantee and error probability.
Base case. At the leaf nodes of the recursion tree, we directly compute the total sensitivities with respect to S' A. Recall that we denote the depth of the tree by D. Then, because S' A is, by definition, p-approximate subspace embedding for A, we have, for each i that indexes a leaf node,
((D.)) P G A (M().)) with probability of 1.
(B.22)
Induction proof. Consider the subtree T, truncated (but including) the nodes level j in the recursion tree T corresponding to Algorithm's 5. Assume the induction hypothesis that all the leaf nodes
in the subtree T, satisfy Definition B.4 with parameters (p (D + 1- -) [BD-i+1_ 8-1
That
is, at the matrix of each such node i, we have an estimate &..) of its total sensitivity G (M such that
(.) (D(+11) G5'A (M(.))))) with probability at least 1 8 (BD-jt1 - B
B-1 (8.(5
(B.23)
Then, each of the leaf nodes of the subtree T (i.e., the subtree that stops one level above T₁) satisfies the assumption in Lemma B.8 since its child nodes are either leaf nodes of the entire recursion tree T (for which this statement has been shown in Incquality (B.22)) if they are not
24

---

leaf nodes of the recursion tree, they satisfy Definition B.4 with parameters p. p and
ő- 8 [=--1+1- E
B-1 Therefore, Lemma B.8 implies that each of the leaf nodes of 7 satisfies
Definition B.4 with approximation parameter p and error probability 8 defined as follows:
- - (=p((D++1))) = p- and 8 B- 26 BBD-++1 B-1 (8---11) )) = 8-
B (8D-j+2_ j+2 B
B B-1
In other words, the outputs (())1,)) of Algorithm 5 the matrices M(-1,k) in the leaf nodes of Tj each satisfy:
B ((88)
(BDD)+2 - B
F)-1K) (D(2(D42-)) 6 ((((M--1A)) with probability at least 8.
(B.24)
The base case in Inequality (B.22) and the satisfaction of the induction hypothesis in Inequality (B.24) together finish the proof of the induction hypothesis.
Finishing the proof of correctness. is at the level j = 1, we may plug in j = 1 in Inequality (B.23) and obtain that the output (1,1) of Algorithm at the root node satisfies
light of the above conclusion, for the root node Tq. which
[DD (BD-1 B -
(1.1) FO(D)) 6 (M(L1) with a probability of at least 5
(B.25)
In Algorithm 5, we choose the parameters 5 = 0.01 and p= ( Plugging these parameters g
into Inequality (B.25), we conclude that the root node satisfies Definition B.4 with 07.0. approximation guarantee. This finishes the proof of correctness of the claim.
Proof of runtime. Let b denote the maximum number of rows of matrix at a leaf node of T. Further, recall from Fact 2.4 that the cost of computing one & sensitivity of an n x d matrix is,
L = Ó ()) + d²)).
(B.26)
Then as stated in Line 1, the algorithm computes the p-approximate total sensitivity at the leaf nodes; Fact 2.4, this incurs a cost of
c(b) = Ö (A) d²).
(B.27)
At any other node, the computational cost is the work done at that level plus the total cost of the recursive calls at that level. Let C(r) denote the runtime of Algorithm 5 when the input is an r x d matrix. Then,
{
r r>b b B.C (r(L+t) log(6-1)) L C(b)
C(r) (m)=
r b b
where b is the number of rows in the base case, and L is the cost of computing leverage scores and equals = Ó(nz)( A) d²)) from Fact 2.1. For the depth of the recursion D (where the top level is 0.01 8D and p O (?), this may be expressed via expanding
D 1), plugging in our choices of 8 the recursion as:
((m) BD C(b) + L. (BD -
(B.28)
B-11
for 2 b. We plug into this expression the values of b and D from Lemma B.6, L from Equa- tion (B.26), C(b) from Equation (B.27), and B = (log(()) to get the claimed runtime.
B.3 Estimating the maximum & sensitivity
Theorem 3.7 (Approximating the Maximum of El Sensitivities). Given matrix A E Brxdd there exists an algorithum, which in time Ó(nnz(A) + d=t1), outputs a positive sealar that satisfies
(A)](() 3 ((valo
25

---

Proof. First, we have
|S |x||1 = ([I
(B.29)
We Equation (B.29) establish the claimed bounds below. First notation. Define x* and ag. as follows:
a x |Ax||
max SERAMEIIA)
(B.30)
Thus, x* is the vector that realizes the maximum sensitivity of the matrix A, and the row a,- the row of A with maximum sensitivity with respect to A. Suppose the matrix S contains the row ap.. Then we have
c x |
c x Ax
la! x AA|| (B.31)
max 01 (c())== = E
max
@(1)
max XER A
O(1)
max XER
XER CC
where the first step is by definition of & sensitivity of S. A with respect to S, the second step is by Equation (B.29), and the third step by noting that the matrix A A subset of the rows of A (which includes a.)). ). By definition of o, (a(.) in Equation (B.30), we have
a x
max XER
1,00,, (A)lle:.
(B.32)
Ax
Then, combining Equation (B.31) and Incquality (B.32) gives the guarantee in this case. In the other case, suppose a,- is not included in S. Then observe that the upper bound from the preceding inequalities still holds. For the lower bound, we observe that
| St
(ss A)l
X
Ax ((1) max xERd
||S Ax
max XER", E8
A
max XER-
||
|Ax
(B.33)
where the the second step is by choosing specific vector in the numerator and the third step uses |S||x|| Ax|]lx) We further have,
| Ax||1
120
AAx 1la
la,.x*' Vxl|Ax*
max XER
2
| 112
vd la lo, (a)
(B.34)
where the first step is by choosing x x", the second step is by the distortion guarantee of Loo. subspace embedding [35], and the final step is by definition of az (a.-) Combining Inequality (C.10)
and Inequality (C.11) gives the claimed lower bound on ||o (5:AASS.. ASloo The runtime follows from the computational cost and row dimension of S from [35] and the cost of computing & sensitivities with respect to S per Fact 2.4.
C Omitted proofs: lp sensitivities
We recall the following notation that we use for stating our results in this setting.
Definition C.1 (Notation for b Results). We introduce the notation LP(m, d, p) to denote the cost of approximating one Ep sensitivity of an m x d matrix up to accuracy of given constant factor.
C.1 Estimating all tp sensitivities
We generalize Algorithm to the case 2 1 below.
Theorem C.2 (Approximating all lp sensitivities). Given a full full-rank matrix A E R an approxi- mation factor < a & n, and a scalar p 2 1, let (a(a)) ) be the tth Ep sensitivity. Then there exists an algorithu that returns a vector ō E 1S such that with high probability, for each E [n], we have
(a.)(a))) < o, < O(a (a())) + =H (A)).
(C.1)
Our algorithm rons in time δ (nnz A) n _ LP( d,p)) (cf. Definition 1.2).
26

---

Algorithm 6 Approximating <p-Sensitivities: Row-wise Approximation
Inputs: Matrix A E Raxd approximation factor a E (1, nl,, scalar p 2 1 Output: Vector o R that satisfies, for each E [n], with probability 0.9, that
o pp(ai) s So; S 0a o (a.) + =" (A))
1: Compute, for A, an lp subspace embedding S E IR (ddduua(t.)/3))
2: Partition A into 2 blocks ......, Bnfa each comprising o randomly selected rows.
3: for the fth block B,, with & E do
4: Sample a-dimensional independent Rademacher vectors ri (0) (f) ( 5: For each jE(100), compute the row vectors r B, E Idd
1100
6: end for
7: Let P E IR. 12 xd be the matrix all vectors computed in Line 5. Compute o S,A (P) using
Definition 1.2
8: for each i E [T]] do
9: Denote by J the of row indices in P that a, is mapped to in Line 5
10: Set ã; max A S (p;))
11: end for
12: Return
Proof. Suppose the 7th row of A falls into the bucket BN. Suppose, further, that the from B, mapped to those in P with row indices in the set J. Then, Algorithm 6 returns a vector sensitivity estimates E IR" with the 7th coordinate defined o; maxjeJ (C()). For any Rd, have
SSpAA| ~ Ax
(C.2)
la x|
We Equation (C.2) establish the claimed bounds below. Let arg maXxeRs
Txx be the
vector that realizes the üth lp sensitivity of A. Further, let the jth row of be r Be Then, with probability of at least 1/2, we have
S. = max XER
Ir B
F
B
O(1)
F B,x 1 2 O(1)
a x Ax 1
((1)or,(a))
Ax
AS
Ax
(C.3)
where the first step is by definition of o (P,), the second is by evaluating the function being maximized at a specific choice of x, the third step is by Equation (C.2), and the final step is by definition of x* and (a.). For the fourth step, we use the fact that 1rk B2x+1 a, x
Ej k. (Bx*), that are -1 or +1 with equal probability. By union bound over B, is mapped to, we establish the claimed lower bound in Incquality (C.3) with probability at least 0.9. To show upper bound o S,A, (p;), we observe that
a, x" PP with probability of at least 1/2 since the vector r
(f) has coordinates independent rows that block
p
r B
B P Ax B
T a x P
max xE
|S,AA||
5Qp-1
max xcRd
S O(a
Σ ca
max xERd
|Ax|E
O(ap-1)
Σ (,))),
jayEB
(C.4)
where the second step is by Equation (C.2) and opening up |Bxx| [ in terms of the rows in Bg, and the final step is by definition of Br. To see the first step, we observe that
( k B < |Bx] ||
(2) 1 < B. "/1
B "-1
where 1 + 1 = 1, the first step is by Hölder's inequality, the second is because each entry of ry is p q
either +1 or - and (2) Ra, and the third is because ạ p 1. Because B, is group of rows
27

---

selected uniformly at random out of rows and contains the row we have:
) } Eo
E
E
opla,)
jzti
Therefore, Markov inequality gives us that with probability of at least 0.9, have
(p,) & 10aa"'lppp0 p((a,) - (A()).
(C.5)
The median trick then establishes Inequality (C.4) and Inequality (C.5) with high probability. To establish the runtime, note that first construct the subspace embedding S,A and then compute the Ep sensitivities of P (with (n/aa) rows) with respect to S,A.A The size of the subspace embedding Sp..A is vp/2 x d) (see [75, Table 1J). This completes the runtime claim.
C.2 Estimating the maximum of sensitivities
In Section 3.2 and Appendix B.3, we showed our result for estimating the maximum & sensitivity, In this section, show how to extend this result to all p 2 1. Algorithm 7 generalizes Algorithm 3.
Algorithm 7 Approximating the Maximum of l Sensitivviiss
Input: Matrix A E Raad and 2 (with # 2) Output: Scalar l2z that satisfies that
(A)ll S C' C.KC../.... (|o)))) (3
1: Compute, for A, an too subspace embedding S A E OROld (d)) such that S A is a subset of the rows of A [35]
Roldixd
p E (1,2)
2: Compute, for A, an bp subspace embedding S,AC
Rold
p> 2
3: Return d 2 SpA (S
Theorem C.3 (Approximating the maximum of ép sensitivities). Given matrix A E guxd and p 2 1 (with? & 2), there exists algorithm, which oulputs positive scalar 8 that satisfies
[[Io,(())ll= S (1p/2110 S (A))))
The runtime of the algorithm is Ó(nzz A) + d=-1) for p E [1,2) and ((nzz(A) + dP/2 LP(O(a)13) y,ddp"",ddd.)..... d)) for > 2.
Proof. First, have
|S + (3).
(C.6)
We use Equation (C.6) establish the claimed bounds below. First we some notation. Define
and a, as follows:
a x P
1 arg
(C.7)
XERMHEIIA
| E
Thus, x" is the vector that realizes the maximum sensitivity of the matrix A, and the row a.- the row of A with maximum sensitivity with respect to A. Suppose the matrix A contains the row a... Then we have
max ia..S
Sp (c.)= S
=
XER
max C.ES
S
C x P
Ax|
= O(1)
max CLES
c x P A
O(1)
max XER
a x
| 7 (C.8)
P
FFor p = 2, the result of [23] gives constant factor approximation to all leverage scores in mnz + d time.
28

---

where the first step is by definition of Ep sensitivity of 2 with respect to SpA, the second step by Equation (C.6), and the third step by noting that the matrix x is subset of the rows of A (which includes a.+). By definition of o (a:)) in Equation (C.7), we have
a x
max RER
= 10 (A)
(C.9)
Ax
Then, combining Equation (C.8) and Equation (C.9) gives the guarantee in this case. In the other case, suppose a,. is not included in A. Then we observe that the upper bound from the preceding inequalities still holds. For the lower bound, we observe that
1x,x/P |
SS & |S
20 Ax
| (S =
max CEEE
max XER
((1)
max XER"
|x
(C.10)
where the the second step is by choosing specific vector in the numerator and the third step IS A = -x] []) We further have,
Ax
|S Ax"
1
AA 1
a.xx de/2 Ax do/2 11 (A)l (C.11)
1
max XER
Ax
|A 1B
dP 1Ax
where the first step is by choosing = x", the second step is by the guarantee of Eoo subspace embedding, and the final step is by definition of o, (a,-). Combining Inequality (C.10) and Inequal-
ity (C.11) gives the claimed lower bound on || aS Aooo The runtime follows from the cost of computing S A from [35] and the cost of computing and size of S, from [23, Figure I] O
D Lower Bounds
Theorem D.1 (lp Regression Reduces to Ep Sensitivities). Suppose that we are given an algorithm
A, which for any matrix A' E Tin' xd' and accuracy parameter e' E (0, 1), computes (1 + E')a (A,) in time T(n', d', mnz( A'),e' Then, there exists algorithm that lakes A E Raxdd and b R" as imputs and computes (1 8) minyeRs ||Ay- - bbl + XP in time T (n+ d+ 1, nnz A) nnz(())nnz(bb nb) 1, E) for any X > 0.
A 0 -
b
Proof. Given A E Raxd and b E Rd, consider the matrix A'
E (n-1))dd Then
the n + 1-th lp sensitivity of A' is
Kea x) A''||
op (2(8/411)
max xER
max x
|A' 1 minyer--1 {y - 2 bb
where the second step is by the scale-invariance of the definition of sensitivity with respect the variable of optimization. By scale-invariance we see that for any e >> O, and miny Iy -cay
d² miny AA_.y -4 lb.
Therefore, note that rearranging gives minyER= ||Ay bil a (0 XP, which is how
A"
can estimate the cost of the regression. Then computing o o, (a/41)) to multiplicative accuracy of e
gives the value of minyer Ayy - blβ up to error e( (( s(>² minyeze-1 ||Ay - bbll)) giving our final bounds.
2"
Using the same technique as in Theorem D.1, we can similarly extend the reduction to a multiple regression task. In particular, we show that computing the values of a family of certain regularized leave-one-out regression problems for matrix may be obtained by simply computing leverage scores of an associated matrix; therefore, this observation demonstrates that finding fast algorithms for sen- sitivities is as hard as multiple regression tasks. Specifically, for some matrix A IRaxd we denote A.-; E xnx as the submatrix of A with its i-th column, denoted as A removed. We show that approximate sensitivity calculations can solve minyeRs AA,, + AAgllp 2ll approximately for all i,
29

---

Lemma D.2 (Regularized Leave-One-Out Ep Multiregression Reduces to Ep Sensitivities). Sup- given algorithm A, which for any matrix A' E R and accuracy parameter Given a matrix A E R
pose that E (0, 1), computes (1 + ,')),((A))))))=) time T(n', ,,', .nnn(AA),
with n 2 ler OPT; == minyeaa-1 ||A A and y; == arg mig - + A.ill for all the i E [d]. Then, there exists an algorithm that takes A Raxd and computes s1OPT + (1 + (Iy) in time T(n d, d, unz(A)e)) f)r any O.
A
Proof. Given A E Raxd consider the matrix A' =
AI E Ri++d)xd formed by vertically
appending a scaled identity matrix to A, with 0. By the definition sensitivities, we have the following.
X a t
of (a
XER
|
XER x
max
|
max XER*
A
minyeRd-1 {1 Aelly 1b
Recall that y; == arg min y AA + Asllp: and since the power is monotone transform,
OPT; == y A Aillp = AA, A.lpp. By icale-inv we see that for any
<> e> cy; = arg minlyeer- AA:.y -y + .s. [1 and miny |A..yy + cA. = eOPT, Therefore,
by plugging y = 2 X note that
min {1 + Xyyl +
A.
Iy 1
However, we also note that by non-negativity.
min yERd-1 {1 + llyly| |AA.
+ POPT,
Combining this all together gives the claimed approximation guarantee by noting that
AP
OPT, + X
XP
OPT,
a a
Therefore, our algorithm is to simply return (p where is the e-approximate estimate of the sensitivity and bound follows. The runtime guarantees follow immediately from assumption that the runtime of calculating the é, sensitivities of A'.
Corollary D.3. Assuming that leave-one-our lp multi-regression with (dd) instances takes poly(d) nnz A) + poly(1/s) time [44, 45], computing Ep sensitivities of an n x d matrix A costs at least (oly(d) mnz(A) + poly(1/E).
E
Additional Experiments
Here include additional experiments other similar datasets.
Total Sensitivity Upper Bound
Brute-Force
Approximation
Brute-Force Runtime (s)
Approximate Runtime (s)
11
4.7
3.9
1940
303
1.5
11
8.3
10.3
1980
280
2.5
20
10.3
11.8
1970
326
36.4
6.9
7.3
1970
371
Table 2: Runtime comparison computing sensitivities the datasct, which shape (617, 11). We include the thenretical upper hound for the total sensitivities using lewis weights calculations.
30

---

Total Sensitivity Upper Bound
Brute-Foree
Approximation
Brute-Force Runtime (8)
Approximate Runtime (S)
11
4.1
S.2
540
154
1.5
11
10.1
6.7
560
201
20
12.2
390
117
36.4
8.8
354
136
Table 3: Runtime comparison for computing total sensitivities for the concrete dataset, which has matrix shape (101, 11). We include the theoretical upper bound for the total sensitivities using lewis weights calculations.
Sensitivity Approximations for Fires Dataset
p=
p = 1.5
additvee bound
additive bound max
max
Absolute Log Hallo
6
5
Absoolute
5
15
20
25
30
35
40
10
20
30
35
40
Approximation (alpha)
Approximation Factor (alpha) p 3
2.5
additive bound average
max
12
Absollute Log Ratto
additive bound
average
15
40
Approximation Factor (alpha)
Approximation (alpha)
Figure 2: Average absolute log ratios for all fp sensitivity approximations for fires.
31


Reviews:
[{"review_id": "p47pZY40K0", "reviewer": "Reviewer_Y5nz", "summary": "Given a matrix $A \\in \\mathbb{R}^{n \\times d}$ and $p \\in (0, \\infty)$, the $l_p$ sensitivity of a vector $a \\in \\mathbb{R}^d$ with respect to $A$ is defined as $\\sigma_p(a) := \\max_x |a^\\top x|^p / |Ax|_p^p$. It is known that by sampling each row $a_i$ of $A$ with probability proportional to (an upper bound on) $\\sigma_p(a_i)$, we can obtain a coreset of $A$ with respect to the $l_p$ loss. Therefore, estimating $\\sigma_p(a_i)$ quickly is desirable. This work presents three fast approximation algorithms for $l_p$ sensitivities. The main focus is on the case where $p = 1$, so I will state the results specifically for that case:\n\n- An algorithm that provides an estimate $\\tilde{\\sigma}$ such that $\\sigma_1(a_i) \\leq \\tilde{\\sigma} \\leq \\sigma_1(a_i) + \\alpha/n \\mathfrak{S}_1(A)$, where $\\mathfrak{S}_1(A)$ is the sum of $l_1$ sensitivities. The algorithm's running time is $O(n/\\alpha (nnz(A) + d^\\omega) + n)$, where $\\omega$ is the matrix multiplication exponent.\n- An algorithm that provides an estimate $\\tilde{\\sigma}$ such that $\\mathfrak{S}_1(A) \\leq \\tilde{\\sigma} \\leq (1 + \\gamma)\\mathfrak{S}_1(A)$. The running time of this algorithm is $O(\\sqrt{d} (nnz(A) + d^\\omega))$, significantly faster than the naive bound of $O(n (nnz(A) + d^\\omega))$.\n- An algorithm that provides an estimate $\\tilde{\\sigma}$ such that $\\max_i \\sigma_1(a_i) \\leq \\tilde{\\sigma} \\leq \\sqrt{d} \\max_i \\sigma_1(a_i)$. The running time for this algorithm is $O(d(nnz(A) + d^\\omega))$, also faster than the naive bound.\n\nThe authors showed that $l_p$ regression reduces to $l_p$ sensitivity calculation, implying that designing a fast algorithm for $l_p$ sensitivity requires a fast algorithm for $l_p$ regression. Experimental results are provided to validate the theoretical bounds.", "questions": "- This is just a comment, but \"additive error of $O(\\alpha^p)$\" should be \"additive error of $\\alpha^p/n \\cdot \\mathfrak{S}_p(A)$\".\n- What is the meaning of the \"$\\rho$-factor subspace embedding\" in Algorithm 2? Should it be \"$\\rho$-approximate subspace embedding\"?\n", "limitations": "The limitations are not clearly mentioned.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The algorithm effectively utilizes $l_p$ subspace embedding to enhance its speed. Additionally, the recursive approach employed for computing total sensitivity is intriguing.", "weaknesses": "- There are reservations about the usefulness of Theorem 1.2. It is possible to estimate sensitivities in $O(n (nnz(A) + d^\\omega))$ time, as indicated in Fact 2.4. This means that, to achieve a significant speedup, $\\alpha$ needs to be a function of $n$, which then results in a substantial additive error.\n- The usefulness of Theorem 1.3 is unclear. Although the sum provides a bound on the coreset's size, it is not evident when one would solely desire the bound without the coreset itself.\n- Similarly, the usefulness of Theorem 1.4 is not clear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are deeply grateful to the reviewer for their careful reading of our manuscript and for raising very thoughtful questions, which we respond to below. \n\n-----------------\n\n### ”Usefulness of Theorem 1.2”: Runtime ###\n**We would like to clarify that our actual runtime is significantly better than what we (mistakenly) first stated in the main submission.** Our sub-optimal runtimes existed only in our theorems, and the algorithms and analyses originally presented indeed give the correct runtimes (e.g., cf. Appendix of the Supplementary file, Lines 551 - 575). Specifically, our correct main result (for the case $p=1$, cf. Lines 73-74 of the Supplementary Material) is that we can approximate all $\\ell_1$ sensitivities at a cost of $$ O\\left(\\text{nnz}(A) + \\frac{n}{\\alpha} \\cdot d^\\omega\\right), $$ with an approximation guarantee $\\sigma_1(a_i) \\leq \\widetilde{\\sigma}_1(a_i)\\leq \\frac{\\alpha}{n} \\mathfrak{S}_1(A)$ for all $i\\in [n]$. As the reviewer noted in their review, naively computing all sensitivities would cost $O(n \\cdot (\\text{nnz}(A) + d^\\omega))$. Thus, **our obtained runtime greatly improves over the trivial and (to our best knowledge) only known result**. Our proof sketch for this runtime uses the construction of $\\ell_1$ subspace embeddings [4] and reduction to solving a $d\\times d$ linear program (and applying [5] to solve it). Please see the top-level response for the details. \n\n-----------------\n\n### \"Usefulness of Theorem 1.2\": Sample Complexity ### \nFor the case $p > 2$, by using Theorem $1.5$ of [1], we have that $$\\text{the sample complexity with approximate sensitivities} = O\\left(\\alpha^{2p} \\cdot \\mathfrak{S}_p^{2-2/p}(A)\\right),$$ as opposed to $$\\text{the sample complexity with true sensitivities} = O\\left(\\mathfrak{S}_p^{2-2/p}(A)\\right),$$ and $$\\text{the sample complexity with Lewis weights} = O\\left(d^{p/2}\\right).$$ Assume that $p>2$ is large and the total sensitivity $\\mathfrak{S}_p(A)$ is small (say, $\\mathfrak{S}_p(A) = d$). Further assume, as a toy example, $n = d^{10}$ and $\\alpha = n^{\\frac{1}{10 p}} = d^{\\frac{1}{p}}$. Then, our approximate sensitivities give a sample complexity of $O(d^4)$, true sensitivities give a sample complexity of $O(d^2)$, and Lewis weights sampling gives $O(d^{p/2})$. Thus, our approximate sensitivities preserve the regression approximation guarantee, while **increasing the total sample complexity by only a small $\\text{poly}(d)$ compared to that with true sensitivities (while still being much smaller than that with Lewis weights) and incurring a much smaller cost.**\n\n-----------------\n\n### ”Usefulness of Theorem $1.3$”. ### \n Please refer to our top-level response for the full answer; Briefly, one potential application of Theorem $1.3$ is to run the corresponding algorithm to determine the total sensitivity and, only if this quantity is substantially smaller than expected, we may use sensitivity sampling; this is based on the result by Woodruff and Yasuda [1] that in cases with low total sensitivity, sensitivity sampling significantly outperforms Lewis weights sampling. \n\n-----------------\n\n### ”Usefulness of Theorem 1.4”. ###\n The maximum sensitivity **captures the importance of the most important datapoint** and finds applications in, e.g., experiment design to detect the most important features and in reweighting matrices for low coherence, which was studied for $p=2$ in [6]. Additionally, it captures the maximum extent to which a datapoint can influence the objective function, thus finding applications in differential privacy [3]. \n\n-----------------\n\n### Clarifying notation. ###\n Yes, $\\rho$-factor subspace embedding means a $\\rho$-approximate subspace embedding, i.e., a matrix $S$ such that for all $x$, we have $\\|SAx\\|_p^p \\in (1\\pm \\rho) \\|Ax\\|_p^p$. We will include this in the notation section and propagate the change wherever applicable. \n\n-----------------\n\n### Conclusion. ### \nWe again thank the reviewer for their time and effort in bringing up very pertinent questions; we will incorporate these answers in our manuscript. \n\nPlease let us know if any questions remain, and if we answered all the questions, we’d like to respectfully request that the reviewer re-consider their score of our submission. We are happy to answer any further questions. \n\n-----------------\n\n### References ### \n\n[1] “Sharper Bounds for $\\ell_p$ Sensitivity Sampling”, David P. Woodruff and Taisuke Yasuda, ICML 2023\n\n[2] “Uniform Sampling for Matrix Approximation”, Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, Aaron Sidford, ITCS 2015\n\n[3] “The Algorithmic Foundations of Differential Privacy”, Cynthia Dwork and Aaron Roth, Foundations and Trends in Theoretical Computer Science 2014\n\n[4] \"$\\ell_p$ Row Sampling by Lewis Weights\", Michael B. Cohen and Richard Peng, STOC 2015 \n\n[5] \"Solving Linear Programs in the Current Matrix Multiplication Time\", Michael B. Cohen, Yin Tat Lee, and Zhao Song, STOC 2019\n\n[6] \"Dimensionality Reduction for Tukey Regression\", Kenneth L. Clarkson, Ruosong Wang, and David P. Woodruff, ICML 2019 \n"}, {"review_id": "4SeHJzy2qL", "reviewer": "Reviewer_T77r", "summary": "This paper proposes a randomized algorithm for efficiently approximating the $\\ell_p$ sensitivities, with a constant approximation parameter guaranteed.  ", "questions": "I am a bit confused about the definition of $\\boldsymbol{S}$ in line 191, where the size of $\\boldsymbol{S}$ is $r$ by $d$, should it be $r$ by $n$? Also, in lines 189 and 193 it is said $\\boldsymbol{S}$ to be 'diagonal' while this seems to contradict its definition, please let me know if I miss anything.", "limitations": "NA", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "This paper presents several novel randomized algorithms for approximating $\\ell_p$ sensitivities and related statistics, based on two key ideas:\n\n(1) Using subspace embeddings to efficiently approximate $\\|Ax\\|_p$ that avoids computing $Ax$ on the whole dataset\n\n(2) Randomly hashing the dataset into small subsets and computing the sensitivities for each subset separately\n\nGiven the above, this work provides efficient approximations of the sensitivities of all data samples, the total sensitivity, and the maximum one. The resulting algorithms are good contributions to the problem of estimating sensitivities that fill in the blank of computing $\\ell_p$ sensitivities efficiently. \n\nThe theoretical guarantees appear to be solid and correct, which is further verified by several experiments. \n\n\n", "weaknesses": " ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for their time spent reviewing our submission and are encouraged by their positive assessment of the motivation of our work, our theory, and our experiments. \n\n---------------\n\nIn response to the reviewer’s question: Yes, thank you for pointing out the typo in Line 191; the size of $S$ should be $r \\times n$. We apologize for the confusion about $S$: one can think of $S$ as being a sparse diagonal matrix, and in the resulting matrix $SA$, delete all the zero rows. We hope this clarifies the structure of $S$ and will incorporate these fixes into our manuscript. \n\n--------------- \n\nWe again thank the reviewer for their time and effort. Please let us know if there are any further questions that we could help clarify!\n"}, {"review_id": "rmDJPfzSPQ", "reviewer": "Reviewer_Vy4U", "summary": "The paper gives algorithms to compute, approximately, the individual sensitivity scores, total sensitivity and maximum sensitivity, for $\\ell_p $ norms. The approximation for the individual sensitivity scores is additive while a relative error approximation is obtained for both total and maximum sensitivity. Calculating exact sensitivity scores often being a computationally expensive problem, the authors provide faster algorithms that use techniques of multiplication with Rademacher vectors and hashing. The authors also empirically validate their results on real datasets.", "questions": "Please see weaknesses section.", "limitations": "Please see weaknesses section.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Strengths of the paper\n1) Methods of constructing coresets for various ML problems rely heavily on good approximation to sensitivity scores. As such the paper is important and will be of interest to the community.\n2) Empirical Results for various values of $p$ which are not usually seen in literature.", "weaknesses": "To me the main weakness appears to be with the writing and clarity of the paper (may be because of space constraints) and also to an extent comparison and discussion with some related works. Here I list out some of the questions/ suggestions that I have: \n1) The authors mention Lewis weights as crude over approximation to sensitivity scores for $\\ell_p$. However, another way to sample rows for $\\ell_p$ subspace embeddings is using the row norms of a well-conditioned basis. The authors have not discussed much from the area. It would be interesting to compare if for some well-conditioned basis, the row norms correspond to some approximation of sensitivity scores.\n2) In all the algorithms the authors rely on constructing SA which is an $\\ell_1$ subspace embedding. There is good amount of work on $\\ell_1$ subspace embeddings using 1-stable Cauchy random variables, exponential random variables etc. It has not been compared with. It would be useful to give a table that compares these methods in terms of time, approximation factors and no. of rows required etc. also how is SA calculated in this work and what is the time taken for it? Please clarify. \n3) Also, algorithm 2 requires computation of leverage scores of A. Are they calculated exactly in which case the time required will be $O(nd^2)$ and if they are also approximated, how does the approximation factor figure in your guarantees?\n4) The implications of the lower bound are not clarified. Please elaborate.\n5) It would be useful to give some motivation as to why maximum of the sensitivity scores is important.\n\nOverall, the paper uses some known techniques in Randomized numerical linear algebra literature to calculate the approximations to $\\ell_p$ sensitivities. However improved writing, better clarification of exact contributions and comparison with existing literature will strengthen it. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are extremely thankful to the reviewer for raising a very thorough set of questions and comments and will incorporate the clarifications for each of them in our manuscript. We are encouraged by the reviewer’s assessment of our work as important and interesting for the community. We address below the reviewer's questions and concerns. \n\n------------------------\n\n### How is $SA$ computed? ### \n We use the results of Cohen-Peng [2] to construct subspace embeddings. These incur a cost of only $O(\\text{nnz}(A))$, which is the best possible cost (since that is the cost of merely reading the input data matrix), and so we do not compare with other methods. However, we acknowledge that there exist many methods for this task (such as the ones noted by the reviewer) and will incorporate them in our \"Related Work\" section to provide a more complete picture of the landscape on subspace embeddings. \n\n----------------------\n\n### Calculation of leverage scores. ###\nAs was shown in [1], leverage scores can be computed up to a constant factor approximation in time $O(\\text{nnz}(A) + d^{\\omega})$ (we are ignoring polylogarithmic factors here). This is the accuracy to which we compute leverage scores in our submission; indeed, since the error accumulates multiplicatively across different steps of the algorithm, a high-accuracy algorithm does not serve any benefit over this constant-accuracy one. \n\n-------------------\n\n### “Implication of lower bound”. ###\nOur hardness results imply that in general one cannot compute all sensitivities as quickly as leverage scores unless there is a major breakthrough in regression. We will make this point clearer in our write-up. \n\n---------------------\n\n### “Usefulness of the maximum sensitivity”. ###\nAs detailed in our top-level response, computing the maximum sensitivity is useful in applications ranging from experiment design to differential privacy [3] to more generic pre-processing tasks such as reducing the coherence of the matrix, which was studied for $p=2$ in [4]. \n\n-------------------------------\n\n### “Row norms of a well-conditioned basis” ### \nIt is indeed true that the Lewis basis is a well-conditioned basis, but as was shown in [5], there are many advantages that sensitivity sampling offers over Lewis weight sampling, such as a much lower sample complexity for $p>2$ or when the total sensitivity is small. Additionally, while other well-conditioned bases exist and can be computed in $\\text{poly}(d)$ time, they result in worse $\\text{poly}(d)$ factor distortion, as was shown in, e.g., [4].   \n\n-------------------------\n\n### Conclusion ### \n\nWe again thank the reviewer for their time and effort in bringing up these questions; we will incorporate these answers in our manuscript. \n\nPlease let us know if any questions remain, and if we answered all the questions, we’d like to respectfully request that the reviewer re-consider their score of our submission. We are happy to answer any further questions. \n\n-------------------------\n\n### References. ### \n\n[1] “Uniform Sampling for Matrix Approximation”, Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, Aaron Sidford, ITCS 2015\n\n[2] \"$\\ell_p$ Row Sampling by Lewis Weights\", Michael B. Cohen and Richard Peng, STOC 2015 \n\n[3] “The Algorithmic Foundations of Differential Privacy”, Cynthia Dwork and Aaron Roth, Foundations and Trends in Theoretical Computer Science 2014\n\n[4] \"Dimensionality Reduction for Tukey Regression\", Kenneth L. Clarkson, Ruosong Wang, and David P. Woodruff, ICML 2019\n\n[5] \"Sharper Bounds for $\\ell_p$ Sensitivity Sampling\", David P. Woodruff and Taisuke Yasuda, ICML 2023\n"}, {"review_id": "ijtJFRHWRh", "reviewer": "Reviewer_RJDW", "summary": "The authors propose randomized algorithms to approximate $\\ell_p$ sensitivity functions for $p \\in [1,\\infty)$, which extend leverage scores beyond the $\\ell_2$ norm. The functions they consider are: 1. Estimating all sensitivities, 2. Estimating total sensitivity, and 3. Estimating maximum sensitivity. They provide different types of approximations for each task.\n\nFor task 1, they give an additive error, constant factor approximation. For task 2, they give a relative error, $(1+O(\\gamma))$ approximation ($\\gamma \\in (0,1)$). And for task 3, they give a constant factor relative error approximation.\n\nThey demonstrate their algorithm initially for $\\ell_1$ and then generalize it to all $\\ell_p$ norms with $p > 1$. They also prove a hardness result by reducing $\\ell_p$ regression to $\\ell_p$ sensitivity estimation. Additionally, they implement their algorithm to estimate all sensitivities on 2 existing datasets and compare the average and maximum approximation ratios with the theoretical results.\n\nThe main techniques used in the algorithms involve hashing using Rademacher combinations, computing sensitivities with respect to subspace embeddings, splitting matrix rows based on leverage score intervals, and utilizing existing results for $\\ell_\\infty$ subspace embeddings. ", "questions": "I would like to understand the comparison between approximate $\\ell_p$ sensitivities and Lewis weights. From what I understand, it is known that Lewis weights cover sensitivities, at least for $p \\in [1,2]$. The authors mentioned that Lewis weights are a crude approximation to sensitivities, although Lewis weights are used as a subroutine to obtain subspace embedding in their proposed algorithms. I am curious to see the benefits that an additive error approximation to sensitivities have over Lewis weight sampling. ", "limitations": "It would be beneficial if the authors consider adding a section addressing the open problems and limitations.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem of sensitivity estimation is useful for regression problems and has received limited attention for general $p$. Approximating total sensitivity is particularly significant in obtaining the sample complexity of learning arbitrary functions. \n\nThe algorithms presented are an interesting combination of known results from sensitivity sampling framework and RNLA. Specifically, their result for approximating total sensitivity is interesting because the computational complexity does not depend polynomially on the number of rows of the matrix, which is usually very large. \n\nFinally, the analysis of the proposed algorithms is largely clear, helping to understand the guarantees of their algorithms. Overall, this paper is a significant contribution to the field of sensitivity sampling and dimensionality reduction. ", "weaknesses": "\n- Algorithm 4 does not appear to be significantly novel, except for integrating generalized sensitivity to an existing $\\ell_\\infty$ subspace embedding technique, but the task is coherent with the other tasks mentioned.\n- While the motivation for total sensitivity is clear, the authors have not motivated the problem of estimating maximum sensitivity. However, this might be because I am not familiar with it. \n- I noticed some typos, missing definitions.. In Algorithm 3, it is unclear where the vector of leverage scores $\\tau(C)$ is being used. Additionally, $\\omega$ is not defined, which I am assuming to be the matrix multiplication exponent. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for their time, effort, and feedback. We are encouraged that they found the problem we study important and our algorithm and techniques novel. We are also very grateful for all the weaknesses and typos pointed out and questions raised, and we will clarify all these points in our manuscript. \n\n-----------------------\n\n### “Usefulness of the maximum sensitivity”. ### \nThe maximum sensitivity **captures the importance of the most important datapoint** and finds applications in, e.g., experiment design to detect the most important features and in reweighting matrices for low coherence [4]. Additionally, it **captures the maximum extent to which a datapoint can influence the objective function**, thus finding applications in differential privacy [1]. \n\n---------------------------\n\n### “Where is $\\tau(C)$ used in Algorithm 3?”. ### \nThe leverage scores of $M$ (which is just a submatrix of $C$) are used critically in Lines 10-14 to split the rows into buckets. This step is actually quite subtle: The sensitivity of the $j^{th}$ row of $M$ with respect to $C$, i.e., $\\sigma_1^{C}(M[j])$, is multiplicatively approximately the sensitivity of that row with respect to $SA$, which in turn is multiplicatively approximately the sensitivity of that row with respect to $A$. Further, the sensitivity of $M[j]$ with respect to $C$ can be sandwiched between appropriately scaled leverage scores of $C$. Therefore, we combine these two facts to obtain upper and lower bounds on the sensitivity of that row with respect to $A$. We formalize this notion in Lines 685-699 of the Supplementary material but will clarify the role of $\\tau(C)$ better in the main text.  \n\n--------------------------\n\n### “Using additive approximate sensitivities instead of Lewis weights”. ### \nFor the case $p > 2$, by using Theorem $1.5$ of [3], we have that $$\\text{the sample complexity with approximate sensitivities} = O\\left(\\alpha^{2p} \\cdot \\mathfrak{S}_p^{2-2/p}(A)\\right),$$ as opposed to $$\\text{the sample complexity with true sensitivities} = O\\left(\\mathfrak{S}_p^{2-2/p}(A)\\right),$$ and $$\\text{the sample complexity with Lewis weights} = O\\left(d^{p/2}\\right).$$ Assume that $p>2$ is large and the total sensitivity $\\mathfrak{S}_p(A)$ is small (say, $\\mathfrak{S}_p(A) = d$). Further assume, as a toy example, $n = d^{10}$ and $\\alpha = n^{\\frac{1}{10 p}} = d^{\\frac{1}{p}}$. Then, our approximate sensitivities give a sample complexity of $O(d^4)$, true sensitivities give a sample complexity of $O(d^2)$, and Lewis weights sampling gives $O(d^{p/2})$. Thus, our approximate sensitivities preserve the regression approximation guarantee, while **increasing the total sample complexity by only a small $\\text{poly}(d)$ compared to that with true sensitivities (while still being much smaller than that with Lewis weights) and incurring a much smaller cost.** A detailed discussion of the connection of Lewis weights to well-conditioned bases can be found in Section $3.3$ of [4]. \n\n-----------------\n\n### Conclusion ###\n\nWe once again thank the reviewer for their time, effort, and very thoughtful questions. We are happy to provide any further clarification required.\n\n------------------\n\n### References ### \n\n[1] “The Algorithmic Foundations of Differential Privacy”, Cynthia Dwork and Aaron Roth, Foundations and Trends in Theoretical Computer Science 2014\n\n[2] “Uniform Sampling for Matrix Approximation”, Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, Aaron Sidford, ITCS 2015\n\n[3] “Sharper Bounds for $\\ell_p$ Sensitivity Sampling”, David P. Woodruff and Taisuke Yasuda, ICML 2023\n\n[4] \"Dimensionality Reduction for Tukey Regression\", Kenneth L. Clarkson, Ruosong Wang, and David P. Woodruff, ICML 2019, https://arxiv.org/abs/1905.05376\n"}]

Author rebuttal:
# Top-Level Response #

We thank all the reviewers for their time, effort, and suggestions. Here, we restate some of our key contributions and answer common questions. We also reply to each reviewer individually. 

---

## Motivating the problem ##
### Why sensitivities? ###
A common preprocessing step for $\ell_p$ regression ($\min_{x} \|Ax\|_p^p$) involves constructing an $\ell_p$ subspace embedding of matrix $A$. This is usually done via a sampling matrix $S$ that ensures $\|S A x\|_p^p = (1\pm \epsilon) \|Ax\|_p^p$ for all vectors $x$. This reduces our problem's data dimension from the number of rows of $A$ to that of $SA$. 

A recent result [1] shows that **a sampling matrix $S$ built using the $\ell_p$ sensitivities of $A$ is much smaller than one built using its $\ell_p$ Lewis weights in many important cases**, e.g., when $p > 2$, or when the total sensitivity is small. Consequently, sensitivity sampling offers significant advantages over Lewis weight sampling, which also extends to many ML applications beyond regression [1, Section 1.1]. 

However, until our work, **there were no known efficient algorithms to approximate $\ell_p$ sensitivities for general $p>0$.** Our paper initiates a systematic study of algorithms for this task.

### Why total sensitivity? (asked by R1). ###
As alluded to above, when the total sensitivity is small, the sample complexity of sensitivity sampling is much lower than that of Lewis weight sampling [1]. Hence, a quick approximation to the total sensitivity can be used as **a fast test for whether or not to proceed with sensitivity sampling (involving the costly task of calculating all sensitivities).**  

Further, for $\ell_2$ sensitivities, the total sensitivity is the rank, which may be used in rank estimation subroutines [6]. The total $\ell_p$ sensitivity is analogous to $\ell_p$ rank, which may be similarly useful.

### Why maximum sensitivity? (asked by R1, R3, and R4). ###
 The maximum sensitivity **captures the importance of the most important datapoint** and is used in, e.g., experiment design and in reweighting matrices for low coherence [7]. Additionally, it **captures the maximum extent to which a datapoint can influence the objective function** and is used in differential privacy [5]. 

---

## Clarification on our results. ##
We first clarify that some statements of our results in the main submission were correct but sub-optimally stated — we have corrected these in our full version (cf. Lines 510 - 519, Lines 72 - 73, Lines 83 - 84, and Line 89 of the Supplementary file). **These sub-optimal statements existed only in our theorems, and the algorithms and analyses originally presented indeed give the correct cost (e.g., cf. Lines 551 - 575 of the Supplementary file).** 

### Computing all sensitivities. ###
Our correct main result for $p=1$ (cf. Lines 73-74 of the Supplementary Material) is that we can approximate all $\ell_1$ sensitivities (with additive error $\frac{\alpha}{n}\mathfrak{S}_1(A)$) at a cost of $ O(nnz(A) + (n/\alpha) \cdot d^\omega). $ As **R1** noted, the naive total cost is $O(n (nnz(A) + d^\omega))$. Thus, **our obtained result greatly improves over the naive and (to our best knowledge) only known result**. 

### Proof sketch. ###
 As a reminder, the proof sketch for this cost is: Constructing an $\ell_1$ subspace embedding $SA$ costs $O(nnz(A))$ [3]  (asked by **R3**); plus, computing $n/\alpha$ sensitivities with respect to the $d \times d$ matrix $SA$ costs $(n/\alpha) \cdot d^\omega$ (by reducing one sensitivity computation w.r.t. $SA$ to solving a $d\times d$ linear program and using the current fastest LP solver [4]).
 
### Sampling with approximate sensitivities (asked by R1 and R4). ###
For $p > 2$, Theorem $1.5$ of [1] implies $$\text{the sample complexity with approximate sensitivities} = O\left(\alpha^{2p} \cdot \mathfrak{S}_p^{2-2/p}(A)\right),$$ as opposed to $$\text{the sample complexity with true sensitivities} = O\left(\mathfrak{S}_p^{2-2/p}(A)\right),$$ and $$\text{the sample complexity with Lewis weights} = O\left(d^{p/2}\right).$$ Assume a large $p>2$ and small total sensitivity, say, $\mathfrak{S}_p(A) = d$. Further assume, as a toy example, $n = d^{10}$ and $\alpha = n^{\frac{1}{10 p}} = d^{\frac{1}{p}}$. Then, sample complexity with approximate sensitivities is $O(d^4)$, with true sensitivities is $O(d^2)$, and with Lewis weights is $O(d^{p/2})$. Thus, our approximate sensitivities preserve the regression approximation guarantee, while **increasing the total sample complexity by only a small $\text{poly}(d)$ compared to that with true sensitivities (while still being much smaller than that with Lewis weights) and incurring a much smaller cost.**

### Computing the total and maximum sensitivities. ###
For the total and maximum $\ell_1$ sensitivities, our runtimes are, respectively, $ O(nnz(A) + d^{\omega + 1/2})$ and $O(nnz(A) + d^{\omega+1}),$ with **no polynomial dependence on $n$** and significantly better than the naive cost.

---
## References ##
[1] “Sharper Bounds for $\ell_p$ Sensitivity Sampling”, David P. Woodruff and Taisuke Yasuda, ICML 2023

[2] “Uniform Sampling for Matrix Approximation”, Michael B. Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, Aaron Sidford, ITCS 2015

[3] "$\ell_p$ Row Sampling by Lewis Weights", Michael B. Cohen and Richard Peng, STOC 2015 

[4] "Solving Linear Programs in the Current Matrix Multiplication Time", Michael B. Cohen, Yin Tat Lee, and Zhao Song, STOC 2019

[5] “The Algorithmic Foundations of Differential Privacy”, Cynthia Dwork and Aaron Roth, Foundations and Trends in Theoretical Computer Science 2014

[6] “Rank Estimation For (Approximately) Low-Rank Matrices”, Niloofar Bayat, Cody Morrin, Yuheng Wang, and Vishal Misra, ACM SIGMETRICS Performance Evaluation Review 2022

[7] "Dimensionality Reduction for Tukey Regression", Kenneth L. Clarkson, Ruosong Wang, and David P. Woodruff, ICML 2019


Meta-review: This paper attacks the problem of computing l_p sensitivities using tools from approximation algorithms and randomized algorithms. This is related to other linear algebra problems which have attracted a lot of interest like estimating leverage scores and lewis weights. After the author-reviewer discussion, some important points of confusion/typos regarding the main results in the original draft were clarified, and a positive consensus was formed of the overall technical contribution. Based on this, I recommend acceptance of this work. Please make the promised edits in the next version of the paper.
Decision: Accept (poster)