An Explicit Sampling Dependent Spectral Error Bound for Column Subset Selection

Introduction

where $P_{C}=CC^{\dagger}$ denotes the projection onto the column space of $C$ with $C^{\dagger}$ being the pseudo inverse of $C$ and $\xi=2$ or $F$ denotes the spectral norm or the Frobenius norm. In this paper, we are particularly interested in the spectral norm reconstruction with respect to a target rank $k$.

To the best of our knowledge, this is the first such kind of error bound for CSS. Compared with existing error bounds, the sampling dependent error bound brings us several benefits: (i) it allows us to better understand the tradeoff in the spectral error of reconstruction due to sampling probabilities, complementary to a recent result on the tradeoff from a statistical perspective (Ma et al., 2014) for least square regression; (ii) it implies that a distribution with sampling probabilities proportional to the square root of the SLS is always better than the uniform sampling, and is potentially better than that proportional to the SLS when they are skewed; (iii) it motivates an optimization approach by solving a constrained optimization problem related to the error bound to attain better performance. In addition to the theoretical analysis, we also develop an efficient bisection search algorithm to solve the constrained optimization problem for finding better sampling probabilities.

By combining our analysis with recent developments for spectral norm reconstruction of CSS (Boutsidis et al., 2011), we also establish the same error bound for an exact rank-$k$ approximation, i.e.,

where $\Pi^{2}_{C,k}(A)$ is the best approximation to $A$ within the column space of $C$ that has rank at most $k$.

The remainder of the paper is organized as follows. We review some closely related work in Section 2, and present the main result in Section 4 with some preliminaries in Section 3. We conduct some empirical studies in Section 5 and present the detailed analysis in Section 6. Finally, conclusion is made.

Related Work

In this section, we review some previous work on CSS, low-rank matrix approximation, and other closely related work on randomized algorithms for matrices. We focus our discussion on the spectral norm reconstruction.

with a running time $O(mnk\log n)$. The same error bound was also achieved by using volume sampling (Deshpande & Rademacher, 2010). The running time of volume sampling based algorithms can be made close to linear to the size of the target matrix. Boutsidis et al. (2009) proposed a two-stage algorithm for selecting exactly $k$ columns and provided error bounds for both the spectral norm and the Frobenius norm, where in the firs stage $\Theta(k\log k)$ columns are sampled based on a distribution related to the SLS and more if for the spectral norm reconstruction and in the second stage $k$ columns are selected based on the rank revealing QR factorization. The spectral error bound in this work that holds with a constant probability 0.8 is following:

The time complexity of their algorithm (for the spectral norm reconstruction) is given by $O(\min(mn^{2},m^{2}n))$ since it requires SVD of the target matrix for computing the sampling probabilities.

Although our sampling dependent error bound is not directly comparable to these results, our analysis exhibits that the derived error bound could be much better than that in (2). When the SLS are nonuniform, our new sampling distributions could lead to a better result than (3). Most importantly, the sampling probabilities in our algorithm are only related to the SLS and that can be computed more efficiently (e.g., exactly in $O(T_{V_{k}})$ or approximately in $O(mn\log n)$ (Drineas et al., 2012)). In simulations, we observe that the new sampling distributions could yield even better spectral norm reconstruction than the deterministic selection criterion in (Boutsidis et al., 2011), especially when the SLS are nonuniform.

There are much more work on studying the Frobenius norm reconstruction of CSS (Drineas et al., 2006a; Guruswami & Sinop, 2012; Boutsidis et al., 2011; Drineas et al., 2008; Boutsidis et al., 2009). For more references, we refer the reader to the survey (Mahoney, 2011). It remains an interesting question to establish sampling dependent error bounds for other randomized matrix algorithms.

Preliminaries

Finally, we let $\mathbf{s}^{*}=(s^{*}_{1},\ldots,s^{*}_{n})$ denote the SLS of $A$ relative to the best rank-$k$ approximation to $A$ (Mahoney, 2011), i.e., $s_{i}^{*}=\|\mathbf{v}_{i}\|_{2}^{2}$. It is not difficult to show that $\sum_{i=1}^{n}s_{i}^{*}=k$.

Main Result

Before presenting our main result, we first characterize scores in $\mathbf{s}$ by two quantities as follows:

Both quantities compare $\mathbf{s}$ to the SLS $\mathbf{s}^{*}$. With $c(\mathbf{s})$ and $q(\mathbf{s})$, we are ready to present our main theorem regarding the spectral error bound.

where $\sigma_{k+1}=\|A-A_{k}\|_{2}$ is the $(k+1)$th singular value of $A$.

Remark: Clearly, the spectral error bound and the successful probability in Theorem 1 depend on the quantities $c(\mathbf{s})$ and $q(\mathbf{s})$. In the subsection below, we study the two quantities to facilitate the understanding of the result in Theorem 1.

The result in Theorem 1 implies that the smaller the quantities $c(\mathbf{s})$ and $q(\mathbf{s})$, the better the error bound. Therefore, we first study when $c(\mathbf{s})$ and $q(\mathbf{s})$ achieve their minimum values. The key results are presented in the following two lemmas with their proofs deferred to the supplement.

The set of scores in $\mathbf{s}$ that minimize $q(\mathbf{s})$ is given by $s_{i}\propto\sqrt{s^{*}_{i}}$, i.e., $s_{i}=\frac{k\sqrt{s^{*}_{i}}}{\sum_{i=1}^{n}\sqrt{s_{i}^{*}}}$.

Remark: The sampling distribution with probabilities that are proportional to the square root of $s_{i}^{*},i\in[n]$ falls in between the uniform sampling and the leverage-based sampling.

$c(\mathbf{s})\geq 1,\forall\mathbf{s}$ such that $\sum_{i=1}^{m}s_{i}=k$. The set of scores in $\mathbf{s}$ that minimize $c(\mathbf{s})$ is given by $s_{i}=s_{i}^{*}$, and the minimum value of $c(\mathbf{s})$ is $1$.

Next, we discuss three special samplings with $\mathbf{s}$ (i) proportional to the square root of the SLS, i.e., $s_{i}\propto\sqrt{s_{i}^{*}}$ (referred to as square-root leverage-based sampling or sqL-sampling for short), (ii) equal to the SLS, i.e., $s_{i}=s_{i}^{*}$ (referred to as leverage-based sampling or L-sampling for short), and (iii) equal to uniform scalars $s_{i}=k/n$ (referred to as uniform sampling or U-sampling for short). Firstly, if $s_{i}\propto\sqrt{s_{i}^{*}}$ , $q(\mathbf{s})$ achieves its minimum value and we have the two quantities written as

Secondly, if $s_{i}\propto s_{i}^{*}$, then $c(\mathbf{s})$ achieves its minimum value and we have the two quantities written as

Lastly, we consider the uniform sampling $s_{i}=\frac{k}{n}$ . Then the two quantities become

Similarly, if $\mathbf{s}_{*}$ is flat, $q_{U}=\sqrt{\frac{n}{k}}$ and $c_{U}=1$. Moreover, it is interesting to compare the two quantities for the sqL-sampling in (7) and for the uniform sampling in (9).

From the above discussions, we can see that when $\mathbf{s}_{*}$ is a flat vector, there is no difference between the three sampling scores for $\mathbf{s}$. The difference comes from when $\mathbf{s}_{*}$ tends to be skewed. In this case, $s_{i}\propto\sqrt{s_{i}^{*}}$ works almost for sure better than uniform distribution and could also be potentially better than $s_{i}\propto s_{i}^{*}$ according to the sampling dependent error bound in Theorem 1. A similar tradeoff between the L-sampling and U-sampling but with a different taste was observed in (Ma et al., 2014), where they showed that for least square approximation by CSS leveraging-based least square estimator could have a large variance when there exist very small SLS. Nonetheless, our bound here exhibits more insights, especially on the sqL-sampling. More importantly, the sampling dependent bound renders the flexibility in choosing the sampling scores by adjusting them according to the distribution of the SLS. In next subsection, we present an optimization approach to find better sampling scores. In Figure 1, we give a quick view of different sampling strategies.

2 Optimizing the error bound

As indicated by the result in Theorem 1, in order to achieve a good performance, we need to make a balance between $c(\mathbf{s})$ an $q(\mathbf{s})$, where $c(\mathbf{s})$ affects not only the error bound but also the successful probability. To address this issue, we propose a constrained optimization approach. More specifically, to ensure that the failure probability is no more than $3\delta$, we impose the following constraint on $c(\mathbf{s})$

Then we cast the problem into minimizing $q(\mathbf{s})$ under the constraint in (10), i.e.,

It is easy to verify that the optimization problem in (11) is convex. Next, we develop an efficient bisection search algorithm to solve the above problem with a linear convergence rate. To this end, we introduce a slack variable $t$ and rewrite the optimization problem in (11) as

We now find the optimal solution by performing bisection search on $t$. Let $t_{\max}$ and $t_{\min}$ be the upper and lower bounds for $t$. We set $t=(t_{\min}+t_{\max})/2$ and decide the feasibility of $t$ by simply computing the quantity

Evidently, $t$ is a feasible solution if $f(t)\leq k$ and is not if $f(t)>k$. Hence, we will update $t_{\max}=t$ if $f(t)\leq k$ and $t_{\min}=t$ if $f(t)>k$. To run the bisection algorithm, we need to decide initial $t_{\min}$ and $t_{\max}$. We can set $t_{\min}=0$. To compute $t_{\max}$, we make an explicit construction of $\mathbf{s}$ by distributing the $(1-\gamma^{-1})$ share of the largest element of $\mathbf{s}_{*}$ to the rest of the list. More specifically, let $j$ be the index for the largest entry in $\mathbf{s}^{*}$. We set $s_{j}=\|\mathbf{s}^{*}\|_{\infty}\gamma^{-1}$ and $s_{i}=s_{i}^{*}+(1-\gamma^{-1})\|\mathbf{s}^{*}\|_{\infty}/(n-1)$ for $i\neq j$. Evidently, this solution satisfies the constraints $s_{i}^{*}\leq\gamma s_{i},i\in[n]$ for $\gamma\geq 1$. With this construction, we can show that

Therefore, we set initial $t_{\max}$ to the value in R.H.S of the above inequality. Given the optimal value of $t=t_{*}$ we compute the optimal value of $s_{i}$ by $s_{i}=\frac{s_{i}^{*}}{\min(\gamma,t_{*}\sqrt{s_{i}^{*}})}.$ The corresponding sampling distribution clearly lies between L-sampling and sqL-sampling. In particular, when $\gamma=1$ the resulting sampling distribution is L-sampling due to Lemma 2 and when $\gamma\rightarrow\infty$ the resulting sampling distribution approaches sqL-sampling.

3 Subsequent Applications

Next, we discuss two subsequent applications of CSS, one for low rank approximation and one for least square approximation.

where $\epsilon(\mathbf{s})$ is given in Theorem 1.

It is worth pointing out that in this case the SLS $\mathbf{s}^{*}=(s^{*}_{1},\ldots,s_{m}^{*})$ are computed based on the the left singular vectors $U$ of $A$ by $s_{i}^{*}=\|U_{i*}\|_{2}^{2}$, where $U_{i*}$ is the $i$-th row of $U$. One might be interested to see whether we can apply our analysis to derive a sampling dependent error bound for the approximation error $\|\mathbf{x}_{opt}-\widehat{\mathbf{x}}_{opt}\|_{2}$ similar to previous bounds of the form $\|\mathbf{x}_{opt}-\widehat{\mathbf{x}}_{opt}\|_{2}\leq\frac{\epsilon}{\sigma_{min}(A)}\|A\mathbf{x}_{top}-\mathbf{b}\|_{2}$. Unfortunately, naively combining our analysis with previous analysis is a worse case analysis, and consequentially yields a worse bound. The reason will become clear in our later discussions. However, the statistical analysis in (Ma et al., 2014) does indicate that $\widehat{\mathbf{x}}_{opt}$ by using sqL-sampling could have smaller variance than that using L-sampling.

Numerical Experiments

Before delving into the detailed analysis, we present some experimental results. We consider synthetic data with the data matrix $A$ generated from one of the three different classes of distributions introduced below, allowing the SLS vary from nearly uniform to very nonuniform.

Nearly uniform SLS (GA). Columns of $A$ are generated from a multivariate normal distribution $\mathcal{N}(\mathbf{1}_{m},\Sigma)$, where $\Sigma_{ij}=2*0.5^{|i-j|}$. This data is referred to as GA data.

Moderately nonuniform SLS ($T_{3}$). Columns of $A$ are generated from a multivariate $t$-distribution with $3$ degree of freedom and covariance matrix $\Sigma$ as before. This data is referred to as $T_{3}$ data.

Very nonuniform SLS ($T_{1}$). Columns of $A$ are generated from a multivariate $t$-distribution with $1$ degree of freedom and covariance matrix $\Sigma$ as before. This data is referred to as $T_{1}$ data.

These distributions have been used in (Ma et al., 2014) to generate synthetic data for empirical evaluations.

More results including relative error versus varying size $n$ of the target matrix, performance on a real data set and the Frobenius norm reconstruction error can be found in supplement.

Analysis

In this section, we present major analysis of Theorem 1 and Theorem 2 with detailed proofs included in supplement. The key to our analysis is the following Theorem.

Let $Y,\Omega_{1},\Omega_{2}$ be defined in (5). Assume that $\Omega_{1}$ has full row rank. We have

The first inequality was proved in (Halko et al., 2011) (Theorem 9.1) and the second inequality is credited to (Boutsidis et al., 2011) (Lemma 3.2) In fact, the first inequality is implied by the second inequality. . Previous work on the spectral norm analysis also start from a similar inequality as above. They bound the second term by using $\|\Sigma_{2}\Omega_{2}\Omega_{1}^{\dagger}\|_{2}\leq\|\Sigma_{2}\Omega_{2}\|_{2}\|\Omega_{1}^{\dagger}\|_{2}$ and then bound the two terms separately. However, we will first write $\left\|\Sigma_{2}\Omega_{2}\Omega_{1}^{\dagger}\right\|_{2}=\|\Sigma_{2}\Omega_{2}\Omega_{1}^{\top}(\Omega_{1}\Omega_{1}^{\top})^{-1}\|_{2}$ using the fact $\Omega_{1}$ has full row rank, and then bound $\|(\Omega_{1}\Omega_{1}^{\top})^{-1}\|_{2}$ and $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$ separately. To this end, we will apply the Matrix Chernoff bound as stated in Theorem 4 to bound $\|(\Omega_{1}\Omega_{1}^{\top})^{-1}\|_{2}$ and apply the matrix Bernstein inequality as stated in Theorem 5 to bound $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$.

Following immediately from Theorem 3, we have

where the last inequality uses the fact $\sqrt{a^{2}+b^{2}}\leq a+b$. Below we bound $\lambda_{\min}(\Omega_{1}\Omega_{1}^{\top})$ from below and bound $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$ from above.

We will utilize Theorem 4 to bound $\lambda_{\min}(\Omega_{1}\Omega_{1}^{\top})$. Define $X_{i}=\mathbf{v}_{i}\mathbf{v}_{i}^{\top}/s_{i}$. It is easy to verify that

We will utilize Theorem 5 to bound $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$. Define $Z_{j}=\mathbf{u}_{i_{j}}\mathbf{v}_{i_{j}}^{\top}/s_{i_{j}}$. Then

We can complete the proof of Theorem 1 by combining the bounds for $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$ and $\lambda^{-1}_{\min}(\Omega_{1}\Omega_{1}^{\top})$ and by setting $\delta=1/2$ in Theorem 6 and using union bounds.

Discussions and Open Problems

From the analysis, it is clear that the matrix Bernstein inequality is the key to derive the sampling dependent bound for $\|\Omega_{2}\Omega_{1}^{\top}\|_{2}$. For bounding $\lambda_{\min}(\Omega_{1}\Omega_{1}^{\top})$, similar analysis using matrix Chernoff bound has been exploited before for randomized matrix approximation (Gittens, 2011).

Since Theorem 3 also holds for the Frobenius norm, it might be interested to see whether we can derive a sampling dependent Frobenius norm error bound that depends on $c(\mathbf{s})$ and $q(\mathbf{s})$, which, however, still remains as an open problem for us. Nonetheless, in experiments (included in the supplement) we observe similar phenomena about the performance of L-sampling, U-sampling and sqL-sampling.

Finally, we briefly comment on the analysis for least square approximation using CSS. Previous results (Drineas et al., 2008, 2006b, 2011) were built on the structural conditions that are characterized by two inequalities

The first condition can be guaranteed by Theorem 6 with a high probability. For the second condition, if we adopt a worse case analysis

and bound the first term in R.H.S of the above inequality using Theorem 7, we would end up with a worse bound than existing ones that bound the left term as a whole. Therefore the naive combination can’t yield a good sampling dependent error bound for the approximation error of least square regression.

Conclusions

In this paper, we have presented a sampling dependent spectral error bound for CSS. The error bound brings a new distribution with sampling probabilities proportional to the square root of the statistical leverage scores and exhibits more tradeoffs and insights than existing error bounds for CSS. We also develop a constrained optimization algorithm with an efficient bisection search to find better sampling probabilities for the spectral norm reconstruction. Numerical simulations demonstrate that the new sampling distributions lead to improved performance.

References