I. INTRODUCTION

We consider a classic active learning problem: given a design matrix A ∈ R n×d and query access to entries of an unknown target (measurement) vector b ∈ R n , how can we compute an approximate minimizer of the regression problem min x∈R d Ax -b while querying as few entries of b as possible? This problem arises in applications where labeled data is expensive: viewing a single entry of b might require running a survey, physical experiment, or time-intensive com-Cameron Musco's work on this project was supported in part by NSF Grants 2046235 and 1763618, along with an Adobe Research Grant. Christopher Musco was supported by NSF Grant 2045590. David P. Woodruff and Taisuke Yasuda were supported by ONR grant N00014-18-1-2562 and a Simons Investigator Award. 1 Extended abstract; full version available at https://arxiv.org/abs/2111. 04888.

puter simulation [44], [45]. Concretely, we study the following problem for general vector norms 2 • : Problem I.1. For A ∈ R n×d , b ∈ R n , and accuracy parameter 0 < ≤ 1, find x ∈ R d satisfying:

Ax -b , while reading as few of the entries {b(1), . . . , b(n)} of the target vector b as possible. 3 Notably, the formulation of Problem I.1 makes no assumptions on A and b. For example, we do not assume that there exists a ground truth x and that Ax -b 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 [9], [33], [44].

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 I.1 inherently requires choosing which entries of b to query in a randomized way: an adversary can easily "fool" any deterministic algorithm by concentrating error in Ax -b on the indices of b that will be deterministically queried.

A. Prior Work

Euclidean Norm. Problem I.1 is fully understood when the error is measured in the 2 norm, w 2 = n i=1 |w i | 2 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 x be the minimizer of this sampled problem, which only involves a fraction of the entries in b. I.e., letting S ∈ R m×n be a sampling matrix with m < n rows (S has one non-zero entry per row), set x = arg min x SAx -Sb . When constraints are selected with probability proportional to the statistical leverage scores of A's rows, Problem I.1 can be solved with O(d/ • log d) 2 Our work will also extend to other loss functions of the form n i=1 M ([Ax -b] i ) that are not necessarily norms. 3 In principal, entries of b 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 I.1 and those studied in this paper are non-adaptive. samples, and thus O(d/ • log d) queries to b [21], [46], [53]. 4Using tools from spectral graph sparsification [6], [35], Chen and Price recently improved the leverage score sampling result to O(d/ ), which is optimal [12].

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 A and b are undesirable [16], [17], [30]. Methods for solving Problem I.1 in the 2 norm also yield robust methods for interpolating sparse Fourier functions, bandlimited and multiband functions, and for data-efficient kernel learning [4], [11], [23].

Other Norms. Much less was known about Problem I.1 beyond the 2 norm until recent work of Chen and Dereziński [10], which proves an upper bound of

This result is tight up to the log d factor. A similar result is obtained in [43]. Chen and Dereziński also prove a result for p norms, w p = (

As for the 2 norm, the results for 1 and p are obtained by subsampling rows of the regression problem independently at random. However, instead of sampling with probabilities proportional to the leverage scores, [10], [43] employ a natural generalization of these scores known as the p Lewis weights [18]. 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 b is a {-1, 1} label vector, and the error is measured via the logistic loss, Munteanu et al. [42] show that poly(d, μ, 1/ ) samples suffice, where μ is a complexity measure of A. This bound has recently been tightened to Õ(dμ 2 / 2 ) [40], using Lewis weight sampling. 5 For other loss functions, such as the Tukey loss and Huber's M -estimators for robust regression [27], we are not aware of any known results solving Problem I.1. 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.

B. Our Contributions a) p Active Regression.:

Our first main result is a new algorithm for solving Problem I.1 for the p norm for any 0 < p < ∞foot_2 . While near-optimal bounds are known for p ∈ {1, 2} [10], [12], [43], the problem is far from settled for all other p. Previously, active p regression for p > 2 and 0 < p < 1 had no known nontrivial algorithms with (1 + ) relative error, and the only known approach was to read all n entries of b and solve the problem using offline results. A natural question is whether a sublinear query complexity is possible in these regimes. For p ∈ (1, 2), [10] achieved an algorithm making O(d 2 / 2 •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.

Theorem I.2 (Main Result for Active p Regression). Given 0 < p < ∞, A ∈ R n×d , and query access to b ∈ R n , there is an algorithm that solves Problem I.1 for the p -norm with probability 99/100 which makes m queries in b, where

.

We complement our algorithmic result with various new lower bounds which show the tightness of our algorithm. For p ∈ (0, 2), our dependence on d and in the query complexity are simultaneously tight up to polylogarithmic factors; we show an Ω(d/ 2 ) lower bound for p ∈ (0, 1) and an Ω(d/ ) lower bound for p ∈ (1, 2). For p > 2, our dependence on d is tight due to a lower bound of Ω(d p/2 ) which we show, while our dependence is off by at most factor of due to an Ω( 1-p ) lower bound for the one-dimensional p power means problem in Theorem 3 of [19]. Note that our active regression lower bounds for p ∈ (0, 2) improve this previous power means lower bound.

Notably, we achieve a linear dependence on for p ∈ (1, 2), which is perhaps surprising given that all previous known approaches to dimension reduction for p regression relied on preserving the p norm of all vectors in a subspace up to (1± ) factors [18], which requires Ω(d/ 2 ) dimensions [37]. It also demonstrates a separation in the query complexity for p ≤ 1 and 1 < p < 2, due to a lower bound of Ω(d/ 2 ) for p = 1 [10], [43] as well as for p ∈ (0, 1) which we show.

Note that Theorem I.2 is stated to solve Problem I.1 with constant probability, 99/100. In general, we show how to obtain 1δ probability with dependence on δ that is only polylogarithmic in 1/δ. In fact, we show that any algorithm that simply samples rows of the regression problem and solves the sampled problem must suffer a 1/δ p-1 dependence. Indeed, such a loss is seen in the algorithm of [10] for p ∈ (1, 2). Thus, a success probability boosting routine, as we give in our work, is required to obtain an O(log(1/δ)) dependence.

b) Sensitivity Bounds and Active Regression for General Losses.: We show that our approach to solving Problem I.1 for p norms generalizes to a broad class of loss functions known as M -estimators [15], which take the form

The only properties that we require are that we can (1) compute a constant factor approximation to Problem I.1 (2) the loss function obeys approximate variants of the triangle 1)).

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

Theorem I.3 (Main Result for Sensitivity Bounds). Let A ∈ R n×d and let M be an M -estimator loss with at most degree p growth. There is an algorithm which, with probability at least 99/100, computes M -sensitivity upper bounds which sum to at most O(d 1∨(p/2) log 2 n + τ )foot_3 in time at most Õ(nnz(A) + nd O(1) /τ ).

Our approach to sensitivity bounds only relies on hashing and the computation of p Lewis weights [18], [24], and avoids the computation of Löwner-John ellipsoids. This allows for input sparsity time algorithms, and answers an open question of [52] on avoiding Löwner-John ellipsoids in the computation of sensitivities. Note that our dependence on d matches the sensitivity bounds for the p loss and is thus tight. We also show 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∨(p/2) log n) for these Mestimators; this is in fact tight for the Tukey loss by our lower bound of Ω(d log n). Thus, we obtain the first tight bounds on the sum of sensitivities, for losses other than p . Overall, we make significant progress on generalizing the theory of matrix approximation beyond p losses to handle general Mestimators, which is a direction that has recently received much attention [13]- [15], [26], [51], [52].

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 [10]. Note that prior to our work, no sublinear query complexity was known for any M -estimator regression, besides the 2 and 1 losses.

Furthermore, our new sensitivity bounds imply significant improvements in previous results using sensitivity sampling, beyond active regression, including Orlicz norm subspace embeddings [50] and robust subspace approximation [14]. We believe that our general technique here will find other further applications, and leave it as an open question to do so. c) Subspace Embeddings for Orlicz Norms.: Orlicz norms can be viewed as scale-invariant extensions of Mestimators, and have recently attracted attention as a general class of norms that admit efficient dimensionality reduction results [3], [50]. In particular, [50] apply sensitivity sampling to obtain subspace embeddings for Orlicz norms, which yields a small weighted subset à of rows of a matrix A ∈ R n×d such that Ãx = (1 ± ) Axfoot_4 for all x ∈ R d . However, the number of rows required by [50] 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∨(p/2+1) poly(log n, -1 ) rows, for any Orlicz norm with a polynomial growth bound of degree p. d) Robust Subspace Approximation.: The robust subspace approximation problem generalizes the classical low rank approximation problem of finding a rank k projection X minimizing AX -A F by replacing the Frobenius norm with an extension of M -estimators to matrix norms. [14] 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, -1 )×poly(k, log n, -1 ) instance. We also extend their method beyond quadratic growth, to any degree p polynomial growth.

e) 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 p active regression results. This is attributed to the use of our net argument for general M -estimators, whereas our p active regression results can make use of more sophisticated chaining arguments of [7], [34], [48]. 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:

Definition I.4 (Huber loss [31]). The Huber loss of width τ ≥ 0 is defined as

and the Huber normfoot_6 is defined as y H := n i=1 H(y(i)). The Huber loss is "arguably one of the most widely used M -estimators" [15], owing its popularity to its convexity and differentiability properties of 2 , which allows for efficient algorithms (see, e.g., [41] for algorithms), in combination with its robustness properties of 1 [29]. This makes it widely applicable in practical big data settings (see, e.g., [5] 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., [38], [39] 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 √ 2 poly(log n, -1 )) queries to b, where 4 -2 √ 2 ≈ 1.17157. For this result, we use the chaining techniques of [7], which provides a more flexible alternative to [34], but requires more technical effort to adapt to the active setting.

Theorem I.5 (Main Result for Huber Active Regression). Let

There is an algorithm which, with probability at least 99/100, returns a x satisfying

Our techniques also yield a subspace embedding result, which constructs a weighted subset

Previously, the best known dimension reduction bound for Huber regression, even in the non-active setting, was d 4 due to [15].

Furthermore, this is, to the best of our knowledge, the first example of a loss function other than p which achieves a sensitivity sampling bound of better than d 2 , despite the fact that such results have been sought in many works [13]- [15], [28], [50], [52]. 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 p norms, the arguments of [7], [48] and their subsequent improvements avoid this problem by using a more sophisticated chaining argument. However, these arguments use the structure of p spaces in crucial ways, such as isometric changes of density using Lewis weights [32], 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.

f) Dimension Reduction for Gamma Functions for Faster p Regression.: One particularly important application of sampling-based dimension reduction for loss functions beyond p losses is, perhaps surprisingly, in the design of algorithms for p regression. The work of [8] introduces gamma functions γ 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 p regression. Subsequently, [2] obtained even faster algorithms by using constant factor approximations of γ p regression as a subroutine, in which the γ 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 p regression [1], [28]. In particular, [1] highlighted the open question of designing sparsification methods for γ p functions for p ∈ (1, 2), and [28] designed a sampling algorithm which samples Õ(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 √ 2 poly(log n, -1 )) rows for any p ∈ [1, 2), and improves to O(d poly(log n, -1 )) rows as p → 2 (see Figure 1 for the trade-off curve).

g) Kronecker Product Regression.: Beyond applications in data-efficient regression, Theorem I.2 implies the first sublinear time algorithm for Kronecker product regression in any p norm, where explicitly constructing the vector b is a computational bottleneck. In q-th order Kronecker product regression, one is given matrices A 1 , A 2 , . . . , A q , where A i ∈ R ni×di , as well as a vector b ∈ R n1n2•••nq , and the goal is to solve:

where ⊗ denotes the Kronecker product. Typically q i=1 d i is much less than q i=1 n i , and the goal is to obtain algorithms that do not explicitly form A 1 ⊗ A 2 ⊗ • • • ⊗ A q or b, which is too expensive. Our results yield the first algorithm for Kronecker product regression, for every p ≥ 1, whose running time does not depend on nnz(b), whereas previous results had a linear dependence on nnz(b), which can be as large as q i=1 n i [22]. Theorem I.6. Let q ≥ 1, p ≥ 1 be constant, and > 0.

Kronecker product regression can be solved up to a (1 +

)-factor with constant probability in Õ( 1) p Active Regression: Our algorithm for solving Problem I.1 uses a novel variation on the "sample-and-solve" approach. In particular, we randomly select a row sampling matrix S ∈ R m×n and return x = arg min x SAx -Sb , which only requires querying m entries of b (those that appear in Sb). To get tight bounds for p regression, we select S using p Lewis weight sampling, a generalization of leverage score sampling for 2 .

It can be shown that the p Lewis weights upper bound the p sensitivities of A, a measure of importance for the rows of A. The p sensitivity of the i th row of A is defined as

where [Ax](i) denotes the i th entry of the vector Ax, and captures how large the i th entry of any Ax ∈ span(A) can be, relative to the p norm. A standard scalar Bernstein bound shows that if S samples rows with probabilities that upper bound the sensitivities, then SAx p p = (1 ± ) Ax p p with high probability, for each x ∈ R d . An -net argument can extend this to a for all claim.

a) Prior Approaches to p Active Regression: While the above ideas give an approach for standard p regression, this bound does not suffice for active p regression. To solve Problem I.1, we actually want that S(Ax -b) p p = (1 ± ) Ax-b p p for any x. Will S provide such a guarantee? The main problem, as discussed in [10], [43] is that the translation by b may introduce outliers, i.e., entries with high sensitivity which are not captured by the sensitivity scores of A. As shown by [10], [43], in the case of 1 , the special structure of the loss function provides a solution. Indeed, by the triangle inequality,

where x * is the optimal solution. This fact can be used to show that sampling by the sensitivities of A preserves the differences between the cost of any x and the optimal x * . However, such a proof cannot work for p = 1, in which case we do not have such a nice inequality. p ∈ (1, 2), [10] 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 .

b) 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 x c such that Ax c -b p p ≤ O(1) • min x Ax -b p p using an idea of [20] and replace b by the residual vector b-Ax c , and then (2) conceptually partition the target vector b into two sets of coordinates, the coordinates i ∈ [n] that are small enough to be comparable to the sensitivity s p i (A) and those that are much larger. That is, we consider the coordinates i ∈

for some C > 0, and all other coordinates. For the former set of coordinates, one can check that the Bernstein bound still applies, and S does preserve the norm of Ax -b 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 Ax can both be close to b(i) in its i th entry, and still close to the remainder of b -the i th entry is simply too large in magnitude. In particular, to have [Ax](i) close to b(i), we would require Ax p p to be much larger than b p p , which by our preprocessing step, is on the order of the optimal cost min x Ax -b p . Via the triangle inequality, this implies that Ax must be far from an optimal solution. Thus, we can argue that any near-optimal solution to min x Ax -b p does not need to fit b(i) with |b(i)| p / b p p much larger than s p i (A). We can effectively ignore the contribution of these rows.

Another technical challenge remains: to obtain an optimal dimension dependence, we need a refined -net argument to make for all statements about x ∈ R d . To do so, we adapt the chaining arguments of Bourgain, Lindenstrauss, and Milman [7] and Ledoux and Talagrand [34] to the active regression setting, avoiding the standard -net and union bound argument used by, e.g., [47]. Although both [7] and [34] provide such approaches, we adapt the (slightly) more complex recursive Lewis weight sampling algorithm of [34] in order to obtain tighter dependencies on . The streamlined proof of [34] 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 [7] argument to handle the Huber loss, in which case the proof of [7] 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 b based on sensitivity weight classes. Aside from our new application of [7], [34], we hope that by translating the arguments of [7], [34] 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 that for any p > 1, any algorithm which simply samples reweighted rows and solves the system must have a polynomial dependence on 1/δ in the query complexity, while our algorithm achieves a log 1 δ dependence, by solving residual problems of a constant factor solution. Thus, our two-stage approach is necessary to achieve our δ dependence. Furthermore, the best known analysis of a simple "one-shot" Lewis weight sampling scheme suffers in dependencies for p > 2, where the one-shot approach is only known to give a Õ(d p/2 / 5 ) bound for subspace embeddings [7], [18], whose losses translate to losses for our active regression algorithms as well, while the recursive approach can achieve Õ(d p/2 / 2 ) [34]. While the [34] result is an existential result, we provide an analysis of the [34] proof to turn it into a randomized algorithm with logarithmic dependencies on the failure rate δ, which achieves the best known dependence on d, , and δ, up to logarithmic factors. We further modify this subspace embedding result for active regression, to optimize our dependence. c) Optimized Dependence for p ∈ (1, 2): For p ∈ (1, 2), the above argument gives a bound of Õ(d/ 2 ). While the linear dependence on d is optimal, confirming the conjecture of [10], it has a quadratic dependence on , which is in fact not optimal. We now show how to improve our bound to Õ(d/ ), which requires additional ideas. We first use strong convexity to show that a

Then, by using that x is close to the optimal solution, we show an improved bound on the difference in the objective values of x and x * , i.e., that x actually has an approximation ratio better than (1 + γ). We then iterate this argument until we obtain a (1+ )-approximation using Õ(d/ ) 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. This improves an p power means lower bound of [19], who only showed a lower bound of Ω( 1-p ) queries. Unfortunately, we are unable to port our algorithmic techniques to the p power means problem in high dimensions, due to difficulties in adapting their chaining argument.

2) Sensitivity Bounds: The notion of 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 × d matrix A. Then, the sensitivity of the ith coordinate with respect to the loss function M is defined as

It is well-established that sensitivities provide a general framework for sampling rows of A that approximate A well under the loss function M [25]. While a rich literature exists for p [18], [20], [49], little was known about the approximation of sensitivities for general loss functions until [52], 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 [52] raise the open question of obtaining general sensitivity bounds without this expensive subroutine.

Our approach to new sensitivity bounds significantly generalizes the approach of [13], 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 ∈ [n] has M -sensitivity α ∈ (0, 1], that is,

and let y = Ax witness this supremum, and assume for simplicity that n i=1 M ([Ax](j)) = 1. Note then that there can be at most 1/α entries j ∈ [n] of y that have coordinate value M (y j ) ≥ M (y i ) = α. Then, if we randomly hash the n coordinates into O(1/α) buckets, then with constant probability, coordinate i will be isolated from any other entry with M (y j ) ≥ M (y i ) = α. Now if M is monotonic, then this means that y 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 α fraction of the total M -mass, so entry i carries a constant fraction of the Mmass in its bucket. In this case, it can be shown that entry i must in fact carry a constant fraction of the 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 2 , then the largest entry will have the largest increase in its normalized contribution. This means that row i must have an 2 leverage score of Ω(1), in this hash bucket.

This leads to the following algorithm: (1) hash the n coordinates into O(1/α) buckets (2) compute 2 leverage scores for each bucket (3) assign an M -sensitivity of α for any coordinate that has 2 leverage score Ω(1). In each of the O(1/α) buckets, we will find at most O(d) coordinates with leverage score at least Ω(1), so we assign an M -sensitivity of α to at most O(d/α) coordinates, which has a total sensitivity contribution of O(d). By repeating this for O(log n) guesses of α 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. 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 I-C1. Furthermore, we demonstrate other applications of our sensitivity bound result, showing how to improve Orlicz norm subspace embeddings and robust subspace approximation.

3) Subspace Embeddings and Active Regression for the Huber Loss:

As discussed previously, we tackle the question of leveraging the theory of [7] 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 b 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 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 -net arguments is that there is no analogue of the chaining constructions of [7], [48] 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-ofdensity type theorem for the Huber loss as done in [36], [48]. We show how to overcome these obstacles in the following discussion.

a) A Sharp Huber Inequality.: Our algorithmic framework follows the Huber algorithm of [15], 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, [15] 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 1 and 2 norms:

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 1 and 2 sensitivities. This motivates the idea of sampling the rows of A with probability proportional to the sum of the 1 and 2 Lewis weights, oversampled by a factor of O(n 1/2 ). This is indeed how [15] proceeds.

The recursion n → 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 [7] net bounds, which reduces the row count in one recursive application to roughly n 1/2 d rather than n 1/2 d 2 . 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 [7] nets. The second improvement is that the Huber inequality as proved in [15] is in fact loose by a polynomial factor in n, and can be improved to the following: Lemma I.8 (Huber Inequality version 2). Let y ∈ R n . Then,

This lemma is tight up to constant factors 10 , and gives a recursion of roughly n → n 1/3 d, giving

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

b) Storing Large Huber Sensitivities.: In order to further improve upon this bound, we crucially make use of our improved sensitivity bounds and a generalized version of the above Huber inequality lemma that is parameterized by an upper bound on the size of the entries of y. Lemma I.9 (Huber Inequality version 3). Let y ∈ R n and let

Then, for some constant c > 0, at least one of the following bounds holds: 10 Consider the vector with one coordinate with n 1/3 and (n -1) coordinates with n -1/3 . By directly including the rows of A with Huber sensitivity at least γ, we exactly preserve the Huber norm inside T = [n] \ T for every y. 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/γ, and the sampling bound inside T , which is roughly (1/γ + (γn)

Then, for some c > 0 and β = 3 -2 √ 2 ≈ 0.17157, at least one of the following bounds holds:

In fact, we prove a generalized bound for the 2q loss for any q ∈ (0, 2) in Lemma I.11. The interval p ∈ [1,2] can be discretized in increments of 1 log n , so with O(log n) applications of [7] nets, we can always find a p within an additive 1 log n of the optimal p for every net vector y, which only affects Lemma I.10 by constant factors whenever y has entries bounded by poly(n). By proceeding as previously discussed, we arrive at our final bound of

d) Extensions to 2q Loss.: We generalize our results to the 2q loss for q ∈ (0, 2). As q ranges from 0 to 1 to 2, the 2q interpolates between the Tukey, Huber, and 2 losses up to constant factors, and provides a natural generalization of these loss function.

Lemma I.11 ( 2q Inequality). Let q ∈ (0, 2) and define

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

where β = 1 2/q -1 (2/q + 1) -2 2/q .

For q ∈ [1, 2), the 1/γ 2/q-1 distortion in Equation ( 2) is smaller than the 1/γ factor incurred from keeping Msensitivities at least γ, so we can balance the parameters as 1/γ = (γn) β as before, which leads to a recursion that gives us a bound of n = d 1+β poly((log n)/ ). For q ∈ (0, 1), the 1/γ 2/q-1 distortion in Equation ( 2) is worse than 1/γ, which means we must balance 1/γ 2/q-1 = (γn) β , or γ = n -β/(2/q-1+β) , which gives a worse bound of n = d γ poly((log n)/ ) for

This is better than a d 2 bound as long as q ≥ ( √ 5 -1) 2 /8 ≈ 0.19098.

Fig. 1: Dependence on d for the active regression sample complexity for the 2q loss. Similar bounds apply to subspace embeddings as well.

D. Conclusions and Future Directions

In this work, we study the sample complexity of active linear regression for both the p norm as well as general Mestimator losses.

For the p norm, we provide optimal algorithms and lower bounds for p ∈ (0, 2), with Θ(d/ 2 ) samples for p ∈ (0, 1) and Θ(d/ ) samples for p ∈ (1, 2). For p > 2, we provide an upper bound of Õ(d p/2 / p ), which is optimal in the d dependence and off by a single factor in the dependence, up to polylogarithmic factors. Our algorithms provide the first nontrivial bounds, i.e., sample complexity less than n, for p ∈ (0, 1) ∪ (2, ∞), while for p ∈ (1, 2), we significantly improve upon the Õ(d 2 / 2 ) upper bound of [10] and answer their main open question. We obtain these results via a twostage 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 dependence for p ∈ (1, 2). Our result is the first to achieve a linear dependence on for dimension reduction for p regression for p ∈ (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 [52] and makes significant progress in the general direction of matrix approximation beyond 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 [10].

For the important special case of the Huber loss, we introduce new techniques which bound Huber sensitivities by the sum of p Lewis weights, which allows us to take advantage of chaining arguments for p in order to obtain an active regression algorithm making at most O(d 4-2 √ 2 poly(log n, -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 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 [15].

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 [51] and robust subspace approximation [14]; our techniques for dimension reduction for the Huber loss gives improved bounds for sparsification for γ p functions for applications in fast algorithms for p regression [1], [28]. We believe that our techniques will be applicable 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 p regression for p > 2: our upper bound is Õ(d p/2 / p ), while the lower bound is Ω(d p/2 + 1-p ). Closing this gap would be interesting. Our bounds are also loose by a factor of log 1 δ for all p > 0, while we can get an optimal dependence on δ if we assume knowledge of the optimal value and sacrifice a factor of . A natural question if one can achieve a simultaneously optimal dependence on d, , and δ, 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 Õ(d 4-2 √ 2 poly log n) for constant , while only a trivial lower bound of Ω(d) is known.