Subgradient and Sampling Algorithms for ℓ_{1} Regression

Kenneth L. Clarkson^{∗}

January 4, 2005

Abstract Given an n × d matrix A and an n-vector b, the ℓ_{1} regression problem is to find the vector x minimizing the objective function ∥Ax−b∥_{1}, where ∥y∥_{1} ≡ ∑_{i} |y_{i}| for vector y. This paper gives an algorithm needing O(n log n)d^{O(1)} time in the worst case to obtain an approximate solution, with objective function value within a fixed ratio of optimum. Given ϵ > 0, a solution whose value is within 1 + ϵ of optimum can be obtained either by a deterministic algorithm using an additional O(n)(d/ϵ)^{O(1)} time, or by a Monte Carlo algorithm using an additional O((d/ϵ)^{O(1)}) time. The analysis of the randomized algorithm shows that weighted coresets exist for ℓ_{1} regression. The algorithms use the ellipsoid method, gradient descent, and random sampling.

1 Introduction

Given a set S of n points, the ℓ_{1} regression problem is to fit a hyperplane to the points, minimizing the sum of the vertical distances of the points to the hyperplane. More formally, given an n × d matrix A and an n-vector b, the problem is to find the d-vector x̂ such that the sum of the absolute values of the entries of Ax − b is minimum. The entries of Ax − b are called residuals. That is, x̂ minimizes the objective function ∥Ax − b∥_{1}, yielding the minimum residual vector b̂ ≡ Ax̂ − b. Another description of x̂ is that it gives the linear combination of the columns of A that is closest in ℓ_{1} distance to b.

This problem arises in statistics, as a more robust alternative to least squares regression, for which the sum of the squares of the residuals is minimized. We can also regard the residual as the difference between b_{i} and its prediction a_{i·}x, where a_{i·} is the i’th row of A. (The points of S are of the form [a_{i·}b_{i}], where a_{i·} is the i’th row of A, so “vertical” here would refer to the last coordinate.)

Related work. While ℓ_{1} regression is an instance of linear programming, and so has good algorithms in theory and practice, it is still of interest to further understand its computational properties. Moreover, there is the usual gap between theory and practice: the best provably-good algorithms take O(n) time, and are exact,[YKII88, MT93] but are complicated, and take time at least exponential in d. (The stated dependence is 3^{d^{2}},[MT93] but a dependence that is single-exponential is likely to be achievable with current methods.) In contrast, the algorithms here, while approximate and O(n log n), are polynomial in d and 1/ϵ.

While interior point or simplex methods can be applied, and may be satisfactory in practice, their worst-case dependence on n may be substantial, since the linear programming formulation has 2n constraints and n + d unknowns.[Sch86] (While the dual formulation allows an interior point iteration in nd^{O(1)} time, the number of such iterations can be large in the worst case. Also, the ellipsoid method could be applied using the subgradient, as used here for rounding.)

∗Bell Labs; 600 Mountain Avenue; Murray Hill, New Jersey 07974; clarkson@research.bell-labs.com

Subgradient descent. A key algorithmic element used here is a version of gradient descent: starting at an estimate of the optimum, move in the direction of the negative gradient of the objective function

F (x) ≡ ∥Ax − b∥_{1},

to obtain a new estimate. Each such step has a certain step length, chosen judiciously (but easily). The descent direction at d-vector x used here is

g_{x} ≡ −A^{T} sgn(Ax − b),

where sgn(y), for a vector y, is the vector of equal dimension whose coordinates are the signs of the corresponding coordinates of y. (Note that the corresponding direction for the least squares objective function ∥Ax − b∥_{2} is −A^{T} (Ax − b); setting that function to zero yields the normal equations.)

Of course, the function ∥Ax−b∥_{1} is only piece-wise linear, and so does not have a gradient everywhere; however, the value g_{x} is a subgradient (as stated below in Lemma 3.1), and the iteration just described is an instance of an unconstrained “subgradient method.” Such methods have a long history and rich theory (see, e.g., [Sho85, NB01, Ber01]). However, their analyses seem to typically focus on limiting convergence behavior: as the number of iterations goes to infinity, does the estimate converge to the optimum? Also, most analyses of subgradient methods focus on a fixed sequence of stepsizes, where here the stepsize is computed using the current function value.

Conditioning. Here the results can be stated with bounds in terms of n, d, and ϵ, as in Theorem 3.7, by taking advantage a particular kind of “conditioning” of A for ℓ_{1} regression, related to polytope rounding.[Lov86] Elementary column operations and scaling on A amount to a change of variable for x. Since such changes of variable can be tracked readily, this paper will freely use column operations and scaling. As discussed in Section 2, such operations can be done so as to obtain a new version of the matrix A with the property that, for any x,

∥x∥_{1} ≥ ∥Ax∥_{1} ≥ ∥x∥_{1}/d\sqrt{d}.

A matrix A will this property will be called ℓ_{1} conditioned. This relation is similar to, but weaker than, the strict equality ∥Ax∥_{2} = ∥x∥_{2} that would hold if Gram-Schmidt orthogonalization were done on the columns of A. The conditioning is done by applying the ellipsoid algorithm, as used to find “weak Loewner-John ellipsoids.” The subgradient is used as a separation oracle. To obtain a good time bound for the conditioning, without a dependence on the bit-complexity of the input, we first apply elementary column operations to make the columns of A orthogonal, and make b orthogonal to A. That is, a bit of “folklore” advice for getting an initial estimate for ℓ_{1} regression, which is to begin with the optimum for ℓ_{2} regression, is used here.

The conditioning of A is helpful both for the ℓ_{1} subgradient method of Section 3, and for the sampling procedure discussed in Section 4. For the subgradient method, the conditioning assures that each iteration makes good progress. For the randomized sampling procedure, the conditioning allows a bound on the norm of the optimum vector x̂, which allows a bound on the variance of the objective function random variable considered in the analysis.

Sampling. The sampling algorithm is not much more than picking a random sample of S, and solving the ℓ_{1} regression problem for the resulting sample problem. (This again is not far from folklore for how to get an initial estimate for ℓ_{1} regression.) The sampling is done after some preprocessing, including the conditioning procedure, and the samples are chosen with probabilities that depend on the point coordinates. Moreover, the resulting sample problem has each point weighted. The sample problem has the provable property that the optimum x̂ of the original problem yields the expected F (x̂) value as in the original problem, and any point yielding a bad objective function value in the original problem also yields a bad value in the sample. (This is the outline of the proof of Theorem 4.7.) So the solution to the sample problem will be a good solution to the original.

The approximate hardness of the weighted sample problem is a property analogous to that of “coresets” for other fitting problems.[BHPI02, HPV02, AHPVar, BC02] A novelty here is that while most of the known coreset constructions are for problems where the fit involves minimizing a maximum of distances, here the fit is the sum of residuals. However, the “coreset” here is not apparently as useful, but may be one step toward coreset-based algorithms for problems related to ℓ_{1} regression. The size of the sample is only polynomially dependent on the dimension, as opposed to some coreset constructions with exponential dependence on the dimension. As shown in Theorem 4.7, the bound is O(d^{3}\sqrt{d}/ϵ^{2}) log(d/γϵ), where γ is a failure probability; this is similar in form to bounds for ϵ-approximations.

The probabilities used in the sampling technique is similar in flavor to that used in some recent randomized approximation algorithms for numerical linear algebra.[DK01, DKM] One difference from those algorithms is that the sampling is done independently for each input point, rather than picking each point of the sample independently from among the input points.

Detailed time bounds. The time bounds for the algorithms are, in more detail,

O(nd^{5}(log n log d + log(d/ϵ)/ϵ^{2}))

in the worst case, for the deterministic algorithm, and

O(nd^{5} log n log d) + O(\frac{d^{8}\sqrt{d}}{ϵ^{4}}(log d)^{3} ln(1/ϵ)^{2} ln(1/γ))

for the Monte Carlo algorithm, where the failure probability is γ. (Note the lack of dependence on n in the second term of the latter bound.) These bounds are from Theorem 3.7 and Theorem 4.8, respectively. Although these dependencies are polynomial, they are relatively large, but it is likely that faster performance is possible in practice. This is true particularly if the rounding/conditioning step is skipped; then the main algorithm comprises only the subgradient descent procedure. In that case O(nd^{2}) time is needed, iterated a data-dependent number of times.

Relation to small-d LP. The form of the time bound for the sampling algorithm is roughly similar to those long known for randomized algorithms for LP (linear programming) in small dimension.[Cla95] Such LP problems arise in, for example, ℓ_{∞} regression. The previous algorithms rely on the existence of a “coreset” of size d for the exact problem; that existence follows by basic arguments. Here the coreset existence is much harder to prove, and uses weights, and is only approximate.

Statistical view. Another motivation for studying this problem is that we might regard the regression hyperplane as a statistic, interpreting S itself as a random sample from an underlying population. We might hope that the expected value of the statistic for the sample is equal to the value of the statistic for the whole set, and indeed, an unweighted (unbiased) sample has the appropriate expectation. However, an unweighted sample doesn’t have the same provable properties, for clear reasons: while a single outlying point can change the regression plane by an arbitrarily large amount, that outlier may well not be in the sample. That is, the ℓ_{1} estimator is not that robust, and so unweighted samples can do poorly. Here the sampling method makes outliers more likely to be picked; this reduces the variance due to their presence.

Outline. Before describing the approximation algorithm in Section 3 and the sampling algorithm in Section 4, the conditioning procedure is described. There are some concluding remarks.

2 Rounding and Conditioning

We will apply a standard rounding procedure, using the ellipsoid method[Lov86], to the polytope P ≡ P (A) ≡ {x | ∥Ax∥_{1} ≤ 1}, which will have the effect of making A wellbehaved with respect to the ℓ_{1} norm.

This procedure requires a separation oracle, that for a given point x ∉ P, returns a hyperplane H that separates x and P. Here the subgradient provides such an oracle, since

∥Ay∥_{1} ≥ ∥Ax∥_{1} + (y − x)^{T} A^{T} sgn(Ax),

as shown in Lemma 3.1 below, so that the hyperplane

{y | (y − x)^{T} A^{T} sgn(Ax) = ((∥Ax∥_{1}) + 1)/2}

separates x and P when ∥Ax∥_{1} > 1.

The rounding procedure maintains a series of ellipsoids E_{m}, m = 0, …, where E_{0} is a Euclidean ball B containing P. The procedure uses the separation oracle to maintain the condition that each E_{m} contains P and has provably smaller volume than E_{m−1}. A critical quantity is the ratio of the volume of B to the volume of P; the running time is proportional to the logarithm of this ratio.

We can construct an instance of the rounding problem where the ratio is a function only of n and d, as follows. Use Gram-Schmidt, or an equivalent procedure, to make the columns of A orthogonal to each other. Also, scale the columns of A to have ℓ_{1} norm equal to one. This needs only elementary column operations, and as noted above, amounts to a change of variable. The following lemma applies.

Lemma 2.1 Given an n × d matrix A where the columns of A are orthogonal and scaled to have unit ℓ_{1} norm,

B_{1/\sqrt{d}} ⊂ P (A) ⊂ B\sqrt{n},

where B_{r} denotes the Euclidean ball centered at the origin, with radius r.

Proof: If ∥x∥_{2} ≤ 1/\sqrt{d}, then ∥x∥_{1} ≤ 1, which implies ∥Ax∥_{1} ≤ 1 by the column scaling, giving the first inclusion.

For the second inclusion, for x with ∥Ax∥_{1} ≤ 1 we have

∥x∥_{2}^{2} = ∑_{k} x_{k}^{2}∥a_{·k}∥_{1}^{2}

≤ n ∑_{k} x_{k}^{2}∥a_{·k}∥_{2}^{2}

= n∥Ax∥_{2}^{2} ≤ n∥Ax∥_{1}^{2} ≤ n, which implies x ∈ B\sqrt{n}.

Lemma 2.2 As applied to the preprocessed n × d matrix A, the rounding procedure takes

O(nd^{5} log n)

time, and returns a pair of concentric ellipsoids E and E′ with E′ ⊂ P (A) ⊂ E, and E′ arises from E by shrinking by a factor of 1/d.

Proof: From the previous lemma, the ratio of the volumes of the initial ellipsoid to the final ellipsoid is at most (dn)^{d/2}, and the ellipsoids shrink by a factor of e^{−3/2(d+1)(2d+1)} at each step,[Lov86] so the number of steps is O(d^{3} log(nd)). A factor d^{2}n accounts for the work at each step of checking inclusions, generating separating hyperplanes, and linear algebra to maintain the ellipsoids. Since P is centrally symmetric, the tighter shrinkage term 1/d can be used, instead of the general 1/d\sqrt{d} for such ellipsoids.

The result of this rounding procedure can be more conveniently used here by applying another change of variables to yield the condition that A is “μ-conditioned” for μ = d\sqrt{d}, where this means that for any x,

∥x∥_{1} ≥ ∥Ax∥_{1} ≥ ∥x∥_{1}/μ.  (1)

Theorem 2.3 Given an n × d matrix A, the rounding procedure and elementary column operations on A yield a new version of A that is d\sqrt{d}-conditioned, in time

O(nd^{5} log n).

Proof: Having applied the rounding procedure, there is an affine transformation with matrix τ that maps the containing ellipsoid E to B_{d}, and the contained ellipsoid E′ to B_{1}. That is, if y = τ(x) for some x, and ∥y∥_{1} ≤ 1, then ∥y∥_{2} ≤ 1, and so y ∈ B_{1} ⊂ τP. Moreover, if y ∈ τP, then y ∈ B_{d} and ∥y∥_{2} ≤ d, so ∥y∥_{1} ≤ d\sqrt{d}. Since τP = P (Aτ^{−1}), renaming Aτ^{−1} to be A yields the conditions that ∥y∥_{1} ≤ 1 implies ∥Ay∥_{1} ≤ 1, and ∥Ay∥_{1} ≤ 1 implies ∥y∥_{1} ≤ d\sqrt{d}. These facts imply the inequalities (1) with μ = d\sqrt{d}.

3 An Approximation Algorithm

As shown below in Lemma 3.2, the rounding procedure has the useful effect that it results in an optimization problem for which the negative gradient of the objective function points from the current location toward an optimal vector. This and other properties will be stated and proven in the analysis of the subgradient iteration scheme.

This property of the gradient immediately suggests a descent scheme: take a step along the negative gradient toward an optimal vector, re-evaluate the gradient, and repeat. A complication here is that the although the negative gradient points toward an optimum vector, it is difficult to tell how far to go in that direction without overshooting the best possible improvement. It’s also difficult to know just how much an improvement will occur.

One approach might be to do a search along the line of the gradient, approximately finding the minimum of the objective function when confined to the given line. This is likely to work well in practice, but seems difficult to analyze.

Here we adopt an approach similar to, but more elaborate than, the smallest ball algorithm of Badiou et al.[BC02]. The step length is calculated using the current function value and other inputs, with a scheme that allows lower and upper bounds on the approximation ratio ∥Ax − b∥_{1}/∥Ax̂ − b∥_{1} to be maintained, where again x̂ is an optimal vector.

3.1 The subgradient iteration

The approximation algorithm is simply the subgradient iteration procedure described below, under the assumption that A has already been μ-conditioned with μ = d\sqrt{d}, as discussed in the last section, and that b is orthogonal to the subspace spanned by the column vectors of A. The latter is easy to achieve by using the orthogonalization of A produced while A was being conditioned.

The subgradient iteration procedure uses several readily calculated values, which will be specified below. The procedure is as follows: for each entry in a decreasing sequence of values α̃_{i}, for i = 0, 1, …, perform a minor iteration with α̃_{i} until α̃_{i} ≤ 1 + ϵ, where 1 + ϵ is the desired approximation ratio. As a special case for i = 0, repeatedly perform minor iterations with α̃_{0} until a minor iteration is unsuccessful, as described below.

Each minor iteration using some α̃ starts with the current estimate x of the solution, and iteratively attempts to improve on x as follows. Let x_{0} := x; for k = 0, 1, …, B(α̃_{i}), where B(α̃_{i}) is given below, let

x_{k+1} := x_{k} + ρ(x_{k}, α̃)\frac{g_{x_{k}}}{∥g_{x_{k}}∥_{2}} ,

where ρ(x_{k}, α̃) is defined below, and

g_{x_{k}} ≡ −A^{T} sgn(Ax_{k} − b),

as defined in the introduction.

A minor iteration can quit before k = B(α̃_{i}), if

F (x_{k})/F (x_{0}) = ∥Ax_{k} − b∥_{1}/∥Ax_{0} − b∥_{1} ≤ (1 + 1/α̃)/2.

Say that a minor iteration that quits early for this reason is successful, and let the current estimate x := x_{k}. Otherwise the minor iteration is unsuccessful, and the final value of x_{k} is not used.

This completes the description of the algorithm, except for the definitions of α̃_{i}, B(α̃) and ρ(x_{k}, α̃). These are

α̃_{i} ≡ \begin{cases}2 & i = 0 \\ \frac{α̃_{i−1}}{4} (1 + \sqrt{1 + 8/α̃_{i−1}}) & i > 0. \end{cases}

and

B(α̃) ≡ \frac{ln[2μ\sqrt{d}(α̃ + 1)/(α̃ − 1)]}{κ((α̃ + 1)/2)^{2}/2},

where

κ(α̃) ≡ \frac{α̃ − 1}{μ\sqrt{d}(α̃ + 1)}.

Finally,

ρ(x, α̃) ≡ F (x)(1 − 1/α̃)/\sqrt{d}.

The motivations for these definitions should be clearer from the following analysis. The general idea is that g_{x} has positive dot product with x̂ − x, and the stepsize ρ is chosen so that progress toward x̂ is made, but the step doesn’t overshoot. The values α̃_{i} are intended as lower bounds on the approximation ratio F (x)/F (x̂). If indeed α̃_{i} ≤ F (x)/F (x̂) then, as will be shown, the minor iteration will be successful, and the approximation ratio will be reduced. The failure of a minor iteration amounts to a proof that α̃_{i} > F (x)/F (x̂), and so the lower bound α̃_{i} can be reduced.

The algorithm and analysis are less elegant than they might be, because for a given step, it hasn’t been proven that the objective function value is reduced, only that the ℓ_{2} distance to an optimum vector is reduced, and that reduction in distance is not directly known. Hence certain conditions must be inferred from the success or failure of a sufficient number of steps, that is, a minor iteration.

3.2 Analysis of the subgradient iteration

First, a theorem regarding the gradient of the ℓ_{1} regression objective function. Some notation will be needed: for vector c, again, let sgn c denote a vector with the same dimensions as c, but with a +1 entry when where c has a positive entry, a −1 when c has a negative entry, and zero where c has a negative entry. With this definition, we can write ∥Ax − b∥_{1} as sgn(Ax − b)^{T} (Ax − b).

When the gradient of ∥Ax − b∥_{1} exists, it is A^{T} sgn(Ax − b): suppose x is such that for any d-vector δ that is short enough, sgn(Ax − b) = sgn(A(x + δ) − b). Then

∥A(x + δ) − b∥_{1}
= sgn(A(x + δ) − b)^{T} (A(x + δ) − b)
= sgn(Ax − b)^{T} (A(x + δ) − b)
= sgn(Ax − b)^{T} (Ax − b) + sgn(Ax − b)^{T} Aδ
= ∥Ax − b∥_{1} + δ^{T} A^{T} sgn(Ax − b),

so the Taylor expansion of

∥A(x + δ) − b∥_{1}

is the above expression, and A^{T} sgn(Ax − b) is the gradient at x.

Indeed, A^{T} sgn(Ax − b) is a subgradient of ∥Ax − b∥_{1}, as shown in the next lemma.

Lemma 3.1 For any d-vectors x and y,

∥Ax − b∥_{1} − ∥Ay − b∥_{1} ≤ (x − y)^{T} A^{T} sgn(Ax − b).

That is, A^{T} sgn(Ax − b) is a subgradient.

Proof: We have

∥Ax − b∥_{1} − ∥Ay − b∥_{1}
= sgn(Ax − b)^{T} (Ax − b) − sgn(Ay − b)^{T} (Ay − b)
≤ sgn(Ax − b)^{T} (Ax − b) − sgn(Ax − b)^{T} (Ay − b)
= sgn(Ax − b)^{T} A[x − y]
= (x − y)^{T} A^{T} sgn(Ax − b),

Next we show that the gradient points toward an optimum vector x̂. As before, we will use the notation F (x) ≡ ∥Ax − b∥_{1}, and g_{x} ≡ −A^{T} sgn(Ax − b).

Lemma 3.2 Given an ℓ_{1} μ-conditioned matrix A and d-vector x, let

α(x, μ) ≡ F (x)/F (x̂),  (2)

and writing simply α when x and μ are understood. Let θ be the angle between g_{x} and x̂−x, so

cos θ ≡ \frac{g_{x}^{T} (x̂ − x)}{∥g_{x}∥_{2}∥x̂ − x∥_{2}} .  (3)

Then

∥x̂ − x∥_{2} cos θ ≥ ρ(x, α) ≡ F (x)(1 − 1/α)/\sqrt{d}  (4)

and

cos θ ≥ κ(α) ≡ \frac{α − 1}{μ\sqrt{d}(α + 1)}.  (5)

Proof: At x, we have

g^{T} (x̂ − x) ≥ F (x) − F (x̂) = F (x)(1 − 1/α),

by the previous lemma and definition of α, and omitting the subscript from g. Also ∥g∥_{2} ≤ \sqrt{d} since the column sums of A have been set to 1. From these considerations,

∥x̂ − x∥_{2} cos θ = \frac{g_{x}^{T} (x̂ − x)}{∥g_{x}∥_{2}} ≥ (F (x) − F (x̂))/\sqrt{d}
= F (x)(1 − 1/α)/\sqrt{d},

the first conclusion of the lemma.

From the triangle inequality, we have

∥Ax̂ − Ax∥_{1} ≤ ∥Ax − b∥_{1} + ∥Ax̂ − b∥_{1} = F (x) + F (x̂),

and so, from the μ-conditioning of A,

\begin{aligned}∥x̂ − x∥_{2} & ≤ ∥x̂ − x∥_{1} \\ & ≤ μ∥Ax̂ − Ax∥_{1} \\ & ≤ μ(F (x) + F (x̂)) \\ & ≤ μ(F (x) + F (x)/α) \\ & = μF (x)(1 + 1/α), \end{aligned}  (6)

and so

cos θ ≥ \frac{F (x)(1 − 1/α)/\sqrt{d}}{∥x̂ − x∥_{2}}

≥ \frac{F (x)(1 − 1/α)/\sqrt{d}}{μF (x)(1 + 1/α)}

= \frac{α − 1}{μ\sqrt{d}(α + 1)}

the other conclusion of the lemma.

We need a quantitative statement about how far to go in the direction g_{x}.

Lemma 3.3 With conditions as in the previous lemma,

∥x + ρ(x, α)\frac{g_{x}}{∥g_{x}∥_{2}} − x̂∥_{2} ≤ (1 − κ(α)^{2}/2)∥x − x̂∥_{2},

where ρ(x, α) and κ are defined above (4,5).

Proof: For convenience in this proof, put x̂ at the origin. We have

∥x + ρ\frac{g_{x}}{∥g_{x}∥_{2}} ∥^{2} = ∥x∥_{2}^{2} + ρ^{2} − 2∥x∥ρ cos θ

and so, writing κ for κ(α) ≡ ρ/∥x∥_{2},

∥x + ρ\frac{g_{x}}{∥g_{x}∥_{2}} ∥^{2}/∥x∥_{2}^{2} = 1 + \frac{ρ^{2}}{∥x∥_{2}^{2}} − 2\frac{ρ}{∥x∥_{2}} cos θ
= 1 + κ^{2} − 2κ cos θ
≤ 1 − κ^{2}

using (4) and (5) of the last lemma. The lemma follows, taking square roots of both sides and using \sqrt{1 − w} ≤ 1 − w/2.

A difficulty in using this result is that generally, F (x̂) is unknown, which implies α is unknown. However, the following lemma implies that we can tell when an estimate α̃ of α is too big.

Lemma 3.4 With the conditions of the previous lemma, suppose α̃ ≤ α ≡ α(x). Let x_{0} ≡ x, and let

x_{k+1} ≡ x_{k} + ρ(x_{k}, α̃)\frac{g_{x_{k}}}{∥g_{x_{k}}∥_{2}} ,

for k ≥ 0. There is a value

B(α̃) ≡ \frac{ln[2μ\sqrt{d}(α̃ + 1)/(α̃ − 1)]}{κ((α̃ + 1)/2)^{2}/2}  (7)

such that for some k ≤ B(α̃),

F (x_{k}) ≤ F (x)(1 + 1/α̃)/2.

That is, if α̃ ≤ α ≡ α(x_{0}), a minor iteration succeeds.

Proof: Suppose the claim is false, so that during the B(α̃) iterations, F (x_{k}) remains above the given bound, implying

F (x_{k}) > F (x_{0})(1 + 1/α̃)/2  (8)
= F (x̂)α(x_{0})(1 + 1/α̃)/2
≥ F (x̂)(α(x_{0}) + 1)/2  (9)

and so α(x_{k}) ≥ (α + 1)/2 for k = 1 … B(α̃), where we abbreviate α(x_{0}) as α. We will show that this assumption implies a contradiction.

Based on this condition, the previous lemma implies that after B(α̃) iterations,

∥x_{k} − x̂∥_{2}

≤ ∥x_{0} − x̂∥_{2} ∏_{1≤k≤B(α̃)}(1 − κ(α(x_{k}))^{2}/2)

≤ (1 − κ((α + 1)/2)^{2}/2)^{B(α̃)}∥x_{0} − x̂∥_{2},  (10)

using the fact that κ() is increasing. Using the triangle inequality and the ℓ_{1} conditioning of A,

F (x_{k}) = ∥Ax_{k} − b∥_{1}
≤ ∥Ax̂ − b∥_{1} + ∥Ax_{k} − Ax̂∥_{1}
≤ F (x̂) + ∥x_{k} − x̂∥_{1}
≤ F (x̂) + ∥x_{k} − x̂∥_{2}\sqrt{d}.

That is, using (10) and (6), and the facts that (α − 1)/(α + 1) and so κ(α) are increasing functions of α,

F (x_{k}) − F (x̂)

≤ ∥x_{k} − x̂∥_{2}\sqrt{d}

≤ (1 − κ((α + 1)/2)^{2}/2)^{B(α̃)}∥x_{0} − x̂∥_{2}\sqrt{d}

≤ exp(−B(α̃)κ((α + 1)/2)^{2}/2)μF (x)(1 + 1/α)\sqrt{d}

≤ \frac{μF (x̂)α(1 + 1/α)\sqrt{d}}{2μ\sqrt{d}(α̃ + 1)/(α̃ − 1)}

≤ F (x̂)(α − 1)/2,

which implies F (x_{k}) ≤ F (x̂)(α + 1)/2, contradicting assumption (9).

The following lemma gives the properties of the sequence of α̃_{i} that will be used in the analysis of the subgradient iteration algorithm.

Lemma 3.5 With α̃_{i} as defined,

α̃_{i+1} = α̃_{i}(1 + 1/α̃_{i+1})/2,

and α̃_{i} ≤ 1 + (3/4)^{i}.

Proof: The first claim requires only the quadratic formula. The second can be proven with induction: it holds for i = 0, and for i > 0, let β denote (3/4)^{i}. Then, using β ≤ 1 and \sqrt{1 + w} ≤ 1 + w/2,

4α̃_{i+1}
= (1 + β)(1 + \sqrt{1 + 8/(1 + β}))
= 1 + β + \sqrt{(1 + β)^{2} + 8(1 + β)}
= 1 + β + 3\sqrt{1 + 10β/9 + β^{2}/9}
≤ 1 + β + 3\sqrt{1 + 11β/9}
≤ 1 + β + 3(1 + 11β/18)
= 4 + 17β/6
≤ 4 + 3β,

and so α̃_{i+1} ≤ 1 + 3β/4 = 1 + (3/4)^{i+1}, and the claim follows inductively.

Lemma 3.6 For μ ≤ d\sqrt{d} and α̃ = 1 + ϵ,

B(α̃) = B(1 + ϵ) = O(d^{4} log(d/ϵ))/ϵ^{2}

as ϵ → 0.

Proof: From (7) and the definition (5) of κ(), we have

B(α̃) ≡ \frac{ln[2μ\sqrt{d}(α̃ + 1)/(α̃ − 1)]}{κ((α̃ + 1)/2)^{2}/2}

≤ \frac{2μ^{2}d(α̃ + 3)^{2} ln[μd\sqrt{d}(α̃ + 1)/(α̃ − 1)]}{(α̃ − 1)^{2}},

from which the lemma follows.

Theorem 3.7 A d-vector x with ℓ_{1} regression value within ϵ of optimal can be found within

O(nd^{5}(log n log d + log(d/ϵ)/ϵ^{2}))

time.

Proof: The preliminary orthogonalization can be done with O(d^{2}) vector operations on b and the column vectors of A, each requiring O(n) time, for O(nd^{2}) time overall.

The rounding procedure, yielding the conditioned version of A, needs O(nd^{5} log n), from Lemma 2.2.

The algorithm requires O(nd) time for each of at most B(α̃_{i}) steps of the i’th minor iteration. By the previous lemma, B(2) = O(d^{4} log d). When the minor iteration with α̃_{i} = 2 is successful, the stopping condition implies that the approximation ratio F (x)/F (x̂) has decreased by a factor of (1 + 1/2)/2 = 3/4. A minor iteration is unsuccessful only when α̃_{i} > α(x), by Lemma 3.4, and so when a minor iteration with α̃_{i} = 2 is unsuccessful, it must be that the approximation ratio α(x) < 2. By beginning with b orthogonal to the columns of A, the initial approximation ratio is at most \sqrt{n}, since

∥b∥_{1} ≤ ∥b∥_{2}\sqrt{n} ≤ ∥b̂∥_{2}\sqrt{n} ≤ ∥b̂∥_{1}\sqrt{n},

where again, b̂ ≡ Ax̂ − b. Therefore, the work for all the minor iterations with α̃_{i} = 2 is

O(nd)B(2) log n = O(nd^{5} log d log n).

For the minor iterations with i > 0 (and so α̃_{i} < 2), note that after the minor iteration is done using α̃_{i}, it holds inductively that α ≤ α̃_{i}: either the iteration fails, which can only happen if α ≤ α̃_{i}, or it succeeds, which implies that the final iterate x_{k} satisfies

F (x_{k}) ≤ F (x_{0})(1 + α̃_{i})/2
≤ F (x_{0})α̃_{i−1}(1 + α̃_{i})/2
= F (x_{0})α̃_{i},

using the condition for success, the inductive assumption, and Lemma 3.5.

Also from that lemma, α̃_{i} ≤ 1 + (3/4)^{i}, and so O(log(1/ϵ)) minor iterations are needed. By the previous lemma, O(d^{4} log(d/(α̃_{i} − 1))/(α̃_{i} − 1)^{2}) steps are needed for each minor iteration, and from that bound the number of steps overall is within a constant factor of the number of steps in the last minor iteration, so the total work for the minor iterations with α̃_{i} < 2 is

O(nd)d^{4} log(d/ϵ)/ϵ^{2}.

Putting this bound together with the bound for α̃_{i} = 2 yields the theorem.

4 A Sampling Algorithm

First, some probabilistic bounds, and then the sampling construction.

4.1 Tail Estimates

Theorem 4.1 Let X_{i}, i = 1 … n, be independent random variables, with bounded EX_{i}^{2} and X_{i} ≥ 0 for all i. For S ≡ ∑_{i} X_{i} and t ≥ 0,

log Prob\{S ≤ ES − t\} ≤ \frac{−t^{2}}{2 ∑_{i} EX_{i}^{2}} .

If there is also some M with X_{i} ≤ EX_{i} + M for all i, then

log Prob\{S ≥ ES + t\} ≤ \frac{−t^{2}}{2tM/3 + 2 ∑_{i} EX_{i}^{2} − [EX_{i}]^{2}}

Proof: The bound on the lower tail seems to be due to Maurer[Mau03] and independently McAllester and Ortiz[MO03]; the former gave a simple proof. The upper tail bound is due to Bernstein.[Ber46]

4.2 The Sampling Algorithm

The first step of the overall algorithm is to apply the rounding algorithm of Section 2 to d\sqrt{d}-condition A. It also is necessary to apply the subgradient algorithm of Section 3 with ϵ = 1 to b, so that ∥b∥_{1} ≤ 2∥Ax̂ − b∥_{1}. (That is, having obtained an approximate solution x̃, replace b by b − Ax̃, and make a change of variable in x to account for the different b.) It will be convenient to then scale b so that ∥b∥_{1} = d.

Following these preprocessing steps, the sampling algorithm is as follows: let

f_{i} ≡ |b_{i}| + ∥a_{i·}∥_{1},

and let

p_{i} ≡ min{1, rf_{i}/2d},

where r is an integer parameter, noting that ∑_{i} f_{i} = 2d.

Let Y and Z be diagonal n × n matrices where

Y_{ii} = \begin{cases}1 & with probability p_{i} \\ 0 & with probability 1 − p_{i}, \end{cases}

and Z_{ii} ≡ Y_{ii}/p_{i}. Now solve the sample problem

min_{x}∥Z(Ax − b)∥_{1}

approximately, and return its output as an approximation to the solution of the original problem. This is the whole of the algorithm.

By construction, we have the expected number of nonzero Y_{ii} to be no more than

E ∑_{i} Y_{ii} = ∑_{i} EY_{ii} = ∑_{i} p_{i} ≤ ∑_{i} rf_{i}/2d = r,

and since EZ_{ii} = EY_{ii}/p_{i} = p_{i}/p_{i} = 1,

E∥Z(Ax − b)∥_{1} = E ∑_{i} |Z_{ii}(b_{i} − a_{i·}x)|

= ∑_{i} EZ_{ii}|b_{i} − a_{i·}x|

= ∥Ax − b∥_{1}.

So the sample problem gives an unbiased estimate of the objective function value of x, for any x. Of course, this does not show that the sample problem has an objective function value near that of x̂, since we don’t know how concentrated these values are around their expectations. However, we can apply the tail bounds, Theorem 4.1, to show that with high probability, ∥Z(Ax̂ − b)∥_{1} is not too large, and for any given x, ∥Z(Ax − b)∥_{1} is not too small. This will require a bound on ∑_{i} E[(Z_{ii}|b_{i} − a_{i·}x|)^{2}], before which, a lemma.

Lemma 4.2 For d\sqrt{d}-conditioned n × d matrix A, and d-vector b with ∥b∥_{1} ≤ 2 min_{x}∥b − Ax∥_{1}, we have

∥x∥_{1} ≤ d\sqrt{d}∥Ax − b∥_{1}.

Proof: Using Lemma 2.3 and the hypothesis,

∥x∥_{1} ≤ d\sqrt{d}∥Ax∥_{1}
≤ d\sqrt{d} (∥b∥_{1} + ∥Ax − b∥_{1})
≤ d\sqrt{d} (2∥Ax̂ − b∥_{1} + ∥Ax − b∥_{1})
≤ 3d\sqrt{d}∥Ax − b∥_{1}.

Lemma 4.3 Under the conditions of the sampling algorithm,

∑_{i} E[(Z_{ii}|b_{i} − a_{i·}x|)^{2}] ≤ \frac{6d^{2}\sqrt{d}}{r}∥Ax − b∥_{1}^{2}.

Proof: Using the definitions, the previous lemma, standard facts, and ignoring those i for which p_{i} = 1 < rf_{i}/2d,

∑_{i}E[(Z_{ii}|b_{i} − a_{i·}x|)^{2}]

= ∑_{i} \frac{|b_{i} − a_{i·}x|^{2}}{p_{i}^{2}} E[Y_{ii}^{2}]

= ∑_{i} \frac{|b_{i} − a_{i·}x|^{2}}{p_{i}}

= ∑_{i} \frac{|b_{i} − a_{i·}x|^{2}}{(|b_{i}| + ∥a_{i·}∥_{1})r/2d}

≤ ∑_{i} \frac{|b_{i} − a_{i·}x|(|b_{i}| + ∥a_{i·}∥_{1}∥x∥_{∞})}{(|b_{i}| + ∥a_{i·}∥_{1})r/2d}

≤ ∑_{i} \frac{|b_{i} − a_{i·}x|(|b_{i}| + ∥a_{i·}∥_{1}) max\{1, ∥x∥_{∞}\}}{(|b_{i}| + ∥a_{i·}∥_{1})r/2d}

= \frac{2d}{r}∥Ax − b∥_{1} max\{1, ∥x∥_{∞}\}

≤ \frac{2d}{r}∥Ax − b∥_{1}(3d\sqrt{d}∥Ax − b∥_{1})

= \frac{6d^{2}\sqrt{d}}{r}∥Ax − b∥_{1}^{2},

and the lemma follows.

Lemma 4.4 For ϵ ≤ 1/2 and a given d-vector x, the probability is at most

exp(−ϵ^{2}r/24d^{2}\sqrt{d})

that

∥Z(Ax − b)∥_{1} ≤ (1 + ϵ)∥b̂∥_{1}

when ∥Ax − b∥_{1} ≥ (1 + 2ϵ)∥b̂∥_{1}.

Proof: The lower tail bound of Theorem 4.1 can be applied, with random variables X_{i} = Z_{ii}|a_{i·}x − b_{i}|, sum S = ∥Z(Ax − b)∥_{1}, and

t = ES − (1 + ϵ)∥b̂∥_{1} = ∥Ax − b∥_{1} − (1 + ϵ)∥b̂∥_{1}.

Using also the previous lemma and ϵ ≤ 1/2, we have

log Prob{∥Z(Ax − b)∥_{1} ≤ (1 + ϵ)∥b̂∥_{1}}

≤ \frac{−(∥Ax − b∥_{1} − ∥b̂∥_{1}(1 + ϵ))^{2}}{6d^{2}\sqrt{d}∥Ax − b∥_{1}^{2}/r}

≤ \frac{−r}{6d^{2}} (1 − \frac{(1 + ϵ)∥b̂∥_{1}}{(1 + 2ϵ)∥b̂∥_{1}} )^{2}

≤ \frac{−rϵ^{2}}{24d^{2}\sqrt{d}},

as claimed.

To apply the upper tail bound of Theorem 4.1, we need an upper bound on each coordinate of b̂.

Lemma 4.5 For all i, Z_{ii}|b̂_{i}| ≤ 6d^{2}\sqrt{d}∥b̂∥_{1}/r.

Proof: From Lemma 4.2, we have

∥x̂∥_{1} ≤ 3d\sqrt{d}∥b̂∥_{1}/r.

Therefore

Z_{ii}|b̂_{i}| = |b̂_{i}|/p_{i} = |b̂_{i}|2d/f_{i}r

≤ \frac{2d}{r} \frac{|b_{i}| + ∥a_{i·}∥_{1}∥x̂∥_{∞}}{|b_{i}| + ∥a_{i·}∥_{1}}

≤ \frac{2d}{r} \frac{|b_{i}| + ∥a_{i·}∥_{1}3d\sqrt{d}∥b̂∥_{1}}{|b_{i}| + ∥a_{i·}∥_{1}}

≤ \frac{6d^{2}\sqrt{d}∥b̂∥_{1}}{r} .

Lemma 4.6 For ϵ ≤ 1, the probability is at most exp(−ϵ^{2}r/16d^{2}\sqrt{d}) that

∥Zb̂∥_{1} ≥ (1 + ϵ)∥b̂∥_{1}.

Proof: The second part of Theorem 4.1 can be applied, with X_{i} = Z_{ii}|b̂_{i}|, S = ∥Zb̂∥_{1}, t = ϵ∥b̂∥_{1}, and from the previous lemma, M = 6d^{2}\sqrt{d}∥b̂∥_{1}/r. Consequently,

log Prob{∥Zb̂∥_{1} > ∥b̂∥_{1} + t}

≤ \frac{−t^{2}}{2tM/3 + 2 ∑_{i} EX_{i}^{2} − [EX_{i}]^{2}}

≤ \frac{−(ϵ∥b̂∥_{1})^{2}}{2(ϵ∥b̂∥_{1})(6d^{2}\sqrt{d}∥b̂∥_{1}/r)/3 + 2^{\frac{6d^{2}\sqrt{d}}{r}}∥b̂∥_{1}^{2}}

= \frac{−ϵ^{2}}{4ϵd^{2}\sqrt{d}/r + 12d^{2}\sqrt{d}/r}

≤ \frac{−ϵ^{2}r}{16d^{2}\sqrt{d}}

for ϵ ≤ 1.

These two lemmas can be combined to prove the following theorem.

Theorem 4.7 For any n × d matrix A and n-vector b, and given ϵ > 0 and γ > 0, there is a value R where

R ≡ \frac{288d^{3}\sqrt{d}}{ϵ^{2}} ln(d/ϵγ)

such that if r ≥ R, then with probability at least 1 − γ, the sample problem with parameter r has optimum solution x_{Z} with

∥Ax_{Z} − b∥_{1} ≤ (1 + ϵ)∥b̂∥_{1},

where ∥b̂∥_{1} = ∥Ax̂ − b∥_{1} is the minimum value of ∥Ax − b∥_{1}.

Proof: We will show that with high probability, x̂ yields an objective function value for the sample problem that is not too big, and that for all vectors, their sample problem objective function value is not too small. Here “too big” means larger than (1 + ϵ/2)∥b̂∥_{1}, and “too small” means smaller than (1+ϵ/2)∥b̂∥_{1}, when the vector x gives objective function value greater than 1 + ϵ times the optimum.

For the latter, it is enough to show that any possible optimum vector x_{Z} gives a value that is not too small. Such optimum vectors come from a finite set of vectors: since ℓ_{1} regression is a linear programming problem, its optima are at vertices of its feasible polytope; here those vertices correspond to hyperplanes determined by d points of the n input points S. There are \binom{n}{d} such hyperplanes. For each such hyperplane, let p_{i_{1}}, p_{i_{2}}, …, p_{i_{d}} denote the points in S that determine it. Let T denote the corresponding set of indices, and let x_{T} denote the corresponding solution.

The probability that x_{T} is a bad solution to the sample problem is bounded by the probability that every index i ∈ T has Y_{ii} = 1, and that the resulting objective function value is too small. The former probability is the product of the probabilities ∏_{i∈T} p_{i}. The latter probability is no more than

exp(−ϵ^{2}r/96d^{2}\sqrt{d})

by Lemma 4.4 (using ϵ/2 as the ϵ of the lemma). Let \binom{N}{d} denote the set of all sets of size d of integers between 1 and n. Using these facts and the bound on r of the theorem statement, the probability that some x_{T}, for T ∈ \binom{N}{d}, has a value that is too small is no more than

∑_{T ∈\binom{N}{d}} exp(−ϵ^{2}R/96d^{2}\sqrt{d}) ∏_{i∈T} p_{i}

≤ exp(−ϵ^{2}R/96d^{2}\sqrt{d}) (∑_{i} p_{i})^{d}

= exp(−ϵ^{2}R/96d^{2}\sqrt{d})R^{d}

≤ γ/2.

Applying Lemma 4.6, the probability that x̂ yields an objective function value for the sample problem that is too large is at most

exp(−ϵ^{2}R/64d^{2}\sqrt{d}) ≤ γ/2,

and so the probability of failure is at most γ, as claimed.

Theorem 4.8 For any n × d matrix A and n-vector b, and given ϵ > 0 and γ > 0, in

O(nd^{5} log n log d) + O(\frac{d^{8}\sqrt{d}}{ϵ^{4}}(log d)^{3} ln(1/ϵ)^{2} ln(1/γ))

time, a vector x^{∗} can be found such that, with probability at least 1 − γ, ∥b − Ax^{∗}∥_{1} ≤ (1 + ϵ) min_{x}∥b − Ax∥_{1}.

Proof: Applying Theorem 3.7 with ϵ = 2, the time needed for preprocessing of A and b is O(nd^{5} log n log d).

The sample size R needed to obtain a sampling approximation less than ϵ/3 is given by the last result; applying Theorem 3.7 again, the time to find a solution with 1 + ϵ/3 of optimum for the sampling problem is

O(Rd^{5}(log R log d + log(d/ϵ)/ϵ^{2}))

= O(\frac{d^{3}\sqrt{d}}{ϵ^{2}}(log d) ln(d/ϵγ)(log(d/ϵ)/ϵ^{2}))

= O(\frac{d^{8}\sqrt{d}}{ϵ^{4}}(log d)^{3} ln(1/ϵ)^{2} ln(1/γ)),

as claimed.

5 Concluding Remarks

The log n term in the runtime might be removable by using a sampling approach to estimate the subgradient, somewhat like the work of Nedić and Bertsekas.[NB01]

Given that the ellipsoid method was used, one might ask why that method could not be simply applied to solve the problem. Perhaps it can be, but it is not obvious that a bit-complexity dependence can be avoided; moreover, the subgradient descent scheme may be of independent interest. It seems likely that the conditioning of A can be done using the regression procedure itself, by making each column “ℓ_{1} independent” of the others, and then the ellipsoid method is not needed.

References

[AHPVar] P. Agarwal, S. Har-Peled, and K. R. Varadarajan. Approximating extent measures of points. J. ACM, to appear.

[BC02] Mihai Bădoiu and K. L. Clarkson. Optimal core-sets for balls. Manuscript,

          available via http://cm.bell-labs.com/who/clarkson/, 2002.

[Ber46] S. Bernstein. The Theory of Probabilities. Gastehizdat Publishing House, 1946.

[Ber01] D. Bertsekas. Convex Analysis and Optimization. Athena Scientific, 2001.

[BHPI02] M. Bădoiu, S. Har-Peled, and P. Indyk. Approximate clustering via core-sets. In Proc. 34th Symp. Theory of Comp., 2002.

[Cla95] K. L. Clarkson. Las Vegas algorithms for linear and integer programming when the dimension is small. Journal of the ACM, 42(2):488–499, 1995.

[DK01] P. Drinneas and R. Kannan. Fast Monte-Carlo algorithms for approximate matrix multiplication. In EEE Symposium on Foundations of Computer Science, 2001.

[DKM] P. Drineas, R. Kannan, and M. W. Mahoney. Fast monte carlo algorithms for matrices I: Approximating matrix multiplication. Manuscript. Available via

          http://cs-www.cs.yale.edu/homes/mmahoney/.

[HPV02] S. Har-Peled and K. R. Varadarajan. Projective clustering in high dimensions using core-sets. In Symp. Comp. Geometry, 2002.

[Lov86] L. Lovász. Algorithmic Theory of Numbers, Graphs, and Convexity. SIAM, 1986.

[Mau03] A. Maurer. A bound on the deviation probability for sums of non-negative random variables. J. Inequalities in Pure and Applied Mathematics, 4, 2003.

[MO03] D. McAllester and L. Ortiz. Concentration inequalities for the missing mass and for histogram rule error. J. Machine Learning Research, 4:895–911, 2003.

[MT93] N. Megiddo and A. Tamir. Linear time algorithms for some separable quadratic programming problems. Operations Research Letters, 13:203–211, March 1993.

[NB01] A. Nedić and D. Bertsekas. Incremental subgradient methods for nondifferentiable optimization. SIAM J. Opt., 12(1):109–138, 2001.

[Sch86] A. Schrijver. Theory of Linear and Integer Programming. Wiley, New York, 1986.

[Sho85] N. Z. Shor. Minimization Methods for Non-differentiable Functions. Springer, 1985.

[YKII88] P. Yamamoto, K. Kato, K. Imai, and H. Imai. Algorithms for vertical and orthogonal L1 linear approximation of points. In Proceedings of the fourth annual symposium on Computational geometry, pages 352–361. ACM Press, 1988.