Fast Regression for Structured Inputs
Authors: Raphael A Meyer
Cameron N Musco
Christopher P Musco
David Woodruff
Samson Zhou
Venue: ICLR 2022 Poster (2022)
Keywords: regression
sublinear time algorithm
structured input

Abstract: We study the $\ell_p$ regression problem, which requires finding $\mathbf{x}\in\mathbb R^{d}$ that minimizes $\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_p$ for a matrix $\mathbf{A}\in\mathbb R^{n \times d}$ and response vector $\mathbf{b}\in\mathbb R^{n}$. There has been recent interest in developing subsampling methods for this problem that can outperform standard techniques when $n$ is very large. However, all known subsampling approaches have run time that depends exponentially on $p$, typically, $d^{\mathcal{O}(p)}$, which can be prohibitively expensive. 

We improve on this work by showing that for a large class of common \emph{structured matrices}, such as combinations of low-rank matrices, sparse matrices, and Vandermonde matrices, there are subsampling based methods for $\ell_p$ regression that depend polynomially on $p$. For example, we give an algorithm for $\ell_p$ regression on Vandermonde matrices that runs in time $\mathcal{O}(n\log^3 n+(dp^2)^{0.5+\omega}\cdot\text{polylog}\,n)$, where $\omega$ is the exponent of matrix multiplication. The polynomial dependence on $p$ crucially allows our algorithms to extend naturally to efficient algorithms for $\ell_\infty$ regression, via approximation of $\ell_\infty$ by $\ell_{\mathcal{O}(\log n)}$. Of practical interest, we also develop a new subsampling algorithm for $\ell_p$ regression for arbitrary matrices, which is simpler than previous approaches for $p \ge 4$.

Published as a conference paper at ICLR 2022
FAST REGRESSION FOR STRUCTURED INPUTS
Raphael A. Meyer New York University tam900@nyu.c
Cameron Musco
Christopher Musco
University of Massachusetts Amherst causcolcs. umass. .dd
New York University
cmuscoênyu. edu
David P. Woodruff
Samson Zhou Camegie Mellon University samsonzhoul gma com
Camegie Mellon University
proodruf@andien cmu. edu
ABSTRACT
We study the lp regression problem, which requires finding x R that mini- mizes | - bllp for matrix A Rnxd and response vector b R". There has been recent interest in developing subsampling methods for this problem that outperform standard techniques when n very large. However, all known subsampling approaches have run time that depends exponentially on p. typically, do(e), which can be prohibitively expensive. We improve on this work by show- ing that for large class of common structured matrices, such as combinations low-rank matrices, sparse matrices, and Vandermonde matrices, there are sub- sampling based methods for l, regression that depend polynomially on p. For ex- ample, we give an algorithm for l regression on Vandermonde matrices that runs in time (n log" + (dp2)o.5++ polylogn). where w is the exponent of matrix multiplication. The polynomial dependence on p crucially allows our algorithms to extend naturally to efficient algorithms for lo regression, via approximation of lo by loolog n)- Of practical interest, we also develop new subsampling al- gorithm for lp regression for arbitrary matrices, which is simpler than previous approaches for 2 4.
INTRODUCTION
Given matrix A Raad and vector b E RR,, the goal of linear regression is find vector x E IRd such that Ax is as close as possible to b. In approximate l, linear regression in particular, we seek to find x E IRd such that, for some approximation parameter > O.
||Ax bl1p I (1 (1++1 E) min ||x bllp-
XER
Here for vector y E R", |ly E // Ep regression is central in statistical data analysis, and has numerous applications in machine learning, data science, and applied mathemat- ics (Friedman et al., 2001; Chatterjee & Hadi, 2006). There are a number of algorithmic approaches to solving fp regression. For example, we can directly apply iterative methods like gradient descent or stochastic gradient descent. Alternatively, we can use iteratively reweighted least squares, which reduces the regression problem to solving poly (d) linear systems (Adil et al., 2019b;a).
Both the above approaches require repeated passes over the matrix A, so while their runtimes are typically linear in nd,
generally on the time to multiply the matrix A by vector, that factor multiplied by other parameters, such as the number of iterations to convergence. An alternative approach, which
lead to faster running time when is large, to apply "sketch-and-solve" methods. This approaches begins with an inexpensive subsampling step, which selects a subset of rows in A to produce a smaller matrix M with poly (d, 1, log n) < rows. M can be written as M SA where S is row sampling and rescaling matrix. The goal is for ||SAx SbIIp to be good approximation A bllp for all Rd If this the then approximate solution to the original é, regression problem can be obtained by solving the subsampled problem, which has smaller size, thus allowing for more efficient computation.

---

Published as a conference paper at ICLR 2022
The standard approach to subsampling for lp regression is to sample rows with probability propor- tional to their so-called Ep Lewis weights (Cohen & Peng, 2015). Unfortunately. for general inputs A, Ep Lewis weight sampling requires O (dmax((_/22) rows, and can be shown that no subsam- pling method can take fewer than d0(p) eLi al., 2021). This means that sampling is only helpful in the limited regime where >>> dp/2. However, there are many applications in which the matrix A has additional structure, which can be leveraged to design more efficient algorithms. For exam- ple, Vandermonde matrices are used in the polynomial regression problem, which has been studied for over 200 years (Gergonne, 1974) and has applications to machine learning (Kalai et al., 2008), applied statistics (MacNeill, 1978), and computer graphics (Pratt, 1987). The goal is to fit signal, which is measured at time points f1.....ty, tn using a degree d polynomial. This problem can be for- mulated as lp regression with Vandermonde feature matrix A, whose 7t row a, is of the form [1, f, (t.))..... ...... (..) Regression problems with Vandermonde matrices also arise in the settings of Fourier-constrained function fitting (Avron et al., 2019) and Toeplitz covariance estimation (EI- dar al., 2020). Notably, (Shi & Woodruff, 2019) leverages the structure of Vandermonde matrices to more quickly build a subsampled matrix M given any Vandermonde matrix A than would be possible for general input. Their method does change how many rows in the matrix M, so the overall algorithm still incurs an exponential dependence in So while the approach is helpful improvement for Vandermonde regression where p small, this leaves an infeasible run- time for important problems like foo regression, which can be approximated by lp regression for p = O (logn).
1.1 OUR CONTRIBUTIONS
We first show that for lp regression on Vandermonde matrices, it is possible to reduce the size of the subsampled matrix M to depend polynomially instead of exponentially on p.
Theorem 1.1. Given E (0, 1) and p 2 Vandermonde matrix A E IRndd and b E R", there exists an algorithm that uses O (n logg n) + 00.5+ poly (2 P...oog n) time to compute a sampling matrix S E R. with O g'd log" vi) rows so that with high probability, for all x e Rdd
- [-=)|lAx bl1p I |SAx Sbll, S (1 +1+ bbbx-
blip,
and then to return vector & E Rd such that ||AR bl1, - (1 + 8) minXER ||Ax - bllp-
The best known previous work required at least Ωn log" n + d=/2 poly(1/s)) time (Shi & Woodruff, 2019). This has an exponential dependence in p. while our algorithm has just a poly- nomial dependence.
Building on Theorem 1.1, we observe that to obtain a (1 + apparoxximttoo to the fundamental problem of lo polynomial regression, it suffices consider l, regression for p O log
Since our results have polynomial dependence on p rather than exponential, we thus obtain the first subsampling guarantees for Vandermonde fox regression.
Theorem 1.2. Given e E (0, 1), Vandermonde matrix A E Ruxd, and E IR, there exists an algorithm that uses O (n log" n) + d005+w poly (,, logn) time to compute a sampling matrix
S Rmxn with O ( [ log" n) such that with high probability, for all x Rd.
- =)= blloo |SAx- - Sbil (1 e)|lAx blloc,
and then to return a vector & E Rdd such that |A bll $ (1 + E) minxeR* || - bll
To the best of our knowledge, this is the first known dimensionality reduction for Eo regression with provable guarantees for any (nontrivial) input matrix. We summarize these results in Table 1.
Our second contribution is to show that improved sampling bounds for f, regression can be extended to broad class of inputs, beyond Vandemonde matrices. We introduce the following definition to capture the "true dimension" of regression problems for structured input matrices.
Definition 1.3 (Rank of Regression Problem). Given an integer p 2 1 and matrix A E R suppose there exists matrix M E Rextt and a fixed function f Rdd - RR so that for all x E Rd.
( x) I = f fxx)
Then call the minimal such t the rank of the E, regression problem.
2

---

Published as a conference paper at ICLR 2022
Rows Sampled, too Regression Rows Sampled, Ep Regression d' poly (log n, 1) d² poly (log n,
Reference
n
(Avron et al., 2013)
(Shi & Woodruff, 2019) Our Results
dp² poly (log n, () (Theorem 1.1)
d dpoly (log n, !) (Theorem 1.2)
Table I: Sample complexity for regression on Vandermonde matrices
Theorems 1.1 and 1.2 rely on the following key structural property that we prove about the p-fold tensor product of rows of the Vandermonde matrix.
Lemma 1.4. For integer p 2 L, the rank of the Ep regression problem on Vandermonde matrix A E Inxdd is O (dp).
Lemma 1.4 implies that the f loss function for row of Vandermonde matrix can potentially be expressed as linear combination of O (dp) variables even though the entries of the measurement vector can be arbitrary. By comparison, the lp loss function for row a; of a general matrix A is (a, x) bp can only be expressed as linear combination of O (pdP) variables, correspond- ing to each of the d' k-wise products of coordinates of for each & [p]. As a corollary of Lemma 1.4, Theorems 1.1 and 1.2 obtain a small coreset for &, regression (as well as loe regression) on Vandermonde matrix, which can thus be used preconditioner for lp regression.
We generalize Lemma 1.4 to similar guarantees for lp regression on a matrix A that the sum of low-rank matrix K and an s-sparse matrix S that has at most non-zero entries per row. Thus using the notion of the rank of the regression problem for such a matrix A, we obtain the following guarantee:
Theorem 1,5. Given E E (0,1) and p 2 1, a rank k matrix K E Raxd and a s-sparse matrix E Raxd so that A == K + S, and b E RR" there exists an algorithm that uses (nk"~1) + n poly (20, "", sh, ,, logn) time to compute sampling matrix T E mmxx containing m = O pd ((++)) log" (pm) rows so that with high probability. for all x E Rd,
(1 -) bllp S TAxx Tblpp S (1 e)
and then to return a vector & E Ra such that || bllp S (1 e) minxeR= ||Ax blip- Fur- ther, if the low-rank factorization of K is given explicitly, then this runtime can be improved to poly (2P, d kP, s², ,, log
Similarly, we obtain efficient guarantees for Ep regression on a matrix A that is the sum of a Van- dermonde matrix V and sparse matrix S that has most non-zero entries per row.
Theorem 1.6. Given E E (0, 1) and p 2 1, a Vandermonde matrix E Rxd and an s-sparse matrix S E Raxd such that A
V + S. and b E R", there exists algorithm that uses O log" + +poly (p²d, d", sP, llg n, {)) time to compute sampling matrix T E Rmxn with m O 10 (170412"((tD)) log" (pm) rows, that with high probability for all x E Rd,,
11 S (1 + Ax - bllp- (1 -1-6)l e)= bllp S | Tbllp
and then to return a vector E Rdd such that |AA blIp S (1 + E) min ||Ax bllp-
Surprisingly, our methods even yield practical algorithms for general matrices with absolutely no structure. Although the rank of the lp regression problem for general matrices is O (dP), we still obtain the optimal sample complexity (see (Li al., 2021) for a lower bound) of roughly O (dr/2). Furthermore, an advantage of our approach is that we only need to perform Eg Lewis weight sam- pling for q [1.4] 4) and there are known efficient iterative methods for computing the ly Lewis weights for q [1, 4]. e.g., see Section 2.1 of this paper or Section 3 in (Cohen & Peng. 2015). By contrast, previous methods relied on computing general lp weights for p > 4, which, prior to the re- cent work of (Fazel et al., 2021), required solving a large convex program, e.g., through semidefinite programming, see Section 4 in (Cohen & Peng, 2015).
Finally, we experimentally validate theory synthetic data. In particular, we consider matri- ces and response vectors that are motivated by existing lower bounds for subsampling and sketch- ing methods. For Vandermonde regression, experiments demonstrate that the number of rows
3

---

Published as a conference paper at ICLR 2022
needed is polynomial in p. reinforcing how structured matrices outperform the worst-case bound. For unstructured matrix Ep regression, we demonstrate that our lg Lewis Weight subsampling scheme is effective and accurate.
1.2 OVERVIEW OF OUR TECHNIQUES
Our algorithmic contributions rely on two key observations, which we describe below. For the sake of presentation, we assume is an integer in this overview, though we handle arbitrary real values of p the subsequent algorithms and analyses.
Reduced rank of the Ep regression problem on structured inputs. The first main ingredient is the simple yet powerful observation that the rank of the lp regression problem on structured inputs such as Vandermonde matrices A E Rnxd does not need to be the O (d") rank of the ép regression problem on general matrices. For a general matrix, xa) for an integer can be rewritten as a linear combination Σ Outy; of O (d²) terms, where each coefficient o; is a p-wise product of entries of the row vector E IRd and similarly each is p-wise product of entries of the vector x E Rd. However, when A is Vandermonde matrix, then the j-th entry of a; is simply A/i²1. Thus each coefficient in (a., x) can be written as a linear combination of 1, A. (Ar2)²..... (A, )((d-)) and so we can express (a., (a. as linear combination of O (dp) variables rather than O (dP) variables.
We can similarly show that the rank of the Ep regression problem on rank k-matrix A is O (kP) by writing each row a, as a linear combination of basis vectors Vi...., Vk- Then, using the Hadamard Product mixed paonudctrruuttttct property, we can rewrite (a, as a linear combination of p-wise products of the variables x)..x...... (V), i..., a linear combination of O (kP) variables
rather than O (d)) variables. It follows that the rank of the lp regression problem on a matrix A whose rows have at most s non-zero entries is O sd) by noting that (1) there are (9) = O (d")
sparsity patterns and that (2) for a fixed sparsity pattern, the rank of the induced matrix is at most which induces a linear combination of O (s²) variables for the fp regression problem. These decomposition techniques can also be generalized to show that the rank of the l, regression problem is low matrices A such that A K V, A = K + S or A = V + S, where K is low-rank matrix, V is a Vandermonde matrix, and S is a sparse matrix.
However, we remark that although the observation that the rank of the Ep regression problem on structured inputs can be low, this itself does not yield an algorithm for lp regression. This is because the loss function is (a,, x) byl" rather than (a. x)" and there can be n possible different values of b; across all i E [n).
Rounding and truncating the measurement vector: novel algorithmic technique. Thus, the second main ingredient is manipulating the measurement vector b E R" so that the loss function can utilize the low-rank property of the lp regression problem on structured inputs. A natural approach would be to round the entries of b. say to the nearest power of (1 + e). Unfortunately, such a rounding approach would roughly preserve |[b|lp up to a multiplicative (1 + E) factor, but it would not preserve ||Ax- blpp- For example, suppose (a,, x) = b = N for some arbitrarily large value N. Then (a., x) - by 12 0, but if by were rounded to (1 + e) N, we would have (a,, x) b,l" = ep which can be arbitrarily large.
The lesson from this counterexample that when b, is significantly larger than (a,, x) b, any rounding technique can be arbitrarily bad because ||b||p can be significantly larger than OPT minxeR" ||Ax bllp- On the other hand, if ||bl|p a constant factor approximation to OPT, then the previous counterexample cannot happen because either (1) is large relative to ||b|p and (a, x) cannot be too close to b, so that rounding 6, will not significantly affect the difference (a, - b, or (2) b, is small relative to ||bllp and any rounding of b, will not affect the contribution of the i-th row of A in the overall loss. Thus our task is reduced to manipulating the input so that |lb||p O (OPT).
To that end, we note that if || argument, we have that the residual vector b' = b- Aă Rd satisfies [|bl)p O (OPT). Namely, we show that to find such a vector x, it suffices by triangle inequality to find a subspace embedding for A, i.e., matrix M E R with < n such that
bllp (OPT) for vector x Rd, then by triangle inequality
(1 A) < | $ (1 Cx]]
4

---

Published as a conference paper at ICLR 2022
for some constant C E (0, 1). Typically, such matrix M can be found by sampling the rows of A according to their Ep Lewis weights to generate matrix with m O (dmax(1.p/2)) rows. However, using our observation that the rank of the Ep regression problem on structured inputs can low, sample O (dp) rows of A with probabilities proportional to their lg Lewis weights, where q- = for integer such that 2 S < 2 Crucially, we can find such vector X without reading the coordinates b since we only require to read the rows of A to perform Lewis weight sampling in this phase. Thus can efficiently find residual vector b' such that ||b'lpp O (OPT).
Partitioning the matrix into groups. We then round the entries of to obtain a vector b" with most 8= 0 log unique values, truncating all entries that are less in magnitude than poly(n) c
to zero instead. We now solve the lp regression problem min || b"llp by partitioning the rows of A into l groups GI...., Gg based on their corresponding values of b" The main point is that all rows in each group G3 have the same entry lg in b". Thus can again observe that (a., x)
- 4 can be written as linear combination of O (p"d) variables and therefore use lg Lewis
weight sampling (for the same q) to reduce each group G₁. down to O (r'd log d) rows. Since there
groups, then there are roughly O (p²d) total rows that have been sampled across all the groups It follows from the decomposition of the ép loss function across each group that the matrix T formed by these rows is a coreset for the lp regression problem. Therefore, by solving the lp regression problem on T, which has significantly smaller dimension than the original input matrix A, we obtain vector X with the desired property that
are f 0 log
AR - blbp S (s)) min | XER
b11p- 4Ax-bllp-
Practical b regression for arbitrary matrices. We now describe a practical procedure for ép regression general matrices that avoids the necessity for convex programming to approximate the E, Lewis weights for p > 4. We first pick an integer r such that P 27 E [2.4). By the same structural argument as in Lemma 1.4, we can tensor product each row a; of A with itself 2" times, thus obtaining an extended matrix of size n x d', where d' = d² independent of the entries in b. We then fp/2- Lewis weight sample on the extended matrix, using the iterative method in Figure since 2 < 4 rather than solving convex program for bp Lewis weight sampling for p > 4. Since P 2
E [2, 4), it follows from Theorem 2.3 that Lewis weight sampling requires roughly (d') rows. Thus we obtain a matrix with O ((dP") = O (dm/2) rows, from which we can compute a residual vector b' such that ||b'll, = (OPT), where OPT == minxeR= ||Ax bllp-
We then round and truncate the entries of b' using the subroutine RoundTrunc to obtain vector b", which allows us to partition the rows of A and the entries of b" into loga
groups. Due to the constant values of b" in each group, we can again 6p/2- Lewis weight sample on an extended matrix using an iterative method and finally solve the lp regression problem on the subsequent rows that have been sampled. We give our algorithm in full in Algorithm 3.
2
Lp REGRESSION FOR VANDERMONDE MATRICES
2.1 PRELIMINARIES ON LEWIS WEIGHT SAMPLING
There are number of known sampling distributions for dimensionality reduction for the E, regres- sion problem, such as the L, leverage scores, e.g., (Cormode et al., 2018; Dasgupta et al., 2008) and
the fp sensitivities, e.g., (Clarkson et al., 2019; Braverman et al., 2020; 2021; Musco et al., 2021); I); in this paper we focus on the &p Lewis weights, e.g., (Cohen & Peng, 2015; Durfee et al., 2018; Chen & Derezinski, 2021; Parulekar et al., 2021).
Definition 2.1 (C, Lewis Weights, (Cohen & Peng, 2015)). Given a matrix A RRaxd and p 2 1, the Ep Lewis weights wi( (A), (A),,.,,,,, (A) are the unique quantities that satisfy
(w,(A))2/p a, WAT--12/ A) a; for all i E [n), where W E Raxx is the diagonal matrix with W = w.,(A) for all i [m].
5

---

Published as a conference paper at ICLR 2022
Definition 2.2 (Lewis Weight Sampling). For an input matrix A B and samples, let the sampling matrix S E Rmnn be generated by independently setting each row of to be the i-th 1/p with probability (A)
standard basis vector multiplied by (m))((A)
Theorem 2.3 (L, Subspace Embedding from Lewis Weight Sampling, (Cohen & Peng, 2015)). For any A Raxd, let the matrix S Rmxx be generated from Lewis weight sampling with m =0 (u0as(1.//2) log(d/8) log(2/e) ) where 2 for E [1.2] and ? 5 for p > 2. Then with
probability at least - 6, we have that simultaneously for all x E Rd4
(1 -1-88 - < |SAx|| S (1 Ax|
Since the Lewis weights are implicitly defined, may not be clear how to compute them exactly. In fact, it suffices to compute constant factor approximations 2015) show that for < 4, there exists a simple iterative approach to compute a constant factor approximation the l Lewis weights in input sparsity time, which we present Figure 1.
the Lewis weights. (Cohen & Peng,
w = Lewislterate (A, p. w)
(1) For 1 I,...,
(a) Let 73 be constant factor approximation to a. A1
(b) f (T,)
(c) Return
w Approx ,p, 3.3,7 ()
(1) Intialize wy = for all i E [n]]
(2) For = 1, I.....TT
(a) w < S.wistteratteA.p.... B.w)
Return
Fig. I: Iterative algorithm for approximate &p Lewis weights with p 4.
At a high level, the correctness of Figure 1 follows from Banach's fixed point theorem and the fact that the subroutine Lewislterate a contraction mapping, because 22/p 1| < for p < (Cohen & Peng, 2015). However, for p 2 4, Figure no longer works. Instead, prior to the recent work of (Fazel et al., 2021), approximating é, Lewis weights for >> seems to require solving the convex program
Q = argmax det (M),
subject to Σ(a Ma)/2 <d d,
= 0.
M
and setting w; (a. Q/) Unfortunately, this often infeasible in practice, and we could not obtain empirical results using it. Therefore, nice advantage of our algorithms, both for structured and unstructured matrices, that we only use Eg Lewis weight sampling for q I [1, 4), even p 2 whereas previous algorithms required using f, Lewis weight sampling for 4.
2.2 ALGORITHM AND ANALYSIS
We first describe the general framework of our algorithm for efficient lp regression, so that given A E Raad and b Rd the goal is to approximately compute OPT == min Ax bllp- Recall that in order to apply our structural results, we first require a measurement vector b" with small number of distinct entries. We obtain b" by first finding a constant factor approximation. ie., x E Rd such that | blI < 0 (OPT). We can E, Lewis weight sample the rows of to do this, but the time to solve the subsequent polynomial regression problem would have exponential dependence on p. due to Theorem 2.3. Instead, we use our structural properties to implicitly create matrix M from A with fewer than d² columns and then f,/2- Lewis weight sample the rows of

---

Published as a conference paper at ICLR 2022
A. where is the integer that satisfies 2" p 2 and then we solve the ép regression problem on the sampled rows of A and entries of to obtain X. We set b' to be the residual vector Ax so that ||b'l I O (OPT). We then use the procedure RoundTrunc, i.c., as in Algorithm which sets the entries of b' that are the smallest magnitude to zero, and rounds the remaining entries of b' the nearest power of (1 E).
Algorithm RoundTrunee Round and truncate coordinates of input vector
Input: Vector b E IR",, accuracy parameter e > 0
Output: Vector
R that is rounded and truncated version of b
1: M f maXign bb:
2: for 1 to i do
3: Ps f |log [b,1]
4: I; + (1 + = sign(b,)
5:
if' [ S no
then
6:
7:
return
Since ||b'lpp = O (OPT). it then follows by triangle inequality that it suffices to approximately solve the Lp regression problem the vector b" instead. We partition the rows of A into groups GI...., Ge based on the values of b" and perform lp/2- Lewis weight sampling again on an implicit matrix that create from the rows of each group. Here we leverage our theory about the rank of the b regression problem for structured inputs. It then suffices to solve the En regression problem on the matrix formed by the sampled rows across all the groups
Algorithm Faster regression for Vandermonde matrices
Input: Vandermonde matrix A Raxd, measurement vector b, accuracy parameter Output: & E IR with | - bllp s (1 + O (E)) minxeR ||Ax bll, bllp
1: r f- [log p_
22" S <p< 2rt1
2: Extend A to Vandermonde matrix M with dimensions n x d', where d' = 1)+
3: Use fp/2--Lewis weight sampling on M to find a set S O (d' log d') indices [n] =
{1,2....,,1} and corresponding rescaling factors.
4: Let A' be the corresponding submatrix of A with indices in S, and scaled accordingly.
5: Compute < min ||A'x bllp-
6: b' + A b" Round TruncTruunc(1 Trunc(b'), - O Log
7: Partition the rows of A into groups G Ge, each containing all rows with the same value
of b"
8: Let Gg be the submatrix of G1 extended d" d= 22r 1) columns.
9: Let 4 be the coordinate of b" corresponding to G for each k E [4]
10: Use fp/2--Lewis weight sampling on [Gu: tx] to find set Sy of O [ log d") indices [n]
and rescaling factors.
11: Let Ta be the corresponding sampling and rescaling matrix for Sy
12: T : :T:
Compute < (1 + e) min TAx Tb/
14: return x
We show in the supplementary material that Algorithm satisfies the guarantees of Theorem 1.1. Moreover, because Theorem 1.1 has polynomial dependence on rather than exponential depen- dence, we obtain the first known sublinear size coreset for the important problem of fox regression. We use the following structural property.
Jogm) Then ||x||>0 S | S + +))
Lemma 2.4. Let E IR" and
Lemma 2.4 implies that to solve Exx regression, we can instead solve En regression for p
log m
Then Theorem 1.2 follows from the fact that even for Ωt logn the matrix T

---

Published as a conference paper at ICLR 2022
in Algorithm 2 satisfies
- []= Tb I =||Ax blI, I +1+ Tb"l,,
for all E RR Hence, can solve the lp regression problem the smaller matrix T to solve the 3o regression on A.
The results of Theorem 1.1 can be further extended matrices with block Vandermonde structure. Corollary 2.5. Given E (0,1) and p 2 E TRxxdy and TKé suppose A =AAll... [A,) for Vandermonde matrices A1...., Ag Raxd, Then there exists an algorithm that, with high probability, returns a vector & R. such that
|A - blpp S $1 (1+
min ||x - bllns XEE
using O (T()) log poly((dp)) logn, 1/E)) time, where T(A) the runtime of mul- tiplying the malrix A by arbitrary vector. For = 2, this can be further optimized to O (nd=2/2-1 polyy(dpp) log 1/e)) time, where (n²*) is the time multiply an x n matrix with an n x n² matrix, so that W2 E [3,4]
We obtain similar algorithms for noisy low-rank matrices and noisy Vandermonde matrices, i.e., Theorem 1.5 and Theorem 1.6. We defer presentation and justification for these results to the sup- plementary material.
Algorithm 3 Faster Ep regression for general matrices
Input: Matrix A E Raxd measurement vector b, accuracy parameter & > 0
Output: & E IR with ||A bll, - (1 E) minxers ||Ax bllp. parameter p 2 4
1: <- [logp] -
2 S p +
2: Extend A to matrix M with dimension
O (df) that each in M, is the fold tensor product of a; reshaped into d2 length vector.
3: Use fl/2i-Lewis weight sampling M to find set S of O (12/2 log d) indices in [n] and rescaling factors.
4: Let A' be the corresponding submatrix of A with indices in S and scaled accordingly.
5: Compute x < 5 m AA'x bllp-
6: b' + b - Ax, b" + Round Trunc (buunccbb, O log
7: Partition the rows of A into groups GI...., Ge. each containing all rows with the same value of b"
8: Let Gy be the corresponding submatrix and tę be the coordinate of corresponding to G2 for each A [!].
9: Use fp/2--Lewis weight sampling on [Gu; 4x] to find a set Sy of O ( and rescaling factors.
( d/2 log d) indices in [n]
10: Let Tx be the corresponding rows with indices in Sy and scaled accordingly.
11: T (Tii.. [T]
12: Compute < (1 + E) minxeRs Txx Tbl
13:
return &
Finally, we describe in Algorithm 3 a practical approach for tp regression for arbitrary matrices without requiring any structural assumptions.
Theorem 2.6. Algorithm outputs a vector & such that |AA bl1p 5 (1 + +) min | bllp- The runtime of Algorithu T Tmm, dd/2) poly (dm/2 logn, n, 1). where T(n, 22/2) is the time to multiply matrix of size x dol2 by vector of length do/2.
Theorem 2.6 achieves the same optimal sample complexity as previous É, Lewis weight algorithms of roughly / However, the main advantage of Algorithm 3 is that it performs Lewis weight sampling for some q E 4), which is quite efficient because we can use an iterative method rather than solving convex program.
8

---

Published as a conference paper at ICLR 2022
10
Layarr Lewis Weights
Relative Ervor (Comptrioa))
10
Delative Error (Compiricat)
Uuilorm
109
10a
107
10-1
10a
2p/2
Lewie Weights
10-3
Uniform
10a
15
500
1,000
1,500
2,000
2,500
3,000
Norm Parameter (p)
Subsample Complexiiy (m)
(a) Relative Error versus p
(b) Relative Error versus
Fig. Comparison relative error (Cempirical) versus ép parameter p and sample complexity m on Vandermonde data. We ran both experiments 30 times and plot the median, 25th quartile, and 75th quartile for each value of and m.
EMPIRICAL VERIFICATION
We provide empirical evidence validate core statistical claim about Vandermonde regression: that the relative eгтоr ε achieved by fp/2" Lewis Weight subsampling is polynomial in p, and not exponential in p. More precisely, we compute the gap between the error achieved from exact lp Vandermonde regression and the error achieved from the subsampled regression:
VX
1p
|| blpp
Cempirical
Vx - blIp
Our theory tells us that (Eempirical) 3
dp"
where m is the total number of subsampled
rows. The prior work unstructured matrices instead suggests (Cempirical) < Ö (
(Shi &
Woodruff, 2019). So, to visually distinguish these two settings, we look at the logarithm of both sides:
log(Eempiricat) s O (In(p) + In( 2
OR log(Eempiricat) S O (pIn(d) - In
In particular, work suggests logarithmic dependence on p, while the prior work suggests
linear dependence on p.
To validate our theory, we plot log(Eempiricat) versus p and m on synthetic data. Specifically, generate n 25, 000 i.i.d. N NO, 1) times samples to form a Vandermonde matrix with = 20 columns, then compute the polynomial qt) 10 each time sample, add N 10a") additive noise, and save the corresponding values in b. We then compute Eempérical for this l, regression problem. Notably, in order to compute x, we omit the rounding procedure in our code, since the rounding is designed for worst-case inputs. Instead, we simply compute & by solving ming // | bl1p where V and b are computed by sampling and rescaling V and with the p//2- Lewis Weights.
Figure shows the result of these experiments, which were run in Julia 1.6.1, on Windows 10 with an Intel 17-7700K CPU and 16Gb RAM. In Figure 2a, we fix m = 1000 and vary E [2, 25). The trendline of Lewis Weight sampling clearly better fits a logarithmic model, as opposed linear model. This reinforces our analysis by showing that the dependence on p notably sub-exponential, beating the known bounds for Ep subsampling on unstructured matrices.
As benchmark, we compare our Lewis Weight sampling method to uniform sampling. The noise in b is large enough that most rows A have little information about the underlying polynomial q(t). Lewis weight sampling takes avoids these rows, while uniform sampling does not, explaining why uniform sampling is much weaker in Figure 2a. Further, in Figure 2b, we fix p = 6 and vary E [100, 3000), showing that Lewis Weight sampling outperforms uniform sampling across both m and p. In Appendix C, we demonstrate similar results for unstructured matrix regression, validating the analysis of Theorem 2.6.

---

Published as a conference paper at ICLR 2022
ACKNOWLEDGEMENTS
Cameron Musco was supported by NSF grants 2046233 and 1763618, and an Adobe Research grant. David P. Woodruff and Samson Zhou were supported by National Institute of Health grant 5401 HG 10798-2 and a Simons Investigator Award. Christopher Musco and Raphael Meyer were supported by NSF grant 2045590 and DOE Award DE-SC0022266.
REFERENCES
Deeksha Adil, Rasmus Kyng, Richard Peng, and Sushant Sachdeva. Iterative refinement for 1p- norm regression. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 1405-1424, 2019a. 1,
Deeksha Adil, Richard Peng, and Sushant Sachdeva. Fast, provably convergent IRLS algorithm for p-norm linear regression. In Advances in Neural Information Processing Systems 32, NeurIPS,
pp. 14166-14177, 2019b. 1, 15
Haim Avron, Vikas Sindhwani, and David P. Woodruff. Sketching structured matrices for faster nonlinear regression. In 27th Annual Conference on Neural Information Processing Systems. Proceedings, pp. 2994-3002, 2013.
Haim Avron, Michael Kapralov, Cameron Musco, Christopher Musco, Ameya Velingker, and Amir Zandieh. A universal sampling method for reconstructing signals with simple fourier transforms. In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing. STOC, pp.
1051-1063, 2019. 2
Vladimir Braverman, Petros Drineas, Cameron Musco, Christopher Musco, Jalaj Upadhyay, David P. Woodruff, and Samson Zhou. Near optimal linear algebra in the online and sliding window models. In 6/st IEEE Annual Symposium Foundations of Computer Science, FOCS, pp. 517-528, 2020. 5,
Vladimir Braverman, Avinatan Hassidim, Yossi Matias, Mariano Schain, Sandeep Silwal, and Sam- son Zhou. Adversarial robustness of streaming algorithms through importance sampling. CoRR, abs/2106. 14952, 2021. 5, 22
Samprit Chatterjee and Ali S Hadi. Regression analysis by example, volume 607. John Wiley and Sons, 2006.
Xue Chen and Michal Derezinski. Query complexity of least absolute deviation regression via robust uniform convergence. In Conference on Learning Theory, COLT, pp. 1144-1179, 2021. 5, 22
Rachit Chhaya, Jayesh Choudhari, Anirban Dasgupta, and Supratim Shit. Streaming coresets for symmetric tensor factorization. In Proceedings of the 37th International Conference Machine Learning, ICML, pp. 1855-1865, 2020. 22
Kenneth L. Clarkson and David Woodruff. Low rank approximation and regression in input sparsity time. In Symposium on Theory of Computing Conference, STOC, pp. 81-90, 2013. 22
Kenneth L. Clarkson, Petros Drineas, Malik Magdon-Ismail, Michael W. Mahoney, Xiangrui Meng, and David P. Woodruff. The fast cauchy transform and faster robust linear regression. SIAM J.
Comput., 45(3):763-88) 810, 2016. 22
Kenneth L. Clarkson, Ruosong Wang, and David P. Woodruff. Dimensionality reduction for tukey regression. In Proceedings of the 36th International Conference on Machine Learning, ICML, pp. 1262 1261, 2019. 5,
Michael B. Cohen and Richard Peng. Ip row sampling by lewis weights. In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC, pp. 183-192, 2015. 2, 15, 21, 22
Graham Cormode, Charlie Dickens, and David P. Woodruff. Leveraging well-conditioned bases:
Streaming and distributed summaries minkowski p-norms. Proceedings of the 35th Interna-
fional Conference on Machine Learning, ICML. pp. 1058-1056, 2018. 5, 22
10

---

Published as a conference paper at ICLR 2022
Anirban Dasgupta, Petros Drineas, Boulos Harb, Ravi Kumar, and Michael W. Mahoney. Sampling algorithms and coresets for lp regression. In Proceedings of the Nineteenth Annual ACM-SIAM Symposium Discrete Algorithms, SODA, pp. 932-9411 2008. 5,
David Durfee, Kevin A. Lai, and Saurabh Sawlani. 1 regression using lewis weights preconditioning and stochastic gradient descent. In Conference On Learning Theory, COLT, pp. 1626-1656, 2018.
5, 22
Yonina C. Eldar, Jerry Li, Cameron Musco, and Christopher Musco. Sample efficient toeplitz co- variance estimation, In Proceedings of the 2020 ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 378 1778-997.00 2020. 2
Maryam Fazel, Yin Tat Lee, Swati Padmanabhan, and Aaron Sidford. Computing lewis weights to high precision. CoRR, abs/21 10.15533, 2021. 3, 6
Jerome Friedman, Trevor Hastie, Robert Tibshirani, et al. The elements of statistical learning. Springer series in statistics New York, 2001.
JD Gergonne. The application of the method of least squares the interpolation of sequences. Historia Mathematica, 1(4):439-44) 447, 1974. 2
I Gohberg and V Olshevsky. Complexity of multiplication with vectors for structured matrices. Linear Algebra and Its Applications, 202:163-192, 1994. 15
Jerrad Hampton and Alireza Doostan. Coherence motivated sampling and convergence analysis of least squares polynomial chaos regression. Computer Methods in Applied Mechanics and Engineering, 290:73-97, 2015. 20
Adam Tauman Kalai, Adam R Klivans, Yishay Mansour, and Rocco A Servedio. Agnostically learning halfspaces. SIAM Journal on Computing, 3766)1177--1805, 2008. 2, 20
Daniel Kane, Sushrut Karmalkar, and Eric Price. Robust polynomial regression up to the information theoretic limit. In 58th IEEE Annual Symposium Foundations of Computer Science, FOCS, pp. 391-402, 2017. 20, 21
Kiran S. Kedlaya and Christopher Umans. Fast polynomial factorization and modular composition.
SIAM J. Comput., 4066) 2011. 17
Yi Li, Ruosong Wang, and David P. Woodruff. Tight bounds for the subspace sketch problem with applications. SIAM J. Comput., 1874)31 2021. 2, 22
lan B MacNeill. Properties of sequences of partial sums of polynomial regression residuals with applications to tests for change of regression unknown times. The Annals of Statistics, pp- 422-433, 1978. 2
Xiangrui Meng and Michael W. Mahoney. Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression. In Symposium Theory of Computing Confer- ence, STOC, pp. 91-100. ACM, 2013. 22
Raphael Meyer, Cameron Musco, Christopher Musco, David Woodruff, and Samson Zhou. Cheyb- shev sampling is universal for 1p polynomial regression, 2021. 20, 21, 22
Cameron Musco, Christopher Musco, David P. Woodruff, and Taisuke Yasuda. Active sampling for linear regression beyond the l2 norm. CoRR, abs/21 11.04888, 2021. 5, 2121.22 22
Jelani Nelson and Huy L. Nguyen. OSNAP: faster numerical linear algebra algorithms via sparser subspace embeddings. In 54th Annual IEEE Symposium Foundations of Computer Science, FOCS, pp. 117-126, 2013. 22
Michael Nüsken and Martin Ziegler. Fast multipoint evaluation of bivariate polynomials. In Algo- rithms ESA, 12th Annual European Symposium, Proceedings, pp. 544-555, 2004. 17
Aditya Parulekar, Advait Parulekar, and Eric Price. I regression with lewis weights subsampling. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, APPROXIANNOM, pp. 49:1499:999 149:21,, 2021. 5,

---

Published as a conference paper at ICLR 2022
Vaughan Pratt. Direct least- squares fitting of algebraic surfaces. ACM SIGGRAPH computer graph- ics, 21(4):145-133 1987. 2, 20
Christopher De Sa, Albert Gu, Rohan Puttagunta, Christopher Ré, and Atri Rudra. A two-pronged progress in structured dense matrix vector multiplication. In Proceedings of the Twenty-Ninth
Anmual ACM-SIAM Symposium on Discrete Algorithms, SODA, pp. 1060-1079, 2018. 17
Tamás Sarlós. Improved approximation algorithms for large matrices via random projections. In 47th Annual IEEE Symposium Foundations of Computer Science (FOCS), Proceedings, pp- 153-152. IEEE Computer Society, 2006, 22
Xiaofei Shi and David P. Woodruff. Sublinear time numerical linear algebra for structured matrices. In The Third AAAI Conference Artificial Intelligence, pp. 4918-4925, 2019. 2, 9, 13
Christian Sohler and David P. Woodruff. Subspace embeddings for the 1 -norr with applications. In Proceedings of the 43rd ACM Symposium on Theory of Computing, STOC, pp. 755-764, 201 22
Ruosong Wang and David P. Woodruff. Tight bounds for lp oblivious subspace embeddings. In Proceedings of the Thirtieth Annual ACM-SIAN Symposium on Discrete Algorithms, SODA, pp- 1825-1843, 2019, 22
David P. Woodruff and Qin Zhang. Subspace embeddings and lp-regression using exponential ran- dom variables. In COLT 2013 The 26th Annual Conference on Learning Theory, pp. 546-567, 2013. 22
BROADER IMPACT
This work is focused improving the runtime of algorithms for lp regression structured inputs, and can allow users to more efficiently solve regression problems on large datasets. The problem we study is abstracted away from specific ethical and unethical applications of Ep regression. So, with the assumption that the research and data science communities will focus on ethical applications while avoiding unethical ones, we believe that the benefits of algorithmic results outweigh the negative risks.
A MISSING PROOFS FROM SECTION 2
In this section, we give the missing proofs from Section 2.
A. t, REGESSION ON VANDERMONDE MATRICES
We first prove simple statement that shows "good" solution to the lp regression problem subspace embedding SA is also a "good" solution to the original input matrix A. Lemma A.1. Let S be sampling and rescaling matrix such that
11 SAx|, - |||x|| <
13
SAXI
E |SAX = ||Ax|
12
12
for all IRd. Let OPT minxerd ||Ax blIp and let IRd be any vector for which ||SAx Sblpp < 5OPT. Then with probability at least 0.79,
Ax blIp I 12OPT.
Proof. Let x be a minimizer of ||Ax b/Ip so that OPT = ||A - blIp- We will prove the claim by contrapositive, we first suppose that |A b/l, 12OPT. Then by the triangle inequality,
SAA Sb/Ip |SA(x x*) ) - |SAA Sb|p.-
Since
11 |SAx|1 12
|Ax|| < 12 SAx|| for all Rd then
SAA Sb/p 2 12
A - x' - ||AAA Sbl1p-
12

---

Published as a conference paper at ICLR 2022
By the triangle inequality,
11
|| |AAx SbI1p
12 (l bblp ||Ax' - bblpp) |SAA -Sbllz-
Note that by Jensen's inequality, SX Sb|lp] S OPT, that by Markov's inequality.
1
Pr Pr[|SAx Sb/pp 2 50PT - 5
Thus with probability at least 0.8,
11
SAA Sbl1
(|) - A ||xx bblp)) -5PP..
12
Thus if ||Ax blIp S 12OPT, then
11 (12OPT - 50PT > 50PT,
SAAx Sblp 2
12
as desired.
Before justifying the correctness of Algorithm 2, we first recall the following algorithm for efficient bp subspace embeddings.
Theorem A.2. (Shi & Woodruff, 2019) Given E (0, and p E [1,2], Vandermonde ma- trix A E Raxe and E RRd lei T(A) be the time takes to perform matrix-vector mul- tiplication, ie., compute Av for arbitrary v E IR.. There exists an algorithm that uses O (A) + d9 pool(1/E)) time, where q = w for p E [1, 4) and q = p/2 + for some fixed constant C > 0 p 4, and with high probability, returns x E Rd such that
||A - bl1p I (1 + E) min XERd
|Ax- - bllg-
We now justify the correctness of Algorithm 2.
Lemma A.3. Given (0. 1) and p 2 Vandermonde matrix A Rexd, and b E Rdd then with high probability, Algorithm returns a vector & E R$ such that
A b/Ip - (1 + O (s)) min ||x blIp-
xER
Proof. Consider Algorithm and let be integer that 2" Sp < 2
Note that
Σ
(a x) = Σ
(( x) 2
tE(m)
iE(n)
Since A is Vandermonde, then
2"
2"
d 1
(dd-1))+1
((a, x))="
Σ
Σ
d Uj
Cijj
where each yy is a fixed function of the coordinates of Notably, the fixed function is the same across all E [n]. Hence, the far subspace embedding problem on an input Vandermonde matrix A E n can be reshaped as a constrained l subspace embedding problem on an input Van- dermonde matrix of size n x (2 (d 1) + 1). Thus, for lp regression with p E [2",20 12",2++1), we have
2 (d 1)+1
D/2
Σ (a,,)) x)
Σ
Vj
ieIn]
v mni
which is a constrained 2//2" regression problem on an input Vandermonde matrix M of size n x (2" (d-1) 1).
By Theorem A.2, we can use fp/2- Lewis weight sampling to find a matrix M' such that
My 1/2
< MMyl
p/2"
p/2
<
My
p/2" p/2
13

---

Published as a conference paper at ICLR 2022
for all y (4-11)+ with high probability. Note that by the above argument, if we take the matrix A' corresponding to the scaled rows A that are sampled by M', then we also have
11
13
Ax 1
A
S
Ax Ax
12
12
and thus
11
13
| 1 A X I
Ax|p
12
12
for all x E Rd+1 with high probability. Thus by Lemma A. we can find a vector such that
Ax blIp - 12 min Ax
bl[p-
xER
We use the subroutine Rrund to create the vector b" which is the vector with all entries of b' rounded to the nearest power of (1 + E), starting at the maximum entry b' absolute value, and stopping after we are poly(nj times that, and replacing all remaining entries with O. By the triangle inequality, have
|A - blll < Ax b'll, + Ilb' - bbll, < |AA blll + 12:OPT S (1+ 12:)] Ax - -b"llp.
for any x E IRd+1
Note that b" has discretized the values of b' into o log possible values. We partition the of A into groups Gi...., Ge, based the corresponding values of b". Suppose that for group Gk, the corresponding values of b" are all 44. Then we have
AA b"I
ΣE
Σ iEGe
K x) telp.
Since A is Vandermonde, for each E Gk.
d
2"
d
2"
2d-
((a,,)) - tx)²² = -t
tk
1.2
23
=
Σ
1=1
j=1
for some fixed values tx,1 ... that be computed from tk. where again each Yj be different function of the coordinates of x. In particular, there can be 2" choices for the exponent of
4. For a fixed choice of the exponent of lx. the exponent of az2 can range from 0 to 2" (d 1).
Hence, the for regression problem submatrix of input Vandermonde matrix A Rn with the same fixed measurement values, i.e., the corresponding coordinates of b" are all the same, can be reshaped to constrained regression problem
input Vandermonde matrix with 2 (d- +1 columns times (22r (d 1) + 1) x (2 (d 1) + 1) diagonal matrix.
Hence, by invoking Theorem A.2 sample rows of M corresponding their f//2- Lewis weights in the submatrix induced by the rows Gg, we obtain a sampling matrix Tg such that with high probability.
D/2"
(d-11)++1
(1-5) T |
P/2"
Σ EEG2
Σ Ø=1
S (.
T&v 1/2"
p/2"
0,2
where Vy corresponds to the vector b" that set to zero outside of coordinates whose values are f.- Note that by the above argument, then also have
(1 - TaAx T& S
Σ (a, - tel = (1+ Tave
REG.
with high probability. Summing over all k, we have
(-=) Σ
Σ
(a,, - tx <1++) Σ Ax Tavell
k
IEG2
with high probability.
14

---

Published as a conference paper at ICLR 2022
Let T be the sampling matrix that T = T1 .... Tr, that Σx T&Ax TaVaslB
TAx- Tb" Ilp Observe that ||xx b" b"1p Σx Eiec, (a,,)) - -x1l. Hence,
TAX Tb < <|AAx b" 1 S (1 + |T Tb"l
and thus
(1 - 8) TAX Tb"l - A b"1 || S (1+ +) TAX Tb"Ilp- -Ab--b'II.
Thus we can compute a vector RR such that
| blpp - (1 +
O (E))
min XER
| - Abll-
Before analyzing the time complexity of Algorithm 2, we first recall the following algorithms for Vandermonde matrix-vector multiplication and approximate lp regression.
Theorem A.4 (Vandermonde Matrix- Multiplication Runtime, e.g., Table 1 in (Gohberg & Olshevsky, 1994)). The runtime of computing Ax for a Vandermonde matrix A E Raxd and a (n log" n).
vector x E Rd for d S n is
Theorem A.5 (Approximate lp Regression Runtime). (Adil et al., 2019b) Given A E Rund and 2 12, 00), there exists an algorithm that makes O (Vnlog ") calls to a linear system solver and computes a vector such that
AA - bllp S (1 + e) min ||A - bllp-
xeRd
We now analyze the runtime of Algorithm 2.
Lemma A.6. Algorithm runs in O (n log n) + d.ww poly ( p. n) time.
Proof. Observe that Algorithm 2 has three main bottlenecks for runtime. Since we only need to Lewis weight sample from the extended matrices, we do not need explicitly form them, which would otherwise require (Rddp) time just list to entries. Hence the first bottleneck is performing the Lewis weight sampling procedure on the extended matrices. The second bottleneck is solving the lp regression problem the final subsampled matrix. The only remaining procedure is rounding and truncating the coordinates of b E R" to form a vector b" e R" using the procedure Round Trunc and then forming the groups GI.... G₁, which clearly takes O (n) arithmetic operations combined. We thus analyze each of the three main runtime bottlenecks.
First observe that the extended matrix M is Vandermonde matrix with O (dp²) columns. (Cohen & Peng, 2015) show that O (log n) matrix-vector multiplication operations can be done to compute approximate Lewis weights for the purposes of l, Lewis weight sampling. By Theorem A.4, each matrix-vector multiplication uses time O (n log² n). Hence computing the extended matrix M uses O (n log" n) time. Similarly, the extended matrix for each group GR is the product of Vander- monde matrix with O (dp²) columns and diagonal matrix. Thus by Theorem A.4, the extended matrices for all the groups G3 can be formed using O (n log n) time in total.
Since each Vandermonde matrix has O (dp?) rows, then observe that each group G2 samples O (d= logd) rows and there are = O (2 log n2) such groups k E [d]. Thus the resulting sub- sampled matrix has O dp² 2s
log221 rows for d S n. To approximately solve the l, regression
problem, Theorem A.5 notes that for 2 and subsampled matrix of size O ( log" d).
require only O ( pvd logn :)) calls to linear system solver. Moreover, on itera-
23/2
tion of the l, regression algorithm of (Adil et al., 2019a) used in Theorem A.5, the linear for a diagonal matrix Ry. Each lin-
system solves the equation X + (A TA, RR A ear system solve can be done in d² poly (1.p, log time. Hence, the total time to approxi- mately solve the Ep regression problem is d0.5+w poly (2 lognn. Therefore, the total runtime is O (n log' n) + 20.5+w poly ( logn).
15

---

Published as a conference paper at ICLR 2022
Lemma 2.4. Let x E I" and log 12
Then ||x||>0 & | S (1
+ [x||o
Proof. For any vector x E Rd, we have
1/8
| = max
Σ
las|
S n/
max [2a|- sen
Ken
Since (1 + 8)3/eps for all > O, then (1 + +)P > n for = 9 Therefore,
1/
max EEn
|ca| S (1 +
max | = (1 + x) En
Corollary 2.5. Given E (0, and p 2 I, E Ruxdy and RR, suppose A= [Aill... [A,] for Vandermonde matrices A..... Ag Raxd, Then there exists an algorithm that, with high probability, returns a vector & E Rdd such that
|AA - blbp S (1
min |Ax bblpp
XERd
using O (dp)) log n pply((dp)) log 1/e)) time, where T(A) is the runtime of mul- this can be further optimized to O (nd=a/2-1) 22-1 + poly((dp)² (ogg))),ooono n, 1/=)) time, where O (n²*) is the time to multiply an n X matrix
tiplying the matrix A by an arbitrary vector. For p
with an x n² matrix, so that w2 E [3.4].
Proof. Recall that key part in the proof of Theorem 1.1 was first the Egp regression on Van- dermonde matrix with dimension Rixd to ε regression Vandermonde matrix with dimension (2 where 2 S p <<p 2r+1, For a matrix A with block Vandermonde structure, we can
similarly write
2"
2"
9
((a,,x))
Cijzz
a (k- 1)+
(2 (d-1))++1
II
Pisk 2,2 (k- 1ya9i.
E
where again each yj is fixed function of the coordinates of Thus we can reshape the fz- regression problem on a matrix A with dimension nnx with block Vandermonde structure to an l regression problem a matrix À with dimension (2x(2"(d-1))*+ Moreover, we can further reshape Á into the concatenation of (dp)--1 Vandermonde matrices, where each Vandermonde matrix has columns that are geometrically growing a,2 but are multiplied by all (dp)" products
II2 kd where Pk E [dp].
matrix-vector multiplication each of the (dp) Vandermonde matrices. Thus by Theorem A.2, can & Lewis weight sample from the rows of the reshaped Å, using O (T(A)ddp)- log + (dp)"a) time. We can similarly write ((a,, x - )² for each tę among the discretized values of the updated b vector as a sum of (2" - 1))9
We
terms that are products of powers of the bases a,,2, 0i,24ds ... and variables V/. as in the proof of Theorem 1.1.
Thus we can partition the for regression problem into & = O logn instances of constrained
l regression problem (dp)e-1 Vandermonde matrices, each with at most 2n 1))² columns. To approximately solve the Ép regression problem, we can thus sample rows by their 6//2r Lewis weights, as Theorem 1.1. Since there are up (dp)--1 Vandermonde matrices, each with at most 2" (d 1))" + columns, then by Theorem A.2. the total time required is O (T(A)(d))*-1 log (dp)wa. poly(logn, 1/e)).
[
16

---

Published as a conference paper at ICLR 2022
We remark that for the special case of 2, (Sa et al., 2018) noted an efficient bivariate matrix multiplication algorithm of (Nüsken & Ziegler, 2004; Kedlaya & Umans, 2011).
Theorem A.7. (Nüsken & Ziegler, 2004; Kedlaya & Umans, 1; Sa el al., 2018) Given a q-variate polynomial f(X1..... .,,,,. such that each variable has degree at most d and N d9 distinct points æ<i)- (a(i))..... ((()g) for [N], there exists an algorithm that uses (d=z(e-1)/2+1)
time to ouiput the vector (J(x(1)), ..... J [xx(N) (N),, where O (1122) is the time multiply an matrix with an n x n² matrix, so that W2 E [3, 4/.
The case of matrix-vector product for g 2 corresponds to the evaluation of d2 points in The- orem A.7. Thus we need to repeat the algorithm in Theorem A.7 a total of d² times to handle all v rows in the input matrix. Since each instance of the algorithm uses O (d'z\/2+1) time, the total time for the matrix-vector product is O (nd"a/2-1)), rather than the naive O (nd) time (recall that W2 E ().
Corollary A.8. Given E (0, and p 2 1, A E Rax2d, and b Rd, suppose A = [A1]A2] for Vandermonde matrices A1, A2 E Rnxd returns vector & RRd such that
Then there exists an algorithm that with high probability
|A 18--bii bblp S (1 +
min |A
blip-
xER
using O (nd=2/2211 + poly((dp)², log n, 1/s)), where O (n²²) is the time to multiply an n x n matrix with an n x n² matrix.
A.2 l REGESSION ON OTHER STRUCTURED INPUT
Algorithm Faster regression for noisy low-rank matrices
Input: Rank l matrix K Raxd matrix S E Raxd with at most & non-zero entries per row, such that A K K S, measurement vector b, accuracy parameter E >> 0
Output: X E IR with ||A - bll, < (1 + O (e)) millgess ||x - bbl,,
1: < [log
2 S < 2
2: Extend A to a matrix M with dimension n x O (d' + + s)²"" so that each entry Mr,j is the
coefficient of the j-th term in the tensor decomposition of a, 622"
3: Use fp/z--Lewis weight sampling on M to find a set S of O (d' log d') indices in [n] and rescal-
ing factors, for d' O (d (k + s)²")
4: Let A' be the corresponding submatrix of A with indices in S and scaled accordingly. -blp-
5: Compute & min || ||A'x-
6: b' b AX, b" +- RoundTTrunc(()
O Log
7: Partition the rows of A into groups GI...., Ge, each containing all rows with the same value
of b"
8: Let Gj be the corresponding submatrix and l be the coordinate of b" corresponding G& for each E [f].
9: Use Leyiis weight sampling on [Gxital fx] to find a set Sy of O ( (% log d") indices in [n]
and rescaling factors, where d" = O (pd + s)²).
10: Let TĄ be the corresponding sampling and rescaling matrix for St
11: T .T....... :T.. Tx]
12: Compute < (1 + minxeR" |Axx TbblIp-
13:
return &
Lemma A.9. Given E 10, and p 2 1, a matrix A Raadd such that A K for rank k matrix K and s-sparse matrix S, and b E IRd, there exists algorithm that with high probability, returns a vector & E gd such that
A blp < (1+ min | blIp.
XER2
17

---

Published as a conference paper at ICLR 2022
Proof. The proof is similar to Lemma A.3. We once again let r be integer so that 2" S p < 2 and observe that
Σ
(a
x)I =
Σ
pp/2"
Σ
((k, + ts, x,, x))2"p12)
iEIN
iE(m)
iE(n]
Since K is low-rank matrix, then for all k;, we can write
k;
fixed set of basis vectors Vi,...:; Vi E Rd. Hence we have
1//22"
&
Σ Ka
x)I
S
+
ai, V
EE(n]
iE(n)
j=1
Since has sparsity 8, then we can further write
2/22"
& e Caisvsi
Σ (a. x) AE|m
sE(ni
j=1
where e1, eI,..., ed denote the elementary vectors. By the Hadamard Product-Kronecke Product mixed-produce property, we have
Σ KEINJ
Σ sE(m]
(y1 1118...0 (22))
where each Yk Va...Vi...... O,sA& S1,10 er.) for i E [2"]. Thus for fixed set of elementary vectors Gij:---- Gi there are (k s)2" possible values for the tensor product (y1 e ... ...).... Since there are (c) choices for the elementary vectors , ei,, then there are at most (Cdd (k + ))²") possible values for the tensor product for an absolute constant C' >> 0. Therefore, the lp subspace embedding problem A iad can be reshaped
constrained fp/2- subspace embedding problem on an input matrix of size n x (Cd' (I + s)²") Hence for É, regression with p E [2", 11, have
2d(d-1)+1
0/2"
xa,
Σ EE|m]
Σ
M
which is a constrained f//2" regression problem on a matrix M of size n x (Cd" + s)²") whose entries can be determined from the decomposition of each row of A.
Using 4//2" Lewis weight sampling, Theorem A.2 implies that we can find a matrix M' such that
12
My 9/2" p/2"
<
M'y
p/2r
<
13
12
|M |
p/2"
p/2"
for all y (CCd""(++))")) with high probability. By Lemma A.1, thus compute vector x such that
Ax blIp < 12 min Ax
Axpblpp-
xeRd
We again set b' b'= Aă and define OPT minxeRd ||x
blI
so that
||b 1|
||A blIp < 12OPT.
Let b' b Ax and OPT mingeR AA b bbI,, so that
|b'
||Ax bIIp S 12OPT.
18

---

Published as a conference paper at ICLR 2022
Let b" be the vector with all entries of b' rounded to the nearest power of (1 + e). starting the maximum entry of b' in absolute value, and stopping after we are poly(n) of that and replacing all remaining entries with O. By the triangle inequality, we have
||Ax blll < || Ax b'lb, ||b' b"ll 2 ||Ax b'll 12:OPT S (1+ + 12:)] Ax bbll for any R++1
Since the coordinates of b" can have e = O log at possible distinct values, we can partition the rows of A into groups, Gi...., Ge. based on the corresponding values of b". Let im be the corresponding value of b" for all rows in group Gms that
|Ax - b"l
Σ
(ai, x) Lmm..
iEG2
By the above argument, we have for each i Gi.
2"
((a; x) - tmm²² = -mm +
Bi Bi
Vi Σ 3=1
= Σ
Bistmalii.
where (1) B are entries of matrix B with (k+s) columns that can be computed from A, (2) 6m,1: 1, fm.2 ... are fixed values that can be computed from my and (3) each y, fixed function of the coordinates of x. Notably. B is the matrix formed by the concatenation of the coefficients of the decomposition of the a-fold tensor product of the row into the a-fold tensor products of the low- rank and sparse basis elements, for each defined in this proof is only the decomposition for a = p. Hence, the for regression problem on submatrix of A Rn the same coordinate of b", can be reshaped constrained & regression problem on a matrix with 2 C'd (k: + s)² columns.
0,... P. By comparison, the matrix M previously
The remainder of the proof follows from the same grouping argument Theorem 1.1. We apply Theorem A.2 by sampling rows of B corresponding to their 6//2" Lewis weights in the submatrix induced by the rows of Gm, we obtain matrix T, such that with high probability.
p/2
(k+s))
1/22 (T. Tm Vm
2/2
S
Σ iEG..
Σ
BistmjUii
S (T vm p/2
where Vm the vector restricted to the coordinates of b" that are equal to bm. Conditioning on the above inequality holding, it follows that
(1 - TT Tm vnl S Σ Ka x)
tml l2 S (1 TT Ax Tm vmll
iEG...
Therefore by summing over all m E [4], we have that with high probability,
(1 - E) Σ
Tm Vm 1 <
(a..x) - t s (1+5)
Σ Tm Tm
m
M
iEGm
For T = Tı o.... Té, we have Σn Ax
Tx b"l Since ||Ax- - -b1
Σm Eiec... ((,, t.m IP, then
(1-e)l Tb"l < ||A- S TAx-Tb"> Tbblpp
Therefore for p 221, 1,
(1-e)] - T S |Ax- - bb < (1+5)
[ TAx- - Tb
Because l = 0 [log and O = -0 (T(A)) we can compute a vector & E Rd such
that
- blIp < (1+ O1+ min | - Abbpp
19

---

Published as a conference paper at ICLR 2022
Lemma A.10. Given the low-rank factorization n poly (2", d", k,, s", loogg) n) time.
of K,
Algorithm
Proof. We analyze the runtime of Algorithm 4. First note that we can compute the extended matrix in time O (n((k + d) s*)) to perform the fp/2- Lewis weight sampling, where we recall that is the unique integer such that 2" p matrix, we require matrix-vector multiplication, which requires time (nk) for a low-rank matrix and time O (ns) for matrix whose rows have at most nonzero entries.
p < 2 To perform Lewis weight sampling on the extended
After the first iteration of Lewis weight sampling, we can round and truncate the coordinates of b R" form vector b" E R" using the procedure Round and then forming the groups GI.....Ga
formed, Ge. which clearly takes O (n) arithmetic operations combined. Once the groups compute the extended matrix in time O (n((k + 1)"dd )")) and perform fp/22 Lewis weight sampling each group, which takes O (nk + ns) total time across all groups. To approximately solve the resulting E, regression problem formed by the subsampled rows, we require poly (d", hP, sP, 2) time. Therefore, the total runtime n poly (2P, d", kP, sP, 2.g. log
Theorem 1.5 then follows from Lemma A.9 and Lemma A.I0. Theorem 1.6 is achieved through similar analysis for a noisy Vandermonde matrix. In particular, it follows from the same proof structure as Lemma A.3 by showing for 2" I p 2"+1, the fz- subspace embedding problem on A V + S for a given Vandermonde matrix V E Rnnx and sparse matrix S E Raxd can be reshaped as a constrained & subspace embedding problem on an input Vandermonde matrix of size d' where d' = (s²)
B APPLICATIONS TO POLYNOMIAL REGRESSION
In the polynomial regression problem, the goal is to find degree d polynomial é such that
(1((1) f(01BB < (1 + +))
min gideg(a)Ed
(1 [())))
where 8 an accuracy parameter given as input and ||-1 is the b, norm the pd power, lIa1l; 5_ ddt. The polynomial regression problem is fundamental problem in statisties, computational mathematies, machine learning, and more. The problem has been studied early as the 19th century with the work of Legendre and Gauss on least squares polynomial regression and has applications in learning half-spaces (Kalai al., 2008), solving parametric PDEs (Hampton & Doostan, 2015), and surface reconstruction (Pratt, 1987).
Given the flexibility choose query locations Ti,..., , problem as an active learning or experimental design problem. Thus we would like to minimize the number of queries 8, function of the approximation degree d, the norm p. and the accuracy parameter E, to find ğ. Observe that d+ queries are obviously necessary. but also that d+1 1 queries suffice when can be exactly by degree d polynomial, by using direct interpolation. In general, however, in the case when ming:deg(g)sd (t) [t] * O we require > d + 1 queries. Our Vandermonde Ep regression results can be used to give the first result showing that for all p 2 1, exp(O(p)) poly (t) queries suffice obtain a (1 + E) approxiaaaio to the best polynomial fit.
T.. we can consider the polynomial regression
We require the following structural theorem reducing the L, polynomial regression problem to problem of solving l, regression on Vandermonde matrices:
Theorem B.I. (Kane al., 2017; Meyer al., 202/) Suppose 81,--.,8og are drawn uniformly
from [- 1, 1, 1]. Let A E Rdtx(d+1) be the associated Vandermonde matrix, so that Aij s Let
70g = ( (zxtzr d³) and let E R. be the evaluations of f. so that b, = f(x)). Then with probability I1 the sketched solution & = argmin, || - bll, satisfies
1o:
|Px - fllE < (1++) min
|P - fll..
XER0+1
Theorem B.1 states that we can uniformly sample exp(O(p)) poly (d. 4) points from [-1, 1]. We can then form lp regression problem by using the evaluation each of the polynomial bases
20

---

Published as a conference paper at ICLR 2022
at the sampled points form the design matrix A and querying the underlying signal f at the sampled points to form the measurement vector b. Theorem B.1 says that the optimal solution & = the best fit degree d polynomial. We can naively
argming || bllp is a (1 +appoxximiiinionnn approximately solve the lp regression on the Vandermonde matrix A and the measurement vector b by standard Pp regression techniques such Lewis weight sampling, which would result in total query complexity exp(O(p))O (مotي d²). However, we can instead note that Lemma 1.4 implies we can instead solve an bg regression problem on a Vandermonde matrix with (dd)) columns for E [1. Crucially for E [1,2], there exist active by regression algorithms that only require reading ((d) poly (!) entries of b:
Theorem B.2. (Musco et al., 2021) Given p E [1, and input matrix A E Ruxd there exists an algorithm that reads Õ(d) poly (1) entries of b and with probability at least 0.99, outputs & such that
Ax bl1p - (1 + min ||A
bllp-
Hence by Theorem B.2, we can approximately solve the bg regression on a Vandermonde matrix with dp rows by reading (ddp) poly (t) entries of b. By Lemma 1.4, the approximate solution will also be a (1 (1-+=)-appp -)aaproxximation to the optimal lp regression on the Vandermonde matrix A. By Theorem B.1, the approximate solution will also form the coefficient vector of polynomial that is (1 + -aprroxaaatiiiiinnnnn to the polynomial bp regression problem. Therefore, we obtain the following guarantees for the polynomial regression problem:
Theorem B.3. For any degree d and norm 2 1, there exists an algorithm that queries f at = dp poly (log(dp). 2) points and outputs degree d polynomial @(t) such that
llā(1) - 1())B S (1 +e)
min [1((t) [()) qideg(a)<a
with probability at least z
The previous-best algorithm (Kane et al., 2017; Meyer et al., 2021) for ép polynomial regression sampled (dd) points uniform from [-1, 1] and then used standard lp regression algorithms that required reading the signal at all O(d)) sampled points, for a total query complexity of OdP). By comparison, Theorem B.3 only has linear dependency in d due the structural property of Lemma 1.4. Since ((d) queries are clearly necessary, our result settles the dependency of in the query complexity for l, polynomial regression for all d and all p 2 1.
C ADDITIONAL EXPERIMENTS
In this section, we experimentally verify that fp/2" Lewis Weight sampling works for unstructured matrices. We take similar approach as in Section to verify that lp/2" Lewis Weight sampling is correct for ép regression on unstructured matrices. For this test, we fix = 25, 000, d = 10, and p 6, while varying
[L. 1000]. We let A G, 0 G2 where G, E R 100x6 and G2 E 2210000x1 are i.i.d. N N(0, 1) matrices. To generate b, sample vector entries are (0, 100") and remaining 4 entries are (0, 1), and let b = Ax + z where Z E 25,000 a (0. 1) iid vector.
R200 whose first 6
We generate this matrix A and response vector b just once and run Ep/2r Lewis Weight sampling many times, so the variance in the plot comes only from the random sampling algorithms. Note that we again omit the rounding procedure on b. Figure 3 shows the result of this test, and we clearly see that the error shrinks quickly in m for our algorithm. This approach is much more practical than the prior Lewis Weight approximation method for unstructured matrices when p > 4. That approach required solving a non-linearly constrained SDP O (log n) times (Cohen & Peng, 2015), while our method requires only Gaussian sketch matrix and the standard Lewis Weight Iteration, which converges very quickly.
Since b is so large on its first 100 entries, it is important for any subsampling algorithm to sample at least 6 of the first 100 rows. Uniform sampling picks none of these rows until m ze 25000
100. = 250,
which is why uniform sampling fails to converge to good solution for small m. Lewis weight sampling instead gives much higher priority to the first 100 rows, avoiding any issue. This is why the the gap between Lewis Weight sampling and uniform sampling is so large for this experiment.
21

---

Published as a conference paper at ICLR 2022
10
Relative Error (Campiriaat)
10
100
10-1
1///2r Lewis Weights
Uniform
10"
10
10a
10
Subsample Complexity (m²)
Fig. 3: Empirical Relative Error (E Sempirical) versus subsample complexity m. We ran the experiment 50 times and plot the median, 25th quartile, and 75th quartile for each value of m.
D ADDITIONAL RELATED WORK
As previously mentioned, subspace embeddings are common tools used to approximately solve Ep regression. Given input matrix A E Buxd a subspace embedding is a matrix M Rmmxd with m < n such that
(1 e) Ax < ||Mx|lp S (1 + [[++)]|Ax]
for all x E Rdd Thus given an instance of é, regression, where the goal is to minimize ||Ax b set B = [A; b]. compute a subspace embedding for B, and then solve constrained 4p regression problem on the smaller matrix B.
A subspace embedding of matrix A Rxxd can be formed by sampling rows of A with proba- bilities proportional to their E, leverage scores, e.g., (Cormode et al., 2018; Dasgupta et al., 2008)
their ly. sensitivities, e.g., (Clarkson al., 2019; Braverman al., 2020; 2021; Musco al., 2021); I); or their ép Lewis weights, e.g., (Cohen & Peng, 2015; Durfee et al., 2018; Chhaya et al., 2020; Chen & Derezinski, 2021; Parulekar al., 2021; Meyer al., 2021). In any of these cases, the result- ing matrix M will contain subset of rows of A that are rescaled by function of their sampling probability, as give unbiased estimate of the actual E, mass.
In addition to sampling methods, sketching is a common approach for subspace embeddings. In these cases, the subspace embedding M is formed by setting M RAA for some (often random) matrix R E Rmx= The advantage of sketching over sampling is that sometimes R can be computed oblivious to the structure of A, whereas the sampling probabilities for each of the above distributions (ep leverage scores, lp sensitivities, and ep Lewis weights) are data dependent and thus require pass over the matrix A. The sketching matrix can be generated from a family of random matrices whose entries Cauchy random variables for Mahoney, 2013; Clarkson et al., 2016), or sub-Gaussian random variables for p = 2, e.g., (Sarlós, 2006; Nelson & Nguyen, 2013; Clarkson & Woodruff, 2013). More generally, exponential random variables can be used for 2 2, though the number of rows m in the resulting sketch matrix now has polynomial dependency in (Woodruff & Zhang, 2013). A line of recent works has studied the tradeoffs between oblivious linear sketches, sampling- based algorithms, and other sketches for Ep subspace embeddings (Wang & Woodruff, 2019; Li al., 2021). this paper, we focus on sampling-be algorithms for l, regression due to preservation of structure when the input matrix A itself has structure.
1, e.g., (Sohler & Woodruff, 2011; Meng &
22


Reviews:
[{"review_id": "OBLej5CqwHV", "reviewer": "Reviewer_JQBd", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths\n1. Author gives algorithms that return samples whose sizes are significantly smaller than existing sample sizes for both Lp and L_\\infty regression.\n2. They define a very useful property called rank for regression, which is extensively used in their analysis.\n3. The algorithms, theorems and lemmas are very well written.\n\nWekaness\n1. Missing related work. Authors should discuss how subspace embedding has been used to get a quick approximation of Lp regression (Clarkson., et al. \"The fast cauchy transform and faster robust linear regression.\", Chhaya, et al. \"Streaming coresets for symmetric tensor factorization.\", Woodruff and Zhang. \"Subspace embeddings and\\ell_p-regression using exponential random variables.\" etc). \n2. ||Ax-b||_p^p can be written as ||By||_p^p, for B = [A, b] (augmenting b as last column) and y = [x, -1] (setting last index as -1). Now by doing this the structure of B might slightly differ from the structure of A. So discussing the possibility of solving the actual problem by preserving the subspace of B would be a better start to motivate towards \"Rounding and truncating measurement vector\" (section 1.2). In particular, I am also curious that how much your guarantee improves by using the rounding and grouping technique instead of just relying on the structure of B (slightly differed from A). \n3. In the abstract it was mentioned that you get a sublinear time algorithm for L_\\infty regression but the running time of Theorem 1.2 is only O(n). Can you clarify, what do you mean by a sublinear time algorithm. \n4. Weak empirical evaluation. The author must give compare the performance (relative error and running time) from their sampling technique with other naive sampling techniques such as uniform, standard lewis weights sampling, or using Cauchy random variables (Clarkson et. al). \n\nMore comment\nIn theorem 1.1 and 1.2 the \\eps for S should be in (0,1/3) in order to get (1+\\delta) approximation to the optimal solution where \\delta \\in (0,1).\n\nUpdate:\nI appreciate for giving a satisfactory reply. I am willing to raise my score by 1.", "summary_of_the_paper": "The paper creates a coreset for solving Lp regression on structured input. For these cases, the paper shows that the coreset size is only polynomial in p. The author defines a notion of rank to Lp regression problem and uses a rounding and grouping technique on the response vector to improve their sample size. They also give a sublinear time algorithm for L_\\infty regression. ", "summary_of_the_review": "In this paper, the author shows good results for an Lp regression on structured input, but the motivation towards this technique is not well presented. It also lacks a discussion on related work and a decent experimental comparison with existing techniques. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "InYm0kPHk96", "reviewer": "Reviewer_2tK5", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Here are the strengths of the paper:\n\n(+) Theoretically, the results in the paper are solid and in particular, the result for $p=\\infty$ is new and is essentially the first non-trivial sublinear time result for $L_\\infty$ regression problem. I really like the $L_\\infty$ regression result.\n\n(+) The paper also presents a result for the general $A$ case (where the improvement is about a quadratic improvement).\n\n(+) The paper is pretty well-written (some minor typos at the end of this section)\n\nMost of my complaints are more from the more practical aspects of the results in the paper:\n\n(-) The improvements in this paper only take effect for $d>n^{2/p}$, which for fixed $p$ is still pretty large. Of course the improvement is much more dramatic for $p=\\infty$.\n\n(-) The paper presents some motivation for $L_p$ regression: it would be useful to clarify which of those need $p>4$. In other words, what are the practical applications for the case of $p>4$ and $d>n^{2/p}$ (for the specific structured $A$ considered in this paper)? I suspect this can handled by motivating the $L_\\infty$ regression problem on Vandermonde matrix.\n\n(-) It would be better to if in the discussion of techniques in Section 1.2, there is more clarity on how these techniques differ from existing techniques in the area. E.g. in the rounding and truncating step the only \"baseline\" provided is the simple rounding to nearest power of $(1+\\epsilon)$-- is there no \"better\" technical baseline to compare the techniques presented in the paper? If the rounding and truncating techniques really does not have any other baseline, then this part should probably be emphasized more.\n\n(-) Since the contribution of the paper is more on the theoretical side, it would have been good to see more proof details in the main part of the paper.\n\n(*) This is not really a concern but the paper did not have an ethics statement. While I can understand that the authors think that since the results in this paper are theoretical there are not ethical issues. While it is correct that there are no direct ethical issues, I would encourage the authors to think about potential mis-uses of the the algorithm by someone who potentially uses the results in this paper. Listing those in an ethical statement could be of use to others.\n\n\nDetailed Comments for the authors\n----------------------------------\n\n[Pg. 2, line below Thm 1.2] Is the claim on $L_\\infty$ regression for any $A$? Please clarify.\n\n[Pg. 3, Thm 1.5+1.6] What are the best known results for these matrices prior to this work?\n\n[Pg. 3, 3 lines below Thm 1.6] \"complexity complexity\"\n\n[Pg. 6, first display equation] $\\mathbf{s}$ should be $\\mathbf{x}$.\n\nUpdate after rebuttal\n------------------------\n\nI like the new result about query complexity of $L_p$ regression problem that the authors have added. Along with the result on $L_\\infty$ for (noisy) Vandemonde matrix, I think the paper is above the accept bar, so am changing my score to a 8.", "summary_of_the_paper": "This paper considers the problem of $L_p$-regression. I.e. given an $n\\times d$ matrix $A$ and a vector $b\\in\\mathbb{R}^n$, compute an approximation $\\tilde{x}\\in\\mathbb{R}^d$ such that ||A\\tilde{x}-b||_p <= (1+\\epsilon)\\min_{x\\in\\mathbb{R}^d} ||Ax-b||_p. In particular, the focus of this paper is on the case of large $p$ and when the matrix $A$ has structure (specifically when A is the Vandermonde matrix or matrices that are sum of sparse and Vandermonde/low rank matrices) and the goal is to perform the task in sub/near-linear time.\n\n\nThis paper focuses on the well studied method of sub-sampling the rows of $A$ (and re-scaling it) to reduce the problem size where more expensive iterative methods can be applied.\n\nFor simplicity we only talk about the results for the case of A being Vandermonde matrix. The results in this paper bring down the previous best know runtime bound of $O(n\\log^2{n})+d^{p/2+C}\\mathrm{poly}(1/\\epsilon)$ (for $p>4$ and some fixed constant $C$) to $O(n\\log^2{n})+p^3d^{3/2}\\mathrm{poly}(1/\\epsilon,\\log{n})$. Thus, for values $d >n^{2/p}$, this presents an improvement. The improvement is especially stark when we take $p=\\infty$, for which case this paper presents the _first_ sublinear-time $L_\\infty$ regression.\n\nTechnically the results in the paper follow from roughly three steps. The first step is to argue that for the structured matrices considered in this paper the rank of the regression problem is much smaller than the bound of $d^{\\Omega(p)}$ for arbitrary $A$. This gives the desired result if $b=\\mathbf{0}$ but to handle general $b$ needs more work. In particular, if all the entries of $b$ were roughly the same then the rank result is enough. Thus, it seems like the bulk of the technical work is in showing how to round and truncate $b$ and then in the third step solve the regression problem in \"groups\" where in each group, the sub-vector of $b$ is essentially \"fixed.\"\n\nThere are some experimental results presented that validate the results in the paper on Vandermonde matrix as well as results for general matrices (where the paper represents a quadratic improvement over the best-known results).", "summary_of_the_review": "Theoretically, the paper seems nice but it could be made stronger by providing motivations for the specific problem being studied in the paper. For the latter reason the paper is not a clear accept for me.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "-jwYgaciYzZ", "reviewer": "Reviewer_KU5L", "summary": "", "questions": "", "limitations": "", "rating": 10, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I find the new algorithmic ideas quite interesting. Much of the recent work on randomized numerical linear algebra, especially those involving p-norms, deal with worst case matrices. It's quite surprising to me that using additional structural information, it's possible to perturb entries of the matrices to simplify the sampling processes.\n\nA significant contribution (which I neglected in my initial read) is that the poly(p) * d bound shown here also work for the p > 2 case. To me this is highly surprising: all previous results require d^{p/2} rows, and such bound is tight for general matrices. This is very intriguing: it clearly demonstrates the value of the approach of looking at specialized matrix classes, and will likely motivate further investigations.\n\nA subtle issue about the initial experiments is that the gains of biased sampling can only be demonstrated on highly non-uniform matrices. This has been addressed in the revision via additional experiments. Here the incoherence of Vandermonde matrices also provides further motivations for studying such instances.", "summary_of_the_paper": "This paper studies ways of subsampling tall-and-dense p-norm regression involving structured Vandermonde matrices. It shows that in this setting with additional structure, sampling by Lewis weights produces poly(p) * d row sized samples for all values of p. Theoretically, this has two major advantages: it avoids more complicated computations of Lewis weights for higher values of p, and the sample size required is significantly better than the d^{p/2} bound for general matrices.\n\nThe algorithm is empirically evaluated on 25k-by-10 synthetic matrices, where it was observed that about 1k rows suffice to give sketches that incur at most 10% relative error. Furthermore, additional experiments provide cases where this type of sampling do significantly better than sampling uniformly.", "summary_of_the_review": "This paper takes an interesting, practically well motivated, study of sketching p-norm regression. The connections with polynomial regression is quite powerful, and I'm only slowly getting a sense of its value.\n\nBoth its theoretical and practical results are quite surprising and interesting to those working on optimization/randomized numerical algorithms.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Dear Authors,

The paper was received nicely and discussed during the rebuttal period. There is consensus among the reviewers that the paper should be accepted:

- The new result about query complexity of regression problem that the authors have added. Along with the result on 
 for (noisy) Vandemonde matrix, these make the paper lie above the accept bar.
- The authors have providing satisfying clarifications during the rebuttal that convinced reviewers to increase further their scores.

The current consensus is that the paper deserves publication.

Best AC
Decision: Accept (Poster)