Optimal CUR Matrix Decompositions

Christos Boutsidis; David P. Woodruff

Source-linked AI summary

Optimal CUR Matrix Decompositions

Christos Boutsidis, David P. Woodruff

arXiv:1405.7910v2May 30, 2014cs.DScs.LGmath.NA

TL;DR

The paper asks how to construct relative-error CUR decompositions with few selected rows and columns, optimal rank, input-sparsity-time performance, and deterministic guarantees. It develops randomized and deterministic algorithms using column and row subset selection, approximate SVDs, low-rank approximation within subspaces, and adaptive sampling. The algorithms achieve c = O(k/ε), r = O(k/ε), and rank(U) = k, with optimality up to constant factors and input-sparsity-time or polynomial-time implementations.

Problem
The paper addresses whether (1 + ε)-error CUR decompositions can simultaneously use optimal rows and columns, rank k, input-sparsity time, and deterministic polynomial time.
Method
The paper combines approximate SVD, leverage-score and BSS sampling, subspace low-rank approximation, and adaptive sampling to instantiate a CUR proto-algorithm.
Results
c = O(k/ε), r = O(k/ε), and rank(U) = k are achieved for relative-error CUR, with a matching lower bound up to constant factors.
Takeaways & Limitations
The resulting algorithms provide simultaneously optimal row, column, and U-rank bounds, with randomized input-sparsity-time and deterministic polynomial-time variants.
Takeaways & Limitations
The approach uses approximate rather than exact rank-k SVD factors because computing the exact choice is costly.

Abstract

from arXiv · show

The CUR decomposition of an $m \times n$ matrix $A$ finds an $m \times c$ matrix $C$ with a subset of $c < n$ columns of $A,$ together with an $r \times n$ matrix $R$ with a subset of $r < m$ rows of $A,$ as well as a $c \times r$ low-rank matrix $U$ such that the matrix $C U R$ approximates the matrix $A,$ that is, $ || A - CUR ||_F^2 \le (1+ε) || A - A_k||_F^2$, where $||.||_F$ denotes the Frobenius norm and $A_k$ is the best $m \times n$ matrix of rank $k$ constructed via the SVD. We present input-sparsity-time and deterministic algorithms for constructing such a CUR decomposition where $c=O(k/ε)$ and $r=O(k/ε)$ and rank$(U) = k$. Up to constant factors, our algorithms are simultaneously optimal in $c, r,$ and rank$(U)$.

1. Introduction.

The paper develops CUR decompositions using actual rows and columns while targeting relative-error approximation, optimal dimensions, rank, speed, and deterministic construction. It answers four open questions with randomized and deterministic algorithms supported by a matching lower bound.

1. Introduction.: CUR approximates A using selected columns C, selected rows R, and a matrix U, preserving actual data rows and columns for feature selection and interpretation.This contrasts with SVD factors, whose singular vectors are not necessarily actual rows or columns of A.
1. Introduction.: Four open questions concern optimal numbers of columns and rows, optimal rank(U), input-sparsity-time construction, and deterministic polynomial-time algorithms.These questions are posed for (1 + ε)-error CUR decompositions.
1. Introduction.: Two novel algorithms settle these questions: one randomized input-sparsity-time method and one deterministic polynomial-time method.Both target relative-error CUR with c = O(k/ε), r = O(k/ε), and rank(U) = k.
1.2. Summary of contributions.: The randomized algorithm computes the CUR factors in O(nnz(A) log n + (m + n) · poly(log n, k, 1/ε)) time.Its running time is proportional to the number of nonzero entries of A up to additional polylogarithmic and parameter-dependent terms.
1.2. Summary of contributions.: O(k/ε) columns and O(k/ε) rows are achieved with rank(U) = k, with a matching lower bound establishing optimality up to constant factors.The paper identifies this as the first relative-error CUR result with O(k/ε) rows and columns.
1.2. Summary of contributions.: The deterministic algorithm runs in O(mn^3k/ε) time and constructs c = O(k/ε) columns, r = O(k/ε) rows, and rank(U) = k.The paper presents this as addressing the open problem of deterministic polynomial-time relative-error CUR.

2. Related work.

Prior CUR methods range from weak or additive guarantees to relative-error algorithms with substantial sampling and computational costs. This paper’s algorithms target fewer rows and columns, optimal rank, and improved runtimes.

Related factorization problems: Unlike CUR, column-only factorizations A = CX + E have optimal algorithms in both spectral and Frobenius norms.The paper uses a sampling method from this line of work to select O(k) columns and rows for optimal CUR.
Prior CUR algorithms: Relative-error CUR methods improve on earlier weak or additive guarantees but require differing numbers of sampled columns and rows.The earliest methods sample by Euclidean length; later methods use subspace sampling and leverage scores.
Comparison table: Table 2.1 compares CUR algorithms by the number of columns c, rows r, and rank u of U.It summarizes the prior algorithms alongside the three algorithms developed in this work.
Prior CUR algorithms: The first relative-error CUR algorithm requires c = O(k log(k/ε^2) log δ^-1) columns and r = O(c log(c/ε^2) log δ^-1) rows, with failure probability δ.Its running time is O(mn min{m, n}) because leverage scores come from the rank-k SVD.
Comparison with prior work: Earlier relative-error algorithms use adaptive applications of the same perturbation framework, but their implementations lead to more columns, rows, and larger rank(U).The paper attributes these differences to using leverage-score sampling without the later BSS and adaptive-sampling improvements.

3. Column subset selection tools from existing literature.

This section introduces the SVD, pseudoinverse, and a perturbation lemma used to analyze column subset selection. The framework relates sampled columns to low-rank approximation quality.

Perturbation framework: Lemma 3.1 bounds a sampled-column approximation when A = AZZ^T + E and rank(Z^TS) = k.With C = AS, the lemma is used extensively in analyzing CUR algorithms.
SVD background: The SVD supplies the best rank-k approximation A_k by minimizing Frobenius error over matrices of rank at most k.The residual A_ρ-k is the difference between A and A_k.
SVD background: The top k singular vectors and singular values form the factors of A_k, while the remaining singular components form the residual.The singular values are ordered nonincreasingly.
Pseudoinverse: The Moore-Penrose pseudoinverse is the unique matrix satisfying four algebraic properties and has reciprocal nonzero singular values.These properties support later manipulations involving sampled-factor pseudoinverses.
Pseudoinverse: For products, (AB)† = B†A† holds when either factor has orthonormal columns or the factors share full column rank.The listed sufficient conditions are ATA = I_n, B^TB = I_ℓ, or rank(A) = rank(B) = n.

3.3. Fast approximate low-rank matrix approximations.

The paper replaces costly SVD computations with approximate low-rank factorization algorithms offering deterministic, randomized, or input-sparsity-time implementations.

Motivation: Approximate low-rank factorizations are used because they can be considerably faster than the SVD while providing comparable approximation quality.The paper omits algorithmic details and focuses on approximation bounds and running times.
Deterministic approximation: A deterministic method computes an orthonormal Z ∈ R^n×k for relative-error low-rank approximation using sub-cubic arithmetic operations.The cited prior work provides the deterministic procedure used in the paper’s toolkit.
Deterministic approximation: The deterministic construction first obtains Q_k and then orthonormalizes the columns of Q_k^T, with the additional step taking O(nk^2) time.This yields the desired orthonormal factorization described in the proof.
Randomized approximation: A randomized algorithm similarly computes an orthonormal Z ∈ R^n×k to accelerate SVD-based approximation.The paper identifies it as a standard randomized method for speeding up the SVD.
Input-sparsity-time approximation: An input-sparsity-time method computes an orthonormal factor with probability at least 0.99 and runtime proportional to nnz(A) plus lower-order terms.Its factorization uses orthonormal L and W with diagonal D.

3.4. Column subset selection techniques.

The section assembles deterministic sparsification, randomized row sampling, and adaptive sampling tools for selecting columns and rows while controlling approximation error.

Deterministic sparsification: BSS spectral-Frobenius sparsification selects at most r weighted vectors from a decomposition of the identity while preserving relevant matrix structure.The weights are represented through a sampling matrix S and can be computed algorithmically.
Randomized sampling: Random sampling with replacement independently selects row indices according to probabilities p_i and rescales each sampled row by 1/√(r p_i).The procedure constructs sampling and rescaling matrices Ω and D.
Adaptive sampling: Adaptive sampling chooses additional columns using residual information from already selected columns, improving constant-factor approximations to relative-error guarantees.The column-sampling implementation requires O(c1mn + c2 log c2) arithmetic operations.
Row-column symmetry: Applying the adaptive framework to A and its transpose transfers the same sampling principles between row and column selection.The paper explicitly switches between sampled rows R_1, R_2, R and columns C_1, C_2, C.

3.5. Low-rank approximations within a subspace.

This section develops methods for computing the best or nearly best rank-k approximation to A within the column span of V, including an input-sparsity-time construction.

Exact subspace approximation: O(m c^2) arithmetic operations compute the QR decomposition of V.
Exact subspace approximation: O(n c^2) arithmetic operations compute the rank-k SVD needed by the exact method.
Exact subspace approximation: The exact method constructs a rank-k approximation lying in the column span of V.It uses a QR decomposition of V followed by an SVD-based construction.
Input-sparsity-time approximation: The exact method can be costly for input-sparsity-time algorithms, motivating an approximate construction based on a sparse subspace embedding.The approximate method targets a rank-k matrix nearly as good as the best rank-k matrix in span(V).
Input-sparsity-time approximation: With probability at least 0.99, the sketched construction achieves the stated near-optimal rank-k approximation within span(V).
Input-sparsity-time approximation: O(nnz(A)+mcξ) time computes Y, Ψ, and ∆ in the approximate construction.Here ξ = O(c^2/ε^2), and ∆ contains the top k left singular vectors of the sketched matrix Ξ.

3.6. Sparse subspace embeddings.

Sparse subspace embeddings reduce dimension while preserving sparsity and geometric information, enabling fast matrix computations with high-probability norm guarantees.

Definition and construction: Sparse subspace embeddings are linear transformations that preserve both input sparsity and geometry during dimensionality reduction.
Definition and construction: W is formed by randomly hashing columns into ξ rows and applying independent random sign flips.The resulting sparse matrix has one nonzero per column, and computing WA takes O(nnz(A)) time.
Norm preservation: ξ = Ω(ρ^2 ε^-2) gives a sparse subspace embedding for a rank-ρ matrix with probability at least 0.99.
Approximation guarantees: For general matrices, ξ = Ω(ε^-2) supports the corresponding sparse embedding guarantee with probability at least 0.99.
Approximation guarantees: The embedding also supplies a high-probability guarantee for the residual of a rank-constrained approximation.

3.7. Johnson Lindestrauss transform.

The Johnson–Lindenstrauss transform reduces the dimension of matrix columns while approximately preserving their Euclidean geometry, with a fast construction for sparse inputs.

Transform and guarantee: The Johnson–Lindenstrauss transform is presented as a standard dimension-reduction method that preserves the geometry of a set of points.
Transform and guarantee: With probability at least 1 − n^-β, every column of B and its transformed counterpart have the stated approximate norm relationship.
Transform and guarantee: S is constructed with independent entries taking values ±1/√s with equal probability, and the transformed matrix is ˜B = SB.
Implementation: O(nnz(B) log n) arithmetic operations construct the transformed matrix ˜B.

3.8. Generalized rank-constrained matrix approximations.

This section characterizes the minimum-Frobenius-norm rank-constrained solution for a generalized matrix approximation problem involving C and R.

Optimal solution: The solution Uopt has rank at most k and minimizes the Frobenius norm among all feasible solutions.
Prior result: The result is attributed to prior work, with a spectral-norm version cited separately.

4. New column subset selection tools.

This section develops new subset-selection and derandomization tools for CUR algorithms, including deterministic adaptive sampling and input-sparsity-time sparsification. These tools support selecting columns and rows while preserving approximation guarantees.

New tools: The paper combines previous subset-selection methods with novel tools for constructing optimal CUR algorithms.The results are also presented as potentially independently interesting.
Deterministic adaptive sampling: 4n-discretization produces a distribution q satisfying qi ≥ pi/4 for every item while enabling pairwise-independent hash-based sampling.The discrete probabilities are integer multiples of 1/(4n), allowing samples to be generated from a pairwise-independent hash family.
Deterministic adaptive sampling: Pairwise independence suffices to control the variance of adaptive sampling averages during derandomization.The variance decreases by a factor of 1/c2 when xj averages c2 pairwise-independent samples.
Deterministic adaptive sampling: The resulting deterministic adaptive-column procedure augments existing columns with c2 columns and supports a residual-based approximation guarantee.The procedure is denoted AdaptiveColsD(A, V, c2).
Input-sparsity-time sparsification: The input-sparsity-time tools include fast versions of BSS sampling and adaptive sampling for constructing CUR factors.The paper combines BSS ideas with sparse subspace embeddings and develops fast adaptive row and column procedures.
Input-sparsity-time sparsification: Input-sparsity-time dual-set spectral-Frobenius sparsification applies a sparse subspace embedding before selecting at most r weighted vectors.With probability at least 0.98, the resulting sampling matrix satisfies the stated sparsification guarantee.

5. Linear-time randomized CUR.

This section presents a simple randomized CUR algorithm that constructs sampled columns, rows, and a rank-k middle factor. Its approximation analysis establishes a relative-error guarantee with constant success probability.

Algorithm: The algorithm is designed for simplicity and may be faster on dense matrices than the alternative randomized algorithm.The paper states that the alternative may be faster especially on sparse matrices, while this algorithm can be faster for dense inputs depending on parameters and dimensions.
Algorithm: Algorithm 2 returns c = O(k/ε) columns, r = O(k/ε) rows, and a matrix U with rank at most k.The inputs are A, a rank parameter k, and an accuracy parameter ε.
Running time: The algorithm’s total arithmetic running time is given by O(n^2k/ε + m^2k/ε + mk^2/ε^2 + k^3/ε^3 + k^4 ln k + k…).The displayed bound is the stated total asymptotic running time, though the supplied passage truncates its final terms.
Guarantee: Theorem 5.1 states that the resulting C, U, and R satisfy the analyzed approximation bound with probability at least 0.2.The proof combines intermediate column, row, and rank-k approximation results.
Column selection: The construction first obtains O(k) columns that give a constant-factor low-rank approximation, then adaptively samples O(k/ε) additional columns for relative error.The adaptive step produces c = O(k) + O(k/ε) columns and uses the best rank-k approximation in span(C).
Row selection: The same adaptive strategy is applied to the transpose to sample rows after constructing the column matrix.The sampled rows depend on the preceding column-sampling process.

6. Input-Sparsity-Time CUR.

This section replaces non-sparse components with sparse embeddings, approximate SVD, and fast sampling procedures to obtain an input-sparsity-time CUR algorithm. The resulting algorithm retains O(k/ε) rows and columns and rank at most k.

Algorithm: Algorithm 3 returns c = O(k/ε) columns, r = O(k/ε) rows, and U with rank at most k.It takes A, k, and ε as input and follows the paper’s CUR proto-algorithm.
Construction: The algorithm first forms O(k) columns, adaptively adds O(k/ε) columns, then applies analogous row sampling and constructs U of rank k.The column and row counts are explicitly assembled from the initial and adaptive sample sizes.
Input-sparsity implementation: The construction uses a sparse subspace embedding W with ξ = Ω(k^2ε^-2) and computes WA in O(nnz(A)) time.The embedding supports sparse implementations of the sampling stages.
Algorithm: Algorithm 3 uses SparseSVD, sparse BSS sampling, and sparse adaptive row and column sampling in place of standard procedures.These substitutions are the central changes from Algorithm 2.
Running time: The total running time is O(nnz(A) log n + (m + n) · poly(log n, k, 1/ε)).This is the stated overall asymptotic running time.

7. Deterministic CUR.

The deterministic algorithm constructs an optimal relative-error rank-k CUR decomposition in polynomial time, selecting O(k/ε) columns and rows and constructing U with rank at most k.

Algorithm and guarantees: The deterministic polynomial-time algorithm returns C with O(k/ε) columns, R with O(k/ε) rows, and U with rank at most k.It selects columns, selects rows, and then constructs an intersection matrix with optimal rank.
Algorithm and guarantees: The algorithm follows a three-step pipeline: deterministic column selection, deterministic row selection, and construction of a rank-k intersection matrix U.Its component procedures, including DeterministicSVD, BssSampling, and adaptive column and row routines, are all deterministic.
Running time: The total asymptotic running time is O(mn^3k/ε).The dominant costs arise in constructing the adaptive column and row sets.
Intermediate approximation results: The first sampled column set preserves the rank of the relevant singular-vector subspace, with rank(Z1^T S1) = k.This rank preservation supports the subsequent constant-factor column-based approximation.
Intermediate approximation results: The constructed C provides a relative-error column-based low-rank approximation, and a rank-k matrix in span(C) achieves a similar bound.The proof combines the column approximation lemmas with adaptive sampling and best-subspace computation.
Intermediate approximation results: The row-selection analysis mirrors the column-selection analysis, using analogous rank-preservation and approximation results for R1.The proof applies the corresponding column arguments to the transpose of A.

8. Lower bound.

The lower-bound construction shows that relative-error CUR decompositions require, up to constants, Ω(k/ε) columns and rows and rank Ω(k) for U.

Main lower bound: Relative-error CUR is impossible unless C has Ω(k/ε) columns, R has Ω(k/ε) rows, and U has rank Ω(k).The proof uses a symmetric hard instance and separately establishes column, row, and intersection-rank requirements.
Contradiction argument: If c = o(k/ε), the columns of CUR lie in the span of C, contradicting the column lower bound.The factorization therefore cannot achieve relative error with too few selected columns.
Contradiction argument: If r = o(k/ε), the rows of CUR lie in the row space of R, contradicting the row lower bound.This yields the analogous necessity for the number of selected rows.
Intersection-rank lower bound: The hard instance also rules out sufficiently low-rank U, since rank(U) < k/2 forces rank(CUR) < k/2 and incurs excessive error.The construction contains repeated diagonal blocks that make this rank deficiency incompatible with ε < 1/3.
Column lower bound: A symmetric hard matrix admits no (1+ε)-approximation from o(k/ε) selected columns.This is formalized by Lemma 8.2 for column-based approximations.
Row lower bound: By symmetry, no subset of o(k/ε) rows of the hard matrix spans a (1+ε)-approximation.Corollary 8.3 transfers the column lower bound to rows.