Linear Time Approximation Algorithm for Column Subset Selection with Local Search

Authors: yuanbin  zou
ziyun  huang
jinhui  xu
jianxin  wang
qilong  feng
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 16:09:20

Abstract: The Column Subset Selection (CSS) problem has been widely studied in dimensionality reduction and feature selection. The goal of the CSS problem is to output a submatrix S, consisting of k columns from an n×d input matrix A that minimizes the residual error ‖A-SS^\dagger A‖_F^2, where S^\dagger is the Moore-Penrose inverse matrix of S. Many previous approximation algorithms have non-linear running times in both n and d, while the existing linear-time algorithms have a relatively larger approximation ratios. Additionally, the local search algorithms in existing results for solving the CSS problem are heuristic. To achieve linear running time while maintaining better approximation using a local search strategy, we propose a local search-based approximation algorithm for the CSS problem with exactly k columns selected. A key challenge in achieving linear running time with the local search strategy is how to avoid exhaustive enumerations of candidate columns for constructing swap pairs in each local search step. To address this issue, we propose a two-step mixed sampling method that reduces the number of enumerations for swap pair construction from O(dk) to k in linear time. Although the two-step mixed sampling method reduces the search space of local search strategy, bounding the residual error after swaps is a non-trivial task. To estimate the changes in residual error after swaps, we propose a matched swap pair construction method to bound the approximation loss, ensuring a constant probability of loss reduction in each local search step. In expectation, these techniques enable us to obtain the local search algorithm for the CSS problem with theoretical guarantees, where a 53(k+1)-approximate solution can be obtained in linear running time O(ndk^4\log k). Empirical experiments show that our proposed algorithm achieves better quality and time compared to previous algorithms on both small and large datasets. Moreover, it is at least 10 times faster than state-of-the-art algorithms across all large-scale datasets.

Main paper: 
 introduction: 
 In machine learning, handling high-dimensional datasets often requires the use of dimensionality reduction techniques, among which Singular Value Decomposition (SVD) is one of the most commonly utilized techniques in practice. The goal of SVD is to minimize the Frobenius norm of the error, aiming to achieve a low-rank approximation of a matrix with theoretical guarantees. An alternative way is to select a small subset of columns from the matrix as representations to well approximate the given matrix, which is known as the Column Subset Selection (CSS) problem. The CSS problem has been widely studied in machine learning for dimensionality reduction with improved interpretability. As pointed out in [4, 27, 17, 13, 20] , the number of selected columns k is much smaller than both n and d, i.e., k ≪ min{n, d}. We consider the following CSS problem. Definition 1.1. Given a matrix A ∈ R n×d and a positive integer k, the goal of CSS problem is to select k columns of A forming a matrix S ∈ R n×k that minimizes the residual error∥A -SS † A∥ 2 F ,where S † represents the Moore-Penrose inverse matrix of S, and∥A∥ 2 F = n i=1 d j=1 A 2 ij denotes the square of Frobenius norm.The CSS problem is known to be UG-hard [9] . For the CSS problem, several heuristic algorithms [2, 22, 23] were proposed. However, a main concern for these heuristics is the lack of theoretical guarantees on both the running time and approximation error. Several (1 + ϵ)-approximation bi-criteria algorithms have been proposed for the CSS problem [19, 4, 13, 15] . Many of these algorithms achieve running time of O(nd • poly(k)), which is referred to as linear in both n and d [26, 11] . Although these algorithms achieve (1 + ϵ)-approximation, they require selecting more than k columns. For the CSS problem with exactly k columns selected, Boutsidis et al. [5] proposed an O(k 2 log k)-approximation algorithm with O(min{n 2 d, nd 2 }) running time, using leverage score sampling and QR decomposition methods. Deshpande and Rademacher [13] proposed a volume sampling algorithm, which yields a (k + 1)-approximation in time O(kdn 3 log n). Guruswami and Sinop [19] gave an improved approximation algorithm with O(n 2 dk) running time using fast volume sampling method, achieving the same (k + 1)-approximation. The running time of these algorithms has at least a quadratic dependence on n or d. Deshpande and Vempala [15] presented a linear-time algorithm with (k + 1)!-approximation using adaptive sampling method.Although the algorithm using the adaptive sampling method [15] achieves the linear running time, its approximation ratio is considerably larger than other algorithms. Moreover, as shown in [23] , local search can improve the quality of the solution for the CSS problem. However, a potential limitation of the local search algorithm proposed in [23] is the lack of theoretical guarantees on the number of local search steps required for reaching a convergence.To apply the local search strategy to handle the CSS problem, the running time of each iteration is quadratically dependent on d, making it impractical for large-scale datasets. More specifically, the single-swap strategy in local search enumerates O(kd) swaps to improve the current solution, where the swap pair is constructed between the given matrix and the set of selected columns during the local search step. This method results in an O(nd 2 k 2 ) running time in each iteration, making it difficult to maintain linear running time in both n and d. Secondly, to the best of our knowledge, there is no available result that provides an approximation guarantee for solving the CSS problem using the local search strategy.In this paper, we propose a local search algorithm for solving the CSS problem with running time linear in both n and d. The key challenge for the local search algorithm is to avoid the O(nd 2 k 2 ) running time caused by enumerating all possible swap pairs. To overcome this challenge, we propose a two-step mixed sampling method that selects a candidate column with a specific probability for swapping, reducing the enumeration of swap pairs from O(dk) to k. By applying the two-step mixed sampling method, the running time of each local search step is reduced from O(nd 2 k 2 ) to O(ndk 2 ).Although the two-step mixed sampling method accelerates the local search process, bounding the residual error of the solution after swaps is a non-trivial task. Specifically, it is challenging to theoretically estimate the improvement in the approximation loss of the current solution after swaps.To address this issue, we propose a matched swap-pair construction method. This method identifies matched column pairs between the current solution and the optimal solution. Based on these matched pairs, we guarantee that the two-step mixed sampling method can find a column. With constant probability, this column reduces the residual error of the updated solution by a multiplicative factor of 1 -Θ(1/k) in each local search step. With this approach, the expected number of iterations can be bounded by O(k2 log k). Therefore, we obtain a local search-based approximation algorithm with O(ndk 4 log k) running time. The main contributions of this paper are summarized as follows.• For the CSS problem, we propose a new algorithm that uses local search with a two-step mixed sampling method. This method avoids the quadratic dependence of d on running time by reducing the number of swap pair enumerations during the local search process. Additionally, we propose a matched swap-pair construction method to bound the improvement of residual error during swaps. With these techniques, we achieve a 53(k + 1)approximation with exactly k columns selected, where the running time of our proposed algorithm is O(ndk 4 log k).• Numerical experiments show that our algorithm performs better in terms of quality on both small and large datasets compared to previous algorithms that selects exactly k columns, and is at least 10 times faster than the state-of-the-art algorithms on all large datasets.Table 1 : Comparison of existing results for the CSS Problem with exactly k columns selected, where n is the number of rows in the given matrix, d is the number of columns, and k is the number of selected columns.References Approximation Ratio Method Running Time [15] (k + 1)! adaptive sampling O(ndk 2 ) [5] k 2 log k leverage score sampling + QR decomposition O(min{n 2 d, nd 2 }) [13] k + 1 volume sampling O(kn 3 d log n) [19] k + 1 fast volume sampling O(n 2 dk) This paper 53(k + 1) local search + two-step mixed sampling O(ndk 4 log k)Within the framework of rank-revealing QR factorization (RRQR) [6] , several poly(k, d)approximation results [3, 7, 18] have been proposed. These results achieve running time of O(nd 2 ) while selecting exactly k columns to solve the CSS problem with Frobenius norm error. Deshpande et al. [14] gave the lower bound of (k + 1)-approximation for the problem. Furthermore, Boutsidis et al. [5] proposed a randomized algorithm with O(min{n 2 d, nd 2 }) running time and O(k 2 log k)approximation. More precisely, Deshpande and Rademacher [13] provided a deterministic (k + 1)approximation with O(dn 3 log nk) running time. To improve the running time, Guruswami and Sinop [19] proposed an O(n 2 dk) time randomized algorithm with (k + 1)-approximation. Additionally, there are many bi-criteria algorithms that relax the number of selected columns. Volume sampling methods have been widely applied to the CSS problem. Deshpande and Vempala [15] utilized these methods to achieve a PTAS, selecting O(k/ϵ 2 + k 2 log k) columns. Boutsidis et al. [4] proposed a linear-time algorithm with a (1 + ϵ)-approximation, requiring O(k/ϵ) columns selected. Guruswami and Sinop [19] developed a deterministic algorithm that also achieves a (1+ϵ)-approximation with O(k/ϵ) columns selected. Civril and Magdon-Ismail [10] gave improved bounds for obtaining a PTAS using k-leverage score sampling and SVD, by selecting O(k log k/ϵ 2 ) 2 columns. Altschuler et al. [1] developed a distributed greedy algorithm for the objective ∥SS † A∥ 2 F . Wang and Singh [28] studied the CSS problem in the missing-data case. Several bi-criteria algorithms have been proposed for the CSS problems with ℓ 1 norm [12, 21, 25, 8] . In both offline and online settings, Woodruff and Yasuda [29] provided several bi-criteria algorithms for the CSS problem with the ℓ p norm. if t is equal to 1 thenInitialize an n × d zero matrix D, and set each diagonal entryD ii = ∥A-SS † A∥ F (52 min{n,d}(k+1)!) 1/2 .A ← A + D. ∈ I such that f (A ′ , A ′ I\{q}∪{p} ) < f (A ′ , S) then 6:Find an index q ∈ I that minimizes f (A ′ , A ′ I\{q}∪{p} ).7:I = I\{q} ∪ {p}. 8: end if 9: return A ′ I . 
preliminaries: 
 For any positive integer n, let [n] denote the set {1, 2, . . . , n}. Given a matrix A ∈ R n×d , let A ij be the element in the i-th row and the j-th column of A, and define the Frobenius norm of A as∥A∥ 2 F = n i=1 d j=1 A 2 ij .Denote A :j as the j-th column of A, and A i: as the i-th row of A. Let A ⊤ be the transpose of A and A † be the Moore-Penrose inverse of A. Given an n × d matrix A, let I be the set of column indices from A, and let A I denote the n × |I| submatrix of A consisting of the columns corresponding to the indices in I. For a matrix A, the linear span of its column vectors is denoted as span (A). For any two n × d matrices A and B, ∥AB∥ F ≤ ∥A∥ F ∥B∥ 2 and ∥AB∥ F ≤ ∥A∥ 2 ∥B∥ F . Given any matrix A ∈ R n×d , the singular value decomposition (SVD) of A can be written as A = n i=1 σ i u i v ⊤ i ,A k = k i=1 σ i u i v ⊤i as the best rank-k approximation to A under the Frobenius norm. The spectral norm of A, denoted by ∥A∥ 2 , is defined as the largest singular value of A, i.e., ∥A∥ 2 = σ max (A). Given a solution S to the CSS problem on matrix A, we define the residual error of S as f (A, S) = ∥A -SS † A∥ 2 F . The LSCSS algorithm mainly comprises local search and two-stage mixed sampling components. The high-level idea behind our proposed local search is to identify a swap pair that minimizes the residual error in each iteration. The swap pair consists of a column from the input matrix to swap in and a column from the current solution to swap out. By repeating this process, the algorithm produces an updated solution with better quality. Moreover, the two-stage mixed sampling method involves two steps for obtaining a candidate column from the input matrix. Firstly, a set of column indices is constructed by sampling each column with probability proportional to its residual error for the current solution. Then, a column is uniformly selected from the set of candidate indices as the final column to swap in. To ensure that the input matrix for the local search process is full rank, we construct a new matrix A ′ by adding a small perturbation matrix D to the original matrix A during the initialization, where D is full-rank and has non-zero values only on its diagonal. The full-rank property of A ′ is used in subsequent analysis.The LSCSS algorithm begins by obtaining an initial solution S with exactly k columns and constructing a full-rank matrix A ′ during the initialization (steps 1-12 of Algorithm 1), which achieves a (k + 1)!-approximate solution on A ′ . We start by initializing an index set I and setting the matrix E = A. A new column index is added to I by sampling each column index i from [d] with probability proportional top i = ∥E :i ∥ 2 2 /∥E∥ 2 F . Then, E is updated as E = A -A I A † I A.Repeating this process k times, we obtain an initial solution S = A I . To construct a full-rank matrix, we construct an n×d zero matrix D and compute the parameter α = ∥A-SS † A∥ F /(52 min{n, d}(k +1)!) -1/2 using the initial solution S and A. Each diagonal entry D ii is set to α. Since rank(A+D) = rank(D), we construct the full-rank matrix A ′ by adding the full-rank matrix D to the input matrix A. To solve the CSS problem on A ′ , we execute steps 3-6 of Algorithm 1 to obtain the solution S = A ′ I . The detailed process described in steps 1-12 of Algorithm 1 requires O(ndk 2 ) time.The local search performed in steps 13-15 of Algorithm 1 plays a crucial role in LSCSS, involving two main steps. Firstly, we compute the matrix E = A ′ -SS † A ′ for the current solution S. Then, a set C of 10k column indices is constructed by sampling each column index i from [d] with probability p i = ∥E :i ∥ 2 2 /∥E∥ 2 F . Next, a column index p is uniformly selected from C, referred to as the swap-in column index. Let I denote the set of column indices of S in A ′ . Subsequently, if there exists an index q ∈ I such that f (A ′ , A ′ I\{q}∪{p} ) < f (A ′ , S), we choose q as the swap-out column index and update the set of indices to I = I\{q} ∪ {p}. Finally, Algorithm 2 returns the solution S = A ′ I . After repeating this process T = O(k 2 log k) times, Algorithm 1 returns the final solution S for the input matrix A.In the following, we explain in more detail how our proposed local search algorithm achieves a 53(k + 1)-approximation for the original matrix A. Given an initial solution, the main idea for analyzing the approximation ratio of our algorithm is to bound the improvement on residual error during the swaps in the local search step. To achieve this bound, we propose a matched swap pair construction method that guarantees an improvement in the current solution by swapping one column (Lemma 3.6 and Lemma 3.7). By carefully analyzing the improvement, we show that with constant probability the approximation loss of the current solution can be reduced by a multiplicative factor of 1 -Θ(1/k) in each iteration of the local search algorithm (Lemma 3.8). This implies that after O(k 2 log k) iterations, we have∥A ′ -SS † A ′ ∥ 2 F ≤ 26(k + 1)∥A ′ -A ′ k ∥ 2 F (Theorem 3.9).Finally, by analyzing the change in the residual error caused by removing matrix D from A ′ , we obtain∥A ′ -SS † A ′ ∥ 2 F ≤ 53(k + 1)∥A -A k ∥ 2F in expectation (Lemma 3.10). We assume that the matrix A has been normalized such that ∥A∥ 2 F = 1/4. Otherwise, we can normalize each element A ij in A as A ij = Aij 2∥A∥ F . Next, we consider a single iteration of Algorithm 2. We assume that the current solution has a high residual error (larger than 25(k + 1)∥A ′ -A ′ k ∥ 2 F before executing Algorithm 2 on A ′ . Otherwise, the initial solution S is a 25(k + 1)-approximation for the input matrix A ′ .Let S * = {s * 1 , . . . , s * k } be the optimal solution with exactly k selected for A ′ , and let S = {s 1 , . . . , s k } be the current solution. We define ϕ(A ′ , S * , S, s * ) = arg min s∈S f (A ′ , S * \{s * } ∪ {s}) as a mapping function that finds s from S such that the residual error f (A ′ , S * \{s * } ∪ {s}) is minimized. Thus, we say that s * is captured by ϕ(A ′ , S * , S, s * ). Each column s * ∈ S * is captured by exactly one column from S. Let I denote the set of column indices of S in matrix A ′ . We denote L as the set of columns indices in S that do not capture any optimal columns. We denote H as the set of indices where each column in S captures exactly one optimal column.The main idea behind the matched swap pair construction method is to analyze the change in residual error caused by swapping an index from set H (or L) with an index from a sampled column, using a two-step mixed sampling approach for the current solution S. For the column s h (where h ∈ H) in S, s h captures exactly the column s * h of the optimal solution S * , serving as the candidate column for s * h . If the residual error of swapping s h to replace s * h is large, we prove that with constant probability, sampling a new column can reduce the residual error and update s h . Similarly, for the column s l (where l ∈ L) in S, s l does not match any optimal column. We also show that, with constant probability, sampling a column from the input matrix A ′ can reduce the residual error for columns in set L. To analyze the improvement in residual error during swaps, we focus on a single swap process, evaluating both the increase in residual error from removing a column s from S and the decrease in residual error from inserting a new column. We give the following definition to measure the change resulting from removing a column. Definition 3.1. Let A ′ ∈ R n×d be a full-rank matrix, and let S be a solution on A ′ . Let I be the set of column indices of S. The change in residual error by removing the column i from I is defined asτ (A ′ , S, I\{i}) = f (A ′ , A ′ I\{i} ) -f (A ′ , S).To bound τ (A ′ , S, I\i) of solution S on the matrix A ′ , we provide the theoretical guarantee in the following lemma. (Detailed proof of Lemma 3.2 is given in Appendix A.1) Lemma 3.2. Let A ′ ∈ R n×d be a full-rank matrix, and let S be a solution on A ′ . Let I be the set of the column indices in S. For i ∈ I, we haveτ (A ′ , S, I\{i}) ≤ ∥A ′ I A ′ † I A ′ ∥ 2 F .To further analyze the bound on τ (A ′ , S, I\{i}), we decompose the projection matrix A ′ I A ′ † I A ′ and show that the expected upper bound of τ (A ′ , S, I\{i}) is proportional to ∥A ′ ∥ 2 F . (Detailed proof of Lemma 3.3 is given in Appendix A.1) Lemma 3.3. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let I be the set of column indices of S for the CSS problem on A ′ . In expectation, the following inequality holds:∥A ′ I A ′ † I A ′ ∥ 2 F ≤ k 2 d 2 ∥A ′ ∥ 2 F .In the following, we theoretically bound the residual error resulting from adding a candidate column index p to the set I of column indices in S, where p is chosen using the two-step mixed sampling method. (Detailed proof of Lemma 3.4 is given in Appendix A.1) Lemma 3.4. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be a solution with the set I of column indices in S. Let E = A ′ -SS † A ′ . The column index p is obtained by executing steps 2-3 of Algorithm 2. In expectation, the following inequality holdsf (A ′ , A ′ I∪{p} ) ≤ f k (A ′ , opt) + 1 10 f (A ′ , S),where f k (A ′ , opt) denotes the best rank-k solution.According to the aforementioned mapping function ϕ(•), we obtain the subset H from the set I of column indices and the set R = I\H. By using the matched swap pair construction method, there are two cases for the residual error of the current solution:1. For the set H, where h∈H f(A ′ , A ′ I\{h} ) > 21 50 i∈I f (A ′ , A ′ I\{i} ). 2. For the set R = I\H, where r∈R f (A ′ , A ′ I\{r} ) ≥ 29 50 i∈I f (A ′ , A ′ I\{i} ).By Lemma 3.2 and Lemma 3.4, we define the good columns s i for i ∈ I with respect to S as follows.Definition 3.5. Let A ′ ∈ R n×d be a full-rank matrix, and let k be a positive integer. Let S * be the optimal solution with exactly k columns selected, and let I * be the set of column indices in S * . Let S be any solution with exactly k columns selected, and let I be the set of column indices in S. A column index i ∈ I is called good iff (A ′ , A ′ I\{i} ) -τ (A ′ , S, I \ {i}) -τ (A ′ , A ′ I∪{p} , (I ∪ {p})\{q}) - 11 10 f (A ′ , A ′ I * \{i * } ) + 1 10 f (A ′ , S) > 1 100k f (A ′ , S).where i * ∈ I * is the column index mapped from i ∈ I by the function ϕ(•).Definition 3.5 estimates the gain from replacing the column s h with a new column obtained using the two-step sampling method. Next, we argue that if case (1) happens, the sum of residual errors for the good columns is large. (Detailed proof of Lemma 3.6 is given in Appendix A.1) Lemma 3.6. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be the solution to the CSS problem on A ′ . Let I be the set of column indices in S.If 50 h∈H f (A ′ , A ′ I\{h} ) ≥ 21 i∈I f (A ′ , A ′ I\{i} ) and f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt), we have h∈H,h is good f (A ′ , A ′ I\{h} ) ≥ 1 125 i∈I f (A ′ , A ′ I\{i} ).Since R = I\H, it holds that L ⊆ R. The index set R contains two subsets: L and R\L, where the indices in L do not capture any optimal columns according to the mapping function ϕ(•) and the indices in R\L capture at least two columns. Similar to case (1), we argue that if case (2) occurs, the sum of residual errors for the good columns is large. (Detailed proof of Lemma 3.7 is given in Appendix A.1) Lemma 3.7. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be a solution for the CSS problem on matrix A ′ . Let I be the set of column indices in S.If r∈R f (A ′ , A ′ I\{r} ) ≥ 29/50 i∈I f (A ′ , A ′ I\{i} ) and f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt), we have r∈R,r is good f (A ′ , A ′ I\{r} ) ≥ 1 125 i∈I f (A ′ , A ′ I\{i} ).In the following, we prove that if the residual error of the current solution is larger than 25(k + 1)f k (A ′ ,′ . Suppose that f (A ′ , S) ≥ 25(k + 1) • f k (A ′ , opt). Then, with probability at least 1/1375, Algorithm 2 returns a new solution S ′ with f (A ′ , S ′ ) ≤ (1 -1/(100k))f (A ′ , S).Subsequently, we prove that the LSCSS algorithm achieves a 26(k + 1)-approximation for A ′ after O(k 2 log k) iterations.Theorem 3.9. Let A ′ ∈ R n×d be the input matrix obtained in step 12 of Algorithm 1, let k be a positive integer, and let S be the solution returned after executing Algorithm 2 T = O(k 2 log k) times. Then, it holds thatE[∥A ′ -SS † A ′ ∥ 2 F ] ≤ 26(k + 1)∥A ′ -A ′ k ∥ 2 F , where A ′ k is the best rank-k approximation of A ′ for the CSS problem. The running time of Algorithm 1 is O(ndk 4 log k).Proof. Let Ŝ denote the submatrix consisting of k columns obtained in step 12 of Algorithm 1. For the initial solution Ŝ, Deshpande and Vempala [15] provide an approximation ratio of (k + 1)!. Before executing steps 13-15 of Algorithm 1, the residual error of the initial solution Ŝ is larger than25(k + 1)∥A ′ -A ′ k ∥ 2 F .According to Lemma 3.8, with probability 1/1375, we can reduce the residual error by a multiplicative factor of (1 -1/100k).Let T = O(k 2 log k). We define a random process P with initial residual error ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F of the solution Ŝ such that for T iterations of Algorithm 2, it reduces the value of ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F by at least (1 -1 100k) with probability 1/1375, and it increases the final value of∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F by 25(k + 1)∥A ′ -A ′ k ∥ 2 F . It is obvious that E[∥A ′ -SS † A ′ ∥ 2 F ] ≤ E[∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F ]. Then, we have In the following, we analyze the change in residual error caused by replacing the input matrix A with A ′ = A + D, which leads to the final solution of Algorithm 1 achieving a 53(k + 1)-approximation. (Detailed proof of Lemma 3.10 is given in Appendix A.1) Lemma 3.10. Let A ∈ R n×d be an input matrix, and let k be a positive integer. Define D as an n × d matrix with elementsE[P] =25(k + 1)∥A ′ -A ′ k ∥ 2 F + ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F • T i=0 T i i 1375 1374 1375 T -i 1 - 1 100k i ≤∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F • 1 - 1 137500k 137500k log(k+1)! + 25(k + 1)∥A ′ -A ′ k ∥ 2 F ≤ ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F (k + 1)! + 25(k + 1)∥A ′ -A ′ k ∥ 2 F . This implies that E[∥A ′ -SS † A ′ ∥ 2 F | Ŝ] ≤ ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F (k+1)! + 25(k + 1)∥A ′ -A ′ k ∥ 2 F . Thus, we obtain E[∥A ′ -SS † A ′ ∥ 2 F ] = Ŝ E[∥A ′ -SS † A ′ ∥ 2 F | Ŝ]P r( Ŝ) ≤ Ŝ P r( Ŝ) ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F (k + 1)! + 25(k + 1)∥A ′ -A ′ k ∥ 2 F ≤ E[∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F ] (k + 1)! + 25(k + 1)∥A ′ -A ′ k ∥ 2 F . Since ∥A ′ -Ŝ Ŝ † A ′ ∥ 2 F ≤ (k + 1)!∥A ′ -A ′ k ∥ 2 F in expectation, we have E[∥A ′ -SS † A ′ ∥ 2 F ] ≤ 26(k + 1)∥A ′ -A ′ k ∥ 2 F . Running Time Analysis. In LSCSSD ij = ∥A-S1S † 1 A∥ F (52 min{n,d}(k+1)!) 1/2 , if i = j 0, otherwise ,where S 1 is obtained by executing the first round of steps 3-6 in Algorithm 1. LetA ′ = A + D. The solution S 2 returned by executing Algorithm 2 for T = O(k 2 log k) iterations satisfies E[∥A ′ -S 2 S † 2 A ′ ∥ 2 F ] ≤ 53(k + 1)∥A -A k ∥ 2 F ,where A k is the best rank-k approximation for A.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [NA] Justification: This paper has no limitation.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach. For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
theory assumptions and proofs: 
 Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: All the theorems, formulas, and proofs are clearly stated in the paper. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems. • The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
experiments: 
 In this section, we compare our algorithm for the CSS problem with the previous ones. For hardware, all the experiments are conducted on a machine with 72 Intel Xeon Gold 6230 CPUs and 2TB memory.Datasets. In this paper, we evaluate the performance of our algorithms on a total of 22 real-world datasets. In previous studies [23, 1] , the CSS problem typically involves datasets with no more than 100,000 rows and 20,000 columns. We include the 14 smaller datasets listed in Table 5 (Appendix A.2). To extend the evaluation to larger datasets, we include 8 additional datasets detailed in Table 2 .Six datasets contain between 40,000 and 480,000 columns, and two contain 400,000 and 8 million rows, respectively. All datasets can be found on the website 345 .Algorithms and parameters. In our experimental evaluation, we consider the following five distinct algorithms:• TwoStage. This is a two-stage algorithm from [5] that combines leverage score sampling and rank-revealing QR factorization.• Greedy. This is an algorithm in [16, 1] , which uses greedy algorithm to generate solution.• VolumeSampling. This is an algorithm in [13] , which uses volume sampling method.• ILS. This is an algorithm in [23] , which uses heuristic local search method.• LSCSS. This is our algorithm given in Algorithm 1, which uses the two-step mixed sampling and local search methods.Methodology We use the error ratio to evaluate the effectiveness of various algorithms, as defined in [23] . The error ratio is given by the formula ∥A -SS † A∥ 2 F /∥A -A k ∥ 2 F, where it quantifies the discrepancy between the selected columns and the optimal rank-k matrix approximation. A smaller error ratio indicates better algorithm performance. Following [23] , we test the TwoStage, VolumeSampling, ILS, and LSCSS algorithms on each dataset 10 times to calculate the average error ratio and running time. Since the Greedy algorithm is deterministic, it is tested only once per dataset. Experimental setup. For the CSS problem with the Frobenius norm, we run the TwoStage, Greedy, ILS, and LSCSS algorithms on both the 8 large datasets and 14 small datasets, providing the average results for each method. The ILS and our LSCSS algorithm are based on local search method. For fair comparison, we set the number of iterations to be 2k for ILS and LSCSS. Since the VolumeSampling requires O(dkn 3 log n) runtime and O(n 2 + d 2 ) memory, it cannot handle the 8 large datasets because the algorithm requires more than 48 hours of runtime and over 2TB of memory. However, the other four algorithms generally produce a solution within 48 hours and with less than 2TB of memory. Thus, we only include VolumeSampling in the comparison on the 14 smaller datasets and exclude its results from Tables 3 and 4 .Results for the CSS problem. Table 3 shows the comparison of running time for varying values of k, where the time is measured by seconds. LSCSS is at least 10 times faster than other algorithms across all datasets and at least 15 times faster than the TwoStage and Greedy algorithms. Our algorithm successfully outputs a feasible solution within 5 hours on all datasets, whereas other algorithms fail to do so within 48 hours or need more than 2TB memory. The comparison of error ratios, reported as mean±std with the best results highlighted in bold, is presented in Table 4 in the Appendix A.2. The LSCSS algorithm achieves the best error ratios on almost all datasets.Moreover, we compare the running time and error ratio of five algorithms with varying values of k on 14 small datasets (Appendix A.2). The experimental results show that the LSCSS algorithm outperforms other algorithms in terms of quality and is at least 2 times faster than Greedy, Volume-Sampling and ILS algorithms on all small datasets.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: The information required to reproduce the main experimental results is provided in the paper. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
conclusion: 
 In this paper, we propose a linear-time approximation algorithm for the CSS problem using local search and two-step mixed sampling methods. Experimental results demonstrate that our framework outperforms previous algorithms for solving the CSS problem with exactly k columns selected. An interesting future direction is how to design multi-swap local search approximation algorithms for handling the CSS problem.A Appendix / supplemental material A.1 Missing Proofs Lemma 3.2. Let A ′ ∈ R n×d be a full-rank matrix, and let S be a solution on A ′ . Let I be the set of column indices in S. For i ∈ I, we haveτ (A ′ , S, I\{i}) ≤ ∥A ′ I A ′ † I A ′ ∥ 2 F .Proof. Given a solution S, we demonstrate an equivalent transformation to simplify the expression of the residual error for S:∥A ′ -SS † A ′ ∥ 2 F = tr((A ′ -SS † A ′ ) ⊤ (A ′ -SS † A ′ )) = tr(A ′ ⊤ A ′ -A ′ ⊤ SS † A ′ -A ′ ⊤ (SS † ) ⊤ A ′ + A ′ ⊤ (SS † ) ⊤ SS † A ′ ) = tr(A ′ ⊤ A ′ -2A ′ ⊤ SS † A ′ + A ′ ⊤ SS † A ′ ) = tr(A ′ ⊤ A ′ ) -tr(A ′ ⊤ SS † A ′ ) = tr(A ′ ⊤ A ′ ) -tr(A ′ ⊤ (SS † ) ⊤ SS † A ′ ) = tr(A ′ ⊤ A ′ ) -∥SS † A ′ ∥ 2 F ,where the third equality follows from S † SS † = S † and (SS † ) ⊤ = SS † .By the above equality, we havef (A ′ , S ′ ) -f (A ′ , S) = ∥A ′ -S I\{i} S † I\{i} A ′ ∥ 2 F -∥A ′ -S I S † I A ′ ∥ 2 F = ∥S I S † I A ′ ∥ 2 F -∥S I\{i} S † I\{i} A ′ ∥ 2 F ≤ ∥S I S † I A ′ ∥ 2 F .Lemma 3.3. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let I be the set of column indices of S for the CSS problem on A ′ . In expectation, the following inequality holds:∥A ′ I A ′ † I A ′ ∥ 2 F ≤ k 2 d 2 ∥A ′ ∥ 2 F .Proof. We begin by seeking to bound the Frobenius norm of the matrix product A ′ I A ′ † I A ′ . Using the submultiplicative property of the Frobenius norm, we get∥A ′ I A ′ † I A ′ ∥ 2 F ≤ ∥A ′ I A ′ † I ∥ 2 F • ∥A ′ ∥ 2 F .Next, by the norm inequality, we have∥A ′ I A ′ † I ∥ 2 F ≤ ∥A ′ † I ∥ 2 2 • ∥A ′ I ∥ 2 F .Since the matrix A is full rank, which is ensured by steps 7-10 of Algorithm 1. Therefore, ∥A ′ † I ∥ 2 2 can be reformulated as the maximum eigenvalue of A ′ ⊤ I A ′ I . Specifically, we have∥A ′ † I ∥ 2 2 = 1 λ min (A ′ ⊤ I A ′ I ) = λ max (A ′ ⊤ I A ′ I ).Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The datasets can be found via the link provided in the paper, and the codes are available upon request via email. Guidelines:• The answer NA means that paper does not include experiments requiring code. • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: For each experiment, the paper provides details on the computer resources required to reproduce the results.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: We present a linear time local search algorithm for the CSS problem with exactly k columns selected, which we expect to be beneficial for handling other dimensionality reduction tasks. The primary purpose is to offer algorithmic insights for obtaining linear time algorithm, where we do not foresee any societal consequences specifically related to the proposed method.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: This paper poses no such risks about safeguards. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [NA] Justification: This paper poses no such risks about safeguards. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: ′ I A ′ † I A ′ ∥ 2 F , we get ∥A ′ I A ′ † I A ′ ∥ 2 F ≤ λ max (A ′ ⊤ I A ′ I ) • ∥A ′ I ∥ 2 F • ∥A ′ ∥ 2 F .Since the maximum eigenvalue λ max (A ′ ⊤ I A ′ I ) is bounded above by ∥A ′ I ∥ 2 F , we obtain∥A ′ I A ′ † I A ′ ∥ 2 F ≤ ∥A ′ I ∥ 4 F • ∥A ′ ∥ 2 F .Finally, by the expectationE[∥A ′ I ∥ 2 F ] = k d ∥A ′ ∥ 2 F and ∥A∥ 2 F = 1 4 , we have ∥A ′ I ∥ 4 F ≤ k d 2 ∥A ′ ∥ 4 F = k 2 d 2 (∥A + D∥ 2 F ) 2 ≤ k 2 d 2 (2∥A∥ 2 F + 2∥D∥ 2 F ) 2 = k 2 d 2 (2∥A∥ 2 F + 2 min{n,d} i=1 D 2 ii ) 2 = k 2 d 2 (2∥A∥ 2 F + 2 min{n, d} ∥A -SS † A∥ 2 F 52 min{n, d}(k + 1)! ) 2 ≤ k 2 d 2 (2∥A∥ 2 F + 2 ∥A∥ 2 F 52(k + 1)! ) 2 ≤ k 2 d 2 (4∥A∥ 2 F ) 2 = k 2 d 2 .where the first equality results from A ′ = A+D, and the second inequality follows from the property∥A + B∥ 2 F ≤ 2(∥A∥ 2 F + ∥B∥ 2 F )for any two matrices A, B ∈ R n×d . Combining the above inequalities, we conclude that∥A ′ I A ′ † I A ′ ∥ 2 F ≤ k 2 d 2 ∥A ′ ∥ 2 F .Lemma 3.4. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be a solution with the set I of column indices in S. Let E = A ′ -SS † A ′ . The column index p is obtained by executing steps 2-3 of Algorithm 2. In expectation, the following inequality holdsf (A ′ , A ′ I∪{p} ) ≤ f k (A ′ , opt) + 1 10 f (A ′ , S),where f k (A ′ , opt) denotes the best rank-k solution.Proof.Let p j = ∥E:j ∥ 2 F ∥E∥ 2 F be the probability that each column index j ∈ [d] is sampled. To bound f (A, A T ∪{p} ), we define a random variable X l i for i ∈ [r], j ∈ [d],and l ∈ [z] as follows:X l i = v ij p j E :j with probability p j .We denote X i as the random variable that is randomly picked from the set C = {X 1 i , X 2 i , . . . , X l i }. Therefore, the expected value of X i is given byE(X i ) = z l=1 Pr(l-th is picked) • E X l i = 1 z z l=1 E(X l i ) = 1 z z l=1 d j=1 p j • v ij p j E :j = 1 z z l=1 d j=1 v ij E :j = 1 z z l=1 Ev i = Ev i ,where the second equality follows from the linearity of expectation, and the last second equality results from the linear combination of v i and E.Letw i = (SS † A)v i + X i for i ∈ [r]. Therefore, we have E(w i ) = σ i u i .Thus, we have the equalityw i -σ i v i = X i -Ev i .By calculating its second central moment, we obtainE(∥w i -σ i u i ∥ 2 2 ) = E(∥X i -Ev i ∥ 2 2 ) = E(∥X i ∥ 2 2 ) -2E(X i ) • Ev i + ∥Ev i ∥ 2 2 = E(∥X i ∥ 2 2 ) -∥Ev i ∥ 2 2 . Next, we seek to bound E(∥X i ∥ 2 2 ), E(∥X i ∥ 2 2 ) = E ∥ z l=1 Pr(l-th is picked) • X l i ∥ 2 = z l=1 (Pr(l-th is picked)) 2 E(∥X l i ∥ 2 2 ) + 1≤l1<l2≤z Pr(l 1 -th and l 2 -th are picked)E(X l1 i • X l2 i ) = 1 z 2 z l=1 E(∥X l i ∥ 2 2 ) + 2 z 2 1≤l1<l2≤z E(X l1 i • X l2 i ) = 1 z 2 z l=1 E(∥X l i ∥ 2 2 ) + z -1 z ∥Ev i ∥ 2 2 .The term of 2z 2 1≤l1<l2≤z E(X l1 i • X l2 i ) = z-1 z ∥Ev i ∥ 2 2follows from the independence of X l1 i and X l2 i . Therefore, we conclude thatE(∥w i -σ i u i ∥ 2 2 ) = 1 z 2 z l=1 E(∥X l i ∥ 2 2 ) - 1 z ∥Ev i ∥ 2 2 .For the first term in the above equation, we haveE(∥X l i ∥ 2 2 ) = d j=1 p j ∥E :j v ij ∥ 2 2 p 2 j = d j=1 ∥E∥ 2 F ∥E :j ∥ 2 F ∥E :j v ij ∥ 2 2 = ∥E∥ 2 F • d i=1 ∥E :j v ij ∥ 2 F ∥E :j ∥ 2 F ≤ ∥E∥ 2 F • ∥Ev i ∥ 2 F ∥E∥ 2 F = ∥Ev i ∥ 2 F .Combining the above inequality, we haveE(∥w i -σ i u i ∥ 2 2 ) ≤ 1 z ∥Ev i ∥ 2 F . Let y i = w i /σ i for i ∈ [k] and let matrix F = ( k i=1 y i u T i )A ′ . Therefore, we have ∥A ′ -A ′ I∪{p} A ′ † I∪{p} A ′ ∥ 2 F ≤ ∥A ′ -F ∥ 2 F. By decomposing F along the right singular vectors {v 1 , . . . , v d }, we have:E(∥A ′ -A ′ I∪{p} A ′ † I∪{p} A ′ ∥ 2 F ) ≤ E(∥A ′ -F ∥ 2 F ) = d i=1 E(∥(A ′ -F )v i ∥ 2 F ) ≤ k i=1 E(∥σ i u i -w i ∥ 2 2 ) + d i=k+1 σ 2 i ≤ f k (A ′ , opt) + k z f (A ′ , A ′ I ),where the last second inequality uses thatf k (A ′ , opt) = d i=k+1 σ 2 i .Setting the z = 10k, we obtain thatf (A ′ , A ′ I∪{p} ) ≤ f k (A ′ , opt) + 1 10 f (A ′ , S).Lemma 3.6. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be the solution to the CSS problem on A ′ . Let I be the set of column indices in S.If 50 h∈H f (A ′ , A ′ I\{h} ) ≥ 21 i∈I f (A ′ , A ′ I\{i} ) and f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt), we have h∈H,h is good f (A ′ , A ′ I\{h} ) ≥ 1 125 i∈I f (A ′ , A ′ I\{i} ).Proof. First, we prove that the residual error is non-increasing. Let S k be a subset of k columns from A ′ . Let p ∈ A ′ be a column not in S k , and define S k+1 = S k ∪ {p}. The residual error function is defined as:f (A ′ , S) = ∥A ′ -SS † A ′ ∥ 2 F = ∥A ′ ∥ 2 F -∥SS † A ′ ∥ 2 F .Our goal is to show that f (A ′ , S k+1 ) ≤ f (A ′ , S k ), i.e., the residual error is non-increasing when a column is added to S k .Since ∥A ′ ∥ 2 F is constant, we only need to prove that ∥S k+1 S † k+1 A ′ ∥ 2 F ≥ ∥S k S † k A ′ ∥ 2 F . For any column u of A ′ , denote ∆(u) = ∥S k+1 S † k+1 u∥ 2 2 -∥S k S † k u∥ 2 2 .If p ∈ span(S k ), then span(S k+1 ) = span(S k ), and thus ∆(u) = 0. Otherwise, if p / ∈ span(S k ), define the orthogonal component of p with respect to span(S k ) asp ⊥ = p -S k S † k p, and normalize it as p ′ = p ⊥ ∥p ⊥ ∥2 .For any u, the projection onto span(S k+1 ) isS k+1 S † k+1 u = S k S † k u + ⟨u, p ′ ⟩p ′ , where ⟨u, p ′ ⟩ denotes the inner product. Therefore, ∥S k+1 S † k+1 u∥ 2 2 = ∥S k S † k u∥ 2 2 + |⟨u, p ′ ⟩| 2 holds, implying that ∆(u) = |⟨u, p ′ ⟩| 2 ≥ 0. Summing over all columns u of A ′ , we obtain ∥S k+1 S † k+1 A ′ ∥ 2 F - ∥S k S † k A ′ ∥ 2 F = u ∆(u) ≥ 0.Thus, the residual error function satisfiesf (A ′ , S k+1 ) = ∥A ′ ∥ 2 F -∥S k+1 S † k+1 A ′ ∥ 2 F ≤ ∥A ′ ∥ 2 F -∥S k S † k A ′ ∥ 2 F = f (A ′ , S k). This proves that adding a column to S k does not increase the residual error. By Lemma 3.2, we haveτ (A ′ , A ′ I∪{p} , I ∪ {p}\{h}) = f (A ′ , A ′ I∪{p}\{h} ) -f (A ′ , A ′ I∪{p} ) ≤ ∥A ′ I∪{p} A ′ † I∪{p} A ′ ∥ 2 F .Similarly, we haveτ (A ′ , S, I\{h}) = f (A ′ , A ′ I\{h} ) -f (A ′ , S) ≤ ∥SS † A ′ ∥ 2 F .Thus, by Lemma 3.3, we obtain τ(A ′ , A ′ I∪{p} , I ∪ {p}\{h}) + τ (A ′ , S, I\{h}) ≤ (k+1) 2 d 2 ∥A∥ 2 F + k 2 d 2 ∥A∥ 2 F ≤ 2(k+1) 2 d 2 ∥A∥ 2 F , where h ∈ H. Next, we have h∈H f (A ′ , S\{h}) ≥ 21 50 i∈I f (A ′ , A ′ I\{i} ). By Definition 3.5 and Lemma 3.2, we haveh∈H,h is not good f (A ′ , A ′ I\{h} ) ≤ h∈H τ (A ′ , S, I\{h}) + τ (A ′ , A ′ I∪{p} , I ∪ {p}\{h}) + 11k 10 (f (A ′ , A ′ I * \{h * } ) + 1 10 f (A ′ , S)) + 1 100 f (A ′ , S) ≤ 2k(k + 1) 2 d 2 ∥A ′ ∥ 2 F + 11k 10 ( k 2 d 2 ∥A ′ ∥ 2 F + f (A ′ , S * )) + 11k + 1 100 f (A ′ , S) ≤ 31k 3 10d 2 ∥A ′ ∥ 2 F + 11k 10 f (A ′ , S * ) + 11k + 1 100 f (A ′ , S) ≤ 31k 3 10d 2 ∥A ′ ∥ 2 F + 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S) ≤ 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S) + 31k 3 10d 2 ∥A∥ 2 F f (A ′ , opt) f k (A ′ , opt) ≤ 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S) + 31k 3 5d 2 f k (A ′ , opt),where the fourth inequality follows from the fact that S * is the optimal solution with exactly k columns selected, and f (A ′ , S * ) ≤ (k + 1)f k (A ′ , opt).Using f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt) and the non-increasing property of residual error, we obtain that h∈H,his not goodf (A ′ , A ′ I\{h} ) ≤ 103k 250 • f (A ′ , S). By kf (A ′ , S) ≤ i∈I f (A ′ , A ′ I\{i} ), h∈H,his not good f (A ′ , A ′ I\{h} ) ≤ 103 250 • i∈I f (A ′ , A ′ I\{i} ).Thus, we haveh∈H,h is good f (A ′ , A ′ I\{h} ) ≥ 1 125 • i∈I f (A ′ , A ′ I\{i} ).Lemma 3.7. Let A ′ ∈ R n×d be a full-rank matrix, k be a positive integer, and let S be a solution for the CSS problem on matrix A ′ . Let I be the set of column indices in S.If r∈R f (A ′ , A ′ I\{r} ) ≥ 29/50 i∈I f (A ′ , A ′ I\{i} ) and f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt), we have r∈R,r is good f (A ′ , A ′ I\{r} ) ≥ 1 125 i∈I f (A ′ , A ′ I\{i} ).′ , A ′ I∪{p} , I ∪ {p}\{l}) + τ (A ′ , S, I\{l}) ≤ (k+1) 2 d 2 ∥A∥ 2 F + k 2 d 2 ∥A∥ 2 F ≤ 2(k+1) 2 d 2 ∥A∥ 2 F , where l ∈ L. We have r∈R f (A ′ , A ′ I\{r} ) ≥ 29/50 i∈I f (A ′ , A ′ I\{i} ). Note that |R| ≤ 2|L|. By Definition 3.5 and Lemma 3.2, we have r∈R,r is not goodf (A ′ , A ′ I\{r} ) ≤ 2|L| min l∈L τ (A ′ , S, I\{l}) + τ (A ′ , A ′ I∪{p} , I ∪ {p}\{l}) + 11k 10 (f (A ′ , A ′ I * \{l * } ) + 1 10 f (A ′ , S)) + 1 100 f (A ′ , S) ≤ 2 l∈L τ (A ′ , S, I\{l}) + τ (A ′ , A ′ I∪{p} , I ∪ {p}\{l}) + 11k 10 (f (A ′ , A ′ I * \{l * } ) + 1 10 f (A ′ , S)) + 1 100 f (A ′ , S) ≤ 4k(k + 1) 2 d 2 ∥A ′ ∥ 2 F + 11k 10 ( k 2 d 2 ∥A ′ ∥ 2 F + f (A ′ , S * )) + 11k + 1 100 f (A ′ , S) ≤ 51k(k + 1) 2 10d 2 ∥A ′ ∥ 2 F + 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S) ≤ 51k 3 10d 2 ∥A ′ ∥ 2 F f k (A ′ , opt) f k (A ′ , opt) + 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S) ≤ 51k 3 5d 2 f k (A ′ , opt) + 11k(k + 1) 10 f k (A ′ , opt) + 11k + 1 100 f (A ′ , S).Using f (A ′ , S) ≥ 25(k + 1)f k (A ′ , opt) and the non-increasing property of the residual error (as in the proof of Lemma 3.6), we obtain r∈R,r is not goodf (A ′ , A ′ I\{l} ) ≤ 143k 250 f (A ′ , S) . By kf (A ′ , S) ≤ i∈I f (A ′ , A ′ I\{i} ), we have r∈R,r is not good f (A ′ , A ′ I\{l} ) ≤ 143 250 i∈I f (A ′ , A ′ I\{i} ).r∈R f (A ′ , A ′ I\{r} ) ≥ 29/50 i∈I f (A ′ , A ′ I\{i} ), we have r∈R,r is good f (A ′ , A ′ I\{r} ) ≥ 1 125 i∈I f (A ′ , A ′ I\{i} ).Lemma 3.8. Let A ′ ∈ R n×d and S ∈ R k×d be the input matrices for Algorithm 2, where k is a positive integer and S is the solution to the CSS problem on A ′ . Suppose that f (A ′ , S) ≥ 25(k + 1) • f k (A ′ , opt). Then, with probability at least 1/1375, Algorithm 2 returns a new solutionS ′ with f (A ′ , S ′ ) ≤ (1 -1/(100k))f (A ′ , S).Proof. Let I denote the set of column indices in S. Let S * be the optimal solution with exactly k columns selected, and let I * be set of the columns indices in S * . Let p denote the column index obtained by the steps 2-3 of Algorithm 2. Let p * ∈ I * be the column corresponding to index p, as determined by the mapping function ϕ(•). According to Lemma 3.6, the probability of sampling a good column index is at least(1/125) h∈H,h is good f (A ′ , A ′ I\{h} ) i∈I f (A ′ , A ′ I\{i} ) ≥ 1 125.Let h represent an arbitrary good column index. By Lemma 3.4 and the non-increasing property of residual error f (•) (as shown in the proof of Lemma 3.6), we haveE[f (A ′ , A ′ I∪{p} )] ≤ f k (A ′ , opt) + 1 10 f (A ′ , S) ≤ f (A ′ , A ′ I * \{h * } ) + 1 10 f (A ′ , S).Hence, by the Markov inequality, we obtainPr[f (A ′ , A ′ I∪{p} ) ≤ 11 10 (f (A ′ , A ′ I * \{h * } ) + 1 10 f (A ′ , S))|h ∈ H] ≥ 1 11 .Combining the above inequalities, the probability that the sampled column index p replaces some good column index h is at least1 125 • 1 11 > 1 1375.Similarly, according to Lemma 3.7, the probability of sampling a good column index is at least(1/125) r∈R,r is good f (A ′ , A ′ I\{r} ) i∈I f (A ′ , A ′ I\{i} ) ≥ 1 125.By Lemma 3.4 and the Markov inequality, we obtainPr[f (A ′ , A ′ I∪{p} ) ≤ 11 10 (f (A ′ , A ′ I * \{r * } ) + 1 10 f (A ′ , S))|r ∈ R] ≥ 1 11 .Thus, the probability that the sampled column index p for some good column index r is at least1 125 • 1 11 > 1 1375.′ , A ′ I∪{p} , I ∪ {p}\{q}) = f (A ′ , A ′ I∪{p}\{q} ) -f (A ′ , A ′ I∪{p} ), and τ (A ′ , S, I\{q}) = f (A ′ , A ′ I\{q} ) -f (A ′ , S). Now, for the column index q ∈ H ∪ R, we can upper bound the updated residual error f (A ′ , S ′ ) as follows.f (A ′ , S ′ ) ≤ f (A ′ , A ′ I\{q}∪{p} ) = f (A ′ , S) -f (A ′ , S) -f (A ′ , A ′ I\{q}∪{p} ) = f (A ′ , S) -f (A ′ , S) -f (A ′ , A ′ I∪{p} ) + f (A ′ , A ′ I∪{p} ) -f (A ′ , A ′ (I∪{p})\{q} ) = f (A ′ , S) -f (A ′ , S) -τ (A ′ , A ′ I∪{p} , (I ∪ {p})\{q}) -f (A ′ , A ′ I∪{p} ) = f (A ′ , S) -(f (A ′ , I\{q}) -τ (A ′ , S, I\{q}) -τ (A ′ , A ′ I∪{p} , (I ∪ {p})\{q}) -f (A ′ , A ′ I∪{p} )) ≤ f (A ′ , S) -(f (A ′ , I\{q}) -τ (A ′ , S, I\{q}) -τ (A ′ , A ′ I∪{p} , (I ∪ {p})\{q}) - 11 10 (f (A ′ , A ′ I * \{q * } ) + 1 10 f (A ′ , S))) ≤ f (A ′ , S) - 1 100k f (A ′ , S).Thus, by combining the two cases, we obtain that with probability at least 1/1375,f (A ′ , S ′ ) ≤ (1 -1/(100k))f (A, S ′ ).Lemma 3.10. Let A ∈ R n×d be an input matrix, and let k be a positive integer. Define D as an n × d matrix with elementsD ij = ∥A-S1S † 1 A∥ F (52 min{n,d}(k+1)!) 1/2 , if i = j 0, otherwise ,where S 1 is obtained by executing the first round of steps 3-6 in Algorithm 1. Let A ′ = A + D. The solution S 2 returned by executing Algorithm 2 for T = O(k 2 log k) iterations satisfiesE[∥A ′ -S 2 S † 2 A ′ ∥ 2 F ] ≤ 53(k + 1)∥A -A k ∥ 2 F, where A k is the best rank-k approximation for A.k be the best rank-k approximation of A ′ , and A k be the best rank-k approximation of A. For the initial solution S 1 , Deshpande and Vempala [15] provide an approximation ratio of (k + 1)!. By Lemma 3.9, we obtainE[∥A ′ -S 2 S † 2 A ′ ∥ 2 F ] ≤ 26(k + 1)∥A ′ -A ′ k ∥ 2 F ≤ 26(k + 1)∥A ′ -A k ∥ 2 F = 26(k + 1)∥A ′ -A + A -A k ∥ 2 F ≤ 52(k + 1) ∥A ′ -A∥ 2 F + ∥A -A k ∥ 2 F ≤ 52(k + 1) ∥D∥ 2 F + ∥A -A k ∥ 2 F = 52(k + 1)   ∥A -A k ∥ 2 F + min{n,d} i=1 D 2 ii   ≤ 52(k + 1) ∥A -A k ∥ 2 F + min{n, d} • ∥A -S 1 S † 1 A∥ 2 F 52 min{n, d} • (k + 1)! = 52(k + 1) ∥A -A k ∥ 2 F + ∥A -S 1 S † 1 A∥ 2 F 52(k + 1)! ≤ 52(k + 1) ∥A -A k ∥ 2 F + ∥A -A k ∥ 2 F 52 = 53(k + 1)∥A -A k ∥ 2 F , where the second inequality holds because ∥A ′ -A ′ k ∥ 2 F ≤ ∥A ′ -B∥ 2F for any n × d matrix B with rank(B) = k, and the third inequality results from the triangle inequality ∥A + B∥ 2 F ≤ 2(∥A∥ 2 F + ∥B∥ 2 F ) for any two matrices A, B ∈ R n×d .In this section, we compare our algorithm with four algorithms (TwoStage, Greedy, VolumeSampling, and ILS) introduced in Section 4 using 14 small real-world datasets. Following [23, 1] , these datasets are listed in Table 5 and all datasets can be found on the website 6 .We use the same experimental settings as in Section 4 to run these algorithms on 14 datasets. The running time results are presented in Tables 6 and 7 . LSCSS is at least 2 times faster than the Greedy, VolumeSampling, and ILS algorithms. The error ratio, reported as mean ± std, with the best results highlighted in bold, is shown in Tables 8 and 9 . LSCSS outperforms all other algorithms in terms of quality across all datasets.In this section, we present experimental results comparing the performance of our proposed LSCSS algorithm with the QR with Column Pivoting (QRP) algorithm [24] for varying k. The QRP algorithm uses column pivoting to improve the traditional QR decomposition.To compare our proposed algorithm with the QRP algorithm, we conducted experiments on five datasets: CMHS, ELD, Gas, FAds, and TGas (listed in Table 2 ). Due to the memory requirement of the QRP algorithm exceeding 2TB for the remaining datasets in Table 2 , those datasets were excluded from our comparison. For the QRP algorithm, we followed the procedure outlined in [24] , first obtaining the permutation matrix P that satisfies AP = QR. Then, we selected the top k elements from the diagonal of P as the indices C of the k columns and computed the error ratio of solution A C . The detailed results in Table 10 show that our algorithm achieves lower error ratios across all datasets and is faster than the QRP algorithm. Justification: The main claims made in the abstract reflect the papers primary contributions and scope. In this paper, an approximation algorithm with running time linear in both n and d is proposed for the Column Subset Selection problem, which selects exactly k columns.• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Justification: The research presented in the paper complies with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).Õ(n) denotes the asymptotic complexity, ignoring polylogarithmic factors, i.e., terms of the form poly(log n).https://archive.ics.uci.edu/datasetshttps://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html.http://www.cs.columbia.edu/CAVE/software/softlib/coil-20.php., http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html., http://www.sheffield.ac.uk/eee/research/iel/research/face., http://www.cs.nyu.edu/~roweis/data.html., http://vision.ucsd.edu/~leekc/ExtYaleDatabase/ExtYaleB.html., http://www.iro.umontreal.ca/~lisa/twiki/bin/view.cgi/Public/PublicDatasets., https://archive.ics.uci.edu/datasets, https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/, and http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html. 


Decision: This paper provides a new computationally efficient local search method for column subset selection, a central problem in machine learning and randomized numerical linear algebra. It provides an O(k) factor approximation using exactly k column, improving on prior work that either returned more columns (i.e., gave a bicriteria approximation) or obtained a much worse approximation factor. All reviewers felt the paper made progress on this longstanding problem, so should be accepted to the conference. We recommend that the reviewers strengthen their experimental evaluation by  including the comparison to pivoted QR that they ran during the rebuttal period. We would also recommend comparing against the recent “Randomly pivoted Cholesky” method by Chen et al.

Reviews: review: paper_summary: paper_summary: The paper addresses the classical CSS problem, and describes significant
improvements over the current state of the art.
The main idea is to run the randomized
adaptive sampling algorithm for selecting k columns,
and follow it by several iterations of swapping (local search),
that further increase the accuracy. The swapping replaces one of the selected
columns with another unselected column when such swapping reduces the error.
strength_weakness: strength_weakness: * There is no experimental comparison with the QRP (QR with Pivoting) method, 
  which is known to show very good results in practice.
  It is significantly faster than your algorithm, as it takes:

  4knd - 2 k^2(d+n) + 4 k^3/3

* The running time of the proposed algorithm is 

O(nd k^3 log k)

I do not consider it fast. There are other algorithms that run much faster
(eg Frieze, Kannan and Vempala, JACM 2004, and the above mentioned QRP).

I am also not that impressed by the relative bound of 110(k+1),
as there are other algorithms with 100 times better bounds (specifically k+1).
What is most surprising here is your experimental results.


Some notes:

line 41. CSS is known to be NP-hard. 
See "Column subset selection is NP-complete" by Shitov, 2021.

line 42. The comment is incorrect wrt reference 2.
questions: questions: * Why is there no comparison with the QRP?
* How good is your accuracy compared to the optimal of worst case of (k+1)?
limitations: limitations: N/A
paper_summary: paper_summary: This paper studies column subset selection. Given a matrix A n*d, how to select a matrix A_S of k columns, which preserves the substance in A? That is, reconstruction error ||A - (A_S A_S#) A||_F should be minimized, where A_S# is the psuedo-inverse of A_S. There are many approximation algorithms for this problem, and it is widely studied in theoretical ML/CS literature. There is a plethora of running-time vs approximation factor guarantee in the literature.

This paper tries to match the SoTA approx. factor with faster running time. Indeed, SoTA algorithm has (k+1)-factor approximation with running time O(n d^2 k). If we wanted linear in n and d, best known algorithm had approx. factor (k+1)! with O(n d k^2) running time. This paper's main contribution is 110 (k+1)-approximation with running time O( n d k^3).

Technically, the main innovation is to show that a constrained version of local search works. Local search is a good algorithm where we initialize a decent solution say the (k+1)! approx. solution, and iteratively swap out a column in the current solution with one not in. This step is expensive, having k*d possible combinations and we need to see which is the best! Multiplying this with # of rounds will yield large running times not-linear in d. 

This paper overcomes this by showng that not all swaps are needed to be tried. Indeed, they sample O(10k) possible columns based on some form of random sampling (by looking at the residual importance of various columns), and further uniformly sample one of them at random to swap in. Now, the complexity of testing swaps is only O(k)! The challenge is now to show that the algorithm will not get stuck in bad local optima.

For this authors present a set of "test swaps" based on optimal set of columns and algo set of columns, and show that we will end up sampling one such pair with decent chance, if the current soln is far from optimal.

Finally they run experiments comparing with prior art and show better running time AND quality. Overall solid piece of work and worthy of NeurIPS. The randomized local search idea is good, and the backing with empirical evals is also welcome.
strength_weakness: strength_weakness: The problem is important, though it would be nice for authors to present why "subset" selection is significant as opposed to top-k-SVD. It is a bit unfair to expect a full motivation because this is standard problem, but a bit more explanation for the importance of the problem would be nice.
questions: questions: Are there practical use-cases of the CSS (where the interpretabiltiy is crucial? as opposed to SVD). 
In typical use-cases, do we sample rows (from n like a coreset) or columns (from d like dimension reduction)? 

The obj. function f doesn't appear to be defined where it is first used.
Can your algo give better guarantees when allowed to pick more than k columns? How does this compare to prior SoTA in this regime?

Where is the randomness captured in Thm 3.9? Is it in expectation? WHP? The proof has some details but the formal statement should also have it.
limitations: limitations: Yes
paper_summary: paper_summary: The Column Subset Selection (CSS) problem aims to select a sub-matrix with $k$ columns from a matrix to minimize the residual error. Previous algorithms often have non-linear running times or higher approximation ratios. This paper proposes a local search-based algorithm with linear running time, utilizing a two-step mixed sampling method to reduce enumerations and a matched swap pair construction to maintain approximation quality. Empirical experiments demonstrate superior performance in quality and speed compared to existing methods, particularly in large datasets.
strength_weakness: strength_weakness: - Line 41: What is UG?

- Line 132: Is Lemma 2.1 necessary for the main text?

- Algorithms 1 and 2: It would be better to give the full name of the algorithms.

- Section 3 does not clearly discuss the algorithmic distinction beyond ILS and Greedy.

- The theoretical results of the algorithm are mixed with the algorithmic details, making the paper poorly structured. Additionally, the assumptions are not explicitly listed and lack sufficient discussion.

- Table 4: The standard error of Greedy hasn’t been reported. Besides, it would be better to put Table 4 into the main text and discuss why LSCSS achieves superior approximation compared to the competitors.

- Lines 320 and 321 are repetitive.
questions: questions: - Algorithm 2 (line 2): How would the numerical performance change if $10k$ columns were replaced with a different number of columns (e.g., $2k$ or even a constant [1])?

- Using the objective function to guide which column to swap out shares a similar spirit with recent advanced algorithms for best-subset selection under linear regression for sparsity-constrained optimization [2, 3]. It would be interesting to discuss these works. 

- How should $k$ be selected in practice?

- Table 9: The error ratio of TwoStage on the ComCri dataset is weird. It is exceptionally large when $k=50$ but suddenly reduces to a comparable result when $k=100$.

### Reference

[1] Bahmani, Sohail, Bhiksha Raj, and Petros T. Boufounos. "Greedy sparsity-constrained optimization." The Journal of Machine Learning Research 14.1 (2013): 807-841.

[2] Wang, Zezhi, et al. "Sparsity-Constraint Optimization via Splicing Iteration." arXiv preprint arXiv:2406.12017 (2024).

[3] Zhu, Junxian, et al. "A polynomial algorithm for best-subset selection problem." Proceedings of the National Academy of Sciences 117.52 (2020): 33117-33123.
limitations: limitations: The authors have discussed the limitations in the end of this paper.
paper_summary: paper_summary: This paper proposes a new algorithm for the column subset selection problem which combines a local search-type strategy with adaptive sampling techniques to obtain an algorithm running in time linear in $nd$ (i.e. the size of the input matrix) for constant $k$. The resulting solution selects exactly $k$ columns and achieves an approximation ratio at most $O(k)$ against the optimal rank $k$ solution given by the SVD, which is comparable to the best known approximation ratio of $(k+1)$. The main idea is that each iteration of the local search can be made to run in linear time by reducing the number of candidates considered via adaptive sampling.
strength_weakness: strength_weakness: The analysis of the algorithm is highly reminiscent of the adaptive sampling analyses from works such as

https://faculty.cc.gatech.edu/~vempala/papers/relative.pdf 

https://dl.acm.org/doi/10.1145/1250790.1250884 

This includes, for example, the idea of relating the cost difference from the optimum to a bound on the success probability of adaptive sampling (Lemma 3.8) and then boosting with repetition. The only new point here seems to be to use similar ideas to throw away “bad” columns (i.e. use local search) so that the solution size stays at $k$ rather than growing with each iteration. However, this is still a nice and new idea and requires delicate work to make this go through.
questions: questions: Questions

- For the problem of row sampling for $\ell_p$ subspace approximation, obtaining tight bounds on the number of rows required for a (1+eps) approximation is an open problem under active research (see, e.g., https://dl.acm.org/doi/10.1145/1250790.1250884). Can the combination of local search and adaptive sampling be used to obtain state of the art results for this problem?

Minor comments

- It would be helpful to move the main theorem Theorem 3.9 earlier so that the formal main result can be seen earlier. In particular, the notion of approximation ratio should be clarified before Table 1, since otherwise it is confusing whether the approximation ratio compares against the best rank k solution or the best subset of k columns.
- The wording of Lemmas 3.7 and 3.8 could be unified for sake of consistency.
- You can remove the section 7 header.
limitations: limitations: The authors have addressed the limitations.
paper_summary: paper_summary: This paper considers the Column Subset Selection (CSS) problem. In this problem the input is an arbitrary $A\in \mathbb{R}^{n\times d}$ and a positive integer $k$ ($k$ is thought to be much smaller than $\min\{n,d\}$). The goal is to output a subset of $k$ columns of $A$ denoted by $S\in \mathbb{R}^{n\times k}$ such that the residual error $\|A- SS^{\dagger}\|_F^2$ is minimized (here $S^{\dagger}$ is the pseudoinverse of $S$). Previously all except one algorithms for this problem ran in time $\Omega(\min\{nd^2,n^2d\})$ time with $O(k^2)$ approximation factors. One algorithm of Deshpande and Vempala did run in linear time, i.e. $O(nd)$ time but had $(k+1)!\sim k^k$ approximation factor. The contribution of this paper is to design an algorithm that achieves linear $O(nd)$ runtime and has $100(k+1)$ approximation ratio.
strength_weakness: strength_weakness: I dont see any major weaknesses. One minor weakness is that the experimental section lacks a discussion regarding the parameter settings and experimental setup for the baseline algorithms.
questions: questions: Can the authors please share the experimental details and parameter settings that are used for implementing the algorithms.
limitations: limitations: Authors have addressed limitations.

Comments: title: Response to rebuttal
comment: Thanks to the authors for the details and for addressing my question. My score remains the same.
Thank you for the rebuttal. I understand the extension to $\ell_p$ subspace approximation would not be obvious. I have increased my score.