Near-optimal Coresets For Least-Squares Regression

Introduction

Linear regression is an important technique in data analysis . Research in the area ranges from numerical techniques to robustness of the prediction error to noise (e.g., using feature selection ). We ask whether it is possible to efficiently identify a small subset of the data that contains all the essential information of a learning problem. Such a subset is called a “coreset”. We show that the answer is yes, for linear regression. Such a coreset is analogous to the support vectors in support vector machines . Such coresets contain the meaningful or important points in the data and can be used to find good approximate solutions to the full problem by solving a (much) smaller problem. When the constraints are complex (e.g., non-convex constraints), solving a much smaller regression problem could be a significant saving .

We present coreset constructions for constrained regression (both simple and multiple response), as well as lower bounds for the size of coresets that achieve certain accuracy. In addition to potential computational savings, a coreset identifies the important core of a machine learning problem and is of considerable interest in applications with huge data where incremental approaches are necessary (for example chunking) and applications where the data is distributed and bandwith is costly (hence communicating only the essential data is imperative ).

Our first contribution is a deterministic, polynomial-time algorithm for constructing a coreset for arbitrarily constrained linear regression. Let $k$ be the “effective dimension” of the data (the rank of the data matrix) and let $\epsilon>0$ be the desired accuracy parameter. Our algorithm constructs a coreset of size $O\left(k/\epsilon^{2}\right)$, which achieves a $\left(1+\epsilon\right)$-relative error performance guarantee. In other words, solving the regression problem on the coreset results in a solution which fits all the data with an error which is at most $\left(1+\epsilon\right)$ worse than the best possible fit to all the data. We extend our results to the setting of multiple response regression using more sophisticated techniques. Our proofs are based on two sparsification tools from linear algebra , which may be of general interest to the machine learning community, and we discuss these in some detail.

Assume the usual setting with $n$ data points $({\mathbf{z}}_{1},y_{1}),\ldots,({\mathbf{z}}_{n},y_{n})$; ${\mathbf{z}}_{i}\in\R^{d}$ are feature vectors (which could have been obtained by applying a non-linear feature transform to raw data) and $y_{i}\in\R$ are targets (responses). The linear regression problem asks to determine a vector ${\mathbf{x}}_{opt}\in{\cal D}\subseteq\R^{d}$ that minimizes

over ${\mathbf{x}}\in{\cal D}$, where $w_{i}\in\R$ are positive weights. So, ${\cal E}({\mathbf{x}}_{opt})\leq{\cal E}({\mathbf{x}})$, for all ${\mathbf{x}}\in{\cal D}$. The domain ${\cal D}$ represents the constraints on the solution, e.g., in non-negative least squares (NNLS) , ${\cal D}=\R^{d}_{+}$, the nonnegative orthant. Our results hold for arbitrary ${\cal D}$.

A coreset of size $r<n$ is a subset of the data points, $({\mathbf{z}}_{i_{1}},y_{i_{1}}),\ldots,({\mathbf{z}}_{i_{r}},y_{i_{r}})$. The coreset regression problem considers the squared error on the coreset with a (possibly) different set of weights $s_{j}>0$,

The algorithm which constructs the coreset should also provide the weights ${\mathbf{s}}_{j}$. For the remainder of the paper, we switch to an equivalent matrix formulation of the problem. (See Appendix for linear algebra background.)

Let ${\mathbf{A}}\in\R^{n\times d}$ be the data matrix whose rows are the weighted data points $\sqrt{w_{i}}{\mathbf{z}}_{i}^{\textsc{T}}$; and let ${\mathbf{b}}\in\R^{n}$ be the similarly weighted target vector, $b_{i}=\sqrt{w_{i}}y_{i}$, where for $i=1,...,n,$ $b_{i}$ denotes the $i$th element of ${\mathbf{b}}\in\R^{n}$. The effective dimension of the data can be measured by the rank of ${\mathbf{A}}$; let $k=\hbox{\rm rank}({\mathbf{A}})$. Our results hold for arbitrary $n>d$, however, in most applications, $n\gg d$ and $\hbox{\rm rank}({\mathbf{A}})\approx d$. We can rewrite the squared error as ${\cal E}({\mathbf{x}})=\mbox{}\|{\mathbf{A}}{\mathbf{x}}-{\mathbf{b}}\|_{2}^{2}$, so,

A coreset of size $r<n$ is a subset ${\mathbf{C}}\in\R^{r\times d}$ of the rows of ${\mathbf{A}}$ and the corresponding elements ${\mathbf{b}}_{c}\in\R^{r}$ of ${\mathbf{b}}$. Let ${\mathbf{D}}\in\R^{r\times r}$ be a positive diagonal matrix for the coreset regression (the weights $s_{j}$ of the coreset regression will depend on ${\mathbf{D}}$). The weighted squared error on the coreset is given by

We say that such a coreset is an $(1+\epsilon)$-coreset if the solution obtained by fitting the coreset data is almost optimal for all the data. Formally,

2 Our contributions

In this section, we discuss our main results for various formulations of linear regression (also summarized in Table 1). In the next section we present the relevant algorithms and proofs.

Our main result for constrained simple regression is Theorem 1, which describes a deterministic polynomial time algorithm that constructs a $(1+\epsilon)$-coreset of size $O\left(k/\epsilon^{2}\right)$. Prior to our work, the best result achieving comparable relative error performance guarantees is Theorem 1 of for constrained regression, and the work of for unconstrained regression. Both of these prior results construct coresets of size $O\left(k\log k/\epsilon^{2}\right)$ and they are randomized, so, with some probability, the fit on all the data can be arbitrarily bad (despite the coreset being a logarithmic factor larger). Our methods have comparable, low order polynomial running times and provide deterministic guarantees. The results in and were achieved using the matrix concentration results in . However, these concentration bounds break unless the coreset size is $\Omega\left(k\log k/\epsilon^{2}\right)$.

2.2 Multi-Objective Regression (Section 3.1)

An important variant of multiple response regression is the so-called multi-objective regression. Let

2.3 Arbitrarily-Constrained Multiple-Response Regression (Section 3.2)

Using the same approach, converting the problem to a single response regression, we construct a $(1+\epsilon)$-coreset for Frobenius-norm arbitrarily-constrained regression in Section 3.2. The coreset size in this case is $O\left(k\omega/\epsilon^{2}\right)$.

2.4 Unconstrained Multiple-Response Regression (Section 4)

In Section 4, we consider coresets for unconstrained multiple-response regression for both the spectral and Frobenius norms. The sizes of the coresets are smaller than the constrained case, and our main results are presented in Theorems 6 and 7. Theorem 6 presents a $(2+\epsilon)$-coreset of size $O((k+\omega)/\epsilon^{2})$ for spectral norm regression, while Theorem 7 presents a $(2+\epsilon)$-coreset of size $O(k/\epsilon^{2})$ for Frobenius norm regression.

2.5 Lower Bounds (Section 5)

Finally, in Section 5, we present lower bounds on coreset sizes. In the single response regression setting, we note that our algorithms need to look at the target vector ${\mathbf{b}}$. We show that this is unavoidable, by arguing that no ${\mathbf{b}}$-agnostic deterministic coreset construction algorithm can construct coresets which are small (Theorem 13). We also present similar results for ${\mathbf{b}}$-agnostic randomized coreset constructions (Theorem 14).

Then, we present lower bounds on the size of coresets for spectral and Frobenius norm multiple response regression that apply in the general, non ${\mathbf{b}}$-agnostic, setting (Theorems 15 and 16).

Constrained Linear Regression

We define constrained linear regression as follows: given ${\mathbf{A}}\in\R^{n\times d}$ of rank $k$, ${\mathbf{b}}\in\R^{n}$, and $\mathcal{D}\subseteq\R^{d}$, we seek ${\mathbf{x}}_{opt}\in{\cal D}$ for which $\mbox{}\|{\mathbf{A}}{\mathbf{x}}_{opt}-{\mathbf{b}}\|_{2}^{2}\leq\mbox{}\|{\mathbf{A}}{\mathbf{x}}-{\mathbf{b}}\|_{2}^{2}$, for all ${\mathbf{x}}\in{\cal D}$ (the domain $\mathcal{D}$ represents the constraints on ${\mathbf{x}}$ and can be arbitrary). To construct a coreset ${\mathbf{C}}\in\R^{r\times d}$ (i.e., ${\mathbf{C}}$ consists of $r$ rows of ${\mathbf{A}}$) and ${\mathbf{b}}_{c}\in\R^{r}$ (i.e., ${\mathbf{b}}_{c}$ consists of $r$ elements of ${\mathbf{b}}$), we introduce sampling and rescaling matrices ${\mathbf{S}}$ and ${\mathbf{D}}$ respectively. More specifically, we define the row-sampling matrix ${\mathbf{S}}\in\R^{r\times n}$ whose rows are basis vectors ${\mathbf{e}}_{i_{1}}^{\textsc{T}},\ldots,{\mathbf{e}}_{i_{r}}^{\textsc{T}}$. Our coreset ${\mathbf{C}}$ is now equal to ${\mathbf{C}}={\mathbf{S}}{\mathbf{A}}$; clearly, ${\mathbf{C}}$ is a matrix whose rows are the rows of ${\mathbf{A}}$ corresponding to indices $i_{1},\ldots,i_{r}$. Similarly, ${\mathbf{b}}_{c}={\mathbf{S}}{\mathbf{b}}$ contains the corresponding elements of the target vector. Next, let ${\mathbf{D}}\in\R^{r\times r}$ be a positive diagonal rescaling matrix and define the ${\mathbf{D}}$-weighted regression problem on the coreset as follows:

In the above, the operator ${\mathbf{D}}{\mathbf{S}}$ first samples and then rescales rows of ${\mathbf{A}}$ and ${\mathbf{b}}$. Theorem 1 is the main result in this section and presents a deterministic algorithm to select a coreset by constructing ${\mathbf{D}}$ and ${\mathbf{S}}$.

The running time of the proposed algorithm is $T\left({\mathbf{U}}_{\left[{\mathbf{A}},{\mathbf{b}}\right]}\right)+O\left(rnk^{2}\right)$, where $T\left({\mathbf{U}}_{\left[{\mathbf{A}},{\mathbf{b}}\right]}\right)$ is the time needed to compute the left singular vectors of the matrix $\left[{\mathbf{A}},{\mathbf{b}}\right]\in\R^{n\times(d+1)}$.

For any $0<\epsilon<1$, we can set $r=k/\epsilon^{2}$ to get an approximation ratio roughly equal to $1+4\epsilon$. This result considerably improves the result in , which needs $r=O(k\log k/\epsilon^{2})$ to achieve the same approximation ratio. Additionally, our bound is deterministic, whereas the bound in fails with constant probability. also requires an SVD computation in the first step, so its running time is comparable to ours.

In order to prove the above theorem, we need a linear algebraic sparsification result from , specifically Theorem 3.1 in , which we restate using our notation (we present the corresponding algorithm below).

We now discuss in more detail the sparsification algorithm of Lemma 2. We present the corresponding algorithm as Algorithm 6. Our notation deviates from the original in ; we employ our own presentation of the corresponding algorithm in . Algorithm 6 is a greedy technique that selects columns one at a time. To describe the algorithm in more detail, it is convenient to view the input matrix as a set of $n$ column vectors,

and let $L({\mathbf{u}},\delta_{L},{\mathbf{A}},{\textsc{l}})$ be defined as

and let $U({\mathbf{u}},\delta_{U},{\mathbf{A}},{\textsc{u}})$ be defined as

The running time of the algorithm is dominated by the search for an index $i_{\tau}$ satisfying

Constrained Multiple-Response Regression

Constrained multiple-response regression in the Frobenius norm can be reduced to simple regression. So, we can apply the results of the previous section to this setting.

Let ${\mathbf{A}}\in\R^{n\times d}$ and ${\mathbf{B}}\in\R^{n\times\omega}$, with $\omega\geq 1$. The objective of multi-objective regression is:

where $[{\mathbf{x}},\ldots,{\mathbf{x}}]\in\R^{d\times\omega}$ contains $\omega$ copies of ${\mathbf{x}}\in\mathcal{D}\subseteq\R^{d}$. Let ${\mathbf{b}}_{avg}={1\over\omega}{\mathbf{B}}\bm{1}_{\omega}$ (here $\bm{1}_{\omega}\in\R^{\omega}$ is a vector of all ones and thus ${\mathbf{b}}_{avg}\in\R^{n}$ is the average of the columns in ${\mathbf{B}}$). Recall that ${\mathbf{A}}\in\R^{n\times d}$, ${\mathbf{B}}\in\R^{n\times\omega}$, and let ${\mathbf{X}}=[{\mathbf{x}},\ldots,{\mathbf{x}}]\in\R^{d\times\omega}$.

The run time of the proposed algorithm is $T\left({\mathbf{U}}_{\left[{\mathbf{A}},{\mathbf{b}}_{avg}\right]}\right)+O\left(n\omega+rnk^{2}\right)$, where $T\left({\mathbf{U}}_{\left[{\mathbf{A}},{\mathbf{b}}_{avg}\right]}\right)$ is the time needed to compute the left singular vectors of the matrix $\left[{\mathbf{A}},{\mathbf{b}}_{avg}\right]\in\R^{n\times(d+1)}$.

We first construct ${\mathbf{D}}$ and ${\mathbf{S}}$ via Theorem 1 applied to ${\mathbf{A}}$ and ${\mathbf{b}}_{avg}$. The running time is $O\left(n\omega\right)$ (the time needed to compute ${\mathbf{b}}_{avg}$) plus the running time of Theorem 1. The result is immediate from the following derivation:

2 Arbitrarily-Constrained Multiple-Response Regression

Multi-objective regression is a special case of constrained multiple-response regression for which we can efficiently obtain the coresets. In the general case, the problem still reduces to simple regression, but the coresets are now larger. The objective of arbitrarily-constrained multiple-response regression is

Since $\R^{d\times\omega}$ is isomorphic to $\R^{d\omega}$, we can view ${\mathbf{X}}\in\R^{d\times\omega}$ as a “stretched out” vector $\hat{\mathbf{X}}\in\R^{d\omega}$; corresponding to the domain ${\cal D}$ is the domain $\hat{{\cal D}}\subseteq\R^{d\omega}$. Similarly, we can stretch out ${\mathbf{B}}\in\R^{n\times\omega}$ to $\hat{\mathbf{B}}\in\R^{n\omega}$. To complete the transformation to simple linear regression, we build a transformed block-diagonal data matrix $\hat{\mathbf{A}}$ from ${\mathbf{A}}$, by repeating $\omega$ copies of ${\mathbf{A}}$ along the diagonal:

So, for the approximation ratio to be $1+O(\epsilon)$, we set $r=O\left(k\omega/\epsilon^{2}\right)$. The running time would involve the time needed to compute the SVD of $[\hat{\mathbf{A}},\hat{\mathbf{B}}]$.

Notice that the coresets are large and somewhat costly to compute and they only work for the Frobenius norm. In the next section, using more sophisticated techniques, we will get smaller coresets for unconstrained regression in both the Frobenius and spectral norms.

Unconstrained Multiple-Response Regression

We can compute ${\mathbf{X}}_{opt}$ via the pseudoinverse of ${\mathbf{A}}$, namely ${\mathbf{X}}_{opt}={\mathbf{A}}^{\dagger}{\mathbf{B}}$. If ${\mathbf{S}}$ and ${\mathbf{D}}$ are sampling and rescaling matrices respectively, then the coreset regression problem is:

Given a matrix ${\mathbf{A}}\in\R^{n\times d}$ with rank $k$, a matrix ${\mathbf{B}}\in\R^{n\times\omega}$, and $r>k$, Algorithm 3 deterministically constructs matrices ${\mathbf{S}}\in\R^{r\times n}$ and ${\mathbf{D}}\in\R^{r\times r}$ such that the solution of the problem of Eqn. (7) satisfies:

The running time of the proposed algorithm is $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)+O\left(rn\left(k^{2}+\omega^{2}\right)\right)$, where $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)$ is the time needed to compute the left singular vectors of ${\mathbf{A}}$.

Since $r>k$, the approximation ratio is $2+O(\sqrt{{\omega}/{r}}+{\omega}/{r}+\sqrt{{k}/{r}})$. So, for $\epsilon>0$ and $r=O((\omega+k)/\epsilon^{2})$ the approximation ratio is $2+\epsilon$. For $r>\omega$, the approximation is $O(1)$, while for $r<\omega$, is asymptotic to $O\left(\omega/r\right)$. We will argue that this is nearly optimal by providing a matching lower bound in Theorem 15.

Given matrix ${\mathbf{A}}\in\R^{n\times d}$ of rank $k$, matrix ${\mathbf{B}}\in\R^{n\times\omega}$, and $r>k,$ Algorithm 4 deterministically constructs a sampling matrix ${\mathbf{S}}\in\R^{r\times n}$ and a rescaling matrix ${\mathbf{D}}\in\R^{r\times r}$ such that the solution of the problem of Eqn. (7) satisfies:

The running time of the proposed algorithm is $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)+O\left(rnk^{2}\right)$, where $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)$ is the time needed to compute the left singular vectors of ${\mathbf{A}}$.

The approximation ratio in the above theorem is $2+O(\sqrt{k/r})$. In Theorem 16, we will give a lower bound for the approximation ratio which is $1+\Omega(k/r)$. We conjecture that our lower bound can be achieved (deterministically), perhaps by a more sophisticated algorithm or analysis.

Finally, we note that the ${\mathbf{B}}$-agnostic randomized construction of achieves a $(1+\epsilon)$ approximation ratio using a significantly larger coreset, $r=O(k\log k/\epsilon^{2})$. Importantly, does not need any access to ${\mathbf{B}}$ in order to construct the coreset, whereas our approach constructs coresets by carefully choosing important data points with respect to the particular target response matrix ${\mathbf{B}}$. We will also discuss ${\mathbf{B}}$-agnostic algorithms in Section 4.2 (Theorem 12) and we will present matching lower bounds in Section 5.

We will make heavy use of facts from Section A in the Appendix. We start with a few simple lemmas.

Let ${\mathbf{E}}={\mathbf{A}}{\mathbf{X}}_{opt}-{\mathbf{B}}\in\R^{n\times\omega}$ be the regression residual. Then, $\hbox{\rm rank}({\mathbf{E}})\leq\min\{\omega,n-k\}$.

Using our notation, ${\mathbf{A}}{\mathbf{X}}_{opt}-{\mathbf{B}}=-\left({\mathbf{I}}_{n}-{\mathbf{U}}_{\mathbf{A}}{\mathbf{U}}_{\mathbf{A}}^{\textsc{T}}\right){\mathbf{B}}=-{\mathbf{U}}_{\mathbf{A}}^{\perp}\left({\mathbf{U}}_{\mathbf{A}}^{\perp}\right)^{\textsc{T}}{\mathbf{B}}$. To conclude notice that $\hbox{\rm rank}({\mathbf{X}}{\mathbf{Y}})\leq\min\{\hbox{\rm rank}({\mathbf{X}}),\hbox{\rm rank}({\mathbf{Y}})\}$ for any matrices ${\mathbf{X}}$ and ${\mathbf{Y}}$.

We now present our main tool for obtaining approximation guarantees for coreset regression.

To simplify notation, let ${\mathbf{W}}={\mathbf{D}}{\mathbf{S}}$. Using the SVD of ${\mathbf{A}}$, ${\mathbf{A}}={\mathbf{U}}_{\mathbf{A}}{\mathbf{\Sigma}}_{\mathbf{A}}{\mathbf{V}}_{\mathbf{A}}^{\textsc{T}}$, we get:

where the last equality follows from properties of the pseudo-inverse and the fact that ${\mathbf{W}}{\mathbf{U}}_{\mathbf{A}}$ is a full-rank matrix (see Lemma 18 in the Appendix). Using ${\mathbf{B}}=\left({\mathbf{U}}_{\mathbf{A}}{\mathbf{U}}_{\mathbf{A}}^{\textsc{T}}+{\mathbf{U}}_{\mathbf{A}}^{\perp}\left({\mathbf{U}}_{\mathbf{A}}^{\perp}\right)^{\textsc{T}}\right){\mathbf{B}}$, we obtain

$(a)$ follows from the assumption that the rank of ${\mathbf{W}}{\mathbf{U}}_{\mathbf{A}}$ is equal to $k$ and thus $\left({\mathbf{W}}{\mathbf{U}}_{\mathbf{A}}\right)^{\dagger}{\mathbf{W}}{\mathbf{U}}_{\mathbf{A}}={\mathbf{I}}_{k}$ and $(b)$ follows by matrix-Pythagoras (Lemma 17). To conclude, we use spectral submultiplicativity on the second term and the fact that ${\mathbf{U}}_{\mathbf{A}}^{\perp}\left({\mathbf{U}}_{\mathbf{A}}^{\perp}\right)^{\textsc{T}}{\mathbf{B}}=-({\mathbf{A}}{\mathbf{X}}_{opt}-{\mathbf{B}})$.

This lemma provides a framework for coreset construction: all we need are sampling and rescaling matrices ${\mathbf{S}}$ and ${\mathbf{D}}$, such that $\hbox{\rm rank}({\mathbf{D}}{\mathbf{S}}{\mathbf{U}}_{\mathbf{A}})=k$ and

is small. The final ingredients for the proofs of Theorems 6 and 7 are two matrix sparsification results that we present in the Appendix.

If ${\mathbf{\Psi}}={\mathbf{I}}_{n}$, the running time of the algorithm reduces to $T_{SVD}\left({\mathbf{Y}}\right)+O\left(rn\rho_{{\mathbf{Y}}}^{2}\right)$. We write $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleSpectralSampling\left({\mathbf{Y}},{\mathbf{\Psi}},r\right)$ to denote such a deterministic procedure.

If ${\mathbf{\Psi}}={\mathbf{I}}_{n}$, the running time of the algorithm reduces to $T_{SVD}\left({\mathbf{Y}}\right)+O\left(rn\rho_{{\mathbf{Y}}}^{2}\right)$. We write $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleFrobeniusSampling\left({\mathbf{Y}},{\mathbf{\Psi}},r\right)$ to denote such a deterministic procedure.

(of Theorem 6) Theorem 6 follows from Lemmas 9 and 10. First, compute the SVD of ${\mathbf{A}}$ to obtain ${\mathbf{U}}_{{\mathbf{A}}}\in\R^{n\times k}$, and let ${\mathbf{E}}={\mathbf{A}}{\mathbf{X}}_{opt}-{\mathbf{B}}={\mathbf{U}}_{{\mathbf{A}}}{\mathbf{U}}_{{\mathbf{A}}}^{\textsc{T}}{\mathbf{B}}-{\mathbf{B}}$. Next, run the algorithm of Lemma 10 to obtain $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleSpectralSampling\left({\mathbf{U}}_{{\mathbf{A}}},{\mathbf{E}},r\right)$. This algorithm runs in time $T_{SVD}\left({\mathbf{E}}\right)+O\left(rn\left(k^{2}+\rho_{{\mathbf{E}}}^{2}\right)\right)$, where $k$ is the rank of ${\mathbf{U}}_{\mathbf{A}}$ and ${\mathbf{A}}$. The total running time of the algorithm is $T({\mathbf{U}}_{\mathbf{A}})+T_{SVD}\left({\mathbf{E}}\right)+O\left(rn\left(k^{2}+\rho_{{\mathbf{E}}}^{2}\right)\right)=T\left({\mathbf{U}}_{\mathbf{A}}\right)+O\left(rn\left(k^{2}+\omega^{2}\right)\right)$.

Lemma 10 guarantees that ${\mathbf{D}}$ and ${\mathbf{S}}$ satisfy the rank assumption of Lemma 9. To conclude the proof, we bound the second term of Lemma 9, using the bounds of Lemma 10 and $\rho_{\mathbf{E}}\leq\min\{\omega,n-k\}\leq\omega$:

(of Theorem 7) The proof is similar to the proof of Theorem 6, using Lemma 11 instead of Lemma 10. Let $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleFrobeniusSampling\left({\mathbf{U}}_{{\mathbf{A}}},{\mathbf{E}},r\right)$ We bound the second term of Lemma 9, using the bounds of Lemma 11:

2 𝐁𝐁{\mathbf{B}}-Agnostic Coreset Construction

All the coreset construction algorithms that we presented so far carefully construct the coreset using knowledge of the response vector. If the algorithm does not need knowledge of ${\mathbf{B}}$ to construct the coreset, and yet can provide an approximation guarantee for every ${\mathbf{B}}$, then the algorithm is ${\mathbf{B}}$-agnostic. A ${\mathbf{B}}$-agnostic coreset construction algorithm is appealing because the coreset, as specified by the sampling and rescaling matrices ${\mathbf{S}}$ and ${\mathbf{D}}$, can be computed off-line and applied to any ${\mathbf{B}}$. We briefly digress to show how our methods can be extended to develop ${\mathbf{B}}$-agnostic coreset constructions.

The running time of the proposed algorithm is $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)+O\left(rnk^{2}\right)$, where $T\left({\mathbf{U}}_{{\mathbf{A}}}\right)$ is the time needed to compute the left singular vectors of ${\mathbf{A}}$.

The proof is similar to the proof of Theorem 6, except we now construct the sampling and rescaling matrices as $\left[{\mathbf{S}},{\mathbf{D}}\right]=MultipleSpectralSampling\left({\mathbf{U}}_{{\mathbf{A}}},{\mathbf{I}}_{n},r\right)$. To bound the second term in Lemma 9, we use

The above bound decreases with $r$ and holds for any ${\mathbf{B}}$, guaranteeing a constant-factor approximation with a constant fraction of the data. The approximation ratio is $O(n/r)$, which seems quite weak. In the next section, we show that this result is indeed tight.

Lower Bounds on Coreset Size

We have just seen a ${\mathbf{B}}$-agnostic coreset construction algorithm with a rather weak worst case guarantee of $O(n/r)$ approximation error. We will now show that no deterministic ${\mathbf{B}}$-agnostic coreset construction algorithm can guarantee a better error (Theorem 13) by providing lower bounds on coreset size as a function of approximation error. These results are also summarized in Table 2.

provides another ${\mathbf{B}}$-agnostic coreset construction algorithm with $r=O(k\log k/\epsilon^{2})$. For a fixed ${\mathbf{B}}$, the method in delivers a probabilistic bound on the approximation error. However, there are target matrices ${\mathbf{B}}$ for which the bound fails by an arbitrarily large amount. The probabilistic algorithms get away with this by brushing all these (possibly large) errors into a low probability event, with respect to random choices made in the algorithm. So, in some sense, these algorithms are not ${\mathbf{B}}$-agnostic, in that they do not construct a coreset which works well for all ${\mathbf{B}}$ with some (say) constant probability. Nevertheless, the fact that they give a constant probability of success for a fixed but unknown ${\mathbf{B}}$ makes these algorithms interesting and useful. We will give a lower bound on the approximation ratio of such algorithms as well, for a given probability of success (Theorem 14). Finally, we will give lower bounds on the size of the coreset for the general (non-agnostic) multiple regression setting (Theorems 15 and 16).

We first present the lower bound for simple regression. Recall that a coreset construction algorithm is ${\mathbf{b}}$-agnostic if it constructs a coreset without knowledge of ${\mathbf{b}}$, and then provides an approximation guarantee for every ${\mathbf{b}}$. We show that no coreset can work for every ${\mathbf{b}}$; therefore a ${\mathbf{b}}$-agnostic coreset will be bad for some vector ${\mathbf{b}}$. In fact, there exists a matrix ${\mathbf{A}}$ such that every coreset has an associated “bad” ${\mathbf{b}}$.

There exists a matrix ${\mathbf{A}}\in\R^{n\times d}$ such that for every coreset ${\mathbf{C}}\in\R^{r\times d}$ of size $r\leq n$, there exists ${\mathbf{b}}\in\R^{n}$ (depending on ${\mathbf{C}}$) for which

Let ${\mathbf{P}}_{\mathbf{A}}$ project onto the columns of ${\mathbf{A}}$ and ${\mathbf{P}}_{{\mathbf{A}}^{(1)}}$ project onto the first column of ${\mathbf{A}}$. The following sequence establishes the result:

We now consider randomized algorithms that construct a coreset without looking at ${\mathbf{b}}$ (e.g. ). These algorithms work for any fixed (but unknown) ${\mathbf{b}}$, and deliver a probabilistic approximation guarantee for any single fixed ${\mathbf{b}}$; in some sense they are ${\mathbf{b}}$-agnostic. By the previous discussion, the returned coreset must fail for some ${\mathbf{b}}$, i.e., the probabilistic guarantee does not hold for all ${\mathbf{b}}$, and, when it fails, it could do so with very bad error. We will now present a lower bound on the approximation accuracy of such existing randomized algorithms for coreset construction, even for a single ${\mathbf{b}}$.

First, we define randomized coreset construction algorithms. Let ${\mathbf{C}}_{1},{\mathbf{C}}_{2},\ldots,{\mathbf{C}}_{\left({{n}\atop{r}}\right)}$ be the $\left({{n}\atop{r}}\right)$ different coresets of size $r$. A randomized algorithm assigns probabilities $p_{1},p_{2},\ldots,p_{\left({{n}\atop{r}}\right)}$ to each coreset, and selects one according to these probabilities. The probabilities $p_{i}$ may depend on ${\mathbf{A}}$. The algorithm is ${\mathbf{b}}$-agnostic if the probabilities $p_{i}$ do not depend on ${\mathbf{b}}$. As usual, let $r$ be the size of the coreset.

2 Lower Bounds for Non-Agnostic Multiple Regression

For both the spectral and the Frobenius norm, we now consider non-agnostic unconstrained multiple regression, and give lower bounds for coresets of size $r>d=\hbox{\rm rank}({\mathbf{A}})$ (for simplicity, we set $\hbox{\rm rank}({\mathbf{A}})=d$). The results are presented in Theorems 15 and 16.

First, we need some results from . Consider the matrix

where ${\mathbf{e}}_{i}\in\R^{\omega}$ are the standard basis vectors. Then, let ${\mathbf{B}}={\mathbf{H}}^{\textsc{T}}\in\R^{(\omega-1)\times\omega}$. Theorem 34 in (with $\alpha=1$) argues the following: given ${\mathbf{B}}$ and any sampling matrix ${\mathbf{S}}\in\R^{r\times(\omega-1)}$ and diagonal rescaling matrix ${\mathbf{D}}\in\R^{r\times r}$, with $\hat{\mathbf{C}}={\mathbf{D}}{\mathbf{S}}{\mathbf{B}}$ (rescaled sampled coreset of ${\mathbf{B}}$), and any $k$ with $1\leq k\leq\omega-1$,

In the above, $\Pi_{\hat{\mathbf{C}},k}({\mathbf{B}})\in\R^{(\omega-1)\times\omega}$ of rank $k$ is the best rank-$k$ approximation to ${\mathbf{B}}$ (in the spectral norm) whose rows lie in the span of all the rows in $\hat{\mathbf{C}}$ (the row-space of $\hat{\mathbf{C}}$); and, ${\mathbf{B}}_{k}\in\R^{(\omega-1)\times\omega}$ of rank $k$ is the best rank-$k$ approximation to ${\mathbf{B}}$ (which could be computed via the truncated SVD of ${\mathbf{B}}$).Actually, ${\mathbf{D}}$ is irrelevant here because the row-space of ${\mathbf{S}}{\mathbf{B}}$ is the same as the row space of ${\mathbf{D}}{\mathbf{S}}{\mathbf{B}}$.

Since $\Pi_{\hat{\mathbf{C}},k}({\mathbf{B}})$ is the best rank-$k$ approximation to ${\mathbf{B}}$ in the row-space of $\hat{\mathbf{C}}$, it follows that

for any $\mathbf{X}$$\in\R^{(\omega-1)\times r}$ with rank at most $k$ (because ${\mathbf{X}}\hat{\mathbf{C}}$ will have rank at most $k$ and is in the row space of $\hat{\mathbf{C}}$). Set ${\mathbf{X}}={\mathbf{U}}_{{\mathbf{B}},k}({\mathbf{D}}{\mathbf{S}}{\mathbf{U}}_{{\mathbf{B}},k})^{\dagger}$, where ${\mathbf{U}}_{{\mathbf{B}},k}\in\R^{(\omega-1)\times k}$ has $k$ columns which are the top-$k$ left singular vectors of ${\mathbf{B}}$. It is easy to verify that ${\mathbf{X}}$ has the correct dimensions and rank at most $k$. Since $\hat{\mathbf{C}}={\mathbf{D}}{\mathbf{S}}{\mathbf{B}}$, we have that

To conclude the proof, observe that ${\mathbf{B}}_{d}={\mathbf{U}}_{{\mathbf{B}},d}{\mathbf{U}}_{{\mathbf{B}},d}^{\textsc{T}}{\mathbf{B}}={\mathbf{A}}{\mathbf{A}}^{\dagger}{\mathbf{B}}={\mathbf{A}}{\mathbf{X}}_{opt}$.

As $\alpha\rightarrow 0$ and $n\rightarrow\infty$ the lower bound is $1+d/r$.

First, we need some results from . For any integer $\gamma>1$ and any integer $k\geq 1$, Theorem 36 in exhibits a matrix ${\mathbf{B}}\in\R^{\gamma k\times(\gamma+1)k}$ such that for any sampling matrix ${\mathbf{S}}\in\R^{r\times\gamma k}$ and diagonal rescaling matrix ${\mathbf{D}}\in\R^{r\times r}$, with $\hat{\mathbf{C}}={\mathbf{D}}{\mathbf{S}}{\mathbf{B}}$ (rescaled sampled coreset of ${\mathbf{B}}$), any $\alpha>0$, and any $r\geq 1$,

The matrix ${\mathbf{B}}$ is constructed as follows. Recall that $\gamma$ is any positive integer with $\gamma>1$. Let ${\mathbf{A}}$ have dimensions $(\gamma+1)\times\gamma$ and be constructed as follows.

where ${\mathbf{e}}_{i}\in\R^{\gamma+1}$ are the standard basis vectors. Now construct ${\mathbf{H}}$ to be block diagonal, with $k$ copies of ${\mathbf{A}}$ along its diagonal; so, the dimensions of ${\mathbf{H}}$ are $(\gamma+1)k\times\gamma k$. Then, ${\mathbf{B}}={\mathbf{H}}^{\textsc{T}}$.

In the above, $\Pi_{\hat{\mathbf{C}},k}({\mathbf{B}})\in\R^{\gamma k\times(\gamma+1)k}$ of rank $k$ is the best rank-$k$ approximation to ${\mathbf{B}}$ (in the Frobenius norm) whose rows lie in the span of all the rows in $\hat{\mathbf{C}}$ (the row-space of $\hat{\mathbf{C}}$); and, ${\mathbf{B}}_{k}\in\R^{\gamma k\times(\gamma+1)k}$ of rank $k$ is the best rank-$k$ approximation to ${\mathbf{B}}$ (which could be computed via the truncated SVD of ${\mathbf{B}}$). Since $\Pi_{\hat{\mathbf{C}},k}({\mathbf{B}})$ is the best rank-$k$ approximation to ${\mathbf{B}}$ in the row-space of $\hat{\mathbf{C}}$, it follows that

for any $\mathbf{X}$$\in\R^{\gamma k\times r}$ with rank at most $k$ (because ${\mathbf{X}}\hat{\mathbf{C}}$ will have rank at most $k$ and is in the row space of $\hat{\mathbf{C}}$). Set ${\mathbf{X}}={\mathbf{U}}_{{\mathbf{B}},k}({\mathbf{D}}{\mathbf{S}}{\mathbf{U}}_{{\mathbf{B}},k})^{\dagger}$, where ${\mathbf{U}}_{{\mathbf{B}},k}\in\R^{\gamma k\times k}$ has $k$ columns which are the top-$k$ left singular vectors of ${\mathbf{B}}$. It is easy to verify that ${\mathbf{X}}$ has the correct dimensions and rank at most $k$. Since $\hat{\mathbf{C}}={\mathbf{D}}{\mathbf{S}}{\mathbf{B}}$, we have that

We now construct the regression problem which proves the lower bound in the theorem. Let ${\mathbf{A}}={\mathbf{U}}_{{\mathbf{B}},d}\in\R^{\gamma d\times d}$ (i.e., we choose $k=d$ in the above discussion), $n=\gamma d$ (i.e. $n$ is a multiple of $d$ in the regression problem), and $\omega=(\gamma+1)d$. ${\mathbf{B}}$ is as we described above. Suppose a coreset construction algorithm gives sampling and rescaling matrices ${\mathbf{S}}$ and ${\mathbf{D}}$, for a coreset of size $r>d$. So, the coreset regression is with ${\mathbf{C}}={\mathbf{D}}{\mathbf{S}}{\mathbf{A}}\in\R^{r\times d}$ and ${\mathbf{D}}{\mathbf{S}}{\mathbf{B}}\in\R^{r\times\omega}$. The solution to the coreset regression is

To conclude the proof, observe that ${\mathbf{B}}_{d}={\mathbf{U}}_{{\mathbf{B}},d}{\mathbf{U}}_{{\mathbf{B}},d}^{\textsc{T}}{\mathbf{B}}={\mathbf{A}}{\mathbf{A}}^{\dagger}{\mathbf{B}}={\mathbf{A}}{\mathbf{X}}_{opt}$ and $\omega=n+d$.

Open problems

An important open problem arises in our work: can we determine the minimum size of a coreset that provides a $(1+\epsilon)$ relative-error guarantee for simple linear regression? We conjecture that $\Omega\left(k/\epsilon\right)$ is a lower bound, which will make our results almost tight. Certainly, coresets of size exactly $k$ cannot be guaranteed: consider two data points $(1,1),(-1,1)$. The optimal regression is zero; however any coreset of size one will give non-zero regression.

Christos Boutsidis acknowledges the support from XDATA program of the Defense Advanced Research Projects Agency (DARPA), administered through Air Force Research Laboratory contract FA8750-12-C-0323. Petros Drineas and Malik Magdon-Ismail have been supported by NSF CCF 1016501, NSF DMS 1008983, and NSF CCF CAREER 824684.

References

Appendix A Linear Algebra Background

The Singular Value Decomposition (SVD) of a matrix ${\mathbf{A}}\in\R^{n\times d}$ of rank $k$ is a decomposition

The singular values $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{k}>0$ are contained in the diagonal matrix ${\mathbf{\Sigma}}_{\mathbf{A}}\in\R^{k\times k}$; ${\mathbf{U}}_{\mathbf{A}}\in\R^{n\times k}$ contains the left singular vectors of ${\mathbf{A}}$; and ${\mathbf{V}}_{{\mathbf{A}}}\in\R^{d\times k}$ contains the right singular vectors. The Moore-Penrose pseudo-inverse of ${\mathbf{A}}$ is ${\mathbf{A}}^{\dagger}={\mathbf{V}}_{{\mathbf{A}}}{\mathbf{\Sigma}}_{\mathbf{A}}^{-1}{\mathbf{U}}_{{\mathbf{A}}}^{\textsc{T}}.$ Given an orthonormal matrix ${\mathbf{U}}_{\mathbf{A}}\in\R^{n\times k}$, the perpendicular matrix ${\mathbf{U}}_{\mathbf{A}}^{\perp}\in\R^{n\times(n-k)}$ to ${\mathbf{U}}_{{\mathbf{A}}}$ satisfies: $({\mathbf{U}}_{\mathbf{A}}^{\perp})^{\textsc{T}}{\mathbf{U}}_{\mathbf{A}}^{\perp}={\mathbf{I}}_{n-k}$, ${\mathbf{U}}_{\mathbf{A}}^{\textsc{T}}{\mathbf{U}}_{\mathbf{A}}^{\perp}=\bm{0}_{k\times(n-k)}$, and ${\mathbf{U}}_{\mathbf{A}}{\mathbf{U}}_{{\mathbf{A}}}^{\textsc{T}}+{\mathbf{U}}_{\mathbf{A}}^{\perp}({\mathbf{U}}_{\mathbf{A}}^{\perp})^{\textsc{T}}={\mathbf{I}}_{n}$. All the singular values of both ${\mathbf{U}}_{\mathbf{A}}$ and ${\mathbf{U}}_{\mathbf{A}}^{\perp}$ are equal to one. Given ${\mathbf{U}}_{{\mathbf{A}}}$, ${\mathbf{U}}_{\mathbf{A}}^{\perp}$ can be computed in deterministic $O\left(n\left(n-k\right)^{2}\right)$ time via the QR factorization.

Let ${\mathbf{X}}$ and ${\mathbf{Y}}$ be two $n\times d$ matrices. If ${\mathbf{X}}{\mathbf{Y}}^{\textsc{T}}=\bm{0}_{n\times n}$ or ${\mathbf{X}}^{\textsc{T}}{\mathbf{Y}}=\bm{0}_{d\times d}$, then

Appendix B Algorithms and Proofs of Lemmas 10 and 11

We now provide all the details of the proofs and the corresponding algorithms of Lemmas 10 and 11. Those results, which have been described in detail in , are slight extensions of two algorithms presented in , which themselves extend the original spectral sparsification result of Batson, Spielman, and Srivastava . More specifically, Lemma 20 below - in some sense - generalizes Lemma 2; indeed, setting ${\mathbf{V}}={\mathbf{Q}}^{\textsc{T}}:={\mathbf{U}}$ in Lemma 20 gives Lemma 2. Lemma 19 below also describes a deterministic algorithm for sampling columns from two matrices but the goal here is to optimize different spectral properties in the sampled matrices.

In this section of the Appendix, we will slightly abuse notation by denoting with $\hat{\mathbf{S}}\in\R^{n\times r}$ a sampling matrix which samples columns - not rows - from matrices. We will later use ${\mathbf{S}}=\hat{\mathbf{S}}^{\textsc{T}}$ to be consistent with the notation used throughout the paper.

We write $[{\mathbf{D}},\hat{\mathbf{S}}]=DeterministicSamplingI({\mathbf{V}}^{\textsc{T}},{\mathbf{B}},r)$ to denote this procedure.

Algorithm 5 is a greedy technique that selects columns one at a time. To describe the algorithm in more detail, it is convenient to view the input matrices as two sets of $n$ vectors,

Given $k$ and $r>k$, introduce the iterator $\tau=0,1,2,...,r-1,$ and define the parameter

For a square symmetric matrix ${\mathbf{A}}\in\R^{k\times k}$ with eigenvalues $\lambda_{1},\ldots,\lambda_{k}$, ${\mathbf{v}}\in\R^{k}$ and ${\textsc{l}}\in\R$, define

and let $L({\mathbf{v}},\delta_{L},{\mathbf{A}},{\textsc{l}})$ be defined as

where ${\textsc{l}}^{\prime}={\textsc{l}}+\delta_{L}={\textsc{l}}+1.$ For a vector ${\mathbf{z}}$ and scalar $\delta>0$, define the function

At each iteration $\tau$, the algorithm selects $i_{\tau}$, $t_{\tau}>0$ for which

The running time of the algorithm is dominated by the search for an index $i_{\tau}$ satisfying

If ${\mathbf{Q}}={\mathbf{I}}_{n}$, it runs in $O(rk^{2}n)$; we write $[{\mathbf{D}},\hat{\mathbf{S}}]=DeterministicSamplingII({\mathbf{V}}^{\textsc{T}},{\mathbf{Q}},r)$ for this procedure.

If ${\mathbf{\Psi}}={\mathbf{I}}_{n}$, the running time of the algorithm reduces to $T_{SVD}\left({\mathbf{Y}}\right)+O\left(rn\rho_{{\mathbf{Y}}}^{2}\right)$. We write $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleSpectralSampling\left({\mathbf{Y}},{\mathbf{\Psi}},r\right)$ to denote such a deterministic procedure.

(a) uses Lemma 18. To obtain the first inequality in the lemma we need to take ${\mathbf{S}}=\hat{\mathbf{S}}^{\textsc{T}}$ and observe that $\mbox{}\|({\mathbf{Y}}^{\textsc{T}}\hat{\mathbf{S}}{\mathbf{D}})^{\dagger}\|_{2}=\mbox{}\|({\mathbf{D}}{\mathbf{S}}{\mathbf{Y}})^{\dagger}\|_{2}$, and $\mbox{}\|({\mathbf{Y}}^{\textsc{T}})^{\dagger}\|_{2}=\mbox{}\|{\mathbf{Y}}^{\dagger}\|_{2}$. We now prove the second inequality in the lemma,

To obtain the second inequality in the lemma we need to take ${\mathbf{S}}=\hat{\mathbf{S}}^{\textsc{T}}$ and use $\mbox{}\|{\mathbf{\Psi}}^{\textsc{T}}\hat{\mathbf{S}}{\mathbf{D}}\|_{2}=\mbox{}\|{\mathbf{D}}{\mathbf{S}}{\mathbf{\Psi}}\|_{2},$ and $\mbox{}\|{\mathbf{\Psi}}^{\textsc{T}}\|_{2}=\mbox{}\|{\mathbf{\Psi}}\|_{2}$.

B.2 Proof of Lemma 11

If ${\mathbf{\Psi}}={\mathbf{I}}_{n}$, the running time of the algorithm reduces to $T_{SVD}\left({\mathbf{Y}}\right)+O\left(rn\rho_{{\mathbf{Y}}}^{2}\right)$. We write $\left[{\mathbf{D}},{\mathbf{S}}\right]=MultipleFrobeniusSampling\left({\mathbf{Y}},{\mathbf{\Psi}},r\right)$ to denote such a deterministic procedure.

which by taking ${\mathbf{S}}=\hat{\mathbf{S}}^{\textsc{T}}$ gives the second inequality in the lemma,

Now we prove the first inequality in the lemma,

(a) uses Lemma 18. To conclude, use $\mbox{}\|({\mathbf{Y}}^{\textsc{T}}{\mathbf{S}}{\mathbf{D}})^{\dagger}\|_{2}=\mbox{}\|({\mathbf{D}}{\mathbf{S}}{\mathbf{Y}})^{\dagger}\|_{2};$ $\mbox{}\|({\mathbf{V}}_{\mathbf{Y}}{\mathbf{\Sigma}}_{\mathbf{Y}})^{\dagger}\|_{2}=\mbox{}\|({\mathbf{Y}}^{\textsc{T}})^{\dagger}\|_{2}=\mbox{}\|({\mathbf{Y}})^{\dagger}\|_{2}$.