Provable Approximations for Constrained $\ell_p$ Regression
Ibrahim Jubran David Cohn Dan Feldman

arabic
plain
[
Provable Approximations for Constrained ℓ_p Regression














equal*


Ibrahim Jubranuni
David Cohnuni
Dan Feldmanuni


uniRobotics and Big Data Lab, Department of Computer Science, University of Haifa, Haifa, Israel

Ibrahim Jubranibrahim.jub@gmail.com
David Cohndavid.cohn.82@gmail.com




Machine Learning, Non-Convex Optimization, Linear Regression, Approximation Algorithms

0.3in
]









  




The ℓ_p linear regression problem is to minimize f(x)=||Ax-b||_p over x∈ℝ^d, where A∈ℝ^n× d, b∈ℝ^n, and p>0.
To avoid overfitting and bound ||x||_2, the constrained ℓ_p regression minimizes f(x) over every unit vector x∈ℝ^d.
This makes the problem non-convex even for the simplest case d=p=2.
Instead, ridge regression is used to minimize the Lagrange form f(x)+λ ||x||_2 over x∈ℝ^d, which yields a convex problem in the price of calibrating the regularization parameter λ>0.

We provide the first provable constant factor approximation algorithm that solves the constrained ℓ_p regression directly, for every constant p,d≥ 1.
Using core-sets, its running time is O(n log n) including extensions for streaming and distributed (big) data.
In polynomial time, it can handle outliers, p∈ (0,1) and minimize f(x) over every x and permutation of rows in A.

Experimental results are also provided, including open source and comparison to existing software.




§ INTRODUCTION
 
One of the fundamental problems in machine learning is ℓ_2 linear regression, where the goal is to fit a hyperplane that minimizes the sum of squared vertical distances to a set of n≫ d input (d+1)-dimensional points (samples, vectors, training data). Formally, the input is an n× d matrix A = (a_1|⋯| a_n)^T and a vector b=(b_1,⋯,b_n)^T in ^n that contains the n labels (heights, or last dimension) of the points. The goal is to minimize the sum ∑_i=1^n(a_i^Tx-b_i)^2 over every d-dimensional vector x of coefficients,

    min_x∈^dAx-b_2.


One disadvantage of these techniques is that overfitting may occur <cit.>. For example, if the entries in b are relatively small, then x=(0,⋯,0) may give an approximated but numerically unstable solution. Moreover, (<ref>) is not robust to outliers, in the sense that e.g. adding a row whose entries are relatively very large would completely corrupt the desired vector x.

This motivates the addition of a constraint c>0 on the norm of x, where c is a constant that may depend on the scale of the input A and b. Without loss of generality, we can assume x_2=c=1, otherwise we divide entries ofP b accordingly. The result is constrained ℓ_2 regression problem,

    min_x∈^d:x_2=1Ax-b_2.

The special case b=(0,⋯,0) can be solved in O(nd^2) time, where the optimum is the smallest singular value of A and x is the largest singular vector of x <cit.>.

A generalization of (<ref>) for a given constant p>0 would be

    min_x∈^d, x_2=1Ax-b_p,

where v_p=(∑_i=1^d |v_i|^p)^1/p for v = (v_1,⋯,v_d)∈^d. Note that for p<1 we obtain a non-standard norm which is a non-convex function over v∈ R^d.

Optimization problem (<ref>) for p=1 can be defined geometrically as follows.
Compute a point x' on the unit sphere that minimizes the weighted sum of distances over n given hyperplanes and n multiplicative weights.
Here, the ith hyperplane is defined by its normal (unit vector) a_i/a_i, its distance from the origin b_i/a_i, and its weight a_i. The weighted distance between x' and the ith hyperplane is defined as a_i· |a_i^T/a_ix' -b_i/a_i| = |a_i^Tx-b_i|.

In the context of machine learning, in linear regression we wish to fit a hyperplane whose unit normal is x', that minimizes the sum of squared vertical distances between the hyperplane at point a_i (predicted value) and b_i (the actual value), over every i∈ [n]. In low-rank approximation (such as SVD / PCA) we wish to fit a hyperplane that passes through the origin and whose unit normal is x', that minimizes the sum of squared Euclidean distances between the data points and the hyperplane.
Our problem is a mixture of these two problems: compute a hyperplane that passes through the origin (as in low-rank approximation) and minimizes sum of squared vertical distances (as in linear regression).

Further generalization of (<ref>) suggests handling data with outliers. For example when one of the rows of A is very noisy, or if an entry of b is unknown. Let k < n be the number of such outliers in our data. In this case, we wish to ignore the largest k distances (fitting errors), i.e., consider only the closest s=n-k points to x. Formally,

    min_x∈^d:x_2=1(Ax-b,n-k)_p,

where (v,s) ∈^s is a vector that consists of the smallest s entries in v∈^n, where s∈[0,n] is an integer.

In some cases, our set of observations is unordered, i.e., we do not know which observation in b matches each point in A.
For example, when half of the points should be assigned to class b_1=..=b_n/2=0 and half of the points to class b_n/2+1..=b_n=1. Here, :1,⋯,n→1,⋯,n denotes a bijective function, called a matching function, and b_ = (b_(1),⋯,b_(n))^T denotes the permutation of the entries in b with respect to . In this case, we need to compute

    min_x,Ax-b__p,

where the minimum is over every unit vector x and matching function :[n]→ [n].





§ RELATED WORK
 
Regression problems are fundamental in statistical data analysis and have numerous applications
in applied mathematics, data mining, and machine learning; see references in <cit.>. Computing the simple (unconstrained) Linear regression in (<ref>) for the case p=2 was known already in the beginning of the previous century <cit.>. Since the ℓ_p norm is a convex function, for p=1 it can be solved using linear programming, and in general for p≥ 1 using convex optimization techniques in time (nd)^O(1) or using recent coreset (data summarization) technique <cit.> in near-linear time.

To avoid overfitting and noise, there is a need to bound the norm of the solution x, which yields problem (<ref>) when p=2.
The constraint x=1 in (<ref>) can be replaced by adding a Lagrange multiplier <cit.> λ>0 to obtain,

    min_λ∈,x∈^dAx-b_2+λx_2-λ.


Unfortunately, (<ref>) and (<ref>) are non-convex problems in quadratic programming as explained in <cit.>, which are also NP-hard if d=n, so there is no hope for running time that is polynomial in d; see Conclusion section. It was proved in <cit.> that the problem is non-convex even if every input point (row) a_i is on the unit circle and d=2.




Similarly, when we are allowed to ignore k outlier as defined in (<ref>), or if we use M-estimators, the problem is no longer convex.


Instead, a common leeway is to “guess" the value of λ in (<ref>), i.e., turn it into an input parameter that is calibrated by the user and is called regularization term <cit.> to obtain a relaxed convex version of the problem,

    min_x∈^dAx-b_2+λx_2.

This problem is the ridge regression which is also called Tikhonov regularization in statistics <cit.>, weight decay in machine learning <cit.>, and constrained linear inversion method <cit.> in optimization. Many heuristics were suggested to calibrate λ automatically in order to remove it such as automatic plug-in estimation, cross-validation, information
criteria optimization, or Markov chain Monte Carlo (MCMC) <cit.> but no provable approximations for the constrained ℓ_p regression (<ref>) are known; see <cit.> and references therein.




Another reminiscent approach is LASSO (least absolute shrinkage and selection operator) <cit.>, which replaces the non-convex constraint x_2=1 in (<ref>) with its convex ℓ_1 inequality x_1≤ t to obtain min_x∈^d:x_1≤ tAx-b_2 for some given parameter t>0.

LASSO is most common technique in regression analysis <cit.> for e.g. variable selection and compressed sensing <cit.> to obtain sparse solutions. As explained in <cit.> these optimization problems can be easily extended to a wide variety of statistical models including generalized linear models, generalized estimating equations, proportional hazards models, and M-estimators. 

Alternatively, the ℓ_2-norm in (<ref>) may be replaced by the ℓ_p-norm for p≥ 1 to obtain the (non-constrained) ℓ_p regression  
    min_x∈^dAx-b_p,

which is convex for the case p≥ 1. Using p∈(0,1] in (<ref>) is especially useful for handling outliers <cit.> which arise in real-world data. However, for the (non-standard) ℓ_p-norm, where p<1, (<ref>) is non-convex.


Adding the constraint x=1 in (<ref>) yields the constrained ℓ_p regression in (<ref>). Only recently, a pair of breakthrough results were suggested for solving (<ref>) if p≠ 2. <cit.> suggested a solution to the constraint ℓ_p regression in (<ref>) for the case p=1. They suggest to convert the constraint x=1 into two inequality constraints x≤ 1 and -x≤ 1. The other result <cit.> suggested a provable constant approximation for the constrained ℓ_p regression problem, in time O(nlog n) for every constant p≥1. However, the result holds only for d=2, and the case d>2 as in our
 paper was left as an open problem. In fact, our main algorithm solves the problem recursively where in the base case d=2 we use the result from <cit.>.

To our knowledge, no existing provable approximation algorithms are known for handling outliers as in (<ref>), unknown matching as in (<ref>), or for (<ref>) for the case p∈ (0,1) and d≥ 3.



  
Coreset for ℓ_p regression in this paper is a small weighted subset of the input that approximates Ax-b for every x∈^d, up to a multiplicative factor of 1±. Solving the constrained ℓ_p regression on such coreset would thus yield an approximation solution to the original (large) data. Such coresets of size independent of n were suggested in <cit.>. In Theorem <ref> we obtain a little smaller coreset by combining <cit.> and the framework from <cit.>. Such coresets can also be maintained for streaming and distributed Big Data in time that is near-logarithmic in n per point as explained e.g. in <cit.>. Applying our main result on this coreset, thus implies its streaming and distributed versions.
We note that this scenario is rare and opposite to the common case: in most coreset related papers, a solution for the problem that takes polynomial time n^O(1) exists and the challenge is to reduces its running time to linear in n, by applying it on a coreset of size that is independent, or at least sub-linear in n. In our case, the coreset exists but not a polynomial time algorithm to apply on the coreset.



§ PAPER OVERVIEW

The rest of the paper is organized as follows. We state our main contributions in Section <ref> and some preliminaries and notations in Section <ref>. Section <ref> suggests an approximation algorithm for solving (<ref>), Section <ref> generalizes this solution of (<ref>) to a wider range of functions, Section <ref> handles the minimization problem in (<ref>), Section <ref> introduces a coreset for the constrained ℓ_p regression, Section <ref> presents our experimental results and Section <ref> concludes our work.



§ OUR CONTRIBUTION
 
Some of the proofs have been placed in the appendix to make the reading of the paper more clear.



  
Constrained ℓ_p regression. We provide the first polynomial time algorithms that approximates, with provable guarantees, the functions in (<ref>), (<ref>) and (<ref>) up to some constant factor that depends only on d and some error parameter ε∈ (0,1). The factor of approximation is (1+ε)4^d-1,4^d-1 and 4^d-1, respectively for (<ref>), (<ref>) and (<ref>). The running time is O(nlogn), n^O(d) and n^O(d), respectively for (<ref>), (<ref>) and (<ref>); see Table <ref>.










  
Coresets. Our main algorithm takes time n^O(d) and is easily generalized for many objective functions via Observation <ref>, Theorem <ref> and Theorem <ref>. It implies the results in the last three rows of Table <ref>.
For the case that p≥ 1 and we wish to minimize Ax-b_p over every unit vector x, we can apply our algorithm on the coreset for ℓ_p regression as explained in the previous section. This reduces the running time to near-linear in n, and enables parallel computation over distribution and streaming data by applying it on the small coreset that is maintained on the main server. This explains the running time and approximation factor for the first three lines of Table <ref>.



  
Our experimental results show that the suggested algorithms perform better, in both accuracy and computation time, compared to the few state of the art methods that can handle this non-convex problem; See Section <ref>.

Table <ref> summarizes the main contributions of this paper.






§ PRELIMINARIES
 
In this section we first give notation and main definitions that are required for the rest of the paper.



  
Notation. 
Let ^n× d be the set of n× d real matrices. In this paper, every vector is a column vector, unless stated otherwise, that is ^d = ^d× 1. We denote by p=p_2=√(p_1^2+…+p_d^2) the length of a point p=(p_1,⋯,p_d) ∈^d, by (p,π) = min_x ∈πp-x_2 the Euclidean distance from p to a subspace π of ^d.
We denote [n]=1,⋯,n for every integer n ≥ 1.
For a bijection function :[n]→ [n] (matching function) and a set Y = (a_1,b_1),⋯,(a_n,b_n)⊆^d× of n pairs, we define Y_ = (a_1,b_(1)),⋯,(a_n,b_(n)).
For every i∈ [d], we denote by e_i ∈^d the ith standard vector, i.e., the ith column of the d× d  identity matrix. For d≥ 2, ^d-1 = x ∈^d |x=1 is the set of all unit vectors in ^d. We denote = ^1 for simplicity.
For every x∈, we define (x) =  1    if  x≥ 0 
 -1    otherwise.
The set (n) denotes the union over every matching function :[n]→[n].
We remind the reader that _x∈ X f(x) is the set (and not a scalar) that contains all the values of x that minimize f(x) over some set X.
We denote by 0⃗(d) = (0,⋯,o)^T ∈^d.



  
Definitions. We first give a brief geometric illustration of the later definitions.
A hyperplane h in ^d that has distance b̂>0 from the origin can be defined by its orthogonal (normal) unit vector â so that h=x∈^d|â^Tx-b̂=0|. More generally, the distance from x∈^d to h is |â^Tx-b̂|. In Definition <ref> below the input is a set h_1,⋯,h_m of m such hyperplanes that are defined by the matrix A = (a_1 |·| a_m)^T ∈^m× d and vector b = (b_1,⋯,b_m)^T∈^m, such that a_i/a_iis the unit normal to the ith hyperplane, and b_i/a_i is its distance from the origin.

In what follows we define a (possibly infinite) set of unit vectors (A,b), which are the unit vectors that are as close as possible to the hyperplane h_m, among all unit vectors that lie on the intersection of h_1,⋯,h_m-1. If m=1, we define this set to be the set of closest vectors on the unit sphere to the hyperplane h_1.
Observe that for every point x∈^d, x ∈ h_i  a_i^Tx=b_i; See Figure <ref> and <ref> for a geometric illustration of the following definition.
 
Let m,d≥ 1 be a pair of integers, A = (a_1 |⋯| a_m)^T ∈^m× d and b = (b_1,⋯,b_m) ∈ [0,∞)^m. We define the set

    (A,b):=
       _x∈^d-1 |a_m^Tx-b_m|,     m=1 
    _x∈^d-1: 
    (a_1 |⋯| a_m-1)^Tx=(b_1,⋯,b_m-1)^T |a_m^Tx-b_m|,    .










In what follows, for every pair of vectors v=(v_1,⋯,v_n) and u=(u_1,⋯,u_n) in ^n we denote v≤ u if v_i ≤ u_i for every i∈ [n]. The function f:^n → [0,∞) is non-decreasing if f(v)≤ f(u) for every v≤ u. For a set I in ^d, and a scalar c∈ we denote I/c=x/c| x∈ I.

The following definition is a generalization of Definition 2.1 in <cit.> from n=1 to n>1 dimensions, and from  to I⊆^n. Intuitively, an r-log-Lipschitz function is a function whose derivative may be large but cannot increase too rapidly (in a rate that depends on r).
 
Let r>0 and let n≥ 1 be an integer. Let I be a subset of ^n, and h:I→[0,∞) be a non-decreasing function.
Then h(x) is r-log-Lipschitz over x∈ I, if for every c ≥ 1 and x∈ I ∩I/c, we have h(c x)≤ c^r h(x). The parameter r is called the log-Lipschitz constant.


The following definition implies that we can partition a function g which is not a log-Lipschitz function into a constant number of log-Lipschitz functions; see Figure <ref> for an illustrative example.

Let g:X→[0,∞) be a continuous function over a set X, and
let (X,dist) be a metric space, i.e. :X^2→[0,∞) is a distance function. Let r> 0. The function g is piecewise r-log-Lipschitz if there is a partition of X into m disjoint subsets X_1,⋯, X_m
such that for every i∈[m]:





  * g has a unique infimum x_i at X_i, i.e., x_i= _x∈ X_ig(x).

  * h_i:[0,max_x∈ X_i(x,x_i)]→[0,∞) is an r-log-Lipschitz function; see Definition <ref>.

  * g(x)=h_i((x_i,x)) for every x∈ X_i.



The set of minima is denoted by M(g)=x_1,⋯,x_m.








§ REGRESSION WITH A GIVEN MATCHING
 
In this section we suggest our main approximation algorithm for solving (<ref>), i.e., when the matching between the rows of A and the entries b = (b_1,⋯,b_n)^T are given.

The following two corollaries lies in the heart of our main result.
 
Let b ≥ 0 and let g:[0,2π) → [0,∞) such that g(α) = |sin(α) - b|. Then g is a piecewise 2 log-Lipschitz function; See Definition <ref>.



 
Let a ∈^2∖0 and b ≥ 0. Let y ∈_x∈^1 |a^T x - b|. Then for every u^*,u' ∈ such that u'-y≤ 2·u^*-y we have

    |a^Tu'-b| ≤ 4· |a^Tu^*-b|.



See proof of Corollary <ref> in the appendix.


In what follows, a non-zero matrix is a matrix of rank at least one (i.e., not all its entries are zero).
Let A = (a_1 |⋯| a_n)^T ∈^n× d and b = (b_1,⋯,b_n) ∈ [0,∞)^n, where we assume that A is a non-zero matrix.
For every i∈ [n], let h_i be a hyperplane whose normal is a_i/a_i_i=1^n and its distances from the origin is b_i/a_i_i=1^n, and let H = h_i_i=1^n denote their union.
Let x^*∈^d be a unit vector. The following lemma generalizes Corollary <ref> to higher dimensions. It states that there is a unit vector x' of minimal distance (h_k,x') to one of the n hyperplanes h_k, that approximates a_i·(h_i,x^*) = |a_i^Tx^*-b_i| for every i∈ [n] up to a multiplicative factor of 4.
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d∖0^d× d be a non-zero matrix such that n ≥ d-1 ≥ 1 points, and let b = (b_1,⋯,b_n)^T ∈ [0,∞)^n and x^* ∈^d-1. Then there exists j∈ [n] where a_j≠ 0 and x_j ∈_x∈^d-1 |a_j^Tx-b_j| such that for every i∈ [n]

    |a_i^Tx_j-b_i| ≤ 4· |a_i^Tx^*-b_i|.



See proof of Lemma <ref> in the appendix.


The following lemma states that there is a set h_1,⋯,h_m⊆ H of m≤ d-1 hyperplanes from H, and a unit vector x' in the intersection of the first m-1 hyperplanes that is closest to the last hyperplane, x'∈_x∈⋂_i∈[m-1]h_i |a_m^Tx-b_m|, that is closer to every one of the n hyperplane in H, up to a factor of 4^d-1 than x^*, i.e., |a_i^Tx'-b_i| ≤ 4^d-1·|a_i^Tx^*-b_i| for every i∈ [n].
The proof is based on applying Lemma <ref> recursively m times. In the following lemma we use opt from Definition <ref>.
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a non-zero matrix of n ≥ d-1 ≥ 1 rows, let b = (b_1,⋯,b_n)^T ∈ [0,∞)^n, and let x^* ∈^d-1.
Then there is a set i_1,⋯,i_r⊆ [n] of r ∈ [d-1] indices such that for X = ((a_i_1|⋯| a_i_r), (b_i_1,⋯,b_i_r)) and every i∈ [n], there is x'∈ X that satisfies

    |a_i^Tx'-b_i| ≤ 4^d-1· |a_i^Tx^*-b_i|.

Moreover, (<ref>) holds for every x'∈ X if |X| = ∞.


See proof of Lemma <ref> in the appendix.





 §.§ Computing the set 
 
In this section we give a suggested simple implementation for computing the set (A,b) given a matrix A∈^m× d and b∈ [0,∞)^m; see Algorithm <ref>.
A call to (A,b) returns (A,b) if |(A,b)| ∈ O(1). Otherwise, it returns some x∈(A,b). In other words, the algorithm always returns a set of finite size.






 §.§ Geometric interpretation and intuition behind Algorithm <ref>. 
Algorithm <ref> takes a set h_1,⋯,h_m of m hyperplanes in ^d as input, each represented by its normal and
its distance from the origin, and computes a point x∈^d-1∩⋂_i∈[m-1] h_i that minimizes its distance to h_m, i.e., (x,h_m). In other words, among all vectors
that lie simultaneously on h_1,⋯,h_m-1, we intend to find the unit vector x which minimizes its distance to h_m. There can either be 0,1, 2 or infinite such points.
In an informal high-level overview, the algorithm basically starts from some unit vector x∈^d, and rotates it until it either intersects h_1, or minimizes its distance to h_1 without intersecting it.
If they intersect, then we rotate x while maintaining that x∈ h_1, until x either intersects h_2 or minimizes its distance to h_2. If x intersects h_2, we rotate x in the intersection h_1∩ h_2 until it intersects h_3 and so on.
If at some iteration we observe that all the remaining subspaces are parallel, then the set (A,b) is of infinite size, so we return one element from it.
We stop this process when at some step k, we can not intersect h_k under the constraint that x is a unit vector and x∈⋂_i∈[k-1]h_i. If k<m, then (A,b) is empty. If k=m we return the vector x that is as close as possible to h_m.

More formally, at each step of the algorithm we do the following. In Line <ref> we check weather the intersection h_1∩^d-1 contains at most 1 point. This happens when the distance b_1/a_1 of the
hyperplane h_1 to the origin is bigger than or equal to 1. This is the simple case in which we terminate. The interesting case is when the distance of the hyperplane to the origin is less than 1.
In this case, if d=2, we compute and return the two possible intersection points in h_1∩^d-1. If d≥ 3, then |h_1∩^d-1|=∞. Observe that h_1∩^d-1 is simply
a sphere ' ⊆ h_1 of dimension d-2, but is not a unit sphere. We rotate the coordinates system such that the hyperplane containing ' is orthogonal to the dth standard vector e_d. We observe that
in order to minimize the distance from a point in ' to h_2, we need to minimize its distance to h_1∩ h_2 which is simply a (d-2)-subspace contained in h_1.
We thus project ' and every (d-2)-subspace in h_1∩ h_i | i∈ [m] ∖1 onto the hyperplane H orthogonal to e_d and passes through the origin, obtaining a sphere ” of dimension d-2, and m-1
(d-2)-subspaces, all contained in H. We scale the system such that ” becomes a unit sphere. We then continue recursively.




  
Overview of Algorithm <ref>. The input is a matrix A and a vector b. The output is a set of unit vectors that satisfies Theorem <ref>. We assume without loss of generality that the entries of b are non-negative, otherwise we change the corresponding signs of rows in A in Lines <ref>-<ref>.
In Lines <ref>-<ref> we run exhaustive search over every possible set of at most d-1 indices, which corresponds to a set S of unit vectors, to find the set X from Lemma <ref>.
The algorithm then returns the union of these sets in Line <ref>.
See Algorithm <ref> for a possible implementation for Line <ref>. Notice that Algorithm <ref> does not return a set of infinite size. In such a case, it returns one element from the infinite set.






What follows is the main theorem of Algorithm <ref>.
The following theorem states that the output X of Algorithm <ref> contains the desired solution.
This is by Lemma <ref> that ensures that one of the sets S contains the desired solution. 
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix of n ≥ d-1 ≥ 1 rows, and let b = (b_1,⋯,b_n)^T ∈^n.
Let X ⊆^d-1 be an output of a call to (A,b); see Algorithm <ref>. Then for every x^* ∈^d-1 there exists a unit vector x' ∈ X such that for every i∈[n],

    |a_i^Tx'-b_i| ≤ 4^d-1· |a_i^Tx^*-b_i|.

Moreover, the set X can be computed in n^O(d) time and its size is |X| ∈ n^O(d).


See proof of Theorem <ref> in the appendix.





 §.§ Generalization
 
We now prove that the output of Algorithm <ref> contains approximations for a large family of optimization functions. Note that each function may be optimized by a different candidate in X.
This family of functions includes squared distances, M-estimators and handling k≥ 1 outliers. It is defined precisely via the following cost function, where the approximation error depends on the parameters r and s.

Let Y = y_1,⋯,y_n be a finite input set of elements and let Q be a set of queries.
Let D:X× Q→ [0,∞) be a function.
Let :[0,∞)→ [0,∞) be an r-log-Lipschitz function and f:[0,∞)^n→ [0,∞) be an s-log-Lipschitz function.
Let q∈ Q. We define

    (Y,q)=f( (D(y_1,q)),⋯, (D(y_n,q))).


In what follows, Y = (a_1,b_1),⋯,(a_n,b_n)⊆^d ×, Q = ^d-1 will be the set of all unit vectors in ^d, and D((a,b),x) = |a^Tx-b|.
Table <ref> gives some examples of different cost functions that satisfy the requirements of Definition <ref>.




The following observation states that if we find a query q∈ Q that approximates the function D for every input element, then it also approximates the function  as defined in Definition <ref>.
 [Observation 5 in <cit.>] 
Let (Y,q)=f( (D(y_1,q)),⋯, (D(y_n,q))) be defined as in Definition <ref>.
Let q^*,q' ∈ Q and let c≥ 1. If D(y_i,q')≤ c· D(y_i,q^*) for every i∈ [n], then

    (Y,q')≤ c^rs·(Y,q^*).




The optimal solution in the following theorem can be computed by taking the optimal solution for the corresponding cost function at hand among the output set X of candidates from Algorithm <ref>.
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix of n ≥ d-1 ≥ 1 rows, and let b = (b_1,⋯,b_n)^T ∈^n. Let , s, r be as defined in Definition <ref> for Y = (a_i,b_i)| i∈ [n] and D((a,b̂),x) = |a^Tx-b̂| for every a∈^d, b̂∈ and x∈^d. Then in n^O(d) time we can compute a unit vector x'∈^d-1 such that

    (Y,x') ≤ 4^(d-1)rs·min_x∈^d-1(Y,x).



See proof of Theorem <ref> in the appendix.





§ HANDLING UNKNOWN MATCHING
 
In this section, we tackle the constrained ℓ_p regression problem with unknown matching between the rows of A and entries of b. That is, where the minimum of (<ref>) is not only over every nit vector x, but also over every permutation of the entries of b (or rows in A). Thus, the problem now at hand is to solve (<ref>).
The main result is derived from the following simple observation. In Lemma <ref> we proved there is a set of r≤ d-1 matching subsets Â and b̂ of A and b respectively, such that (Â,b̂) contains the desired unit vector x'. This result still holds when the matching is unknown, though we dont know whose the subset b̂ in b that corresponds to Â. Hence, for every such Â, all we have to do is apply another exhaustive search over subset of b in order to find the subset b̂ that best matches Â.

Algorithm <ref> uses the following definition of optimal matching between an input set of pairs Y, a query q, and a cost function.
 
Let (n) denote the union over every matching function (bijection function) :[n]→[n] and let Q be a set. For every i∈ [n], let y_i = (a_i,b_i) be a pair of elements, and let Y = y_1,⋯,y_n be their union.
Consider a function  as defined in Definition <ref> for f(v) = v_1 and let q ∈ Q.
We denote by (Y,q,) the matching function that minimizes (Y_,q) over every permutation ∈(n). Formally,

    (Y,q,) ∈_∈(n)(Y_,q).







  
Overview of Algorithm <ref>. In Lines <ref>-<ref>, we iterate over every pair of sets I=i_1,⋯,i_r and ℓ=ℓ_1,⋯,ℓ_r of r indices, for r≤ d-1. We assume that a_i_k matches b_ℓ_k for every k∈ [r]. We change the signs of a_i_k and b_ℓ_k as in Algorithm <ref>. We make a call to Algorithm <ref> to compute and store in S the set (a_i_1'|⋯| a_i_r'),(b_ℓ_1',⋯,b_ℓ_r'). Notice that we do not use infinite sets in Algorithm <ref>. We then add S to X;
In Lines <ref>-<ref> we compute for every x∈ X the matching (Y,x,) that minimizes the given cost function , and store the pair in S. In Line <ref> we return the pair in S whose cost is minimal.




The following theorem states that Algorithm <ref> indeed returns an approximation to the optimal cost, even when the matching between rows of A and b is not given.
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix containing n ≥ d-1 points in its rows, let b = (b_1,⋯,b_n)^T ∈^n, and let Y = (a_i,b_i)| i∈ [n]. Let D((a,b̂),x) = |a^Tx-b̂| for every a∈^d, b̂∈ and x∈^d. Let r be a scalar and  be a function that satisfy Definition <ref> for D and f(v) = v_1. Let (x̃,) be a pair of unit vector and permutation (matching function) which is the output of a call to (A,b,); see Algorithm <ref>. Then it holds that

    (Y_,x̃) ≤ 4^(d-1)r·min_x,(Y_,x),

where the minimum is over every unit vector x∈^d-1 and ∈(n). Moreover, (x̃,) is computed in n^O(d) time.


See proof of Theorem <ref> in the appendix.





§ CORESET FOR LINEAR REGRESSION
 
In this section we provide a coreset for the constrained ℓ_p regression, which is a small improvement and generalization of previous results <cit.>



 §.§ Improvements via coreset

An ε-coreset, which is a compression scheme for the data, is suggested in Theorem <ref>. Streaming and distribution in near-logarithmic update time and space, including support for deletion of points in near-logarithmic time in n (using linear space) is also supported by our algorithm when the matching between the rows of the matrix A and the entries of b is given; see Section <ref>. This is due to the fact that our coresets are composable, i.e., can be merged and reduced over time and computed independently on different subsets. This is now a standard technique, we refer the reader to  <cit.> for further details.



 
Let d ≥ 2 be a constant integer. Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix containing n ≥ d-1 points in its rows, let b = (b_1,⋯,b_n)^T ∈^n, and let w = (w_1,⋯,w_n) ∈ [0,∞)^n.
Let ε, δ∈ (0,1) and let z∈ [1,∞).
Then in O(n logn) time we can compute a weights vector u = (u_1,⋯,u_n)∈[0,∞)^n that satisfies the following pair of properties.



  * With probability at least 1-δ, for every x∈^d-1 it holds that

    (1-ε) ·∑_i∈ [n] w_i ·  |a_i^Tx-b_i|^z  ≤∑_i∈ [n]u_i· |a_i^Tx-b_i|^z
    ≤   (1+ε) ·∑_i∈ [n] w_i · |a_i^Tx-b_i|^z.


  * The weights vector u has O(log1/δ/ε^2) non-zero entries.



See proof of Theorem <ref> in the appendix.





§ EXPERIMENTAL RESULTS
 
To demonstrate the correctness and robustness of our algorithms, we implemented them in Matlab and compared them to some commercial software for non-convex optimization. A discussion is provided after each test, and an overall discussion is provided in Section <ref>. Open code is provided <cit.>.

Hardware. All the following tests were conducted using Matlab R2019a and Maple 2018 on a Lenovo W541 laptop with an Intel i7-4710MQ CPU @ 2.50GHZ and 8GB RAM.


Let N(μ,σ) denote a Gaussian distribution with mean μ and standard deviation σ and let U(r) denote a uniform distribution over [0,r].















 §.§ Constrained ℓ_p regression
 

The dataset. We conducted 3 tests with different values of p and d. The data for each test was a matrix A∈^n× d of n points, and a vector b∈^n, where each entry of A and b was sampled from U(200), i.e. uniformly at random from [0,200].

Objective. The test aims to solve problem (<ref>).
The goal was to minimize the constrained ℓ_p regression function f_p = Ax-b_p for different values of p≥ 1.

Algorithms. We compare the following 5 algorithms, which aim to minimize f_p over every unit vector x.

: This algorithm computes X = (A,b), and returns min_x∈ X f_p.

: Maple <cit.> provides a function named GlobalOptima from the DirectSearch optimization package version 2 <cit.>. This function takes as input an objective function, a constraint, and an integer which specifies the number of initial ”simulated points“, where both the accuracy and the computation time increase with larger number of simulated points. The function aims to compute the global minima of non-linear multivariate functions under given constraints. We call this function from Matlab with the objective function f_p, the constraint x=1, and the default number (30) of initial simulated points.

: same as  but with number=100 of initial simulated points.

: The Yalmip library <cit.> for Matlab provides a function named optimize, which takes as input an objective function to minimize, a constraint, and a solver program, and aims to optimize the given function under the given constraint using the given solver.
We run this function with the objective function f_p, the constraint x=1, and the bmibnb solver <cit.> which aims to solve non-convex problems.

: same as the previous algorithm, but with the Baron solver <cit.> which is a widely used commercial software for global optimization.

The results are shown in Figures <ref> and <ref>.

Discussion. Figure <ref> presents the computation time required for the suggested methods with p≥ 1. We conducted two tests: p=1 and p=3.5, both for increasing values of n.

As the graph for p=1 shows,  managed to compute the output in a fraction of the time it took other methods to compute their output, except for , which was a close second. In this test, the cost (error) of all methods was similar.

The graph for p=3.5 presents similar results. However, this time  took time which is more than 10 folds larger than the other methods. Hence it is presented as a thick black line at the top of the graph. Also for this case the cost (error) of all methods was similar.

Figure <ref> presents both the cost (error) and the computation time of the suggested methods for p=0.1 with increasing values of n. In this test, we were only able to run  and , as the other methods either did not support p<1 or did not return an output in a reasonable amount of time.

As the time comparison graph shows,  required time which is 2 orders of magnitude smaller than the time it took  and .

The costs graph shows that at n=10, the cost of  and  start to grow exponentially faster than the cost of . The gap between the cost of and the cost of the other methods continues to grow with n.





 §.§ Robustness to outliers
 
The test in this subsection aims to solve Problem (<ref>).


The dataset. We generated a matrix A ∈^80× 3 where each entry was sampled randomly from U(200), and a unit vector x∈^3 was randomly generated. We then define b = Ax . Gaussian noise drawn from N(0,30) was then added to each of the entries of A and b respectively.
We then added noise sampled randomly from U(20000) to k· A entries of A and their corresponding entries of b, where k∈ [0,1/2].



  
Objective. This test aims to solve (<ref>).
The goal of the experiment was to minimize the f_O(A,b,x) = ∑_i∈ [n](Ax-b,3n/4)_1 over every unit vector x.






  
Algorithms. We compared the following 4 algorithms, which aim to minimize f_O(A,b,x) over every unit vector x.

 above with objective function f_O(A,b,x).

: A random sample consensus (RANSAC) scheme <cit.>, where at each iteration we sample uniformly at random corresponding subsets Â, and b̂ of size 3 from A and b respectively, called  with objective function Âx-b̂_1 and constraint x=1 to obtain some unit vector x̂∈^3, computed the number of pairs (a_i,b_i) that satisfy |a_i^Tx̂-b_i| ≤ 200 (inlier points) and pick the unit vector with maximal number of inliers denoted by A^* and b^*. Final call is made to  with objective function A^*x-b^*_1 and constrain x=1 that returns its output.

: As  but calls  instead of .

: As  but calls  instead of .

The results are shown in Figure <ref>.

Discussion. Figure <ref> presents both the cost (error) and the computation time of the suggested methods, while increasing the value of k (the fraction of outliers in the data).

As the time comparison graph shows,  managed to compute the output faster than the other methods.

The costs graph demonstrates the robustness in practice of our algorithm to outliers. The cost of  was roughly constant as the fraction k of outliers grew until k=0.25. At k>0.25, the cost of  had a sharp increase. This is not surprising since the cost function sums over 75% of the points with the smallest cost. In other words, when k>0.25, 75% of the data with the smallest cost must contain points who are outliers.
The other methods showed an immediate increase in cost as soon as outliers were introduced to the data (k>0).



 §.§ Unknown correspondences between A and b
 
The test in this subsection aims to solve (<ref>).

The dataset is a matrix A ∈^n× 3 where each entry was sampled randomly from U(200), and Ax=b for some random unit vector x∈^d. The rows of A have then been shuffled with a random permutation. We then added small Gaussian noise drawn from N(0,10) to every entry of A and b respectively.



Objective. Let Y=(a_i,b_i)_i=1^n. The goal of the experiment was to minimize f_M(Y_,x) = ∑_i∈ [n] |a_i^Tx-b_(i)| over every unit vector x∈^3 and permutation .



  
Algorithms. We compared the 5 following algorithms.

: Run (A,b,f_M); See Algorithm <ref>.

: An iterative closest point (ICP) scheme where we alternate between the following two steps until there is no sufficient change in the cost function: (1) Call  with objective function A'x-b_1 and constraint x=1 to obtain a unit vector x', where A'=A at the first iteration, and (2) Compute the optimal permutation (a shuffling of A') that minimizes the differences between A'x' and b, and apply the shuffling to A'. This is done using the Hungarian method <cit.>.

: Similar to  but calls  instead of .

: Similar to  but calls  instead of .

: Similar to  but calls  instead of .

The results are shown in Figure <ref>.

Discussion. As the costs graph in Figure <ref> shows,  solved (<ref>) with higher accuracy than all other methods.




  §.§.§ Unknown correspondences between A and b, without noise. 

We conducted another test similar to the previous test, for the case where the correspondences are unknown, but without adding noise to the data. In other words, a matrix A ∈^n× 3 was generated, where each entry was sampled randomly from U(200), and a corresponding vector b∈^n was computed such that Ax=b for some random unit vector x∈^d. The rows of A have then been shuffled with a random permutation. We did not add noise to the data. The objective function and the algorithms tested are similar to the previous test.

The results are shown in Figure <ref>.

Discussion. As expected due to the constant factor approximation guaranteed of all our algorithms, the cost of  was constant at 0. Which means was able to recover the unit vector and the matching function correctly, while the other methods could not.






 §.§ Overall discussion
 


  
Variance. As all the figures in Section <ref> show, the variance of our algorithms was smaller and more stable than the variance of the other methods. This happens since our algorithms are guaranteed to compute an approximated ("good") result. In other words, the output cost of our algorithms can not have a dignificant change, since it is guaranteed to be with in small range of the optimal objective function value. Hence the small variance. The running time of our algorithms is deterministic. Hence the computation time variance is small as well.



  
The constants behind our approximation factors. 
Our theorems guarantee that the output cost of our algorithms is always within range of the optimal value of the objective function used, up to a constant factor of 4^d-1∈ O(4^d). This bound on the approximation constant is a worst case analysis bound. In practice however, this constant is much smaller.



  
Decreasing the computation time of our algorithms.
Observe that since our algorithms are embarrassingly parallel, using a computer with M cores would have reduced the running time of our algorithms by a factor of M.




§ CONCLUSION AND OPEN PROBLEMS
 
We proposed the first provable polynomial time approximation algorithms for the constrained ℓ_p regression problem, for any constant p>0 and d≥1, including versions for handling outliers, and unknown order of rows in A. Using coresets, the running time is near linear in some cases. Experimental results show that our algorithms outperform existing commercial solvers. Open problems: (i) Running time that is polynomial in d is hopeless since the problem is NP-hard, but additive approximations may be obtained via projection on random subspaces or PCA, (ii) (1+)-approximations, (iii) near-linear time algorithms for the unknown matching case, and streaming version for these cases.



icml2019















§ REGRESSION WITH GIVEN MATCHING


The following corollary states that if we double the distance on a unit sphere to a point y, which is the closest point on the unit sphere to a line ℓ, then the distance to ℓ will grow by at most a multiplicative factor of 4.
 
Let a ∈^2∖0 and b ≥ 0. Let y ∈_x∈^1 |a^T x - b|. Then for every u^*,u' ∈ such that u'-y≤ 2·u^*-y we have

    |a^Tu'-b| ≤ 4· |a^Tu^*-b|.



Let g:[0,2π) → [0,∞) such that g(α) = |sin(α) - b/a|.
Let M(g) = _α∈ [0,2π)g(α) be the set of minima of g.

By replacing b with b/a in Corollary <ref>, g is a piecewise 2-log-Lipschitz function. Hence, by Definition <ref>, there is a partition X_1,⋯,X_m of [0,2π), a set of 2-log-Lipschitz functions h_1,⋯,h_m, and a set M(g) = α^*_1,⋯,α^*_m of minima such that for every i∈ [m] and α_1,α_2 ∈ X_i where |α_2-x_i| ≤ 2· |α_1-_i|
we have that

    g(α_2) = h_i(|α_2-x_i|) ≤ 4· h_i(|α_1-x_i| =4· g(α_1),

where the first derivation holds by property (iii) of Definition <ref>, the second derivation holds by combining that h_i is a 2-lop-Lipschitz function with |α_2-x_i| ≤ 2· |α_1-x_i|, and the last derivation holds by Property (iii) of Definition <ref>.





Without loss of generality, assume that a/a = (1,0)^T. Otherwise rotate the coordinates system.
Let x:→^2 such that x(α) = (sinα,-cosα)^T.
For every α∈ [0,2π) we now have that

    |a^Tx(α) - b|    = a|a^T/ax(α)-b/a|
       = a|sinα-b/a| = a· g(α).


Let α_y ∈ [0,2π) such that y = x(α_y). We now have that

    a g(α_y)     = |a^Ty-b| =min_α∈[0,2π)|a^Tx(α)-b|
       = min_α∈[0,2π)a g(α),

where the first derivation is by combining (<ref>) and the definitions of y and α_y, the second derivation is by the definition of x(α), and the last derivation is by (<ref>).
Hence, α_y ∈ M(g).

Let u',u^* ∈ and α',α^* ∈ [0,2π) such that u' = x(α') and u^*=x(α^*). Since u'-y≤ 2·u^*-y by the assumption of the Corollary, we get that |α'-α_y| ≤ 2 |α^*-α_y| by simple linear algebra.
Due to the last inequality and since α_y ∈ M(g), we can substitute α_2 = α' and α_1 = α^* in (<ref>) to obtain

    g(α') ≤ 4 g(α^*).

Corollary <ref> now holds as

    |a^Tu'-b|     = |a^Tx(α')-b|=a· g(α')
       ≤ 4a· g(α^*) = 4|a^Tx(α^*)-b|
       =4|a^Tu^*-b|,

where the first derivation is by the definition of u', the second derivation is by substituting α=α' in (<ref>), the third derivation is by (<ref>), the fourth derivation is by substituting α=α^* in (<ref>) and the last derivation is by the definition of α^*.




Let A = (a_1 |⋯| a_n)^T ∈^n× d∖0^d× d be a non-zero matrix such that n ≥ d-1 ≥ 1 points, and let b = (b_1,⋯,b_n)^T ∈ [0,∞)^n and x^* ∈^d-1. Then there exists j∈ [n] where a_j≠ 0 and x_j ∈_x∈^d-1 |a_j^Tx-b_j| such that for every i∈ [n]

    |a_i^Tx_j-b_i| ≤ 4· |a_i^Tx^*-b_i|.



For every i∈ [n], if a_i = 0, then for every x ∈^d-1 it holds that |a_i^Tx-b_i| = |b_i| ≤ 4 |b_i| = 4 |a_i^Tx^*-b_i|.
Therefore, we ignore vectors of zero length and assume that A does not contain any zero row.

For every i∈ [n], let

    M_i = _x∈^d-1|a_i^T/a_i x -b_i/a_i|,


    x' ∈_x ∈⋃_i∈ [n] M_ix^* - x,


    j ∈ [n]  such that x' ∈ M_j,

i.e., M_i is the set of all unit vectors that minimize |a_i^Tx-b_i|, x' is the unit vector that is closest to x^* among all vectors in ⋃_i∈ [n] M_i, and j is the index of the set M_j such that x' ∈ M_j.

Put i∈ [n]. We prove that

    |a_i^Tx'-b_i| ≤ 4· |a_i^Tx^*-b_i|.

This would prove the lemma since x' ∈ M_j = _x∈^d-1 |a_j^Tx-b_j|.

Indeed, assume without loss of generality that a_i/a_i = (1,0,⋯,0)^T and that x^* = (x^*_1,x^*_2,0,⋯,0)^T, x^*_2 ≥ 0. Otherwise, rotate the coordinate system as follows: rotate a_i/a_i until it coincides with the x-axis. Then, rotate the system around the x-axis (i.e., without changing a_i/a_i), until x^* intersects the xy-plane at the halfspace with positive y values.

We prove (<ref>) via the following case analysis: (i) b_i ≤a_i and (ii) b_i>a_i.

Case (i): b_i ≤a_i. Let y = (y_1,y_2,⋯,y_d)^T := (b_i/a_i,√(1-b_i^2/a_i^2),0,⋯,0)^T. Observe that y ∈ M_i since

    |a_i^Ty-b_i|     = a_i|a_i^T/a_iy - b_i/a_i|
        = a_i|b_i/a_i - b_i/a_i| = 0.


Identify x' = (x'_d1,⋯,x'_d) and x^* = (x^*_1,⋯,x^*_d).
Let v' = (x'_1, √(1-x'_1^2))^T, v^* = (x^*_1,√(1-x^*_1^2))^T and y' = (y_1,y_2)^T = (b_i/a_i,√(1-b_i^2/a_i^2))^T.
Since v'_1 = x'_1 and v' = x'=1 we obtain

    v'_2^2     =v'^2-v_1^2=v'^2-x'_1^2
       =x'^2-x'_1^2 =∑_j=2^d x'_j^2 ≥x'_2^2

We now have that

    x'-x^*^2     = x'^2 + x^*^2 -2x'^Tx^* 
        = v'^2 + v^*^2 -2(x'_1· x^*_1 + x'_2· x^*_2) 
        = v'^2 + v^*^2 -2(v'_1· v^*_1 + x'_2· x^*_2) 
       ≥v'^2 + v^*^2 -2(v'_1· v^*_1 + v'_2· x^*_2) 
        = v'-v^*^2 ,

where (<ref>) holds since only the first two entries of x^* are non-zero, (<ref>) is by the definitions of v' and v^*, and (<ref>) holds by combining the squared root of (<ref>) with the assumption that x^*_2 ≥ 0.

Hence, we obtain that

    v'-v^*   ≤x'-x^*
       ≤y-x^*
        = y'-v^*,

where (<ref>) is by taking the squared root of (<ref>), (<ref>) holds by combining the definition of x' and the fact that y ∈ M_i ⊆⋃_l∈ [n] M_l, and (<ref>) holds by the definitions of y' and v^*.
Hence,

    v'-y'≤v'-v^* + v^*-y'≤ 2·v^*-y',

where the first inequality is the triangle inequality, and the second inequality is by (<ref>).

Since |(1,0)y'-b_i/a_i| = 0 we have y' ∈_x∈ |(1,0)x - b_i/a_i|.
Therefore, substituting a = (1,0)^T, b=b_i/a_i, y=y', u^* = v^* and u' = v' in Corollary <ref> yields that

    |(1,0)v'-b_i/a_i| ≤ 4·|(1,0)v^*-b_i/a_i|.


Hence,

    |a_i^Tx'-b_i|     = a_i·|a_i^T/a_ix'-b_i/a_i|
        = a_i·|(1,0)v'-b_i/a_i|
       ≤ 4a_i·|(1,0)v^* -b_i/a_i|
    = 4· |a_i^Tx^* - b_i|,

where in the first equality we multiply and divide by a_i, the second and last equality are by the assumption a_i/a_i=(1,…,0)^T, and the inequality is by (<ref>).

Case (ii): b_i>a_i. Let a = a_i/a_i and b = b_i/a_i.
In this case, we have

    |a_i^Tx'-b_i|     = a_i·|b-a_i^T/a_ix'|
       = a_i·|b-1+1-a_i^T/a_ix'|
       = a_i·(b-1 + |1-a^Tx'|),

where the last derivation holds since b-1>0 by the assumption of Case (ii), and 1≥ x'_1 = a_i^T/a_ix' = a^Tx' since x' is a unit vector.

Since for b_i=1 and a_i = a the condition b_i ≤a_i of Case (i) holds, we obtain by Case (i) that

    |1-a^Tx'| = |a^Tx'-b_i| ≤ 4|a^Tx^*-b_i| = 4 |a^Tx^*-1|.


This proves Case (ii) as

    |a_i^Tx'-b_i|     = a_i·(b-1 + |1-a^Tx'|) 
       ≤a_i·(b-1 + 4·|1-a^Tx^*|) 
       ≤ 4a_i·(b-1 + |1-a^Tx^*|) 
        = 4a_i·|b-a^Tx^*| 
        = 4· |a_i^Tx^*-b_i| ,

where (<ref>) is by (<ref>), (<ref>) is by (<ref>), (<ref>) holds since b-1 ≥ 0, and (<ref>) holds since a^Tx^* ≤ 1 ≤ b.


As explained in Section <ref>, the optimization problem min_x∈^d-1Ax-b_1 can be interpreted as computing a point x' on the unit sphere that minimizes the sum of distances to n given hyperplanes H = h_1,⋯,h_n in ^d.
The following lemma suggests that there are m≤ d-1 hyperplanes from H and a unit vector x' in the intersection in the intersection of the first m-1 hyperplanes, that is closest to the last hyperplane, that is closer to every one of the n hyperplane in H, up to a factor of 4^d-1 than x^*, i.e., |a_i^Tx'-b_i| ≤ 4^d-1·|a_i^Tx^*-b_i| for every i∈ [n].

To find such a vector x', we begin with an initial set of candidate solutions ^d-1, i.e., all unit vectors in ^d.
The proof of the following lemma consists of at most d-1 steps, each step adds c≥ 1 more constraints on the candidate unit vectors for the solutions, and adds another factor of 4 to the final approximation factor. This induction terminates when the number of added constraints is d-1 (so there is a finite set of candidates for x').
Each step suggests that we can rotate x', until it minimizes its distances |a_j^Tx'-b_j| to one of the hyperplanes (the jth hyperplane in the proof), without increasing each of the other n-1 distances by more than a multiplicative factor of 4. This follows from Lemma <ref>.
Afterwards, we reduce the problem to an instance of the same optimization problem but with less free parameters.
 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a non-zero matrix of n ≥ d-1 ≥ 1 rows, let b = (b_1,⋯,b_n)^T ∈ [0,∞)^n, and let x^* ∈^d-1.
Then there is a set i_1,⋯,i_r⊆ [n] of r ∈ [d-1] indices such that for X = ((a_i_1|⋯| a_i_r), (b_i_1,⋯,b_i_r)) and every i∈ [n], there is x'∈ X that satisfies

    |a_i^Tx'-b_i| ≤ 4^d-1· |a_i^Tx^*-b_i|.

Moreover, (<ref>) holds for every x'∈ X if |X| = ∞.


For every i∈ [n], if a_i = 0, then for every x ∈^d-1 we have |a_i^Tx-b_i| = |b_i| ≤ 4· |b_i| = 4· |a_i^Tx^*-b_i| and (<ref>) trivially holds.
Therefore, in the rest of the proof we assume a_i≠ 0 for every i∈ [n].


The proof is by induction on the dimension d.

Base case for d=2: Substituting A, b, and d=2 in Lemma <ref> yields that there is j∈ [n] and x_j∈_x∈|a_j^Tx-b_j| that satisfy, for every i∈ [n],

    |a_i^Tx_j-b_i| ≤ 4· |a_i^Tx^*-b_i|.

For every a∈^2 and b ≥ 0, we have that

    _x∈|a^Tx-b|
       = a/a,     b ≥a
    x∈|a^T/ax = b/a,    ;

see Figure <ref> for a geometric illustration of the solutions for _x∈|a^Tx-b|, in both cases where b ≥a and b < a.
Substituting m=1, a_1 = a_j and b_1=b_j in Definition <ref> implies x_j ∈_x∈ |a_j^Tx-b_j| = (x_j,b_j). By (<ref>), |(x_j,b_j)| = |_x∈ |a_j^Tx-b_j|| ≤ 2.

Hence, X = (x_j,b_j) and x' = x_j satisfies (<ref>) in Lemma <ref>.

Case d≥ 3: Inductively assume that Lemma <ref> holds for d'=d-1.
For every i∈ [n], let _i = (a_i,b_i); See Definition <ref>. Put i∈ [n].
If a_i≤ b_i, then

    _i     = _x∈^d-1|a_i^Tx-b_i| = _x∈^d-1(b_i - a_i^Tx)
       = _x∈^d-1 a_i^Tx = a_i/a_i.

If a_i > b_i, then min_x∈^d-1|a_i^Tx-b_i| = 0. In this case,

    _i     = _x∈^d-1|a_i^Tx-b_i| = x ∈^d-1| a_i^Tx=b_i
       = x ∈^d-1|a_i/a_i^Tx=b_i/a_i.


Combining (<ref>) and (<ref>) yields that

    _i = a_i/a_i,    a_i≤ b_i 
    x∈^d-1|a_i/a_i^Tx=b_i/a_i,    a_i > b_i.
             .

Therefore,

    |_i| = 
               1,    a_i≤ b_i 
    ∞,    a_i > b_i
             .


Substituting A,b and x^* in Lemma <ref> yields that there is an index j ∈ [n] and a corresponding unit vector x_j∈_x∈^d-1 |a_j^Tx-b_j| = _j such that for every i∈ [n]

    |a_i^Tx_j-b_i| ≤ 4· |a_i^Tx^*-b_i|.

We continue with the following case analysis: (a): a_j≤ b_j and (b): a_j > b_j.

Case (a): a_j≤ b_j. Similarly to (<ref>), we have that

    _j = a_j/a_j,    a_j≤ b_j 
    x∈^d-1|a_j/a_j^Tx=b_j/a_j,    a_j > b_j.
             .

Combining the previous definition of _j with the assumption of Case (a) yields

    |_j| = |a_j/a_j| = 1.

By combining (<ref>), (<ref>) and the fact that a_j/a_j∈_j = (a_j,b_j), we obtain that Lemma <ref> holds for the case
that a_j≤ b_j and every d≥ 3, by letting r=1, i_1 = j and x' = a_j/a_j.

Case (b): a_j > b_j. Put i∈ [n].
Assume without loss of generality that a_j/a_j = e_d = (0,⋯,0,1)∈^d. Otherwise, rotate the coordinates system.
Hence,

    _j     = x∈^d-1|a_j^T/a_jx = b_j/a_j
        = x∈^d-1| e_d^Tx = b_j/a_j
        = x = (x^1,⋯,x^d)∈^d-1| x^d = b_j/a_j,

where the first derivation holds by substituting i=j in (<ref>), and the second derivation holds by the assumption that a_j/a_j = e_d.

For every m ∈ [n], identify the entries of a_m by (a_m^1,⋯,a_m^d), let

    b_m' = -(a_m^d · b_j/a_j - b_m),

and

    a_m' = √(1-(b_j/a_j)^2)· (a_m^1,⋯,a_m^d-1)^T ∈^d-1.

For every x = (x^1,⋯,x^d) ∈_j, let

    v(x) = (x^1,⋯,x^d-1)^T/(x^1,⋯,x^d-1)∈^d-2,

and observe that (x^1,⋯,x^d-1)≠ 0 since x=1 and |x^d| = |b_j|/a_j < 1.
For every x∈_j we have that

    |a_i^Tx-b_i|
       = | (a_i^1,⋯,a_i^d-1)(x^1,⋯,x^d-1)^T + a_i^d· x^d - b_i| 
        = | (a_i^1,⋯,a_i^d-1)(x^1,⋯,x^d-1)^T + a_i^d ·b_j/a_j - b_i| 
        = | (a_i^1,⋯,a_i^d-1)(x^1,⋯,x^d-1)^T - b_i'| 
        = | √(1-(b_j/a_j)^2)· (a_i^1,⋯,a_i^d-1)(x^1,⋯,x^d-1)^T/√(1-(b_j/a_j)^2) - b_i' | 
        = | √(1-(b_j/a_j)^2)· (a_i^1,⋯,a_i^d-1)(x^1,⋯,x^d-1)^T/(x^1,⋯,x^d-1) - b_i' | 
        = | a_i'^T v(x) -b_i' | ,

where (<ref>) is by (<ref>), (<ref>) is by the definition of b_i', (<ref>) holds since x is a unit vector, and (<ref>) is by the definition of a_i'.

We continue the proof of Case (b) with the following subcase analysis: b(i) a_i'= 0⃗ for every i∈ [n]∖j, and b(ii) a_l'≠0⃗(d-1) for some index l∈ [n]∖j.

Subcase b(i): a_i'= 0⃗(d-1) for every i∈ [n]∖j.
Put x' ∈_j and i∈ [n]∖j. In this case, for every x̂∈^d-2 we have that

    |a_i'x̂-b_i'| = |b_i'|.

We thus obtain that

    |a_i^Tx'-b_i| = |a_i'^Tv(x')-b_i'| = |b_i'|,

where the first equality is by (<ref>) and the last equality is by (<ref>).

Combining (<ref>) and the fact that x_j ∈_j yields

    |a_i^Tx'-b_i| = |a_i^Tx_j-b_i|.

Combining the last equality with (<ref>) yields

    |a_i^Tx'-b_i| ≤ 4· |a_i^Tx^*-b_i|.


Combining that the last equality holds for every x'∈_j and i∈ [n], with the fact that |_j| = ∞ by (<ref>) proves Subcase b(i) of Lemma <ref> for X = _j.



Subcase b(ii): a_l' ≠0⃗(d-1) for some l∈ [n]∖j.
We have that

    |a_i'^Tv(x_j)-b_i'| = |a_i^Tx_j-b_i| ≤ 4 · |a_i^Tx^*-b_i|,

where the equality follows by substituting x=x_j in (<ref>), and the inequality holds by (<ref>).

Let A' = (a_1' |⋯| a_j-1' | a_j+1' |⋯,a_n')^T∈^(n-1)×(d-1), b' = (b_1',⋯,b_j-1',b_j+1',⋯,b_n')^T ∈^n-1.
Applying the inductive assumption on A',b' and v(x_j) yields that there exists a set i_1,⋯,i_m⊆ [n]∖j of
m ∈ [d-2] indices such that for ' = ((a_i_1'|⋯| a_i_m'),(b_i_1',⋯,b_i_m')) and every i∈ [n]∖j there is x̂∈' that satisfies

    |a_i'x̂-b_i'| ≤ 4^d-2· |a_i'v(x_j)-b_i'|.

Moreover, (<ref>) holds for every x̂∈ X' if |'|=∞.

Let

    = ((a_j| a_i_1|⋯| a_i_m),(b_j,b_i_1,⋯,b_i_m)).


We continue to prove Lemma <ref> for Subcase b(ii) using another case analysis. Subcase b(ii,1): |'| ∈ O(1) and there is x̂∈' that satisfies (<ref>), and Subcase b(ii,2): |'| = ∞ and every x̂∈' satisfies (<ref>).

Subcase b(ii,1): |'| ∈ O(1) and there is x̂∈' that satisfies (<ref>). We prove there is x'∈ that satisfies (<ref>) and that || ∈ O(1).

Let

    x' = (√(1-(b_j/a_j)^2)·x̂^T |b_j/a_j)^T ∈^d.

Observe that

    v(x') = x̂∈',

by the definitions of x' and x̂ respectively.
Combining the definition of x' and (<ref>) yields that x' ∈_j ⊆^d-1. 
Therefore,

    |a_i^Tx'-b_i|     = |a_i'^Tv(x')-b_i'| 
        = |a_i'^Tx̂-b_i'| 
       ≤ 4^d-2· |a_i'^Tv(x_j)-b_i'| 
       ≤ 4^d-1· |a_i^Tx^*-b_i|.

where (<ref>) holds by substituting x=x' in (<ref>), (<ref>) is by (<ref>), (<ref>) is by (<ref>),  and (<ref>) holds by (<ref>).

We now present the following observation that will be used afterwards in the proof.
 


    = x' ∈_j | v(x') ∈'.




    =((a_j| a_i_1|⋯| a_i_m),(b_j,b_i_1,⋯,b_i_m))
       =_x'∈^d-1: 
    a_j^Tx'=b_j |a_i_m^Tx'-b_i_m|,     m=1 
    _x'∈^d-1: 
    (a_j | a_i_1|⋯| a_i_m-1)^Tx'=
    (b_j, b_i_1,⋯,b_i_m-1)^T |a_i_m^Tx'-b_i_m|,    
       =_x'∈^d-1: 
    a_j^Tx'=b_j |a_i_m^Tx'-b_i_m|,     m=1 
    _x'∈^d-1: 
    a_j^Tx'=b_j,
    
    (a_i_1|⋯| a_i_m-1)^Tx'=
    (b_i_1,⋯,b_i_m-1)^T |a_i_m^Tx'-b_i_m|,    
       =_x'∈^d-1: 
    a_j^T/a_jx'=b_j/a_j |a_i_m^Tx'-b_i_m|,     m=1 
    _x'∈^d-1: 
    a_j^T/a_jx'=b_j/a_j,
    
    (a_i_1|⋯| a_i_m-1)^Tx'=
    (b_i_1,⋯,b_i_m-1)^T |a_i_m^Tx'-b_i_m|,    
       =_x'∈^d-1: 
    
    x' ∈_j |a_i_m^Tx'-b_i_m|,     m=1 
    _x'∈^d-1: 
    
    x' ∈_j,
    
    (a_i_1|⋯| a_i_m-1)^Tx'=
    (b_i_1,⋯,b_i_m-1)^T |a_i_m^Tx'-b_i_m|,    
       =_x'∈^d-1: 
    
    x' ∈_j |a_i_m'^Tv(x')-b_i_m'|,     m=1 
    _x'∈^d-1: 
    
    x' ∈_j,
    
    (a_i_1' |⋯| a_i_m-1')^Tv(x')=
    (b_i_1',⋯,b_i_m-1')^T |a_i_m'^Tv(x')-b_i_m'|,    
       =x' ∈_j | v(x') ∈((a_i_1'|⋯| a_i_m'),(b_i_1',⋯,b_i_m'))
       =x' ∈_j | v(x') ∈'

where (<ref>) holds by (<ref>),  (<ref>) holds by  (<ref>) and since there are at least 2 entries in (b_j,b_i_1,⋯,b_i_m),  (<ref>) holds since it's the same as (<ref>) just written differently. (<ref>) holds since a_j > b_j≥ 0 and therefore there's no division by 0. (<ref>) holds by  (<ref>). (<ref>) holds by (<ref>). (<ref>) holds by  (<ref>) and (<ref>) holds by the definition of '.


Combining the fact that x' ∈_j and (<ref>) with Observation <ref> yields that

    x' ∈.

It also holds that

    || = |x' ∈_j | v(x') ∈'|
        = |x'=(x^1,⋯,x^d) | x^d=b_j/a_j and  v(x') ∈'|
       ≤ |'| ∈ O(1),

where the last derivation is by the assumption of Subcase b(ii,1).

Hence, Lemma <ref> holds for Subcase b(ii,1) with X = and x' ∈ X.

Subcase b(ii,2): |'| = ∞ and every x̂∈' satisfies (<ref>). We prove that || = ∞ and that every x'∈ satisfies (<ref>).




Similarly to (<ref>), we have that || ≤ |'| = ∞.
For every x' ∈, by Observation <ref> we have that x' ∈_j and v(x') ∈'. Therefore, for every x' ∈ we have that

    |a_i^Tx'-b_i| = |a_i'^Tv(x')-b_i'|   ≤ 4^d-2· |a_i'^Tv(x_j)-b_i'|
       ≤ 4^d-1· |a_i^Tx^*-b_i|.

where the first equality holds by substituting x=x' in (<ref>), the first inequality holds by substituting x̂ = v(x') ∈' in (<ref>), and the last inequality is by (<ref>).

Hence, Lemma <ref> holds for Subcase b(ii,2) with X =.

Combining Subcase b(ii,1) and Subcase b(ii,2) prove Lemma <ref> for Subcase (ii). Combining Subcase (i) and Subcase (ii) proves  Lemma <ref> holds for Case (b) for every d≥ 3.
Lemma <ref> now holds for every d≥ 3 by combining Case (a) and Case (b).



 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix of n ≥ d-1 ≥ 1 rows, and let b = (b_1,⋯,b_n)^T ∈^n.
Let X ⊆^d-1 be the output of a call to (A,b); see Algorithm <ref>. Then for every x^* ∈^d-1 there is a unit vector x' ∈^d-1 such that for every i∈[n],

    |a_i^Tx'-b_i| ≤ 4^d-1· |a_i^Tx^*-b_i|.

Moreover, the set X can be computed in n^O(d) time and its size is |X| ∈ n^O(d).


We use the variables as in Algorithm <ref> for a call to (A,b).

Put x^* ∈^d-1 and i∈ [n].
In Lines <ref>–<ref> of Algorithm <ref>, we define

    b_i' = |b_i|,
         a_i' = (b_i)· a_i.

Observe that for every a∈^d, b̂∈ and every x∈^d

    |a^Tx-b̂| = |-a^Tx-(-b̂)|.


Substituting A with A'=a_1',⋯,a_n', and b with b'=b_1',⋯,b_n' and plugging x^* in Lemma <ref> yields that there is a set j_1,⋯,j_r⊆ [n] of m ∈ [d-1] indices, such that for X' = ((a_j_1'|⋯| a_j_m'), (b_j_1',⋯,b_j_m')) and every i∈ [n], there is x̂∈ X' that satisfies

    |a_i'^Tx̂-b_i'| ≤ 4^d-1· |a_i'^Tx^*-b_i'|.

Moreover, if |X'| = ∞ then every x̂∈ X' satisfies (<ref>).

Let x' be x̂ if |X'| ∈ O(1), and x' ∈ X' be an arbitrary element if |X'| = ∞. Hence, x' satisfies (<ref>).
We now have that

    |a_i^Tx'-b_i| = |a_i'^Tx'-b_i'|    ≤ 4^d-1· |a_i'^Tx^*-b_i'|
       = 4^d-1· |a_i^Tx^*-b_i|,

where the first and last equalities hold by combining the definitions of a_i' and b_i' with (<ref>), and the inequality holds since x' satisfies (<ref>).

It is left to prove that x' is in the output set X.
In Lines <ref>–<ref> of Algorithm <ref>, we iterate over every r ∈ [d-1] and every subset i_1,⋯,i_r∈ [n] of r indices and compute the set S = ((a_i_1'|⋯| a_i_r'),(b_i_1',⋯,b_i_r')) using Algorithm <ref>.

Therefore, when r=m, and i_1=j_1,⋯,i_r = j_r, the call to Algorithm <ref> in Line <ref> is guaranteed to compute a set S that satisfies

    S =
         X',     |X'|∈ O(1) 
    x̂∈ X',    .
       .

We then add S to the output set X. By the definition of x' and the output guarantee of Algorithm <ref>, x' ∈ S. Hence, the output set X is guaranteed to contain an element that satisfies (<ref>).

The running time of Algorithm <ref> is dominated by the nested d·nd∈ n^O(d) ”for“ loops, each iteration taking O(d) time, for a total of n^O(d) running time.
In each iteration, we add O(1) unit vectors to X. Hence, |X| ∈ n^O(d).


 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix of n ≥ d-1 ≥ 1 rows, and let b = (b_1,⋯,b_n)^T ∈^n. Let , s, r be as defined in Definition <ref> for Y = (a_i,b_i)| i∈ [n] and D((a,b̂),x) = |a^Tx-b̂| for every a∈^d, b̂∈ and x∈^d. Then in n^O(d) time we can compute a unit vector x'∈^d-1 such that

    (Y,x') ≤ 4^(d-1)rs·min_x∈^d-1(Y,x).



Let x^* ∈_x∈^d-1(Y,x).
Plugging A,b and x^* in Theorem <ref> yields that in n^O(d) time we can compute a unit vector x'∈^d-1 such that for every i∈ [n]

    D((a_i,b_i),x')=|a_i^Tx'-b_i|    ≤ 4^d-1· |a_i^Tx^*-b_i|
       = 4^d-1· D((a_i,b_i),x^*).

Theorem <ref> now holds by plugging Y, D, , r,s, q'=x', q^*=x^* and c=4^d-1 in Observation <ref> as

    (Y,x')    ≤ 4^(d-1)rs(Y,x^*)
       = 4^(d-1)rsmin_x∈^d-1(Y,x).






§ HANDLING UNKNOWN MATCHING


 
Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix containing n ≥ d-1 ≥ 1 rows, let b = (b_1,⋯,b_n)^T ∈^n, and let D((a,b̂),x) = |a^Tx-b̂| for every a∈^d, b̂∈ and x∈^d. Let , r be as defined in Definition <ref> for Y = (a_i,b_i)| i∈ [n], D and f(v) = v_1. Let (x̃,) be a pair of unit vector and permutation (matching function) which is the output of a call to (A,b,); see Algorithm <ref>. Then it holds that

    (Y_,x̃) ≤ 4^(d-1)r·min_x,(Y_,x),

where the minimum is over every unit vector x∈^d-1 and ∈(n). Moreover, (x̃,) is computed in n^O(d) time.


Put (x^*,^*) ∈_x,(Y_,x), where the minimum is over every unit vector x∈^d-1 and permutation ∈(n), and put i∈ [n].

In Lines <ref>-<ref> of Algorithm <ref>, we define

    b_i' = |b_i|,
         a_i' = (b_i)· a_i.

For every a∈^d, b̂∈ and every x∈^d-1, we have

    |a^Tx-b̂| = |-a^Tx-(-b̂)|.


Substituting A with A'=a_1',⋯,a_n', b with b'=b_^*(1)',⋯,b_^*(n)' and plugging x^* in Lemma <ref> yields that there is a set j_1,⋯,j_m⊆ [n] of m ∈ [d-1] indices, such that for X' = ((a_j_1'|⋯| a_j_m'), (b_^*(j_1)',⋯,b_^*(j_m)')) and every i∈ [n], there is x̂∈ X' that satisfies

    |a_i'^Tx̂-b_M^*(i)'| ≤ 4^d-1· |a_i'^Tx^*-b_M^*(i)'|.

Moreover, if |X'| = ∞ then every x̂∈ X' satisfies (<ref>).








Let x' be x̂ if |X'| ∈ O(1), and x' ∈ X' be an arbitrary element if |X'| = ∞. Hence, x' satisfies (<ref>).
We now have

    |a_i^Tx'-b_^*(i)|     = |a_i'^Tx'-b_^*(i)'|
       ≤ 4^d-1· |a_i'^Tx^*-b_^*(i)'|
       = 4^d-1· |a_i^Tx^*-b_^*(i)|.

where the first and last equalities hold by combining the definitions of a_i' and b_i' with (<ref>), and the inequality holds since x' satisfies (<ref>).

Since (<ref>) holds for every i∈ [n], we can substitute Y = Y_^* = (a_i,b_^*(i))_i=1^n, q' = x', q^*=x^* and c=4^d-1 in Observation <ref> to obtain that

    (Y_^*,x') ≤ 4^r(d-1)·(Y_^*,x^*).


It is left to prove that x' is in the output set X.
In Lines <ref>-<ref> of Algorithm <ref>, we iterate over every r ∈ [d-1] and every subsets i_1,⋯,i_r⊆ [n] and ℓ_1,⋯,ℓ_r⊆ [n] of r indices, we compute the set S = ((a_i_1'|⋯| a_i_r'),(b_ℓ_1',⋯,b_ℓ_r')) using Algorithm <ref>, we then add S to X.

Therefore, when r=m, i_1=j_1,⋯,i_r = j_r and ℓ_1 = ^*(j_1),⋯,ℓ_r = ^*(j_r), the call to Algorithm <ref> in Line <ref> is guaranteed to compute a set S that satisfies

    S =
         X',     |X'|∈ O(1) 
    x̂∈ X',    .
       .

We then add S to the output set X. By the definition of x' and the output guarantee of Algorithm <ref>, x' ∈ S. Hence, the output set X is guaranteed to contain an element that satisfies (<ref>).

Afterwards, in Line <ref>, we compute the optimal permutation (Y,x,) for every unit vector x∈ X.
To compute the optimal matching one could compute the pairwise cost matrix Z ∈^n× n, where the entry (i,j) contains ((a_i,b_j),x) (i.e., the cost of pairing a_i with b_j when using x), and then apply the Hungarian Method <cit.>.

Since there is x' ∈ X that satisfies (<ref>), we get that (x',(Y,x',)) ∈ S. Let M' = (Y,x,).
By Definition <ref> of optimal matching,

    (Y_M',x') ≤min_∈(n)(Y_,x') ≤(Y_M^*,x').

In Line <ref>, we pick the pair (x̃,) ∈_(x,)∈ S(Y_,x). Therefore, the main claim of Theorem <ref> holds as

    (Y_,x̃) ≤(Y_',x')   ≤(Y_M^*,x')
       ≤ 4^r(d-1)·(Y_^*,x^*),

where the first inequality is by the definition of (x̃,), the second inequality is by (<ref>) and the last inequality is by (<ref>).

Since the call to Algorithm <ref> in Line <ref> returns at most 2 unit vectors, the size of the set X is proportional to the number of iterations of the nested for loops, which is d·nd^2 ∈ n^O(d).
The running time of Algorithm <ref> is dominated by Line <ref>, which runs an optimal matching algorithm for every unit vector in X.
The Hungarian method for optimal matching takes O(n^3) time. Hence, the total running time is n^O(d).





§ CORESET FOR LINEAR REGRESSION



Let M be an n× m matrix of rank ρ. Let z∈ [1,∞), and α,β≥ 1. An n×ρ matrix U is an (α,β,z)-conditioned basis
for M if the column vectors of U span the column space of M, and additionally U satisfies that: (1) ∑i,j |a_i,j|^z ≤α^z, (2) for all u∈^ρ,
u_z'≤βUu_z, where ·_z' is the dual norm for ·_z (i.e. 1/z + 1/z'=1).


 
Let M = (m_1|⋯| m_n)^T be an n× m matrix of rank ρ. Let z∈ [1,∞). Let U = (u_1|⋯| u_n)^T be an (α,β,z)-conditioned basis for M. For every vector u∈^m, the following inequality
holds: |m_i^Tu|^z ≤ (u_i_z^z ·β^z)Mu_z^z.


The following is a restatement of Lemma 16 in <cit.>
 
Let P = p_1,⋯,p_n⊆^d be a set of n ≥ d points. Let = (x,b) ∈^d×|x=1 be the set of all hyperplanes in ^d and let z ∈ [1,∞).
Let s:P → [0,∞) where s(p_i) := sup_(x,b)∈|p_i^Tx-b|^z/∑_j∈ [n]|p_j^Tx-b|^z is the sensitivity of a point p_i ∈ P.
Then the total sensitivity ∑_i∈ [n] s(p_i) of P is O(d^1+z/2) if z∈ [1,2), O(d) for z=2 and O(d^z) for z>2.



 
Let P = p_1,⋯,p_n⊆^d be a set of n ≥ d points and let z∈ [1,∞). Then a function s:P → [0,∞) can be computed in O(nd^5 logn)
time such that for every i∈ [n],

    s(p_i) ≥sup_(x,b) ∈|p_i^T x-b|^z/∑_j ∈ [n]|p_j^Tx-b|^z,

where the sup is over every (x,b) ∈ such that ∑_j ∈ [n]|p_j^Tx-b|^z ≠ 0.


For every i∈ [n], let M_i = (p_i^T | 1)^T, and let M = (m_1 |⋯| m_n)^T ∈^n× (d+1). Let U = (u_1|⋯| u_n)^T be a (α,β,z)-conditioned basis for M.
Put i∈ [n]. Then for every unit vector x∈^d and b ≥ 0

    |p_i^T x-b|^z = |(p_i^T | 1) (x^T| b)^T|^z = |m_i^T (x^T| b)^T|^z,

and

    ∑_j ∈ [n]|p_j^Tx-b|^z     = ∑_j ∈ [n]|(p_j^T | 1) (x^T| b)^T|^z
        = ∑_j ∈ [n]|m_j^T (x^T| b)^T|^z
       = M (x^T| b)^T_z^z.


It thus holds that

    sup_(x,b) ∈|p_i^T x-b|^z/∑_j ∈ [n]|p_j^Tx-b|^z    = sup_(x,b) ∈|m_i^T (x^T| b)^T|^z/M (x^T| b)^T_z^z
       ≤u_i_z^z ·β^z,

where the first equality is by (<ref>) and (<ref>), and the inequality holds by plugging M, U and u = (x^T| b)^T in Lemma <ref>.
Hence, by letting s(p_i) = u_i_z^z ·β^z, we get that

    s(p_i) ≥sup_(x,b) ∈|p_i^T x-b|^z/∑_j ∈ [n]|p_j^Tx-b|^z.


The time it takes to compute the function s is dominated by the time to compute the matrix U. By Theorem 3.1 in <cit.>,
it takes O(nd^5logn) time to compute the (α,β,z)-conditioned basis matrix U for M since M is an n × (d+1) matrix of rank at most d+1.
Hence, it takes O(nd^5logn) time to compute the function s.


 
Let P be a finite set, and let w:P→ [0,∞). Let Q be a function that maps every set S⊆ P to a corresponding set Q(S), such that Q(T) ⊆ Q(S) for every T⊆ S. Let f:P× Q(P) → be a cost function. The tuple (P,w,Q,f) is called a query space.


 
For a query space (P,w,Q,f), S⊆ P, q∈ Q(S) and r∈ [0,∞) we define (q,r) = p∈ P | w(p)· f(p,q) ≤ r.
The VC-dimension of (P,w,Q,f) is the smallest integer d_VC such that for every S⊆ P we have

    | (q,r) | q∈ Q(S), r∈ [0,∞)| ≤ |S|^d_VC.



 
Let (P,w,Q,f) be a query space; see Definition <ref>. Let s:P→ [0,∞) such that

    sup_qw(p)f(p,q)/∑_p∈ P w(p)f(p,q)≤ s(p),

for every p∈ P and q∈ Q(P) such that the denominator is non-zero. Let t = ∑_p∈ P s(p) and Let d_VC be the VC-dimension of query space (P,w,Q,f); See Definition <ref>. Let c ≥ 1 be a sufficiently large constant and let ε, δ∈ (0,1). Let S be a random sample of

    |S| ≥ct/ε^2(d_VClogt+log1/δ)

points from P, such that p is sampled with probability s(p)/t for every p∈ P. Let u(p) = t· w(p)/s(p)|S| for every p∈ S. Then, with probability at least 1-δ, for every q∈ Q it holds that

    (1-ε)∑_p∈ P w(p)· f(p,q)    ≤∑_p∈ S u(p)· f(p,q)
       ≤ (1+ε)∑_p∈ P w(p)· f(p,q).




Let d ≥ 2 be a constant integer. Let A = (a_1 |⋯| a_n)^T ∈^n× d be a matrix containing n ≥ d-1 points in its rows, let b = (b_1,⋯,b_n)^T ∈^n, and let w = (w_1,⋯,w_n) ∈ [0,∞)^n.
Let ε, δ∈ (0,1) and let z∈ [1,∞).
Then in O(n logn) time we can compute a weights vector u = (u_1,⋯,u_n)∈[0,∞)^n that satisfies the following pair of properties.



  * With probability at least 1-δ, for every x∈^d-1 it holds that

    (1-ε) ·∑_i∈ [n] w_i ·    |a_i^Tx-b_i|^z  ≤∑_i∈ [n]u_i· |a_i^Tx-b_i|^z
       ≤ (1+ε) ·∑_i∈ [n] w_i · |a_i^Tx-b_i|^z.


  * The weights vector u has O(log1/δ/ε^2) non-zero entries.



Put x∈^d-1 and i∈ [n].
Let W ∈^n× n be a diagonal matri whose diagonal elements are (w_1^1/z,⋯,w_n^1/z).
Let A' = W A ∈^n× d, b' = Wb ∈^n and P = (p_1|⋯| p_n)^T= [A' | b'] ∈^n× (d+1). For every unit vector x∈^d it holds that

    ∑_i∈ [n] w_i · |a_i^Tx-b_i|^z = A'x-b'_z^z = P (x^T | 1)^T_z^z.


Plugging P in Corollary <ref> yields that in O(nd^5logn) time we can compute a function s:P→ [0,∞) such that

    s(p_i) ≥sup_(y,b) ∈^d×|p_i^T y-b|^z/P y-b_z^z.


Therefore, we have that

    sup_x ∈^d-1   w_i · |a_i^Tx-b_i|^z/∑_j∈ [n] w_j · |a_j^Tx-b_j|^z
    = sup_x ∈^d-1|p_i^T (x^T | 1)^T|^z/P (x^T | 1)^T_z^z
       ≤sup_y ∈^d√(2) |p_i^T y|^z/√(2)P y_z^z
       ≤sup_(y,b) ∈^d×|p_i^T y-b|^z/P y-b_z^z
       ≤ s(p_i).

where (<ref>) holds by the definition of P, (<ref>) holds since (x^T | 1)^T/√(2)|x = 1⊆^d,
(<ref>) holds since for every y ∈^d there is (a,b) ∈^d× such that a=y, and (<ref>) is by (<ref>).

The VC-dimension of the corresponding query space is bounded by d_VC∈ O(d) by <cit.>; See Definition <ref>.

Let t = ∑_i∈ [n] s(p_i) and Y = (a_i,b_i)_i=1^n.
By Lemma <ref> we have that t ∈ O(d^1+z/2) if z∈ [1,2), O(d) if z=2 and O(d^z) if z>2.

Let Z ⊆ Y be a random sample of

    |Z| ∈ O(t/ε^2(dlogt+log1/δ))
    =d^O(1)/ε^2log1/δ

pairs from Y, where (a_i,b_i) ∈ Y is sampled with probability s(p_i)/t and let u = (u_1,⋯,u_n) where

    u_i = t· w_i/s(p_i)|Z|,     if  (a_i,b_i)∈ Z
    
    0,     otherwise.

By substituting P=Y, Q(·)≡^d-1, f((a,b),x)=w_i· |a^Tx-b|^z for every (a,b)∈ Y and x∈^d-1, t ∈ d^O(1) and d_VC=O(d), in Theorem <ref>, Property (i) of Theorem <ref> holds as

    (1-ε) ·∑_i∈ [n]    w_i · |a_i^Tx-b_i|^z ≤∑_i∈ [n]u_i· |a_i^Tx-b_i|^z
       ≤ (1+ε) ·∑_i∈ [n] w_i · |a_i^Tx-b_i|^z.


Furthermore, Property (i) of Theorem <ref> holds since the number of non-zero entries of u is equal to |Z| ∈d^O(1)/ε^2log1/δ = O(log1/δ/ε^2).

The time needed to compute u is bounded by the computation time of s, which is bounded by O(nlogn) since d is a constant.














































