Two-way Spectrum Pursuit for $\boldsymbol{CUR}$ Decomposition and Its Application in Joint Column/Row Subset Selection

Matrix factorization provides a concise representation of data. Despite desirable uniqueness conditions and computational simplicity of the well-known singular value decomposition (SVD), it comes with some fundamental shortcomings. The intrinsic structure of data is not inherited to the singular components. Moreover, SVD implies orthogonality on the components which is irrelevant to the underlying structure of the original data. This enforced structure makes the bases, a.k.a. singular vectors, hard to interpret [9]. On the other hand, it is shown that borrowing bases from the actual samples of a dataset provides a robust representation, which can be employed in interesting applications where sampling is their key factor [20]. This problem is studied under the literature of column subset selection problem (CSSP) [5, 2] and CUR decomposition [16, 3]. A general problem for CSSP and CUR can be written in the following form:

where, $\!\small{\boldsymbol{X}}\!\!\in\!\!\mathbb{R}^{N\times M}\!$ is the data matrix containing $M$ data points in an $N$-dimensional space. Here, $\pi_{\mathbb{S}^{c}}^{c}(.)$ and $\pi_{\mathbb{S}^{r}}^{r}(.)$ indicate column space projection and row space projection, respectively. These operators project all columns (rows) to a low-dimensional subspace spanned by selected columns (rows) of matrix $\boldsymbol{X}$ indexed by the set $\mathbb{S}^{c}$ ($\mathbb{S}^{r}$). The chronological order in applying $\pi_{\mathbb{S}^{c}}^{c}(.)$ and $\pi_{\mathbb{S}^{r}}^{r}(.)$ does not affect the problem since these operators are linear. Moreover, substituting $\pi_{\mathbb{S}^{r}}^{r}$ with the identity projection simplifies the problem to CSSP. Fig. 1 illustrates the structure of CUR matrix decomposition as a self-representative approach.

A versatile metric for evaluating the performance of a data subset selection algorithm can be defined by the approximation error resulted from the projection of the entire data to the span of selected rows/columns. How close to the optimal selection an algorithm can reach, is determined by comparing its approximation error to the best low-rank approximation error specified by the spectral decomposition. Recently, we proposed a fast and accurate algorithm for solving CSSP which is called spectrum pursuit (SP) [12].

Inspired by SP, we propose a new algorithm to address the more general case of the CUR matrix decomposition. The main contributions of our paper are summarized as follows:

• •

  A novel algorithm for CUR decomposition, referred to as two-way spectrum pursuit (TWSP), is proposed. TWSP provides an accurate solution for CUR decomposition.

• •

  TWSP enjoys a linear complexity w.r.t. the number of columns and the number of rows of a matrix.

• •

  TWSP is a parameter-free algorithm that only requires the number of desired columns and rows for selection. Thus, TWSP does not require any parameter fine-tuning.

• •

  The TWSP algorithm is put to the test and investigated in a set of synthetic and real experiments.

• •

  The role of the core matrix $\boldsymbol{U}$ in $\boldsymbol{CUR}$ decomposition is illustrated which shows the connection between selected columns and rows. Based on analysis of $\boldsymbol{U}$, an interesting application for joint sensor/channel selection is presented.

Notations: Throughout this paper, vectors, and matrices are denoted by bold lowercase, and bold uppercase letters, respectively. Moreover sets are indicated by blackboard bold characters. The $m^{\text{th}}$ column of matrix $\boldsymbol{X}$ is denoted by $\boldsymbol{X}(:,m)$ and the $n^{\text{th}}$ row of matrix $\boldsymbol{X}$ is denoted by $\boldsymbol{X}(n,:)$.

The rest of paper is organized as follows. Sec. 2 studies the problem of column subset selection as the conventional one-way selection scheme. In Sec. 3, our main work is presented for joint column and row subset selection. Experimental results are exhibited in Sec. 4.

Selecting the most diverse subset of data in an optimal sense is studied vastly [5, 14, 13]. However, these methods do not guarantee that the unselected columns are well represented by the selected ones. Further, outliers are selected with a high probability using such algorithms due to their diversity [12]. A more effective approach is selecting some representatives which are able to approximate the rest of data accurately [7] as defined as a special case of (1). This is an NP-hard problem [18] and there are several efforts for solving this problem [1, 5, 17, 21]. There are computationally expensive approaches based on convex relaxation [7, 6] that are not computationally feasible for large datasets since their complexity is of order $O(M^{3})$, where $M$ is the number of original columns. Recently, we proposed a new algorithm for solving CSSP with a linear complexity which is called Spectrum pursuit (SP) [12]. The SP algorithm finds $K$ columns of $\boldsymbol{X}$ such that their span is close to that of the best rank-$K$ approximation of $\boldsymbol{X}$. SP is an iterative approach where at each iteration one sample selection is optimized such that the ensemble of selected samples describes the whole dataset more accurately. SP finds representatives such that the column space is spanned accurately via consecutive rank-$1$ approximations as theoretically analyzed in [11]. In the present paper, we extend the SP algorithm for selecting columns and rows jointly such that their outer product can represent the whole matrix accurately. A naive approach is applying SP algorithm on the matrix of interest to select a subset of columns and applying SP on its transpose in order to find a subset of rows. However, this approach is not efficient and we will compare it with our proposed approach which is optimized through a joint representation of selected columns and rows.

The introduced joint column/row subset selection in (1) can be written as a CUR decomposition in the following form in which factor matrices must be drawn from actual columns/rows of the original matrix as

In this problem, $\mathbb{X}_{c}$, and $\mathbb{X}_{r}$ indicate the set of normalized columns and rows of matrix $\boldsymbol{X}$, respectively. Here, $\boldsymbol{c}_{k}$ and $\boldsymbol{r}^{T}_{k}$ denote the $k^{th}$ column and the $k^{th}$ row of $\boldsymbol{C}$ and $\boldsymbol{R}$, respectively. In other words, $\mathbb{X}_{c}\!=\!\{\boldsymbol{X}(:,m)/\|\boldsymbol{X}(:,m)\|\}$ for all columns and $\mathbb{X}_{r}\!=\!\{\boldsymbol{X}(n,:)/\|\boldsymbol{X}(n,:)\|\}$ for all rows. Please note that replacing constraints in (2) with orthogonality constraint on $\boldsymbol{c}_{k}$’s and on $\boldsymbol{r}_{k}$’s results in the truncated singular value decomposition (SVD) with $K$ most significant components. In this case, $\boldsymbol{C}$ and $\boldsymbol{R}$ contain the first $K$ left singular vectors and the first $K$ right singular vectors of $\boldsymbol{X}$, respectively. Moreover, the core matrix will be diagonal and the entries will be singular values with diagonal entries as singular values. However, the underlying constraints in (2) turn the problem into a joint subset of row and column selection problem instead of matrix low-rank approximation problem.

To solve this complicated problem, we split it into two consecutive problems for optimization on the $k^{\text{th}}$ selected column/row. Our optimization approach is alternative, i.e., a random subset of columns and rows are picked. Then, one column or row is considered to be replaced with a more efficient one at each iteration. Since scaling a vector does not change its span, without loss of generality we assume that the column or the row subject of the optimization lie on the unit sphere. At each iteration, a rank-$1$ component is optimized characterized by $\boldsymbol{c}\boldsymbol{g}^{T}$ or $\boldsymbol{h}\boldsymbol{r}^{T}$ given by

Matrix $\boldsymbol{C}_{k}$ is the set of selected columns except the $k^{\text{th}}$ one and $\boldsymbol{R}_{k}$ is the set of selected rows except the $k^{\text{th}}$ row. The first subproblem can be solved easily w.r.t. $\boldsymbol{c}$ using singular value decomposition. In other words, $\boldsymbol{cg}^{T}$ and $\boldsymbol{hr}^{T}$ are the best rank-$1$ approximations of the residual $\boldsymbol{E}^{c}$ and $\boldsymbol{E}^{r}$, respectively. The obtained $\boldsymbol{c}$/$\boldsymbol{r}$ is the best column/row that can be added to the pool of selected columns/rows. However, the obtained vector is not available in the given dataset as a column or row since it is a singular vector which is a function of all columns/rows. The following step re-imposes the underlying constraints at each iteration,

Here, $\mathbb{S}_{k}^{c}$ indicates a singleton that contains $k^{\text{th}}$ selected column and $\mathbb{S}_{k}^{r}$ corresponds to the $k^{\text{th}}$ selected row. In each iteration, one column or one row is the subject of optimization. The impact of the latest estimation for that column/row on the representation is neglected. Then, an optimized replacement is found. At each iteration of TWSP, sub-problems in (3) are solved and their solutions are matched to the accessible column samples and row samples through matching equations in (4). The new selected column or row is stored in $\mathbb{S}_{k}^{c}$ or $\mathbb{S}_{k}^{r}$, respectively. At each iteration we need to compute only the first singular vector and there are fast methods to do so [4]. The pair of (3a) and (4a) optimizes and matches a column. Similarly, performing (3b) and (4b) provides us with an optimized row. However, we do not update both of them per each iteration. In fact, we need to perform a column update or a row update in each iteration. It should be determined that which update (column or row) is more efficient in each iteration of the algorithm. To this aim, first, we choose a random previously selected column and a random previously selected row. Then, we find the best possible replacement column and the best possible replacement row. Accordingly, we choose whichever that minimizes the cost function more. The best modified column-wise subset is denoted by $\tilde{\mathbb{S}}^{c}$ and $\tilde{\mathbb{S}}^{r}$ denotes the best row-wise modified subset. Alg. 1 indicates the steps of TWSP algorithm. Here, $\dagger$ refers to the Moore–Penrose pseudo-inverse operator. Iterations can be terminated either once CUR decomposition error is saturated or when a maximum number of iterations is reached.

The proposed TWSP provides the selected columns and selected rows in order to form matrix $\boldsymbol{C}$ and $\boldsymbol{R}$ in CUR decomposition. It is straightforward to estimate the core matrix $\boldsymbol{U}$. Mathematically,

This matrix is a two-way compressed replica of the whole dataset and it contains valuable information in practice as will be discussed in Sec. 4.2. In general, the number of selected columns may differ from the desired number of rows. Here, $K_{1}$ refers to the number of columns and $K_{2}$ points to the number of rows. It is worthwhile to mention that the complexity order of TWSP is bottle-necked by computational burden for two pseudo-inverses and two singular vectors computation. Thus, the complexity can be expressed as $\small{O(NK_{1}^{2}\!+\!MK_{2}^{2}\!+\!MN)}$ per iteration and the algorithm needs $\small{O(\text{max}(K_{1},K_{2}))}$ iterations. Moreover, operator rnd($K$) refers to an integer random number less than or equal to $K$. In the next section, we evaluate the performance of our proposed algorithm.

The CSSP and CUR decomposition problems are NP-hard, i.e., a combinatorial search is required to find the best columns and rows. The proposed TWSP algorithm minimizes the main cost function (1) in practice. However, there is no theoretical guarantee for convergence of the proposed TWSP algorithm Alg. 1. In order to improve the convergence behavior of TWSP, we evade updating both columns and rows in each iteration. Rather, we prioritize updating a row or a column depending on which exhibits a smaller projection error, i.e., which is a better minimizer for the cost function per each iteration. The implementation steps of TWSP are summarized in Alg. 1. As seen in Alg. 1, the method adds up a revision feature to remedy the greediness associated with consecutive selection by removing samples and optimizing to select a sample for replacement and reiterating this procedure.

In order to evaluate TWSP on machine-learning tasks in terms of CUR decomposition accuracy, we apply the proposed TWSP on synthetic data as well as three real applications.

Two-way spectrum pursuit is proposed as an accurate and efficient algorithm for solving CUR decomposition. TWSP can be employed for joint selection of columns and rows such that their outer products is able to reconstruct the whole matrix as accurate as possible. Some applications of the proposed algorithm are presented to show the efficacy of the proposed method. However, they are not limited to the mentioned applications. Moreover, the proposed algorithm can be extended to $n$-way spectrum pursuit for efficient tensor subset selection.