# Computing Approximate $\ell_{p}$ Sensitivities††thanks: Authors listed
alphabetically. A preliminary version of this paper was first published in
NeurIPS 2023.

Swati Padmanabhan
Massachusetts Institute of Technology. Email<EMAIL_ADDRESS>Most of the work
done while a Ph.D. student at UW Seattle and as a student researcher with
Google. David P. Woodruff
Carnegie Mellon University. Email<EMAIL_ADDRESS>Qiuyi (Richard) Zhang
Google Research. Email<EMAIL_ADDRESS>
###### 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 $\ell_{p}$ sensitivities, which we
show is equivalent to approximate $\ell_{p}$ regression, are known for only
the $\ell_{2}$ setting, in which they are termed leverage scores.

In this work, we provide efficient algorithms for approximating $\ell_{p}$
sensitivities and related 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 $O(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 to the total sensitivity at the cost of roughly $O(\sqrt{d})$
sensitivity computations. Furthermore, we estimate the maximum $\ell_{1}$
sensitivity, up to a $\sqrt{d}$ factor, using $O(d)$ sensitivity computations.
We generalize all these results to $\ell_{p}$ norms for $p>1$. Lastly, we
experimentally show that for a 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.

##  1 Introduction

Many modern large-scale machine learning datasets comprise tall matrices
$\mathbf{A}\in\mathbb{R}^{n\times d}$ with $d$ features and $n\gg d$ training
examples, often making it computationally expensive to run them through even
basic algorithmic primitives. Therefore, many applications spanning dimension
reduction, privacy, and fairness aim to estimate the importance of a given
datapoint as a first step. Such importance estimates afford us a principled
approach to focus our attention on a 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, 2, 3, 4]
such as [5, 6, 7, 8, 9, 10, 11, 12, 13, 14] developed for the “finite-sum
minimization” problem

$\begin{array}[]{ll}\mbox{minimize}_{\boldsymbol{\theta}\in\mathcal{X}}&\sum_{i=1}^{n}f_{i}(\boldsymbol{\theta}).\end{array}$

where, in the context of empirical risk minimization, each
$f_{i}:\mathcal{X}\mapsto\mathbb{R}_{\geq 0}$ in Section 1 measures the loss
incurred by the $i$-th data point from the training set, and in generalized
linear models, each $f_{i}$ represents a link function applied to a linear
predictor evaluated at the $i$-th data point. The aforementioned algorithms
randomly sample a function $f_{i}$ and make progress via gradient estimation
techniques.

However, uniform sampling can lead to significant information loss when the
number of important training examples is 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, 18, 19, 20,
21]. One such technique is based on _sensitivity scores_ , defined [15] for
each $i\in[n]$ as:

$\boldsymbol{\sigma}_{i}\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\max_{\mathbf{x}\in\mathcal{X}}\frac{f_{i}(\mathbf{x})}{\sum_{j=1}^{n}f_{j}(\mathbf{x})}.$

Sampling each $f_{i}$ independently with probability $p_{i}\propto\sigma_{i}$
and scaling the sampled row with $1/p_{i}$ preserves the objective in Section
1 in expectation for every $\mathbf{x}\in\mathcal{X}$. Further, one gets a
$(1\pm\varepsilon)$-approximation to this objective by sampling
$\widetilde{O}(\mathfrak{S}d\varepsilon^{-2})$ functions [22, Theorem $1.1$],
where $\mathfrak{S}=\sum_{i=1}^{n}\boldsymbol{\sigma}_{i}$ is called the total
sensitivity, and $d$ is an associated VC dimension.

Numerous advantages of sensitivity sampling over $\ell_{p}$ Lewis weight
sampling [23] were highlighted in the recent work [24], prominent among which
is that a sampling matrix built using the $\ell_{p}$ sensitivities of
$\mathbf{A}$ is significantly smaller than one built using its $\ell_{p}$
Lewis weights in many important cases, e.g., when the total sensitivity
$\mathfrak{S}$ is small for $p>2$ (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
$\epsilon$-approximators for shape-fitting functions, also called strong
coresets [15, 16, 26, 27, 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]. This is in contrast to the progress on fast computation of _leverage
scores_ and _$\ell_{p}$ Lewis weights_, two closely related quantities
extensively used in importance sampling (cf. Section 2). As a brief example,
for $f_{i}(\mathbf{x})=(\mathbf{a}_{i}^{\top}\mathbf{x})^{2}$, when the
sensitivity score of $f_{i}$ is exactly the leverage score of
$\mathbf{a}_{i}$, the fastest algorithm [36] can compute all sensitivities in
runtime on the order of $O(\log(n))$ independent sensitivity calculations, as
opposed to the $O(n)$ sensitivity calculations one would naively require.
Similarly, there exist fast algorithms [23, 37] for computing even high-
accuracy $\ell_{p}$ Lewis weights for all $p>0$, at the cost of computing
$\textrm{poly}(p)$ leverage scores.

###  1.1 Our contributions

In this work, we initiate a systematic study of algorithms for approximately
computing the $\ell_{p}$ sensitivities of a matrix, generally up to a constant
factor. We first define $\ell_{p}$ sensitivities:

######  Definition 1.1 ($\ell_{p}$ sensitivities [15]).

Given a full-rank matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with $n\geq d$
and a scalar $p\in(0,\infty)$, let $\mathbf{a}_{i}$ denote the $i$’th row of
matrix $\mathbf{A}$. Then the vector of $\ell_{p}$ sensitivities of
$\mathbf{A}$ is defined as111From hereon, for notational conciseness, we omit
stating $\mathbf{A}\mathbf{x}\neq\mathbf{0}$ in the problem constraint. (cf.
Section 2 for the notation)

$\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\max_{\mathbf{x}\in\mathbb{R}^{d},\mathbf{A}\mathbf{x}\neq\mathbf{0}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|^{p}}{\|\mathbf{A}\mathbf{x}\|_{p}^{p}},\text{
for all $i\in[n]$},$

and the total $\ell_{p}$ sensitivity is defined as
$\mathfrak{S}_{p}(\mathbf{A}){\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}}\sum_{i\in[n]}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})$.

These $\ell_{p}$ sensitivities yield sampling probabilities for $\ell_{p}$
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 [38], sparse recovery [39],
graph-based semi-supervised learning [40, 41, 42], and data clustering [43].
However, algorithms for approximating $\ell_{p}$ sensitivities were known for
only $p=2$ (when they are called “leverage scores”). As is typical in large-
scale machine learning, we assume $n\gg d$. We now state our contributions.

(1)

Approximation algorithms. We design algorithms for three prototypical
computational tasks associated with $\ell_{p}$ sensitivities of a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$. Specifically, for an approximation
parameter $\alpha>1$, we can simultaneously estimate all $\ell_{p}$
sensitivities up to an additive error of
$O\left(\frac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A})\right)$ via
$O(n/\alpha)$ individual sensitivity calculations, reducing our runtime by
sacrificing some accuracy. We limit most of our results in the main text to
$\ell_{1}$ sensitivities, but extend these to all $p\geq 1$ in Appendix C. To
state our results for the $\ell_{p}$ setting, we use
$\mathop{\mathbf{nnz}}(\mathbf{A})$ to denote the number of non-zero entries
of the matrix $\mathbf{A}$, $\omega$ for the matrix multiplication constant,
and we introduce the following piece of notation.

######  Definition 1.2 (Notation for $\ell_{p}$ results).

We introduce the notation $\mathsf{LP}(m,d,p)$ to denote the cost of
approximating one $\ell_{p}$ sensitivity of an $m\times d$ matrix up to an
accuracy of a given constant factor.

######  Theorem 1.3 (Estimating all sensitivities; informal version of Theorem
3.1 and Theorem C.1).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and $1\leq\alpha\ll n$,
there exists an algorithm (Algorithm 1) which returns a vector
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{\geq 0}$ such that, with
high probability, for each $i\in[n]$, we have

$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})+\tfrac{\alpha}{n}\mathfrak{S}_{1}(\mathbf{A})).$

Our algorithm’s runtime is
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\frac{n}{\alpha}\cdot
d^{\omega}\right)$. More generally, for $p\geq 1$, we can compute (via
Algorithm 6) an estimate
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{\geq 0}$ satisfying
$\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\alpha^{p-1}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})+\frac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A}))$
with high probability for each $i\in[n]$, at a cost of
$O\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\tfrac{n}{\alpha}\cdot\mathsf{LP}(O(d^{\max(1,p/2)}),d,p)\right)$.

Two closely related properties arising widely in algorithms research are the
_total sensitivity_ $\mathfrak{S}_{1}(\mathbf{A})$ and the _maximum
sensitivity_ $\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}$. 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 $\ell_{1}$ 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.

While the total sensitivity helps us 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 [44]. Additionally, it captures the
maximum extent to which a datapoint can influence the objective function and
is therefore used in differential privacy [45]. While one could naively
estimate the total and maximum sensitivities by computing all $n$
sensitivities using Theorem 1.3, we give faster algorithms that have no
polynomial dependence on $n$, which, in the assumed regime of $n\gg d$, is an
improvement on the runtime.

######  Theorem 1.4 (Estimating the total sensitivity; informal version of
Theorem 3.5).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, a scalar $p\geq 1$, and
a small constant $\gamma<1$, there exists an algorithm (Algorithm 2), which
returns a positive scalar $\widehat{s}$ such that, with a probability of at
least $0.99$, we have

$\mathfrak{S}_{p}(\mathbf{A})\leq\widehat{s}\leq(1+O(\gamma))\mathfrak{S}_{p}(\mathbf{A}).$

Our algorithm’s runtime is
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\frac{1}{\gamma^{2}}\cdot
d^{|p/2-1|}\mathsf{LP}(O(d^{\max(1,p/2)},d,p)\right)$.

In Algorithm 4 (presented in Section B.1), we give an alternate method to
estimate $\mathfrak{S}_{1}(\mathbf{A})$ based on recursive leverage score
sampling, without $\ell_{p}$ Lewis weight calculations, which we believe to be
of independent interest. Our Algorithm 2 utilizes approximation properties of
$\ell_{p}$ 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.2).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, there exists an
algorithm (Algorithm 3) which returns a scalar
$\widetilde{\boldsymbol{\sigma}}>0$ such that

$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}\leq\widetilde{\boldsymbol{\sigma}}\leq
O(\sqrt{d}\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}).$

Our algorithm’s runtime is
$\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+d^{\omega+1})$. Furthermore,
via Algorithm 7, this result holds for all $p\in[1,2]$ and for $p>2$, this
guarantee costs
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+d^{p/2}\cdot\mathsf{LP}(O(d^{p/2}),d,p)\right)$.

(2)

Hardness results. In the other direction, while it is known that $\ell_{p}$
sensitivity calculation reduces to $\ell_{p}$ regression, we show the
converse. Specifically, we establish hardness of computing $\ell_{p}$
sensitivities in terms of the cost of solving multiple corresponding
$\ell_{p}$ regression problems, by showing that a
$(1+\varepsilon)$-approximation to $\ell_{p}$-sensitivities solves $\ell_{p}$
multi-regression up to an accuracy factor of $1+\varepsilon$.

######  Theorem 1.6 (Leave-one-out $\ell_{p}$ multiregression reduces to
$\ell_{p}$ sensitivities; informal222Please see Lemma D.2 for the formal
version with a small additive error.).

Suppose that we are given an sensitivity approximation algorithm
$\mathcal{A}$, which for any matrix
$\mathbf{A}^{\prime}\in\mathbb{R}^{n^{\prime}\times d^{\prime}}$ and accuracy
parameter $\varepsilon^{\prime}\in(0,1)$, computes
$(1\pm\varepsilon^{\prime})\boldsymbol{\sigma}_{p}(\mathbf{A}^{\prime})$ in
time
$\mathcal{T}(n^{\prime},d^{\prime},\mathop{\mathbf{nnz}}(\mathbf{A}^{\prime}),\varepsilon^{\prime})$.
Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ with $n\geq d$, let
$\textrm{OPT}_{i}:=\min_{\mathbf{y}\in\mathbb{R}^{d-1}}\|\mathbf{A}_{:-i}\mathbf{y}+\mathbf{A}_{:i}\|_{p}^{p}$
and
$\mathbf{y}_{i}^{\star}:=\arg\min_{\mathbf{y}\in\mathbb{R}^{d-1}}\|\mathbf{A}_{:-i}\mathbf{y}+\mathbf{A}_{:i}\|_{p}$
for all the $i\in[d]$. Then, there exists an algorithm that takes
$\mathbf{A}\in\mathbb{R}^{n\times d}$ and computes
$(1\pm\varepsilon)\textrm{OPT}_{i}$ for all $i$ in time
$\mathcal{T}(n+d,d,\mathop{\mathbf{nnz}}(\mathbf{A}),\varepsilon)$.

For $p\neq 2$, approximating all $\ell_{p}$ sensitivities to a high accuracy
would therefore require improving $\ell_{p}$ multi-regression algorithms with
$\Omega(d)$ instances, for which the current best high-accuracy algorithms
take $\text{poly}(d)\mathop{\mathbf{nnz}}(\mathbf{A})$ time [46, 47]. More
concretely, our hardness results imply that in general one cannot compute all
sensitivities as quickly as leverage scores unless there is a major
breakthrough in multi-response regression. In order to show this, we introduce
a reduction algorithm that efficiently regresses columns of $\mathbf{A}$
against linear combinations of other columns of $\mathbf{A}$, by simply adding
a row to $\mathbf{A}$ capturing the desired linear combination and then
computing sensitivities, which will reveal the cost of the corresponding
regression problem. We can solve $\Theta(d)$ of these problems by augmenting
the matrix to increase the rows and columns by at most a 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 [48] 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$-$5$x better. We also found these sensitivities to be easy to
approximate, with the accuracy-runtime tradeoff much better than our theory
suggests, with up to $40$x 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 $4$x.

###  1.2 Related work

Introduced in the pioneering work of [15], sensitivity sampling falls under
the broad framework of “importance sampling”. It has found use in constructing
$\epsilon$-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, 27, 28], clustering [16, 29, 30], logistic regression
[31, 32], least squares regression [33, 34, 49], principal component analysis
[50], reinforcement learning [51, 52], and pruning of deep neural networks
[53, 54, 55]. As mentioned earlier, [24] made significant progress in
sharpening bounds on sample complexities obtained with sensitivity sampling.
Sampling algorithms for $\ell_{p}$ subspace embeddings have also been
extensively studied in functional analysis [56, 57, 58, 59, 60, 61] as well as
in theoretical computer science [23, 36, 46, 47].

Another notable line of work in importance sampling uses Lewis weights to
determine sampling probabilities. As we explain in LABEL:{sec:notation}, 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 $\ell_{p}$ regression. They have
subsequently been used in row sampling algorithms for data pre-processing [33,
62, 63, 36, 23], computing dimension-free strong coresets for $k$-median and
subspace approximation [64], and fast tensor factorization in the streaming
model [65]. Lewis weights are also used for $\ell_{1}$ regression in: [66] for
stochastic gradient descent pre-conditioning, [67] for quantile regression,
and [68] to provide algorithms for linear algebraic problems in the sliding
window model. Lewis weights have also become widely used in convex geometry
[69], randomized numerical linear algebra [23, 44, 70, 65, 71], and machine
learning [72, 73, 74].

##  2 Notation and preliminaries

We use boldface uppercase and lowercase letters for matrices and vectors
respectively. We denote the $i^{\mathrm{th}}$ row vector of a matrix
$\mathbf{A}$ by $\mathbf{a}_{i}$. Given a matrix $\mathbf{M}$, we use
$\mathrm{Tr}(\mathbf{M})$ for its trace, $\mathrm{rank}(\mathbf{M})$ for its
rank, $\mathbf{M}^{\dagger}$ for its Moore-Penrose pseudoinverse, and
$|\mathbf{M}|$ for the number of its rows. When $x$ is an integer, we use
$[x]$ for the set of integers $1,2,\dotsc,x$. For two positive scalars $x$ and
$y$ and a scalar $\alpha\in(0,1)$, we use $x\approx_{\alpha}y$ to denote
$(1-\alpha)y\leq x\leq(1+\alpha)y$; we sometimes use $x\approx y$ to mean
$(1-c)y\leq x\leq(1+c)y$ for some appropriate universal constant $c$. We use
$\widetilde{O}$ to hide dependence on $n^{o(1)}$. We acknowledge that there is
some notation overloading between $p$ (when used to denote the scalar in, for
example, $\ell_{p}$ norms) and $p_{i}$ (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:
$\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i})$ to denote the
$\ell_{p}$ sensitivity of $\mathbf{a}_{i}$ with respect to $\mathbf{B}$; when
computing the sensitivity with respect to the same matrix that the row is
drawn from, we omit the matrix superscript:

$\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i}):=\max_{\mathbf{x}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|^{p}}{\|\mathbf{B}\mathbf{x}\|_{p}^{p}}\text{
and
}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}):=\max_{\mathbf{x}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|^{p}}{\|\mathbf{A}\mathbf{x}\|_{p}^{p}}.$

When referring to the vector of $\ell_{p}$ sensitivities of all the rows of a
matrix, we omit the subscript in the argument: so,
$\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{A})$ is the vector of $\ell_{p}$
sensitivities of rows of the matrix $\mathbf{A}$, each computed with respect
to the matrix $\mathbf{B}$ (as in Section 2), and analogously,
$\boldsymbol{\sigma}_{p}(\mathbf{A})$ is the vector of $\ell_{p}$
sensitivities of the rows of matrix $\mathbf{A}$, each computed with respect
to itself. Similarly, the total sensitivity is the sum of $\ell_{p}$
sensitivities of the rows of $\mathbf{A}$ with respect to matrix $\mathbf{B}$,
for which we use the following notation:

$\mathfrak{S}_{p}^{\mathbf{B}}(\mathbf{A}):=\sum_{i\in[|\mathbf{A}|]}\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i})\text{
and
}\mathfrak{S}_{p}(\mathbf{A}):=\sum_{i\in[|\mathbf{A}|]}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}),$

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 $\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i})$ could be
infinite, by appending $\mathbf{a}_{i}$ to $\mathbf{B}$, the generalized
sensitivity becomes once again contained in $[0,1]$. Specifically, a simple
rearrangement of the definition gives us the identity
$\boldsymbol{\sigma}_{p}^{\mathbf{B}\cup\mathbf{a}_{i}}(\mathbf{a}_{i})=1/(1+1/\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i}))$.
We now define leverage scores, which are a special type of sensitivities and
also amenable to efficient computation [36], which makes them an ideal
“building block” in algorithms for approximating sensitivities.

##### Leverage scores are $\ell_{2}$ sensitivities.

Leverage scores are a key tool in theoretical computer science, machine
learning, and statistics that quantify the importance of any row of a matrix
in composing its rowspace. Formally, for a matrix
$\mathbf{A}\in\mathbb{R}^{n\times d}$, one may define the $i^{\mathrm{th}}$
leverage score as
$\tau_{i}(\mathbf{A})=\min_{\mathbf{x}\in\mathbb{R}^{n}:\mathbf{A}^{\top}\mathbf{x}=\mathbf{a}_{i}}\|\mathbf{x}\|_{2}^{2}$
[36]. This characterization captures the notion that if a row $\mathbf{a}_{i}$
is orthogonal to all other rows of $\mathbf{A}$, then its leverage score is
$1$, and if all the rows are the same, then their leverage scores are all
$d/n$. The previous convex program may be solved to see that the leverage
score of the $i^{\mathrm{th}}$ row $\mathbf{a}_{i}$ of $\mathbf{A}$ has the
closed-form solution
$\tau_{i}(\mathbf{A})\stackrel{{\scriptstyle\mathrm{{\scriptscriptstyle
def}}}}{{=}}\mathbf{a}_{i}^{\top}(\mathbf{A}^{\top}\mathbf{A})^{\dagger}\mathbf{a}_{i}$.
The $i^{\mathrm{th}}$ leverage score is also the $i^{\mathrm{th}}$ diagonal
entry of the orthogonal projection matrix
$\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{\dagger}\mathbf{A}^{\top}$. Since
the eigenvalues of an orthogonal projection matrix are zero or one, it implies
$\tau_{i}(\mathbf{A})=\mathbf{1}^{\top}_{i}\mathbf{A}(\mathbf{A}^{\top}\mathbf{A})^{\dagger}\mathbf{A}^{\top}\mathbf{1}_{i}\leq
1$ for all $i\in[n]$, and the sum of all leverage scores satisfies
$\sum_{i=1}^{n}\tau_{i}(\mathbf{A})\leq d$ (see [75]).

It turns out that
$\tau_{i}(\mathbf{A})=\boldsymbol{\sigma}_{2}(\mathbf{a}_{i})$ (see [76]).
This may be verified by applying a change of basis to the variable
$\mathbf{y}=\mathbf{A}\mathbf{x}$ in the definition of sensitivity and working
out that the maximizer is the vector parallel to
$\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{a}_{i}$. This also gives an explicit
formula for the total $\ell_{2}$ sensitivity:
$\mathfrak{S}_{2}(\mathbf{A})=\textrm{rank}(\mathbf{A})$.

The fastest algorithm for constant factor approximation to leverage scores is
by [36], as stated next.

######  Fact 2.1 ([36, Lemma 7]).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, we can compute constant-
factor approximations to all its leverage scores in time
$O(\mathop{\mathbf{nnz}}(\mathbf{A})\log(n)+d^{\omega}\log(d))$.

In this paper, the above is the result we use to compute leverage scores (to a
constant accuracy). 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. We now state the connections between leverage scores and
$\ell_{p}$ sensitivities for $p\neq 2$.

Leverage scores approximate $\ell_{1}$ sensitivities up to a distortion factor
of $\sqrt{n}$, as seen from the below known fact we prove, for completeness,
in Appendix A.

######  Fact 2.2 (Crude sensitivity approximation via leverage scores).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, let
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})$ denote the $i^{\mathrm{th}}$
$\ell_{1}$ sensitivity of $\mathbf{A}$ with respect to $\mathbf{A}$, and let
$\tau_{i}(\mathbf{A})$ denote its $i^{\mathrm{th}}$ leverage score with
respect to $\mathbf{A}$. Then we have
$\sqrt{\frac{\tau_{i}(\mathbf{A})}{n}}\leq\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\leq\sqrt{\tau_{i}(\mathbf{A})}.$

In general, $\ell_{p}$ sensitivities of rows of a matrix suffer a distortion
factor of at most $n^{|1/2-1/p|}$ from corresponding leverage scores, the
proof of which follows from Hölder’s inequality in the same as Fact 2.2. The
$\ell_{p}$ generalizations of leverages scores are called $\ell_{p}$ Lewis
weights [56, 23], which may be defined as the unique333The existence and
uniqueness of such a $\mathbf{w}$ was shown by D.R.Lewis, after whom the
weights are named [56] vector of weights satisfying the recurrence
$\mathbf{w}_{i}=\tau_{i}(\mathbf{W}^{1/2-1/p}\mathbf{A})$. Lewis weights have
notably been at the core of recent breakthroughs in fast algorithms for linear
programs, and we refer the reader to [77] for these details and alternate
definitions and fundamental properties of Lewis weights.

Unlike when $p=2$, $\ell_{p}$ Lewis weights do not equal $\ell_{p}$
sensitivities, and from an approximation perspective, they provide only a one-
sided bound: $\boldsymbol{\sigma}_{p}^{\mathbf{A}}(\mathbf{a}_{i})\leq
d^{\,\max(0,p/2-1)}\,\mathbf{w}^{\mathbf{A}}_{p}(\mathbf{a}_{i})$ [76, 78].
While we have algorithms for efficiently computing leverage scores [36] and
Lewis weights [23, 37], extensions for utilizing these ideas for fast and
accurate sensitivity approximation algorithms for all $p$ 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 [79].

##### Subspace embeddings.

We use sampling-based constant-approximation $\ell_{p}$ subspace embeddings
that we denote by $\mathbf{S}_{p}$. When the type of embedding is clear from
context, we omit the subscript.

######  Definition 2.3 ($\ell_{p}$ subspace embedding).

Given $\mathbf{A}\in\mathbb{R}^{n\times d}$ (typically $n\gg d$), we call
$\mathbf{S}\in\mathbb{R}^{r\times n}$ (typically $d\leq r\ll n$) an
$\epsilon$-approximate $\ell_{p}$ subspace embedding of $\mathbf{A}$ if
$\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{p}{\approx_{\varepsilon}}\|\mathbf{A}\mathbf{x}\|_{p}$
$\forall\mathbf{x}\in\mathbb{R}^{d}$.

When the diagonal matrix $\mathbf{S}$ is defined as
$\mathbf{S}_{ii}=1/\sqrt{p_{i}}$ with probability
$p_{i}=\Omega(\tau_{i}(\mathbf{A}))$ (and $0$ otherwise), then with high
probability $\mathbf{S}\mathbf{A}$ is an $\ell_{2}$ subspace embedding for
$\mathbf{A}$ [33, 80]; further, $\mathbf{S}$ has at most
$\widetilde{O}(d\varepsilon^{-2})$ non-zero entries. For general $p$, we use
Lewis weights to construct $\ell_{p}$ subspace embeddings fast [23, Theorem
7.1]. This fast construction has been made possible by the reduction of Lewis
weight computation to a few (polynomial in $p$) leverage score computations
for $p<4$ [23] as well as $p\geq 4$ [37]. This efficient construction is the
reason for our choice of this subspace embedding. We note that there has been
an extensive amount of literature on designing $\ell_{1}$ subspace embeddings
for example using Cauchy [81] and exponential [82] random variables.

For $p\leq 2$, it is known that the expected sample complexity of rows is
$\sum_{i=1}^{n}p_{i}=O(d)$ [34], and for $p>2$, this is known to be at most
$O(d^{p/2})$ [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 case applications,
particularly those seen in practice. This gap is captured by the total
sensitivity $\mathfrak{S}_{p}(\mathbf{A})$, which has been used to perform
more sample efficient subspace embeddings [68, 27, 28], as sampling
$\widetilde{O}(\mathfrak{S}_{p}(\mathbf{A})d/\epsilon^{2})$ rows suffice to
provide an $\ell_{p}$ subspace embedding.

##### Accurate sensitivity calculations.

While leverage scores give a crude approximation, constant-factor
approximations of $\ell_{p}$ 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 $\ell_{p}$ regression since by
scale invariance, we have
$\frac{1}{\boldsymbol{\sigma}_{p}^{\mathbf{B}}(\mathbf{a}_{i})}=\min_{\mathbf{a}_{i}^{\top}\mathbf{x}=1}{\|\mathbf{B}\mathbf{x}\|_{p}^{p}}$.

This reduction allows us to use recent developments in fast approximate
algorithms for $\ell_{p}$ regression [46, 25, 47], which utilize subspace
embeddings to first reduce the number of rows of the matrix to
$O(\text{poly}(d))$. For the specific case of $p=1$, which is the focus of the
main body of our paper, we may reduce to approximate $\ell_{1}$ regression on
a $O(d)\times d$ matrix, which is a form of linear programming.

######  Fact 2.4.

Given an $n\times d$ matrix, the cost of computing $k$ of its $\ell_{1}$
sensitivities is $\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+k\cdot
d^{\omega})$.

##  3 Approximating functions of sensitivities

We first present a constant-probability algorithm approximating the $\ell_{1}$
sensitivities in Algorithm 1. Our algorithm is a natural one: Since computing
the $\ell_{1}$ sensitivities of all the rows simultaneously is computationally
expensive, we instead hash $\alpha$-sized blocks of random rows into smaller
(constant-sized) buckets and compute the sensitivities of the smaller,
$O(n/\alpha)$-sized matrix $\mathbf{P}$ so generated, computing the
sensitivities with respect to $\mathbf{S}\mathbf{A}$, the $\ell_{1}$ subspace
embedding of $\mathbf{A}$. Running this algorithm multiple times gives our
high-probability guarantee via the standard median trick.

Algorithm 1 Approximating $\ell_{1}$-Sensitivities: Row-wise Approximation

1:Input: Matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, approximation factor
$\alpha\in(1,n)$

2:Output: Vector $\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{>0}$ that
satisfies, for each $i\in[n]$, with probability $0.9$, that
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i}){}\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\boldsymbol{\sigma}_{1}(\mathbf{a}_{i}){}+\mathfrak{S}_{1}(\mathbf{A})\tfrac{\alpha}{n})$

3:Compute for $\mathbf{A}$ an $\ell_{1}$ subspace embedding
$\mathbf{S}\mathbf{A}\in O(d)\times d$ (cf. Definition 2.3)

4:Partition $\mathbf{A}$ into $\tfrac{n}{\alpha}$ blocks
$\mathbf{B}_{1},\dotsc,\mathbf{B}_{n/\alpha}$ each comprising $\alpha$
randomly selected rows

5:for the $\ell^{\mathrm{th}}$ block $\mathbf{B}_{\ell}$, with
$\ell\in[\tfrac{n}{\alpha}]$ do

6: Sample $\alpha$-dimensional independent Rademacher vectors
$\mathbf{r}_{1}^{(\ell)},\dotsc,\mathbf{r}_{100}^{(\ell)}$

7: For each $j{\in}[100]$, compute the row vectors444We are omitting the
transposes on $\mathbf{r}_{i}^{\ell}$ to avoid a cluttered notation.
$\mathbf{r}_{j}^{(\ell)}\mathbf{B}_{\ell}\in\mathbb{R}^{d}$

8:end for

9:Let $\mathbf{P}\in\mathbb{R}^{100\tfrac{n}{\alpha}\times d}$ be the matrix
of all vectors from 7. Compute
$\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{P})$ using Fact 2.4

10:for each $i\in[n]$ do

11: Denote by $J$ the set of row indices in $\mathbf{P}$ that $\mathbf{a}_{i}$
is mapped to in 7

12: Set $\widetilde{\boldsymbol{\sigma}}_{i}=\max_{j\in
J}(\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j}))$

13:end for

14:Return $\widetilde{\boldsymbol{\sigma}}$

######  Theorem 3.1.

Given a full-rank matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and an
approximation factor $1<\alpha\ll n$, let
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})$ be the $\ell_{1}$ sensitivity of the
$i^{\mathrm{th}}$ row of $\mathbf{A}$. Then there exists an algorithm that, in
time
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\frac{n}{\alpha}\cdot
d^{\omega}\right)$, returns
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{\geq 0}$ such that with
high probability, for each $i\in[n]$,

$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})+\tfrac{\alpha}{n}\mathfrak{S}_{1}(\mathbf{A})).$

###### Proof.

We achieve the guarantee of this theorem via Algorithm 1. Note that in 4 of
Algorithm 1, the rows of $\mathbf{A}$ are partitioned into randomly created
$n/\alpha$ blocks. Suppose $\mathbf{a}_{i}$, the $i^{\mathrm{th}}$ row of
$\mathbf{A}$, falls into the block $\mathbf{B}_{\ell}$. Further, in 11, the
rows from $\mathbf{B}_{\ell}$ are mapped to those in $\mathbf{P}$ with row
indices in the set $J$, from which we compute the $i^{\mathrm{th}}$ $\ell_{1}$
sensitivity estimate as $\widetilde{\boldsymbol{\sigma}}_{i}=\max_{j\in
J}\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})$. Because of
$\mathbf{S}\mathbf{A}$ being an $\ell_{1}$ subspace embedding for
$\mathbf{A}$, we have that for all $\mathbf{x}\in\mathbb{R}^{d}$,

$\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}=\Theta(\|\mathbf{A}\mathbf{x}\|_{1}).$

We use Section 3 below to establish the claimed bounds. Let
$\mathbf{x}^{\star}=\arg\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}$
be the vector that realizes the $i^{\mathrm{th}}$ $\ell_{1}$ sensitivity of
$\mathbf{A}$. Further, let the $j^{\mathrm{th}}$ row of $\mathbf{P}$ be
$\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}$. Then, with a probability of at
least $1/2$, we have

$\displaystyle\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})$
$\displaystyle=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}|}{\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}}\geq\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}^{\star}|}{\|\mathbf{S}\mathbf{A}\mathbf{x}^{\star}\|_{1}}=\Theta(1)\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}^{\star}|}{\|\mathbf{A}\mathbf{x}^{\star}\|_{1}}\geq\Theta(1)\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}|}{\|\mathbf{A}\mathbf{x}^{\star}\|_{1}}=\Theta(1)\boldsymbol{\sigma}_{1}(\mathbf{a}_{i}),$
(3.3)

where the first step is by definition of
$\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})$, the second
is by evaluating the function being maximized at a specific choice of
$\mathbf{x}$, the third step is by applying Section 3, and the final step is
by definition of $\mathbf{x}^{\star}$ and
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})$. For the fourth step, we use the
fact that
$|\mathbf{r}_{k}^{\ell}\mathbf{B}_{\ell}\mathbf{x}^{\star}|=|\mathbf{r}_{k,i}^{\ell}\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}+\sum_{j\neq
i}\mathbf{r}_{k,j}^{\ell}(\mathbf{B}_{\ell}\mathbf{x}^{\star})_{j}|\geq|\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}|$
with a probability of at least $1/2$ since the vector $\mathbf{r}_{k}^{\ell}$
has coordinates that are $+1$ or $-1$ with equal probability. By a union bound
over the $|J|$ independent rows that block $\mathbf{B}_{\ell}$ is mapped to,
we establish the claimed lower bound in Equation 3.3 with a probability of at
least $0.9$. To show an upper bound on
$\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})$, we observe
that

$\displaystyle\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})$
$\displaystyle=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}|}{\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}}\leq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{B}_{\ell}\mathbf{x}\|_{1}}{\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}}=\Theta(1)\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{B}_{\ell}\mathbf{x}\|_{1}}{\|\mathbf{A}\mathbf{x}\|_{1}}\leq\Theta(1)\sum_{j:\mathbf{a}_{j}\in\mathbf{B}_{\ell}}\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{j}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}},$
(3.4)

where the second step is by Hölder inequality and the fact that the entries of
$\mathbf{r}_{k}^{(\ell)}$ are all bounded between $-1$ and $+1$, and the third
step is by Section 3. The final term of the above chain of inequalities may be
rewritten as
$\Theta(1)\sum_{j:\mathbf{a}_{j}\in\mathbf{B}_{\ell}}\boldsymbol{\sigma}_{1}(\mathbf{a}_{j})$,
by the definition of $\mathbf{B}_{\ell}$. Because $\mathbf{B}_{\ell}$ is a
group of $\alpha$ rows selected uniformly at random out of $n$ rows and
contains the row $\mathbf{a}_{i}$, we have:

$\mathbb{E}\left\\{\sum_{\begin{subarray}{c}j:\mathbf{a}_{j}\in\mathbf{B}_{\ell},\\\
j\neq
i\end{subarray}}\boldsymbol{\sigma}_{1}(\mathbf{a}_{j})\right\\}=\frac{\alpha-1}{n-1}\sum_{j\neq
i}\boldsymbol{\sigma}_{1}(\mathbf{a}_{j}).$

Therefore, Markov inequality gives us that with a probability of at least
$0.9$, we have

$\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{p}_{j})\leq\boldsymbol{\sigma}_{1}(\mathbf{a}_{j})+\mathfrak{S}_{1}(\mathbf{A})O(\tfrac{\alpha}{n}).$

The median trick then establishes the bounds in Equation 3.3 and Section 3
with high probability. The claimed runtime is obtained by summing the cost of
constructing $\mathbf{S}\mathbf{A}$ and $O(n/\alpha)$ computations of
$\ell_{1}$ sensitivities with respect to $\mathbf{S}\mathbf{A}$, for which we
may use Fact 2.4. ∎

As seen in Section C.1, our techniques described above also generalize to the
$p\geq 1$ case. Specifically, in Theorem C.1, we show that reducing $\ell_{p}$
sensitivity calculations by an $\alpha$ factor gives an approximation
guarantee of the form
$O(\alpha^{p-1}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}))+\tfrac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A}))$.

######  Remark 3.2.

The sensitivities returned by our algorithm are approximate, with relative and
additive error, but are useful because they preserve $\ell_{p}$ regression
approximation guarantees while increasing total sample complexity by only a
small $\textrm{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 example, consider the case $p>2$. Note that by Theorem $1.5$ of
[24], the sample complexity using approximate sensitivities is
$\widetilde{O}(\alpha^{2p-2}\mathfrak{S}_{p}^{2-2/p}(\mathbf{A}))$, using the
$\alpha^{p}$ factor blowup via Theorem C.1, that using true sensitivities is
$\widetilde{O}(\mathfrak{S}_{p}^{2-2/p}(\mathbf{A}))$, and that using
$\ell_{p}$ Lewis weights is $O(d^{p/2})$. Suppose also that the total
sensitivity $\mathfrak{S}_{p}(\mathbf{A})=d$. Further assume that $n=d^{10}$
and $\alpha=n^{\frac{1}{10p}}=d^{1/p}$; then the dimension-dependence of
sample complexities given by our approximate sensitivities, true
sensitivities, and Lewis weights are, respectively,
$\widetilde{O}(d^{4}),\widetilde{O}(d^{2})$, and $\widetilde{O}(d^{p/2})$
(ignoring small factors in the exponent that do not affect the asymptotic
analysis). Thus, our approximate sensitivities provide a tradeoff in
computational cost and sample complexity between using _true_ sensitivities
and Lewis weights.

###  3.1 Estimating the sum of $\ell_{p}$ sensitivities

In this section, we present a one-shot algorithm that provides a constant-
approximation to the total $\ell_{p}$ sensitivity for all $p\geq 1$. This
algorithm is based on importance sampling using $\ell_{p}$ Lewis weights. For
our desired guarantees, we crucially need the approximation bounds of
sensitivities via $\ell_{p}$ Lewis weights, as provided by Fact 3.3, efficient
computation of $\ell_{p}$ Lewis weights as provided by [23] for $p<4$ and by
[37] for $p\geq 4,$ and the fact that a random $\ell_{p}$ sampling matrix
$\mathbf{S}$ (cf. Fact 3.4) yields $\mathbf{S}\mathbf{A}$, a constant-
approximation subspace embedding to $\mathbf{A}$ [24, Theorem $1.8$].

Algorithm 2 Approximating the Sum of $\ell_{p}$-Sensitivities for $p\geq 1$

1:Input: Matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, approximation factor
$\gamma\in(0.01,1)$, and a scalar $p\geq 1$

2:Output: A positive scalar that satisfies, with probability $0.99$, that
$\mathfrak{S}_{p}(\mathbf{A})\leq\widehat{s}\leq(1+O(\gamma))\mathfrak{S}_{p}(\mathbf{A})$

3:Compute all $\ell_{p}$ Lewis weights of $\mathbf{A}$, denoting the
$i^{\mathrm{th}}$ weight by $\mathbf{w}_{p}(\mathbf{a}_{i})$

4:Define the sampling vector $v\in\mathbb{R}^{n}_{\geq 0}$ by
$v_{i}=\frac{\mathbf{w}_{p}(\mathbf{a}_{i})}{d}$ for all $i\in[n]$

5:Sample $m=O(d^{|1-p/2|})$ rows with replacement, picking the
$i^{\mathrm{th}}$ row with a probability of $v_{i}$

6:Construct a random $\ell_{p}$ sampling matrix $\mathbf{S}_{p}\mathbf{A}$
with probabilities $\\{v_{i}\\}_{i=1}^{n}$ (cf. Fact 3.4)

7:For each sampled row $i_{j}$, $j\in[m]$, compute
$r_{j}=\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i_{j}})}{v_{i_{j}}}$

8:Return $\tfrac{1}{m}\sum_{j=1}^{m}r_{j}$

######  Fact 3.3 (Lewis Weights bound Sensitivities, [44, 83]).

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and $p\geq 1$, the
$\ell_{p}$ sensitivity $\boldsymbol{\sigma}_{p}^{\mathbf{A}}(\mathbf{a}_{i})$
and $\ell_{p}$ Lewis weight $\mathbf{w}_{p}(\mathbf{a}_{i})$ of the row
$\mathbf{a}_{i}$ satisfy, for all $i\in[n]$,

$\boldsymbol{\sigma}_{p}^{\mathbf{A}}(\mathbf{a}_{i})\leq
d^{\,\max(0,p/2-1)}\,\mathbf{w}_{p}(\mathbf{a}_{i}).$

######  Fact 3.4 (Sampling via Lewis Weights [23, 35]).

Given $\mathbf{A}\in\mathbb{R}^{n\times d}$ and $p>0$. Consider a random
diagonal $\ell_{p}$ sampling matrix $\mathbf{S}\in\mathbb{R}^{n\times n}$ with
sampling probabilities $\\{p_{i}\\}$ proportional to the $\ell_{p}$ Lewis
weight of $\mathbf{A}$, i.e., for each $i\in[n]$, the $i^{\mathrm{th}}$
diagonal entry is independently set to be

$\mathbf{S}_{i,i}=\left\\{\begin{array}[]{ll}1/p_{i}^{1/p}&\mbox{ }\text{ with
probability }p_{i}\\\ 0&\mbox{ }\text{ otherwise}\end{array}\right..$

Then, with high probability, $\mathbf{S}$ with
$O(\varepsilon^{-2}d^{\max(1,p/2)}(\log d)^{2}\log(d/\varepsilon))$ rows is an
$\ell_{p}$ subspace embedding for $\mathbf{A}$ (cf. Definition 2.3).

######  Theorem 3.5.

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and an approximation
factor $\gamma\in(0,1)$, there exists an algorithm, which returns a positive
scalar $\widehat{s}$ such that, with a probability $0.99$, we have

$\mathfrak{S}_{p}(\mathbf{A})\leq\widehat{s}\leq(1+O(\gamma))\mathfrak{S}_{p}(\mathbf{A}).$

Our algorithm’s runtime is
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\frac{1}{\gamma^{2}}\cdot
d^{|p/2-1|}\mathsf{LP}(O(d^{\max(1,p/2)},d,p)\right)$.

###### Proof.

We achieve our theorem’s guarantee via Algorithm 2. Without loss of
generality, we may assume $\mathbf{A}$ to be full rank. Then, its Lewis
weights satisfy $\sum_{i=1}^{n}\mathbf{w}_{p}(\mathbf{a}_{i})=d$. Per 4 of
Algorithm 2, our sampling distribution is chosen to be
$v_{i}=\frac{\mathbf{w}_{p}(\mathbf{a}_{i})}{d}$ for all $i\in[n]$. We sample
from $\mathbf{A}$ rows with replacement, with row $i$ picked with a
probability of $v_{i}$. From the definition of $v_{i}$ in 4 and $r_{j}$ in 7
and the fact that $\mathbf{S}_{p}\mathbf{A}$ is a constant factor $\ell_{p}$
subspace embedding of $\mathbf{A}$, our algorithm’s output satisfies the
following unbiasedness condition:

$\mathbb{E}\left(\frac{1}{m}\sum_{j\in[m]}r_{j}\right)=\mathbb{E}\left(\frac{1}{m}\sum_{j=1}^{m}\sum_{i_{j}=1}^{n}\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i_{j}})}{v_{i_{j}}}\cdot
v_{i_{j}}\right)=\mathfrak{S}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{A})\leq
2\mathfrak{S}_{p}(\mathbf{A}).$

By independence, we also have the following bound on the variance of
$\frac{1}{m}\sum_{j\in[m]}r_{j}$:

$\textrm{Var}\left(\frac{1}{m}\sum_{j\in[m]}r_{j}\right)\leq\frac{1}{m}\sum_{i=1}^{n}\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})^{2}}{v_{i}}=\frac{d}{m}\cdot\sum_{i=1}^{n}\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})^{2}}{\mathbf{w}_{p}(\mathbf{a}_{i})},$

with the final step holding by the choice of $v_{i}$ in 4. In the case that
$p\geq 2$, we have

$\textrm{Var}\left(\frac{1}{m}\sum_{j\in[m]}r_{j}\right)\leq\frac{d}{m}\cdot\sum_{i=1}^{n}\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})^{2}}{\mathbf{w}_{p}(\mathbf{a}_{i})}\leq\frac{d}{m}\cdot\sum_{i\in[n]}\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})\cdot
d^{p/2-1}\leq\frac{2d^{p/2}}{m}\mathfrak{S}_{p}(\mathbf{A}),$

where the second inequality uses Fact 3.3. Therefore, applying Chebyshev’s
inequality on $\frac{1}{m}\sum_{j\in[m]}r_{j}$ (as defined in 7) with
$m=O\left(\frac{d^{p/2}}{\mathfrak{S}_{p}(\mathbf{A})\gamma^{2}}\right)$ gives
us a $\gamma$-multiplicative accuracy in approximating
$\mathfrak{S}_{p}(\mathbf{A})$ (with the desired constant probability). For
$p\geq 2$, we additionally have [24, Theorem $1.7$] the lower bound
$\mathfrak{S}_{p}(\mathbf{A})\geq d$, which when plugged into the value of $m$
gives the claimed sample complexity. A similar argument may be applied for the
case $p<2$; specifically, we have that

$\textrm{Var}\left(\frac{1}{m}\sum_{j\in[m]}r_{j}\right)\leq\frac{d}{{m}}\cdot\sum_{i=1}^{n}\frac{\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})^{2}}{\mathbf{w}_{p}(\mathbf{a}_{i})}\leq\frac{d}{{m}}\cdot\sum_{i\in[n]}\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{a}_{i})\leq\frac{2d}{m}\mathfrak{S}_{p}(\mathbf{A}),$

where the second step is by
$\mathbf{w}_{p}(\mathbf{a}_{i})\geq\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})$
from Fact 3.3. For $p<2$, we also have $\mathfrak{S}_{p}(\mathbf{A})\geq
d^{p/2}$ from [24, Theorem $1.7$]. Applying Chebyshev’s inequality with this
fact gives a sample complexity of $m=O\left(d^{1-p/2}/\gamma^{2}\right)$. This
completes the correctness guarantee.

##### Runtime.

We first compute all $\ell_{p}$ Lewis weights of
$\mathbf{A}\in\mathbb{R}^{n\times d}$ up to a constant multiplicative
accuracy, the cost of which is
$O\left(\frac{1}{1-|1-p/2|}\log(\log(n))\right)$ leverage score computations
for $p<4$ [23] and $O(p^{3}\log(np))$ leverage score computations for $p\geq
4$ [37]. Next, we compute $m=O(d^{|1-p/2|})$ sensitivities with respect to
$\mathbf{S}_{p}\mathbf{A}$. From [23, 35], $\mathbf{S}_{p}$ has, with high
probability, $O(d^{p/2}\log d)$ rows when $p>2$ and $O(d\log d)$ rows when
$p\leq 2$. Summing the cost of computing these $m$ sensitivities and that of
computing the Lewis weights gives the claimed runtime.

∎

######  Remark 3.6.

We present in Section B.1 an alternate algorithm, Algorithm 4, for estimating
the total $\ell_{p}$ sensitivity for $p=1$. Algorithm 4 uses recursive
computations of leverage scores, in contrast to Algorithm 2 which uses Lewis
weights in a one-shot manner, and may be of independent interest.

###  3.2 Approximating the maximum $\ell_{1}$ sensitivity

In this section, we present an algorithm that approximates
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}=\max_{i\in[|\mathbf{A}|]}\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})$,
the maximum of the $\ell_{1}$ sensitivities of the rows of a matrix
$\mathbf{A}$. As alluded to in Section 1, a first approach to this problem
would be to simply estimate all $n$ sensitivities and compute their maximum.
To do better than this, one idea, inspired by the random hashing approach of
Algorithm 1, is as follows.

If the matrix has a large number of high-sensitivity rows, then, intuitively,
the appropriately scaled maximum sensitivity of a uniformly sampled subset of
these rows should approximate
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}$. Specifically, assume that
the matrix has at least $\alpha$ rows of sensitivity at least
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}/\sqrt{\alpha}$; uniformly
sample $n/\alpha$ rows and return $\widetilde{\boldsymbol{\sigma}}_{a}$, the
(appropriately scaled) maximum of
$\boldsymbol{\sigma}_{1}^{\mathbf{S}_{1}\mathbf{A}}(\mathbf{a}_{i})$, for the
sampled rows $i$. Then,
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}\leq\widetilde{\boldsymbol{\sigma}}_{a}\leq
O(\sqrt{\alpha}\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty})$ with a
constant probability. Here the upper bound guarantee is without any condition
on the number of rows with large $\ell_{1}$ 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 $\alpha$
rows each (just like in Algorithm 1). Then,
$\widetilde{\boldsymbol{\sigma}}_{b}$, the scaled maximum of the $\ell_{1}$
sensitivities of these rows satisfies
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}\leq\widetilde{\boldsymbol{\sigma}}_{b}\leq
O(\sqrt{\alpha}\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty})$. Here, it is
the _lower_ bound that comes for free (i.e., without any condition on the
number of rows with large $\ell_{1}$ sensitivities).

Despite the above strategies working for each case, there is no way to combine
these approaches without knowledge of whether the input matrix has a enough
“high-sensitivity” rows or not. We therefore avoid this approach and instead
present Algorithm 3, where we make use of $\mathbf{S}_{\infty}\mathbf{A}$, an
$\ell_{\infty}$ subspace embedding of $\mathbf{A}$. Our algorithm hinges on
the recent development [76, Theorem 1.3] on efficient construction of such an
embedding with simply a subset of $\widetilde{O}(d)$ rows of $\mathbf{A}$.

Algorithm 3 Approximating the Maximum of $\ell_{1}$-Sensitivities

1:Input: Matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$

2:Output: Scalar $\widehat{s}\in\mathbb{R}_{\geq 0}$ that satisfies
$\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}\leq\widehat{s}\leq
C\cdot\sqrt{d}\cdot\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}$

3:Compute, for $\mathbf{A}$, an $\ell_{\infty}$ subspace embedding
$\mathbf{S}_{\infty}\mathbf{A}\in\mathbb{R}^{O(d\log^{2}(d))\times d}$ such
that $\mathbf{S}_{\infty}\mathbf{A}$ is a subset of the rows of $\mathbf{A}$
[76, Theorem $1.3$]

4:Compute, for $\mathbf{A}$, an $\ell_{1}$ subspace embedding
$\mathbf{S}_{1}\mathbf{A}\in\mathbb{R}^{O(d)\times d}$

5:Return
$\sqrt{d}\|\boldsymbol{\sigma}_{1}^{\mathbf{S}_{1}\mathbf{A}}(\mathbf{S}_{\infty}\mathbf{A})\|_{\infty}$

######  Theorem 3.7.

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, there exists an
algorithm, which in time
$\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+d^{\omega+1})$, outputs a
positive scalar $\widehat{s}$ that satisfies

$\Omega(\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty})\leq\widehat{s}\leq
O(\sqrt{d}\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}).$

###### Proof.

We achieve our guarantee via Algorithm 3. First, by the $\ell_{1}$ subspace
embedding property of $\mathbf{S}_{1}\mathbf{A}$, we have for all
$\mathbf{x}\in\mathbb{R}^{d}$ that

$\|\mathbf{S}_{1}\mathbf{A}\mathbf{x}\|_{1}=\Theta(\|\mathbf{A}\mathbf{x}\|_{1}).$

Next, we set some notation. Define $\mathbf{x}^{\star}$ and
$\mathbf{a}_{i^{\star}}$ as follows:

$\mathbf{x}^{\star},i^{\star}=\arg\max_{\mathbf{x}\in\mathbb{R}^{d},i\in[|\mathbf{A}|]}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}$

Thus, $\mathbf{x}^{\star}$ is the vector that realizes the maximum sensitivity
of the matrix $\mathbf{A}$, and the row $\mathbf{a}_{i^{\star}}$ is the row of
$\mathbf{A}$ with maximum sensitivity with respect to $\mathbf{A}$. Suppose
the matrix $\mathbf{S}_{\infty}\mathbf{A}$ — which is a subset of the rows of
$\mathbf{A}$ — contains the row $\mathbf{a}_{i^{\star}}$. Then we have

$\displaystyle\max_{i:\mathbf{c}_{i}\in\mathbf{S}_{\infty}\mathbf{A}}\boldsymbol{\sigma}_{1}(\mathbf{c}_{i})$
$\displaystyle=\max_{\mathbf{x}\in\mathbb{R}^{d},\mathbf{c}_{k}\in\mathbf{S}_{\infty}\mathbf{A}}\frac{|\mathbf{c}_{k}^{\top}\mathbf{x}|}{\|\mathbf{S}_{1}\mathbf{A}\mathbf{x}\|_{1}}=\Theta(1)\max_{\mathbf{x}\in\mathbb{R}^{d},\mathbf{c}_{k}\in\mathbf{S}_{\infty}\mathbf{A}}\frac{|\mathbf{c}_{k}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}=\Theta(1)\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i^{\star}}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}},$
(3.8)

where the first step is by definition of $\ell_{1}$ sensitivity of
$\mathbf{S}_{\infty}\mathbf{A}$ with respect to $\mathbf{S}_{1}\mathbf{A}$,
the second step is by Section 3.2, and the third step by our assumption in
this case that $\mathbf{a}_{i^{\star}}\in\mathbf{S}_{\infty}\mathbf{A}$. By
definition of $\boldsymbol{\sigma}_{1}(\mathbf{a}_{i^{\star}})$ in Section
3.2, we have

$\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i^{\star}}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}=\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty}.$

Then, combining Equation 3.8 and Section 3.2 gives the guarantee in this case.
In the other case, suppose $\mathbf{a}_{i^{\star}}$ is not included in
$\mathbf{S}_{\infty}\mathbf{A}$. Then we observe that the upper bound from the
preceding inequalities still holds. For the lower bound, we observe that

$\|\boldsymbol{\sigma}_{1}^{\mathbf{S}_{1}\mathbf{A}}(\mathbf{S}_{\infty}\mathbf{A})\|_{\infty}=\max_{\mathbf{x}\in\mathbb{R}^{d},\mathbf{c}_{j}\in\mathbf{S}_{\infty}\mathbf{A}}\frac{\mathbf{c}_{j}^{\top}\mathbf{x}}{\|\mathbf{S}_{1}\mathbf{A}\mathbf{x}\|_{1}}\geq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{S}_{\infty}\mathbf{A}\mathbf{x}\|_{\infty}}{\|\mathbf{S}_{1}\mathbf{A}\mathbf{x}\|_{1}}=\Theta(1)\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{S}_{\infty}\mathbf{A}\mathbf{x}\|_{\infty}}{\|\mathbf{A}\mathbf{x}\|_{1}},$

where the the second step is by choosing a specific vector in the numerator
and the third step uses Section 3.2. We further have,

$\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{S}_{\infty}\mathbf{A}\mathbf{x}\|_{\infty}}{\|\mathbf{A}\mathbf{x}\|_{1}}\geq\frac{\|\mathbf{S}_{\infty}\mathbf{A}\mathbf{x}^{\star}\|_{\infty}}{\|\mathbf{A}\mathbf{x}^{\star}\|_{1}}\geq\frac{\|\mathbf{A}\mathbf{x}^{\star}\|_{\infty}}{\sqrt{d}\|\mathbf{A}\mathbf{x}^{\star}\|_{1}}\geq\frac{|\mathbf{a}_{i^{\star}}^{\top}\mathbf{x}^{\star}|}{\sqrt{d}\|\mathbf{A}\mathbf{x}^{\star}\|_{1}}=\frac{1}{\sqrt{d}}\|\boldsymbol{\sigma}_{1}(\mathbf{A})\|_{\infty},$

where the first step is by choosing $\mathbf{x}=\mathbf{x}^{\star}$, the
second step is by the distortion guarantee of $\ell_{\infty}$ subspace
embedding [35], and the final step is by definition of
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i^{\star}})$. Combining Section C.2 and
Section C.2 gives the claimed lower bound on
$\|\boldsymbol{\sigma}_{1}^{\mathbf{S}_{1}\mathbf{A}}(\mathbf{S}_{\infty}\mathbf{A})\|_{\infty}$.
The runtime follows from the computational cost and row dimension of
$\mathbf{S}_{\infty}\mathbf{A}$ from [35] and the cost of computing $\ell_{1}$
sensitivities with respect to $\mathbf{S}_{1}\mathbf{A}$ as per Fact 2.4. ∎

##  4 Numerical experiments

We demonstrate our fast sensitivity approximations on multiple real-world
datasets in the UCI Machine Learning Dataset Repository [48], such as the wine
and fires datasets, for different $p$ and varying approximation parameters
$\alpha$. We focus on the wine dataset here, with full details in Appendix E.

##### Experimental Setup.

For each dataset and each $p$, we first apply the Lewis weight subspace
embedding to compute a smaller matrix $\mathbf{S}\mathbf{A}$. Then, we run our
sensitivity computation Algorithm 1 and measure the average and maximum
absolute log ratio $|\log(\sigma_{\text{approx}}/\sigma_{\text{true}})|$ to
capture the relative error of the sensitivity approximation. Note that this is
stricter than our guarantees, which give only an $O(\alpha)$ 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(\alpha^{p})$, which is
the worst case additive error approximation although it does provide relative
error guarantees with respect to 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 easy-to-approximate
sensitivities with total $\ell_{p}$ sensitivity about $2$-$5$ times lower than
the theoretical upper bound. Figure 1 shows that the average absolute log
ratios for approximating the $\ell_{p}$ sensitivities are much smaller than
those suggested by the theoretical upper bound, even for large $\alpha$.
Specifically, when $\alpha=40$, we incur only a $16$x accuracy deterioration
in exchange for a $40$x faster algorithm. This is significantly better than
the worst-case accuracy guarantee which for $p=3$ would be
$\alpha^{p}=40^{3}$.

More importantly, we find that empirical total sensitivities are much lower
than the theoretical 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 to an overestimate within a factor of
$1.3$, in less than a quarter of the time of the brute-force algorithm. Our
observation generalizes to other real-world datasets (see Appendix E).

Figure 1: Average absolute log ratios for all $\ell_{p}$ sensitivity
approximations for wine, with the theoretical additive bound (dashed).

$p$ |  Total Sensitivity Upper Bound | Brute-Force/True | Approximation | Brute-Force Runtime (s) | Approximate Runtime (s)
---|---|---|---|---|---
1 | 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 $\ell_{p}$ sensitivity sampling; and our NeurIPS
reviewers for their time, effort, and many constructive suggestions. Swati
gratefully acknowledges funding for her student researcher position from
Google Research and research 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

  * Robbins and Monro [1951] Herbert Robbins and Sutton Monro.  A stochastic approximation method.  _The annals of mathematical statistics_ , 1951.
  * Bottou and Cun [2003] Léon Bottou and Yann Cun.  Large scale online learning.  _Advances in neural information processing systems_ , 16, 2003.
  * Zhang [2004] Tong Zhang.  Solving large scale linear prediction problems using stochastic gradient descent algorithms.  In _Proceedings of the twenty-first international conference on Machine learning_ , page 116, 2004.
  * Bottou [2012] Léon Bottou.  Stochastic gradient descent tricks.  In _Neural networks: Tricks of the trade_ , pages 421–436. Springer, 2012.
  * Chang and Lin [2011] 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://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * Roux et al. [2012] Nicolas Roux, Mark Schmidt, and Francis Bach.  A stochastic gradient method with an exponential convergence rate for finite training sets.  _Advances in neural information processing systems_ , 25, 2012.
  * Shalev-Shwartz and Zhang [2013] Shai Shalev-Shwartz and Tong Zhang.  Stochastic dual coordinate ascent methods for regularized loss minimization.  _Journal of Machine Learning Research_ , 14(2), 2013.
  * Johnson and Zhang [2013] Rie Johnson and Tong Zhang.  Accelerating stochastic gradient descent using predictive variance reduction.  _Advances in neural information processing systems_ , 26, 2013.
  * Mahdavi et al. [2013] Mehrdad Mahdavi, Lijun Zhang, and Rong Jin.  Mixed optimization for smooth functions.  _Advances in neural information processing systems_ , 26, 2013.
  * Defazio et al. [2014] Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.  Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.  _Advances in neural information processing systems_ , 27, 2014.
  * Mairal [2015] Julien Mairal.  Incremental majorization-minimization optimization with application to large-scale machine learning.  _SIAM Journal on Optimization_ , 25(2), 2015.
  * Allen-Zhu and Yuan [2016] Zeyuan Allen-Zhu and Yang Yuan.  Improved svrg for non-strongly-convex or sum-of-non-convex objectives.  In _International conference on machine learning_ , 2016.
  * Hazan and Luo [2016] Elad Hazan and Haipeng Luo.  Variance-reduced and projection-free stochastic optimization.  In _International Conference on Machine Learning_ , 2016.
  * Schmidt et al. [2017] Mark Schmidt, Nicolas Le Roux, and Francis Bach.  Minimizing finite sums with the stochastic average gradient.  _Mathematical Programming_ , 162(1), 2017.
  * Langberg and Schulman [2010] Michael Langberg and Leonard J Schulman.  Universal $\varepsilon$-approximators for integrals.  In _Proceedings of the twenty-first annual ACM-SIAM symposium on Discrete Algorithms_ , 2010.
  * Feldman and Langberg [2011] 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.
  * Needell et al. [2014] Deanna Needell, Rachel Ward, and Nati Srebro.  Stochastic gradient descent, weighted sampling, and the randomized kaczmarz algorithm.  In 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://proceedings.neurips.cc/paper_files/paper/2014/file/f29c21d4897f78948b91f03172341b7b-Paper.pdf.
  * Zhao and Zhang [2015] Peilin Zhao and Tong Zhang.  Stochastic optimization with importance sampling for regularized loss minimization.  In _Proceedings of the 32nd International Conference on Machine Learning_ , pages 1–9, 2015.
  * Katharopoulos and Fleuret [2017] Angelos Katharopoulos and François Fleuret.  Biased importance sampling for deep neural network training.  _arXiv preprint arXiv:1706.00043_ , 2017.
  * Johnson and Guestrin [2018] Tyler B. Johnson and Carlos Guestrin.  Training deep models faster with robust, approximate importance sampling.  In _Advances in Neural Information Processing Systems 31_ , 2018.
  * Mindermann et al. [2022] 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.  Prioritized training on points that are learnable, worth learning, and not yet learnt.  In _International Conference on Machine Learning_ , 2022.
  * Braverman et al. [2016] Vladimir Braverman, Dan Feldman, Harry Lang, Adiel Statman, and Samson Zhou.  New frameworks for offline and streaming coreset constructions.  _arXiv preprint arXiv:1612.00889_ , 2016.
  * Cohen and Peng [2015] 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.
  * Woodruff and Yasuda [2023a] David Woodruff and Taisuke Yasuda.  Sharper bounds for $\ell_{p}$ sensitivity sampling.  In _International Conference on Machine Learning_ , pages 37238–37272. PMLR, 2023a.
  * Meyer et al. [2022] 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.
  * Varadarajan and Xiao [2012] 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.
  * Braverman et al. [2021] Vladimir Braverman, Dan Feldman, Harry Lang, Adiel Statman, and Samson Zhou.  Efficient coreset constructions via sensitivity sampling.  _Proceedings of Machine Learning Research_ , 157(2021), 2021.
  * Tukan et al. [2022a] 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_ , 2022a.
  * Bachem et al. [2015] Olivier Bachem, Mario Lucic, and Andreas Krause.  Coresets for nonparametric estimation-the case of dp-means.  In _International Conference on Machine Learning_ , 2015.
  * Lucic et al. [2016] 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.
  * Huggins et al. [2016] Jonathan Huggins, Trevor Campbell, and Tamara Broderick.  Coresets for scalable bayesian logistic regression.  _Advances in Neural Information Processing Systems_ , 29, 2016.
  * Munteanu et al. [2018] Alexander Munteanu, Chris Schwiegelshohn, Christian Sohler, and David Woodruff.  On coresets for logistic regression.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Drineas et al. [2006] Petros Drineas, Michael W Mahoney, and Shan Muthukrishnan.  Sampling algorithms for $\ell_{2}$ regression and applications.  In _Proceedings of the seventeenth annual ACM-SIAM symposium on Discrete algorithm_ , 2006.
  * Mahoney et al. [2011] Michael W Mahoney et al.  Randomized algorithms for matrices and data.  _Foundations and Trends® in Machine Learning_ , 3(2), 2011.
  * Woodruff and Yasuda [2023b] David P Woodruff and Taisuke Yasuda.  Online lewis weight sampling.  In _Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , 2023b.
  * Cohen et al. [2015] 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.
  * Fazel et al. [2022] 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 on Discrete Algorithms (SODA)_ , 2022.
  * Chierichetti et al. [2017] Flavio Chierichetti, Sreenivas Gollapudi, Ravi Kumar, Silvio Lattanzi, Rina Panigrahy, and David P Woodruff.  Algorithms for $\ell_{p}$ low-rank approximation.  In _International Conference on Machine Learning_ , 2017.
  * Candes and Tao [2005] Emmanuel J Candes and Terence Tao.  Decoding by linear programming.  _IEEE transactions on information theory_ , 51(12), 2005.
  * Alamgir and Luxburg [2011] Morteza Alamgir and Ulrike Luxburg.  Phase transition in the family of $p$-resistances.  _Advances in neural information processing systems_ , 24, 2011.
  * Calder [2019] Jeff Calder.  Consistency of lipschitz learning with infinite unlabeled data and finite labeled data.  _SIAM Journal on Mathematics of Data Science_ , 1(4), 2019.
  * Flores et al. [2022] Mauricio Flores, Jeff Calder, and Gilad Lerman.  Analysis and algorithms for $\ell_{p}$-based semi-supervised learning on graphs.  _Applied and Computational Harmonic Analysis_ , 60, 2022.
  * Elmoataz et al. [2015] Abderrahim Elmoataz, Matthieu Toutain, and Daniel Tenbrinck.  On the $p$-laplacian and $\infty$-laplacian on graphs with applications in image and data processing.  _SIAM Journal on Imaging Sciences_ , 8(4), 2015.
  * Clarkson et al. [2019] Kenneth Clarkson, Ruosong Wang, and David Woodruff.  Dimensionality reduction for tukey regression.  In _International Conference on Machine Learning_ , 2019.
  * Dwork et al. [2014] Cynthia Dwork, Aaron Roth, et al.  The algorithmic foundations of differential privacy.  _Foundations and Trends® in Theoretical Computer Science_ , 9(3–4):211–407, 2014.
  * Adil et al. [2019] 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.
  * Jambulapati et al. [2022] Arun Jambulapati, Yang P Liu, and Aaron Sidford.  Improved iteration complexities for overconstrained p-norm regression.  In _Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing_ , 2022.
  * Asuncion and Newman [2007] Arthur Asuncion and David Newman.  Uci machine learning repository, 2007.
  * Raj et al. [2020] Anant Raj, Cameron Musco, and Lester Mackey.  Importance sampling via local sensitivity.  In _International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Maalouf et al. [2020] 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–2061, 2020.
  * Wang et al. [2020] 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.
  * Kong et al. [2021] Dingwen Kong, Ruslan Salakhutdinov, Ruosong Wang, and Lin F Yang.  Online sub-sampling for reinforcement learning with general function approximation.  _arXiv preprint arXiv:2106.07203_ , 2021.
  * Liebenwein et al. [2019] 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.
  * Tukan et al. [2022b] Murad Tukan, Loay Mualem, and Alaa Maalouf.  Pruning neural networks via coresets and convex geometry: Towards no assumptions.  _Advances in Neural Information Processing Systems_ , 2022b.
  * Mussay et al. [2021] 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–7841, 2021.
  * Lewis [1978] D Lewis.  Finite dimensional subspaces of $\ell_{p}$.  _Studia Mathematica_ , 1978.
  * Schechtman [1987] Gideon Schechtman.  More on embedding subspaces of $\ell_{p}$ in $\ell^{n}_{r}$.  _Compositio Mathematica_ , 1987.
  * Bourgain et al. [1989] Jean Bourgain, Joram Lindenstrauss, and Vitali Milman.  Approximation of zonoids by zonotopes.  1989\.
  * Talagrand [1990] Michel Talagrand.  _Proceedings of the American Mathematical Society_ , 1990.
  * Ledoux and Talagrand [1991] Michel Ledoux and Michel Talagrand.  _Probability in Banach Spaces: isoperimetry and processes_ , volume 23.  Springer Science & Business Media, 1991.
  * Schechtman and Zvavitch [2001] Gideon Schechtman and Artem Zvavitch.  Embedding subspaces of $\ell_{p}$ into $\ell^{n}_{p}$, $0<p<1$.  _Mathematische Nachrichten_ , 2001.
  * Drineas et al. [2012] 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.
  * Li et al. [2013] Mu Li, Gary L Miller, and Richard Peng.  Iterative row sampling.  In _2013 IEEE 54th Annual Symposium on Foundations of Computer Science_ , 2013.
  * Sohler and Woodruff [2018] Christian Sohler and David P Woodruff.  Strong coresets for k-median and subspace approximation: Goodbye dimension.  In _2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS)_ , 2018.
  * Chhaya et al. [2020] Rachit Chhaya, Jayesh Choudhari, Anirban Dasgupta, and Supratim Shit.  Streaming coresets for symmetric tensor factorization.  In _International Conference on Machine Learning_ , 2020.
  * Durfee et al. [2018] David Durfee, Kevin A Lai, and Saurabh Sawlani.  $\ell_{1}$ regression using lewis weights preconditioning and stochastic gradient descent.  In _Conference On Learning Theory_ , 2018.
  * Li et al. [2020a] 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_ , 2020a.
  * Braverman et al. [2020] 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.
  * Laddha et al. [2020] 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.
  * Li et al. [2020b] Yi Li, Ruosong Wang, Lin Yang, and Hanrui Zhang.  Nearly linear row sampling algorithm for quantile regression.  In _International Conference on Machine Learning_ , 2020b.
  * Mahankali and Woodruff [2021] Arvind V Mahankali and David P Woodruff.  Optimal $\ell_{1}$ column subset selection and a fast ptas for low rank approximation.  In _Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA)_ , 2021.
  * Mai et al. [2021] Tung Mai, Cameron Musco, and Anup Rao.  Coresets for classification–simplified and strengthened.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Chen and Derezinski [2021] Xue Chen and Michal Derezinski.  Query complexity of least absolute deviation regression via robust uniform convergence.  In _Conference on Learning Theory_ , 2021.
  * Parulekar et al. [2021] Aditya Parulekar, Advait Parulekar, and Eric Price.  L1 regression with lewis weights subsampling.  _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques_ , 2021.
  * Foster [1953] Frederic G Foster.  On the stochastic matrices associated with certain queuing processes.  _The Annals of Mathematical Statistics_ , 24(3), 1953.
  * Woodruff and Yasuda [2022] David P Woodruff and Taisuke Yasuda.  High-dimensional geometric streaming in polynomial space.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 732–743, 2022.
  * Lee and Sidford [2019] Yin Tat Lee and Aaron Sidford.  Solving linear programs with sqrt (rank) linear system solves.  _arXiv preprint arXiv:1910.08033_ , 2019.
  * Musco et al. [2021] David P Woodruff Cameron Musco, Christopher Musco, David P Woodruff, and Taisuke Yasuda.  Active sampling for linear regression beyond the l2 norm.  _CoRR_ , 2021.
  * Woodruff et al. [2014] David P Woodruff et al.  Sketching as a tool for numerical linear algebra.  _Foundations and Trends® in Theoretical Computer Science_ , 10(1–2), 2014.
  * Rudelson and Vershynin [2007] Mark Rudelson and Roman Vershynin.  Sampling from large matrices: An approach through geometric functional analysis.  _Journal of the ACM (JACM)_ , 2007.
  * Sohler and Woodruff [2011] Christian Sohler and David P Woodruff.  Subspace embeddings for the $\ell\\_1$-norm with applications.  In _Proceedings of the forty-third annual ACM symposium on Theory of computing_ , pages 755–764, 2011.
  * Woodruff and Zhang [2013] David Woodruff and Qin Zhang.  Subspace embeddings and $\ell\\_p$-regression using exponential random variables.  In _Conference on Learning Theory_ , pages 546–567. PMLR, 2013.
  * Musco et al. [2022] Cameron Musco, Christopher Musco, David P Woodruff, and Taisuke Yasuda.  Active linear regression for $\ell_{p}$ norms and beyond.  In _2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS)_ , 2022.
  * Cohen et al. [2021] 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.
  * Vershynin [2020] Roman Vershynin.  _High-dimensional probability: An introduction with applications in data science_.  Cambridge university press, 2020.

##  Appendix A Omitted proofs: general technical results

See Fact 2.2

###### Proof.

The proof of this claim follows from a simple application of standard norm
inequalities, and we present one here for completeness. For any
$\mathbf{u}\in\mathbb{R}^{n}$, we have
$\|\mathbf{u}\|_{2}\leq\|\mathbf{u}\|_{1}\leq\sqrt{n}\|\mathbf{u}\|_{2}$.
Therefore, for any $\mathbf{x}\in\mathbb{R}^{d}$, we have

$\frac{|\mathbf{x}^{\top}\mathbf{a}_{i}|}{\sqrt{n}\|\mathbf{A}\mathbf{x}\|_{2}}\leq\frac{|\mathbf{x}^{\top}\mathbf{a}_{i}|}{\|\mathbf{A}\mathbf{x}\|_{1}}\leq\frac{|\mathbf{x}^{\top}\mathbf{a}_{i}|}{\|\mathbf{A}\mathbf{x}\|_{2}}.$

Let $\mathbf{x}_{\tau}$ be the vector that realizes the $i^{\mathrm{th}}$
leverage score, and $\mathbf{x}_{\sigma}$ be the vector that realizes the
$i^{\mathrm{th}}$ $\ell_{1}$ sensitivity (i.e. they are the maximizers of
sensitivity). Then, we may conclude from Appendix A that

$\sqrt{\frac{|\mathbf{x}_{\tau}^{\top}\mathbf{a}_{i}|^{2}}{n\|\mathbf{A}\mathbf{x}_{\tau}\|_{2}^{2}}}\leq\frac{|\mathbf{x}_{\tau}^{\top}\mathbf{a}_{i}|}{\|\mathbf{A}\mathbf{x}_{\tau}\|_{1}}\leq\frac{|\mathbf{x}_{\sigma}^{\top}\mathbf{a}_{i}|}{\|\mathbf{A}\mathbf{x}_{\sigma}\|_{1}}\leq\sqrt{\frac{|\mathbf{x}_{\sigma}^{\top}\mathbf{a}_{i}|^{2}}{\|\mathbf{A}\mathbf{x}_{\sigma}\|_{2}^{2}}}\leq\sqrt{\frac{|\mathbf{x}_{\tau}^{\top}\mathbf{a}_{i}|^{2}}{\|\mathbf{A}\mathbf{x}_{\tau}\|_{2}^{2}}},$

which gives the claimed result. ∎

See Fact 2.4

###### Proof.

Given an $n\times d$ matrix $\mathbf{A}$, we first compute
$\mathbf{S}\mathbf{A}$, its $\ell_{1}$ subspace embedding of size $O(d)\times
d$, and compute the $i^{\mathrm{th}}$ sensitivity of $\mathbf{A}$ as follows.

$\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{i})=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\Theta(\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1})}\approx_{O(1)}\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{a}_{i}),$

where the second step is by the definition of $\ell_{1}$ subspace embedding
(cf. Definition 2.3). To compute
$\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{a}_{i})$, we need to
compute
$\min_{\mathbf{a}_{i}^{\top}\mathbf{x}=1}\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}$,
which, by introducing a new variable for each of the rows of
$\mathbf{S}\mathbf{A}\mathbf{x}$, can be transformed into a linear program
with $O(d)$ variables and $d$ constraints. To see the claim on runtime, the
cost of computing the subspace embedding is
$\mathop{\mathbf{nnz}}(\mathbf{A})$. Having computed this once, we can then
solve $k$ linear programs, each at cost $d^{\omega}$ (which is the current
fastest LP solver [84]). Putting together these two costs gives the claimed
runtime. ∎

######  Lemma A.1.

Given positive numbers $a_{1},a_{2},\dotsc,a_{m}$, let
$r:=\tfrac{\max_{i\in[m]}a_{i}}{\min_{i\in[m]}a_{i}}$ and
$A^{(\text{true})}:=\sum_{i\in[m]}a_{i}$. Suppose we sample, uniformly at
random with replacement, a set $S\in[m]$ of these numbers, and let
$A^{(\text{est})}:=\tfrac{m}{|S|}\cdot\sum_{i:a_{i}\in S}a_{i}$. Then, if
$|S|\geq
10\left(\frac{r(1+\gamma)}{\gamma^{2}}\log\left(\frac{1}{\delta}\right)\right)$
for a large enough absolute constant $C$, we can ensure with at least a
probability of $1-\delta$, that
$|A^{(\text{est})}-A^{(\text{true})}|\leq\gamma A^{(\text{true})}$.

###### Proof.

By construction, the expected value of a sample $a_{i}$ in $S$ is
$\tfrac{1}{m}A^{(\text{true})}$. Denoting by $\mathcal{P}$ the uniform
distribution over $a_{1},a_{2},\dotsc,a_{m}$, and let $\sigma(\mathcal{P})$ be
the standard deviation of this distribution $\mathcal{P}$. Then, we can apply
Bernstein’s inequality [85, Theorem $2.8.1$] to see that the absolute error
$|A^{(\text{est})}-A^{(\text{true})}|$ satisfies the following guarantee for
all $t>0$:

$\displaystyle\Pr\left\\{|A^{(\text{est})}-A^{(\text{true})}|\geq t\right\\}$
$\displaystyle=\Pr\left\\{\left|\sum_{i:a_{i}\in
S}\left(a_{i}-\tfrac{1}{m}A^{(\text{true})}\right)\right|\geq
t\cdot\tfrac{|S|}{m}\right\\}$ $\displaystyle\leq
e^{-\left\\{\tfrac{\tfrac{1}{2}t^{2}\cdot\tfrac{|S|^{2}}{m^{2}}}{|S|\cdot\sigma(\mathcal{P})^{2}+\tfrac{1}{3}t\cdot\frac{|S|}{m}\cdot\max_{i\in[m]}a_{i}}\right\\}}$
$\displaystyle=e^{-\left\\{\tfrac{\tfrac{1}{2}t^{2}\cdot|S|}{m^{2}\cdot\sigma(\mathcal{P})^{2}+\tfrac{1}{3}\cdot
m\cdot t\cdot\max_{i\in[m]}a_{i}}\right\\}}.$ (A.2)

Since each of the $a_{i}$s is positive, we note that

$\sigma(\mathcal{P})^{2}=\frac{1}{m}\sum_{i\in[m]}\left(a_{i}-\frac{A^{\text{(true)}}}{m}\right)^{2}\leq\frac{1}{m}\sum_{i\in[m]}a_{i}^{2}\leq\frac{1}{m}{\cdot\max_{i\in[m]}a_{i}\cdot\sum_{i\in[m]}a_{i}}=\frac{1}{m}\cdot\max_{i\in[m]}a_{i}\cdot
A^{(\text{true})}.$

Combining Equation A.2 and Appendix A for $t=\gamma A^{\text{(true)}}$ implies
that

$\Pr\left\\{|A^{(\text{est})}-A^{(\text{true})}|\geq\gamma
A^{(\text{true})}\right\\}\leq
e^{-\left\\{\tfrac{\tfrac{1}{2}\gamma^{2}A^{(\text{true})}|S|}{m(1+\tfrac{1}{3}\cdot\gamma)\cdot\max_{i\in[m]}a_{i}}\right\\}}\leq\delta$

holds for the choice of $|S|\geq
10\left(\frac{m}{\gamma^{2}}(1+\gamma)\frac{\max_{i\in[m]}a_{i}}{A^{(\text{true})}}\log\left(\frac{1}{\delta}\right)\right)$
yields the claimed error guarantee. ∎

##  Appendix B Omitted proofs: $\ell_{1}$ sensitivities

###  B.1 Estimating the sum of $\ell_{1}$ sensitivities

In this section, we present a randomized algorithm that approximates the total
sensitivity $\mathfrak{S}_{1}(\mathbf{A})$, up to a $\gamma$-multiplicative
approximation for some $\gamma\in(0,1)$ that is not too small. This algorithm
is significantly less general than Algorithm 2 (which holds for all $p\geq 1$)
but it uses extremely simple facts about leverage scores and what we think is
a new recursion technique in this context, which we believe could be of
potential independent interest.

######  Theorem B.1.

Given a matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and an approximation
factor $\gamma\in(0,1)$ such that
$\gamma\geq\Omega\left(\tfrac{1}{n^{3}}\right)$, there exists an algorithm,
which in time
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\tfrac{d^{\omega}}{\gamma^{2}}\cdot\max\left(\sqrt{d},\tfrac{1}{\gamma^{2}}\right)\right)$,
returns a positive scalar $\widehat{s}$ such that, with a probability of
$0.99$, we have

$\mathfrak{S}_{1}(\mathbf{A})\leq\widehat{s}\leq(1+O(\gamma))\mathfrak{S}_{1}(\mathbf{A}).$

Algorithm 4 Approximating the Sum of $\ell_{1}$-Sensitivities

1:Input: Matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$ and approximation factor
$\gamma\in(0,1)$

2:Output: A positive scalar that satisfies, with probability $0.99$, that
$\mathfrak{S}_{1}(\mathbf{A})\leq\widehat{s}\leq(1+O(\gamma))\mathfrak{S}_{1}(\mathbf{A})$

3:Compute $\mathbf{S}\mathbf{A}\in\mathbb{R}^{O(d)\times d}$, a
$1/2$-approximate $\ell_{1}$ subspace embedding of $\mathbf{A}$ (cf.
Definition 2.3)

4:Set $\rho=O({\gamma}/{\log(\log(n+d))})$. Compute
$\mathbf{S}^{\prime}\mathbf{A}$, a $\rho$-approximate $\ell_{1}$ subspace
embedding of $\mathbf{A}$

5:Construct a matrix $\widehat{\mathbf{A}}$ comprising only those rows of
$\mathbf{A}$ with leverage scores at least $1/n^{10}$

6:Let $s$ be the output of Algorithm 5 with inputs $\widehat{\mathbf{A}}$,
$\mathbf{S}\mathbf{A}$, $\mathbf{S}^{\prime}\mathbf{A}$, $n$, and $\rho$

7:Return
$\widehat{s}:=(1+\gamma)\left(s+\frac{1}{n^{5}}(|\mathbf{A}|-|\widehat{\mathbf{A}}|)\right)$

Algorithm 5 Subroutine for Approximating the Sum of $\ell_{1}$-Sensitivities

1:Input: Matrices $\mathbf{M}\in\mathbb{R}^{r\times d}$,
$\mathbf{S}\mathbf{A}\in\mathbb{R}^{O(d)\times d}$,
$\mathbf{S}^{\prime}\mathbf{A}\in\mathbb{R}^{O(d)\times d}$, scalars $n$ and
$\rho\in(0,1)$

2:Output: A positive scalar $\widehat{s}$ that satisfies, with probability
$0.99$, that
$\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})\leq\widehat{s}\leq(1+O(\rho))\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}).$

3:Set $\mathcal{B}=\Theta(\log n)$, $D=O(\log(\log(n+d)))$,
$\delta=\frac{0.01}{\mathcal{B}^{D}}$ and
$b=\widetilde{O}\left(\frac{D^{4}}{\gamma^{2}}\max\left(\frac{D^{4}}{\gamma^{2}},\sqrt{d}\right)\right)$

4:if the number of rows of $\mathbf{M}$ is at most $b$ then

5: Let $\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})$ be the
sum of $\ell_{1}$ sensitivities of the rows of ${\mathbf{M}}$ with respect to
$\mathbf{S}^{\prime}\mathbf{A}$

6: Using Fact 2.4 on $\mathbf{M}$, compute $\widehat{s}$ such that
$\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})\leq\widehat{s}\leq(1+O(\rho))\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})$

7: Return $\widehat{s}$

8:else

9: Construct the matrix $\mathbf{C}:=\begin{bmatrix}\mathbf{M}\\\
\mathbf{S}\mathbf{A}\end{bmatrix}$. Set
$r=O(\sqrt{|\mathbf{C}|}(1+\rho)\rho^{-2}\log(\delta^{-1}))$.

10: Compute $\tau(\mathbf{C})$, the vector of leverage scores of $\mathbf{C}$

11: Divide $\left[\frac{1}{n^{20}},1\right]$ into $\mathcal{B}$ geometrically
decreasing sub-intervals $[1/2,1]$, $[1/4,1/2],\dotsc$, etc.

12: for each bucket $i\in[\mathcal{B}]$ created in the previous step do

13: Construct matrix $\mathbf{M}_{i}$ with those rows of $\mathbf{M}$ whose
leverage scores fall in the $i^{\mathrm{th}}$ bucket

14: Construct matrix $\widetilde{\mathbf{M}}_{i}$ composed of $r$ uniformly
sampled (with replacement) rows of $\mathbf{M}_{i}$

15: Compute $\widetilde{s}_{i}$, the output of Algorithm 5 with inputs
$\widetilde{\mathbf{M}}_{i}$, $\mathbf{S}\mathbf{A}$,
$\mathbf{S}^{\prime}\mathbf{A}$, $n$, and $\rho$.

16: end for

17: Return
$\widehat{s}:=(1+\rho)\sum_{i\in[\mathcal{B}]}\frac{|\mathbf{M}_{i}|}{|\widetilde{\mathbf{M}}_{i}|}\widetilde{s}_{i}$

18:end if

####  B.1.1 Proof sketch.

We achieve our guarantee via Algorithm 4 and Algorithm 5, which are based on
the following key principles. The first is Fact 2.2: the $i^{\mathrm{th}}$
sensitivity and $i^{\mathrm{th}}$ leverage score of
$\mathbf{A}\in\mathbb{R}^{n\times d}$ satisfy
$\sqrt{\frac{\tau_{i}(\mathbf{A})}{n}}\leq\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\leq\sqrt{\tau_{i}(\mathbf{A})}$.
The second idea, derived from Bernstein’s inequality (and formalized in Lemma
A.1), is: the sum of a set of positive numbers with values bounded between $a$
and $b$ can be approximated to a $\gamma$-multiplicative factor using
$\widetilde{O}((b/a)\gamma^{-2})$ uniformly sampled (and appropriately scaled)
samples. The third principle is that
$\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\approx\boldsymbol{\sigma}_{1}^{\mathbf{S}\mathbf{A}}(\mathbf{a}_{i})$
when $\mathbf{S}\mathbf{A}$ is a constant-approximation $\ell_{1}$ subspace
embedding of $\mathbf{A}$.

Based on these ideas, we first divide the rows of $\mathbf{A}$ into
$\mathcal{B}=\Theta(\log n)$ submatrices based on which of the $\mathcal{B}$
buckets $[1/2,1],[1/4,1/2],[1/8,1/4]$, etc. their leverage scores (computed
with respect to $\mathbf{A}$) fall into. Since each of these $\mathcal{B}$
submatrices comprises rows with leverage scores in the range $[a,2a]$ for some
$a$, per Fact 2.2, we have
$\sqrt{\frac{a}{n}}\leq\boldsymbol{\sigma}_{1}(\mathbf{a}_{i})\leq\sqrt{2a}$.
Therefore, by Lemma A.1, the total sensitivity of this submatrix may be
approximated by the appropriately scaled total sensitivity of $O(\sqrt{n})$ 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 [36]). Computing the total sensitivity of an
$\widetilde{O}(\sqrt{n})$-sized matrix thus costs
$\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+\sqrt{n}\cdot d^{\omega})$
under this scheme.

Applying the idea sketched above recursively further reduces the cost.
Suppose, after dividing the input matrix $\mathbf{A}$ into $\mathcal{B}$
submatrices $\mathbf{M}_{i}$ and subsampling $\widetilde{\mathbf{M}}_{i}$ from
these matrices, we try to recurse on each of these $\mathcal{B}$ submatrices
$\widetilde{\mathbf{M}}_{i}$. A key requirement for this idea to work is that
the sum of the $\ell_{1}$ sensitivities of the rows of
$\widetilde{\mathbf{M}}_{i}$ in isolation be close enough to the true sum of
$\ell_{1}$ sensitivities of those rows with respect to the original matrix
$\mathbf{A}$. In general, this is not true. However, we can meet our
requirement using the $\ell_{1}$ subspace embedding $\mathbf{S}\mathbf{A}$ of
the (original) matrix $\mathbf{A}$. Specifically, if we vertically concatenate
a matrix $\widetilde{\mathbf{M}}_{i}$ with $\mathbf{S}\mathbf{A}$ and call
this matrix $\mathbf{C}$, then for each row
$\widetilde{\mathbf{M}}_{i}[j]\in\widetilde{\mathbf{M}}_{i}$, we have
$\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\widetilde{\mathbf{M}}_{i}[j])\approx\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\widetilde{\mathbf{M}}_{i}[j])$.
Further, $|\mathbf{C}|\leq O(\sqrt{n})$, and so in the step subsampling the
rows of $\widetilde{\mathbf{M}}_{i}$, we choose $O(n^{1/4})$ samples.
Recursing this exponentially decreases the row dimension of the input matrix
at each level, with $(\log n)^{k}$ matrices at the $k^{\mathrm{th}}$ level of
recursion.

####  B.1.2 Helper lemmas for proving Theorem B.1

For our formal proof, we need Definition B.2, Fact B.3, and Definition B.4.

######  Definition B.2 (Notation for recursion tree of Algorithm 5).

Recall that $\mathcal{B}=\Theta(\log(n))$ as in 3 of Algorithm 5. Then we have
the following notation corresponding to Algorithm 5.

  * We define a rooted $\mathcal{B}$-ary tree $\mathcal{T}$ corresponding to the execution of Algorithm 5.

  * We use $\mathcal{T}_{j}$ to denote the subtree starting from the root node with all the nodes up to and including those at the $j^{\mathrm{th}}$ level. We denote the $i^{\mathrm{th}}$ node at the $j^{\mathrm{th}}$ level by $\mathcal{T}_{(j,i)}$. Thus, every node is uniquely specified by its level in the tree and its index within that level.

  * The root node $\mathcal{T}_{(1,1)}$ at level $1$ corresponds to the first call to Algorithm 5 from Algorithm 4.

  * For the node $\mathcal{T}_{(j,i)}$, we denote by $\mathbf{M}^{(j,i)}$ 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, we use $\widehat{s}_{(j,i)}$ to denote our estimate of $\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(j,i)})$, the sum of $\ell_{1}$ sensitivities of rows of the matrix $\mathbf{M}^{(j,i)}$ in the node $\mathcal{T}_{(j,i)}$, defined with respect to the sketched matrix $\mathbf{S}^{\prime}\mathbf{A}$.

######  Fact B.3.

If event $A$ happens with probability $1-\tau_{1}$ and event $B$ happens with
probability $1-\tau_{2}$ conditioned on event $A$, then the probability of
both events $A$ and $B$ happening is at least $1-\tau_{1}-\tau_{2}$.

######  Definition B.4.

We say that a node $\mathcal{N}$ with matrix $\mathbf{M}$ satisfies the
$(\widehat{\rho},\widehat{\delta})$-approximation guarantee if $\widehat{s}$,
the output of Algorithm 5 on $\mathbf{M}$, satisfies the following guarantee
for the true sum
$\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})$ of $\ell_{1}$
sensitivities of the rows in $\mathbf{M}$:

$\widehat{s}\approx_{\widehat{\rho}}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M})\text{
with a probability at least }1-\widehat{\delta}$

for some approximation parameter $\widehat{\rho}\in(0,1)$ and an error
probability parameter $\widehat{\delta}\in(0,1)$.

######  Lemma B.5 (Partitioning of $\mathbf{M}$).

In Algorithm 5, even though we ignore the range
$\left[0,\frac{1}{n^{20}}\right]$ when creating the submatrices
$\mathbf{M}_{1},\mathbf{M}_{2},\dotsc,\mathbf{M}_{\mathcal{B}}$ from
$\mathbf{M}$, _every_ row in $\mathbf{M}$ is part of exactly one of the
submatrices $\mathbf{M}_{1},\mathbf{M}_{2},\dotsc,\mathbf{M}_{\mathcal{B}}$.
In other words, the matrix $\mathbf{M}$ is partitioned into the matrices
$\mathbf{M}_{1},\mathbf{M}_{2},\dotsc,\mathbf{M}_{\mathcal{B}}$, with no row
missed out.

###### Proof.

Our proof strategy is to show that for all row indices $k\in[|\mathbf{M}|]$,
each of the leverage scores $\tau_{k}(\mathbf{C})$ satisfy
$\tau_{k}(\mathbf{C})\geq\frac{1}{n^{20}}$. Therefore, creating buckets of
rows with the range of leverage score starting at $\frac{1}{n^{20}}$ (instead
of $0$) does not miss out any rows.

To this end, we note that for every row $\mathbf{a}_{k}$ in the input
${\mathbf{M}}$ to Algorithm 5, we have:

$\sqrt{\frac{1}{n^{11}}}\leq\sqrt{\frac{\tau_{k}(\mathbf{A})}{|\mathbf{A}|}}\leq\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{k})=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{a}_{k}|}{\|\mathbf{A}\mathbf{x}\|_{1}}\leq\max_{\mathbf{x}\in\mathbb{R}^{d}}3\frac{|\mathbf{x}^{\top}\mathbf{a}_{k}|}{\|\mathbf{C}\mathbf{x}\|_{1}}=3\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{a}_{k})\leq
3\sqrt{\tau_{k}(\mathbf{C})},$

where the first step is because in 5 of Algorithm 4, we discard those rows
$\mathbf{a}_{k}$ of $\mathbf{A}$ with leverage scores smaller than $1/n^{10}$
and there are $n$ rows in $\mathbf{A}$, so every row $\mathbf{a}_{k}$ in the
input matrix ${\mathbf{M}}$ to Algorithm 5 satisfies
$\tau_{k}(\mathbf{A})\geq\frac{1}{n^{10}}$; the second step is by Fact 2.2
applied to $\mathbf{a}_{k}$; the third step is by the definition of the
$\ell_{1}$ sensitivity of $\mathbf{a}_{k}$ with respect to $\mathbf{A}$; the
fourth step is by observing that
$\|\mathbf{C}\mathbf{x}\|_{1}=\|\mathbf{S}\mathbf{A}\mathbf{x}\|_{1}+\|\mathbf{M}\mathbf{x}\|_{1}\leq
2\|\mathbf{A}\mathbf{x}\|_{1}+\|\mathbf{A}\mathbf{x}\|_{1}=3\|\mathbf{A}\mathbf{x}\|_{1}$
since $\mathbf{S}\mathbf{A}$ is a constant-factor (with the constant assumed
$1/2$) $\ell_{1}$ subspace embedding of $\mathbf{A}$ and because $\mathbf{M}$
is a submatrix of $\mathbf{A}$; the fifth step is by the definition of
$\ell_{1}$ sensitivity of $\mathbf{a}_{k}$ with respect to $\mathbf{C}$; the
final step is by Fact 2.2 applied to $\mathbf{a}_{k}$ within $\mathbf{C}$.
Observing the first and last terms of this inequality chain, we have
$\tau_{k}(\mathbf{C})\geq\frac{1}{3n^{11}}$ for all $k\in[|\mathbf{C}|]$. This
finishes the proof of the claim. ∎

######  Lemma B.6.

For the recursion tree $\mathcal{T}$ corresponding to Algorithm 5, let the
size $b$ of the input matrices at the leaf nodes be
$b=\tfrac{50D^{3}\log(100\mathcal{B})}{\gamma^{2}}\cdot\left(\tfrac{2D^{3}\log(100\mathcal{B})}{\gamma^{2}}+\sqrt{d\log(d)}\right)$.
Then, the depth of the recursion tree is less than
$D:=1+\log(\log(2n+2d\log(d)))$.

###### Proof.

For now, let $\beta$ be arbitrary but less than $\beta^{2}\leq n+d^{\prime}$,
where $d^{\prime}$ is the row dimension of the subspace embedding
$\mathbf{S}\mathbf{A}$. Then, the size of the input matrix to the nodes of the
recursion tree as we increase in depth evolves as

$\beta\sqrt{d^{\prime}+n},\,\beta\sqrt{d^{\prime}+\beta\sqrt{d^{\prime}+n}},\,\beta\sqrt{d^{\prime}+\beta\sqrt{d^{\prime}+\beta\sqrt{d^{\prime}+n}}}.\dotsc$

Thus, this sequence evolves as

$x_{t+1}=\beta\sqrt{d^{\prime}+x_{t}},\text{ with }x_{0}=n+d^{\prime}.$

To see the limit of this sequence, we set $x=\beta\sqrt{d^{\prime}+x},$ and
solve for this quadratic equation as

$x=\frac{\beta^{2}+\sqrt{\beta^{4}+4\beta^{2}d^{\prime}}}{2}\leq\beta(\beta+\sqrt{d^{\prime}}).$

We set the base case to be $b^{\star}:=7\beta(\beta+\sqrt{d^{\prime}})$, a
constant times larger than this asymptotic limit. We now compute the depth
$t^{\star}$ at which the recursion attains the value $b^{\star}$ is attained.
To do so, we consider a sequence $\\{y_{t}\\}$ which overestimates the
sequence $\\{x_{t}\\}$ of matrix sizes for a big range of $t$. We define
$\\{y_{t}\\}$ as

$y_{0}=n+d^{\prime},\text{ and }y_{t+1}=\beta\sqrt{2y_{t}},$

and it can be checked that for $t$ such that $y_{t}\geq d^{\prime}$, we have
$y_{t+1}\geq x_{t+1}$. Moreover, for all $t\geq 1$,
$y_{t}=\beta^{1+1/2+\dots+1/2^{t-1}}(2y_{0})^{1/2^{t}}=\beta^{2-2^{1-t}}(2y_{0})^{1/2^{t}}$.
The limit of this sequence can easily be seen to $\beta^{2}$. We now split our
analysis in two cases based on the limits of these two sequences.

##### Case 1: $\beta^{2}\geq d^{\prime}$.

Since the limit of the sequence $\\{y_{t}\\}$ is $\beta^{2}$, it is easy to
see that $y_{t}$ is always greater $d^{\prime}$ (since $y_{t}$ is
monotonically decreasing as $y_{0}\geq\beta^{2}$). Thus, $x_{t}$ will always
be less than $y_{t}$ in this case. We will now show that $y_{t}$ (and thus
$x_{t}$) will be less than $b^{\star}$ in $t_{0}=2\log\log(2y_{0})$ steps.
Solving the inequality for $t$ as follows:
$y_{t}\leq\beta^{2}(2y_{0})^{1/2^{t}}\leq 4\beta^{2}\leq b^{\star}$, it
suffices to ensure that $(2y_{0})^{1/2^{t}}\leq 4$, which happens if $t\geq
t_{0}$. Thus, $t^{\star}\leq 2\log\log 2y_{0}$.

##### Case $2$: $\beta^{2}<d^{\prime}$.

In this case, the sequence $y_{t}$ will eventually go below $d^{\prime}$ and
thus the inequality $y_{t}\geq x_{t}$ might eventually break down. Let
$t=t^{\prime}$ denote the first time step at which $y_{t}\leq 5d^{\prime}$
(which will be finite in this case). Then, by definition, $y_{t}\geq x_{t}$
for all $t<t^{\prime}$. We will now upper bound the value of $t^{\prime}$ as
follows: $y_{t}\leq\beta^{2}(2y_{0})^{1/2^{t}}\leq
d^{\prime}(2y_{0})^{1/2^{t}}\leq 4d^{\prime}$ ; Thus, it suffices to ensure
that $(2y_{0})^{1/2^{t}}\leq 4$, which happens if $t\geq
t_{0}:=2\log\log(2y_{0})$. Hence, $t^{\prime}\leq t_{0}$. Moreover, the value
of $y_{t^{\prime}}\leq 5d^{\prime},$ implies that $x_{t^{\prime}-1}\leq
y_{t^{\prime}-1}\leq\tfrac{({5}d^{\prime})^{2}}{2\beta^{2}}$. We will now show
that $x_{t^{\prime}+1}\leq b^{\star}$ using the evolution of
$x_{t+1}=\beta\sqrt{d^{\prime}+x_{t}}$ as follows: First,
$x_{t^{\prime}}\leq\beta\sqrt{d^{\prime}}+\beta\sqrt{x_{t^{\prime}-1}}\leq\beta\sqrt{d^{\prime}}+5d^{\prime}\leq
6d^{\prime}$; Then,
$x_{t^{\prime}+1}\leq\beta\sqrt{d^{\prime}}+\beta\sqrt{6d^{\prime}}\leq
7\beta\sqrt{d^{\prime}}\leq b^{\star}$. Thus, $t^{\star}\leq 1+2\log\log
2y_{0}$. Therefore, in both of these cases, we have that after at most
$D=1+\log\log(2n+2d^{\prime})$ depth, the base case of
$b^{\star}=7\beta(\beta+\sqrt{d^{\prime}})$ is achieved at the leaf node.
Recall that the algorithm sets $D=2\log\log(n+d^{\prime})+1$ and
$\rho=\gamma/D$ for $d^{\prime}=d\log(d)$ which implies that
$\beta:=C\frac{(1+\rho)}{\rho^{2}}\log(\delta^{-1})=CD\frac{(D+\gamma)}{\gamma^{2}}\log(100\mathcal{B}^{D})\leq\frac{2D^{3}\log(100\mathcal{B})}{\gamma^{2}}$.
Since Algorithm 5 sets
$b=\tfrac{50D^{3}\log(100\mathcal{B})}{\gamma^{2}}\cdot\left(\tfrac{2D^{3}\log(100\mathcal{B})}{\gamma^{2}}+\sqrt{d\log(d)}\right)$,
which is larger than $b^{\star}$, we may conclude that the choice of $b$
implies that the recursion tree $\mathcal{T}$ of Algorithm 5 has a depth less
than $D=1+\log(\log(2n+2d\log(d)))$.

∎

######  Lemma B.7.

Suppose the root node $\mathcal{T}_{(1,1)}$ of Algorithm 5 satisfies a
$(\widehat{\rho},\widehat{\delta})$-approximation guarantee as defined in
Definition B.4, and denote by $\widehat{s}$ the output of Algorithm 4. Then,
with a probability of at least $1-\widehat{\delta}$, we have

$\frac{1}{1+\gamma}\widehat{s}\approx_{\widehat{\rho}+\rho}\mathfrak{S}_{1}^{\mathbf{A}}(\mathbf{A}).$

As a result of this lemma, in the final proof of correctness of Algorithm 4
(appearing later as proof of Theorem B.1), it suffices to show that
$\widehat{s}$ satisfies a $(\widehat{\rho},\widehat{\delta})$-approximation
guarantee for some appropriate choices of these parameters.

###### Proof of Lemma B.7.

Suppose the root node satisfies Definition B.4 with some
$(\widehat{\rho},\widehat{\delta})$-approximation guarantee. This means that
there exists an $\widetilde{s}$ such that the output of Algorithm 5 on the
input matrix $\widehat{\mathbf{A}}$ satisfies

$\widetilde{s}\approx_{\widehat{\rho}}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\widehat{\mathbf{A}})\text{
with a probability at least }1-\widehat{\delta}.$

Since $\mathbf{S}^{\prime}\mathbf{A}$ is, by design, a $\rho$-approximation
$\ell_{1}$ subspace embedding for $\mathbf{A}$, it means for all
$\mathbf{x}\in\mathbb{R}^{d}$ we have

$(1-\rho)\|\mathbf{A}\mathbf{x}\|_{1}\leq\|\mathbf{S}^{\prime}\mathbf{A}\mathbf{x}\|_{1}\leq(1+\rho)\|\mathbf{A}\mathbf{x}\|_{1}.$

Therefore, for every row $\mathbf{a}_{i}\in\widehat{\mathbf{A}}$, we have by
the definition of $\ell_{1}$ sensitivity and Section B.1.2 that the $\ell_{1}$
sensitivity of $\mathbf{a}_{i}$ computed with respect to
$\mathbf{S}^{\prime}\mathbf{A}$ is a $\rho$-approximation of the $\ell_{1}$
sensitivity of $\mathbf{a}_{i}$ with respect to $\mathbf{A}$:

$\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{a}_{i})=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\|\mathbf{S}^{\prime}\mathbf{A}\mathbf{x}\|_{1}}\approx_{\rho}\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|}{\|\mathbf{A}\mathbf{x}\|_{1}}=\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{i}).$

Therefore, by linearity, this $\rho$-approximation factor transfers over in
connecting the total sensitivity of $\widehat{\mathbf{A}}$ with respect to
$\mathbf{S}^{\prime}\mathbf{A}$ to the total sensitivity of
$\widehat{\mathbf{A}}$ with respect to $\mathbf{A}$:

$\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\widehat{\mathbf{A}})=\sum_{i\in[|\widehat{\mathbf{A}}|]}\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{a}_{i})\approx_{\rho}\sum_{i\in[|\widehat{\mathbf{A}}|]}\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{i})=\mathfrak{S}_{1}^{\mathbf{A}}(\widehat{\mathbf{A}}).$

Therefore, chaining Section B.1.2 and Section B.1.2 yields the following
approximation guarantee:

$\widetilde{s}\approx_{\widehat{\rho}+\rho}\mathfrak{S}_{1}^{\mathbf{A}}(\widehat{\mathbf{A}})\text{
with a probability at least }1-\widehat{\delta}.$

There is no change in the error probability above from Section B.1.2 because
Section B.1.2 holds deterministically. Next, since in 5 of Algorithm 4, we
dropped rows of $\mathbf{A}$ that have leverage scores less than
$\frac{1}{n^{10}}$, for each of the dropped rows $\mathbf{a}_{i}$, we have

$\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{i})\leq\sqrt{\tau_{i}(\mathbf{A})}\leq\frac{1}{n^{5}}.$

Therefore, the quantity we return, $\widehat{s}$, satisfies the following
approximation guarantee:

$\frac{1}{1+\gamma}\widehat{s}:=\widetilde{s}+\frac{1}{n^{5}}(|\mathbf{A}|-|\widehat{\mathbf{A}}|)\geq(1-(\widehat{\rho}+\rho))\mathfrak{S}_{1}^{\mathbf{A}}(\widehat{\mathbf{A}})+\sum_{i:\mathbf{a}_{i}\notin\widehat{\mathbf{A}}}\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{a}_{i})\geq(1-(\widehat{\rho}+\rho))\mathfrak{S}_{1}^{\mathbf{A}}(\mathbf{A}),$

where the first step is by definition of $\widehat{s}$ in Algorithm 4; the
second step is by Section B.1.2; the third step is by applying Section B.1.2
to each of the rows not present in $\widehat{\mathbf{A}}$. In the other
direction,

$\frac{1}{1+\gamma}\widehat{s}=\widetilde{s}+\frac{1}{n^{5}}(|\mathbf{A}|-|\widehat{\mathbf{A}}|)\leq(1+(\widehat{\rho}+\rho))\mathfrak{S}_{1}^{\mathbf{A}}(\widehat{\mathbf{A}})+\frac{1}{n^{4}}\leq(1+(\widehat{\rho}+\rho))\mathfrak{S}_{1}^{\mathbf{A}}(\mathbf{A}),$

where the first step is by definition of $\widehat{s}$ in Algorithm 4; the
second step is by the fact that $|\mathbf{A}|\leq n$; the final step is by
combining the assumption $\rho\geq\frac{1}{n^{4}}$ and the facts that all
sensitivities are non-negative and that the total $\ell_{1}$ sensitivity is at
least one. Combining Section B.1.2 and Section B.1.2 finishes the proof. ∎

######  Lemma B.8 (Single-Level Computation in Recursion Tree).

Consider a node $\mathcal{N}$ in the recursion tree corresponding to Algorithm
5, and let $\mathbf{M}$ be the input matrix for $\mathcal{N}$. Further assume
that $\mathcal{N}$ is not a leaf node. Denote by
$\mathcal{N}_{1},\mathcal{N}_{2},\dotsc,\mathcal{N}_{\mathcal{B}}$ (where
$\mathcal{B}=\Theta(\log(n))$ is the number of recursive calls at each node in
Algorithm 5) the children of node $\mathcal{N}$, with each node
$\mathcal{N}_{i}$ with the input matrix $\widetilde{\mathbf{M}}_{i}$ (as
constructed in Algorithm 5). Suppose each of these child nodes
$\mathcal{N}_{i}$ satisfies a
$(\widehat{\rho},\widehat{\delta})$-approximation guarantee (cf. Definition
B.4). Then $\mathcal{N}$ satisfies a
$(\widehat{\rho}+\rho,(\widehat{\delta}+\delta)\cdot\mathcal{B})$-approximation
guarantee.

###### Proof.

We first reiterate the steps of Algorithm 5 that are crucial to this proof. As
proved in Lemma B.5, Algorithm 5 partitions the rows of the matrix
$\mathbf{M}$ in node $\mathcal{N}$ into $\mathcal{B}$ submatrices
$\mathbf{M}_{1},\mathbf{M}_{2},\dotsc,\mathbf{M}_{\mathcal{B}}$ based on their
leverage scores. Next, for each matrix $\mathbf{M}_{i}$, we uniformly sample
with replacement $O(\sqrt{|\mathbf{C}|}(1+\rho)\rho^{-2}\log(\delta^{-1}))$
rows and term this set of rows as matrix $\widetilde{\mathbf{M}}_{i}$, which
we then pass to Algorithm 5 recursively in node $\mathcal{N}_{i}$ (recall that
we assume that the node $\mathcal{N}$ is not a leaf node). In other words,
each child node $\mathcal{N}_{i}$ has as its input matrix
$\widetilde{\mathbf{M}}_{i}$.

Since we assume that all the child nodes
$\mathcal{N}_{1},\mathcal{N}_{2},\dotsc,\mathcal{N}_{\mathcal{B}}$ satisfy a
$(\widehat{\rho},\widehat{\delta})$-approximation guarantee as in Definition
B.4, it means that, for all $i\in[\mathcal{B}]$, the output $\widehat{s}_{i}$
of Algorithm 5 on node $\mathcal{N}_{i}$ satisfies the following guarantee:

$\widehat{s}_{i}\approx_{\widehat{\rho}}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}({\widetilde{\mathbf{M}}}_{i})\text{
with a probability at least }1-\widehat{\delta}.$

Next, we recall that each ${\mathbf{M}}_{i}$ is composed of only those rows of
$\mathbf{M}$ whose leverage scores (computed with respect to $\mathbf{C}$) lie
in the $i^{\mathrm{th}}$ bucket $[2^{-i},2^{-i+1}]$. Combining this with Fact
2.2, this implies that the $\ell_{1}$ sensitivity of the $j^{\mathrm{th}}$ row
$\mathbf{M}_{i}[j]$ of the matrix $\mathbf{M}_{i}$, computed with respect to
$\mathbf{C}$, satisfies

$\sqrt{\frac{2^{-i}}{|\mathbf{C}|}}\leq\sqrt{\frac{\tau_{j}^{\mathbf{C}}(\mathbf{M}_{i})}{|\mathbf{C}|}}\leq\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{M}_{i}[j])\leq\sqrt{\tau_{j}^{\mathbf{C}}(\mathbf{M}_{i})}\leq\sqrt{2^{-i+1}}.$

Additionally, we claim that

$\frac{1}{3(1+\rho)}\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{M}_{i}[j])\leq\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}_{i}[j])\leq
3(1+\rho)\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{M}_{i}[j]).$

To see this, we observe that

$\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{M}_{i}[j])=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{\|\mathbf{C}\mathbf{x}\|_{1}}\geq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{3\|\mathbf{A}\mathbf{x}\|_{1}}\geq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{3(1+\rho)\|\mathbf{S}^{\prime}\mathbf{A}\mathbf{x}\|_{1}}=\frac{\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}_{i}[j])}{3(1+\rho)},$

where the first step is by the definition of $\ell_{1}$ sensitivity of
$\mathbf{M}_{i}[j]$ (i.e., the $j^{\mathrm{th}}$ row of matrix
$\mathbf{M}_{i}$) with respect to $\mathbf{C}$; the second step is by the
definition of $\mathbf{C}$, the fact that $\mathbf{S}\mathbf{A}$ is a
$1/2$-factor $\ell_{1}$ subspace embedding of $\mathbf{A}$, and because
$\mathbf{M}$ is a subset of rows of $\mathbf{A}$; the third step is because
$\mathbf{S}^{\prime}\mathbf{A}$ is a $\rho$-approximate $\ell_{1}$ subspace
embedding for $\mathbf{A}$, and the final step is by the definition of
$\boldsymbol{\sigma}_{1}^{\mathbf{A}}(\mathbf{M}_{i}[j])$. In the other
direction, we have

$\boldsymbol{\sigma}_{1}^{\mathbf{C}}(\mathbf{M}_{i}[j])=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{\|\mathbf{C}\mathbf{x}\|_{1}}\leq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{0.5\|\mathbf{A}\mathbf{x}\|_{1}}\leq\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{(1+\rho)|\mathbf{x}^{\top}\mathbf{M}_{i}[j]|}{0.5\|\mathbf{S}^{\prime}\mathbf{A}\mathbf{x}\|_{1}}\leq
3(1+\rho)\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}_{i}[j]),$

where the second step is by the definition of $\mathbf{C}$ as a vertical
concatenation of $\mathbf{S}\mathbf{A}$ and $\mathbf{M}$ and further using
that $\mathbf{S}\mathbf{A}$ is a constant factor $\ell_{1}$ subspace embedding
of $\mathbf{A}$ and dropping the non-negative term
$\|\mathbf{M}\mathbf{x}\|_{1}$. Combining Section B.1.2 and Section B.1.2, we
see that for all rows $j\in[|\mathbf{M}_{i}|]$, we have

$\frac{1}{3(1+\rho)}\sqrt{\frac{2^{-i}}{|\mathbf{C}|}}\leq\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}_{i}[j])\leq
3(1+\rho)\sqrt{2^{-i+1}}.$

Since Section B.1.2 shows upper and lower bounds on each sensitivity
$\boldsymbol{\sigma}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}_{i})$ of
the matrix $\mathbf{M}_{i}$, we may use Lemma A.1 to approximate the sum of
these sensivities. In particular, by Lemma A.1, we may construct a matrix
$\widetilde{\mathbf{M}}_{i}$ composed of
$\frac{20\sqrt{|C|}(1+2\rho)\log(\delta^{-1})}{\rho^{2}}$ rows of
$\mathbf{M}_{i}$ sampled uniformly at random with replacement, and we are
guaranteed

$(1+\rho)\frac{|\mathbf{M}_{i}|}{|\widetilde{\mathbf{M}}_{i}|}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}({\widetilde{\mathbf{M}}}_{i})\approx_{{\rho}}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}({{\mathbf{M}}}_{i})\text{
with a probability at least }1-{\delta}.$

Note that the output $\widehat{s}$ of Algorithm 5 when applied to node
$\mathcal{N}$ is defined as

$\widehat{s}:=(1+\rho)\sum_{i\in[\mathcal{B}]}\frac{|\mathbf{M}_{i}|}{|\widetilde{\mathbf{M}}_{i}|}\widehat{s}_{i}.$

Therefore, by union bound over the failure probability (cf. Fact B.3), we can
combine Section B.1.2 and sum Section B.1.2 over all $i\in\mathcal{B}$ child
nodes to conclude that $\widehat{s}$ from Section B.1.2 satisfies a
$(\widehat{\rho}+\rho,(\widehat{\delta}+\delta)\mathcal{B})$-approximation
guarantee of Definition B.4. ∎

####  B.1.3 Proof of Theorem B.1

###### Proof of Theorem B.1.

We first sketch our strategy to prove Theorem B.1’s correctness guarantee.
Essentially, our goal is to establish Definition B.4 for the root node
$\mathcal{T}_{(1,1)}$. That this goal suffices to prove the theorem is
justified in Lemma B.7. We now show that the node $\mathcal{T}_{(1,1)}$
satisfies a $(\widehat{\rho},\widehat{\delta})$-approximation guarantee for
some $\widehat{\rho}$ and $\widehat{\delta}$. We use the method of induction.

##### Correctness guarantee of Theorem B.1.

We know that the leaf nodes all satisfy the $(\rho,0)$-approximation guarantee
in Definition B.4 since at this level, the base case is triggered, where we
compute a $\rho$-approximation estimate of the total $\ell_{1}$ sensitivity
with respect to the subspace embedding $\mathbf{S}^{\prime}\mathbf{A}$. Having
established this approximation guarantee at the leaf nodes, we can use Lemma
B.8 to inductively propagate the property of $(\rho,\delta)$-approximation
guarantee up the recursion tree $\mathcal{T}$. From Lemma B.6, we have that
the recursion tree $\mathcal{T}$ of Algorithm 5 is of depth at most
$D=O(\log\log(n+d))$, and the number of child nodes at any given node is at
most $\mathcal{B}=\Theta(\log(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 $\mathbf{S}^{\prime}\mathbf{A}$. Recall that we
denote the depth of the tree by $D$. Then, because
$\mathbf{S}^{\prime}\mathbf{A}$ is, by definition, a $\rho$-approximate
subspace embedding for $\mathbf{A}$, we have, for each $i$ that indexes a leaf
node,

$\widehat{s}_{(D,i)}\approx_{\rho}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(D,i)})\text{
with a probability of }1.$

##### Induction proof.

Consider the subtree $\mathcal{T}_{j}$ truncated at (but including) the nodes
at level $j$ in the recursion tree $\mathcal{T}$ corresponding to Algorithm 5.
Assume the induction hypothesis that all the leaf nodes in the subtree
$\mathcal{T}_{j}$ satisfy Definition B.4 with parameters
$\left(\rho\cdot(D+1-j),\delta\cdot\left(\frac{\mathcal{B}^{D-j+1}-\mathcal{B}}{\mathcal{B}-1}\right)\right)$.
That is, at the matrix of each such node $i$, we have an estimate
$\widehat{s}_{(j,i)}$ of its total sensitivity
$\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(j,i)})$ such
that

$\widehat{s}_{(j,i)}\approx_{\rho\cdot(D+1-j)}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(j,i)})\text{
with a probability at least
}1-\delta\cdot\left(\frac{\mathcal{B}^{D-j+1}-\mathcal{B}}{\mathcal{B}-1}\right).$

Then, each of the leaf nodes of the subtree $\mathcal{T}_{j-1}$ (i.e., the
subtree that stops one level above $\mathcal{T}_{j}$) satisfies the assumption
in Lemma B.8 since its child nodes are either leaf nodes of the entire
recursion tree $\mathcal{T}$ (for which this statement has been shown in
Section B.1.3) or, if they are not leaf nodes of the recursion tree, they
satisfy Definition B.4 with parameters $\widehat{\rho}=\rho\cdot(D+1-j)$ and
$\widehat{\delta}=\delta\cdot\left(\frac{\mathcal{B}^{D-j+1}-\mathcal{B}}{\mathcal{B}-1}\right)$.
Therefore, Lemma B.8 implies that each of the leaf nodes of
$\mathcal{T}_{j-1}$ satisfies Definition B.4 with approximation parameter
$\widehat{\rho}$ and error probability $\widehat{\delta}$ defined as follows:

$\widehat{\rho}=\rho+\rho\cdot(D+1-j)=\rho\cdot(D+2-j)\text{ and
}\widehat{\delta}=\mathcal{B}\cdot\left(\delta+\delta\cdot\left(\frac{\mathcal{B}^{D-j+1}-\mathcal{B}}{\mathcal{B}-1}\right)\right)=\delta\cdot\left(\frac{\mathcal{B}^{D-j+2}-\mathcal{B}}{\mathcal{B}-1}\right).$

In other words, the outputs $\widehat{s}_{(j-1,k)}$ of Algorithm 5 on the
matrices $\mathbf{M}^{(j-1,k)}$ in the leaf nodes of $\mathcal{T}_{j-1}$ each
satisfy:

$\widehat{s}_{(j-1,k)}\approx_{O\left(\rho\cdot\left(D+2-j\right)\right)}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(j-1,k)})\text{
with a probability at least
}1-\delta\cdot\left(\frac{\mathcal{B}^{D-j+2}-\mathcal{B}}{\mathcal{B}-1}\right).$

The base case in Section B.1.3 and the satisfaction of the induction
hypothesis in Section B.1.3 together finish the proof of the induction
hypothesis.

##### Finishing the proof of correctness.

In light of the above conclusion, for the root node $\mathcal{T}_{(1,1)}$,
which is at the level $j=1$, we may plug in $j=1$ in Section B.1.3 and obtain
that the output $\widehat{s}_{(1,1)}$ of Algorithm 5 at the root node
satisfies

$\widehat{s}_{(1,1)}\approx_{O(D\rho)}\mathfrak{S}_{1}^{\mathbf{S}^{\prime}\mathbf{A}}(\mathbf{M}^{(1,1)})\text{
with a probability of at least
}1-\delta\left(\frac{\mathcal{B}^{D}-\mathcal{B}}{\mathcal{B}-1}\right).$

In Algorithm 5, we choose the parameters $\delta=\frac{0.01}{\mathcal{B}^{D}}$
and $\rho=O\left(\frac{\gamma}{D}\right)$. Plugging these parameters into
Section B.1.3, we conclude that the root node satisfies Definition B.4 with a
$(\gamma,0.01)$-approximation guarantee. This finishes the proof of
correctness of the claim.

##### Proof of runtime.

Let $b$ denote the maximum number of rows of a matrix at a leaf node of
$\mathcal{T}$. Further, recall from Fact 2.4 that the cost of computing one
$\ell_{1}$ sensitivity of an $n\times d$ matrix is,

$L=\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+d^{\omega}).$

Then as stated in 7, the algorithm computes the $\rho$-approximate total
sensitivity at the leaf nodes; Fact 2.4, this incurs a cost of

$\mathcal{C}(b)=\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+b\cdot
d^{\omega}).$

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 $\mathcal{C}(r)$
denote the runtime of Algorithm 5 when the input is an $r\times d$ matrix.
Then,

$\mathcal{C}(r)=\left\\{\begin{array}[]{ll}\mathcal{B}\cdot\mathcal{C}\left(\frac{\sqrt{r}(1+\rho)}{\rho^{2}}\log(\delta^{-1})\right)+L&\mbox{
}r>b\\\ \mathcal{C}(b)&\mbox{ }r\leq b\end{array}\right.,$

where $b$ is the number of rows in the base case, and $L$ is the cost of
computing leverage scores and equals
$L=\widetilde{O}(\mathop{\mathbf{nnz}}(\mathbf{A})+d^{\omega})$ from Fact 2.1.
For the depth of the recursion $D$ (where the top level is $D=1$), plugging in
our choices of $\delta=\frac{0.01}{\mathcal{B}^{D}}$ and
$\rho=O\left(\frac{\gamma}{D}\right)$, this may be expressed as via expanding
the recursion as:

$\displaystyle\mathcal{C}(n)$
$\displaystyle=\mathcal{B}^{D}\cdot\mathcal{C}(b)+L\cdot\left(\frac{\mathcal{B}^{D}-1}{\mathcal{B}-1}\right)$
(B.24)

for $n\geq b$. We plug into this expression the values of $b$ and $D$ from
Lemma B.6, $L$ from Section B.1.3, $\mathcal{C}(b)$ from Section B.1.3, and
$\mathcal{B}=\Theta(\log(n))$ to get the claimed runtime. ∎

##  Appendix C Omitted proofs: $\ell_{p}$ sensitivities

We recall the following notation that we use for stating our results in this
setting. See Definition 1.2

###  C.1 Estimating all $\ell_{p}$ sensitivities

We generalize Algorithm 1 to the case $p\geq 1$ below.

Algorithm 6 Approximating $\ell_{p}$-Sensitivities: Row-wise Approximation

1:Inputs: Matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, approximation factor
$\alpha\in(1,n]$, scalar $p\geq 1$

2:Output: Vector $\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{>0}$ that
satisfies, for each $i\in[n]$, with probability $0.9$, that
$\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\alpha^{p-1}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})+\tfrac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A}))$

3:Compute, for $\mathbf{A}$, an $\ell_{p}$ subspace embedding
$\mathbf{S}_{p}\mathbf{A}\in\mathbb{R}^{\widetilde{O}(d^{\max(1,p/2)})\times
d}$

4:Partition $\mathbf{A}$ into $\tfrac{n}{\alpha}$ blocks
$\mathbf{B}_{1},\dotsc,\mathbf{B}_{n/\alpha}$ each comprising $\alpha$
randomly selected rows.

5:for the $\ell^{\mathrm{th}}$ block $\mathbf{B}_{\ell}$, with
$\ell\in[\tfrac{n}{\alpha}]$ do

6: Sample $\alpha$-dimensional independent Rademacher vectors
$\mathbf{r}_{1}^{(\ell)},\dotsc,\mathbf{r}_{100}^{(\ell)}$

7: For each $j{\in}[100]$, compute the row vectors
$\mathbf{r}_{j}^{(\ell)}\mathbf{B}_{\ell}\in\mathbb{R}^{d}$

8:end for

9:Let $\mathbf{P}\in\mathbb{R}^{100\tfrac{n}{\alpha}\times d}$ be the matrix
of all vectors computed in 7. Compute
$\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{P})$ using
Definition 1.2

10:for each $i\in[n]$ do

11: Denote by $J$ the set of row indices in $\mathbf{P}$ that $\mathbf{a}_{i}$
is mapped to in 7

12: Set $\widetilde{\boldsymbol{\sigma}}_{i}=\max_{j\in
J}(\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{p}_{j}))$

13:end for

14:Return $\widetilde{\boldsymbol{\sigma}}$

######  Theorem C.1 (Approximating all $\ell_{p}$ sensitivities).

Given a full-rank matrix $\mathbf{A}\in\mathbb{R}^{n\times d}$, an
approximation factor $1<\alpha\ll n$, and a scalar $p\geq 1$, let
$\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})$ be the $i^{\mathrm{th}}$ $\ell_{p}$
sensitivity. Then there exists an algorithm that returns a vector
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}_{\geq 0}$ such that with
high probability, for each $i\in[n]$, we have

$\Omega(\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}))\leq\widetilde{\boldsymbol{\sigma}}_{i}\leq
O(\alpha^{p-1}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}))+\tfrac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A})).$

Our algorithm runs in time
$\widetilde{O}\left(\mathop{\mathbf{nnz}}(\mathbf{A})+\tfrac{n}{\alpha}\cdot\mathsf{LP}(d^{\max(1,p/2)},d,p)\right)$
(cf. Definition 1.2).

###### Proof.

Suppose the $i^{\mathrm{th}}$ row of $\mathbf{A}$ falls into the bucket
$\mathbf{B}_{\ell}$. Suppose, further, that the rows from $\mathbf{B}_{\ell}$
are mapped to those in $\mathbf{P}$ with row indices in the set $J$. Then,
Algorithm 6 returns a vector of sensitivity estimates
$\widetilde{\boldsymbol{\sigma}}\in\mathbb{R}^{n}$ with the $i^{\mathrm{th}}$
coordinate defined as $\widetilde{\boldsymbol{\sigma}}_{i}=\max_{j\in
J}\boldsymbol{\sigma}_{p}(\mathbf{c}_{j})$. For any
$\mathbf{x}\in\mathbb{R}^{d}$, we have

$\|\mathbf{S}_{p}\mathbf{A}\mathbf{x}\|_{p}^{p}\approx\|\mathbf{A}\mathbf{x}\|_{p}^{p}.$

We use Section C.1 to establish the claimed bounds below. Let
$\mathbf{x}^{\star}=\arg\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}|^{p}}{\|\mathbf{A}\mathbf{x}\|_{p}^{p}}$
be the vector that realizes the $i^{\mathrm{th}}$ $\ell_{p}$ sensitivity of
$\mathbf{A}$. Further, let the $j^{\mathrm{th}}$ row of $\mathbf{P}$ be
$\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}$. Then, with probability of at least
$1/2$, we have

$\displaystyle\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{p}_{j})$
$\displaystyle=\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}|^{p}}{\|\mathbf{S}_{p}\mathbf{A}\mathbf{x}\|_{p}^{p}}\geq\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}^{\star}|^{p}}{\|\mathbf{S}_{p}\mathbf{A}\mathbf{x}^{\star}\|_{p}^{p}}=\Theta(1)\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}^{\star}|^{p}}{\|\mathbf{A}\mathbf{x}^{\star}\|_{p}^{p}}\geq\Theta(1)\frac{|\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}|^{p}}{\|\mathbf{A}\mathbf{x}^{\star}\|_{p}^{p}}=\Theta(1)\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}),$
(C.3)

where the first step is by definition of
$\boldsymbol{\sigma}_{p}(\mathbf{p}_{j})$, the second is by evaluating the
function being maximized at a specific choice of $\mathbf{x}$, the third step
is by Section C.1, and the final step is by definition of $\mathbf{x}^{\star}$
and $\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})$. For the fourth step, we use the
fact that
$|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}^{\star}|^{p}=|\mathbf{r}_{k,i}^{(\ell)}\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}+\sum_{j\neq
i}\mathbf{r}_{k,j}^{(\ell)}(\mathbf{B}\mathbf{x}^{\star})_{j}|^{p}\geq|\mathbf{a}_{i}^{\top}\mathbf{x}^{\star}|^{p}$
with a probability of at least $1/2$ since the vector
$\mathbf{r}_{k}^{(\ell)}$ has coordinates that are $-1$ or $+1$ with equal
probability. By a union bound over $j\in J$ independent rows that block
$\mathbf{B}_{\ell}$ is mapped to, we establish the claimed lower bound in
Equation C.3 with probability at least $0.9$. To show an upper bound on
$\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{p}_{j})$, we
observe that

$\displaystyle\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}|^{p}}{\|\mathbf{S}_{p}\mathbf{A}\mathbf{x}\|_{p}^{p}}\leq{\alpha^{p-1}}\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{\|\mathbf{B}_{\ell}\mathbf{x}\|_{p}^{p}}{\|\mathbf{S}_{p}\mathbf{A}\mathbf{x}\|_{p}^{p}}\leq\Theta(\alpha^{p-1})\sum_{j:\mathbf{a}_{j}\in\mathbf{B}_{\ell}}\max_{\mathbf{x}\in\mathbb{R}^{d}}\frac{|\mathbf{a}_{j}^{\top}\mathbf{x}|^{p}}{\|\mathbf{A}\mathbf{x}\|_{p}^{p}}=\Theta(\alpha^{p-1})\sum_{j:\mathbf{a}_{j}\in\mathbf{B}_{\ell}}\boldsymbol{\sigma}_{p}(\mathbf{a}_{j}),$
(C.4)

where the second step is by Section C.1 and opening up
$\|\mathbf{B}_{\ell}\mathbf{x}\|_{p}^{p}$ in terms of the rows in
$\mathbf{B}_{\ell}$, and the final step is by definition of
$\mathbf{B}_{\ell}$. To see the first step, we observe that

$|\mathbf{r}_{k}^{(\ell)}\mathbf{B}_{\ell}\mathbf{x}|^{p}\leq\|\mathbf{B}_{\ell}\mathbf{x}\|_{p}^{p}\cdot\|\mathbf{r}_{k}^{(\ell)}\|_{q}^{p}\leq\|\mathbf{B}_{\ell}\mathbf{x}\|_{p}^{p}\cdot\alpha^{p/q}=\|\mathbf{B}_{\ell}\mathbf{x}\|_{p}^{p}\cdot\alpha^{p-1},$

where $\frac{1}{p}+\frac{1}{q}=1$, the first step is by Hölder’s inequality,
the second is because each entry of $\mathbf{r}_{k}^{(\ell)}$ is either $+1$
or $-1$ and $\mathbf{r}_{k}^{(\ell)}\in\mathbb{R}^{\alpha}$, and the third is
because $\frac{p}{q}=p-1$. Because $\mathbf{B}_{\ell}$ is a group of $\alpha$
rows selected uniformly at random out of $n$ rows and contains the row
$\mathbf{a}_{i}$, we have:

$\mathbb{E}\left\\{\sum_{j:\mathbf{a}_{j}\in\mathbf{B}_{\ell},j\neq
i}\boldsymbol{\sigma}_{p}(\mathbf{a}_{j})\right\\}=\frac{\alpha-1}{n-1}\sum_{j\neq
i}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i}).$

Therefore, Markov inequality gives us that with a probability of at least
$0.9$, we have

$\boldsymbol{\sigma}_{p}^{\mathbf{S}_{p}\mathbf{A}}(\mathbf{p}_{j})\leq
O(\alpha^{p-1}\boldsymbol{\sigma}_{p}(\mathbf{a}_{i})+\tfrac{\alpha^{p}}{n}\mathfrak{S}_{p}(\mathbf{A})).$

The median trick then establishes Equation C.4 and Section C.1 with high
probability. To establish the runtime, note that we first construct the
subspace embedding $\mathbf{S}_{p}\mathbf{A}$ and then compute the $\ell_{p}$
sensitivities of $\mathbf{P}$ (with $O(n/\alpha)$ rows) with respect to
$\mathbf{S}_{p}\mathbf{A}$. The size of the subspace embedding
$\mathbf{S}_{p}\mathbf{A}$ is $O(d^{\max(1,p/2)}\times d)$ (see [76, Table
$1$]). This completes the runtime claim. ∎

###  C.2 Estimating the maximum of $\ell_{p}$ sensitivities

In Section 3.2, we showed our result for estimating the maximum $\ell_{1}$
sensitivity. In this section, we show how to extend this result to all $p\geq
1$. Algorithm 7 generalizes Algorithm 3.

Algorithm 7 Approximating the Maximum of $\ell_{p}$-Sensitivities