The Road to Low-Rank Approximation (IV): ID

低秩近似之路（四）：ID 苏剑林 · 2024-10-30 · Mathematics approximation optimal matrix low-rank

The protagonist of this article is ID (Interpolative Decomposition), which can also be understood as a low-rank decomposition with a specific structure, where one side of the decomposition consists of several columns of the original matrix (of course, if you have a preference for rows, choosing rows works just as well). In other words, ID tries to find a number of key columns from a matrix to serve as a "skeleton" (often also called a "sketch") to approximate the original matrix.

Many readers may never have heard of ID — even Wikipedia only offers a few vague sentences about it (link). But in fact, ID has long been built into SciPy just like SVD (see scipy.linalg.interpolative), which is itself indirect evidence of ID's practical value.

Basic Definition

In the first three articles of this series, we introduced the pseudo-inverse, SVD, and CR approximation, all of which can be viewed as ways of finding low-rank approximations with a particular structure:

\begin{equation}\mathop{\text{argmin}}_{\text{rank}(\tilde{\boldsymbol{M}})\leq r}\Vert \tilde{\boldsymbol{M}} - \boldsymbol{M}\Vert_F^2\end{equation}more

where $\boldsymbol{M}\in\mathbb{R}^{n\times m}$. When no further constraints are added, the optimal solution is given by SVD; when we require $\tilde{\boldsymbol{M}}=\boldsymbol{A}\boldsymbol{B}$ and $\boldsymbol{A},\boldsymbol{B}$, with one of them given and we solve for the optimal other half, the optimal solution can be given via the pseudo-inverse; if we require $\boldsymbol{M}=\boldsymbol{X}\boldsymbol{Y}$ and $\tilde{\boldsymbol{M}}=\boldsymbol{X}_{[:, S]}\boldsymbol{Y}_{[S,:]}$, that's the problem CR approximation is concerned with.

CR approximation constructs a low-rank approximation by selecting rows/columns of the original matrix, which makes the approximation result more interpretable and also applicable to some nonlinear scenarios. But the premise of CR approximation is that the matrix $\boldsymbol{M}$ itself arises as the product of two matrices — its original motivation was to reduce the computational cost of matrix multiplication. For scenarios where the matrix $\boldsymbol{M}$ is given directly, an analogous low-rank approximation is instead provided by ID.

Specifically, in ID we have $\tilde{\boldsymbol{M}}=\boldsymbol{C}\boldsymbol{Z}$, where $\boldsymbol{C}=\boldsymbol{M}_{[:,S]}$ consists of several columns of $\boldsymbol{M}$, and $\boldsymbol{Z}$ is arbitrary — that is, we approximate $\boldsymbol{M}$ using some of its own columns as the skeleton:

\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 S\subset\{0,1,\cdots,m-1\},|S|=r,\boldsymbol{Z}\in\mathbb{R}^{r\times m}\end{equation}

According to the results in The Road to Low-Rank Approximation (I): The Pseudo-Inverse, if $\boldsymbol{C}$ has already been determined, then the optimal solution for $\boldsymbol{Z}$ is $\boldsymbol{C}^{\dagger} \boldsymbol{M}$. So the real difficulty of ID lies entirely in the optimization of $S$, i.e., the selection of columns. This is a combinatorial optimization problem, and solving it exactly is NP-hard, so current efforts are mainly focused on finding approximation algorithms that strike a good balance between efficiency and accuracy.

Geometric Meaning

Before we try to solve it, let's first get a better understanding of the geometric meaning of ID, which will help us better understand its use cases and solution strategies. Let's write $\boldsymbol{C}$ in column-vector form as $\boldsymbol{C}=(\boldsymbol{c}_1,\boldsymbol{c}_2,\cdots,\boldsymbol{c}_r)$. Then for any column vector $\boldsymbol{z}=(z_1,z_2,\cdots,z_r)^{\top}$, we have

\begin{equation}\boldsymbol{C}\boldsymbol{z} = \begin{pmatrix}\boldsymbol{c}_1 & \boldsymbol{c}_2 & \cdots & \boldsymbol{c}_r\end{pmatrix}\begin{pmatrix}z_1 \\ z_2 \\ \vdots \\ z_r\end{pmatrix} = \sum_{i=1}^r z_i \boldsymbol{c}_i\end{equation}

So the geometric meaning of $\boldsymbol{C}\boldsymbol{z}$ is a linear combination of the column vectors of $\boldsymbol{C}$. Note that $\boldsymbol{c}_1,\boldsymbol{c}_2,\cdots,\boldsymbol{c}_r$ is selected from $\boldsymbol{M}=(\boldsymbol{m}_1,\boldsymbol{m}_2,\cdots,\boldsymbol{m}_m)$, so ID amounts to choosing a number of columns as (approximate) basis vectors, and expressing all the remaining columns as linear combinations of these basis vectors — this is the meaning of the "I" (Interpolative) in ID.

As we know, "Interpolative" more precisely means "interpolating," and to better highlight this "interpolating" property, some references add the condition $|z_{i,j}| \leq 1$ to the definition of ID (where $z_{i,j}$ is an arbitrary element of the matrix $\boldsymbol{Z}$). Of course this condition is actually rather demanding — guaranteeing it holds strictly is likely also NP-hard, so many references relax it to $|z_{i,j}| \leq 2$, a bound that most approximation algorithms can satisfy in practice. If there are no other requirements and we only care about optimizing the approximation error, this restriction can also simply be dropped.

QR Decomposition

Algorithms for solving ID fall into two broad categories: deterministic algorithms and randomized algorithms. Deterministic algorithms are more computationally expensive but tend to give better approximations, whereas randomized algorithms are more computationally efficient but slightly less accurate. Note that these are all merely approximation algorithms that perform reasonably well in practice, and none of them rule out the possibility of extreme cases where they fail entirely.

The first algorithm regarded as a standard approach is based on QR decomposition — more precisely, column-pivoting QR decomposition (常译作"列主元QR分解" in Chinese, though the author feels it might as well be more freely translated as "column-driven QR decomposition"). It is a deterministic algorithm. Why is ID connected to QR decomposition? We can understand this by starting from how $\boldsymbol{Z}$ is computed.

As mentioned earlier, if $\boldsymbol{C}$ is given, the optimal solution for $\boldsymbol{Z}$ is $\boldsymbol{C}^{\dagger}\boldsymbol{M}$. This answer is of course correct, but not very intuitive. Without loss of generality, assume $\boldsymbol{c}_1,\boldsymbol{c}_2,\cdots,\boldsymbol{c}_r$ are linearly independent. Then, from a geometric point of view, finding the optimal approximation of the form $\boldsymbol{C}\boldsymbol{Z}$ is exactly the same as projecting each column vector of $\boldsymbol{M}$ onto the $r$-dimensional subspace spanned by $\boldsymbol{c}_1,\boldsymbol{c}_2,\cdots,\boldsymbol{c}_r$. To compute this projection, we can first apply Gram-Schmidt orthogonalization to $\boldsymbol{c}_1,\boldsymbol{c}_2,\cdots,\boldsymbol{c}_r$ to turn it into an orthonormal basis, and then projecting onto an orthonormal basis becomes much simpler. This orthogonalization process is, naturally, exactly QR decomposition.

Gram-Schmidt orthogonalization recursively performs the following steps:

\begin{equation}\boldsymbol{q}_1 = \frac{\boldsymbol{c}_1}{\Vert\boldsymbol{c}_1\Vert},\quad \boldsymbol{q}_k = \frac{\hat{\boldsymbol{q}}_k}{\Vert\hat{\boldsymbol{q}}_k\Vert},\quad\hat{\boldsymbol{q}}_k = \boldsymbol{c}_k - \sum_{i=1}^{k-1} (\boldsymbol{c}_k^{\top} \boldsymbol{q}_i)\boldsymbol{q}_i,\quad k = 2,3,\cdots,r\end{equation}

The result expresses $\boldsymbol{C}$ as:

\begin{equation}\boldsymbol{C} = \underbrace{\begin{pmatrix}\boldsymbol{q}_1 & \boldsymbol{q}_2 & \cdots & \boldsymbol{q}_r\end{pmatrix}}_{\boldsymbol{Q}}\underbrace{\begin{pmatrix}R_{1,1} & R_{1,2} & \cdots & R_{1,r} \\ 0 & R_{2,2} & \cdots & R_{2,r} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_{r,r} \\ \end{pmatrix}}_{\boldsymbol{R}}\end{equation}

Given $\boldsymbol{q}_1,\boldsymbol{q}_2,\cdots,\boldsymbol{q}_r$, the optimal approximation and error for the $k$-th column $\boldsymbol{m}_k$ of matrix $\boldsymbol{M}$, projected onto $\boldsymbol{C}$, are respectively

\begin{equation}\sum_{i=1}^r (\boldsymbol{m}_k^{\top} \boldsymbol{q}_i)\boldsymbol{q}_i\qquad\text{and}\qquad \left\Vert\boldsymbol{m}_k - \sum_{i=1}^r (\boldsymbol{m}_k^{\top} \boldsymbol{q}_i)\boldsymbol{q}_i\right\Vert^2\end{equation}

Column-Pivoting QR

Of course, the above results were obtained under the assumption that $\boldsymbol{C}$ is already known. So how do we pick a good set of $r$ columns from $\boldsymbol{M}$ to form $\boldsymbol{C}$? Column-pivoting QR decomposition offers one answer.

Normally, if we perform Gram-Schmidt orthogonalization on $\boldsymbol{m}_1,\boldsymbol{m}_2,\cdots,\boldsymbol{m}_m$, we proceed in order, starting from $\boldsymbol{m}_1$, then $\boldsymbol{m}_2,\boldsymbol{m}_3,\cdots$, and so on. Column-pivoting QR decomposition instead modifies the orthogonalization order based on vector norms, which can be written as

\begin{equation}\begin{gathered} \boldsymbol{q}_1 = \frac{\boldsymbol{m}_{\rho_1}}{\Vert\boldsymbol{m}_{\rho_1}\Vert},\quad \boldsymbol{q}_k = \frac{\hat{\boldsymbol{q}}_k}{\Vert\hat{\boldsymbol{q}}_k\Vert},\quad\hat{\boldsymbol{q}}_k = \boldsymbol{m}_{\rho_k} - \sum_{i=1}^{k-1} (\boldsymbol{m}_{\rho_k}^{\top} \boldsymbol{q}_i)\boldsymbol{q}_i \\ \rho_1 = \mathop{\text{argmax}}_{i\in\{1,2,\cdots,m\}} \Vert \boldsymbol{m}_i\Vert,\quad \rho_k = \mathop{\text{argmax}}_{i\in\{1,2,\cdots,m\}\backslash\{\rho_1,\rho_2,\cdots,\rho_{k-1}\}} \left\Vert \boldsymbol{m}_i - \sum_{j=1}^{k-1} (\boldsymbol{m}_i^{\top} \boldsymbol{q}_j)\boldsymbol{q}_j\right\Vert \end{gathered}\end{equation}

In plain terms, column-pivoting QR decomposition simply chooses, at each step, the remaining column with the largest error to orthonormalize next. Aside from this change in execution order, column-pivoting QR decomposition involves no other differences from ordinary QR decomposition in terms of the actual computation, so its final form can be written as

\begin{equation}\boldsymbol{M}\boldsymbol{\Pi} = \underbrace{\begin{pmatrix}\boldsymbol{q}_1 & \boldsymbol{q}_2 & \cdots & \boldsymbol{q}_m\end{pmatrix}}_{\boldsymbol{Q}}\underbrace{\begin{pmatrix}R_{1,1} & R_{1,2} & \cdots & R_{1,m} \\ 0 & R_{2,2} & \cdots & R_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_{m,m} \\ \end{pmatrix}}_{\boldsymbol{R}}\end{equation}

where $\boldsymbol{\Pi}$ is a column permutation matrix. Because at each step we select the column with the largest error (norm), we can show that for any $k$, the first column of the submatrix $\boldsymbol{R}_{[k-1:,k-1:]}$ has the largest norm, which is no smaller than the norm of any of the remaining columns:

\begin{equation}R_{k,k}^2 \geq \sum_{i=k}^j R_{i,j}^2,\quad \forall j = k,k+1,\cdots,m\end{equation}

From this we can further derive $|R_{1,1}|\geq |R_{2,2}| \geq \cdots\geq |R_{m,m}|$. These properties give us confidence that, if we want a rank-$r$ approximation of $\boldsymbol{M}\boldsymbol{\Pi}$, keeping only the first $r$ rows of $\boldsymbol{R}$ should be a good choice, i.e.,

\begin{equation}\boldsymbol{M}\boldsymbol{\Pi} = \boldsymbol{Q}\boldsymbol{R} \approx \boldsymbol{Q}_{[:,:r]}\boldsymbol{R}_{[:r,:]}=\boldsymbol{Q}_{[:,:r]}\big[\boldsymbol{R}_{[:r,:r]},\boldsymbol{R}_{[:r,r:]}\big]=\boldsymbol{Q}_{[:,:r]}\boldsymbol{R}_{[:r,:r]}\big[\boldsymbol{I}_r,\boldsymbol{R}_{[:r,:r]}^{-1}\boldsymbol{R}_{[:r,r:]}\big]\end{equation}

Note that we previously agreed that slicing takes priority over matrix inversion, so here $\boldsymbol{R}_{[:r,:r]}^{-1}$ means $(\boldsymbol{R}_{[:r,:r]})^{-1}$. It's not hard to see that $\boldsymbol{Q}_{[:,:r]}\boldsymbol{R}_{[:r,:r]}$ is in fact exactly the $r$ columns of matrix $\boldsymbol{M}$, so the equation above actually gives us an ID approximation:

\begin{equation}\boldsymbol{M} \approx \boldsymbol{C}\boldsymbol{Z},\quad \boldsymbol{C}=\boldsymbol{Q}_{[:,:r]}\boldsymbol{R}_{[:r,:r]},\quad \boldsymbol{Z}=\big[\boldsymbol{I}_r,\boldsymbol{R}_{[:r,:r]}^{-1}\boldsymbol{R}_{[:r,r:]}\big]\boldsymbol{\Pi}^{\top}\end{equation}

The above is the ID-solving algorithm based on column-pivoting QR decomposition, and it's also the algorithm built into SciPy (with rand=False). Note that this algorithm cannot guarantee $|z_{i,j}| \leq 1$ or $|z_{i,j}| \leq 2$, but feedback from many references suggests that in practice $|z_{i,j}| > 2$ almost never occurs, so this counts as a fairly good solving algorithm. In addition, SciPy also has column-pivoting QR decomposition built in, which can be enabled by setting pivoting=True in scipy.linalg.qr.

Randomized Solving

Each orthogonalization step of column-pivoting QR decomposition requires scanning through all remaining vectors to find the one with the largest error, which is often unacceptable when $m$ is large. On the other hand, if $n$ is large, the cost of computing norms and inner products also becomes high. This is where randomized algorithms come in — they try to reduce the effective value of $n$ or $m$ to lower the computational complexity.

Let's first look at the idea of reducing $n$, that is, reducing the dimensionality of each column vector of $\boldsymbol{M}$. A common method is random projection, which is exactly analogous to the "JL Lemma" introduced in The Amazing Johnson-Lindenstrauss Lemma: Theory. Specifically, suppose $\boldsymbol{\Omega}\in\mathbb{R}^{d\times n}$ is some random projection matrix (where $d\ll n$), whose entries are drawn i.i.d. from some distribution such as $\mathcal{N}(0,1/n)$. We then perform column-pivoting QR decomposition on the smaller matrix $\boldsymbol{\Omega}\boldsymbol{M}\in\mathbb{R}^{d\times m}$ to determine the positions of the $r$ selected columns. For a more detailed treatment, see Randomized algorithms for pivoting and for computing interpolatory and CUR factorizations.

Based on the author's limited research, SciPy's randomized algorithm for solving ID follows a similar idea, but replaces the randomly sampled matrix with a more structured "Subsampled Randomized Fourier Transform" (SRFT), so that the computation in step $\boldsymbol{\Omega}\boldsymbol{M}$ can be reduced from $\mathcal{O}(mnd)$ to $\mathcal{O}(mn\log d)$. However, the author is not familiar with the implementation details of SRFT or of SciPy, so interested readers may consult Enabling very large-scale matrix computations via randomization, A brief introduction to Randomized Linear Algebra, and other resources for further study.

Another reason for not delving deeply into SRFT and other sophisticated random-projection methods is that the paper Efficient Algorithms for Constructing an Interpolative Decomposition found that a much simpler column-sampling approach often yields better results, and it's also very easy to understand: randomly sample $k > r$ columns from $\boldsymbol{M}$, then use column-pivoting QR decomposition to select $r$ columns from these $k$ columns as $\boldsymbol{C}$, and finally solve for $\boldsymbol{Z}$ based on $\boldsymbol{C}$. This reduces the size of the matrix used for column-pivoting QR decomposition from $n\times m$ to $n\times k$.

Experiments show that this simple idea at most slightly increases the risk of $|z_{i,j}| > 2$ on a few individual tasks, but it has a clear advantage in terms of error:

Comparison of the maximum absolute value of the Z matrix between column sampling (Optim-RID) and SciPy's built-in algorithm (SciPy-RID)Comparison of error between column sampling (Optim-RID) and SciPy's built-in algorithm (SciPy-RID)Comparison of efficiency between column sampling (Optim-RID) and SciPy's built-in algorithm (SciPy-RID)

Improving Accuracy

There's a perhaps surprising result worth noting in the table above: the randomly column-sampled Optim-RID not only outperforms SciPy-RID (also a randomized algorithm) in terms of error, but on some individual tasks it even outperforms the deterministic algorithms SciPy-ID and Optim-ID (which are mathematically equivalent, both being based on the full column-pivoting QR decomposition, differing only in implementation efficiency).

This seemingly counterintuitive phenomenon actually reveals a fact: although column-pivoting QR decomposition can serve as a decent baseline for ID, its ability to choose a good basis may not be much better than random selection at all — the main role of column-pivoting QR decomposition is really just to guarantee, with high probability, that $|z_{i,j}| < 2$ holds. This is actually not hard to understand. Take $r=1$ as an example: in this case column-pivoting QR decomposition simply returns the column with the largest norm. But is the column with the largest norm necessarily a good basis (i.e., one that minimizes the reconstruction error)? Clearly not — a good basis vector should point in the direction that most of the vectors collectively align with, and having the largest norm doesn't capture that at all.

For ID, column-pivoting QR decomposition is essentially a greedy algorithm: it greedily reduces the problem of selecting $r$ columns into a recursive sequence of selecting one column $1$ at a time. When $r=1$, in scenarios where $m$ is not too large or where high precision is required, it is computationally feasible to solve exactly via enumeration:

\begin{equation}\mathop{\text{argmin}}_i \sum_{j=1}^m \left\Vert\boldsymbol{m}_j - \frac{(\boldsymbol{m}_j^{\top} \boldsymbol{m}_i)\boldsymbol{m}_i}{\Vert\boldsymbol{m}_i\Vert^2}\right\Vert^2\end{equation}

That is, we enumerate all choices of $\boldsymbol{m}_i$, project all remaining columns onto $\boldsymbol{m}_i$ to compute the total error, and choose the $\boldsymbol{m}_i$ that minimizes this total error, at a complexity proportional to $m^2$. If we replace the "choose the column with the largest norm" step of column-pivoting QR decomposition, at each step, with "choose the option that minimizes the total error" as above, we can find a better basis and thereby achieve lower reconstruction error (naturally at the cost of higher complexity, and with even less guarantee that $|z_{i,j}| < 2$ holds).

Overall, because exact solving is NP-hard, there is a huge variety of approaches to solving ID, and the ones listed above are just a small sample. Interested readers can dig deeper by searching around keywords like Randomized Linear Algebra and Column Subset Selection. It's worth particularly noting that Randomized Linear Algebra, which aims to accelerate various matrix computations via randomized methods, has itself become a rich subject in its own right. The randomized ID discussed in this article, and the sampling-based CR approximation from the previous article, are both classic examples from this field.

Summary

This article introduced ID (Interpolative Decomposition), which approximates the original matrix by selecting a number of its columns to serve as a "skeleton." It is a low-rank decomposition with a specific structure whose geometric meaning is relatively more intuitive. Its central difficulty lies in the selection of columns, which is fundamentally an NP-hard discrete optimization problem.