The Path to Low-Rank Approximation (Part 3): CR Approximation

Problem Background#

The general formulation for optimal rank-$r$ approximation of a matrix is:

where $\boldsymbol{M},\tilde{\boldsymbol{M}}\in\mathbb{R}^{n\times m},r < \min(n,m)$. In the first two articles of this series, we have discussed two scenarios:

1. If $\tilde{\boldsymbol{M}}$ has no additional constraints, then the optimal solution for $\tilde{\boldsymbol{M}}$ is $\boldsymbol{U}_{[:,:r]}\boldsymbol{\Sigma}_{[:r,:r]}\boldsymbol{V}_{[:,:r]}^{\top}$, where $\boldsymbol{M}=\boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^{\top}$ is the Singular Value Decomposition (SVD) of $\boldsymbol{M}$;

2. If we require $\tilde{\boldsymbol{M}}=\boldsymbol{A}\boldsymbol{B}$ (with $\boldsymbol{A}\in\mathbb{R}^{n\times r},\boldsymbol{B}\in\mathbb{R}^{r\times m}$) and $\boldsymbol{A}$ (or $\boldsymbol{B}$) is already given, then the optimal solution for $\tilde{\boldsymbol{M}}$ is $\boldsymbol{A} \boldsymbol{A}^{\dagger} \boldsymbol{M}$ (or $\boldsymbol{M} \boldsymbol{B}^{\dagger} \boldsymbol{B}$), where ${}^\dagger$ denotes the "pseudo-inverse".

Both results have extensive applications, but neither explicitly introduces structural relationships between $\tilde{\boldsymbol{M}}$ and $\boldsymbol{M}$. This makes it difficult to intuitively understand the connection between $\tilde{\boldsymbol{M}}$ and $\boldsymbol{M}$—in other words, $\tilde{\boldsymbol{M}}$ lacks strong interpretability.

Here $S$ can be understood as a slice (following Python slicing conventions). We refer to $\boldsymbol{X}_{[:,S]}\boldsymbol{Y}_{[S,:]}$ as the "CR approximation" of $\boldsymbol{X}\boldsymbol{Y}$. Note that this slicing operation can also be equivalently described using selection matrices. Suppose the first $1,2,\cdots,r$ columns of $\boldsymbol{X}_{[:,S]}$ correspond to columns $s_1,s_2,\cdots,s_r$ of $\boldsymbol{X}$, respectively. We can then define a selection matrix $\boldsymbol{S}\in\{0,1\}^{l\times r}$:

That is, the $s_j$-th element of the $j$-th column of $\boldsymbol{S}$ equals 1, while all others equal 0. This yields $\boldsymbol{X}_{[:,S]}=\boldsymbol{X}\boldsymbol{S}$ and $\boldsymbol{Y}_{[S,:]}=\boldsymbol{S}^{\top} \boldsymbol{Y}$.

Preliminary Approximation#

If we express $\boldsymbol{X},\boldsymbol{Y}$ respectively as:

where $\boldsymbol{x}_i\in\mathbb{R}^{n\times 1},\boldsymbol{y}_i\in\mathbb{R}^{m\times 1}$ are column vectors, then $\boldsymbol{X}\boldsymbol{Y}$ can be written as:

Finding the optimal CR approximation of $\boldsymbol{X}\boldsymbol{Y}$ can be equivalently expressed as:

A simple case that admits an exact solution occurs when $\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l$ are mutually orthogonal. In this scenario:

Therefore, its minimum value is the sum of the smallest $l-r$ values of $\Vert\boldsymbol{v}_i\Vert^2$. That is, set $\gamma_i = 1$ for the $l-r$ vectors $\boldsymbol{v}_i$ with the smallest norms, and $\gamma_i = 0$ for the remaining ones. When the mutual orthogonality condition is not strictly satisfied, we can still use "selecting the $l-r$ vectors $\boldsymbol{v}_i$ with the smallest norms" as an approximate solution.

Returning to the original CR approximation problem, we have $\Vert\boldsymbol{x}_i\boldsymbol{y}_i^{\top}\Vert_F = \Vert\boldsymbol{x}_i\Vert \Vert \boldsymbol{y}_i\Vert$. Thus, a baseline approach for the optimal CR approximation of $\boldsymbol{X}\boldsymbol{Y}$ is to retain the $r$ column/row pairs with the largest products of column vector norms from $\boldsymbol{X}$ and the corresponding row vector norms from $\boldsymbol{Y}$.

Sampling Perspective#

In some scenarios, we can relax Equation $\eqref{eq:xy-l-k}$ to:

where $\#[\lambda_i\neq 0]$ outputs 1 if $\lambda_i\neq 0$, and 0 otherwise. This relaxed version essentially extends the CR approximation form from $\boldsymbol{C}\boldsymbol{R}$ to $\boldsymbol{C}\boldsymbol{\Lambda}\boldsymbol{R}$, where $\boldsymbol{\Lambda}\in\mathbb{R}^{r\times r}$ is a diagonal matrix. That is, it allows us to adjust the diagonal matrix $\boldsymbol{\Lambda}$ to achieve higher accuracy. Correspondingly, Equation $\eqref{eq:v-l-k}$ becomes:

Up to this point, our theoretical results are independent of the choice of distribution $\boldsymbol{p}$—meaning they hold for any $\boldsymbol{p}$. This provides us with the possibility to select an optimal $\boldsymbol{p}$. How do we measure the quality of $\boldsymbol{p}$? Clearly, we want the error of each sampling estimate to be as small as possible. Therefore, we can use the expected squared error of the sampling estimate:

to compare the effectiveness of different $\boldsymbol{p}$ distributions. Using the inequality of arithmetic and geometric means (AM-GM inequality), we have:

Equality is achieved when $\Vert\boldsymbol{v}_i\Vert^2 / p_i = p_i Z^2$, yielding the optimal $\boldsymbol{p}$:

The corresponding error becomes:

Extended Discussions#

Whether sorting by norm or random sampling according to $p_i\propto \Vert\boldsymbol{v}_i\Vert$, both methods allow us to construct a CR approximation with linear complexity [i.e., $\mathcal{O}(l)$]. This is certainly ideal for real-time computation. However, since both sorting and sampling depend solely on $\Vert\boldsymbol{v}_i\Vert$, the accuracy can only be considered moderate. If we can accept higher complexity, how can we improve the accuracy of CR approximation?

Another question worth considering is: if we allow the CR approximation to be relaxed to $\boldsymbol{C}\boldsymbol{\Lambda}\boldsymbol{R}$, what is the optimal solution for $\boldsymbol{\Lambda}$? If no structure is imposed on $\boldsymbol{\Lambda}$, the answer can be given by the pseudo-inverse:

What if $\boldsymbol{\Lambda}$ must be a diagonal matrix? Then we can first reformulate the problem as: given $\{\boldsymbol{u}_1,\boldsymbol{u}_2,\cdots,\boldsymbol{u}_r\}\subset\{\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l\}$, find:

Denote $\boldsymbol{U} = (\boldsymbol{u}_1,\boldsymbol{u}_2,\cdots,\boldsymbol{u}_r), \boldsymbol{V} = (\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l), \boldsymbol{\lambda}=(\lambda_1,\lambda_2,\cdots,\lambda_r)^{\top}$. The optimization objective can then be written as:

This can also be expressed using the pseudo-inverse to obtain the optimal solution:

The last equality assumes $\boldsymbol{U}^{\top}\boldsymbol{U}$ is invertible, which is typically satisfied. If not, replace $(\boldsymbol{U}^{\top}\boldsymbol{U})^{-1}$ with $(\boldsymbol{U}^{\top}\boldsymbol{U})^{\dagger}$.

The current issue is that directly applying this formula results in excessive computational cost for the original problem, because $\boldsymbol{v}_i = \text{vec}(\boldsymbol{x}_i \boldsymbol{y}_i^{\top})$—meaning $\boldsymbol{v}_i$ is an $mn$-dimensional vector. Consequently, $\boldsymbol{V}$ has size $mn\times l$ and $\boldsymbol{U}$ has size $mn\times r$, which becomes problematic when $m,n$ are large.

Utilizing the relationship $\boldsymbol{v}_i = \text{vec}(\boldsymbol{x}_i \boldsymbol{y}_i^{\top})$ helps us further simplify the above expression. Let $\boldsymbol{u}_i = \text{vec}(\boldsymbol{c}_i \boldsymbol{r}_i^{\top})$. Then:

That is, $\boldsymbol{U}^{\top}\boldsymbol{V}=(\boldsymbol{C}^{\top}\boldsymbol{X})\odot (\boldsymbol{R}\boldsymbol{Y}^{\top})$ and $\boldsymbol{U}^{\top}\boldsymbol{U}=(\boldsymbol{C}^{\top}\boldsymbol{C})\odot (\boldsymbol{R}\boldsymbol{R}^{\top})$, where $\odot$ denotes the Hadamard product. Through these identity transformations, the computational cost of $\boldsymbol{U}^{\top}\boldsymbol{V}$ and $\boldsymbol{U}^{\top}\boldsymbol{U}$ is significantly reduced.

Article Summary#