We study active sampling algorithms for linear regression, which aim to query only a small number of entries of a target vector $\bfb \in \mathbb{R}^n$ and output a near minimizer to $\min_{\bfx \in \mathbb{R}^d} \| \bfA\bfx - \bfb\|$, where $\bfA \in \mathbb{R}^{n \times d}$ is a design matrix and $\| \cdot \|$ is some loss function.

For $\ell_p$ norm regression for any $0<p<\infty$, we give an algorithm based on Lewis weight sampling which outputs a $(1+\epsilon)$-approximate solution using just $\tilde{O}(d/\eps^2)$ queries to $\bfb$ for $p\in(0,1)$, $\tilde{O}(d/\eps)$ queries for $p\in(1,2)$, and $\tilde{O}(d^{p/2}/\eps^{p})$ queries for $p\in(2,\infty)$. For $p\in(0,2)$, our bounds are optimal up to logarithmic factors, thus settling the query complexity for this range of $p$. For $p\in(2,\infty)$, our dependence on $d$ is optimal, while our dependence on $\eps$ is off by at most a single $\eps$ factor, up to logarithmic factors. Our result resolves an open question of Chen and Dereziński, who gave near optimal bounds for the $\ell_1$ norm, but required at least $d^2/\eps^2$ samples for $\ell_p$ regression with $p \in (1,2)$, and gave no bounds for $p\in(2,\infty)$ or $p\in(0,1)$.

We also provide the first total sensitivity upper bound of $O(d^{\max\{1,p/2\}}\log^2 n)$ for loss functions with at most degree $p$ polynomial growth. This improves a recent result of Tukan, Maalouf, and Feldman. By combining this with our techniques for $\ell_p$ regression, we obtain an active regression algorithm making $\tilde O(d^{1+\max\{1,p/2\}}/\poly(\eps))$ queries for such loss functions, including the important cases of the Tukey and Huber losses. This answers another question of Chen and Dereziński. For the Huber loss, we further improve our bound to a sample complexity of $\tilde O(d^{4-2\sqrt 2}/\poly(\eps))$ where $4-2\sqrt 2\approx 1.17157$.

Our sensitivity bounds also give improvements to a variety of previous results using sensitivity sampling, including Orlicz norm subspace embeddings, robust subspace approximation, and dimension reduction for smoothed $p$-norms.

Finally, our active sampling results give the first sublinear time algorithms for Kronecker product regression under every $\ell_p$ norm. Previous results required reading the entire $\bfb$ vector in the kernel feature space.

§ INTRODUCTION

We consider a classic active learning problem: given a design matrix $\bfA\in \R^{n\times d}$ and query access to entries of an unknown target (measurement) vector $\bfb \in \R^n$, how can we compute an approximate minimizer of the regression problem $\min_{\bfx\in\R^d} \|\bfA\bfx - \bfb\|$ while querying as few entries of $\bfb$ as possible? This problem arises in applications where labeled data is expensive: viewing a single entry of $\bfb$ might require running a survey, physical experiment, or time-intensive computer simulation [60, 55].
Concretely, we study the following problem for general vector norms[Our work will also extend to other loss functions of the form $\sum_{i=1}^n M([\bfA\bfx-\bfb]_i)$ that are not necessarily norms.] $\|\cdot \|$:

For $\bfA \in \R^{n\times d}, \bfb\in \R^n$, and accuracy parameter $0 < \epsilon \leq 1$, find $\tilde{\bfx}\in \R^d$ satisfying:

\begin{align*}
\|\bfA\tilde{\bfx} - \bfb\| \leq (1+\epsilon)\cdot \min_{\bfx\in\R^d} \|\bfA\bfx - \bfb\|,
\end{align*}

while reading as few of the entries $\{\bfb(1), \ldots, \bfb(n)\}$ of the target vector $\bfb$ as possible.[In principal, entries of $\bfb$ can be read adaptively – i.e., we can select indices to query based on the results of other queries. However, the benefits of adaptivity appear limited. Most methods for solving Problem <ref> and those studied in this paper are non-adaptive.]

Notably, the formulation of Problem <ref> makes no assumptions on $\bfA$ and $\bfb$. For example, we do not assume that there exists a ground truth $\bar{\bfx}$ and that  $\bfA\bar{\bfx} - \bfb$ is bounded in magnitude, or follows some distribution (e.g., has random Gaussian entries). Under these stronger assumptions, much is known about the problem, which has been studied for decades in the statistics literature on “optimal design of experiments”, as well as in machine learning [41, 55, 13].

In contrast, progress on the assumption-free version of the problem has only come in recent years, thanks to advances in random matrix theory and randomized numerical linear algebra. This is for good reason: solving Problem <ref> inherently requires choosing which entries of $\bfb$ to query in a randomized way: an adversary can easily “fool” any deterministic algorithm by concentrating error in $\bfA\bfx - \bfb$ on the indices of $\bfb$ that will be deterministically queried.

§.§ Prior Work

Euclidean Norm.
Problem <ref> is fully understood when the error is measured in the $\ell_2$ norm, $\|\bfw\|_2 = \left(\sum_{i=1}^n |\bfw_i|^2\right)^{1/2}$ – i.e., for least squares regression. The typical approach is to subsample and reweight rows (i.e., constraints) of the regression problem and to let $\tilde{\bfx}$ be the minimizer of this sampled problem, which only involves a fraction of the entries in $\bfb$. I.e., letting $\bfS \in \R^{m\times n}$ be a sampling matrix with $m < n$ rows ($\bfS$ has one non-zero entry per row), set $\tilde{\bfx} = \argmin_\bfx \|\bfS \bfA\bfx - \bfS\bfb\|$. When constraints are selected with probability proportional to the statistical leverage scores of $\bfA$'s rows, Problem <ref> can be solved with $O(d/\epsilon \cdot \log d)$ samples, and thus $O(d/\epsilon \cdot \log d)$ queries to $\bfb$ [56, 73, 28].[All query complexity bounds in this section are stated for solving Problem <ref> with high constant probability – e.g., probability $99/100$. In later sections we will include an explicit dependence on a failure probability $\delta$.]
Using tools from spectral graph sparsification [10, 45], Chen and Price recently improved the leverage score sampling result to $O(d/\epsilon)$, which is optimal [18].

In practice, methods based on leverage score sampling (also known as “coherence motivated sampling”) have found many applications. They are widely used in high-dimensional function fitting problems arising in the solution of parametric partial differential equations, where even mild assumptions on $\bfA$ and $\bfb$ are undesirable [12, 16, 37]. Methods for solving Problem <ref> in the $\ell_2$ norm also yield robust methods for interpolating sparse Fourier functions, bandlimited and multiband functions, and for data-efficient kernel learning [14, 2, 29].

Other Norms. Much less was known about Problem <ref> beyond the $\ell_2$ norm until recent work of Chen and Dereziński [11], which proves an upper bound of ${O}(d /\epsilon^2 \cdot \log d)$ queries for the $\ell_1$ norm, $\|\bfw\|_1 = \sum_{i=1}^n |\bfw_i|$. This result is tight up to the $\log d$ factor. A similar result is obtained in [53]. Chen and Dereziński also prove a result for $\ell_p$ norms, $\|\bfw\|_p = \left(\sum_{i=1}^n |\bfw_i|^p\right)^{1/p}$, for $p\in (1,2)$, in which they show that $O(d^2/\epsilon^2 \cdot \log d)$ queries suffice to solve Problem <ref>. As for the $\ell_2$ norm, the results for $\ell_1$ and $\ell_p$ are obtained by subsampling rows of the regression problem independently at random. However, instead of sampling with probabilities proportional to the leverage scores, [11, 53] employ a natural generalization of these scores known as the $\ell_p$ Lewis weights [17]. They left open the question of whether a linear in $d$ dependence is possible for $1 < p < 2$, and any bounds at all for $p > 2$.

Beyond norms, if $\bfb$ is a $\{-1,1\}$ label vector, and the error is measured via the logistic loss, Munteanu et al. [50] show that $\poly(d, \mu, 1/\epsilon)$ samples suffice, where $\mu$ is a complexity measure of $\bfA$.
This bound has recently been tightened to $\tilde O(d \mu^2/\epsilon^2)$ [49], using Lewis weight sampling.[Throughout, $\tilde{O}$ is used to suppress polylogarithmic factors in the argument.] For other loss functions, such as the Tukey loss and Huber's $M$-estimators for robust regression [32], we are not aware of any known results solving Problem <ref>. Chen and Dereziński also pose the open question of obtaining active regression bounds for other loss functions, in particular the Tukey and Huber losses, which are important in practice.

§.§ Our Contributions

$\ell_p$ Active Regression.
Our first main result is a new algorithm for solving Problem <ref> for the $\ell_p$ norm for any $0 < p < \infty$[Note that for $p\in(0,1)$, $\norm*{\cdot}_p$ is not a norm, but we refer to it as a norm by a standard abuse of notation.]. While near-optimal bounds are known for $p\in\{1,2\}$ [18, 11, 53], the problem is far from settled for all other $p$. Previously, active $\ell_p$ regression for $p>2$ and $0<p<1$ had no known nontrivial algorithms with $(1+\eps)$ relative error, and the only known approach was to read all $n$ entries of $\bfb$ and solve the problem using offline results. A natural question is whether a sublinear query complexity is possible in these regimes. For $p\in(1,2)$, [11] achieved an algorithm making $O(d^2/\eps^2\cdot \log d)$ queries, thus achieving the first sublinear query complexity. One of their main open questions is whether the dependence on $d$ can be improved to linear or not. Our main result answers all of these questions.

Given $0 < p < \infty$, $\bfA \in \R^{n \times d}$, and query access to $\bfb\in \R^n$, there is an algorithm (Algorithm <ref>) that solves Problem <ref> for the $\ell_p$-norm with probability $99/100$ which makes $m$ queries in $\bfb$, where

\[
m = \begin{dcases}
O\parens*{\frac{d}{\eps^{2}}(\log d)^2(\log(d/\eps))} & p\in (0,1) \\
O\parens*{\frac{d}{\eps}(\log d)^2(\log(d/\eps))} & p\in (1,2) \\
O\parens*{\frac{d^{p/2}}{\eps^{p}}(\log d)^2(\log(d/\eps))^{p-1}} & p\in(2,\infty)
\end{dcases}.
\]

Our main algorithm, Algorithm <ref>, is introduced and analyzed in Section <ref>, culminating in Theorem <ref>. We complement our algorithmic result with various new lower bounds which show the tightness of our algorithm, proven in Section <ref>. For $p\in(0,2)$, our dependence on $d$ and $\eps$ in the query complexity are simultaneously tight up to polylogarithmic factors; we show an $\Omega(d/\eps^2)$ lower bound for $p\in(0,1)$ and an $\Omega(d/\eps)$ lower bound for $p\in(1,2)$. For $p>2$, our dependence on $d$ is tight due to a lower bound of $\Omega(d^{p/2})$ which we show, while our $\eps$ dependence is off by at most factor of $\eps$ due to an $\Omega(\eps^{1-p})$ lower bound for the one-dimensional $\ell_p$ power means problem in Theorem 3 of [19]. Note that our active regression lower bounds for $p\in(0,2)$ improve this previous power means lower bound.

Notably, we achieve a linear dependence on $\eps$ for $p\in(1,2)$, which is perhaps surprising given that all previous known approaches to dimension reduction for $\ell_p$ regression relied on preserving the $\ell_p$ norm of all vectors in a subspace up to $(1\pm\eps)$ factors [17], which requires $\Omega(d/\eps^2)$ dimensions [47]. It also demonstrates a separation in the query complexity for $p \leq 1$ and $1 < p < 2$, due to a lower bound of $\Omega(d/\eps^2)$ for $p = 1$ [11, 53] as well as for $p\in(0,1)$ which we show.

Note that Theorem <ref> is stated to solve Problem <ref> with constant probability, $99/100$. In general, we show how to obtain $1-\delta$ probability with dependence on $\delta$ that is only polylogarithmic in $1/\delta$. In Section <ref>, we show that any algorithm that simply samples rows of the regression problem and solves the sampled problem must suffer a $1/\delta^{p-1}$ dependence. Indeed, such a loss is seen in the algorithm of [11] for $p\in(1,2)$. Thus, a success probability boosting routine, as we give in Section <ref>, is required to obtain a $\poly\log(1/\delta)$ dependence.

Upper and lower bounds for Problem <ref> for various norms and loss functions. New results are highlighted in blue.

For simplicity, we suppress leading constants depending only on $p$, as well as $\poly\log n$ factors for $M$-estimator results. Our results significantly strengthen and generalize prior work, providing the first query complexity result with a tight $d$ dependence for $\ell_p$ norms. We also give the first results for $M$-estimators as well as $\ell_p$ norms for $p>2$ and $p\in(0,1)$, and matching lower bounds in many cases.

Loss Function

2c|Prior Work      2c|Our Work     2c|Lower Bound

$\ell_2$					    $O(d/\epsilon)$    [18]     2c|–     $\Omega(d/\epsilon)$     [18]

$\ell_1$ 						    $\tilde{O}(d/\epsilon^2)$    [11]     2c|–     $\Omega(d/\epsilon^2)$    [11]

$\ell_p$, $p\in (1,2)$     $\tilde{O}(d^2/\epsilon^2)$    [11]     blue!15$\tilde{O}(d/\eps)$    blue!15 (Thm. <ref>)    blue!15$\Omega(d/\eps)$    blue!15 (Thm. <ref>)

$\ell_p$, $p > 2$ 		    2c|–    blue!15$\tilde{O}(d^{\frac{p}{2}}/\eps^p)$     blue!15 (Thm. <ref>)    blue!15$\Omega(d^{\frac{p}{2}}+\eps^{1-p})$    blue!15 (Thm. <ref>)

$\ell_p$, $p \in (0,1)$ 		    2c|–    blue!15$\tilde{O}(d/\eps^2)$     blue!15 (Thm. <ref>)    blue!15$\Omega(d/\eps^2)$    blue!15 (Thm. <ref>)

$M$-estimators		    2c|–    blue!15$\tilde{O}(d^{\frac{p}2+O(1)}/\eps^c)$     blue!15 (Thm. <ref>)    $\Omega(d)$

Huber loss		    2c|–    blue!15$\tilde{O}(d^{4-2\sqrt2}/\eps^c)$     blue!15 (Thm. <ref>)    $\Omega(d)$

Tukey loss		    2c|–    blue!15$\tilde{O}(d^{\frac{p}{2}+O(1)}/\eps^c)$     blue!15 (Thm. <ref>)    $\Omega(d)$

Sensitivity Bounds and Active Regression for General Losses.
We show that our approach to solving Problem <ref> for $\ell_p$ norms generalizes to a broad class of loss functions known as $M$-estimators [22], which take the form $\sum_{i=1}^n M([\bfA\bfx - \bfb]_i)$. The only properties that we require are that we can (1) compute a constant factor approximation to Problem <ref> (2) the loss function obeys approximate variants of the triangle inequality and (3) we can bound the so-called sensitivities of the loss, which bound the fraction of the total loss that can be concentrated at any coordinate $i\in[n]$ (see Equation (<ref>)).

To the best of our knowledge, the only prior result achieving sensitivity bounds for general loss functions is [66]. However, this work makes use of Löwner-John ellipsoids, which leads to practically inefficient algorithms, and loses a factor of $\sqrt d$ in the total sensitivity due to the ellipsoidal rounding. As our second main result, we develop new sensitivity bounds for $M$-estimators that significantly simplify and improve this result.

Let $\bfA\in\mathbb R^{n\times d}$ and let $M$ be an $M$-estimator loss with at most degree $p$ growth. Then, with probability at least $99/100$, Algorithm <ref> computes $M$-sensitivity upper bounds which sum to at most $O(d^{1\lor (p/2)}\log^2 n + \tau)$[Here, $a\lor b$ denotes $\max(a,b)$, and $a\land b$ denotes $\min(a,b)$.] in time at most $\tilde O(\nnz(\bfA) + nd^{C}/\tau)$ for some $C = O(1)$.

Our approach to sensitivity bounds only relies on hashing and the computation of $\ell_p$ Lewis weights [17, 31], and avoids the computation of Löwner-John ellipsoids. This allows for input sparsity time algorithms, and answers an open question of [66] on avoiding Löwner–John ellipsoids in the computation of sensitivities. Note that our dependence on $d$ matches the sensitivity bounds for the $\ell_p$ loss and is thus tight. We also show in Section <ref> that the dependence on $n$ is necessary for loss functions such as the Huber and Tukey losses. Furthermore, our algorithm can be turned into a non-algorithmic proof that the sensitivities sum to at most $O(d^{1\lor(p/2)}\log n)$ for these $M$-estimators as shown in Section <ref>; this is in fact tight for the Tukey loss by our lower bound of $\Omega(d\log n)$ in Section <ref>. Thus, we obtain the first tight bounds on the sum of sensitivities, for losses other than $\ell_p$. Overall, we make significant progress on generalizing the theory of matrix approximation beyond $\ell_p$ losses to handle general $M$-estimators, which is a direction that has recently received much attention [33, 22, 21, 23, 62, 66].

Combined with our active regression techniques, our sensitivity bounds yield active regression algorithms for general loss functions, including the Huber and Tukey losses, answering an open question of Chen and Dereziński [11]. Note that prior to our work, no sublinear query complexity was known for any $M$-estimator regression, besides the $\ell_2$ and $\ell_1$ losses.

Furthermore, our new sensitivity bounds imply significant improvements in previous results using sensitivity sampling, beyond active regression, including Orlicz norm subspace embeddings [61] and robust subspace approximation [21]. We believe that our general technique here will find other further applications, and leave it as an open question to do so.

Our new sensitivity computation algorithm for general losses is given in Algorithm <ref>, and its guarantees are stated and proven in Theorem <ref>. Its application to active regression is given in Theorem <ref>, its applications to Orlicz norm subspace embeddings are discussed in Section <ref>, and its applications to robust subspace approximation are discussed in Section <ref>.

Subspace Embeddings for Orlicz Norms.

Orlicz norms can be viewed as scale-invariant extensions of $M$-estimators, and have recently attracted attention as a general class of norms that admit efficient dimensionality reduction results [4, 61]. In particular, [61] apply sensitivity sampling to obtain subspace embeddings for Orlicz norms, which yields a small weighted subset $\tilde\bfA$ of rows of a matrix $\bfA\in\mathbb R^{n\times d}$ such that $\norm{\tilde\bfA\bfx} = (1\pm\eps)\norm{\bfA\bfx}$[For $a,b\geq 0$, $a \pm b$ denotes a number $c$ such that $a-b \leq c \leq a+b$.] for all $\bfx\in\mathbb R^d$. However, the number of rows required by [61] is a large polynomial in $d$, and is also restricted to Orlicz norms of at most quadratic growth. We show that by applying our new sensitivity bounds, we can obtain subspace embeddings for Orlicz norms with $d^{2\lor(p/2+1)}\poly(\log n,\eps^{-1})$ rows, for any Orlicz norm with a polynomial growth bound of degree $p$.

Robust Subspace Approximation.

The robust subspace approximation problem generalizes the classical low rank approximation problem of finding a rank $k$ projection $\bfX$ minimizing $\norm*{\bfA\bfX - \bfA}_F$ by replacing the Frobenius norm with an extension of $M$-estimators to matrix norms. [21] showed the first dimensionality reduction results for this problem for a general class of $M$-estimators of at most quadratic growth via a recursive sampling scheme using the sensitivity sampling framework. However, due to the use of looser sensitivity bounds, they suffer an undesirable factor of $(\log n)^{O(\log k)}$ in their sample complexities. Our new sensitivity bounds allow us to remove this factor, giving a dimension reduction result into a $\poly(k,\log n,\eps^{-1})\times \poly(k,\log n,\eps^{-1})$ instance. We also extend their method beyond quadratic growth, to any degree $p$ polynomial growth.

Active Regression for the Huber Loss.
Our active regression result for general $M$-estimators discussed above is loose by a factor of $d$ in the sample complexity, compared to our $\ell_p$ active regression results.
This is attributed to the use of our net argument for general $M$-estimators, whereas our $\ell_p$ active regression results can make use of more sophisticated chaining arguments of [6, 46, 63]. A natural question is if this gap can be improved.

We consider the important special case of the Huber loss, which is defined as follows:

The Huber loss of width $\tau\geq 0$ is defined as

\begin{align*}
H(x) &\coloneqq \begin{cases}
x^2 / 2\tau  & \text{if $\abs{x}\leq \tau$} \\
\abs{x} - \tau/2 & \text{otherwise}
\end{cases}
\end{align*}

and the Huber norm[Again, this is a standard abuse of notation, and the Huber norm is not an actual norm.] is defined as $\norm*{\bfy}_H \coloneqq \sqrt{\sum_{i=1}^n H(\bfy(i))}$.

The Huber loss is “arguably one of the most widely used $M$-estimators” [22], owing its popularity to its convexity and differentiability properties of $\ell_2$, which allows for efficient algorithms (see, e.g., [48] for algorithms), in combination with its robustness properties of $\ell_1$ [35]. This makes it widely applicable in practical big data settings (see, e.g., [9] for a list of popular software packages implementing Huber regression as well as references that make use of Huber regression). Variations on Huber regression have also recently been shown to hold theoretical guarantees in the robust statistics literature (see, e.g., [43, 44] and references therein).

For the Huber loss, we show that it is indeed possible to leverage the chaining techniques in order to obtain improved sample complexity bounds for active regression. We show that we can improve beyond the $d^2$ bound obtained by our general $M$-estimator algorithm as applied to the Huber loss, and obtain a sample complexity of $O(d^{4-2\sqrt 2}\poly(\log n,\eps^{-1}))$ queries to $\bfb$, where $4-2\sqrt2 \approx 1.17157$. For this result, we use the chaining techniques of [6], which provides a more flexible alternative to [46], but requires more technical effort to adapt to the active setting.

[Main Result for Huber Active Regression]theoremActiveHuber
Let $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$. Then, with probability at least $99/100$, Algorithm <ref> returns a $\tilde\bfx$ satisfying

\[
\norm{\bfA\tilde \bfx - \bfb}_H \le (1+\epsilon) \cdot \min_\bfx\norm*{\bfA\bfx - \bfb}_H
\]

Furthermore, the algorithm reads at most $d^{4-2\sqrt2}\poly(\log n,\eps^{-1})$ entries of $\bfb$.

Our techniques also yield a subspace embedding result, which constructs a weighted subset $\tilde\bfA$ of $O(d^{4-2\sqrt 2}\poly(\log n,\eps^{-1}))$ rows such that $\norm{\tilde\bfA\bfx}_H = (1\pm\eps)\norm{\tilde\bfA\bfx}_H$ for all $\bfx\in\mathbb R^d$, contained in Theorem <ref>. Previously, the best known dimension reduction bound for Huber regression, even in the non-active setting, was $d^4$ due to [22].

Furthermore, this is, to the best of our knowledge, the first example of a loss function other than $\ell_p$ which achieves a sensitivity sampling bound of better than $d^2$, despite the fact that such results have been sought in many works [22, 21, 61, 23, 66, 34]. The reason for this is that $d^2$ is a natural bound for sensitivity sampling, attributed to one $d$ factor from the sum of sensitivities and one $d$ factor from carrying out a union bound over a net of $\exp(d)$ vectors. For $\ell_p$ norms, the arguments of [6, 63] and their subsequent improvements avoid this problem by using a more sophisticated chaining argument. However, these arguments use the structure of $\ell_p$ spaces in crucial ways, such as isometric changes of density using Lewis weights [40], and do not generalize easily to other loss functions.

It is an interesting open question to determine whether our dimension reduction bound for the Huber loss can be improved all the way down to $d$.
Our results for the Huber loss are found in Section <ref>.

Dimension Reduction for Gamma Functions for Faster $\ell_p$ Regression.

One particularly important application of sampling-based dimension reduction for loss functions beyond $\ell_p$ losses is, perhaps surprisingly, in the design of algorithms for $\ell_p$ regression. The work of [5] introduces gamma functions $\gamma_p$, which are generalizations of the Huber loss which behave quadratically near the origin and like $|x|^p$ away from the origin, in the context of algorithms for $\ell_p$ regression. Subsequently, [3] obtained even faster algorithms by using constant factor approximations of $\gamma_p$ regression as a subroutine, in which the $\gamma_p$ loss is minimized over a subspace. Dimension reduction for this loss function has been a crucial ingredient for recent results in fast algorithms for $\ell_p$ regression [1, 34]. In particular, [1] highlighted the open question of designing sparsification methods for $\gamma_p$ functions for $p\in(1,2)$, and [34] designed a sampling algorithm which samples $\tilde O(d^3)$ rows. By generalizing our dimension reduction techniques for the Huber loss, we obtain an algorithm which samples at most $O(d^{4-2\sqrt 2}\poly(\log n,\eps^{-1}))$ rows for any $p\in[1,2)$, and improves to $O(d\poly(\log n,\eps^{-1}))$ rows as $p\to 2$ (see Figure <ref> for the trade-off curve). We give a further discussion in Section <ref>.

Kronecker Product Regression. Beyond applications in data-efficient regression, Theorem <ref> implies the first sublinear time algorithm for Kronecker product regression in any $\ell_p$ norm, where explicitly constructing the vector $\bfb$ is a computational bottleneck. We detail this result in Section <ref>.
In $q$-th order Kronecker product regression, one is
given matrices $\bfA_1, \bfA_2,\ldots, \bfA_q$, where $\bfA_i \in \mathbb{R}^{n_i \times d_i}$,
as well as a vector $\bfb \in \mathbb{R}^{n_1 n_2 \cdots n_q}$, and
the goal is to solve:
$\min_{\bfx \in \R^{d_1 d_2 \cdots d_q}} \|(\bfA_1 \otimes \bfA_2 \cdots \otimes \bfA_q)\bfx - \bfb\|_p,$
where $\otimes$ denotes the Kronecker product. Typically $\prod_{i=1}^q d_i$ is much less
than $\prod_{i=1}^q n_i$, and the goal is to obtain algorithms that do not explicitly form
$\bfA_1 \otimes \bfA_2 \otimes \cdots \otimes \bfA_q$ or $\bfb$, which is too expensive. Our results yield the first algorithm for Kronecker product regression,
for every $p \geq 1$, whose running time does not depend on $\nnz(\bfb)$, whereas previous results had a linear dependence on $\nnz(\bfb)$, which can be as large as $\prod_{i=1}^q n_i$ [25].

Let $q \geq 1$, $p \geq 1$ be constant, and $\epsilon > 0$. Kronecker product regression
can be solved up to a $(1+\epsilon)$-factor with constant probability in
$\tilde{O}(\sum_{i=1}^q \nnz(\bfA_i) + \poly(\prod_{i=1}^q d_i / \epsilon))$ time.

§.§ Technical Approach

§.§.§ lp Active Regression

Our algorithm for solving Problem <ref> uses a novel variation on the “sample-and-solve” approach. In particular, we randomly select a row sampling matrix $\bfS\in \R^{m\times n}$ and return $\tilde{\bfx} = \argmin_x\|\bfS \bfA\bfx - \bfS\bfb\|$, which only requires querying $m$ entries of $\bfb$ (those that appear in $\bfS\bfb$).
To get tight bounds for $\ell_p$ regression, we select $\bfS$ using $\ell_p$ Lewis weight sampling, a generalization of leverage score sampling for $\ell_2$.

It can be shown (Lemma <ref>) that the $\ell_p$ Lewis weights upper bound the $\ell_p$ sensitivities of $\bfA$, a measure of importance for the rows of $\bfA$. The $\ell_p$ sensitivity of the $i^\text{th}$ row of $\bfA$ is defined as

\[
\bfs^p_i(\bfA) \coloneqq \max_{\bfx \in \R^d\setminus\{0\}} \frac{|[\bfA\bfx](i)|^p}{\norm{\bfA\bfx}_p^p},
\]

where $[\bfA\bfx](i)$ denotes the $i^\text{th}$ entry of the vector $\bfA\bfx$, and captures how large the $i^\text{th}$ entry of any $\bfA\bfx\in\colspan(\bfA)$ can be, relative to the $\ell_p$ norm. A standard scalar Bernstein bound shows that if $\bfS$ samples rows with probabilities that upper bound the sensitivities, then $\|\bfS \bfA\bfx\|_p^p = (1\pm\epsilon)\|\bfA\bfx\|_p^p$ with high probability, for each $\bfx\in\mathbb R^d$. An $\epsilon$-net argument can extend this to a for all claim.

Prior Approaches to $\ell_p$ Active Regression.

While the above ideas give an approach for standard $\ell_p$ regression, this bound does not suffice for active $\ell_p$ regression. To solve Problem <ref>, we actually want that $\|\bfS (\bfA\bfx - \bfb)\|_p^p = (1\pm\epsilon)\|\bfA\bfx-\bfb\|_p^p$ for any $\bfx$. Will $\bfS$ provide such a guarantee? The main problem, as discussed in [11, 53] is that the translation by $\bfb$ may introduce outliers, i.e., entries with high sensitivity which are not captured by the sensitivity scores of $\bfA$. As shown by [11, 53], in the case of $\ell_1$, the special structure of the loss function provides a solution. Indeed, by the triangle inequality,

\[
\abs*{\parens*{\abs*{[\bfA\bfx-\bfb](i)} - \abs*{[\bfA\bfx^*-\bfb](i)}}} \leq \abs*{[\bfA\bfx-\bfA\bfx^*](i)} = \abs*{[\bfA(\bfx-\bfx^*)](i)}
\]

where $\bfx^*$ is the optimal solution. This fact can be used to show that sampling by the sensitivities of $\bfA$ preserves the differences between the cost of any $\bfx$ and the optimal $\bfx^*$. However, such a proof cannot work for $p\neq 1$, in which case we do not have such a nice inequality. For $p\in(1,2)$, [11] take the approach of bounding the residual error terms from the above approach by using a Taylor approximation, but this leads to a sample complexity of at least $d^2$.

Our Solution: Partitions by Sensitivity.

Instead of relying on the technique of “cancelling out the outliers”, we take a conceptually different approach. We proceed in two stages, where we (1) first find a constant factor solution $\bfx_c$ such that $\norm*{\bfA\bfx_c-\bfb}_p^p \leq O(1)\cdot\min_\bfx\norm*{\bfA\bfx-\bfb}_p^p$ using an idea of [24] and replace $\bfb$ by the residual vector $\bfb - \bfA\bfx_c$, and then (2) conceptually partition the target vector $\bfb$ into two sets of coordinates, the coordinates $i\in[n]$ that are small enough to be comparable to the sensitivity $\bfs_i^p(\bfA)$ and those that are much larger. That is, we consider the coordinates $i\in[n]$ such that $|\bfb(i)|^p / \|\bfb\|_p^p \leq C\cdot \bfs^p_i(\bfA)$ for some $C>0$, and all other coordinates. For the former set of coordinates, one can check that the Bernstein bound still applies, and $\bfS$ does preserve the norm of $\|\bfA\bfx - \bfb\|_p^p$, when restricted to these coordinates. On the other hand, for the latter set of coordinates, we show that no vector of the form $\bfA\bfx$ can both be close to $\bfb(i)$ in its $i^\text{th}$ entry, and still close to the remainder of $\bfb$ – the $i^\text{th}$ entry is simply too large in magnitude. In particular, to have $[\bfA\bfx](i)$ close to $\bfb(i)$, we would require $\norm{\bfA\bfx}_p^p$ to be much larger than $\norm{\bfb}_p^p$, which by our preprocessing step, is on the order of the optimal cost $\min_\bfx \|\bfA\bfx - \bfb\|_p$. Via the triangle inequality, this implies that $\bfA\bfx$ must be far from an optimal solution. Thus, we can argue that any near-optimal solution to $\min_\bfx \|\bfA\bfx - \bfb\|_p$ does not need to fit $\bfb(i)$ with $|\bfb(i)|^p / \|\bfb\|_p^p$ much larger than $\bfs^p_i(\bfA)$. We can effectively ignore the contribution of these rows.

Another technical challenge remains: to obtain an optimal dimension dependence, we need a refined $\epsilon$-net argument to make for all statements about $\bfx\in \R^d$. To do so, we adapt the chaining arguments of Bourgain, Lindenstrauss, and Milman [6] and Ledoux and Talagrand [46] to the active regression setting, avoiding the standard $\eps$-net and union bound argument used by, e.g., [57]. Although both [6] and [46] provide such approaches, we adapt the (slightly) more complex recursive Lewis weight sampling algorithm of [46] in order to obtain tighter dependencies on $\eps$. The streamlined proof of [46] also adapts nicely to the active regression setting with minimal changes to the original argument. We note here that we will later also need to adapt the much more involved [6] argument to handle the Huber loss, in which case the proof of [6] allows for more fine-grained control over bounding the sensitivity sampling algorithm, but requires a more complex argument based on carefully partitioning the coordinates of the target vector $\bfb$ based on sensitivity weight classes.

Aside from our new application of [6, 46], we hope that by translating the arguments of [6, 46] to the language of theoretical computer science and matrix approximation, they will find further applications to randomized algorithm design.

We note that our algorithm is quite a bit more involved than a simple scheme of sampling proportionally to Lewis weights and solving. This is for good reasons. Not only is it not clear that such an approach works at all, we show in Theorem <ref> that for any $p>1$, any algorithm which simply samples reweighted rows and solves the system must have a polynomial dependence on $1/\delta$ in the query complexity, while our algorithm achieves a $\log\frac1\delta$ dependence, by solving residual problems of a constant factor solution. Thus, our two-stage approach is necessary to achieve our $\delta$ dependence. Furthermore, the best known analysis of a simple “one-shot” Lewis weight sampling scheme suffers in $\eps$ dependencies for $p>2$, where the one-shot approach is only known to give a $\tilde O(d^{p/2}/\eps^5)$ bound for subspace embeddings [6, 17], whose losses translate to losses for our active regression algorithms as well, while the recursive approach can achieve $\tilde O(d^{p/2}/\eps^2)$ [46]. While the [46] result is an existential result, we provide an analysis of the [46] proof in Theorem <ref> to turn it into a randomized algorithm with logarithmic dependencies on the failure rate $\delta$, which achieves the best known dependence on $d$, $\eps$, and $\delta$, up to logarithmic factors. We further modify this subspace embedding result for active regression, to optimize our $\eps$ dependence.

Optimized $\eps$ Dependence for $p\in(1, 2)$.

For $p\in(1,2)$, the above argument gives a bound of $\tilde O(d/\eps^2)$. While the linear dependence on $d$ is optimal, confirming the conjecture of [11], it has a quadratic dependence on $\eps$, which is in fact not optimal. We now show how to improve our bound to $\tilde O(d/\eps)$, which requires additional ideas. We first use strong convexity to show that a $(1+\gamma)$-approximate solution $\hat\bfx\in\mathbb R^d$ satisfying $\norm*{\bfA\hat\bfx-\bfb}_p^p \leq (1+\gamma)\norm*{\bfA\bfx^*-\bfb}_p^p$, for the optimal solution $\bfx^*\in\mathbb R^d$, in fact satisfies $\norm*{\bfA\hat\bfx-\bfA\bfx^*}_p \leq O(\sqrt\gamma)\norm*{\bfA\bfx^*-\bfb}_p$. Then, by using that $\hat\bfx$ is close to the optimal solution, we show an improved bound on the difference in the objective values of $\hat\bfx$ and $\bfx^*$, i.e., that $\hat\bfx$ actually has an approximation ratio better than $(1+\gamma)$. We then iterate this argument until we obtain a $(1+\eps)$-approximation using $\tilde O(d/\eps)$ queries, at which point we can no longer get improvements. The chaining argument used in this proof, while similar to the previous proofs, has a different geometry than the previous chaining arguments, and requires additional ideas.

Our upper bound is tight up to polylogarithmic factors due to a lower bound we show in Theorem <ref>. This improves an $\ell_p$ power means lower bound of [19], who only showed a lower bound of $\Omega(\eps^{1-p})$ queries. Unfortunately, we are unable to port our algorithmic techniques to the $\ell_p$ power means problem in high dimensions, due to difficulties in adapting their chaining argument.

§.§.§ Sensitivity Bounds

The notion of $\ell_p$ sensitivities, as discussed above, naturally generalizes to loss functions that take the form of coordinate-wise sums. Consider a loss function $M$ and an $n\times d$ matrix $\bfA$. Then, the sensitivity of the $i$th coordinate with respect to the loss function $M$ is defined as

\begin{equation}\label{eq:sensitivity-def}
\bfs_i^M(\bfA) \coloneqq \sup_{\bfx\in\mathbb R^d\setminus\{0\}}\frac{M([\bfA\bfx](i))}{\sum_{j=1}^n M([\bfA\bfx](j))}.
\end{equation}

It is well-established that sensitivities provide a general framework for sampling rows of $\bfA$ that approximate $\bfA$ well under the loss function $M$ [30]. While a rich literature exists for $\ell_p$ [24, 58, 17], little was known about the approximation of sensitivities for general loss functions until [66], which used Löwner-John ellipsoids to obtain sensitivity bounds for a general family of near-convex losses. However, the computation of Löwner-John ellipsoids has running time that is a large polynomial in $n$ and $d$, and is impractical for large datasets, and [66] raise the open question of obtaining general sensitivity bounds without this expensive subroutine.

Our approach to new sensitivity bounds significantly generalizes the approach of [23], whose algorithm can be seen as a way to use hashing and Lewis weights to compute sensitivities for the Tukey loss, but heavily uses the properties of the Tukey loss in their analysis.

Suppose that a coordinate $i\in[n]$ has $M$-sensitivity $\alpha\in(0,1]$, that is,

\[
\bfs_i^M(\bfA) = \sup_{\bfx\in\mathbb R^d\setminus\{0\}}\frac{M([\bfA\bfx](i))}{\sum_{j=1}^n M([\bfA\bfx](j))} = \alpha,
\]

and let $\bfy = \bfA\bfx$ witness this supremum, and assume for simplicity that $\sum_{i=1}^n M([\bfA\bfx](j)) = 1$. Note then that there can be at most $1/\alpha$ entries $j\in[n]$ of $\bfy$ that have coordinate value $M(\bfy_j)\geq M(\bfy_i) = \alpha$. Then, if we randomly hash the $n$ coordinates into $O(1/\alpha)$ buckets, then with constant probability, coordinate $i$ will be isolated from any other entry with $M(\bfy_j)\geq M(\bfy_i) = \alpha$. Now if $M$ is monotonic, then this means that $\bfy_i$ is the largest coordinate in its hash bucket. Furthermore, the sum of the $M$-mass of all of the other coordinates in $i$'s hash bucket is only an $\alpha$ fraction of the total $M$-mass, so entry $i$ carries a constant fraction of the $M$-mass in its bucket. In this case, it can be shown that entry $i$ must in fact carry a constant fraction of the $\ell_2$ mass inside its hash bucket, if $M$ is a function of at most quadratic growth. This is because when we switch the error metric from $M$ to $\ell_2$, then the largest entry will have the largest increase in its normalized contribution (see Lemma <ref>). This means that row $i$ must have an $\ell_2$ leverage score of $\Omega(1)$, in this hash bucket.

This leads to the following algorithm: (1) hash the $n$ coordinates into $O(1/\alpha)$ buckets (2) compute $\ell_2$ leverage scores for each bucket (3) assign an $M$-sensitivity of $\alpha$ for any coordinate that has $\ell_2$ leverage score $\Omega(1)$. In each of the $O(1/\alpha)$ buckets, we will find at most $O(d)$ coordinates with leverage score at least $\Omega(1)$, so we assign an $M$-sensitivity of $\alpha$ to at most $O(d/\alpha)$ coordinates, which has a total sensitivity contribution of $O(d)$. By repeating this for $O(\log n)$ guesses of $\alpha$ in powers of $2$, this gives a total sensitivity bound of $O(d\log n)$. The constant probability events in the hashing process can be boosted to probability $1-1/\poly(n)$ by repeating the procedure $O(\log n)$ times, which increases the total sensitivity to roughly $O(d\log^2 n)$.

By sampling according to these sensitivities and applying a union bound over a net, we obtain the first active regression algorithms for general loss functions in Theorem <ref>. Note that this result is made possible by a combination of both our new sensitivity bounds for $M$-estimators and our new active regression techniques as discussed in Section <ref>. Furthermore, we demonstrate other applications of our sensitivity bound result, showing how to improve Orlicz norm subspace embeddings in Section <ref> and robust subspace approximation in Section <ref>.

§.§.§ Subspace Embeddings and Active Regression for the Huber Loss

As discussed previously, we tackle the question of leveraging the theory of [6] nets in order to obtain sample complexities for the Huber loss beyond $d^2$. Our algorithmic framework for active regression is based on the earlier idea of partitioning the entries of $\bfb$ by sensitivity and then applying sensitivity sampling, so we focus on the problem of preserving the Huber norm using an improved sensitivity sampling technique. Note that unlike the $\ell_p$ losses, the Huber loss is not scale-invariant. Furthermore, perhaps the largest obstacle in designing row sampling algorithms for the Huber loss going beyond standard $\eps$-net arguments is that there is no analogue of the chaining constructions of [6, 63] for the Huber loss. This can also be attributed to the fact that the Huber loss is not scale-invariant, which precludes an isometric change-of-density type theorem for the Huber loss as done in [42, 63]. We show how to overcome these obstacles in the following discussion.

A Sharp Huber Inequality.

Our algorithmic framework follows the Huber algorithm of [22], which is a recursive sampling algorithm which reduces the number of rows from $n$ to roughly $n^{1/2}d^2$ in each recursive application of the algorithm. To show this result, [22] first show in their Lemma 2.1 that the Huber norm is within a factor of $O(n^{1/2})$ of the smaller of the $\ell_1$ and $\ell_2$ norms:

Let $\bfy\in\mathbb R^n$. Then,

\[
\norm*{\bfy}_H^2 = \sum_{i=1}^n H(\bfy_i) \geq \Omega(n^{-1/2})\min\braces{\norm*{\bfy}_1, \norm*{\bfy}_2^2}
\]

It can be shown that the above lemma implies that the Huber sensitivities are within a factor of $O(n^{1/2})$ of the sum of the $\ell_1$ and $\ell_2$ sensitivities.

This motivates the idea of sampling the rows of $\bfA$ with probability proportional to the sum of the $\ell_1$ and $\ell_2$ Lewis weights, oversampled by a factor of $O(n^{1/2})$. This is indeed how [22] proceeds.

The recursion $n\to n^{1/2}d^2$ solves to a final row count of around $d^4$, which is quadratically worse than our general loss function result of $d^2$ using our new sensitivity upper bounds and our general framework. To improve this further, first note that two improvements can be made to the above argument. First, by using the Huber inequality in a different way, we can use it in conjunction with the [6] net bounds, which reduces the row count in one recursive application to roughly $n^{1/2}d$ rather than $n^{1/2}d^2$ (see Lemma <ref>). This reduces the overall row count to $d^2$ after solving for the recursion, but this still does not beat our general purpose sensitivity sampling algorithm, despite the use of the [6] nets. The second improvement is that the Huber inequality as proved in [22] is in fact loose by a polynomial factor in $n$, and can be improved to the following:

Let $\bfy\in\mathbb R^n$. Then,

\[
\norm*{\bfy}_H^2 = \sum_{i=1}^n H(\bfy_i) \geq \Omega(n^{-1/3})\min\braces{\norm*{\bfy}_1, \norm*{\bfy}_2^2}
\]

This lemma is tight up to constant factors[Consider the vector with one coordinate with $n^{1/3}$ and $(n-1)$ coordinates with $n^{-1/3}$.], and gives a recursion of roughly $n\to n^{1/3}d$, giving

\[
O(d^{3/2}\poly(\eps^{-1},\log n))
\]

rows, which shaves a factor of approximately $\sqrt d$ over the naïve Bernstein bound over a net.

Storing Large Huber Sensitivities.

In order to further improve upon this bound, we crucially make use of our improved sensitivity bounds from Section <ref> and a generalized version of the above Huber inequality lemma that is parameterized by an upper bound on the size of the entries of $\bfy$.

Let $\bfy\in\mathbb R^n$ and let $0 < \gamma \leq 1$. Let

\[
T \supseteq \braces*{i \in [n] :H(\bfy_i) \leq \gamma\norm*{\bfy}_H^2}.
\]

Then, for some constant $c>0$, at least one of the following bounds holds:

\begin{align*}
\norm*{\bfy\mid_T}_H^2 = \sum_{i\in T} H(\bfy_i) &\geq c\frac{1}{(\gamma n)^{1/3}}\min\braces{\norm*{\bfy\mid_T}_1, \norm*{\bfy\mid_T}_2^2} \\
\norm*{\bfy}_H^2 = \sum_{i=1}^n H(\bfy_i) &\geq c\gamma\min\braces{\norm*{\bfy}_1, \norm*{\bfy}_2^2}.
\end{align*}

By directly including the rows of $\bfA$ with Huber sensitivity at least $\gamma$, we exactly preserve the Huber norm inside $\overline T = [n]\setminus T$ for every $\bfy$. On the remaining coordinates inside $T$, we then have an improved Huber inequality, which implies an improved sampling bound. By balancing the number of rows which we directly include, which is roughly $d/\gamma$, and the sampling bound inside $T$, which is roughly $(1/\gamma + (\gamma n)^{1/3})d$, we obtain a bound of roughly $n^{1/4}d$ rows by choosing $\gamma = n^{-1/4}$ at each step. By recursively applying this result, we obtain an improved sampling bound of

\[
O(d^{4/3}\poly(\eps^{-1},\log n)).
\]

Comparing to Every $p\in[1,2]$.

Finally, to achieve our final optimization, we further drive down the ratio between $\norm*{\bfy}_H^2$ and $\norm*{\bfy}_p^p$ by choosing the best $p\in[1,2]$ for each $\bfy$:

Let $\bfy\in\mathbb R^n$, $\alpha \in [0,1/2]$, $\gamma = n^{-\alpha}$ with $2/n \leq \gamma \leq 1$. Let

\[
T \supseteq \braces*{i\in[n] : H(\bfy_i) \leq \gamma \norm*{\bfy}_H^2}.
\]

Then, for some $c>0$ and $\beta = 3-2\sqrt 2\approx 0.17157$, at least one of the following bounds holds:

\begin{align*}
\norm*{\bfy\mid_T}_H^2 &\geq c\frac1{(\gamma n)^{\beta}}\min_{p\in[1,2]}\norm*{\bfy\mid_T}_p^p \\
\norm*{\bfy}_H^2 &\geq c\gamma\min_{p\in\{1,2\}}\norm*{\bfy}_p^p
\end{align*}

In fact, we prove a generalized bound for the $\ell_2$-$\ell_q$ loss for any $q\in(0,2)$ in Lemma <ref>. The interval $p\in[1,2]$ can be discretized in increments of $\frac1{\log n}$, so with $O(\log n)$ applications of [6] nets, we can always find a $p$ within an additive $\frac1{\log n}$ of the optimal $p$ for every net vector $\bfy$, which only affects Lemma <ref> by constant factors whenever $\bfy$ has entries bounded by $\poly(n)$. This is formalized in Corollary <ref>. By proceeding as previously discussed, we arrive at our final bound of

\[
O(d^{4-2\sqrt 2}\poly(\eps^{-1},\log n)).
\]

Extensions to $\ell_2$-$\ell_q$ Loss.

We generalize our results to the $\ell_2$-$\ell_q$ loss for $q\in(0,2)$. As $q$ ranges from $0$ to $1$ to $2$, the $\ell_2$-$\ell_q$ interpolates between the Tukey, Huber, and $\ell_2$ losses up to constant factors, and provides a natural generalization of these loss function.

Let $q\in(0,2)$ and define

\[
M(x) = \begin{cases}
\abs*{x}^2 & \text{if $\abs{x} \leq 1$} \\
\abs*{x}^q & \text{if $\abs{x} > 1$}
\end{cases}.
\]

Let $\bfy\in\mathbb R^n$ and let $\alpha \in [0,q/2]$ and $\gamma = n^{-\alpha}$ with $2/n \leq \gamma \leq 1$. Let

\[
T \supseteq \braces*{i\in[n] : M(\bfy_i) \leq \gamma \norm*{\bfy}_M^2}.
\]

Then, for some constant $c>0$, at least one of the following bounds holds:

\begin{equation}\label{eq:l2-lq-whole-vector}
\begin{aligned}
\norm*{\bfy}_M^2 &\geq c\gamma^{2/q - 1}\min_{p\in\{q,2\}}\norm*{\bfy}_p^p \\
\norm*{\bfy\mid_T}_M^2 &\geq c\frac1{(\gamma n)^{\beta}}\min_{p\in[q,2]}\norm*{\bfy\mid_T}_p^p
\end{aligned}
\end{equation}

\[
\beta = \frac1{2/q - 1}\bracks*{(2/q + 1) - 2\sqrt{2/q}}.
\]

A proof of this bound can be found in Appendix <ref>.

For $q\in[1,2)$, the $1/\gamma^{2/q-1}$ distortion in Equation (<ref>) is smaller than the $1/\gamma$ factor incurred from keeping $M$-sensitivities at least $\gamma$, so we can balance the parameters as $1/\gamma = (\gamma n)^\beta$ as before, which leads to a recursion that gives us a bound of $n = d^{1+\beta} \poly((\log n)/\eps)$. For $q\in(0,1)$, the $1/\gamma^{2/q-1}$ distortion in Equation (<ref>) is worse than $1/\gamma$, which means we must balance $1/\gamma^{2/q-1} = (\gamma n)^{\beta}$, or $\gamma = n^{-\beta/(2/q-1+\beta)}$, which gives a worse bound of $n = d^{\gamma}\poly((\log n)/\eps)$ for

\[
\gamma = \frac{2/q-1+\beta}{2/q-1+\beta(2-2/q)}.
\]

This is better than a $d^2$ bound as long as $q \geq (\sqrt 5-1)^2/8 \approx 0.19098$.

[->] (-0.5,0) – (2.2,0) node[right] $q$;
[->] (0,-0.5) – (0,2.2 - );
[domain=1:1.995,smooth,variable=,thick,line width=0.9mm,samples=100,blue] plot (,1 + ((1/(2/-1) * ((2/+1)-2 * sqrt(2/)))) - );
[domain=0.191:1,smooth,variable=,thick,line width=0.9mm,samples=100,red] plot (,(2/- 1 + ((1/(2/-1) * ((2/+1)-2 * sqrt(2/))))) / (2/- 1 + ((1/(2/-1) * ((2/+1)-2 * sqrt(2/)))) * (2 - 2/)) - );
[thick,line width=0.9mm] (0,2- )–(0.2,2- );

[below] (a) at (0.6,1.1) $d^{\gamma}$;
[below] (a) at (1.5,0.9) $d^{1+\beta}$;

[gray, dashed] (-0.5,1 - )–(2.2,1 - );
[left] (a) at (-0.5,1 - ) $d$;

[gray, dashed] (-0.5,2 - )–(2.2,2 - );
[left] (a) at (-0.5,2 - ) $d^2$;

[gray, dashed] (-0.5,1.172 - )–(2.2,1.172 - );
[left] (a) at (-0.5,1.172 - ) $d^{4-2\sqrt 2}$;

[gray, dashed] (1,-0.5)–(1,1.2 - );
[below] (a) at (1,-0.5) $1$;

[gray, dashed] (2,-0.5)–(2,1 - );
[below] (a) at (2,-0.5) $2$;

[gray, dashed] (0.191,-0.5)–(0.191,2 - );
[below] (a) at (0.191,-0.5) $\frac{(\sqrt 5-1)^2}{8}\approx 0.191$;

Dependence on $d$ for the active regression sample complexity for the $\ell_2$-$\ell_q$ loss. Similar bounds apply to subspace embeddings as well.

§.§ Conclusions and Future Directions

In this work, we study the sample complexity of active linear regression for both the $\ell_p$ norm as well as general $M$-estimator losses.

For the $\ell_p$ norm, we provide optimal algorithms and lower bounds for $p\in(0, 2)$, with $\tilde\Theta(d/\eps^2)$ samples for $p\in(0,1)$ and $\tilde\Theta(d/\eps)$ samples for $p\in(1,2)$. For $p>2$, we provide an upper bound of $\tilde O(d^{p/2}/\eps^p)$, which is optimal in the $d$ dependence and off by a single $\eps$ factor in the $\eps$ dependence, up to polylogarithmic factors. Our algorithms provide the first nontrivial bounds, i.e., sample complexity less than $n$, for $p\in(0,1)\cup(2,\infty)$, while for $p\in(1,2)$, we significantly improve upon the $\tilde O(d^2/\eps^2)$ upper bound of [11] and answer their main open question. We obtain these results via a two-stage algorithm and a novel sensitivity partitioning technique for every $p$, as well as an iterative improvement argument via strong convexity and Lewis bases to improve the $\eps$ dependence for $p\in(1,2)$. Our result is the first to achieve a linear dependence on $\eps$ for dimension reduction for $\ell_p$ regression for $p\in(1,2)$.

Next, we obtain a new sensitivity bound which achieves optimal total sensitivity bounds for $M$-estimators of at most polynomial growth, which runs in input sparsity time and avoids the use of Löwner–John ellipsoids. This answers an open question of [66] and makes significant progress in the general direction of matrix approximation beyond $\ell_p$ losses. By combining this with our new active regression techniques, we obtain active regression algorithms for general $M$-estimator losses, including the Tukey and Huber losses, which answers an open question of [11].

For the important special case of the Huber loss, we introduce new techniques which bound Huber sensitivities by the sum of $\ell_p$ Lewis weights, which allows us to take advantage of chaining arguments for $\ell_p$ in order to obtain an active regression algorithm making at most $O(d^{4-2\sqrt 2}\poly(\log n,\eps^{-1}))$ queries. Our techniques also give subspace embeddings with the same number of rows. This is the first dimension reduction result for losses other than $\ell_p$ to approximate a $d$-dimensional subspace with fewer than $d^2$ dimensions. This improves over a previous bound of $d^4$ for the Huber loss, which held only for subspace embeddings, and not active regression, in [22].

Finally, our results and techniques give many applications in a wide variety of related problems. Our lower bounds for active regression give improved lower bounds for the sublinear power means problem [19]; our new sensitivity bounding techniques sharpen and generalize previous results on Orlicz norm subspace embeddings [62] and robust subspace approximation [21]; our techniques for dimension reduction for the Huber loss gives improved bounds for sparsification for $\gamma_p$ functions for applications in fast algorithms for $\ell_p$ regression [34, 1]. We believe that our techniques will be applicable much further, and hope to see more uses in future work.

We conclude with questions that are still left open by our work. Perhaps the most pressing is to resolve the query complexity of active $\ell_p$ regression for $p>2$: our upper bound is $\tilde O(d^{p/2}/\eps^p)$, while the lower bound is $\Omega(d^{p/2}+\eps^{1-p})$. Closing this gap would be interesting. Our bounds are also loose by a factor of $\log\frac1\delta$ for all $p>0$, while we can get an optimal dependence on $\delta$ if we assume knowledge of $\OPT$ (Section <ref>) and sacrifice a factor of $\eps$. A natural question if one can achieve a simultaneously optimal dependence on $d$, $\eps$, and $\delta$, up to logarithmic factors, and without assumptions. Another gap to close is the query complexity of Huber regression, or more generally $M$-estimator regression, even for just the $d$ dependence: our upper bound is $\tilde O(d^{4-2\sqrt 2}\poly\log n)$ for constant $\eps$, while only a trivial lower bound of $\Omega(d)$ is known.

§ PRELIMINARIES

Throughout, we assume that $p$ is a fixed constant, and thus do not include constants depending on $p$ (e.g. $p$, $2^p$) in our big-O notation. We do include precise constants involving $p$ that appear in the exponents of other terms like $d$ and $1/\epsilon$. For simplicity, we also assume that our design matrix $\bfA\in \R^{n\times d}$ is overdetermined ($n > d$) and has full-rank (i.e., rank $d$). This is without loss of generality, as if $\bfA$ had rank $r < d$ we could replace it with an $n \times r$ basis for its column span and solve the regression problem with that basis as the design matrix.

Throughout, we will use the following fact, which follows from the convexity of $x^p$ for $x \ge 0$.

For any $p \ge 1$ and any $a,b \in \R$, $|a+b|^p \le 2^{p-1} (|a|^p + |b|^p) = O(|a|^p + |b|^p)$.

We denote the maximum of two real numbers $a,b$ by $a\lor b$ and the minimum as $a\land b$. We access the coordinates of vectors as $\bfv_i$ or $\bfv(i)$, which means the $i$th coordinate of $\bfv\in\mathbb R^n$.

§.§ Sensitivities and Lewis Weights

We use two central and closely related notions of matrix row importance for $\ell_p$ linear regression.

For any $\bfA \in \R^{n \times d}$ and $0 < p < \infty$ define the $\ell_p$ sensitivity of row $i$ for $i\in[n]$ as:

\begin{align*}
\bfs^p_i(\bfA) \coloneqq \max_{\bfx \in \R^d\setminus\{0\}} \frac{|[\bfA\bfx](i)|^p}{\norm{\bfA\bfx}_p^p}.
\end{align*}

$[\bfA\bfx](i)$ denotes the $i^\text{th}$ entry of the vector $\bfA\bfx$. The $\ell_p$ sensitivity measures how large an entry of a vector in the column span of $\bfA$ can be compared to the total $\ell_p$ norm of that vector. When $p = 2$, $\bfs^2_i(A)$ is equivalent to the statistical leverage score of row $i$.

A related, but not equivalent, generalization of the leverage scores, are the Lewis weights:

For any $\bfA \in \R^{n \times d}$ and $0<p<\infty$ the $\ell_p$ Lewis weights $\{\bfw^p_1(\bfA), \ldots, \bfw^p_n(\bfA) \}$ are the unique set of weights such that, if we let $\bfW \in \R^{n \times n}$ be the diagonal matrix with $\bfW(i,i) = \bfw_i^p(\bfA)$, then for all $i$,

\begin{align*}
\bfw_i^p(\bfA) = \bfs^2_i (\bfW^{1/2-1/p} \bfA).
\end{align*}

Here $\bfs^2_i (\bfW^{1/2-1/p} \bfA)$ is the $i^\text{th}$ leverage score of $\bfW^{1/2-1/p} \bfA$. Note that the leverage scores of a full-rank matrix always sum to $d$, so we have $\sum_{i=1}^n \bfw^p_i(\bfA) = d$.

Note that [17] only (explicitly) define $\ell_p$ Lewis weights for $p\geq 1$. However, the definition in fact makes sense for the range $p\in(0,1)$ as well, as shown by [63].

From Definition <ref> we can see that for $p = 2$, we have $\bfs^2_i(\bfA) = \bfw^2_i(\bfA)$, so the Lewis weights also correspond to the statistical leverage scores. While not equivalent for other values of $p$, the Lewis weights do upper bound the sensitivities, up to a polynomial in $d$ factor.

For $\bfA \in \R^{n \times d}$ and $0<p<\infty$,
$$\bfs_i^p(\bfA) \le d^{\max(0,p/2-1)} \cdot \bfw_i^p(\bfA).$$

The result is stated in Lemma 3.8 of [23] for $p\geq 1$. See Lemma <ref> for a proof of the case when $0 < p < 1$.

Lemma <ref> immediately yields a bound on the total sum of sensitivities.

For $\bfA \in \R^{n \times d}$ and $0<p<\infty$, we have $ \sum_{i=1}^n \bfs^p_i(\bfA) \le d^{\max(1,p/2)}$.

By Definition <ref>, the Lewis weights are the leverage scores of $\bfW^{1/2-1/p} \bfA$. It is well known that $\sum_{i=1}^n \bfw_i^p(\bfA) = \rank(\bfW^{1/2-1/p} \bfA) \le d$. The lemma then follows by applying Lemma <ref>.

Finally, we recall that Lewis weights can be efficiently approximated:

Let $\bfA\in\mathbb R^{n\times d}$. There is an algorithm which computes upper bounds $\tilde\bfw_i^p(\bfA) \geq \bfw_i^p(\bfA)$ to the $\ell_p$ Lewis weights of $\bfA$ such that

\[
\sum_{i=1}^n \tilde\bfw_i^p(\bfA) = O(d)
\]

with the following running times:

* for $p = 2$: $O(\nnz(\bfA) + d^\omega\log^2 d)$

* for $0 < p < 4$: $O(\nnz(\bfA)\log n + d^{\omega})$

* for $p \geq 2$: $O(\nnz(\bfA)\log n + d^{p/2 + O(1)})$

Here, $\omega \approx 2.37286$ is the current exponent of fast matrix multiplication.

§.§ Lewis Weight Sampling and lp Subspace Embedding

Using Lemma <ref> and a standard Bernstein bound, one can show that sampling $m$ rows of $\bfA$ for

\[
m = \tilde O(d^{\max(1,p/2)}/\epsilon^2)
\]

according to their sensitivities (and appropriately reweighting to keep the expectation correct) preserves $\norm{\bfA\bfx}_p$ to $(1 \pm \epsilon)$ relative error for any fixed $\bfx \in \R^d$ with high probability. Via an $\epsilon$-net argument, this can be shown to hold for all $\bfx$ simultaneously, if $\tilde O(d^{\max(1,p/2)+1}/\epsilon^2)$ rows are sampled, giving a method for $\ell_p$ subspace embedding. This argument holds when sampling by any set of upper bounds on the sensitivities that sum to $O(d^{\max(1,p/2)})$, including, e.g., the Lewis weights scaled up by a $d^{\max(0,p/2-1)}$ factor by Lemma <ref>. Note that this result gives a dimension bound that is independent of $n$.

Importantly however, the Lewis weights can be combined with more refined chaining arguments [6, 17] to shave a $d$ factor from the above bound, giving $\ell_p$ subspace embedding bounds into $\tilde O(d^{\max(1,p/2)}/\poly(\eps)))$ dimensions, which achieves a tight dependence on $d$, as shown by [47]. These bounds have been subsequently sharpened by [64, 65, 74, 63, 46], which use Lewis weights in a slightly different manner to achieve a bound of $\tilde O(d^{\max(1,p/2)}/\eps^2)$ dimensions. For $p\leq 2$, this type of analysis can in fact be turned into an analysis for the algorithm for sampling directly by the Lewis weights, via a reduction found in Lemma 7.4 of [17]. This proof crucially uses the monotonicity properties of Lewis weights for $p\leq 2$, and does not work for $p > 2$. To simplify the presentation and work with similar techniques simultaneously for all $0 < p < \infty$, we stick with analyzing the “recursive” form of Lewis weight sampling rather than the “one shot” form of Lewis weight sampling.

§.§.§ Chaining and the Iteration Procedure

We now describe the improved chaining arguments of [64, 65, 74, 63, 46], since we build on these to obtain our active regression algorithms. In all of these works, the algorithms proceed in multiple rounds of reducing the number of rows by a constant factor, rather than reducing the rows to $\tilde O(d^{\max(1,p/2)}/\eps^2)$ rows in one shot. Roughly, the idea is to sample each row $i\in[n]$ of $\bfA$ with probability $1/2$ and scale the result by $2$. This amounts to analyzing the quantity

\begin{align*}
\Lambda &\coloneqq \sup_{\norm*{\bfA\bfx}_p\leq 1} \abs*{\norm*{\bfS\bfA\bfx}_p^p - \norm*{\bfA\bfx}_p^p} \\
&= \sup_{\norm*{\bfA\bfx}_p \leq 1}\abs*{\sum_{i=1}^n (1+\sigma_i)\abs*{[\bfA\bfx](i)}^p - \sum_{i=1}^n \abs*{[\bfA\bfx](i)}^p} \\
&= \sup_{\norm*{\bfA\bfx}_p \leq 1}\abs*{\sum_{i=1}^n \sigma_i\abs*{[\bfA\bfx](i)}^p}
\end{align*}

where $\sigma_i\in\{\pm1\}$ are independent Rademacher variables and $\bfS\in\mathbb R^{n\times n}$ is a diagonal sampling matrix with $\bfS_{i,i} = (1+\sigma_i)$. The following is known:

Suppose that $\bfA\in\mathbb R^{n\times d}$ has Lewis weights bounded by $cd/n$ for some constant $c > 0$. Then, there exists a constant $C>0$ such that the following holds:

* if $2 \leq p < \infty$, then

\[
\E_{\bfsigma}[\Lambda] \leq \bracks*{C p^2 \frac{d^{p/2}}{n}(\log d)^2\log n}^{1/2}
\]

* if $1 < p \leq 2$, then

\[
\E_{\bfsigma}[\Lambda] \leq \bracks*{C\frac{d}{n}(\log d)^2 \max\braces*{\frac1{p-1}, \log n}}^{1/2}
\]

* if $0 < p < 1$, then

\[
\E_{\bfsigma}[\Lambda] \leq \bracks*{\frac{C}{p}\frac{d}{n}(\log d)^3}^{1/2}
\]

If $p$ is a constant, then we may concisely write the result of Theorem <ref> (giving up some log factors) as

\[
\E_{\bfsigma}[\Lambda] \leq \bracks*{C\frac{d^{\max(1,p/2)}}{n}(\log d)^2\log n}^{1/2}.
\]

The bounds of Theorem <ref> are proven using Dudley's entropy integral inequality, which comes with the following tail bound, which gives sub-Gaussian tails on the above quantity.

Let $(X_t)_{t\in T}$ be a Gaussian process with pseudo-metric $d_X(s,t)\coloneqq \norm*{X_s - X_t}_2$. Let $E(T, d_X, u)$ denote the minimal number of $d_X$-balls of radius $u$ required to cover $T$. Then, for every $u\geq 0$, we have that

\[
\Pr\braces*{\sup_{t\in T}X_t \geq C\bracks*{\int_0^\infty \sqrt{\log E(T, d_X, u)}~du + z\cdot \diam(T)}} \leq 2\exp(-z^2)
\]

This results in the following:

\begin{equation}\label{eq:lt-sz-tail-bound}
\Pr\braces*{\Lambda \geq C\bracks*{\frac{d^{\max(1,p/2)}}{n}}^{1/2}\bracks*{(\log d)\sqrt{\log n}+z}} \leq 2\exp(-z^2).
\end{equation}

The necessary diameter calculations are carried out in Lemma <ref> for $1 < p < \infty$ and in [63]for $0 < p < 1$.

Furthermore, the assumption that the Lewis weights are bounded by $O(d/n)$ can be enforced by a standard procedure of “splitting rows” (see, e.g., Remark 2.2 of [63]):

Let $\bfA\in\mathbb R^{n\times d}$ and let $0 < p < \infty$. Let $\tilde\bfw_i^p(\bfA) \geq \bfw_i^p(\bfA)$ be Lewis weight upper bounds such that

\[
\sum_{i=1}^n \tilde\bfw_i^p(\bfA) \leq C_1\cdot d.
\]

Let $C_2 > 0$. Then, there exists an $n'\times d$ matrix with $n \leq n' \leq (1+2C_1/C_2)n$ such that

\[
\norm*{\bfA\bfx}_p = \norm*{\bfA'\bfx}_p
\]

for all $\bfx\in\mathbb R^d$, and furthermore, $\bfw_i^p(\bfA') \leq C_2d/n$ for all $i\in[n']$.

Let $\bfA'$ be the $(n+k-1)\times d$ matrix obtained by replacing the first row $\bfa_1$ by $k$ copies of $\bfa_1 / k^{1/p}$. Let $\bfw_i$ denote the Lewis weights of $\bfA$. Then,

\[
\sum_{i=1}^n \bfw_i^{1-2/p}\bfa_i\bfa_i^\top = k\cdot (\bfw_1/k)^{1-2/p} \frac{\bfa_1}{k^{1/p}}\frac{\bfa_1^\top}{k^{1/p}} +  \sum_{i=2}^n \bfw_i^{1-2/p}\bfa_i\bfa_i^\top
\]

\[
\bracks*{\frac{\bfa_i^\top}{k^{1/p}}(\bfA^\top\bfW^{1-2/p}\bfA)^{-1}\frac{\bfa_i}{k^{1/p}}}^{p/2} = \frac{\bfw_i}{k}
\]

so $\bfw_1/k$ is the Lewis weight for all $k$ copies of $\bfa_1$ in $\bfA'$. Furthermore, $\norm*{\bfA\bfx}_p = \norm*{\bfA'\bfx}_p$ for every $\bfx\in\mathbb R^d$. Now suppose that for any row $i\in[n]$ with $\tilde\bfw_i^p(\bfA) \geq C_2 d/n$, we replace the row with $\ceil*{\tilde\bfw_i^p(\bfA)/(C_2 d/n)}$ scaled copies of $\bfa_i$. Then, we add at most

\begin{align*}
\sum_{i : \tilde\bfw_i^p(\bfA) \geq C_2d/n} \ceil*{\frac{\tilde\bfw_i^p(\bfA)}{C_2d/n}} &\leq \sum_{i : \tilde\bfw_i^p(\bfA) \geq C_2d/n} \frac{\tilde\bfw_i^p(\bfA)}{C_2d/n} + 1 \\
&\leq \frac{\sum_{i:\tilde\bfw_i^p(\bfA) \geq C_2d/n}\tilde\bfw_i^p(\bfA)}{C_2d/n} + \abs*{\braces*{i:\tilde\bfw_i^p(\bfA) \geq C_2d/n}} \\
&\leq 2\frac{C_1 d}{C_2 d/n} = \frac{2C_1}{C_2}n
\end{align*}

new rows, and furthermore, all Lewis weights of the new matrix $\bfA'$ are at most $C_2 d/n$.

Given approximations to the Lewis weights that sum to $C_1 d$ (using, e.g., Theorem <ref>), we may choose $C_2$ to be a large enough constant so that $2C_1/C_2 \leq 1/3$. Then, after splitting rows and sampling half of the rows, we only have $(4/3) \cdot (1/2) = 2/3$ of the rows remaining in expectation, and thus at most $3/4$ of the rows with probability at least $1 - \delta$ as long as the expected number of rows is at least $\Theta(\log(1/\delta))$, by Chernoff bounds. We condition on this event. Now let $n_i$ denote the number of rows remaining after $i$ rounds of splitting and sampling, until we have at most $m \coloneqq O(\eps^{-2}d^{\max(1,p/2)}[(\log d)^2\log n + \log(1/\delta)])$ rows remaining. Furthermore, let $\Lambda_i$ denote the value of $\Lambda$ after the $i$th round. Then, by applying Equation (<ref>) with $z$ set to $\sqrt{\log((\log n)/\delta)}$, we can union bound over the at most $O(\log n)$ iterations to conclude that with probability at least $1-\delta$, we simultaneously have

\[
\Lambda_i \leq C\bracks*{\frac{d^{\max(1,p/2)}}{n_i}\bracks*{(\log d)^2\log n+\log\frac1\delta}}^{1/2}
\]

for all iterations $i$. Then, the total distortion of the subspace embedding over all of the iterations is at most

\[
\prod_{i} (1+\Lambda_i) \leq \exp\parens*{\sum_i \Lambda_i} \leq \exp\parens*{O(\eps) \cdot \sum_{j=0}^\infty (\sqrt{3/4})^i} = \exp(O(\eps)) \leq 1 + O(\eps).
\]

This gives the following theorem:

Let $0 < p < \infty$. There is a randomized algorithm which constructs a sampling matrix $\bfS\in\mathbb R^{r\times n}$ such that with probability at least $1-\delta$, we have

\[
r = O\parens*{\frac{d^{\max(1,p/2)}}{\eps^2}\bracks*{(\log d)^2\log n + \log\frac1\delta}}
\]

and $\norm*{\bfS\bfA\bfx}_p = (1\pm \eps)\norm*{\bfA\bfx}_p$ for every $\bfx\in\mathbb R^d$.

Although Theorem <ref> has a $\log n$ dependence, this can be replaced with a $\log(d/\eps)$ dependence by first applying the subspace embedding result using the simpler union bound over a net, which has no dependence on $n$.

We give the pseudocode for the splitting and sampling procedure below in Algorithm <ref>.

$\ell_p$ sampling
input: Matrix $\bfA \in \R^{n \times d}$, constants $C_1, C_2 > 0$.

output: Diagonal sampling matrix $\bfS$.
[1]
Let $\tilde\bfw_i^p(\bfA)\geq\bfw_i^p(\bfA)$ be $\ell_p$ Lewis weight upper bounds summing to $C_1 d$ (Theorem <ref>).
Obtain $\bfA'$ by replacing each row $i\in[n]$ with $k = \ceil*{\tilde\bfw_i^p(\bfA)/(C_2d/n)}$ copies of $\bfa_i/k^{1/p}$.
Sample each row of $\bfA'$ with probability $1/2$ and scale by $2^{1/p}$.
$\bfS$, the matrix corresponding to the above procedure.

§ UPPER BOUND FOR LP REGRESSION

At a high level, our active regression algorithm follows the $\ell_p$ subspace embedding algorithm of Theorem <ref> by using Lewis weights to split rows and sample half of the rows.

§.§ Constant Factor Approximation

Analyzing the first stage is simple: as mentioned in Section <ref>, constant factor approximation algorithms for Problem <ref> can be obtained from existing subspace embedding results. Formally, we study Algorithm <ref>, which is a simple “sample-and-solve” approach to Problem <ref>. We analyze the method using the $\ell_p$ subspace embedding of Theorem <ref>, along with Markov's inequality.

Constant factor $\ell_p$ regression
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb \in \R^n$.

output: Approximate solution $\tilde \bfx \in \R^d$ to $\min_\bfx\|\bfA\bfx - \bfb\|_p$.
[1]
Let $\bfS \in \R^{m \times n}$ be an $1/2$-approximate $\ell_p$ subspace embedding for $\bfA$ (Theorem <ref>).
$\tilde \bfx$ with $\norm{\bfS\bfA\tilde \bfx - \bfS\bfb}_p \le (1+\gamma) \cdot \min_{\bfx \in \R^d} \norm{\bfS\bfA\bfx - \bfS\bfb}_p$ for $\gamma \geq 0$.

Running Algorithm <ref> only requires querying $m$ entries of $\bfb$ in order to construct the vector $\bfS\bfb$. Also note that in Line <ref> of the algorithm, we would have $\gamma = 0$ if an exact minimizer of the subsampled regression problem $\min_\bfx\norm{\bfS\bfA \bfx - \bfS\bfb}_p^p$ was obtained. To allow for the use of approximation algorithms in implementing Line <ref>, we state the method for a general $\gamma \geq 0$.

We first give an algorithm which works with constant probability, and then show how to boost the probability to $1-\delta$ for any $\delta \in (0,1)$ incuring an $O(\log(1/\delta))$ factor overhead in our sample complexity.

For $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $0 < p < \infty$, let $\OPT = \min_{\bfx \in \R^d} \norm{\bfA\bfx - \bfb}_p$. For any $\delta \in (0,1]$, if $\tilde\bfx$ is the output of Algorithm <ref>, then with probability at least $1-\delta$,

\[
\norm{\bfA\tilde \bfx - \bfb}_p \le 2^{2\max\{0,1/p-1\}+1+1/p}(3+\gamma)/\delta^{1/p} \cdot \OPT.
\]

When $\delta$ is constant (e.g., $\delta = 1/100$) and $(1+\gamma)$ is constant (e.g., $\gamma = 0$) then $\norm{\bfA\tilde\bfx - \bfb}_p \le C\cdot \OPT$ for constant $C$.

Let $\bfx^* = \argmin_{\bfx \in \R^d} \norm{\bfA\bfx - \bfb}_p$. By triangle inequality for $p\geq1$ or subadditivity and approximate triangle inequality (Fact <ref>) for $p\in(0,1)$,

\[
\norm{\bfA\tilde\bfx - \bfb}_p \le 2^{\max\{0,1/p - 1\}}\parens*{\norm{\bfA\bfx^* - \bfb}_p + \norm{\bfA\tilde\bfx - \bfA\bfx^*}_p} = 2^{\max\{0,1/p - 1\}}\parens*{\OPT + \norm{\bfA\tilde \bfx - \bfA \bfx^*}_p}.
\]

Applying the subspace embedding property of Theorem <ref> with $\epsilon =1/2$ and failure probability $\delta/2$, we conclude that, with probability at least $1-\delta/2$,

\begin{align*}
\norm{\bfA\tilde \bfx - \bfb}_p \le 2^{\max\{0,1/p - 1\}}\parens*{\OPT + 2 \norm{\bfS\bfA\tilde \bfx - \bfS\bfA \bfx^*}_p}.
\end{align*}

By similar reasoning, we have $\parens*{\norm{\bfS\bfA\tilde \bfx - \bfS\bfA \bfx^*}_p} \le 2^{\max\{0,1/p - 1\}}\parens*{\norm{\bfS\bfA\tilde \bfx - \bfS\bfb}_p + \norm{\bfS\bfA \bfx^* - \bfS\bfb}_p}$. We know that $\norm{\bfS\bfA\tilde \bfx - \bfS\bfb}_p \le (1+\gamma) \cdot \min_{\bfx \in \R^d} \norm{\bfS\bfA\bfx - \bfS\bfb}_p \le (1+\gamma) \cdot \norm{\bfS\bfA \bfx^* - \bfS\bfb}_p$, so we conclude that

\begin{align*}
\norm{\bfS\bfA\tilde \bfx - \bfS\bfA \bfx^*}_p \le 2^{\max\{0,1/p - 1\}}(2+\gamma)\norm{\bfS\bfA \bfx^* - \bfS\bfb}_p.
\end{align*}

Finally, by Markov's inequality, since $\E[\norm{\bfS\bfA \bfx^* - \bfS\bfb}_p^p] = \OPT^p$, with probability $\ge 1-\delta/2$, $\norm{\bfS\bfA \bfx^* - \bfS\bfb}_p^p \le \OPT^p/(\delta/2)$ and so $\norm{\bfS\bfA \bfx^* - \bfS\bfb}_p \le \OPT/(\delta/2)^{1/p}$.
Combining all these bounds we have that with probability $1-\delta$,

\begin{align*}
\norm{\bfA\tilde \bfx - \bfb}_p &\le 2^{\max\{0,1/p - 1\}}\parens*{\OPT + 2 \cdot 2^{\max\{0,1/p - 1\}}(2+\gamma) \cdot 2^{1/p} \OPT/\delta^{1/p}} \\
&\le 2^{2\max\{0,1/p-1\}+1+1/p}(3+\gamma)\OPT/\delta^{1/p}.\qedhere
\end{align*}

§.§.§ Probability Boosting for Constant Factor Approximation

We now show a boosting step for our constant factor approximation algorithm (Algorithm <ref>), described in Algorithm <ref>. If we repeat the constant factor approximation algorithm with success probability $99/100$ for a total of $O(\log(1/\delta))$ times, then via a standard Chernoff bound, with probability at least $1-\delta$, at least $9/10$ of the computed $\bfx_c$ will satisfy the guarantee of Theorem <ref> – i.e., that $\norm{\bfA\bfx_c-\bfb}_p = O(\OPT)$. Thus, we just need to identify one of these good solutions, which Algorithm <ref> does, deterministically, and without reading any entries of $\bfb$. The approach simply computes pairwise distances between solutions and returns any solution with a relatively low distance to at least $1/2$ of the other solutions. For later use, we state the result in terms of a general error measure $\norm*{\cdot}$ which satisfies an approximate triangle inequality (for example, $\norm*{\cdot}_p$ for $p\in(0,1)$ satisfies an approximate triangle inequality with constant $2^{1/p - 1}$ by Fact <ref>).

Constant factor $\norm*{\cdot}$ regression – Boosted Success Probability
input: $\ell$ candidate solutions $\bfx_1,\ldots,\bfx_\ell$ with at least $9/10\cdot \ell$ satisfying $\norm{\bfA\bfx_i-\bfb} \le \alpha \min_{\bfx} \norm{\bfA\bfx-\bfb}$.

output: Approximate solution $\tilde \bfx \in \R^d$ to $\min_{\bfx} \norm{\bfA\bfx-\bfb}$.
[1]

Let $\bfd \in \R^{\ell^2}$ contain all pairwise distances $\norm{\bfA\bfx_{i}-\bfA\bfx _{j}}$ (over ordered pairs $(i,j)$) sorted in increasing order. Let $\tau = \bfd(\floor*{\ell^2 \cdot 8/10})$ be the $80^{th}$ percentile distance.
Return any $\bfx_{i}$ such that $\norm{\bfA\bfx _{i}-\bfA\bfx _{j}} \le \tau$ for at least $1/2 \cdot \ell$ vectors $\bfx_{j}$.

Consider $\bfA \in \R^{n \times d}$, $\bfb \in \R^{n}$, and an error measure $\norm*{\cdot}$ which satisfies an approximate triangle inequality, that is, there exists a constant $\kappa\geq1$ such that $\norm*{\bfy_1+\bfy_2} \leq \kappa(\norm*{\bfy_1} + \norm*{\bfy_2})$ for any two vectors $\bfy_1,\bfy_2\in\R^n$. Let $\OPT = \min_{\bfx \in \R^d} \norm{\bfA\bfx-\bfb}$. Given a set of solution vectors $\bfx_1,\ldots,\bfx_\ell \in \R^d$ where $\norm{\bfA\bfx_i - \bfb} \le \alpha \cdot \OPT$ for at least $9/10 \cdot \ell $ of the vectors, Algorithm <ref> identifies $\bfx_i$ with $\norm{\bfA\bfx_i - \bfb} \le (\kappa\alpha + 2\kappa^3(\alpha+1)) \cdot \OPT$, without querying any entries of $\bfb$.

Let $\bfx^* = \argmin_{\bfx \in \R^d} \norm{\bfA\bfx-\bfb}$. Call $\bfx_i$ good if $\norm{\bfA\bfx_i - \bfb} \le \alpha \cdot \OPT$. By approximate triangle inequality , for any good $\bfx_i$,

\begin{align*}
\norm{\bfA\bfx_i-\bfA\bfx^*} \le \kappa(\norm{\bfA\bfx_i-\bfb}_p + \norm{\bfA\bfx^*-\bfb}_p) = \kappa(\alpha+1) \cdot \OPT.
\end{align*}

Thus, again via approximate triangle inequality, for any good $\bfx_i,\bfx_j$,

\[
\norm{\bfA\bfx_i-\bfA\bfx_j} \le \kappa(\kappa(\alpha+1) \cdot \OPT + \kappa(\alpha+1) \cdot \OPT) = 2\kappa^2(\alpha+1) \cdot \OPT.
\]

Thus, for the pairwise distance vector $\bfd \in \R^{\ell^2}$ computed in line 1 of Algorithm <ref>, at least $(9/10)^2 \cdot \ell^2 \ge 8/10 \cdot \ell^2$ of the distances will be upper bounded by $2\kappa^2(\alpha+1) \cdot \OPT$. Thus, the threshold $\tau$ computed in Line <ref>, which is the $80^{th}$ percentile of the distances, gives a lower bound $\tau \le 2\kappa^2(\alpha+1) \cdot \OPT$. In Line <ref>, we return any $\bfx_i$ with $\norm{\bfA\bfx_i - \bfA\bfx_j} \le \tau$ for at least $1/2 \cdot \ell$ vectors $\bfx_j$. First observe that at least one such $\bfx_i$ must exist. Otherwise, at most $1/2\cdot \ell^2$ of the pairwise distances would lie below $\tau$.

Additionally, observe that since at least $9/10 \cdot \ell$ of the $\bfx_i$ are good, if $\bfx_i$ is returned, it must have $\norm{\bfA\bfx_i - \bfA\bfx_j} \le \tau \le 2\kappa^2(\alpha+1) \cdot \OPT$ for at least one good $\bfx_j$. Since this good $\bfx_j$ has $\norm{\bfA\bfx_j - \bfb} \le \alpha \cdot \OPT$, by approximate triangle inequality, the returned $\bfx_i$ must then satisfy

\[
\norm{\bfA\bfx_i -\bfb} \le \kappa(\alpha + 2\kappa^2(\alpha+1)) = (\kappa\alpha + 2\kappa^3(\alpha+1))\cdot \OPT.
\]

§.§ Relative Error Approximation

We next show how to achieve a relative error $(1+\epsilon)$ solution to Problem <ref>. The relative error algorithm simply solves a Lewis weight sampled regression problem multiple times and outputs a boosted solution as done in Theorem <ref>. Pseudocode is given in Algorithm <ref> below.

Recursive relative error $\ell_p$ regression
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb \in \R^n$, number of measurements $m$.

output: Approximate solution $\tilde \bfx \in \R^d$ to $\min_\bfx\|\bfA\bfx - \bfb\|_p$.
[1]
$n \leq m$
Compute $\bar\bfx$ with $\norm{\bfA\bar\bfx - \bfb}_p \le (1+\gamma) \cdot \min_{\bfx \in \R^d} \norm{\bfA\bfx - \bfb}_p$ for $\gamma \geq 0$.
Let $\bfS \in \R^{m \times n}$ be generated via Algorithm <ref>.
Run Algorithm <ref> with inputs $\bfS\bfA$ and $\bfS\bfb$ to get $\tilde\bfx$.
$\tilde \bfx$

High probability relative error $\ell_p$ regression
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb \in \R^n$, number of measurements $m$.

output: Approximate solution $\tilde \bfx \in \R^d$ to $\min_\bfx\|\bfA\bfx - \bfb\|_p$.
[1]
Run Algorithm <ref> with inputs $\bfA$ and $\bfb$ independently $\ell$ times for $\ell = O(\log\frac1\delta)$.
Run Algorithm <ref> (and <ref>) to find $\bfx_c$ such that $\norm{\bfA\bfx_c-\bfb}_p \leq \alpha \min_\bfx\norm{\bfA\bfx-\bfb}_p$ for $\alpha > 0$.
Discard the $1/10$ fraction of trials with the largest $\norm*{\bfS(\bfb - \bfA\bfx_c)}_p^p$.
Run Algorithm <ref> on the remaining candidate solutions.

Let $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $0 < p < \infty$. Let

\[
m = O(1)\frac{d^{1\lor (p/2)}}{\eps^{2\lor p}}\bracks*{(\log d)^2\log n + \log\frac1\delta}\log\frac1\delta.
\]

Then, with probability at least $1-\delta$, Algorithm <ref> returns $\tilde\bfx\in\mathbb R^d$ such that

\[
\norm*{\bfA\tilde\bfx-\bfb}_p \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx-\bfb}_p
\]

and nonadaptively reads at most $m$ entries of $\bfb$.

Our main result, Theorem <ref>, follows directly from Theorem <ref>, up to a $\log n$ factor that is replaced by a $\log(d/\eps)$ factor by results in Section <ref>, as well as an optimization in the $\eps$ dependence for $p\in(1,2)$ as done in Section <ref>.
The proof of Theorem <ref> proceeds in 2 steps. We first show in Lemma <ref> that any entries of the residual $\bfz = \bfb-\bfA\bfx^*$ whose contribution to $\norm{\bfz}_p^p$ is significantly larger than the corresponding Lewis weight in $\bfA$ can be effectively ignored, since no $\bfx$ achieving small error can accurately fit such entries. We then show in Lemma <ref> that if we ignore these entries, sampling by the $\ell_p$ Lewis weights approximately preserves $\norm{\bfA\bfx-\bfz}_p$ for all $\bfx \in \R^d$. The result follows by combining Lemmas <ref> and <ref>.

Let $\bfA\in\mathbb R^{n\times d}$ and $\bfb\in\mathbb R^n$, and let $0 < p < \infty$. Let $\bfz = \bfb - \bfA\bfx^*\in \R^n$. Let

\[
\mathcal B \coloneqq \braces*{i\in[n] : \frac{|\bfz(i)|^p}{\OPT^p} \ge \frac{d^{\max(0,p/2-1)} \cdot \bfw_{i}^p(\bfA)}{\epsilon^p}}
\]

Let $\bar \bfz \in \R^n$ be equal to $\bfz$ but with all entries in $\mathcal{B}$ set to $0$. Then for all $\bfx \in \R^d$ with $\norm{\bfA\bfx}_p = O(\OPT)$,

$$\left | \norm{\bfA\bfx-\bfz}_p^p - \norm{\bfA\bfx-\bar \bfz}_p^p - \norm*{\bfz-\bar\bfz}_p^p \right | = O( \epsilon) \cdot \OPT^p.$$

For any $\bfx \in \R^d$ and any $i \in \mathcal{B}$, we have that

\begin{align*}
|[\bfA\bfx](i)|^p &\le \norm{\bfA\bfx}_p^p \cdot \bfs_i^p(\bfA) \\
&\le \norm{\bfA\bfx}_p^p \cdot d^{\max(0,p/2-1)} \cdot \bfw^p_i(\bfA) && \text{Lemma \ref{lem:lwBound}} \\
&\le \norm{\bfA\bfx}_p^p \cdot \epsilon^p \cdot \frac{|\bfz(i)|^p}{\OPT^p} = O(\epsilon^p) \cdot |\bfz(i)|^p,
\end{align*}

where we use the assumption that $\norm{\bfA\bfx}_p^p = O(\OPT^p)$.
From the above, we have that

\begin{align*}
|[\bfA\bfx](i) - \bfz(i)|^p - |[\bfA\bfx](i) - \bar \bfz(i)|^p &= |[\bfA\bfx](i) - \bfz(i)|^p - |[\bfA\bfx](i)|^p = (1 \pm O(\epsilon)) \cdot |\bfz(i)|^p
\end{align*}

for $i\in\mathcal B$.
Since $\norm{\bfz-\bar \bfz}_p^p \le \norm{\bfz}_p^p = \OPT^p$,

\begin{align*}
\left | \norm{\bfA\bfx-\bfz}_p^p - \norm{\bfA\bfx-\bar \bfz}_p^p -\norm{\bfz}_p^p \right | = O( \epsilon) \cdot \norm{\bfz-\bar \bfz}_p^p = O(\epsilon) \cdot \OPT^p.
\end{align*}

Consider the setting of Lemma <ref>. Let $\bfS \in \R^{m \times n}$ be formed as in Line <ref> of Algorithm <ref>. With probability at least $99/100$, $\norm{\bfS\bfz}_p^p = O(\OPT^p)$ and further, for all $\bfx \in \R^d$ with $\norm{\bfA\bfx}_p = O(\OPT)$,

\[
\left | \norm{\bfS\bfA\bfx-\bfS\bfz}_p^p - \norm{\bfS\bfA\bfx-\bfS\bar \bfz}_p^p - \norm{\bfS(\bfz-\bar \bfz)}_p^p \right | = O(\eps) \cdot \OPT^p.
\]

The proof follows that of Lemma <ref>. We have, for all $\bfx \in \R^d$ and any $i \in \mathcal{B}$ sampled by $\bfS$,

\begin{align*}
|[\bfS\bfA\bfx](i) - \bfS\bfz(i)|^p - |[\bfS\bfA\bfx](i) - \bfS\bar \bfz(i)|^p &= |[\bfS\bfA\bfx](i) - \bfS\bfz(i)|^p - |[\bfS\bfA\bfx](i)|^p \\
&= (1 \pm O(\epsilon)) \cdot |\bfS\bfz(i)|^p,
\end{align*}

where we have used that $\bfS$ is a subspace embedding so $\norm*{\bfS\bfA\bfx}_p = O(\norm*{\bfA\bfx}_p) = O(\OPT)$. Now, $\E[\norm{\bfS \bfz}_p^p] = \norm{\bfz}_p^p = O(\OPT^p)$. Then by Markov's inequality, with probability at least $99/100$, $\norm{\bfS(\bfz-\bar \bfz)}_p^p \le \norm{\bfS \bfz}_p^p = O(\OPT^p)$. We then have

\begin{align*}
\left | \norm{\bfS\bfA\bfx-\bfS\bfz}_p^p - \norm{\bfS\bfA\bfx-\bfS\bar \bfz}_p^p - \norm{\bfS(\bfz-\bar \bfz)}_p^p\right | = O(\epsilon) \cdot \norm{\bfS(\bfz-\bar \bfz)}_p^p = O(\eps)\cdot \OPT^p.
\end{align*}

While Lemma <ref> only succeeds with constant probability, we may boost it by repeating the procedure $\ell = O(\log\frac1\delta)$ times, sorting the trials by $\norm*{\bfS(\bfb-\bfA\bfx_c)}_p^p$ for a constant factor solution $\bfx_c$, and then discarding the top $1/10$ of trials (see Line <ref>). By Chernoff bounds, we will then select a trial satisfying

\[
\norm*{\bfS\bfz}_p \leq \norm*{\bfS\bfA(\bfx^*-\bfx_c)}_p + \norm*{\bfS(\bfb-\bfA\bfx_c)}_p = O(\OPT)
\]

with probability at least $1-\delta$.

We now state the following lemma, which shows that Lewis weight sampling as in Line <ref> of Algorithm <ref> will approximately preserve the $\ell_p$ norms of $\norm*{\bfA\bfx-\bar\bfz}_p$, provided that $\norm*{\bfA\bfx}_p = O(\OPT)$.

Consider the setting of Lemma <ref>. Let $\bfS \in \R^{m \times n}$ be formed as in Line <ref> of Algorithm <ref>. Let

\[
\Lambda \coloneqq \sup_{\norm*{\bfA\bfx}_p \leq O(\OPT)}\abs*{\norm*{\bfS\bfA\bfx-\bfS\bar\bfz}_p^p - \norm*{\bfA\bfx-\bar\bfz}_p^p}.
\]

\[
\Pr\braces*{\Lambda \geq \bracks*{C\frac1n \frac{d^{1\lor (p/2)}}{\eps^{0\lor(p-2)}}}^{1/2} \bracks*{(\log d)\sqrt{\log n} + z} \OPT^p} \leq 2\exp(-z^2)
\]

We defer the proof of Lemma <ref> to Section <ref> and first show how to use it to prove Theorem <ref>.

Let $\bfS$ be the sketching matrix generated as a result of the $O(\log n)$ rounds of recursion. As in the proof of Theorem <ref>, the algorithm makes at most $O(\log n)$ levels of recursion with probability at least $1-\delta$, and generates a $\bfS$ such that

\[
\Pr\braces*{\sup_{\norm*{\bfA\bfx}_p \leq O(\OPT)}\abs*{\norm*{\bfS\bfA\bfx-\bfS\bar\bfz}_p^p - \norm*{\bfA\bfx-\bar\bfz}_p^p} \geq \eps \OPT^p} \leq \frac{\delta}{\log\frac1\delta}
\]

Consider $\ell = O(\log\frac1\delta)$ independent repetitions of Algorithm <ref>. By a union bound over all $\ell$ repetitions, the above guarantee holds simultaneously for all repetitions. Note also that with probability at least $1-\delta$, we only keep trials such that $\norm*{\bfS\bfz}_p = O(\OPT)$ in Line <ref>. Furthermore, with probability at least $1-\delta$, the boosting algorithm of Theorem <ref> selects a run outputting a constant factor approximation. Let $\bar\bfx$ be the solution that is output by the algorithm. Since $\bar\bfx$ is a constant factor approximation, we have the guarantee that

\[
\norm*{\bfA\bar\bfx}_p \leq \norm*{\bfA\bar\bfx-\bfz}_p + \norm*{\bfz}_p = O(\OPT).
\]

Similarly, $\norm*{\bfA\bfx^*}_p = O(\OPT)$, and thus $\norm*{\bfA(\bar\bfx-\bfx^*)}_p = O(\OPT)$. Then,

\begin{align*}
\norm*{\bfA\bar\bfx-\bfb}_p^p &= \norm*{\bfA(\bar\bfx-\bfx^*)-\bfz}_p^p \\
&\leq \norm*{\bfA(\bar\bfx-\bfx^*)-\bar\bfz}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:clip}} \\
&\leq \norm*{\bfS\bfA(\bar\bfx-\bfx^*)-\bfS\bar\bfz}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:additive2}} \\
&\leq \norm*{\bfS\bfA(\bar\bfx-\bfx^*)-\bfS\bfz}_p^p - \norm*{\bfS(\bfz-\bar\bfz)}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:clip2}} \\
&\leq \norm*{\bfS\bfA\bfx-\bfS\bfb}_p^p - \norm*{\bfS(\bfz-\bar\bfz)}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p \\
&\leq \norm*{\bfS\bfA\bfx^*-\bfS\bfb}_p^p - \norm*{\bfS(\bfz-\bar\bfz)}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{near optimality of $\bar\bfx$} \\
&\leq \norm*{\bfS\bfA\bfx^*-\bfS\bar\bfz}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:clip2}} \\
&\leq \norm*{\bfA\bfx^*-\bar\bfz}_p^p + \norm*{\bar\bfz - \bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:additive2}} \\
&\leq \norm*{\bfA\bfx^*-\bfz}_p^p + O(\eps)\cdot\OPT^p && \text{Lemma \ref{lem:clip}} \qedhere
\end{align*}

§.§ Proof of Lemma <ref>

We return to proving Lemma <ref>, which shows that Lewis weight sampling using the Lewis weights of $\bfA$ can preserve the norm of $\bfA\bfx-\bar\bfz$. To do this, we will modify the proofs of the Lewis weight sampling subspace embeddings of [46, 63]. We return to the framework of Section <ref> which considers sampling each row with probability $1/2$ under the assumption that $\bfA$ has $\ell_p$ Lewis weights bounded by $O(d/n)$. However, this time, we must preserve the norm of $\bfA\bfx - \bar\bfz$ rather than $\bfA\bfx$. We thus consider the quantity

\[
\Lambda \coloneqq \sup_{\bfy\in T}\abs*{\sum_{i=1}^n \sigma_i\abs*{\bfy(i)}^p} = \sup_{\bfy\in T}\abs*{\sum_{i=1}^n \bfw_i \sigma_i\abs*{[\bfW^{-1/p}\bfy](i)}^p}
\]

where $T \coloneqq \braces*{\bfA\bfx-\bar\bfz : \norm*{\bfA\bfx}_p \leq O(\OPT)}$ and $\bfw_i$ are the $\ell_p$ Lewis weights of $\bfA$ which are bounded by $O(d/n)$. We now handle this by considering three different parameter regimes: $p\in (2, \infty)$, $p\in(1,2)$, and $p\in(0,1)$. Note that we exclude $p = 1$ and $p = 2$, as they have already been handled by previous work [18, 53, 11]. Throughout this section, we assume that $\OPT = O(1)$ by scaling.

We first adapt the work of [46] to handle the cases of $p\in (1,2)$ and $p\in (2,\infty)$. The following is an adaptation of Lemma 15.17 of [46], which is essentially the statement that Lewis weights of $\bfA$ uniformly bound the sensitivities over the set $T$, which is in turn a combination of Lemma <ref> and the definition of $\mathcal B$ in Lemma <ref>.

Let $p\geq 1$ and $\bfy\in T$. Then,

\[
\norm*{\bfW^{-1/p}\bfy}_\infty \leq O\parens*{\frac{d^{0\lor(1/2-1/p)}}{\eps}}\OPT
\]

Write $\bfy = \bfA\bfx-\bar\bfz$ for $\bfx\in\mathbb R^d$. For any coordinate $i\in[n]$, we have

\begin{align*}
\abs*{[\bfW^{-1/p}\bfy](i)}^p &= \bfw_i^{-1}\abs*{[\bfA\bfx-\bar\bfz](i)}^p \\
&\leq O(2^p)\bfw_i^{-1}(\abs*{[\bfA\bfx](i)}^p+\abs*{[\bar\bfz](i)}^p) \\
&\leq O(2^p)\bfw_i^{-1}\cdot \bfw_i d^{0\lor(p/2-1)}\bracks*{\norm*{\bfA\bfx}_p^p+\frac1{\eps^p}\norm*{\bar\bfz}_p^p} \\
&= O(2^p) \frac{d^{0\lor(p/2-1)}}{\eps^p}\OPT^p
\end{align*}

Taking $p$th roots yields the desired result.

We now separately handle indices $i\in[n]$ such that the Lewis weights $\bfw_i$ are less than some $\poly(\eps/n)$.

Define the set $J\subseteq[n]$ by

\[
J \coloneqq \braces*{i\in[n] : \bfw_i \geq \frac{\eps^p d}{n^2}}.
\]

For indices not contained in $J$, we use the following simple bound:

We have with probability $1$ that

\[
\sup_{\bfy\in T}\abs*{\sum_{i\notin J} \bfw_i \sigma_i \abs*{[\bfW^{-1/p}\bfy](i)}^p} \leq O\parens*{\frac{d^{1\lor (p/2)}}{n}}\OPT^p
\]

For any $\bfy\in T$, we have that

\begin{align*}
\abs*{\sum_{i\notin J} \bfw_i \sigma_i \abs*{[\bfW^{-1/p}\bfy](i)}^p} &\leq \sum_{i\notin J}\bfw_i\norm*{\bfW^{-1/p}\bfy}_\infty^p && \text{triangle inequality} \\
&\leq \sum_{i\notin J}\bfw_i O\parens*{\frac{d^{0\lor(p/2-1)}}{\eps^p}}\OPT^p && \text{Lemma \ref{lem:sensitivity-bound-T}} \\
&\leq \sum_{i\notin J}\frac{\eps^p d}{n^2} \cdot  O\parens*{\frac{d^{0\lor(p/2-1)}}{\eps^p}}\OPT^p && \text{$i\notin J$} \\
&\leq O\parens*{\frac{d^{1\lor(p/2)}}{n}}\OPT^p.\qedhere
\end{align*}

For the indices contained in $J$, we use a comparison theorem for Rademacher and Gaussian averages, and then use Dudley's entropy integral (Theorem <ref>) to bound the corresponding Gaussian process.

Let $\{X_t\}_{t\in T}$ be a Gaussian process. Define the pseudo-metric on $T$ via

\[
d_X(s,t) \coloneqq \norm*{X_s - X_t}_2.
\]

Let $E(T, d_X, u)$ denote the minimal number of $d_X$-balls of radius $u$ required to cover $T$. Then,

\[
\E\sup_{t\in T} X_t \leq 24 \int_0^\infty \sqrt{\log E(T, B_{X}, u)} ~du .
\]

We first use that $\bfw_i \leq O(d/n)$ to bound

\begin{equation}\label{eq:J-bound}
\E\sup_{\bfy\in T}\abs*{\sum_{i\in J} \bfw_i \sigma_i \abs*{[\bfW^{-1/p}\bfy](i)}^p} \leq \frac{O(1)}{\sqrt{n}}\cdot \E\sup_{\bfy\in T}\abs*{\sum_{i\in J} (d\bfw_i)^{1/2} \sigma_i \abs*{[\bfW^{-1/p}\bfy](i)}^p}
\end{equation}

We now study the latter quantity, with Rademacher variables $\sigma_i$ replaced by standard Gaussian variables $g_i$.

Let $\{X_\bfy\}_{\bfy\in T}$ be the Gaussian process defined by

\[
X_\bfy \coloneqq \sum_{i\in J}(d\bfw_i)^{1/2}g_i\abs*{[\bfW^{-1/p}\bfy](i)}^p.
\]

Note that the pseudo-metric $d_X$ associated with this Gaussian process (see Theorem <ref>) is

\begin{align*}
d_X(\bfy,\bfy') &= \norm*{X_\bfy - X_{\bfy'}}_2 \\
&= \sqrt{\E\bracks*{\sum_{i\in J}(d\bfw_i)^{1/2}g_i\abs*{[\bfW^{-1/p}\bfy](i)}^p - \sum_{i\in J}(d\bfw_i)^{1/2}g_i\abs*{[\bfW^{-1/p}\bfy'](i)}^p}^2} \\
&= \sqrt{\E\bracks*{\sum_{i\in J}(d\bfw_i)^{1/2}g_i(\abs*{[\bfW^{-1/p}\bfy](i)}^p - \abs*{[\bfW^{-1/p}\bfy'](i)}^p)}^2} \\
&= O(1)\sqrt{\sum_{i\in J}d\bfw_i(\abs*{[\bfW^{-1/p}\bfy](i)}^p - \abs*{[\bfW^{-1/p}\bfy'](i)}^p)^2}
\end{align*}

For tail bounds, we can use the same Rademacher contraction principle <cit.> and Rademacher average comparison <cit.> theorems along with a standard comparison lemma due to Panchenko <cit.>.

§.§.§ Entropy Bounds

To apply Dudley's entropy integral theorem, we now require estimates on the metric entropy $E(T, d_X, u)$. For easier entropy calculations, we bound $d_X$ by a more convenient metric, adapting Equations (15.18) and (15.19) in [46] to encorporate the change in Lemma <ref>.

For $\bfy,\bfy'\in T$, we have the following:

* if $2\leq p < \infty$, then

\[
d_X(\bfy,\bfy') \leq O(p\sqrt d) O\parens*{\frac{d^{1/2-1/p}}{\eps}}^{p/2-1} \norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty
\]

* if $1 < p < 2$, then

\[
d_X(\bfy,\bfy') \leq O(\sqrt d)\norm*{(\bfW^{-1/p}\bfy - \bfW^{-1/p}\bfy')\vert_J}_\infty^{p/2}
\]

We first handle $p\geq 2$. For $a,b\geq 0$, we have by convexity for $p > 1$ that

\[
a^p - b^p \leq p(a^{p-1}+b^{p-1})\abs*{a-b}.
\]

\begin{align*}
d_X(\bfy,\bfy')^2 &\leq O(d)\sum_{i\in J}\bfw_i\parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^p - \abs*{[\bfW^{-1/p}\bfy'](i)}^p}^2 \\
&\leq O(p^2d)\norm*{(\bfW^{-1/p}\bfy - \bfW^{-1/p}\bfy')\vert_J}_\infty^2\sum_{i\in J}\bfw_i\max\braces*{\abs*{[\bfW^{-1/p}\bfy](i)}, \abs*{[\bfW^{-1/p}\bfy'](i)}}^{2(p-1)}.
\end{align*}

We then use Lemma <ref> to bound

\begin{align*}
&\sum_{i\in J}\bfw_i\max\braces*{\abs*{[\bfW^{-1/p}\bfy](i)}, \abs*{[\bfW^{-1/p}\bfy'](i)}}^{2(p-1)} \\
\leq ~ &\max\braces*{\norm*{\bfW^{-1/p}\bfy}_\infty, \norm*{\bfW^{-1/p}\bfy'}_\infty}^{p-2} \sum_{i\in J}\bfw_i \parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^p + \abs*{[\bfW^{-1/p}\bfy'](i)}^p} \\
\leq ~ &O\parens*{\frac{d^{1/2-1/p}}{\eps}\OPT}^{p-2}(\norm*{\bfy}_p^p + \norm*{\bfy'}_p^p)\leq O\parens*{\frac{d^{1/2-1/p}}{\eps}}^{p-2}
\end{align*}

so combining the bounds and taking square roots gives the result.

For $1 < p < 2$, we have using the subadditivity of $|\cdot|^{p/2}$ that

\begin{align*}
d_X(\bfy,\bfy')^2 &\leq O(d)\sum_{i\in J}\bfw_i\parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^p - \abs*{[\bfW^{-1/p}\bfy'](i)}^p}^2 \\
&= O(d)\sum_{i\in J}\bfw_i\parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^{p/2} - \abs*{[\bfW^{-1/p}\bfy'](i)}^{p/2}}^2 \\
&\hspace{6.4em}\parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^{p/2} + \abs*{[\bfW^{-1/p}\bfy'](i)}^{p/2}}^2 \\
&\leq O(d)\sum_{i\in J}\bfw_i\abs*{[\bfW^{-1/p}\bfy - \bfW^{-1/p}\bfy'](i)}^p \parens*{\abs*{[\bfW^{-1/p}\bfy](i)}^{p} + \abs*{[\bfW^{-1/p}\bfy'](i)}^{p}} \\
&\leq O(d)\norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty^p (\norm*{\bfy}_p^p + \norm*{\bfy'}_p^p) \\
&\leq O(d)\norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty^p.\qedhere
\end{align*}

Because $\bfw_i$ is bounded below for $i\in J$, we may further bound

\begin{equation}\label{eq:dx-bound-wq}
\norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty \leq O(1)\norm*{\bfW^{-1/p}(\bfy-\bfy')}_{\bfw,q}
\end{equation}

for $q = O(\log n)$, where

\[
\norm*{\bfv}_{\bfw,q} \coloneqq \bracks*{\sum_{i=1}^n \frac{\bfw_i}{d}\abs*{\bfv(i)}^q}^{1/q}.
\]

Indeed, if $i\in J$, then

\[
\frac1{\poly(n)}\sum_{i\in J}\abs*{\bfv(i)}^q \leq \sum_{i=1}^n\frac{\bfw_i}{d}\abs*{\bfv(i)}^q = \norm*{\bfv}_{\bfw,q}^q
\]

so $\norm*{\bfv\vert_J}_\infty \leq O(1)\norm*{\bfv}_{\bfw,q}$.

Now by the results of Lemma <ref> and Equation (<ref>), we can replace $d_X$ balls by $\norm*{\cdot}_{\bfw,q}$ balls of an appropriate radius. Furthermore, we may note that the set $T$ is simply a translation of $\ell_p$ ball of radius $O(\OPT) = O(1)$ by $\bar\bfz$, so it suffices to cover the unit $\ell_p$ ball $B_p$ by $\norm*{\cdot}_{\bfw,q}$ balls. These are exactly the types of entropy estimates given by [6, 46]. Indeed, we have the following result:

Let $1 < p < \infty$ and $q = O(\log n)$. Let $E$ be the subspace spanned by the columns of $\bfW^{-1/p}\bfA$ and let $B_{\bfw,p}\subseteq E$ denote the unit $\norm*{\cdot}_{\bfw,p}$-ball in this subspace. Then, for some universal constant $C>0$, the following holds:

* if $2 \leq p < \infty$, then

\[
\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u) \leq C\frac{d\log n}{u^2}.
\]

* if $1 < p \leq 2$, then

\[
\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u) \leq C\frac{d\log n}{u^p}.
\]

For $u$ tending to $0$, a standard volume argument gives a better bound:

Consider the setting of Lemma <ref>. Then, for some universal constant $C>0$ we have that

\[
\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u) \leq C d\log \frac{d}{u}.
\]

Diameter Bounds.

Finally, we give diameter bounds for the application of the Dudley tail bound of Theorem <ref>.

The $d_X$-diameter of $T$ is bounded by $[O(d^{p/2}/\eps^{p-2})]^{1/2}$ for $2 \leq p < \infty$ and $O(\sqrt d)$ for $1 < p < 2$.

Let $\bfy = \bfA\bfx-\bar\bfz$ and $\bfy' = \bfA\bfx'-\bar\bfz$ be two points in $T$. Then,

\begin{align*}
\norm*{[\bfW^{-1/p}(\bfy-\bfy')]\vert_J}_\infty^p &= \norm*{[\bfW^{-1/p}\bfA(\bfx-\bfx')]\vert_J}_\infty^p \\
&\leq \norm*{\bfW^{-1/p}\bfA(\bfx-\bfx')}_\infty^p \\
&\leq \max_{i=1}^n O(d^{0\lor (p/2-1)})\bfw_i\cdot \bfw_i^{-1}\norm*{\bfA(\bfx-\bfx')}_p^p \\
&= O(d^{0\lor (p/2-1)})\norm*{\bfy-\bfy'}_p^p.
\end{align*}

It follows from the above calculation and Lemma <ref> that the $d_X$-diameter of the set $T$ is bounded by

\begin{align*}
\sup_{\bfy,\bfy'\in T}d_X(\bfy,\bfy') &\leq \sup_{\bfy,\bfy'\in T}O(p\sqrt d) O\parens*{\frac{d^{1/2-1/p}}{\eps}}^{p/2-1} \norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty \leq \bracks*{O(1)\frac{d^{p/2}}{\eps^{p-2}}}^{1/2}
\end{align*}

for $2\leq p < \infty$ and

\begin{align*}
\sup_{\bfy,\bfy'\in T}d_X(\bfy,\bfy') &\leq \sup_{\bfy,\bfy'\in T}O(\sqrt d) \norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_J}_\infty^{p/2} \leq O(\sqrt d)
\end{align*}

for $1 < p < 2$.

§.§.§ Proof for 2 < p < inf

We will calculate Dudley's entropy integral. As previously noted, we may replace $T$ by the set

\[
B_p \coloneqq \braces*{\bfA\bfx : \norm*{\bfA\bfx}_p \leq \OPT} = \braces*{\bfA\bfx : \norm*{\bfW^{-1/p}\bfA\bfx}_{\bfw,p} \leq \frac{\OPT}{d^{1/p}}} = \frac{\Theta(1)}{d^{1/p}}\bfW^{1/p}(B_{\bfw,p})
\]

since translations do not change cover numbers. Furthermore, we have by Lemma <ref> and Equation (<ref>) that for

\[
\alpha = O(p\sqrt d) O\parens*{\frac{d^{1/2-1/p}}{\eps}}^{p/2-1},
\]

we have that

\begin{align*}
\log E(B_p, d_X, t) &\leq \log E(\Theta(d^{-1/p})B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, t/\alpha) \\
&= \log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, \Theta(d^{1/p}t/\alpha))
\end{align*}

since linear transformations do not change cover numbers. Then, by Theorem <ref>, we have that

\begin{align*}
&\E\sup_{\bfy\in T}\abs*{(d\bfw_i)^{1/2} g_i \abs*{[\bfW^{-1/p}\bfy](i)}^p} \\
\leq ~&24\int_0^\infty \sqrt{\log E(T, d_X, u)}~du \\
\leq ~&\frac{O(\alpha)}{d^{1/p}}\int_0^\infty \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u)}~du \\
= ~&\frac{O(\alpha)}{d^{1/p}}\bracks*{\int_0^1 \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u)}~du + \int_1^\infty \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u)}~du} \\
\leq ~&\frac{O(\alpha)}{d^{1/p}}\bracks*{\int_0^1 \sqrt{d\log\frac{d}{u}}~du + \int_1^{\poly(d)} \frac{\sqrt{d\log n}}{u}~du} \tag{Lemmas \ref{lem:lt2011-entropy-estimate} and \ref{lem:simple-entropy-estimate}} \\
\leq ~&O(1)\frac{(d^{p/2})^{1/2}(\log d)\sqrt{\log n}}{\eps^{p/2-1}}.
\end{align*}

Combining the above bound with Equation (<ref>) and Lemma <ref> yields the result that

\[
\E\Lambda \leq \bracks*{O(1)\frac1n \frac{d^{p/2}(\log d)^2\log n}{\eps^{p-2}}}^{1/2}.
\]

Furthermore, applying Theorem <ref> and Lemma <ref> yields the result that

\[
\Pr\braces*{\Lambda \geq C\frac1{\sqrt n}\bracks*{\frac{d^{p/2}}{\eps^{p-2}}}^{1/2}\bracks*{(\log d)\sqrt{\log n}+z}} \leq 2\exp(-z^2).\qedhere
\]

§.§.§ Proof for 1 < p < 2

Our proof is very similar to that of $2 < p < \infty$. By combining the reasoning of the proof for $2 < p < \infty$ with Lemma <ref> and Equation (<ref>), we obtain that

\begin{align*}
\log E(T, d_X, t) &\leq \log E(\Theta(d^{-1/p})B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, (t/\sqrt d)^{2/p}) \\
&= \log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, \Theta(d^{1/p})(t/\sqrt d)^{2/p}) \\
&= \log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, \Theta(t^{2/p})).
\end{align*}

Then by Theorem <ref>, we have that

\begin{align*}
&\E\sup_{\bfy\in T}\abs*{(d\bfw_i)^{1/2} g_i \abs*{[\bfW^{-1/p}\bfy](i)}^p} \\
\leq ~&24\int_0^\infty \sqrt{\log E(T, d_X, u)}~du \\
\leq ~&O(1)\int_0^\infty \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u^{2/p})}~du \\
= ~&O(1)\bracks*{\int_0^1 \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u^{2/p})}~du + \int_1^\infty \sqrt{\log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, u^{2/p})}~du} \\
\leq ~&O(1)\bracks*{\int_0^1 \sqrt{d\log\frac{d}{u}}~du + \int_1^{\poly(d)} \frac{\sqrt{d\log n}}{u}~du} \tag{Lemmas \ref{lem:lt2011-entropy-estimate} and \ref{lem:simple-entropy-estimate}} \\
\leq ~&O(1)\sqrt d(\log d)\sqrt{\log n}.
\end{align*}

Combining the above bound with Equation (<ref>) and Lemma <ref> yields the result that

\[
\E\Lambda \leq \bracks*{O(1)\frac1n d(\log d)^2\log n}^{1/2}.
\]

Furthermore, applying Theorem <ref> and Lemma <ref> yields the result that

\[
\Pr\braces*{\Lambda \geq C\frac1{\sqrt n}\sqrt d\bracks*{(\log d)\sqrt{\log n}+z}} \leq 2\exp(-z^2).\qedhere
\]

§.§.§ Proof for 0 < p < 1

Note that in the proof of $1 < p < 2$, essentially the only thing that changed between the active setting and the proof of the subspace embedding was the translation by a vector $\bar\bfz$ with norm at most $O(\OPT) = O(1)$, which did not affect the entropy integral bound. In the case of $0 < p < 1$, the proof of [63] is structured similarly and can be easily checked to have essentially the same straightforward modifications.

§.§ Sharp Dependence on eps for 1 < p < 2

We will now optimize our dependence on $\eps$ for $1 < p < 2$.

§.§.§ Closeness of Near-Optimal Solutions

We first show that near-optimal solutions $\bfx\in\mathbb R^d$ are close to the optimal solution $\bfx^*$, using strong convexity with respect to the $\norm*{\cdot}_p$ norm.

Let $1 < p \leq 2$. Let $\bfA\in\mathbb R^{n\times d}$ and $\bfb\in\mathbb R^n$. Then, for any $\bfx\in\mathbb R^d$ such that

\[
\norm*{\bfA\bfx-\bfb}_p \leq (1+\gamma) \OPT
\]

with $\gamma\in(0,1)$, we have that

\[
\norm*{\bfA\bfx-\bfA\bfx^*}_p \leq O(\sqrt\gamma) \OPT,
\]

where $\bfx^* \coloneqq \arg\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx-\bfb}_p$

For $p = 2$, this follows from the Pythagorean theorem. Otherwise, note first that the KKT conditions require that

\[
\angle*{(\bfA\bfx^*-\bfb)^{p-1},\bfA\bfx} = 0
\]

for all $\bfx\in\mathbb R^d$, where the power to $p-1$ denotes the signed entrywise power. Then by the strong convexity of $\norm*{\cdot}_p^2$ with respect to $\norm*{\cdot}_p$ <cit.>, we have that

\begin{align*}
\norm*{\bfA\bfx-\bfb}_p^2 &\geq \norm*{\bfA\bfx^*-\bfb}_p^2 + 2\norm*{\bfA\bfx^*-\bfb}_p^{2-p}\angle*{(\bfA\bfx^*-\bfb)^{p-1},\bfA\bfx-\bfA\bfx^*} + \frac{p-1}2\norm*{\bfA\bfx-\bfA\bfx^*}_p^2 \\
&\geq \norm*{\bfA\bfx^*-\bfb}_p^2 + \frac{p-1}2\norm*{\bfA\bfx-\bfA\bfx^*}_p^2
\end{align*}

which rearranges to

\[
\norm*{\bfA\bfx-\bfA\bfx^*}_p \leq \sqrt{\frac{\norm*{\bfA\bfx-\bfb}_p^2 - \norm*{\bfA\bfx^*-\bfb}_p^2}{(p-1)/2}} \leq \sqrt{\frac{O(\gamma)\OPT^2}{(p-1)/2}} \leq O(\sqrt\gamma)\OPT.\qedhere
\]

§.§.§ Cost Difference

By Theorem <ref>, if we first obtain a $(1+\gamma)$-approximate solution, then this solution must be within a distance of $O(\sqrt\gamma)\OPT$ from the optimum in the $\ell_p$ norm. Thus, we may in fact just analyze the distortion of the sampling process when restricted to a ball of radius $O(\sqrt\gamma)\OPT$. In this case, we will show that we obtain improved approximations.

We now consider sampling each row independently with probability $1/2$. This can be expressed as multiplying each row by $(1+\sigma_i)$ for independent Rademacher variables $\sigma_i\in\{\pm1\}$ for $i\in[n]$.w Suppose that we find $\hat \bfx\coloneqq \arg\min_{\bfx\in\mathbb R}\hat f(\bfx)$ for

\begin{equation}\label{eq:def-f-hat}
\hat f(\bfx) \coloneqq \sum_{i=1}^n (1+\sigma_i) \abs*{[\bfA\bfx-\bfb](i)}^p.
\end{equation}

\begin{align*}
f(\hat \bfx) - f(\bfx^*) &= [f(\hat \bfx) - \hat f(\hat \bfx)] + [\hat f(\hat \bfx) - \hat f(\bfx^*)] + [\hat f(\bfx^*) - f(\bfx^*)] \\
&\leq [f(\hat \bfx) - \hat f(\hat \bfx)] - [f(\bfx^*) - \hat f(\bfx^*)]
\end{align*}

where we have used that $\hat f(\hat \bfx) \leq \hat f(\bfx^*)$. Thus, the difference in the quality of the two solutions $\hat \bfx$ and $\bfx^*$ depends on the difference in $f-\hat f$ between $\hat \bfx$ and $\bfx^*$. We thus define the difference function

\[
\bar f(\bfx) \coloneqq f(\bfx) - \hat f(\bfx)
\]

for our analyses. Note then that

\[
\bar f(\bfx) = \sum_{i=1}^n \sigma_i \abs*{[\bfA\bfx-\bfb](i)}^p
\]

up to switching the signs on the Rademacher variables $\sigma_i$. The problem then is to bound $\abs*{\bar f(\bfx_1) - \bar f(\bfx_2)}$ for $\bfx_1, \bfx_2$ such that $\norm*{\bfA\bfx_1 - \bfA\bfx_2}_p \leq O(\sqrt\gamma)\OPT$. We will do so by bounding

\[
\abs*{\bar f(\bfx)-\bar f(0)} = \abs*{\sum_{i=1}^n \sigma_i \parens*{\abs*{[\bfA\bfx-\bfb](i)}^p - \abs*{\bfb(i)}^p}}
\]

for $\bfx$ with $\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT$. We thus define

\[
\Lambda \coloneqq \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT} \abs*{\sum_{i=1}^n \sigma_i \parens*{\abs*{\bfA\bfx-\bfb}^p - \abs*{\bfb(i)}^p}}.
\]

We will now bound this quantity.

§.§.§ Chaining Argument

For a subset $S\subseteq[n]$, define

\[
\Lambda_S \coloneqq \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT} \abs*{\sum_{i\in S} \sigma_i \parens*{\abs*{\bfA\bfx-\bfb}^p - \abs*{\bfb(i)}^p}}.
\]

Via a Gaussian comparison theorem (see, e.g., [17, 19]), one can show that

\begin{equation}\label{eq:sym-gp}
\begin{aligned}
\E_\Omega\Lambda_S &\leq \E_\Omega \E_\sigma \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\abs*{\sum_{i\in S}\sigma_i \cdot\parens*{\abs*{[\bfA\bfx-\bfb](i)}^p - \abs*{\bfb(i)}^p}} \\
&= \sqrt{2\pi}\E_\Omega \E_g \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\abs*{\sum_{i\in S}g_i \cdot \parens*{\abs*{[\bfA\bfx-\bfb](i)}^p - \abs*{\bfb(i)}^p}},
\end{aligned}
\end{equation}

where $\sigma_i$ are independent Rademacher variables and $g_i$ are independent standard Gaussian variables.

By Lemma <ref>, we may assume that $\norm*{\bfW^{-1/p}\bfb}_\infty \leq O(1/\eps)\OPT$. Now as done previously, we first handle the coordinates whose Lewis weights are very small:

Suppose that $\norm*{\bfW^{-1/p}\bfb}_\infty \leq O(1/\eps)\OPT$. Let

\[
R \coloneqq \braces*{i\in[n] : \bfw_i^p(\bfA) < \sqrt{\frac{\gamma d}{n}}\frac{\eps^{p}}{n}}.
\]

Then, $\Lambda_R \leq O(\sqrt{\gamma d/n})\OPT^p$ almost surely.

Using (<ref>), we have that

\begin{align*}
\Lambda_R &= \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\abs*{\sum_{i\in R}\sigma_i \cdot\parens*{\abs*{\bfA\bfx-\bfb}^p - \abs*{\bfb(i)}^p}} \\
&\leq \sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\sum_{i\in R} \abs*{\bfA\bfx-\bfb}^p + \abs*{\bfb(i)}^p \\
&\leq O(1)\sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\sum_{i\in R} \bfw_i^p(\bfA)\bracks*{\norm*{\bfA\bfx}_p^p + \frac{\sqrt\gamma^p}{\eps^p}\OPT^p} \\
&\leq O(1)\cdot n \cdot \sqrt{\frac{\gamma d}{n}}\frac{\eps^{p}}{n} \cdot \frac{1}{\eps^p}\OPT^p \\
&\leq O(\sqrt{\gamma d/n})\OPT^p.\qedhere
\end{align*}

For the remaining coordinates, we partition into two sets: those with large coordinates of $\bfb$ and those with small coordinates. Define the sets

\begin{align*}
L_\bfb &\coloneqq \braces*{i\in [n] \setminus R : \abs*{\bfb(i)} \geq C\sqrt\gamma\cdot \bfw_i^p(\bfA)^{1/p}\OPT} \\
S_\bfb &\coloneqq \braces*{i\in [n] \setminus R : \abs*{\bfb(i)} < C\sqrt\gamma\cdot \bfw_i^p(\bfA)^{1/p}\OPT}
\end{align*}

where $R$ is as defined in Lemma <ref> and $C>0$ is a sufficiently large constant. We then bound $\E\Lambda_{[n]\setminus R}$ as

\begin{equation}\label{lem:lambda-n-r-bound}
\begin{aligned}
\E\Lambda_{[n]\setminus R} &\leq \E\sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\abs*{\sum_{i\in [n]\setminus R}\sigma_i \cdot\parens*{\abs*{\bfA\bfx-\bfb}^p - \abs*{\bfb(i)}^p}} \\
&\leq \E\sup_{\norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}\frac1{\sqrt n}\abs*{\sum_{i\in [n]\setminus R}(d\bfw_i)^{1/2}\sigma_i \cdot\parens*{\abs*{\bfW^{-1/p}(\bfA\bfx-\bfb)}^p - \abs*{\bfW^{-1/p}\bfb(i)}^p}}
\end{aligned}
\end{equation}

as done previously, and bound the latter term by bounding the corresponding terms for $L_\bfb$ and $S_\bfb$. We now bound the corresponding Gaussian processes given by

\[
X_\bfy \coloneqq C_S\sum_{i\in S} (d\bfw_i)^{1/2} g_i \cdot \parens*{\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p - \abs*{\bfW^{-1/p}\bfb(i)}^p}
\]

for $\bfy\in \braces*{\bfA\bfx : \bfx\in\mathbb R^d, \norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}$, where $S$ can be either $L_\bfb$ or $S_\bfb$, and $C_S$ is a normalizing constant associated with $S$. Our first task is to bound the associated metric $d_X$.

Let $C_S \coloneqq 1/ \sqrt\gamma^p\OPT^p$. Consider the metric

\begin{align*}
d_X(\bfy,\bfy') &\coloneqq \norm*{X_\bfy-X_{\bfy'}}_2 \\
&= \frac{1}{\sqrt\gamma^{p}\OPT^p}\bracks*{\sum_{i\in S_\bfb}(d\bfw_i)\parens*{\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p - \abs*{\bfW^{-1/p}(\bfy'-\bfb(i))}^p}^2}^{1/2}
\end{align*}

for $\bfy,\bfy'\in \braces*{\bfA\bfx : \bfx\in\mathbb R^d, \norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}$. Then,

\[
d_X(\bfy,\bfy') \leq O(\sqrt d)\frac{\norm*{(\bfW^{-1/p}\bfy-\bfW^{-1/p}\bfy')\vert_{S_\bfb}}_\infty^{p/2}}{\sqrt\gamma^{p/2}\OPT^{p/2}}
\]

The proof follows from a straightforward adaptation of Lemma <ref>.

Bounding the metric for $L_\bfb$ requires a different approach than the metric used in [46].

Let $C_S \coloneqq 1/\sqrt\gamma\OPT^p$. Consider the metric

\begin{align*}
d_X(\bfy,\bfy') &\coloneqq \norm*{X_\bfy-X_{\bfy'}}_2 \\
&= \frac{1}{\sqrt\gamma\OPT^p}\bracks*{\sum_{i\in L_\bfb}(d\bfw_i)\parens*{\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p - \abs*{\bfW^{-1/p}(\bfy'-\bfb(i))}^p}^2}^{1/2}
\end{align*}

for $\bfy,\bfy'\in \braces*{\bfA\bfx : \bfx\in\mathbb R^d, \norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}$. Then,

\[
d_X(\bfy,\bfy') \leq O(\sqrt d)\frac{\norm*{(\bfW^{-1/p}(\bfy-\bfy'))\vert_{L_\bfb}}_\infty^{p/2}}{\sqrt\gamma^{p/2}\OPT^{p/2}}
\]

Note that

\begin{align*}
\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p &= \parens*{\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)} \pm \abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}}^p \\
&= \abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^p \pm O(1)\abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^{p-1}
\end{align*}

since $i\in L_\bfb$. Thus,

\begin{align*}
\abs*{\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p - \abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^p} &\leq O(1)\abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^{p-1}
\end{align*}

so we have that

\begin{align*}
&\sum_{i\in L_\bfb}(d\bfw_i)\parens*{\abs*{[\bfW^{-1/p}(\bfy-\bfb)](i)}^p - \abs*{\bfW^{-1/p}(\bfy'-\bfb(i))}^p}^2 \\
\leq~&O\parens*{d}\norm*{(\bfW^{-1/p}(\bfy-\bfy'))\vert_{L_\bfb}}_\infty^p \sum_{i\in L_\bfb}\bfw_i \abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}^{2-p}\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^{2p-2}.
\end{align*}

Now note that $\frac{p}{2-p} > 1$ and $\frac{p}{2p-2} > 1$ are Hölder conjugates for $1 < p < 2$. Then, Hölder's inequality on the sum yields

\begin{align*}
&\sum_{i\in L_\bfb}\bfw_i \abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}^{2-p}\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^{2p-2} \\
\leq~&\parens*{\sum_{i\in L_\bfb}\bfw_i \abs*{[\bfW^{-1/p}(\bfy-\bfy')](i)}^{p}}^{\frac{2-p}{p}}\parens*{\sum_{i\in L_\bfb}\bfw_i\abs*{[\bfW^{-1/p}(\bfy'-\bfb)](i)}^{p}}^{1 - \frac{2-p}{p}} \\
\leq~&\norm*{\bfy-\bfy'}_p^{2-p} \norm*{\bfy'-\bfb}_p^{2p-2} \leq O(\sqrt\gamma^{2-p})\OPT^p.
\end{align*}

Combining these bounds yields the claim.

Now given the bounds on the metrics, we can finish by using the entropy bounds from Section <ref> and then calculating Dudley's entropy integral as done in Section <ref>.

We have that

\[
\E\Lambda \leq O(1) \bracks*{\gamma\frac{d}{n}(\log d)^2(\log n)}^{1/2}\OPT^p
\]

\[
\Pr\braces*{\Lambda \geq C\sqrt{\frac{\gamma d}{n}}\bracks*{(\log d)\sqrt{\log n}+z}\OPT^p} \leq 2\exp(-z^2).
\]

Our proof is very similar to that of Lemma <ref> for $1 < p < 2$. Let

\[
T = \braces*{\bfA\bfx : \bfx\in\mathbb R^d, \norm*{\bfA\bfx}_p \leq O(\sqrt\gamma)\OPT}.
\]

Then, as reasoned in Lemma <ref> for $1 < p < 2$, for the Gaussian process $X_\bfy$ associated with $S_\bfb$ (see Lemma <ref>),

\begin{align*}
\log E(T, d_X, t) &\leq \log E(\Theta(\sqrt\gamma d^{-1/p}\OPT)B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, (t/\sqrt d)^{2/p}\sqrt\gamma\OPT) \\
&= \log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, \Theta(t^{2/p})).
\end{align*}

Similarly, for the Gaussian process $X_\bfy$ associated with $L_\bfb$ (see Lemma <ref>),

\[
\log E(T, d_X, t) \leq \log E(B_{\bfw,p}, \norm*{\cdot}_{\bfw,q}, \Theta(t^{2/p})).
\]

Then by Theorem <ref>, following the calculation in Lemma <ref> for $1 < p < 2$, we have that

\begin{align*}
\E\sup_{\bfy\in T}\abs*{X_\bfy} \leq O(1)\sqrt d(\log d)\sqrt{\log n}
\end{align*}

for both of the Gaussian processes. We combine this with (<ref>) to obtain

\begin{align*}
\E\Lambda_{[n]\setminus R} &\leq O(1)\frac1{\sqrt n}\bracks*{\frac1{C_{L_\bfb}}\sqrt d(\log d)\sqrt{\log n} + \frac1{C_{S_\bfb}}\sqrt d(\log d)\sqrt{\log n}} \\
&\leq O(1) \frac1n\sqrt d(\log d)\sqrt{\log n}[\sqrt\gamma\OPT^p + \sqrt\gamma^p\OPT^p] \\
&\leq O(1) \bracks*{\gamma\frac{d}{n}(\log d)^2(\log n)}^{1/2}\OPT^p
\end{align*}

where $C_{L_\bfb} = 1/\sqrt\gamma\OPT^p$ and $C_{S_\bfb} = 1/\sqrt\gamma^p\OPT^p$ are the $C_S$ constants for $S = L_\bfb$ and $S = S_\bfb$ (see Lemmas <ref> and <ref>). Finally, together with Lemma <ref>, we conclude that

\[
\E\Lambda \leq O(1) \bracks*{\gamma\frac{d}{n}(\log d)^2(\log n)}^{1/2}\OPT^p.
\]

In a similar manner as Lemma <ref>, the $d_X$-diameter of $T$ is $O(\sqrt d)$ so Theorem <ref> gives us that

\[
\Pr\braces*{\Lambda \geq C\sqrt{\frac{\gamma d}{n}}\bracks*{(\log d)\sqrt{\log n}+z}\OPT^p} \leq 2\exp(-z^2).\qedhere
\]

§.§.§ Iterative Size Reduction

We will now assemble the previous lemmas into an iterative size reduction argument.

Because Lemma <ref> differs from Lemma <ref> only by a factor of $\sqrt\gamma$, one can show as done in Theorem <ref> that $m$ entries need to be read to produce a $(1+\eps)$-approximation, for

\[
m = O(1)\frac{\gamma d}{\eps^{2}}\bracks*{(\log d)^2\log n + \log\frac1\delta},
\]

as long as it can be shown that reading $m$ entries gave a $(1+\gamma)$-approximation. Now let

\[
C = O(1)\bracks*{(\log d)^2\log n + \log\frac1\delta}
\]

and suppose that reading $Cd/\eps^\beta$ is sufficient to obtain a $(1+\eps)$-approximation (for instance, Theorem <ref> applied directly gives that $\beta\leq 2$). Then, if we set $\gamma = \eps^{2/(1+\beta)}$, then making $Cd/\gamma^\beta$ gives a $(1+\gamma)$-approximation. Furthermore, note that

\[
C\frac{\gamma d}{\eps^2} = Cd \eps^{\frac{2}{1+\beta} - 2} = Cd \eps^{\frac{2-2(1+\beta)}{1+\beta}} = Cd \eps^{-\frac{2\beta}{1+\beta}} = C\frac{d}{\gamma^\beta}
\]

so reading $Cd/\gamma^\beta = C\gamma d/\eps^2$ entries is in fact sufficient to obtain a $(1+\eps)$-approximation for this setting of $\gamma$. We may now iterate this argument. Consider the sequence $\beta_i$ given by

\[
\beta_1 = 2, \qquad \beta_{i+1} \to \frac{2\beta_i}{1+\beta_i}.
\]

It can be checked that the solution is

\[
\beta_i = \frac{2^i}{2^i - 1} = 1 + \frac1{2^i-1},
\]

so applying this argument $O(\log\log\frac1\eps)$ times yields that $\beta_i \leq 1 + O(1/\log(\frac1\eps))$ which means that reading only $O(1)Cd/\eps$ entries suffices.

We have shown the following:

Let $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $1 < p < 2$. Let

\[
m = O(1)\frac{d}{\eps}\bracks*{(\log d)^2\log n + \log\frac1\delta}\log\frac1\delta.
\]

Then, with probability at least $1-\delta$, Algorithm <ref> returns $\tilde\bfx\in\mathbb R^d$ such that

\[
\norm*{\bfA\tilde\bfx-\bfb}_p \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx-\bfb}_p
\]

and reads at most $m$ entries of $\bfb$.

§.§ Removing the log n

As shown in Section <ref>, Algorithm <ref> gives a relative error approximation to $\min \norm{\bfA\bfx-\bfb}_p$ using $\eps^{-2} d^{\max(1,p/2)} \poly(\log d, \log n)$ samples from $\bfb$. However, the dependence on $\log n$ is undesirable, and we now show how to replace the $\log n$ dependence with a $\log(d/\eps)$ dependence.

§.§.§ Bounding n

Due to the use of the sophisticated chaining argument, Theorem <ref> has an undesirable dependence on $\log n$. By replacing this net construction with a standard $\epsilon$-net, we can still show that Algorithm <ref> solves $\ell_p$ regression to $(1+\epsilon)$ error, albeit using $\tilde{O}(d^{\max(1,p/2)+1}/\poly(\epsilon))$ samples (without dependence on $\log n$). This bound is loose by a $d$ factor as compared to Theorem <ref> due to the simpler net construction.
It can, however, be used as a preprocessing step for Theorem <ref>, to ensure that $n$ is only $\poly(d/\eps)$ and so that the $\log n$ dependence is removed — after sampling by Lewis weights, we are left with a problem where $n = \poly(d,1/\epsilon)$. We can then apply Theorem <ref> to further subsample this problem, achieving $\tilde{O}(d^{\max(1,p/2)}/\poly(\epsilon))$ sample complexity. See Algorithm <ref> for pseudocode for the full approach.

We use a more standard Lewis weight sampling algorithm for this result:

Let $\bfA\in\mathbb R^{n\times d}$. Let $\tilde\bfw_i^p(\bfA) \geq \bfw_i^p(\bfA)$ be Lewis weight upper bounds, and let $m$ be a parameter. Define probabilities

\[
\bfp_i \coloneqq \min\{1, m\cdot d^{0\lor (p/2-1)}\tilde\bfw_i^p(\bfA)\}
\]

Then, $\bfS$ is a Lewis weight sampling matrix if for each $i\in[n]$, $\bfS(i,i) = \bfp_i^{-1/p}$ independently with probability $\bfp_i$ and $0$ otherwise.

Relative error $\ell_p$ regression – No Assumptions
input: Matrix $\bfA \in \R^{n \times d}$, measurement vector $\bfb \in \R^n$.

output: Approximate solution $\tilde \bfx \in \R^d$ to $\min_\bfx\|\bfA\bfx - \bfb\|_p$.
[1]
Let $\bfS_1 \in \R^{m_1 \times n}$ be generated via Def. <ref> for $m_1=O \left (d^{\max(1,p/2)+1} \cdot \frac{\log(1/\epsilon\delta)}{\epsilon^{2+p}}\right )$.
Run Algorithm <ref> with inputs $\bfS_1\bfA$ and $\bfS_1\bfz$ to obtain $\tilde\bfx$.
$\tilde \bfx$

Let $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $0 < p < \infty$. Let

\[
m = \begin{dcases}
O(1)\frac{d}{\eps^{2}}\bracks*{(\log d)^2\log(d/\eps) + \log\frac1\delta}\log\frac1\delta & p\in(0,1) \\
O(1)\frac{d}{\eps}\bracks*{(\log d)^2\log(d/\eps) + \log\frac1\delta}\log\frac1\delta & p\in(1,2) \\
O(1)\frac{d^{p/2}}{\eps^{p}}\bracks*{(\log d)^2\log(d/\eps) + \log\frac1\delta}\log\frac1\delta & p\in(2,\infty)
\end{dcases}
\]

Then, with probability at least $1-\delta$, Algorithm <ref> returns $\tilde\bfx\in\mathbb R^d$ such that

\[
\norm*{\bfA\tilde\bfx-\bfb}_p \leq (1+\eps)\min_{\bfx\in\mathbb R^d}\norm*{\bfA\bfx-\bfb}_p
\]

and nonadaptively reads at most $m$ entries of $\bfb$.

Theorem <ref> will follow from an analog of Lemma <ref>, when $m=O \left (d^{\max(1,p/2)+1} \cdot \frac{\log(1/\epsilon\delta)}{\epsilon^{2+p}}\right )$. In particular:

Let $\bfA \in \R^{n \times d}$, $\bfb \in \R^n$, and $0 < p < \infty$. Let $\bfz = \bfb - \bfA\bfx^*\in \R^n$, with $\norm{\bfz}_p = \OPT$. Let

\[
\mathcal B\coloneqq \braces*{i\in[n] : \frac{|\bfz(i)|^p}{\OPT^p} \ge \frac{d^{\max(0,p/2-1)} \cdot \bfw_i^p(\bfA)}{\epsilon^p}}
\]

Let $\bar \bfz \in \R^n$ be equal to $\bfz$ but with all entries in $\mathcal{B}$ set to $0$. Let $\bfS_1 \in \R^{m_1 \times n}$ be formed as in Line <ref> of Algorithm <ref>. With probability at least $1-\delta$, for all $\bfx \in \R^d$ with $\norm{\bfA\bfx}_p = O(\OPT)$,

\begin{align*}
\left | \left \| \bfS_1\bfA \bfx -\bfS_1 \bar \bfz \right \|_p^p - \norm{\bfA\bfx - \bar \bfz}_p^p \right | = O(\epsilon) \cdot \OPT^p.
\end{align*}

We first use Lemma <ref> to prove Theorem <ref>. Then we prove the lemma.

By restricting to the $m_1$ rows sampled by Line <ref> of Algorithm <ref>, we have a $\poly(d/\eps)\times d$ matrix $\bfS_1\bfA$ such that

\begin{align*}
\left | \left \| \bfS_1\bfA \bfx -\bfS_1 \bar \bfz \right \|_p^p - \norm{\bfA\bfx - \bar \bfz}_p^p \right | = O(\epsilon) \cdot \OPT^p
\end{align*}

by Lemma <ref>. In combination with Lemmas <ref> and <ref>, Theorem <ref> and Lemma <ref> give that with probability at least $1-O(\delta)$, there is some fixed $C\geq0$ with $C = O(\OPT^p)$ such that for all $\bfx \in \R^d$ with $\norm{\bfA\bfx}_p = O(\OPT)$,

\begin{align*}
|\norm{\bfS_1\bfA\bfx-\bfS_1\bfz}_p^p - \norm{\bfA\bfx-\bfz}_p^p - C| = O(\epsilon) \cdot \OPT^p.
\end{align*}

Following the proof of Theorem <ref>, we then have that if $\bar \bfx \in \R^d$ satisfies $\norm{\bfS_1 \bfA \bar \bfx -\bfS_1 \bfz}_p \le (1+\epsilon/2) \cdot \min_{\bfx \in \R^d} \norm{\bfS_1 \bfA \bfx -\bfS_1 \bfz}_p$ then $\norm{\bfA \bar \bfx - \bfz}_{p} \le (1+\epsilon) \cdot \OPT$. We can verify that $\bar \bfx$ as computed in Line <ref> of Algorithm <ref> does indeed satisfy this near optimality condition with probability $1-O(\delta)$ by Theorem <ref>.  Overall, with probability at least $1-O(\delta)$, we have that $\bar \bfx$ as computed in Line <ref> of Algorithm <ref> satisfies $\norm{\bfA \bar \bfx - \bfz}_{p} \le (1+\epsilon) \cdot \OPT$.

For simplicity, we assume via scaling that $\norm{\bfz}_p \le 1$ and $\OPT = \Theta(1)$. By a standard argument (see, e.g., Lemma 2.5 of [6]), it suffices to show that with high probability, for all $\bfy' \in \mathcal{N}_\epsilon$,

\begin{align}\label{eq:netError}
\left | \left \|\bfS_1\bfy' -\bfS_1 \bar \bfz \right \|_p^p - \norm{\bfy ' - \bar \bfz}_p^p \right | \le \epsilon,
\end{align}

where $\mathcal{N}_\epsilon$ is an $\epsilon$-net in the $p$-norm over $\{\bfA\bfx: \norm{\bfA\bfx}_p \le 1\}$. By a standard volume argument (see, e.g., Lemma 2.4 of [6]), it is known that one can construct this net such that $\log |\mathcal{N}_\epsilon| = O(d\log(1/\epsilon))$.

We prove (<ref>) via a Bernstein inequality and union bound. We have $\E[\norm{\bfS _1\bfy' -\bfS_1 \bar \bfz}_p^p] = \norm{\bfy ' - \bar \bfz}_p^p = O(1)$. Additionally, by definition, $|\bar \bfz(i)|^p \le \frac{d^{\max(0,p/2-1)}\cdot \bfw_i^p(\bfA)}{\epsilon^p}$ for all $i$. Similarly, by Lemma <ref> which bounds the $\ell_p$ sensitivities by the Lewis weights, $ |\bfy '(i)|^p \le d^{\max(0,p/2-1)}\cdot \bfw_i^p(\bfA) $. Overall,

\begin{align*}
|\bfy'(i) - \bar \bfz(i)|^p = O \left ( \frac{d^{\max(0,p/2-1)}\cdot \bfw_i^p(\bfA)}{\epsilon^p}\right)
\end{align*}

By the construction of $\bfS_1$ (Definition <ref>), this gives

$$|[\bfS_1\bfy'-\bfS_1\bar \bfz](i)|^p \le \frac{d}{m \cdot \bfw_i^p(\bfA)} \cdot \frac{d^{\max(0,p/2-1)}\cdot \bfw_i^p(\bfA)}{\epsilon^p} =O \left (\frac{d^{\max(1,p/2)}}{m \cdot \epsilon^p} \right ),$$

and so applying a Bernstein bound

\begin{align*}
\Pr \left [|\norm{\bfS _1\bfy' -\bfS_1\bar \bfz}_p^p - \norm{\bfy '-\bar \bfz}_p^p |\le \epsilon \right ] & \le 2 \exp \left (- \Omega \left ( \frac{\epsilon^2}{{d^{\max(0,p/2-1)}} / { m \cdot \epsilon^p}} \right ) \right ) \\
&= 2 \exp \left (- \Omega \left ( \frac{\epsilon^{2+p} \cdot m}{d^{\max(0,p/2-1)}} \right ) \right ).
\end{align*}

\[
m = O \left (\frac{d^{\max(1,p/2)}}{\epsilon^{2+p}} \cdot \log \frac{|\mathcal{N}_\epsilon|}{\delta} \right ) \leq O \parens*{\frac{d^{\max(1,p/2)+1}}{\epsilon^{2+p}} \log\frac1{\eps\delta}}
\]

and applying a union bound gives that this bound holds with high probability for all $\bfy ' \in \mathcal{N}_\epsilon$. This completes the proof.

§.§ Optimal Dependence on delta with Additive Error

In this section, we show that with a knowledge of some overestimate $E$ of $\OPT$, we can obtain an optimal dependence on $\delta$ with an additive error of $\eps E$. Our approach is as follows. We first use Theorem <ref> to compute a constant factor approximation $\bfx_c$. Now with the knowledge of $E$ as well as $\bfx_c$, we can explicitly clip the large entries of $\bfb$ as done in Lemma <ref>, where $\bfz$ can be replaced by $\bfz = \bfb - \bfA\bfx_c$ for a $\bfx_c$ which can be found with high probability by Theorem <ref>. In this case, we can explicitly approximate the cost of this clipped vector by sampling, which means we can collect $O(\log\frac1\delta)$ runs of a constant probability solution, and output the best solution among $O(\log\frac1\delta)$ many repetitions. This achieves an optimal dependence on $\delta$ up to constant factors due to a lower bound of <cit.>, as this hard instance has an $\OPT$ that is fixed up to a $\Theta(1)$ factor, with high probability.

Note first that in Lemma <ref>, we may replace $\OPT$ by $E$, $\bfz = \bfb-\bfA\bfx^*$ by $\bfz = \bfb-\bfA\bfx_c$, and $\bar\bfz$ by the explicitly clipped $\bar\bfz$ using $E$ instead of $\OPT$, to obtain the following:
Let $\bfA\in\mathbb R^{n\times d}$ and $\bfb\in\mathbb R^n$, and let $0 < p < \infty$. Let $\bfz = \bfb - \bfA\bfx_c\in \R^n$ satisfy $\norm*{\bfz}_p = O(\OPT)$. Let $E\geq\OPT$ and let

\[
\mathcal B \coloneqq \braces*{i\in[n] : \frac{|\bfz(i)|^p}{E^p} \ge \frac{d^{\max(0,p/2-1)} \cdot \bfw_{i}^p(\bfA)}{\epsilon^p}}
\]

Let $\bar \bfz \in \R^n$ be equal to $\bfz$ but with all entries in $\mathcal{B}$ set to $0$. Then for all $\bfx \in \R^d$ with $\norm{\bfA\bfx}_p = O(E)$,

$$\left | \norm{\bfA\bfx-\bfz}_p^p - \norm{\bfA\bfx-\bar \bfz}_p^p - \norm*{\bfz-\bar\bfz}_p^p \right | = O( \epsilon) \cdot E^p.$$

With this modification, the following version of Lemma <ref> holds with the same proof:

Consider the setting of Lemma <ref>. Let $\bfS \in \R^{m \times n}$ be formed as in Line <ref> of Algorithm <ref>. Let

\[
\Lambda \coloneqq \sup_{\norm*{\bfA\bfx}_p \leq O(E)}\abs*{\norm*{\bfS\bfA\bfx-\bfS\bar\bfz}_p^p - \norm*{\bfA\bfx-\bar\bfz}_p^p}.
\]

\[
\Pr\braces*{\Lambda \geq \bracks*{C\frac1n \frac{d^{1\lor (p/2)}}{\eps^{0\lor(p-2)}}}^{1/2} \bracks*{(\log d)\sqrt{\log n} + z} E^p} \leq 2\exp(-z^2)
\]

We then obtain the following theorem, which has an improved dependence on $\delta$.

Let $E\geq \OPT = \min_\bfx\norm*{\bfA\bfx-\bfb}_p$ be known. There is an algorithm which, with probability at least $1-\delta$, reads only $M$ entries of $\bfb$ and outputs $\tilde\bfx$ satisfying

\[
\norm*{\bfA\tilde\bfx-\bfb}_p \leq \OPT + \eps E
\]

\[
m = \begin{dcases}
O(1)\frac{d}{\eps^{2}}(\log d)^2\log(d/\eps) \log\frac1\delta & p\in(0,2) \\
O(1)\frac{d^{p/2}}{\eps^{p}}(\log d)^2\log(d/\eps) \log\frac1\delta & p\in(2,\infty)
\end{dcases}
\]

We first prepare $\bfz$ as required in Lemma <ref> using Theorem <ref>. This requires reading only $O(d(\log d)^3\log\frac1\delta)$ entries. Then by using our knowledge of $E$ and $\bfz$, we may replace $\bfb$ by the clipped residual $\bar\bfz$ at a cost of only an additive $\eps\cdot E$ error by Lemma <ref>. Next, we sample $\ell = O(\log\frac1\delta)$ independent trials of $\bfS^{(i)}$ for $i\in[\ell]$ each with a constant probability success guarantee and solve $\min_\bfx\norm{\bfS^{(i)}\bfA\bfx-\bfS^{(i)}\bar\bfz}_p$ to obtain $\ell$ candidate solutions $\tilde\bfx^{(i)}$. Note that each $\bfS^{(i)}$ only requires reading

\[
O(1)\frac{d^{1\lor(p/2)}}{\eps^{2\lor p}}(\log d)^2\log(d/\eps)
\]

entries of $\bfb$, so we only require $O(m)$ entries altogether. Then by Chernoff bounds, with probability at least $1-\delta$, at least $9/10$ of these solutions satisfy

\[
\norm*{\bfA\tilde\bfx^{(i)}-\bfb}_p \leq \OPT + \eps\cdot E
\]

and satisfy $\norm{\bfA\tilde\bfx^{(i)}}_p \leq O(E)$, and we may discard any solution which fails to have $\norm{\bfA\tilde\bfx^{(i)}}_p \leq O(E)$ by using our knowledge of $E$. Finally, we sample a single high probability sketch $\bfS$ which satisfies the guarantee of Lemma <ref> with probability $1-\delta$. This single sketch requires only reading $O(m)$ entries of $\bfb$. Then, up to an $\eps E$ additive error we may approximately evaluate the cost of each of the $\ell$ solutions $\tilde\bfx$ such that $\norm{\bfA\tilde\bfx^{(i)}}_p \leq O(E)$ and select the best candidate using $\bfS$, by using $\norm*{\bfS\bfA\tilde\bfx^{(i)}-\bfS\bar\bfz}_p$ as the proxy cost. We conclude by adjusting $\delta$ and $\eps$ by constant factors.

§ SENSITIVITY BOUNDS FOR M-ESTIMATORS

In this section, we present our new sensitivity bounds. In Section <ref>, we collect basic definitions and lemmas concerning $M$-estimators. Section <ref> develops basic notions for sensitivity sampling for $M$-estimators. In Section <ref>, we describe our efficient algorithm for computing sensitivities for a broad class of $M$-estimators. In Section <ref>, we show that a variation on our efficient algorithm can be used to show an existential bound of $O(d^{\max\{1,p_M/2\}}\log n)$ total sensitivity for the same class of $M$-estimators. Finally, in Section <ref>, we show that the Tukey loss has a total sensitivity of $\Omega(d\log n)$, and that the Huber loss has a total sensitivity of $\Omega(d\log\log n)$.

§.§ Geometry of M-norms

In this section, we define $M$-norms and collect some of their geometric properties. This is a slight generalization of Section 4.1 of [21] which allows for a broader class of $M$-norms (namely with a relaxed polynomial lower bound condition). With applications to active regression in mind, we also slightly generalize the results to handle translations by a single vector $\bfb$, which can be taken to be $0$ to retrieve the original results.

Let $M:\mathbb R_{\geq0}\to\mathbb R_{\geq0}$ be increasing. If there exist constants $p>0$ and $c_U\geq 1$ such that for all $y>x$,

\[
\frac{M(y)}{M(x)} \leq c_U \parens*{\frac{y}{x}}^p,
\]

then we say that $M$ is polynomially bounded above with degree $p$ and constant $c_U$. Similarly, if there exists constants $q>0$ and $c_L\geq 1$ such that for all $y>x$,

\[
\frac{M(y)}{M(x)} \geq c_L \parens*{\frac{y}{x}}^q,
\]

then we say that $M$ is polynomially bounded below with degree $q$ and constant $c_L$.

As noted in [21], it can be shown that convex functions are polynomially bounded below with degree $1$.

Throughout this work, we will consider the constants $p, q, c_U, c_L$ in Definition <ref> to be absolute constants that don't depend on other parameters under consideration.

We define the $M$-norm as follows. Note that despite our abuse of notation and terminology, the $M$-norm need not be an actual norm.

Let $M:\mathbb R_{\geq0}\to\mathbb R_{\geq0}$ be such that

* $M(0) = 0$

* $M$ is nondecreasing

* $M$ is polynomially bounded above with degree $p_M$ and constant $c_U$ (see Definition <ref>)

Let $\bfw\in\mathbb R^n$ be a set of weights such that

\[
\bfw_i \geq 1
\]

for all $i\in[n]$. Then, we define the $M$-norm of a vector $\bfx\in\mathbb R^n$ as

\[
\norm*{\bfx}_{M,\bfw} \coloneqq \bracks*{\sum_{i=1}^n \bfw_i M(\abs*{\bfx_i})}^{1/p_M}.
\]

If $\bfw$ is the vector of all ones, we simple write $\norm*{\bfx}_{M}$ for $\norm*{\bfx}_{M,\bfw}$. If $M(x) = \abs*{x}^p$ for some $p>0$, then we write $\norm*{\bfx}_{p,\bfw}$ for $\norm*{\bfx}_{M,\bfw}$.

Let $\bfA\in\mathbb R^{n\times d}$ and let $\mathcal V = \Span(\bfA)$. Let $M:\mathbb R_{\geq0}\to\mathbb R_{\geq0}$ satisfy the conditions of Definition <ref>, and let $\bfw\geq\mathbf{1}_n$ be a set of weights. Define the ball $\mathcal B^{M,\bfw}_\rho$ of radius $\rho>0$ as

\[