Year: 2021

Abstract:
Lp Regression on Structured Inputs is an important problem in data analysis and machine learning where we find a vector x ∈ R that minimizes ‖Ax − b‖p for a structured matrix A ∈ Rn×d and response vector b ∈ R. Unfortunately, for many common classes of matrices, sampling-based algorithms for approximately solving Lp regression require runtime that is exponential in p, e.g., dO(p), which is prohibitively expensive. We show that for a large class of structured inputs, such as combinations of low-rank matrices, sparse matrices, and Vandermonde matrices, Lp regression can be approximately solved using runtime that is polynomial in p. For example, we show that Lp regression on Vandermonde matrices can be approximately solved using timeO(T (A) log n+(dp) ·polylogn), where T (A) is the time to multiply Ax for an arbitrary vector x ∈ R, and ω is the exponent of matrix multiplication. The polynomial dependence on p also crucially allows our algorithms to extend naturally to sublinear time algorithms for L∞ regression. Of independent interest, we develop a new algorithm for solving Lp regression for arbitrary matrices, which is significantly faster in practice for every p ≥ 4.

heading: 1 INTRODUCTION
text: Given a matrix A ∈ Rn×d and a vector b ∈ Rn, the goal of linear regression is to find a vector x ∈ Rd that best fits A to b. However for d < n, the linear system is overconstrained and thus there may not exist a vector x such that Ax = b. Thus the Lp regression problem seeks to find a vector x̃ ∈ Rd such that
‖Ax̃− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p, for some approximation parameter ε > 0. Here for a vector y ∈ Rn, ‖y‖pp = ∑n i=1 |yi|p. The Lp regression problem is a central question in statistical data analysis that has numerous applications in machine learning, data science, and applied mathematics (Friedman et al., 2001; Chatterjee & Hadi, 2006) and thus has received a considerable amount of attention. A direct way to solve the Lp regression problem is to use an iterative solver (Adil et al., 2019b;a) that makes poly(n) calls to a linear system solver. However, the linear system solver uses time poly(n, d), which can be prohibitively slow for n d. Therefore, a standard approach for Lp regression is to first sample a number of rows of A down to a smaller matrix M that generally has poly ( d, 1ε , log n ) rows by generating a sampling and rescaling matrix S such that M = SA and ‖SAx − Sb‖p is a “good” approximation to ‖Ax − b‖p for all x ∈ Rd. Subsequently, an iterative solver can now be used on M, which serves as a proxy for the original problem but has smaller input size, thus allowing for more efficient computation. Unfortunately, to handle worst-case inputs, M requires sampling O ( dmax(1,p/2) ) rows, so that not
only is the overall runtime exponential in p, but also the sampling is only helpful if n dp/2. Fortunately, there are many applications in which the design matrix A is known to have additional structure that can be leveraged to design more efficient algorithms. For example, Vandermonde matrices are used in the polynomial regression problem, which has been studied for over 200 years (Gergonne, 1974) and recently motivated by applications to machine learning (Kalai et al., 2008), applied statistics (MacNeill, 1978), and computer graphics (Pratt, 1987) . The goal is to fit a continuous signal, which can be measured at time choices t1, . . . , tn to form a measurement vector b ∈ Rn, using the best degree d polynomial whose coefficients can be represented by the vector x ∈ Rd+1. Thus, entry Ai,j of the feature matrix A is the value of the j-th polynomial basis at time ti so that the
problem has been restated as the Lp regression problem with a Vandermonde matrix A, where the entries of each row ai form a sequence 1, ti, (ti)2, . . . , (ti)d−1 for some parameter ti. Regression problems with Vandermonde matrices also arise in the settings of Fourier fitting (Avron et al., 2019) and Toeplitz covariance estimation (Eldar et al., 2020). Shi & Woodruff (2019) show that for Vandermonde matrices A, the time used to compute the smaller subsampled matrix M can be improved upon the naı̈ve approach by using fast Vandermonde matrix-vector multiplication. However, the number of rows of M is still exponential in p so that unfortunately, the runtime of the iterative solver for Lp regression still uses time exponential in p, which suffices for small p but crucially does not obtain any improvements for the important problem of L∞ regression.

heading: 1.1 OUR CONTRIBUTIONS
text: We first show that for Lp regression on Vandermonde matrices and its noisy perturbations, the size of the subsampled matrix M and therefore the runtime of the iterative solver only needs polynomial dependence in p. Theorem 1.1. Given ε ∈ (0, 1) and p ≥ 1, a Vandermonde matrix A ∈ Rn×d, and b ∈ Rd, there exists an algorithm that creates a sampling matrix S ∈ Rm×n with m = O ( p2d ε3 log 2 n )
rows so that with high probability,
(1− ε)‖Ax− b‖p ≤ ‖SAx− Sb‖p ≤ (1 + ε)‖Ax− b‖p, for all x ∈ Rd. The algorithm uses O ( n log2 n ) + p3d3/2 poly ( 1 ε , log n ) total time and with high probability, returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p.
The best known previous work has runtime either (1) O ( n log2 n+ dq · poly(1/ε) ) time (Shi & Woodruff, 2019), where q = ω ∈ [2, 3) is the exponent for matrix multiplication for p ∈ [1, 4] and q = p/2 + C for some fixed constant C > 0 for p > 4, which has an exponential dependence in p, or (2) using poly(n) solves on the n × d linear system Adil et al. (2019b;a). By contrast, our work only has a polynomial dependence on p.
In fact, we show that to obtain a (1 + ε)-approximation to the fundamental problem of L∞ regression, it suffices to consider Lp regression for p = O ( logn ε ) and since our results have polynomial dependence on p rather than exponential, we obtain sublinear size coresets for L∞ regression: Theorem 1.2. Given ε ∈ (0, 1), a Vandermonde matrix A ∈ Rn×d, and b ∈ Rd, there exists an algorithm that uses O ( n log2 n ) + d3/2 poly ( 1 ε , log n ) time that computes a sampling matrix
S ∈ Rm×n with m = d poly ( 1 ε , log n ) such that with high probability,
(1− ε)‖Ax− b‖∞ ≤ ‖SAx− Sb‖∞ ≤ (1 + ε)‖Ax− b‖∞ for all x ∈ Rd. In particular,
min x∈Rd ‖SAx− Sb‖∞ ≤ (1 + ε) min x∈Rd ‖Ax− b‖∞.
To the best of our knowledge, this is the first known dimensionality reduction for L∞ regression with provable guarantees for any (nontrivial) input matrix.
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 ≥ 1 and a matrix A ∈ Rn×d, suppose there exists a matrix M ∈ Rn×t and a fixed function f : Rd → Rt so that for all x ∈ Rd,
|〈ai,x〉|p = |〈mi, f(x)〉|.
Then we call t the rank of the Lp regression problem.
Our result relies 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 ≥ 1, the rank of the Lp regression problem on a Vandermonde matrix A ∈ Rn×d is O (dp).
Lemma 1.4 implies that the Lp loss function for a row of a Vandermonde matrix can potentially be expressed as a linear combination of O (dp) variables even though the entries of the measurement vector b can be arbitrary. By comparison, the Lp loss function for a row ai of a general matrix A is |〈ai,x〉 − bi|pp can only be expressed as a linear combination ofO (pdp) variables, corresponding to each of the dk k-wise product of coordinates of x, for each k ∈ [p]. As a corollary of Lemma 1.4, we obtain a small coreset for Lp regression (as well as L∞ regression) on a Vandermonde matrix, which can thus be used as a preconditioner. We summarize these results in Table 1.
We generalize Lemma 1.4 to similar guarantees for Lp regression on a matrix A that is the sum of a low-rank matrix L and a sparse matrix S that has at most s non-zero entries per row.
Theorem 1.5. Given ε ∈ (0, 1) and p ≥ 1, a matrix A ∈ Rn×d such that A = K + S for a rank k matrix K and a s-sparse matrix S, and b ∈ Rd, there exists an algorithm that uses O (nk + ns) + poly ( ds, kp, sp, 1ε ) time, and with high probability, returns a vector x̂ ∈ Rd with
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p.
Similarly, we obtain efficient guarantees for Lp regression on a matrix A that is the sum of a Vandermonde matrix V and a sparse matrix S that has at most s non-zero entries per row.
Theorem 1.6. Given ε ∈ (0, 1) and p ≥ 1, a matrix A ∈ Rn×d such that A = V + S for a Vandermonde matrix V and an s-sparse matrix S, and b ∈ Rd, there exists an algorithm that usesO ( T (A)) log n+ poly ( p2d, ds, sp, log n, 1ε )) time and with high probability, returns a vector x̂ ∈ Rd such that ‖Ax̂− b‖p ≤ (1 + ε) min
x∈Rd ‖Ax− b‖p.
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 et al. (2021) for a lower bound) of roughly O ( dp/2 ) . Furthermore, the advantage of our approach is that we now only need to perform Lq Lewis weight sampling for q ∈ [1, 4]. There are known efficient iterative methods for computing the Lq Lewis weights for q ∈ [1, 4], e.g., see Section 2.1 or Section 3 in Cohen & Peng (2015). By contrast, previous methods relied on computing general Lp weights for p > 4, which required solving a large convex program, e.g., through semidefinite programming, see Section 4 in Cohen & Peng (2015).
Finally, we experimentally validate our theory on synthetic data. In particular, we consider matrices and response vectors that are motivated by existing lower bounds for subsampling and sketching methods. For Vandermonde regression, our experiments demonstrate that the number of rows needed is subexponential in p, reinforcing how structured matrices outperform the worst-case bound. For unstructured matrix Lp regression, we demonstrate that our Lq Lewis Weight subsampling scheme is effective and accurate. For instance, we can find a near-optimal solution vector by solving the regression problem on just 2% of the input matrix’s rows.

heading: 1.2 OVERVIEW OF OUR TECHNIQUES
text: Our algorithmic contributions rely on two key observations, which we describe below. For the sake of presentation, we assume p is an integer in this overview, though we handle arbitrary real values of p in the subsequent algorithms and analyses.
Reduced rank of the Lp 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 ∈ Rn×d does not need to be the O (dp) rank of the Lp regression problem on general matrices. For a general matrix, 〈ai,x〉p for an integer p can be rewritten as a linear combination ∑ αiyi of O (dp) terms, where each coefficient αi is a p-wise product of entries of the row vector ai ∈ Rd and similarly each yi is a p-wise product of entries of the vector x ∈ Rd. However, when A is a Vandermonde matrix, then the j-th entry of ai is simply A j−1 i,2 . Thus each coefficient in 〈ai,x〉p can be written as a linear combination of 1, Ai,2, (Ai,2)2, . . . , (Ai,2)p(d−1) and so we can express 〈ai,x〉p as a linear combination of O (dp) variables rather than O (dp) variables. We can similarly show that the rank of the Lp regression problem on rank k-matrix A is O (kp) by writing each row ai as a linear combination of basis vectors v1, . . . ,vk. Then, using the Hadamard Product-Kronecker Product mixed-product property, we can rewrite 〈ai,x〉p as a linear combination of p-wise products of the variables 〈v1,x〉, . . ., 〈vk,x〉, i.e., a linear combination ofO (kp) variables rather than O (dp) 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 (dssp) by noting that (1) there are ( d s ) = O (ds) sparsity patterns and that (2) for a fixed sparsity pattern, the rank of the induced matrix is at most s, which induces a linear combination of O (sp) variables for the Lp regression problem. These decomposition techniques can also be generalized to show that the rank of theLp regression problem is low on matrices A such that A = K + V, A = K + S or A = V + S, where K is a 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 Lp 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 |〈ai,x〉 − bi|p rather than 〈ai,x〉p and there can be n possible different values of bi across all i ∈ [n]. Rounding and truncating the measurement vector: a novel algorithmic technique. Thus, the second main ingredient is manipulating the measurement vector b ∈ Rn 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 + ε). Unfortunately, such a rounding approach would roughly preserve ‖b‖p up to a multiplicative (1 + ε) factor, but it would not preserve ‖Ax−b‖p. For example, suppose 〈ai,x〉 = bi = N for some arbitrarily large valueN . Then |〈ai,x〉−bi|p = 0, but if bi were rounded to (1+ε)N , we would have |〈ai,x〉−bi|p = εpNp, which can be arbitrarily large.
The lesson from this counterexample is that when bi is significantly larger than 〈ai,x〉 − bi, any rounding technique can be arbitrarily bad because ‖b‖p can be significantly larger than OPT := minx∈Rd ‖Ax − b‖p. On the other hand, if ‖b‖p is a constant factor approximation to OPT, then the previous counterexample cannot happen because either (1) bi is large relative to ‖b‖p and 〈ai,x〉 cannot be too close to bi so that rounding bi will not significantly affect the difference 〈ai,x〉 − bi or (2) bi is small relative to ‖b‖p and so any rounding of bi 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 ‖b‖p = O (OPT).
To that end, we note that if ‖Ax̃−b‖p = O (OPT) for a vector x̃ ∈ Rd, then by a triangle inequality argument, we have that the residual vector b′ = b−Ax̃ ∈ Rd satisfies ‖b′‖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., a matrix M ∈ Rm×d with m n such that
(1− C)‖Ax‖p ≤ ‖Mx‖p ≤ (1 + C)‖Ax‖p
for some constant C ∈ (0, 1). Typically, such a matrix M can be found by sampling the rows of A according to their Lp Lewis weights to generate a matrix withm = O ( dmax(1,p/2) ) rows. However, using our observation that the rank of the Lp regression problem on structured inputs can be low, we can sample O (dp) rows of A with probabilities proportional to their Lq Lewis weights, where q = p2r for an integer r such that 2
r ≤ p < 2r+1. Crucially, we can find such a vector x̃ without reading the coordinates of b since we only require to read the rows of A to perform Lewis weight sampling in this phase. Thus we can efficiently find a residual vector b′ such that ‖b′‖p = O (OPT).
Partitioning the matrix into groups. We then round the entries of b′ to obtain a vector b′′ with at most ` := O ( logn ε ) unique values, truncating all entries that are less in magnitude than 1poly(n) to zero instead. We now solve the Lp regression problem minx∈Rd ‖Ax− b′′‖p by partitioning the rows of A into ` groups G1, . . . , G` based on their corresponding values of b′′. The main point is that all rows in each group Gk have the same entry tk in b′′. Thus we can again observe that |〈ai,x〉 − tk|p can be written as a linear combination of O ( p2d )
variables and therefore use Lq Lewis weight sampling (for the same q) to reduce each group Gk down to O ( p2d ε2 log d ) rows.
Since there are ` = O (
logn ε
) groups, then there are roughly O ( p2d )
total rows that have been sampled across all the groups. It follows from the decomposition of the Lp 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 a vector x̂ with the desired property that
‖Ax̂− b‖p ≤ (1 +O (ε)) min x∈Rd ‖Ax− b‖p.
Practical Lp regression for arbitrary matrices. We now describe a practical procedure for Lp regression on general matrices that avoids the necessity for convex programming to approximate the Lp Lewis weights for p > 4. We first pick an integer r such that p2r ∈ [2, 4). By the same structural argument as in Lemma 1.4, we can tensor product each row ai of A with itself 2r times, thus obtaining an extended matrix of size n × d′, where d′ = d2r , independent of the entries in b. We then Lp/2r Lewis weight sample on the extended matrix, using the iterative method in Figure 1 since p2r < 4 rather than solving a convex program for Lp Lewis weight sampling for p > 4. Since p 2r ∈ [2, 4), it follows from Theorem 2.3 that Lewis weight sampling requires roughly (d ′)p/2 rows.
Thus we obtain a matrix with O ( (d2 r )p/2 r+1 ) = O ( dp/2 ) rows, from which we can compute a
residual vector b′ such that ‖b′‖p = O (OPT), where OPT := minx∈Rd ‖Ax− b‖p. We then round and truncate the entries of b′ using the subroutine RoundTrunc to obtain a vector b′′, which allows us to partition the rows of A and the entries of b′′ into O ( logn ε ) groups. Due to the constant values of b′′ in each group, we can again Lp/2r 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

heading: 2.1 PRELIMINARIES ON LEWIS WEIGHT SAMPLING
text: There are a number of known sampling distributions for dimensionality reduction for the Lp regression problem, such as the Lp leverage scores, e.g., (Cormode et al., 2018; Dasgupta et al., 2008) and the Lp sensitivities, e.g., (Clarkson et al., 2019; Braverman et al., 2020; Musco et al., 2021a); in this paper we focus on the Lp Lewis weights, e.g., (Cohen & Peng, 2015; Durfee et al., 2018; Chen & Derezinski, 2021; Parulekar et al., 2021).
Definition 2.1 (Lp Lewis Weights, Cohen & Peng (2015)). Given a matrix A ∈ Rn×d and p ≥ 1, the Lp Lewis weights w1(A), . . . , wn(A) are the unique quantities that satisfy
(wi(A)) 2/p = a>i (A >W1−2/pA)−1ai
for all i ∈ [n], where W ∈ Rn×n is the diagonal matrix with Wi,i = wi(A) for all i ∈ [n]. Definition 2.2 (Lewis Weight Sampling). For an input matrix A ∈ Rn×d and m samples, let the sampling matrix S ∈ Rm×n be generated by independently setting each row of S to be the i-th
standard basis vector multiplied by (
d m·wpi (A)
)1/p with probability w p i (A)
d .
Theorem 2.3 (Lp Subspace Embedding from Lewis Weight Sampling, (Cohen & Peng, 2015)). For any A ∈ Rn×d, let the matrix S ∈ Rm×n be generated from Lewis weight sampling with
m = O ( dmax(1,p/2) log(d/δ) log(1/ε)
εγ
) , where γ = 2 for p ∈ [1, 2] and γ = 5 for p > 2. Then with
probability at least 1− δ, we have that simultaneously for all x ∈ Rd,
(1− ε)‖Ax‖p ≤ ‖SAx‖p ≤ (1 + ε)‖Ax‖p.
Since the Lewis weights are implicitly defined, it may not be clear how to compute them exactly. In fact, it suffices to compute constant factor approximations to the Lewis weights. Cohen & Peng (2015) show that for p ≤ 4, there exists a simple iterative approach to compute a constant factor approximation to the Lp Lewis weights in input sparsity time, which we present in Figure 1.
At a high level, the correctness of Figure 1 follows from Banach’s fixed point theorem and the fact that the subroutine LewisIterate is a contraction mapping, because |2/p− 1| < 1 for p < 4 (Cohen & Peng, 2015). However, for p ≥ 4, Figure 1 no longer works. Instead, approximating Lp Lewis weights for p > 4 seems to require solving the convex program
Q = argmax M det(M), subject to ∑ i (a>i Ma) p/2 ≤ d, M 0,
and setting wi = (a>i Qa) p/2. Unfortunately, this is often infeasible in practice, and we could not obtain empirical results using it. Therefore, a nice advantage of our algorithms, both for structured and unstructured matrices, is that we only use Lq Lewis weight sampling for q ≤ [1, 4), even if p ≥ 4, whereas previous algorithms required using Lp Lewis weight sampling for p > 4.

heading: 2.2 ALGORITHM AND ANALYSIS
text: We first describe the general framework of our algorithm for efficient Lp regression, so that given A ∈ Rn×d and b ∈ Rd, the goal is to approximately compute OPT := minx∈Rd ‖Ax − b‖p. Recall that in order to apply our structural results, we first require a measurement vector b′′ with a small number of distinct entries. We obtain b′′ by first finding a constant factor approximation, i.e., x̃ ∈ Rd such that ‖Ax̃ − b‖p ≤ O (OPT). We can Lp Lewis weight sample the rows of A 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 a matrix M from A with fewer than dp columns and then Lp/2r Lewis weight sample the rows of A, where r is the integer that satisfies 2r ≤ p < 2r+1, and then we solve the Lp regression problem on the sampled rows of A and entries of b to obtain x̃. We set b′ to be the residual vector b−Ax̃ so that ‖b′‖p = O (OPT). We then use the procedure RoundTrunc, i.e., as in Algorithm 1, which sets the entries of b′ that are the smallest in magnitude to zero, and rounds the remaining entries of b′ to the nearest power of (1 + ε).
Since ‖b′‖p = O (OPT), it then follows by triangle inequality that it suffices to approximately solve the Lp regression problem on the vector b′′ instead. We partition the rows of A into groups
Algorithm 1 RoundTrunc: Round and truncate coordinates of input vector b Input: Vector b ∈ Rn, accuracy parameter ε > 0 Output: Vector x ∈ Rn that is a rounded and truncated version of b
1: M ← maxi∈n |bi| 2: for i = 1 to i = n do 3: pi ← blog1+ε |bi|c 4: xi ← (1 + ε)pi · sign(bi) 5: if |xi| ≤ Mn5 then 6: xi ← 0 7: return x
G1, . . . , G` based on the values of b′′ and performLp/2r Lewis weight sampling again on an implicit matrix that we create from the rows of each group. Here we leverage our theory about the rank of the Lp regression problem for structured inputs. It then suffices to solve the Lp regression problem on the matrix formed by the sampled rows across all the groups.
Algorithm 2 Faster regression for Vandermonde matrices
Input: Vandermonde matrix A ∈ Rn×d, measurement vector b, accuracy parameter ε > 0 Output: x̂ ∈ R with ‖Ax̂− b‖p ≤ (1 +O (ε)) minx∈Rd ‖Ax− b‖p
1: r ← blog pc . 2r ≤ p < 2r+1. 2: Extend A to a Vandermonde matrix M with dimensions n× (2r(d− 1) + 1). 3: Use Lp/2r -Lewis weight sampling on M to find a set S of O (d log d) indices in [n] = {1, 2, . . . , n}, and corresponding rescaling factors. 4: Let A′ be the corresponding submatrix of A with indices in S, and scaled accordingly. 5: Compute x̃ ≤ 5 minx∈Rd ‖A′x− b‖p. 6: b′ ← b−Ax̃, b′′ ← RoundTrunc(b′), `← O ( logn ε
) 7: Partition the rows of A into groups G1, . . . , G`, each containing all rows with the same value
of b′′ 8: Let Gk be the submatrix of Gk extended to 22r(d− 1) + 1 columns. 9: Let tk be the coordinate of b′′ corresponding to Gk for each k ∈ [`].
10: Use Lp/2r -Lewis weight sampling on [Gk; tk] to find a set S′k ofO ( 1 εd log d ) indices in [n] and
rescaling factors. 11: Let Tk be the corresponding rows with indices in S′k and scaled accordingly. 12: T← [T1; . . . ;Tk]> 13: Compute x̂ ≤ (1 + ε) minx∈Rd ‖Tx− b‖p. 14: return x̂
We show in the supplementary material that Algorithm 2 satisfies the guarantees of Theorem 1.1. Moreover, because Theorem 1.1 has polynomial dependence on p rather than exponential dependence, we obtain the first known sublinear size coreset for the important problem of L∞ regression. We use the following structural property.
Lemma 2.4. Let x ∈ Rn and p = Ω (
logn ε ) . Then ‖x‖∞ ≤ ‖x‖p ≤ (1 + ε)‖x‖∞.
Lemma 2.4 implies that to solve L∞ regression, we can instead solve Lp regression for p = Ω (
logn ε
) . Then Theorem 1.2 follows from the fact that even for p = Ω ( logn ε ) , the matrix T
in Algorithm 2 satisfies
(1− ε)‖Tx− b′′‖p ≤ ‖Ax− b′′‖p ≤ (1 + ε)‖Tx− b′′‖p for all x ∈ Rd. Hence, we can solve the Lp regression problem on the smaller matrix T to solve the L∞ regression on A.
The results of Theorem 1.1 can be further extended to matrices with block Vandermonde structure. Corollary 2.5. Given ε ∈ (0, 1) and p ≥ 1, A ∈ Rn×dq , and b ∈ Rd, suppose A = [A1| . . . |Aq] for Vandermonde matrices A1, . . . ,Aq ∈ Rn×d. Then there exists an algorithm that, with high
probability, returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p,
using O ( T (A)(dp)q−1 log n+ poly((dp)q log n, 1/ε) ) time, where T (A) is the runtime of mul-
tiplying the matrix A by an arbitrary vector. For q = 2, this can be further optimized to O ( ndω2/2−1 + poly((dp)2, log n, 1/ε) ) time, whereO (nω2) is the time to multiply an n×nmatrix
with an n× n2 matrix, so that ω2 ∈ [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 supplementary material.
Algorithm 3 Faster Lp regression for general matrices
Input: Matrix A ∈ Rn×d, measurement vector b, accuracy parameter ε > 0 Output: x̂ ∈ R with ‖Ax̂− b‖p ≤ (1 + ε) minx∈Rd ‖Ax− b‖p, parameter p ≥ 4
1: r ← blog pc − 1 . 2r+1 ≤ p < 2r+2. 2: Extend A to a matrix M with dimension n × O ( d2 r)
so that each row in Mi is the 2r-fold tensor product of ai reshaped into a d2 r
length vector. 3: Use Lp/2r -Lewis weight sampling on M to find a set S of O ( dp/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 minx∈Rd ‖A′x− b‖p. 6: b′ ← b−Ax̃, b′′ ← RoundTrunc(b′), `← O ( logn ε
) 7: Partition the rows of A into groups G1, . . . , G`, each containing all rows with the same value
of b′′ 8: Let Gk be the corresponding submatrix and tk be the coordinate of b′′ corresponding to Gk for
each k ∈ [`]. 9: Use Lp/2r -Lewis weight sampling on [Gk; tk] to find a set S′k of O ( 1 ε5 d p/2 log d )
indices in [n] and rescaling factors.
10: Let Tk be the corresponding rows with indices in S′k and scaled accordingly. 11: T← [T1; . . . ;Tk]> 12: Compute x̂ ≤ (1 + ε) minx∈Rd ‖Tx− b‖p. 13: return x̂
Finally, we describe in Algorithm 3 a practical approach for Lp regression for arbitrary matrices without requiring any structural assumptions.
Theorem 2.6. Algorithm 3 outputs a vector x̂ such that ‖Ax̂−b‖p ≤ (1+ε) min ‖Ax−b‖p. The runtime of Algorithm 3 is T (n, dp/2) + poly ( dp/2, log n, 1ε ) , where T (n, dp) is the time to multiply a matrix of size n× dp by a vector of length dp.
Theorem 2.6 achieves the same optimal sample complexity as previous Lp Lewis weight algorithms of roughly dp/2. However, the main advantage of Algorithm 3 is that it performs Lq Lewis weight sampling for some q ∈ [2, 4), which is quite efficient because we can use an iterative method rather than solving a convex program.

heading: 3 EMPIRICAL VERIFICATION
text: We provide empirical evidence to validate our core statistical claim about Vandermonde regression: that the relative error ε achieved by Lp/2r 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:
εempirical = ‖Vx̂− b‖p − ‖Vx∗ − b‖p
‖Vx∗ − b‖p
Our theory tells us that (εempirical)3 ≤ Õ ( dp2
m
) , where m is the total number of subsampled
rows. The prior work on unstructured matrices instead suggests (εempirical)c ≤ Õ ( dO(p)
m
) (Shi &
Woodruff, 2019). So, to visually distinguish these two settings, we look at the logarithm of both sides:
log(εempirical) ≤ O ( ln(p) + ln( dm ) ) OR log(εempirical) ≤ O ( p ln(d) + ln( 1m ) ) In particular, our work suggests a logarithmic dependence on p, while the prior work suggests a
linear dependence on p.
To validate our theory, we plot log(εempirical) versus p and m on synthetic data. Specifically, we generate n = 25, 000 i.i.d. N(0, 1) times samples to form a Vandermonde matrix V with d = 20 columns, then compute the polynomial q(t) = t10 at each time sample, add N(0, 1010) additive noise, and save the corresponding values in b. We then compute εempirical for this Lp 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 x̂ by solving minx ‖V̂x − b̂‖p where V̂ and b̂ are computed by sampling and rescaling V and b with the Lp/2r Lewis Weights.
Figure 2 shows the result of these experiments, which were run in Julia 1.6.1, on Windows 10 with an Intel i7-7700K CPU and 16Gb RAM. In Figure 2a, we fix m = 1000 and vary p ∈ [2, 25]. The trendline of Lewis Weight sampling clearly better fits a logarithmic model, as opposed to a linear model. This reinforces our analysis by showing that the dependence on p is notably sub-exponential, beating the known bounds for Lp subsampling on unstructured matrices.
As a benchmark, we compare our Lewis Weight sampling method to uniform sampling. The noise in b is large enough that most rows of 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 m ∈ [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.

heading: A MISSING PROOFS FROM SECTION 2
text: In this section, we give the missing proofs from Section 2.
A.1 Lp REGESSION ON VANDERMONDE MATRICES
We first prove a simple statement that shows a “good” solution to the Lp regression problem on a subspace embedding SA is also a “good” solution to the original input matrix A. Lemma A.1. Let S be a sampling and rescaling matrix such that
11 12 ‖SAx‖p ≤ ‖Ax‖p ≤ 13 12 ‖SAx‖p, E
[ ‖SAx‖pp ] = ‖Ax‖pp
for all x ∈ Rd. Let OPT = minx∈Rd ‖Ax− b‖p and let x ∈ Rd be any vector for which ‖SAx− Sb‖p ≤ 5OPT. Then with probability at least 0.79,
‖Ax− b‖p ≤ 12OPT.
Proof. Let x∗ be a minimizer of ‖Ax− b‖p so that OPT = ‖Ax∗ − b‖p. We will prove the claim by contrapositive, so we first suppose that ‖Ax− b‖p ≤ 12OPT. Then by the triangle inequality,
‖SAx− Sb‖p ≥ ‖SA(x− x∗)‖p − ‖SA∗ − Sb‖p.
Since 1112‖SAx‖p ≤ ‖Ax‖p ≤ 13 12‖SAx‖p for all x ∈ R d, then
‖SAx− Sb‖p ≥ 11 12 ‖A(x− x∗)‖p − ‖SA∗ − Sb‖p.
By the triangle inequality,
‖SAx− Sb‖p ≥ 11 12 (‖Ax− b‖p − ‖Ax∗ − b‖p)− ‖SA∗ − Sb‖p.
Note that by Jensen’s inequality, E[‖SAx∗ − Sb‖p] ≤ OPT, so that by Markov’s inequality,
Pr [‖SAx∗ − Sb‖p ≥ 5OPT] ≤ 1
5 .
Thus with probability at least 0.8,
‖SAx− Sb‖p ≥ 11 12 (‖Ax− b‖p − ‖Ax∗ − b‖p)− 5OPT.
Thus if ‖Ax− b‖p ≤ 12OPT, then
‖SAx− Sb‖p ≥ 11
12 (12OPT− OPT)− 5OPT > 5OPT,
as desired.
Before justifying the correctness of Algorithm 2, we first recall the following algorithm for efficient Lp subspace embeddings. Theorem A.2. (Shi & Woodruff, 2019) Given ε ∈ (0, 1) and p ∈ [1, 2], a Vandermonde matrix A ∈ Rn×d, and b ∈ Rd, let T (A) be the time it takes to perform matrix-vector multiplication, i.e., compute Av for an arbitrary v ∈ Rd. There exists an algorithm that uses O (T (A) log n+ dq · poly(1/ε)) time, where q = ω for p ∈ [1, 4) and q = p/2 + C for some fixed constant C > 0 for p > 4, and with high probability, returns x̃ ∈ Rd such that
‖Ax̃− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p.
We now justify the correctness of Algorithm 2. Lemma A.3. Given ε ∈ (0, 1) and p ≥ 1, a Vandermonde matrix A ∈ Rn×d, and b ∈ Rd, then with high probability, Algorithm 2 returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 +O (ε)) min x∈Rd ‖Ax− b‖p.
Proof. Consider Algorithm 2 and let r be an integer so that 2r ≤ p < 2r+1. Note that∑ i∈[n] |〈ai,x〉|p = ∑ i∈[n] |(〈ai,x〉)2 r |p/2 r .
Since A is Vandermonde, then
(〈ai,x〉)2 r =  d∑ j=1 ai,jxj 2 r =  d∑ j=1 aji,1xj 2 r = 2r(d−1)+1∑ j=1 aj−1i,2 yj ,
where each yj is a fixed function of the coordinates of x. Notably, the fixed function is the same across all i ∈ [n]. Hence, the L2r subspace embedding problem on an input Vandermonde matrix A ∈ Rn×d can be reshaped as a constrained L1 subspace embedding problem on an input Vandermonde matrix of size n × (2r(d − 1) + 1). Thus, for Lp regression with p ∈ [2r, 2r+1), we have ∑
i∈[n] |〈ai,x〉|p = ∑ i∈[n] ∣∣∣∣ 2 r(d−1)+1∑ j=1 aj−1i,2 yj ∣∣∣∣p/2r , which is a constrained Lp/2r regression problem on an input Vandermonde matrix M of size n × (2r(d− 1) + 1). By Theorem A.2, we can use Lp/2r Lewis weight sampling to find a matrix M′ such that
1 5 ‖My‖p/2 r p/2r ≤ ‖M ′y‖p/2 r p/2r ≤ 5‖My‖ p/2r p/2r
for all y ∈ R2r(d−1)+1 with high probability. Note that by the above argument, if we take the matrix A′ corresponding to the scaled rows of A that are sampled by M′, then we also have
11 12 ‖Ax‖pp ≤ ‖A′x‖pp ≤ 13 12 ‖Ax‖pp
and thus 11
12 ‖Ax‖p ≤ ‖A′x‖p ≤
13 12 ‖Ax‖p
for all x ∈ Rd+1 with high probability. Thus by Lemma A.1, we can find a vector x̃ such that ‖Ax̃− b‖p ≤ 12 min
x∈Rd ‖Ax− b‖p.
We use the subroutine RoundTrunc to create the vector b′′ which is the vector with all entries of b′ rounded to the nearest power of (1 + ε), starting at the maximum entry of b′ in absolute value, and stopping after we are 1poly(n) times that, and replacing all remaining entries with 0. By the triangle inequality, we have
‖Ax− b′′‖p ≤ ‖Ax− b′‖p + ‖b′ − b′′‖p ≤ ‖Ax− b′‖p + 12εOPT ≤ (1 + 12ε)‖Ax− b′′‖p,
for any x ∈ Rd+1. Note that b′′ has discretized the values of b′ into ` = O (
logn ε
) possible values. We partition the
rows of A into ` groups G1, . . . , G`, based on the corresponding values of b′′. Suppose that for a group Gk, the corresponding values of b′′ are all tk. Then we have
‖Ax− b′′‖pp = ∑ k ∑ i∈Gk |〈ai,x〉 − tk|p .
Since A is Vandermonde, for each i ∈ Gk,
(〈ai,x〉 − tk)2 r = −tk + d∑ j=1 ai,jxj 2 r = −tk + d∑ j=1 aj−1i,2 xj 2 r = 22r(d−1)+1∑ j=1 aj−1i,2 tk,jyj ,
for some fixed values tk,1, tk,2, . . . that can be computed from tk, where again each yj can be a different function of the coordinates of x. In particular, there can be 2r choices for the exponent of tk. For a fixed choice of the exponent of tk, the exponent of ai,2 can range from 0 to 2r(d− 1).
Hence, the L2r regression problem on a submatrix of an input Vandermonde matrix A ∈ Rn×d with the same fixed measurement values, i.e., the corresponding coordinates of b′′ are all the same, can be reshaped to a constrainedL1 regression problem on an input Vandermonde matrix with 22r(d−1)+1 columns times a (22r(d− 1) + 1)× (22r(d− 1) + 1) diagonal matrix. Hence, by invoking Theorem A.2 to sample rows of M corresponding to their Lp/2r Lewis weights in the submatrix induced by the rows of Gk, we obtain a matrix Rk such that with high probability,
(1− ε)‖Rky − tk‖p/2 r p/2r ≤ ∑ i∈Gk 2r(d−1)+1∑ j=1 aj−1i,2 tk,jyj p/2 r ≤ (1 + ε)‖Rky − tk‖p/2 r p/2r .
Note that by the above argument, if we take the matrix Tk corresponding to the scaled rows of A that are sampled by Rk, then we also have
(1− ε)‖Tkx− tk‖pp ≤ ∑ i∈Gk |〈ai,x〉 − tk|p ≤ (1 + ε)‖Tkx− tk‖pp,
with high probability. Summing over all k, we have (1− ε) ∑ k ‖Tkx− tk‖pp ≤ ∑ k ∑ i∈Gk |〈ai,x〉 − tk|p ≤ (1 + ε) ∑ k ‖Tkx− tk‖pp,
with high probability. Let T = T1 ◦ . . . ◦ T`, so that ∑ k ‖Tkx − tk‖pp = ‖Tx − b′′‖pp. Observe that ‖Ax − b′′‖pp =∑
k ∑ i∈Gk |〈ai,x〉 − tk| p. Hence,
(1− ε)‖Tx− b′′‖pp ≤ ‖Ax− b′′‖pp ≤ (1 + ε)‖Tx− b′′‖pp and thus (1− ε)‖Tx− b′′‖p ≤ ‖Ax− b′′‖p ≤ (1 + ε)‖Tx− b′′‖p. Thus we can compute a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 +O (ε)) min x∈Rd ‖Ax− b‖p.
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-Vector Multiplication Runtime, e.g., Table 1 in Gohberg & Olshevsky (1994)). The runtime of computing Ax for a Vandermonde matrix A ∈ Rn×d and a vector x ∈ Rd for d ≤ n is O ( n log2 n ) .
Theorem A.5 (Approximate Lp Regression Runtime). (Adil et al., 2019b) Given A ∈ Rn×d and p ≥ [2,∞), there exists an algorithm that makes O (√ n log nε ) calls to a linear system solver and computes a vector x̃ such that
‖Ax̃− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p.
We now analyze the runtime of Algorithm 2. Lemma A.6. Algorithm 2 runs in O ( n log2 n+ p 3d3/2
ε7/2 log4 n log 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 to explicitly form them, which would otherwise require Ω(ndp2) time just to 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 on the final subsampled matrix. The only remaining procedure is rounding and truncating the coordinates of b ∈ Rn to form a vector b′′ ∈ Rn using the procedure RoundTrunc and then forming the groupsG1, . . . , G`, which clearly takesO (n) arithmetic operations combined. We thus analyze each of the three main runtime bottlenecks.
First observe that the extended matrix M is a Vandermonde matrix with O ( dp2 )
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 Lp Lewis weight sampling. By Theorem A.4, each matrix-vector multiplication uses time O ( n log2 n ) . Hence computing the extended matrix M uses
O ( n log2 n ) time. Similarly, the extended matrix for each group Gk is the product of a Vander-
monde matrix with O ( dp2 )
columns and a diagonal matrix. Thus by Theorem A.4, the extended matrices for all the groups Gk can be formed using O ( n log2 n ) time in total.
Since each Vandermonde matrix has O ( dp2 ) rows, then observe that each group Gk samples
O ( dp2 ε2 log d ) rows and there are ` = O ( 1 ε log n ) such groups k ∈ [`]. Thus the resulting sub-
sampled matrix has O ( dp2 ε3 log 2 n ) rows for d ≤ n. To approximately solve the Lp regression
problem, Theorem A.5 notes that for p ≥ 2 and a subsampled matrix of size O ( dp2 ε3 log 2 d ) , we
require only O ( p √ d
ε3/2 log n log dε
) calls to a linear system solver. Moreover, on an iteration t of the
Lp regression algorithm of Adil et al. (2019a) used in Theorem A.5, the linear system solves the equation xt ← (A>RTA)−1A>Rtb for a diagonal matrix RT . Thus by Theorem A.4, each linear system solve can be done in O ( dp2 ε3 log 4 n ) time for a Vandermonde matrix A. Hence, the total
time to approximately solve the Lp regression problem is O ( p3d3/2 ε7/2 log4 n log nε ) . Therefore, the
total runtime is O ( n log2 n+ p 3d3/2
ε7/2 log4 n log nε
) .
Lemma 2.4. Let x ∈ Rn and p = Ω (
logn ε ) . Then ‖x‖∞ ≤ ‖x‖p ≤ (1 + ε)‖x‖∞.
Proof. For any vector x ∈ Rd, we have
‖x‖∞ = max i∈n |xi| ≤ (∑ i∈ |xi|p )1/p ≤ n1/p ·max i∈n |xi|.
Since (1 + ε)3/eps > e for all ε > 0, then (1 + ε)p > n for p = Ω (
logn ε
) . Therefore,
n1/p ·max i∈n |xi| ≤ (1 + ε) max i∈n |xi| = (1 + ε)‖x‖∞.
Corollary 2.5. Given ε ∈ (0, 1) and p ≥ 1, A ∈ Rn×dq , and b ∈ Rd, suppose A = [A1| . . . |Aq] for Vandermonde matrices A1, . . . ,Aq ∈ Rn×d. Then there exists an algorithm that, with high probability, returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p,
using O ( T (A)(dp)q−1 log n+ poly((dp)q log n, 1/ε) ) time, where T (A) is the runtime of mul-
tiplying the matrix A by an arbitrary vector. For q = 2, this can be further optimized to O ( ndω2/2−1 + poly((dp)2, log n, 1/ε) ) time, whereO (nω2) is the time to multiply an n×nmatrix
with an n× n2 matrix, so that ω2 ∈ [3, 4].
Proof. Recall that a key part in the proof of Theorem 1.1 was to first the L2r regression on a Vandermonde matrix with dimension Rn×d to L1 regression on a Vandermonde matrix with dimension Rn×(2r(d−1)+1), where 2r ≤ p < 2r+1. For a matrix A with block Vandermonde structure, we can similarly write
(〈ai,x〉)2 r =  dq∑ j=1 ai,jxj 2 r =  q∑ k=1 d∑ j=1 aj−1i,2+(k−1)dxj+(k−1)d 2 r
= (2r(d−1))q+1∑ j=1 q∏ k∈[q],∑ pj,k=2r a pj,k i,2+(k−1)dyj ,
where again each yj is a fixed function of the coordinates of x. Thus we can reshape the L2r regression problem on a matrix A with dimension Rn×dq with block Vandermonde structure to an L1 regression problem on a matrix Ã with dimension Rn×(2
r(d−1))q+1. Moreover, we can further reshape Ã into the concatenation of (dp)q−1 Vandermonde matrices, where each Vandermonde matrix has columns that are geometrically growing in ai,2 but are multiplied by all (dp)q−1 products∏q−1 k=1 a pk i,2+kd, where pk ∈ [dp].
We can now use matrix-vector multiplication on each of the (dp)q−1 Vandermonde matrices. Thus by Theorem A.2, we can L1 Lewis weight sample from the rows of the reshaped Ã, using O ( T (A)(dp)q−1 log n+ (dp)ωq ) time. We can similarly write (〈ai,x− tk〉)2 r
for each tk among the discretized values of the updated b vector as a sum of (2r(d − 1))q + 1 terms that are all products of powers of the bases ai,2, ai,2+d, . . . and a variables yj , as in the proof of Theorem 1.1.
Thus we can partition the L2r regression problem into ` = O ( logn ε ) instances of a constrained L1 regression problem on (dp)q−1 Vandermonde matrices, each with at most 2r(d − 1))q + 1 columns. To approximately solve the Lp regression problem, we can thus sample rows by their Lp/2r Lewis weights, as in Theorem 1.1. Since there are up to (dp)q−1 Vandermonde matrices, each with at most 2r(d − 1))q + 1 columns, then by Theorem A.2, the total time required is O ( T (A)(dp)q−1 log n+ (dp)ωq · poly(log n, 1/ε) ) .
We remark that for the special case of q = 2, Sa et al. (2018) noted an efficient bivariate matrix multiplication algorithm of Nüsken & Ziegler (2004); Kedlaya & Umans (2011).
Theorem A.7. (Nüsken & Ziegler, 2004; Kedlaya & Umans, 2011; Sa et al., 2018) Given a q-variate polynomial f(X1, . . . , Xq) such that each variable has degree at most d − 1 and N = dq distinct points x(i) = (x(i)1, . . . , x(i)q) for i ∈ [N ], there exists an algorithm that uses O ( dω2(q−1)/2+1 ) time to output the vector (f(x(1)), . . . , f(x(N))), where O (nω2) is the time to multiply an n × n matrix with an n× n2 matrix, so that ω2 ∈ [3, 4].
The case of a matrix-vector product for q = 2 corresponds to the evaluation of d2 points in Theorem A.7. Thus we need to repeat the algorithm in Theorem A.7 a total of nd2 times to handle all n rows in the input matrix. Since each instance of the algorithm uses O ( dw2/2+1 ) time, the total
time for the matrix-vector product is O ( ndω2/2−1 ) , rather than the naı̈ve O (nd) time (recall that ω2 ∈ [3, 4]).
Corollary A.8. Given ε ∈ (0, 1) and p ≥ 1, A ∈ Rn×2d, and b ∈ Rd, suppose A = [A1|A2] for Vandermonde matrices A1,A2 ∈ Rn×d. Then there exists an algorithm that with high probability returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p, using O ( ndω2/2−1 + poly((dp)2, log n, 1/ε) ) , where O (nω2) time.
A.2 Lp REGESSION ON OTHER STRUCTURED INPUT
Algorithm 4 Faster regression for noisy low-rank matrices
Input: Sparse matrix A ∈ Rn×d, measurement vector b, accuracy parameter ε > 0 Output: x̂ ∈ R with ‖Ax̂− b‖p ≤ (1 +O (ε)) minx∈Rd ‖Ax− b‖p
1: r ← blog pc . 2r ≤ p < 2r+1. 2: Extend A to a matrix M with dimension n × O ( ds(d+ s)2 r) so that each entry Mi,j is the
coefficient of the j-th term in the tensor decomposition of a⊗2 r
i . 3: UseLp/2r -Lewis weight sampling on M to find a set S ofO (d 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 minx∈Rd ‖A′x− b‖p. 6: b′ ← b−Ax̃, b′′ ← RoundTrunc(b′), `← O ( logn ε
) 7: Partition the rows of A into groups G1, . . . , G`, each containing all rows with the same value
of b′′ 8: Let Gk be the corresponding submatrix and tk be the coordinate of b′′ corresponding to Gk for
each k ∈ [`]. 9: Use Lp/2r -Lewis weight sampling on [Gk; tk] to find a set S′k ofO ( 1 εd log d ) indices in [n] and
rescaling factors. 10: Let Tk be the corresponding rows with indices in S′k and scaled accordingly. 11: T← [T1; . . . ;Tk]> 12: Compute x̂ ≤ (1 + ε) minx∈Rd ‖Tx− b‖p. 13: return x̂
Lemma A.9. Given ε ∈ (0, 1) and p ≥ 1, a matrix A ∈ Rn×d such that A = K + S for a rank k matrix K and an s-sparse matrix S, and b ∈ Rd, there exists an algorithm that with high probability, returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p.
Proof. The proof is similar to Lemma A.3. We once again let r be an integer so that 2r ≤ p < 2r+1 and observe that∑
i∈[n] |〈ai,x〉|p = ∑ i∈[n] |(〈ai,x〉)2 r |p/2 r = ∑ i∈[n] |(〈ki + si,x〉)2 r |p/2 r .
Since K is a low-rank matrix, then for all ki, we can write
ki = k∑ j=1 αi,jvj ,
for a fixed set of basis vectors v1, . . . ,vk ∈ Rd. Hence we have
∑ i∈[n] |〈ai,x〉|p = ∑ i∈[n] ∣∣∣∣∣∣ 〈 si + k∑ j=1 αi,jvj ,x 〉2r ∣∣∣∣∣∣ p/2r .
Since S has sparsity s, then we can further write
∑ i∈[n] |〈ai,x〉|p = ∑ i∈[n] ∣∣∣∣∣∣ 〈 s∑ j=1 βi,jeij + k∑ j=1 αi,jvj ,x 〉2r ∣∣∣∣∣∣ p/2r ,
where e1, . . . , ed denote the elementary vectors. By the Hadamard Product-Kronecker Product mixed-product property, we have∑
i∈[n] |〈ai,x〉|p = ∑ i∈[n] ∣∣∣(y1 ⊗ . . .y2r ) x⊗(2r)∣∣∣p/2r , where each yk ∈ {αi,1v1, . . . , αi,kvk, βi,1ei1 , . . . , βi,sei,s} for i ∈ [2r]. Thus for a fixed set of elementary vectors ei1 , . . . , eis , there are (d + s)
2r possible values for the tensor product (y1 ⊗ . . .y2r ). Since there are ( d s ) choices for the elementary vectors ei1 , . . . , eis , then there are at most (Cds(d+s)2 r
) possible values for the tensor product for an absolute constant C > 0. Therefore, the Lp subspace embedding problem on A ∈ Rn×d can be reshaped as a constrained Lp/2r subspace embedding problem on an input matrix of size n × (Cds(d + s)2r ). Hence for Lp regression with p ∈ [2r, 2r+1), we have
∑ i∈[n] |〈ai,x〉|p = ∑ i∈[n] ∣∣∣∣ 2 r(d−1)+1∑ j=1 Mi,jyj ∣∣∣∣p/2r , which is a constrained Lp/2r regression problem on a matrix M of size n× (Cds(d+ s)2 r
) whose entries can be determined from the decomposition of each row of A.
Using Lp/2r Lewis weight sampling, Theorem A.2 implies that we can find a matrix M′ such that
11 12 ‖My‖p/2 r p/2r ≤ ‖M ′y‖p/2 r p/2r ≤ 13 12 ‖My‖p/2 r p/2r
for all y ∈ R(Cds(d+s)2 r ) with high probability. By Lemma A.1, we can thus compute a vector x̃ such that ‖Ax̃− b‖p ≤ 12 min
x∈Rd ‖Ax− b‖p.
We again set b′ = b−Ax̃ and define OPT = minx∈Rd ‖Ax− b‖p, so that
‖b′‖p = ‖Ax̃− b‖p ≤ 12OPT.
Let b′ = b−Ax̃ and OPT = minx∈Rd ‖Ax− b‖p, so that
‖b′‖p = ‖Ax̃− b‖p ≤ 12OPT.
Let b′′ be the vector with all entries of b′ rounded to the nearest power of (1 + ε), starting at the maximum entry of b′ in absolute value, and stopping after we are 1poly(n) of that and replacing all remaining entries with 0. By the triangle inequality, we have
‖Ax− b′′‖p ≤ ‖Ax− b′‖p + ‖b′ − b′′‖p ≤ ‖Ax− b′‖p + 12εOPT ≤ (1 + 12ε)‖Ax− b′′‖p,
for any x ∈ Rd+1. Since the coordinates of b′′ can have ` = O (
logn ε
) possible distinct values, we can partition the
rows of A into ` groups, G1, . . . , G`, based on the corresponding values of b′′. Let tm be the corresponding value of b′′ for all rows in a group Gm, so that
‖Ax− b′′‖pp = ∑ m ∑ i∈Gk |〈ai,x〉 − tm|p.
By the above argument, we have for each i ∈ Gk,
(〈ai,x〉 − tm)2 r = ∣∣∣∣∣∣−tm + 〈 s∑ j=1 βi,jeij + k∑ j=1 αi,jvj ,x 〉∣∣∣∣∣∣ 2r = 2rCds(d+s)2 r∑ j=1 Bi,jtm,jyj ,
where (1)Bi,j are entries of a matrix B with 2rCds(d+s)2 r
columns that can be computed from A, (2) tm,1, tm,2, . . . are fixed values that can be computed from mk, and (3) each yj is a fixed function of the coordinates of x. Notably, B is the matrix formed by the concatenation of the coefficients of
the decomposition of the α-fold tensor product of the row into the α-fold tensor products of the lowrank and sparse basis elements, for each α = 0, . . . , p. By comparison, the matrix M previously defined in this proof is only the decomposition for α = p. Hence, the L2r regression problem on a submatrix of A ∈ Rn×d the same coordinate of b′′, can be reshaped as a constrained L1 regression problem on a matrix with 2rCds(d+ s)2 r columns.
The remainder of the proof follows from the same grouping argument as Theorem 1.1. We apply Theorem A.2 by sampling rows of B corresponding to their Lp/2r Lewis weights in the submatrix induced by the rows of Gm, we obtain a matrix Rm such that with high probability,
(1− ε)‖Rmy − tm‖p/2 r p/2r ≤ ∑ i∈Gm 2rCds(d+s)2r∑ j=1 Bi,jtm,jyj p/2 r ≤ (1 + ε)‖Rmy − tm‖p/2 r p/2r .
Conditioning on the above inequality holding, then by taking the matrix Tm corresponding to the scaled rows of A that are sampled by Rm, it follows that
(1− ε)‖Tmx− tm‖pp ≤ ∑ i∈Gm |〈ai,x〉 − tm|p ≤ (1 + ε)‖Tmx− tm‖pp.
Therefore by summing over all m ∈ [`], we have that with high probability, (1− ε) ∑ m ‖Tmx− tm‖pp ≤ ∑ m ∑ i∈Gm |〈ai,x〉 − tm|p ≤ (1 + ε) ∑ m ‖Tmx− tm‖pp.
For T = T1 ◦ . . . ◦ T`, we have ∑ m ‖Tmx − tm‖pp = ‖Tx − b′′‖pp. Since ‖Ax − b′′‖pp =∑
m ∑ i∈Gm |〈ai,x〉 − tm| p, then
(1− ε)‖Tx− b′′‖pp ≤ ‖Ax− b′′‖pp ≤ (1 + ε)‖Tx− b′′‖pp. Therefore for p ≥ 1,
(1− ε)‖Tx− b′′‖p ≤ ‖Ax− b′′‖p ≤ (1 + ε)‖Tx− b′′‖p. Because ` = O (
logn ε
) and ∑ O (T (Gk)) = O (T (A)), we can compute a vector x̂ ∈ Rd such
that ‖Ax̂− b‖p ≤ (1 +O (ε)) min
x∈Rd ‖Ax− b‖p.
Lemma A.10. The runtime of Algorithm 4 is O (pnpk + pnsp) + poly ( ds, kp, sp, 1ε ) .
Proof. We analyze the runtime of Algorithm 4. First note that we can compute the extended matrix in time O (n((k + 1)pdssp)) to perform the Lp/2r Lewis weight sampling, where we recall that r is the unique integer such that 2r ≤ p < 2r+1. To perform Lewis weight sampling on the extended matrix, we require matrix-vector multiplication, which requires time O (nk) for a low-rank matrix and time O (ns) for a matrix whose rows have at most s nonzero entries. After the first iteration of Lewis weight sampling, we can round and truncate the coordinates of b ∈ Rn to form a vector b′′ ∈ Rn using the procedure RoundTrunc and then forming the groups G1, . . . , G`, which clearly takes O (n) arithmetic operations combined. Once the groups are formed, we can compute the extended matrix in time O (n((k + 1)pdssp)) and perform Lp/2r Lewis weight sampling on each group, which takesO (nk + ns) total time across all groups. To approximately solve the resulting Lp regression problem formed by the subsampled rows, we require poly ( ds, kp, sp, 1ε ) time. Therefore, the total runtime is O (nk + ns) + poly ( ds, kp, sp, 1ε ) .
We obtain similar results for a noisy Vandermonde matrix. Lemma A.11. Given ε ∈ (0, 1) and p ≥ 1, a matrix A ∈ Rn×d such that A = V + S for a Vandermonde matrix V and an s-sparse matrix S, and b ∈ Rd, there exists an algorithm that with high probability, returns a vector x̂ ∈ Rd such that
‖Ax̂− b‖p ≤ (1 + ε) min x∈Rd ‖Ax− b‖p. Lemma A.12. The runtime of Algorithm 4 is O (nk + ns) + poly ( ds, pd, sp, 1ε ) .
Theorem 1.6 then follows from Lemma A.11 and Lemma A.12.

heading: B APPLICATIONS TO POLYNOMIAL REGRESSION
text: In the polynomial regression problem, the goal is to find a degree d polynomial q̂ such that
‖q̂(t)− f(t)‖pp ≤ (1 + ε) · min q:deg(q)≤d ‖q(t)− f(t)‖pp,
where ε > 0 is an accuracy parameter given as input and ‖·‖pp is the Lp norm to the pth power, ‖g‖pp = ∫ 1 −1 |g(t)|
pdt. The polynomial regression problem is a fundamental problem in statistics, computational mathematics, machine learning, and more. The problem has been studied as 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 et al. (2008), solving parametric PDEs Hampton & Doostan (2015), and surface reconstruction Pratt (1987).
Given the flexibility to choose query locations x1, . . . , xs, we can consider the polynomial regression problem as an active learning or experimental design problem. Thus we would like to minimize the number of queries s, as a function of the approximation degree d, the norm p, and the accuracy parameter ε, to find q̂. Observe that d+1 queries are obviously necessary, but also that d+1 queries suffice when f can be exactly fit by a degree d polynomial, by using direct interpolation. In general, however, in the case when minq:deg(q)≤d ‖q(t) − f(t)‖pp 6= 0 we require s > d + 1 queries. Our Vandermonde Lp regression results can be used to give the first result showing that for all p ≥ 1, d exp(O(p)) ·poly ( 1 ε ) queries suffice to obtain a (1 + ε)-approximation to the best polynomial fit.
We require the following structural theorem reducing the Lp polynomial regression problem to a problem of solving Lp regression on Vandermonde matrices: Theorem B.1. Kane et al. (2017); Meyer et al. (2021) Suppose s1, . . . , sn0 are drawn uniformly from [−1, 1]. Let A ∈ Rn0×(d+1) be the associated Vandermonde matrix, so that Ai,j = sj−1i . Let n0 = exp(O(p))Õ ( 1 ε2+2p d 5 )
and let b ∈ Rn0 be the evaluations of f , so that bi = f(si). Then with probability 1112 , the sketched solution x̂ = argminx ‖Ax− b‖p satisfies
‖Px− f‖pp ≤ (1 + ε) min x∈Rd+1 ‖Px− f‖pp.
Theorem B.1 states that we can uniformly sample exp(O(p)) · poly ( d, 1ε ) points from [−1, 1]. We can then form an Lp regression problem by using the evaluation each of the polynomial bases at the sampled points to 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 x̂ = argminx ‖Ax−b‖p is a (1 + ε)-approximation to the best fit degree d polynomial. We can naı̈vely approximately solve the Lp regression on the Vandermonde matrix A and the measurement vector b by standard Lp regression techniques such as Lewis weight sampling, which would result in total query complexity exp(O(p))Õ ( 1 ε2+2p d 5 ) . However, we can instead note that Lemma 1.4 implies we can instead solve an Lq regression problem on a Vandermonde matrix with O(dp) columns for q ∈ [1, 2]. Crucially for q ∈ [1, 2], there exist active Lq regression algorithms that only require reading Õ(d) · poly ( 1 ε ) entries of b: Theorem B.2. Musco et al. (2021b) Given p ∈ [1, 2] and an input matrix A ∈ Rn×d, there exists an algorithm that reads Õ(d) · poly ( 1 ε ) entries of b and with probability at least 0.99, outputs x̃ such that ‖Ax− b‖p ≤ (1 + ε) min
x ‖Ax− b‖p.
Hence by Theorem B.2, we can approximately solve the Lq regression on a Vandermonde matrix with dp rows by reading Õ(dp) · poly ( 1 ε ) entries of b. By Lemma 1.4, the approximate solution will also be a (1+ε)-approximation to the optimal Lp regression on the Vandermonde matrix A. By Theorem B.1, the approximate solution will also form the coefficient vector of a polynomial that is a (1+ε)-approximation to the polynomial Lp regression problem. Therefore, we obtain the following guarantees for the polynomial regression problem: Theorem B.3. For any degree d and norm p ≥ 1, there exists an algorithm that queries f at s = dp poly ( log(dp), 1ε ) points and outputs a degree d polynomial q̂(t) such that
‖q̂(t)− f(t)‖pp ≤ (1 + ε) · min q:deg(q)≤d ‖q(t)− f(t)‖pp,
with probability at least 23 .
The previous-best algorithm Kane et al. (2017); Meyer et al. (2021) for Lp polynomial regression sampled O(d5) points uniform from [−1, 1] and then used standard Lp regression algorithms that required reading the signal at all O(d5) sampled points, for a total query complexity of O(d5). By comparison, Theorem B.3 only has linear dependency in d due to the structural property of Lemma 1.4. Since Ω(d) queries are clearly necessary, our result settles the dependency of d in the query complexity for Lp polynomial regression for all d and all p ≥ 1.

heading: C ADDITIONAL EXPERIMENTS
text: In this section, we experimentally verify that Lp/2r Lewis Weight sampling works for unstructured matrices. We take a similar approach as in Section 3 to verify that Lp/2r Lewis Weight sampling is correct for Lp regression on unstructured matrices. For this test, we fix n = 25, 000, d = 10, and p = 6, while varying m ∈ [100, 3000]. We let A = [ G1 0 0 G2 ] where G1 ∈ R50×4 and G2 ∈ R24,950×6 are i.i.d. N(0, 1) matrices. The first 50 entries of b are sampled from N(0, 1002), while the rest of the entries are N(0, 1).
To compute the Lewis Weights of the expanded matrix M, due to numerical instability, we cannot apply Lewis Weight iteration directly. Instead, we compute 10 steps of Lewis Weight Iteration on the random sketch MG, where G ∈ R(dp/2×100) has i.i.d. N(0, 1) entries. We then subsample A and b with respect to these Lewis Weights, and directly solve Lp regression on the subsampled matrices. 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. For instance, to achieve 2.5% error, takingm = 500, or just 2% of the rows of A suffices. Also note 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 a Gaussian sketch matrix and the standard Lewis Weight Iteration, which converges very quickly.
Since b is so large on its first 50 entries, it is important for any subsampling algorithm to sample at least 4 of the first 50 rows. Uniform sampling picks few (if any) of these rows, which is why uniform sampling fails to converge to a good solution, even for large m. Lewis weight sampling instead gives much higher priority to the first 50 rows, avoiding any issue. This is why the the gap between Lewis Weight sampling and uniform sampling is so large for this experiment.

heading: D ADDITIONAL RELATED WORK
text: As previously mentioned, subspace embeddings are common tools used to approximately solve Lp regression. Given an input matrix A ∈ Rn×d, a subspace embedding is a matrix M ∈ Rm×d with
m n such that (1− ε)‖Ax‖p ≤ ‖Mx‖p ≤ (1 + ε)‖Ax‖p,
for all x ∈ Rd. Thus given an instance of Lp regression, where the goal is to minimize ‖Ax− b‖p, we can set B = [A;b], compute a subspace embedding for B, and then solve a constrained Lp regression problem on the smaller matrix B.
A subspace embedding of a matrix A ∈ Rn×d can be formed by sampling rows of A with probabilities proportional to their Lp leverage scores, e.g., (Cormode et al., 2018; Dasgupta et al., 2008) their Lp sensitivities, e.g., (Clarkson et al., 2019; Braverman et al., 2020; Musco et al., 2021a); or their Lp Lewis weights, e.g., (Cohen & Peng, 2015; Durfee et al., 2018; Chhaya et al., 2020; Chen & Derezinski, 2021; Parulekar et al., 2021; Meyer et al., 2021). In any of these cases, the resulting matrix M will contain a subset of rows of A that are rescaled by a function of their sampling probability, as to give an unbiased estimate of the actual Lp 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 = RA for some (often random) matrix R ∈ Rm×n. 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 (Lp leverage scores, Lp sensitivities, and Lp Lewis weights) are data dependent and thus require a pass over the matrix A. The sketching matrix can be generated from a family of random matrices whose entries are Cauchy random variables for p = 1, e.g., Sohler & Woodruff (2011); Meng & Mahoney (2013); Clarkson et al. (2016), or sub-Gaussian random variables for p = 2, e.g., Sarlós (2006); Nelson & Nguyen (2013); Clarkson & Woodruff (2013). More generally, exponential random variables can be used for p ≥ 2, though the number of rows m in the resulting sketch matrix now has a polynomial dependency in n Woodruff & Zhang (2013). A line of recent works has studied the tradeoffs between oblivious linear sketches, sampling-based algorithms, and other sketches for Lp subspace embeddings Wang & Woodruff (2019); Li et al. (2021). In this paper, we focus on sampling-based algorithms for Lp regression due to preservation of structure when the input matrix A itself has structure.

Review summaries:
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.  The 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.

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. 

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.   This 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.  For 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.  Technically 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."  There 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).