低秩近似之路（三）：CR

Author: [苏剑林]

Link: [https://zhuanlan.zhihu.com/p/3977810830]

最佳排版请看原博客：

低秩近似之路（三）：CR - 科学空间|Scientific Spaces (opens new window)在https://kexue.fm/archives/10407 (opens new window)中，我们证明了SVD可以给出任意矩阵的最优低秩近似。那里的最优近似是无约束的，也就是说SVD给出的结果只管误差上的最小，不在乎矩阵的具体结构，而在很多应用场景中，出于可解释性或者非线性处理等需求，我们往往希望得到具有某些特殊结构的近似分解。

因此，从这篇文章开始，我们将探究一些具有特定结构的低秩近似，而本文将聚焦于其中的CR近似（Column-Row Approximation），它提供了加速矩阵乘法运算的一种简单方案。

问题背景

矩阵的最优$r$ 近似的一般提法是

${\text{argmin}}_{\text{rank}(\tilde{\mathit{M}})\le r}\Vert \tilde{\mathit{M}}-\mathit{M}{\Vert }_F^2$ 其中$\mathit{M},\tilde{\mathit{M}}\in {\mathbb{R}}^{n\times m},r<min(n,m)$。在前两篇文章中，我们已经讨论了两种情况：

1、如果$\tilde{\mathit{M}}$ 再有其他约束，那么$\tilde{\mathit{M}}$ 最优解就是${\mathit{U}}_{[:,:r]}{\mathbf{\Sigma }}_{[:r,:r]}{\mathit{V}}_{[:,:r]}^{\mathrm{\top }}$，其中$\mathit{M}=\mathit{U}\mathbf{\Sigma }{\mathit{V}}^{\mathrm{\top }}$ $\mathit{M}$ 奇异值分解（SVD）；

2、如果约定$\tilde{\mathit{M}}=\mathit{A}\mathit{B}$（$\mathit{A}\in {\mathbb{R}}^{n\times r},\mathit{B}\in {\mathbb{R}}^{r\times m}$），且$\mathit{A}$（或$\mathit{B}$）已经给定，那么$\tilde{\mathit{M}}$ 最优解是$\mathit{A}{\mathit{A}}^{†}\mathit{M}$（或$\mathit{M}{\mathit{B}}^{†}\mathit{B}$），这里的${}^{†}$ “https://kexue.fm/archives/10366 (opens new window)”。

这两个结果都有很广泛的应用，但它们都没有显式地引入$\tilde{\mathit{M}}$ $\mathit{M}$ 构上的关联，这就导致了很难直观地看到$\tilde{\mathit{M}}$ $\mathit{M}$ 关联，换言之$\tilde{\mathit{M}}$ 可解释性不强。 此外，如果目标中包含非线性运算如$\varphi (\mathit{X}\mathit{W})$，通常也不允许我们使用任意实投影矩阵来降维，而是要求使用“选择矩阵（Selective Matrix）”，比如$\varphi (\mathit{X}\mathit{W})\mathit{S}=\varphi (\mathit{X}\mathit{W}\mathit{S})$ 于任意矩阵$\mathit{S}$ 是恒成立的，但对于选择矩阵$\mathit{S}$ 恒成立的。 所以，接下来我们关注选择矩阵约束下的低秩近似。具体来说，我们有$\mathit{X}\in {\mathbb{R}}^{n\times l},\mathit{Y}\in {\mathbb{R}}^{l\times m}$，然后选定$\mathit{M}=\mathit{X}\mathit{Y}$，我们的任务是从$\mathit{X}$ 选出$r$ 、从$\mathit{Y}$ 选出相应的$r$ 来构建$\tilde{\mathit{M}}$，即 ${\text{argmin}}_S\Vert \underset{\mathit{C}}{\underbrace{{\mathit{X}}_{[:,S]}}}\underset{\mathit{R}}{\underbrace{{\mathit{Y}}_{[S,:]}}}-\mathit{X}\mathit{Y}{\Vert }_F^2\text{s.t.}S\subset \{0,1,\cdots ,l-1\},|S|=r$

这里的$S$ 以理解为slice，即按照Python的切片规则来理解，我们称${\mathit{X}}_{[:,S]}{\mathit{Y}}_{[S,:]}$ $\mathit{X}\mathit{Y}$ “CR近似”。注意这种切片结果也可以用选择矩阵来等价描述，假设${\mathit{X}}_{[:,S]}$ 第$1,2,\cdots ,r$ 分别为$\mathit{X}$ 第$s_1,s_2,\cdots ,s_r$ ，那么可以定义选择矩阵$\mathit{S}\in \{0,1{\}}^{l\times r}$：

$S_{i,j}=\{\begin{array}{rlr}& 1,& i=s_j\\ & 0,& i\ne s_j\end{array}$

即$\mathit{S}$ 第$j$ 的第$s_j$ 元素为1，其余都为0，这样一来就有${\mathit{X}}_{[:,S]}=\mathit{X}\mathit{S}$ 及${\mathit{Y}}_{[S,:]}={\mathit{S}}^{\mathrm{\top }}\mathit{Y}$。

初步近似

如果我们将$\mathit{X},\mathit{Y}$ 别表示成

$\mathit{X}=({\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_l),\mathit{Y}=\left(\begin{array}{c}{\mathit{y}}_1^{\mathrm{\top }}\\ {\mathit{y}}_2^{\mathrm{\top }}\\ ⋮\\ {\mathit{y}}_l^{\mathrm{\top }}\end{array}\right)$

其中${\mathit{x}}_i\in {\mathbb{R}}^{n\times 1},{\mathit{y}}_i\in {\mathbb{R}}^{m\times 1}$ 是列向量，那么$\mathit{X}\mathit{Y}$ 以写成

$\mathit{X}\mathit{Y}=\sum _{i=1}^l{\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}$

而找$\mathit{X}\mathit{Y}$ 最优CR近似则可以等价地写成

${\text{argmin}}_{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_l\in \{0,1\}}{\Vert \sum _{i=1}^l{\lambda }_i{\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}-\sum _{i=1}^l{\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}\Vert }_F^2\text{s.t.}\sum _{i=1}^l{\lambda }_i=r$

我们知道，矩阵的$F$ 数实际上就是将矩阵展平成向量来算模长，所以这个优化问题本质上就相当于给定$l$ 向量${\mathit{v}}_1,{\mathit{v}}_2,\cdots ,{\mathit{v}}_l\in {\mathbb{R}}^d$，求

${\text{argmin}}_{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_l\in \{0,1\}}{\Vert \sum _{i=1}^l{\lambda }_i{\mathit{v}}_i-\sum _{i=1}^l{\mathit{v}}_i\Vert }^2\text{s.t.}\sum _{i=1}^l{\lambda }_i=r$

其中${\mathit{v}}_i=\text{vec}({\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }})$，$d=nm$。记${\gamma }_i=1-{\lambda }_i$，那么可以进一步简化成

${\text{argmin}}_{{\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_l\in \{0,1\}}{\Vert \sum _{i=1}^l{\gamma }_i{\mathit{v}}_i\Vert }^2\text{s.t.}\sum _{i=1}^l{\gamma }_i=l-r$

如果笔者没有理解错，这个优化问题的精确求解是NP-Hard的，所以一般情况下只能寻求近似算法。一个可精确求解的简单例子是${\mathit{v}}_1,{\mathit{v}}_2,\cdots ,{\mathit{v}}_l$ 两垂直，此时

${\Vert \sum _{i=1}^l{\gamma }_i{\mathit{v}}_i\Vert }^2=\sum _{i=1}^l{\gamma }_i^2\Vert {\mathit{v}}_i{\Vert }^2$

所以它的最小值就是最小的$l-r$ $\Vert {\mathit{v}}_i{\Vert }^2$ 和，即让模长最小的$l-r$ ${\mathit{v}}_i$ ${\gamma }_i$ 于1，剩下的${\gamma }_i$ 等于0。当两两正交的条件不严格成立时，我们依然可以将选择模长最小的$l-r$ ${\mathit{v}}_i$ 为一个近似解。回到原始的CR近似问题上，我们有$\Vert {\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}{\Vert }_F=\Vert {\mathit{x}}_i\Vert \Vert {\mathit{y}}_i\Vert$，所以$\mathit{X}\mathit{Y}$ 最优CR近似的一个baseline，就是保留$\mathit{X}$ 列向量与$\mathit{Y}$ 应的行向量模长乘积最大的$r$ 列/行向量。

采样视角

有一些场景允许我们将式$\text{(???)}$ 宽为

${\text{argmin}}_{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_l\in \mathbb{R}}{\Vert \sum _{i=1}^l{\lambda }_i{\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}-\sum _{i=1}^l{\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}\Vert }_F^2\text{s.t.}\sum _{i=1}^l\mathrm{\#}[{\lambda }_i\ne 0]=r$

其中$\mathrm{\#}[{\lambda }_i\ne 0]$ 示${\lambda }_i\ne 0$ 输出1，否则输出0。这个宽松版本其实就是将CR近似的形式从$\mathit{C}\mathit{R}$ 展成$\mathit{C}\mathbf{\Lambda }\mathit{R}$，其中$\mathbf{\Lambda }\in {\mathbb{R}}^{r\times r}$ 对角阵，即允许我们调整对角阵$\mathbf{\Lambda }\in {\mathbb{R}}^{r\times r}$ 达到更高的精度。相应地，式$\text{(???)}$ 为

${\text{argmin}}_{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_l\in \mathbb{R}}{\Vert \sum _{i=1}^l{\lambda }_i{\mathit{v}}_i-\sum _{i=1}^l{\mathit{v}}_i\Vert }^2\text{s.t.}\sum _{i=1}^l\mathrm{\#}[{\lambda }_i\ne 0]=r$

这样放宽之后，我们可以从采样视角来看待这个问题。首先我们引入任意$l$ 分布$\mathit{p}=(p_1,p_2,\cdots ,p_l)$，然后我们就可以写出

$\sum _{i=1}^l{\mathit{v}}_i=\sum _{i=1}^l{p}_i\times \frac{{\mathit{v}}_i}{p_i}={\mathbb{E}}_{i\sim \mathit{p}}\left[\frac{{\mathit{v}}_i}{p_i}\right]$

也就是说，${\mathit{v}}_i/p_i$ 数学期望正好是我们要逼近的目标，所以我们可以通过从$\mathit{p}$ 布独立重复采样来构建近似：

$\sum _{i=1}^l{\mathit{v}}_i={\mathbb{E}}_{i\sim \mathit{p}}\left[\frac{{\mathit{v}}_i}{p_i}\right]\approx \frac{1}{r}\sum _{j=1}^r\frac{{\mathit{v}}_{s_j}}{p_{s_j}},s_1,s_2,\cdots ,s_r\sim \mathit{p}$

这意味着当$i$ $s_1,s_2,\cdots ,s_r$ 一时有${\lambda }_i=(r{p}_i{)}^{-1}$，否则${\lambda }_i=0$。可能有读者疑问为什么要独立重复采样，而不是更符合逼近需求的不放回采样呢？无他，纯粹是因为独立重复采样使得后面的分析更简单。

到目前为止，我们的理论结果跟分布$\mathit{p}$ 选择无关，也就是说对于任意$\mathit{p}$ 是成立的，这给我们提供了选择最优$\mathit{p}$ 可能性。那如何衡量$\mathit{p}$ 优劣呢？很显然我们希望每次采样估计的误差越小越好，因此可以用采样估计的误差

${\mathbb{E}}_{i\sim \mathit{p}}\left[{\Vert \frac{{\mathit{v}}_i}{p_i}-\sum _{i=1}^l{\mathit{v}}_i\Vert }^2\right]=\left(\sum _{i=1}^l\frac{\Vert {\mathit{v}}_i{\Vert }^2}{p_i}\right)-{\Vert \sum _{i=1}^l{\mathit{v}}_i\Vert }^2$

来比较不同的$\mathit{p}$ 间的优劣。接着利用均值不等式有

$\sum _{i=1}^l\frac{\Vert {\mathit{v}}_i{\Vert }^2}{p_i}=(\sum _{i=1}^l\frac{\Vert {\mathit{v}}_i{\Vert }^2}{p_i}+p_i{Z}^2)-Z^2\ge (\sum _{i=1}^{l}2\Vert {\mathit{v}}_i\Vert Z)-Z^2$

等号在$\Vert {\mathit{v}}_i{\Vert }^2/p_i=p_i{Z}^2$ 取到，由此可得最优的$\mathit{p}$

$p_i^{\ast }=\frac{\Vert {\mathit{v}}_i\Vert }{Z},Z=\sum _{i=1}^l\Vert {\mathit{v}}_i\Vert$

对应的误差为

${\mathbb{E}}_{i\sim \mathit{p}}\left[{\Vert \frac{{\mathit{v}}_i}{p_i}-\sum _{i=1}^l{\mathit{v}}_i\Vert }^2\right]={(\sum _{i=1}^l\Vert {\mathit{v}}_i\Vert )}^2-{\Vert \sum _{i=1}^l{\mathit{v}}_i\Vert }^2$

最优的$p_i$ 好正比于$\Vert {\mathit{v}}_i\Vert$，所以概率最大的$r$ ${\mathit{v}}_i$ 正是模长最大的$r$ ${\mathit{v}}_i$，这就跟上一节的近似联系起来了。该结果出自2006年的论文https://www.stat.berkeley.edu/~mmahoney/pubs/matrix1_SICOMP.pdf (opens new window)，初衷是加速矩阵乘法，它表明只要按照$p_i\propto \Vert {\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }}{\Vert }_F=\Vert {\mathit{x}}_i\Vert \Vert {\mathit{y}}_i\Vert$ 采样$\mathit{X},\mathit{Y}$ 应的列/行，并乘以$(r{p}_i{)}^{-1/2}$，就可以得到$\mathit{X}\mathit{Y}$ 一个CR近似，从而将乘法复杂度从$\mathcal{O}(lmn)$ 低到$\mathcal{O}(rmn)$。

延伸讨论

不管是按模长排序还是按$p_i\propto \Vert {\mathit{v}}_i\Vert$ 机采样，它们都允许我们在线性复杂度【即$\mathcal{O}(l)$】内构建一个CR近似，这对于实时计算来说当然是很理想的，但由于排序或采样都只依赖于$\Vert {\mathit{v}}_i\Vert$，所以精度只能说一般。如果我们可以接受更高的复杂度，那么如何提高CR近似的精度呢？

我们可以尝试将排序的单位改为$k$ 组。简单起见，假设$k\le l-r$ $l-r$ 一个因数，$l$ 向量${\mathit{v}}_1,{\mathit{v}}_2,\cdots ,{\mathit{v}}_l$ $k$ 的组合数为$C_l^k$，对于每个组合$\{s_1,s_2,\cdots ,s_k\}$ 们都可以算出向量和的模长$\Vert {\mathit{v}}_{s_1}+{\mathit{v}}_{s_2}+\cdots +{\mathit{v}}_{s_k}\Vert$。有了这些数据，我们就可以贪婪地构建$\text{(???)}$ 近似解：

初始化$\mathrm{\Omega }=\{1,2,\cdots ,l\},\mathrm{\Theta }=\{\}$

对于$t=1,2,\cdots ,(l-r)/k$，执行：

$\mathrm{\Theta }=\mathrm{\Theta }\cup \underset{\{s_1,s_2,\cdots ,s_k\}\subset \mathrm{\Omega }}{\text{argmin}}\Vert {\mathit{v}}_{s_1}+{\mathit{v}}_{s_2}+\cdots +{\mathit{v}}_{s_k}\Vert$；

$\mathrm{\Omega }=\mathrm{\Omega }\mathrm{\setminus }\mathrm{\Theta }$；

返回$\mathrm{\Theta }$。

说白了，就是每次都从剩下的向量中挑选和模长最小的$k$ 向量，重复挑选$(l-r)/k$ 即得到$l-r$ 向量，它是按照单个向量模长来排序的自然推广，其复杂度为$\mathcal{O}(C_l^k)$，当$k>1$ $l$ 较大时可能难以承受，这也侧面体现了原问题精确求解的复杂性。

另一个值得思考的问题是如果允许CR近似放宽为$\mathit{C}\mathbf{\Lambda }\mathit{R}$，那么$\mathbf{\Lambda }$ 最优解是什么呢？如果不限定$\mathbf{\Lambda }$ 结构，那么答案可以由伪逆给出

${\mathbf{\Lambda }}^{\ast }={\text{argmin}}_{\mathbf{\Lambda }}\Vert \mathit{C}\mathbf{\Lambda }\mathit{R}-\mathit{X}\mathit{Y}{\Vert }_F^2={\mathit{C}}^{†}\mathit{X}\mathit{Y}{\mathit{R}}^{†}$

如果$\mathbf{\Lambda }$ 须是对角阵呢？那可以先将问题重述为给定$\{{\mathit{u}}_1,{\mathit{u}}_2,\cdots ,{\mathit{u}}_r\}\subset \{{\mathit{v}}_1,{\mathit{v}}_2,\cdots ,{\mathit{v}}_l\}$，求

${\text{argmin}}_{{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_r}{\Vert \sum _{i=1}^r{\lambda }_i{\mathit{u}}_i-\sum _{i=1}^l{\mathit{v}}_i\Vert }^2$

我们记$\mathit{U}=({\mathit{u}}_1,{\mathit{u}}_2,\cdots ,{\mathit{u}}_r),\mathit{V}=({\mathit{v}}_1,{\mathit{v}}_2,\cdots ,{\mathit{v}}_l),\mathit{\lambda }=({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_r{)}^{\mathrm{\top }}$，那么优化目标可以写成

${\text{argmin}}_{\mathit{\lambda }}{\Vert \mathit{U}\mathit{\lambda }-\mathit{V}{\mathbf{1}}_{l\times 1}\Vert }^2$

这同样可以通过伪逆写出最优解

${\mathit{\lambda }}^{\ast }={\mathit{U}}^{†}\mathit{V}{\mathbf{1}}_{l\times 1}=({\mathit{U}}^{\mathrm{\top }}\mathit{U}{)}^{-1}{\mathit{U}}^{\mathrm{\top }}\mathit{V}{\mathbf{1}}_{l\times 1}$

最后一个等号假设了${\mathit{U}}^{\mathrm{\top }}\mathit{U}$ 逆，这通常能满足，如果不满足的话$({\mathit{U}}^{\mathrm{\top }}\mathit{U}{)}^{-1}$ $({\mathit{U}}^{\mathrm{\top }}\mathit{U}{)}^{†}$ 行。

现在的问题是直接套用上式的话对原始问题来说计算量太大，因为${\mathit{v}}_i=\text{vec}({\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }})$，即${\mathit{v}}_i$ $mn$ 向量，所以$\mathit{V}$ 小为$mn\times l$、$\mathit{U}$ 小为$mn\times r$，这在$m,n$ 大时比较难受。利用${\mathit{v}}_i=\text{vec}({\mathit{x}}_i{\mathit{y}}_i^{\mathrm{\top }})$ 帮我们进一步化简上式。不妨设${\mathit{u}}_i=\text{vec}({\mathit{c}}_i{\mathit{r}}_i^{\mathrm{\top }})$，那么

$\begin{array}{rl}({\mathit{U}}^{\mathrm{\top }}\mathit{V}{)}_{i,j}=& ⟨{\mathit{c}}_i{\mathit{r}}_i^{\mathrm{\top }},{\mathit{x}}_j{\mathit{y}}_j^{\mathrm{\top }}{⟩}_F=\text{Tr}({\mathit{r}}_i{\mathit{c}}_i^{\mathrm{\top }}{\mathit{x}}_j{\mathit{y}}_j^{\mathrm{\top }})=({\mathit{c}}_i^{\mathrm{\top }}{\mathit{x}}_j)({\mathit{r}}_i^{\mathrm{\top }}{\mathit{y}}_j)\\ =& [({\mathit{C}}^{\mathrm{\top }}\mathit{X})\otimes (\mathit{R}{\mathit{Y}}^{\mathrm{\top }}){]}_{i,j}\end{array}$

即${\mathit{U}}^{\mathrm{\top }}\mathit{V}=({\mathit{C}}^{\mathrm{\top }}\mathit{X})\otimes (\mathit{R}{\mathit{Y}}^{\mathrm{\top }}),{\mathit{U}}^{\mathrm{\top }}\mathit{U}=({\mathit{C}}^{\mathrm{\top }}\mathit{C})\otimes (\mathit{R}{\mathit{R}}^{\mathrm{\top }})$，这里的$\otimes$ https://en.wikipedia.org/wiki/Hadamard_product_(matrices) (opens new window)，这样恒等变换之后${\mathit{U}}^{\mathrm{\top }}\mathit{V}$ ${\mathit{U}}^{\mathrm{\top }}\mathit{U}$ 计算量就降低了。

文章小结

本文介绍了矩阵乘法的CR近似，这是一种具有特定行列结构的低秩近似，相比由SVD给出的最优低秩近似，CR近似具有更直观的物理意义以及更好的可解释性。