The Road to Low-Rank Approximation (5): CUR
Back once more to the road of low-rank approximation. In Low-Rank Approximation, Part 4: ID, we introduced "Interpolative Decomposition (ID)", which is the process of finding an approximation of the form $\boldsymbol{C}\boldsymbol{Z}$ for a matrix $\boldsymbol{M}\in\mathbb{R}^{n\times m}$, where $\boldsymbol{C}\in\mathbb{R}^{n\times r}$ consists of several columns of the matrix $\boldsymbol{M}$, and $\boldsymbol{Z}\in\mathbb{R}^{r\times m}$ is an arbitrary matrix.
In this article we will introduce CUR decomposition, which shares the same underlying idea as interpolative decomposition: both build an approximation of the original matrix using its own rows and/or columns as a "skeleton". Unlike ID, which uses only rows or only columns, CUR decomposition uses both rows and columns simultaneously.
Basic Definition
This is actually not the first time CUR decomposition has appeared on this site. Back in Nyströmformer: A Linear-Complexity Attention Scheme Based on Matrix Decomposition, we already introduced the Nyström approximation of a matrix, which is in fact CUR decomposition. Later, in Using CUR Decomposition to Accelerate Retrieval for Interactive Similarity Models, we also discussed the application of CUR decomposition to reducing the retrieval complexity of interactive similarity models. more
The reason CUR decomposition supports these applications lies precisely in the "C" and "R" of its name. Specifically, CUR decomposition seeks an approximation of matrix $\boldsymbol{M}\in\mathbb{R}^{n\times m}$ of the following form:
\begin{equation}\mathop{\text{argmin}}_{S_1,S_2,\boldsymbol{\mathcal{U}}}\Vert \underbrace{\boldsymbol{M}_{[:,S_1]}}_{\boldsymbol{\mathcal{C}}}\boldsymbol{\mathcal{U}}\underbrace{\boldsymbol{M}_{[S_2,:]}}_{\boldsymbol{\mathcal{R}}} - \boldsymbol{M}\Vert_F^2\quad\text{s.t.}\quad \left\{\begin{aligned}&S_1\subset\{0,1,\cdots,m-1\},|S_1|=r\\ &S_2\subset\{0,1,\cdots,n-1\},|S_2|=r \\ &\boldsymbol{\mathcal{U}}\in\mathbb{R}^{r\times r} \end{aligned}\right.\end{equation}
To distinguish it from the $\boldsymbol{U}$ of SVD, we use the script $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{R}}$ here. For comparison, the ID introduced in the previous article was
\begin{equation}\mathop{\text{argmin}}_{S,\boldsymbol{Z}}\Vert \underbrace{\boldsymbol{M}_{[:,S]}}_{\boldsymbol{C}}\boldsymbol{Z} - \boldsymbol{M}\Vert_F^2\quad\text{s.t.}\quad \left\{\begin{aligned} &S\subset\{0,1,\cdots,m-1\},|S|=r\\ &\boldsymbol{Z}\in\mathbb{R}^{r\times m} \end{aligned}\right.\end{equation}
while the low-rank approximation found by SVD is
\begin{equation}\mathop{\text{argmin}}_{\boldsymbol{U},\boldsymbol{\Sigma},\boldsymbol{V}}\Vert \boldsymbol{U}_{[:,:r]}\boldsymbol{\Sigma}_{[:r,:r]}\boldsymbol{V}_{[:,:r]}^{\top} - \boldsymbol{M}\Vert_F^2\quad\text{s.t.}\quad \left\{\begin{aligned} &\boldsymbol{U}\in\mathbb{R}^{n\times n}, \boldsymbol{U}^{\top}\boldsymbol{U} = \boldsymbol{I}_n \\ &\boldsymbol{V}\in\mathbb{R}^{m\times n}, \boldsymbol{V}^{\top}\boldsymbol{V} = \boldsymbol{I}_m \\ &\boldsymbol{\Sigma}=\text{diag}(\sigma_1,\cdots,\sigma_{\min(n,m)})\in\mathbb{R}_{\geq 0}^{n\times m} \end{aligned}\right.\end{equation}
In the SVD article we proved that SVD finds the optimal solution for a $r$-rank approximation, but it is itself expensive to compute, and $\boldsymbol{U},\boldsymbol{V}$ doesn't have a very intuitive physical meaning. By contrast, CUR decomposition replaces $\boldsymbol{\mathcal{U}},\boldsymbol{V}$ with the original matrix's own columns $\boldsymbol{\mathcal{C}}$ and rows $\boldsymbol{\mathcal{R}}$; although it is not as accurate as SVD, it is superior in terms of interpretability, storage cost, and computational cost.
On the surface, the left and right matrices $\boldsymbol{U},\boldsymbol{V}$ of the SVD approximation are more complex while the middle matrix $\boldsymbol{\Sigma}$ is simpler, whereas CUR decomposition is the opposite: the left and right matrices $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$ are simpler while the middle matrix $\boldsymbol{\mathcal{U}}$ is more complex.
Choosing U
Clearly, the difficulty of CUR decomposition lies in choosing the rows and columns, because once $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$ is given, the optimal solution for $\boldsymbol{\mathcal{U}}$ can be expressed analytically using the pseudo-inverse:
\begin{equation}\boldsymbol{\mathcal{U}}^* = \boldsymbol{\mathcal{C}}^{\dagger}\boldsymbol{M}\boldsymbol{\mathcal{R}}^{\dagger}\end{equation}
The derivation can be found in the article on pseudo-inverses. In fact this solution is quite intuitive: if $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$ were both invertible matrices, then the solution to the equation $\boldsymbol{\mathcal{C}}\boldsymbol{\mathcal{U}}\boldsymbol{\mathcal{R}}=\boldsymbol{M}$ would naturally be $\boldsymbol{\mathcal{U}}=\boldsymbol{\mathcal{C}}^{-1}\boldsymbol{M}\boldsymbol{\mathcal{R}}^{-1}$, and when they are not invertible we simply replace the inverse ${}^{-1}$ with the pseudo-inverse ${}^{\dagger}$.
Besides this theoretically optimal solution, there is another commonly used, and in some sense more intuitive, choice for CUR decomposition:
\begin{equation}\boldsymbol{\mathcal{U}} = \boldsymbol{M}_{[S_2,S_1]}^{\dagger}\end{equation}
Note that slicing takes precedence over transpose and pseudo-inverse, so this $\boldsymbol{\mathcal{U}}$ is actually the pseudo-inverse of the submatrix formed by the common part of $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$.
How should we understand this choice? By permuting rows and columns, we can move the selected rows and columns to the front, and assuming $\boldsymbol{M}_{[S_2,S_1]}$ is invertible, the result can be written in block-matrix form as
\begin{equation}\underbrace{\begin{pmatrix}\boldsymbol{A} & \boldsymbol{B} \\ \boldsymbol{C} & \boldsymbol{D}\end{pmatrix}}_{\boldsymbol{M}} \approx \underbrace{\begin{pmatrix}\boldsymbol{A} \\ \boldsymbol{C}\end{pmatrix}}_{\boldsymbol{\mathcal{C}}}\,\,\underbrace{\boldsymbol{A}^{-1}}_{\boldsymbol{\mathcal{U}}}\,\,\underbrace{\begin{pmatrix}\boldsymbol{A} & \boldsymbol{B}\end{pmatrix}}_{\boldsymbol{\mathcal{R}}} = \begin{pmatrix}\boldsymbol{A} & \boldsymbol{B} \\ \boldsymbol{C} & \boldsymbol{C}\boldsymbol{A}^{-1}\boldsymbol{B}\end{pmatrix}\label{eq:id-abcd}\end{equation}
We can see that in this case the CUR decomposition exactly reconstructs the selected $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C}$ (or $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$), and uses $\boldsymbol{C}\boldsymbol{A}^{-1}\boldsymbol{B}$ to approximate $\boldsymbol{D}$; in this sense CUR decomposition becomes a kind of "matrix completion" method.
It's worth pointing out that since both $\boldsymbol{\mathcal{U}}$ involve pseudo-inverses, and the definition of the pseudo-inverse does not require a square matrix, the most general form of CUR decomposition does not actually require $\boldsymbol{\mathcal{C}}$/$\boldsymbol{\mathcal{R}}$ to have the same number of columns/rows. If necessary, we can choose different numbers of columns/rows for $\boldsymbol{\mathcal{C}}$/$\boldsymbol{\mathcal{R}}$.
Row and Column Selection
Having resolved $\boldsymbol{\mathcal{U}}$, what remains is mainly the choice of $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$. Since choosing rows and choosing columns is essentially equivalent, we'll use column selection as our example below.
That is, our task below is to select $r$ key columns from the matrix $\boldsymbol{M}$, to serve as its "skeleton" — which could also be called an "outline" or "sketch". This is actually a problem we already explored in the previous two articles (the CR article and the ID article), and the schemes discussed there can also be used to construct the $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{R}}$ of CUR decomposition, including
1. Choosing the $r$ columns with the largest norm;
2. Randomly sampling $r$ columns with probability weighted by norm;
3. Uniformly randomly sampling $r$ columns;
4. Selecting the first $r$ columns via column-pivoted QR decomposition.
Each of these schemes has its own strengths and weaknesses, along with its own applicable scenarios and implicit assumptions. Beyond these, we can also consider some more intuitive approaches. For instance, given that the notion of "key column" seems similar to that of a "cluster center", we could cluster the $n$ column vectors into $k$ clusters, and then select the $k$ vectors closest to the cluster centers. When $n$ is too large, we can first randomly subsample a subset and then run the above selection algorithm on that subset.
In general, column selection is a classic problem in matrix approximation; the relevant English search terms are Randomized Linear Algebra, Column Subset Selection, and so on — a quick search will turn up plenty of material.
Leverage Scores
Naturally, as befits a new article, we should introduce some new methods too, so next we'll present two more approaches to column selection. The first we call "Leverage Scores", which performs column selection via the idea of linear regression.
First, we regard the matrix $\boldsymbol{M}$ as $m$ samples of dimension $n$, and correspondingly we have $m$ vectors of dimension $d$ making up a target matrix $\boldsymbol{Y}$. Our task is to predict $\boldsymbol{Y}$ using $\boldsymbol{M}$, with the simplest linear model, and the optimization objective is least squares:
\begin{equation}\boldsymbol{W}^* = \mathop{\text{argmin}}_{\boldsymbol{W}} \Vert\boldsymbol{Y} - \boldsymbol{W}\boldsymbol{M}\Vert_F^2\label{eq:linear-loss}\end{equation}
We've already solved this objective in the pseudo-inverse article; the answer is $\boldsymbol{W}^* = \boldsymbol{Y}\boldsymbol{M}^{\dagger}$. Assuming $n < m$ and that $\boldsymbol{M}$ has rank $n$, we can further write $\boldsymbol{W}^* = \boldsymbol{Y}\boldsymbol{M}^{\top}(\boldsymbol{M}\boldsymbol{M}^{\top})^{-1}$, giving us
\begin{equation}\hat{\boldsymbol{Y}} = \boldsymbol{W}^*\boldsymbol{M} = \boldsymbol{Y}\boldsymbol{M}^{\top}(\boldsymbol{M}\boldsymbol{M}^{\top})^{-1}\boldsymbol{M} = \boldsymbol{Y}\boldsymbol{H}\end{equation}
Here $\boldsymbol{H}=\boldsymbol{M}^{\top}(\boldsymbol{M}\boldsymbol{M}^{\top})^{-1}\boldsymbol{M}$ is called the "Hat Matrix", supposedly because it turns $\boldsymbol{Y}$ into $\hat{\boldsymbol{Y}}$ — as if putting a hat on $\boldsymbol{Y}$ (i.e., $\hat{}$). Let $\boldsymbol{m}_i$ be the $i$-th column vector of $\boldsymbol{M}$ — here, the $i$-th sample — then we take
\begin{equation}\boldsymbol{H}_{i,i} = \boldsymbol{m}_i^{\top}(\boldsymbol{M}\boldsymbol{M}^{\top})^{-1}\boldsymbol{m}_i\end{equation}
to measure that sample's contribution when predicting $\hat{\boldsymbol{Y}}$; this is the "Leverage Score". We consider that choosing $r$ key columns is equivalent to choosing the $r$ most important samples, so we can select the $r$ columns with the largest $\boldsymbol{H}_{i,i}$.
When $\boldsymbol{M}\boldsymbol{M}^{\top}$ is not invertible, the paper Input Sparsity Time Low-Rank Approximation via Ridge Leverage Score Sampling generalizes Leverage Scores to "Ridge Leverage Scores", essentially by adding a regularization term to the objective $\eqref{eq:linear-loss}$ to make it invertible. But in fact, as we know, the concept of the pseudo-inverse does not require full rank, so we can compute the pseudo-inverse directly via SVD, without needing to introduce an extra regularization term.
Let the SVD of $\boldsymbol{M}$ be $\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}$, then
\begin{equation}\boldsymbol{H} = \boldsymbol{M}^{\dagger}\boldsymbol{M} = (\boldsymbol{V}\boldsymbol{\Sigma}^{\dagger}\boldsymbol{U}^{\top})(\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}) = \boldsymbol{V}(\boldsymbol{\Sigma}^{\dagger}\boldsymbol{\Sigma})\boldsymbol{V}^{\top}\end{equation}
Assuming $\boldsymbol{M}$ has rank $\gamma$ ($\gamma$ need not equal $r$), then according to the rules for computing the pseudo-inverse, $\boldsymbol{\Sigma}^{\dagger}\boldsymbol{\Sigma}$ is a $m\times m$ diagonal matrix whose first $\gamma$ diagonal entries are all 1 and the rest are 0, so
\begin{equation}\boldsymbol{H} = \boldsymbol{V}_{[:,:\gamma]}\boldsymbol{V}_{[:,:\gamma]}^{\top}\quad\Rightarrow\quad \boldsymbol{H}_{i,i} = \Vert\boldsymbol{V}_{[i-1,:\gamma]}\Vert^2\end{equation}
Note that $\boldsymbol{H}_{i,i}$ denotes the entry of $\boldsymbol{H}$ in row $i$, column $i$, with counting starting from 1, but slicing follows Python convention with counting starting from 0, so the final slice is ${}_{[i-1,:\gamma]}$. We now see that to compute the leverage scores of the columns of $\boldsymbol{M}$, we only need to take its SVD $\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{\top}$ and then compute the squared norm of each row of $\boldsymbol{V}_{[:,:\gamma]}$; likewise, to compute the leverage scores of the rows, we only need to compute the squared norm of each row of $\boldsymbol{U}_{[:,:\gamma]}$.
Since $\boldsymbol{V}$ itself is an orthogonal matrix, the following always holds:
\begin{equation}\sum_{i=1}^m \boldsymbol{H}_{i,i} = \sum_{i=1}^m\Vert\boldsymbol{V}_{[i-1,:\gamma]}\Vert^2 = \gamma\end{equation}
So besides selecting the $r$ columns with the largest leverage score, we can also construct a distribution $p_i = \boldsymbol{H}_{i,i} / \gamma$ for random sampling.
Leverage scores are related to the rank $\gamma$ of $\boldsymbol{M}$, and the rank of $\boldsymbol{M}$ equals the number of nonzero singular values of $\boldsymbol{M}$, so it is affected by singular values close to zero. But these smaller singular values contribute little in practice, so in practice $\gamma$ is usually taken to be the number of significant (dominant) singular values $\boldsymbol{M}$. Another issue with leverage scores is that they require SVD upfront; in practice, one often resorts to approximate SVD algorithms. As for the details of approximate SVD algorithms, we'll leave that for another occasion.
The DEIM Method
The other column-selection method we want to introduce is called DEIM, short for Discrete Empirical Interpolation Method. We won't dig into the origin of the name here, but broadly speaking, both leverage scores and DEIM are commonly used column-selection methods for CUR decomposition, and DEIM is more closely tied to CUR decomposition, which is why it has become increasingly popular in recent years.
DEIM's starting point is the identity $\eqref{eq:id-abcd}$: under the $\boldsymbol{\mathcal{C}},\boldsymbol{\mathcal{U}},\boldsymbol{\mathcal{R}}$ of this identity, the error of the CUR decomposition depends on $\Vert \boldsymbol{D} - \boldsymbol{C}\boldsymbol{A}^{-1}\boldsymbol{B}\Vert_F$. When would this expression be small? Intuitively, if $\boldsymbol{A}$ is fairly large while $\boldsymbol{B},\boldsymbol{C},\boldsymbol{D}$ are both fairly small, the whole expression should certainly be small too. But $\boldsymbol{A}$ is a matrix — how do we measure its "size"? The absolute value of the determinant can serve as one reference metric. So one feasible scheme is to choose the rows and columns that maximize the absolute value of the determinant of $\boldsymbol{M}_{[S_2,S_1]}$.
Of course, this scheme is of theoretical value only, since exactly finding the submatrix with the largest determinant in absolute value is also NP-hard. But it does give us a target we can try to approximate with a greedy solution: when $r=1$, finding the largest determinant in absolute value amounts to finding the largest element in absolute value, which is quite tractable, and this can then be applied recursively. DEIM follows this line of thinking, but instead of starting from $\boldsymbol{M}$, it borrows the approach used for leverage scores and starts from the post-SVD $\boldsymbol{V}$.
Leverage scores turn the problem of finding key columns in $\boldsymbol{M}$ into finding key rows in $\boldsymbol{V}_{[:,:\gamma]}$, ranked by squared row norm; DEIM instead tries to find key rows by seeking a CUR approximation for $\boldsymbol{V}$. But wait — we haven't even finished solving the CUR problem for $\boldsymbol{M}$, and now there's another CUR problem for $\boldsymbol{V}$? Doesn't that make things even more complicated? Not really — this one is simpler, because $\boldsymbol{V}$ is the result of SVD on $\boldsymbol{M}$, and it has already been sorted by singular-value magnitude. So we can regard the columns of $\boldsymbol{V}$ as already ranked by importance, meaning the most important $r$ columns must be $\boldsymbol{V}_{[:,:r]}$ itself — we only need to select the rows.
As mentioned above, the solving strategy is a greedy algorithm. The most important column is naturally the first one, $\boldsymbol{V}_{[:,0]}$ — so what's the most important row? We should select the row where, among the entries at the intersection of that row and the first column, the absolute value is largest — in other words, the row containing the largest-magnitude element of the first column. This gives us our initial row and column. Suppose we have already selected $k$ key rows, with index set $S_k$; how do we select the $k+1$-th key row? First, we know that the CUR approximation built from the already-selected $k$ rows and the first $k$ columns is $\boldsymbol{V}_{[:,:k]}\boldsymbol{V}_{[S_k,:k]}^{-1}\boldsymbol{V}_{[S_k,:]}$, and the error in the $k+1$-th column is
\begin{equation}\boldsymbol{V}_{[:,k]} - \boldsymbol{V}_{[:,:k]}\boldsymbol{V}_{[S_k,:k]}^{-1}\boldsymbol{V}_{[S_k,k]}\end{equation}
From equation $\eqref{eq:id-abcd}$ we know that this kind of CUR approximation exactly recovers the selected rows and columns, so in the expression above, the entries corresponding to the already-selected $k$ rows must be zero. Hence the remaining largest-magnitude nonzero element cannot lie in one of the already-selected $k$ rows, and we select the row it belongs to as the $k+1$-th key row.
In short, DEIM exploits the fact that the SVD has already sorted the column vectors of $\boldsymbol{V}$ by importance, turning CUR decomposition into a pure row-search problem, thereby reducing the search space, and then solves it via a greedy algorithm, where at each step the choice is based on which row contains the largest-error element. For a more detailed treatment and proofs, see A DEIM Induced CUR Factorization and CUR Matrix Factorizations: Algorithms, Analysis, Applications.
Summary
This article introduced CUR decomposition, which can be seen as a further extension of the interpolative decomposition (ID) discussed in the previous article. Its distinguishing feature is that it builds a low-rank approximation using both several rows and several columns of the original matrix as its skeleton simultaneously.
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.