《低秩近似之路(二):SVD》中,我们证明了SVD可以给出任意矩阵的最优低秩近似。那里的最优近似是无约束的,也就是说SVD给出的结果只管误差上的最小,不在乎矩阵的具体结构,而在很多应用场景中,出于可解释性或者非线性处理等需求,我们往往希望得到具有某些特殊结构的近似分解。

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

问题背景 #

矩阵的最优r秩近似的一般提法是
\begin{equation}\mathop{\text{argmin}}_{\text{rank}(\tilde{\boldsymbol{M}})\leq r}\Vert \tilde{\boldsymbol{M}} - \boldsymbol{M}\Vert_F^2\label{eq:loss-m2}\end{equation}
其中\boldsymbol{M},\tilde{\boldsymbol{M}}\in\mathbb{R}^{n\times m},r < \min(n,m)。在前两篇文章中,我们已经讨论了两种情况:

1、如果\tilde{\boldsymbol{M}}不再有其他约束,那么\tilde{\boldsymbol{M}}的最优解就是\boldsymbol{U}_{[:,:r]}\boldsymbol{\Sigma}_{[:r,:r]}\boldsymbol{V}_{[:,:r]}^{\top},其中\boldsymbol{M}=\boldsymbol{U}\boldsymbol{\Sigma} \boldsymbol{V}^{\top}\boldsymbol{M}的奇异值分解(SVD);

2、如果约定\tilde{\boldsymbol{M}}=\boldsymbol{A}\boldsymbol{B}\boldsymbol{A}\in\mathbb{R}^{n\times r},\boldsymbol{B}\in\mathbb{R}^{r\times m}),且\boldsymbol{A}(或\boldsymbol{B})已经给定,那么\tilde{\boldsymbol{M}}的最优解是\boldsymbol{A} \boldsymbol{A}^{\dagger} \boldsymbol{M}(或\boldsymbol{M} \boldsymbol{B}^{\dagger} \boldsymbol{B}),这里的{}^\dagger是“伪逆”。

这两个结果都有很广泛的应用,但它们都没有显式地引入\tilde{\boldsymbol{M}}\boldsymbol{M}结构上的关联,这就导致了很难直观地看到\tilde{\boldsymbol{M}}\boldsymbol{M}的关联,换言之\tilde{\boldsymbol{M}}的可解释性不强。

此外,如果目标中包含非线性运算如\phi(\boldsymbol{X}\boldsymbol{W}),通常也不允许我们使用任意实投影矩阵来降维,而是要求使用“选择矩阵(Selective Matrix)”,比如\phi(\boldsymbol{X}\boldsymbol{W})\boldsymbol{S} = \phi(\boldsymbol{X}\boldsymbol{W}\boldsymbol{S})对于任意矩阵\boldsymbol{S}不是恒成立的,但对于选择矩阵\boldsymbol{S}是恒成立的。

所以,接下来我们关注选择矩阵约束下的低秩近似。具体来说,我们有\boldsymbol{X}\in\mathbb{R}^{n\times l},\boldsymbol{Y}\in\mathbb{R}^{l\times m},然后选定\boldsymbol{M}=\boldsymbol{X}\boldsymbol{Y},我们的任务是从\boldsymbol{X}中选出r列、从\boldsymbol{Y}中选出相应的r行来构建\tilde{\boldsymbol{M}},即
\begin{equation}\mathop{\text{argmin}}_S\Vert \underbrace{\boldsymbol{X}_{[:,S]}}_{\boldsymbol{C}}\underbrace{\boldsymbol{Y}_{[S,:]}}_{\boldsymbol{R}} - \boldsymbol{X}\boldsymbol{Y}\Vert_F^2\quad\text{s.t.}\quad S\subset \{0,1,\cdots,l-1\}, |S|=r\end{equation}
这里的S可以理解为slice,即按照Python的切片规则来理解,我们称\boldsymbol{X}_{[:,S]}\boldsymbol{Y}_{[S,:]}\boldsymbol{X}\boldsymbol{Y}的“CR近似”。注意这种切片结果也可以用选择矩阵来等价描述,假设\boldsymbol{X}_{[:,S]}的第1,2,\cdots,r列分别为\boldsymbol{X}的第s_1,s_2,\cdots,s_r列,那么可以定义选择矩阵\boldsymbol{S}\in\{0,1\}^{l\times r}
\begin{equation}S_{i,j}=\left\{\begin{aligned}&1, &i = s_j \\ &0, &i\neq s_j\end{aligned}\right.\end{equation}
\boldsymbol{S}的第j列的第s_j个元素为1,其余都为0,这样一来就有\boldsymbol{X}_{[:,S]}=\boldsymbol{X}\boldsymbol{S}以及\boldsymbol{Y}_{[S,:]}=\boldsymbol{S}^{\top} \boldsymbol{Y}

初步近似 #

如果我们将\boldsymbol{X},\boldsymbol{Y}分别表示成
\begin{equation}\boldsymbol{X} = (\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_l),\quad \boldsymbol{Y}=\begin{pmatrix}\boldsymbol{y}_1^{\top} \\ \boldsymbol{y}_2^{\top} \\ \vdots \\ \boldsymbol{y}_l^{\top}\end{pmatrix}\end{equation}
其中\boldsymbol{x}_i\in\mathbb{R}^{n\times 1},\boldsymbol{y}_i\in\mathbb{R}^{m\times 1}都是列向量,那么\boldsymbol{X}\boldsymbol{Y}可以写成
\begin{equation}\boldsymbol{X}\boldsymbol{Y} = \sum_{i=1}^l \boldsymbol{x}_i\boldsymbol{y}_i^{\top}\end{equation}
而找\boldsymbol{X}\boldsymbol{Y}的最优CR近似则可以等价地写成
\begin{equation}\mathop{\text{argmin}}_{\lambda_1,\lambda_2,\cdots,\lambda_l\in\{0,1\}}\left\Vert\sum_{i=1}^l \lambda_i \boldsymbol{x}_i\boldsymbol{y}_i^{\top} - \sum_{i=1}^l\boldsymbol{x}_i\boldsymbol{y}_i^{\top}\right\Vert_F^2\quad\text{s.t.}\quad \sum_{i=1}^l \lambda_i = r\label{eq:xy-l-k}\end{equation}
我们知道,矩阵的F范数实际上就是将矩阵展平成向量来算模长,所以这个优化问题本质上就相当于给定l个向量\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l\in\mathbb{R}^d,求
\begin{equation}\mathop{\text{argmin}}_{\lambda_1,\lambda_2,\cdots,\lambda_l\in\{0,1\}}\left\Vert\sum_{i=1}^l \lambda_i \boldsymbol{v}_i - \sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2\quad\text{s.t.}\quad \sum_{i=1}^l \lambda_i = r\label{eq:v-l-k}\end{equation}
其中\boldsymbol{v}_i = \text{vec}(\boldsymbol{x}_i \boldsymbol{y}_i^{\top})d=nm。记\gamma_i = 1 - \lambda_i,那么可以进一步简化成
\begin{equation}\mathop{\text{argmin}}_{\gamma_1,\gamma_2,\cdots,\gamma_l\in\{0,1\}}\left\Vert\sum_{i=1}^l \gamma_i \boldsymbol{v}_i\right\Vert^2\quad\text{s.t.}\quad \sum_{i=1}^l \gamma_i = l-r\label{eq:v-l-k-0}\end{equation}
如果笔者没有理解错,这个优化问题的精确求解是NP-Hard的,所以一般情况下只能寻求近似算法。一个可精确求解的简单例子是\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l两两垂直,此时
\begin{equation}\left\Vert\sum_{i=1}^l \gamma_i \boldsymbol{v}_i\right\Vert^2 = \sum_{i=1}^l \gamma_i^2 \Vert\boldsymbol{v}_i\Vert^2\end{equation}
所以它的最小值就是最小的l-r\Vert\boldsymbol{v}_i\Vert^2之和,即让模长最小的l-r\boldsymbol{v}_i\gamma_i等于1,剩下的\gamma_i则等于0。当两两正交的条件不严格成立时,我们依然可以将选择模长最小的l-r\boldsymbol{v}_i作为一个近似解。回到原始的CR近似问题上,我们有\Vert\boldsymbol{x}_i\boldsymbol{y}_i^{\top}\Vert_F = \Vert\boldsymbol{x}_i\Vert \Vert \boldsymbol{y}_i\Vert,所以\boldsymbol{X}\boldsymbol{Y}的最优CR近似的一个baseline,就是保留\boldsymbol{X}的列向量与\boldsymbol{Y}对应的行向量模长乘积最大的r个列/行向量。

采样视角 #

有一些场景允许我们将式\eqref{eq:xy-l-k}放宽为
\begin{equation}\mathop{\text{argmin}}_{\lambda_1,\lambda_2,\cdots,\lambda_l\in\mathbb{R}}\left\Vert\sum_{i=1}^l \lambda_i \boldsymbol{x}_i\boldsymbol{y}_i^{\top} - \sum_{i=1}^l\boldsymbol{x}_i\boldsymbol{y}_i^{\top}\right\Vert_F^2\quad\text{s.t.}\quad \sum_{i=1}^l \#[\lambda_i\neq 0] = r\end{equation}
其中\#[\lambda_i\neq 0]表示\lambda_i\neq 0时输出1,否则输出0。这个宽松版本其实就是将CR近似的形式从\boldsymbol{C}\boldsymbol{R}拓展成\boldsymbol{C}\boldsymbol{\Lambda}\boldsymbol{R},其中\boldsymbol{\Lambda}\in\mathbb{R}^{r\times r}是对角阵,即允许我们调整对角阵\boldsymbol{\Lambda}\in\mathbb{R}^{r\times r}来达到更高的精度。相应地,式\eqref{eq:v-l-k}变为
\begin{equation}\mathop{\text{argmin}}_{\lambda_1,\lambda_2,\cdots,\lambda_l\in\mathbb{R}}\left\Vert\sum_{i=1}^l \lambda_i \boldsymbol{v}_i - \sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2\quad\text{s.t.}\quad \sum_{i=1}^l \#[\lambda_i\neq 0] = r\end{equation}

这样放宽之后,我们可以从采样视角来看待这个问题。首先我们引入任意l元分布\boldsymbol{p}=(p_1,p_2,\cdots,p_l),然后我们就可以写出
\begin{equation}\sum_{i=1}^l\boldsymbol{v}_i = \sum_{i=1}^l p_i\times\frac{\boldsymbol{v}_i}{p_i} = \mathbb{E}_{i\sim \boldsymbol{p}} \left[\frac{\boldsymbol{v}_i}{p_i}\right] \end{equation}
也就是说,\boldsymbol{v}_i/p_i的数学期望正好是我们要逼近的目标,所以我们可以通过从\boldsymbol{p}分布独立重复采样来构建近似:
\begin{equation}\sum_{i=1}^l\boldsymbol{v}_i = \mathbb{E}_{i\sim \boldsymbol{p}} \left[\frac{\boldsymbol{v}_i}{p_i}\right] \approx \frac{1}{r}\sum_{j=1}^r \frac{\boldsymbol{v}_{s_j}}{p_{s_j}},\quad s_1,s_2,\cdots,s_r\sim \boldsymbol{p}\end{equation}
这意味着当is_1,s_2,\cdots,s_r之一时有\lambda_i = (r p_i)^{-1},否则\lambda_i=0。可能有读者疑问为什么要独立重复采样,而不是更符合逼近需求的不放回采样呢?无他,纯粹是因为独立重复采样使得后面的分析更简单。

到目前为止,我们的理论结果跟分布\boldsymbol{p}的选择无关,也就是说对于任意\boldsymbol{p}都是成立的,这给我们提供了选择最优\boldsymbol{p}的可能性。那如何衡量\boldsymbol{p}的优劣呢?很显然我们希望每次采样估计的误差越小越好,因此可以用采样估计的误差
\begin{equation}\mathbb{E}_{i\sim \boldsymbol{p}} \left[\left\Vert\frac{\boldsymbol{v}_i}{p_i} - \sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2\right] = \left(\sum_{i=1}^l \frac{\Vert\boldsymbol{v}_i\Vert^2}{p_i}\right) - \left\Vert\sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2 \end{equation}
来比较不同的\boldsymbol{p}之间的优劣。接着利用均值不等式有
\begin{equation}\sum_{i=1}^l \frac{\Vert\boldsymbol{v}_i\Vert^2}{p_i} = \left(\sum_{i=1}^l \frac{\Vert\boldsymbol{v}_i\Vert^2}{p_i} + p_i Z^2\right) - Z^2\geq \left(\sum_{i=1}^l 2\Vert\boldsymbol{v}_i\Vert Z\right) - Z^2\end{equation}
等号在\Vert\boldsymbol{v}_i\Vert^2 / p_i = p_i Z^2时取到,由此可得最优的\boldsymbol{p}
\begin{equation}p_i^* = \frac{\Vert\boldsymbol{v}_i\Vert}{Z},\quad Z = \sum\limits_{i=1}^l \Vert\boldsymbol{v}_i\Vert\end{equation}
对应的误差为
\begin{equation}\mathbb{E}_{i\sim \boldsymbol{p}} \left[\left\Vert\frac{\boldsymbol{v}_i}{p_i} - \sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2\right] = \left(\sum_{i=1}^l \Vert\boldsymbol{v}_i\Vert\right)^2 - \left\Vert\sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2 \end{equation}
最优的p_i正好正比于\Vert\boldsymbol{v}_i\Vert,所以概率最大的r\boldsymbol{v}_i也正是模长最大的r\boldsymbol{v}_i,这就跟上一节的近似联系起来了。该结果出自2006年的论文《Fast Monte Carlo Algorithms for Matrices I: Approximating Matrix Multiplication》,初衷是加速矩阵乘法,它表明只要按照p_i\propto \Vert \boldsymbol{x}_i\boldsymbol{y}_i^{\top}\Vert_F = \Vert \boldsymbol{x}_i\Vert \Vert\boldsymbol{y}_i\Vert来采样\boldsymbol{X},\boldsymbol{Y}对应的列/行,并乘以(r p_i)^{-1/2},就可以得到\boldsymbol{X}\boldsymbol{Y}的一个CR近似,从而将乘法复杂度从\mathcal{O}(lmn)降低到\mathcal{O}(rmn)

延伸讨论 #

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

我们可以尝试将排序的单位改为k元组。简单起见,假设k \leq l-rl-r的一个因数,l个向量\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_lk个的组合数为C_l^k,对于每个组合\{s_1,s_2,\cdots,s_k\}我们都可以算出向量和的模长\Vert \boldsymbol{v}_{s_1} + \boldsymbol{v}_{s_2} + \cdots + \boldsymbol{v}_{s_k}\Vert。有了这些数据,我们就可以贪婪地构建\eqref{eq:v-l-k-0}的近似解:

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

对于t=1,2,\cdots,(l-r)/k,执行:
     \Theta = \Theta\,\cup\,\mathop{\text{argmin}}\limits_{\{s_1,s_2,\cdots,s_k\}\subset \Omega}\Vert \boldsymbol{v}_{s_1} + \boldsymbol{v}_{s_2} + \cdots + \boldsymbol{v}_{s_k}\Vert
     \Omega = \Omega\,\backslash\,\Theta

返回\Theta

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

另一个值得思考的问题是如果允许CR近似放宽为\boldsymbol{C}\boldsymbol{\Lambda}\boldsymbol{R},那么\boldsymbol{\Lambda}的最优解是什么呢?如果不限定\boldsymbol{\Lambda}的结构,那么答案可以由伪逆给出
\begin{equation}\boldsymbol{\Lambda}^* = \mathop{\text{argmin}}_{\boldsymbol{\Lambda}}\Vert \boldsymbol{C}\boldsymbol{\Lambda}\boldsymbol{R} - \boldsymbol{X}\boldsymbol{Y}\Vert_F^2 = \boldsymbol{C}^{\dagger}\boldsymbol{X}\boldsymbol{Y}\boldsymbol{R}^{\dagger}\end{equation}
如果\boldsymbol{\Lambda}必须是对角阵呢?那可以先将问题重述为给定\{\boldsymbol{u}_1,\boldsymbol{u}_2,\cdots,\boldsymbol{u}_r\}\subset\{\boldsymbol{v}_1,\boldsymbol{v}_2,\cdots,\boldsymbol{v}_l\},求
\begin{equation}\mathop{\text{argmin}}_{\lambda_1,\lambda_2,\cdots,\lambda_r}\left\Vert\sum_{i=1}^r \lambda_i \boldsymbol{u}_i - \sum_{i=1}^l\boldsymbol{v}_i\right\Vert^2\end{equation}
我们记\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},那么优化目标可以写成
\begin{equation}\mathop{\text{argmin}}_{\boldsymbol{\lambda}}\left\Vert\boldsymbol{U}\boldsymbol{\lambda} - \boldsymbol{V}\boldsymbol{1}_{l\times 1}\right\Vert^2\end{equation}
这同样可以通过伪逆写出最优解
\begin{equation}\boldsymbol{\lambda}^* = \boldsymbol{U}^{\dagger}\boldsymbol{V}\boldsymbol{1}_{l\times 1} = (\boldsymbol{U}^{\top}\boldsymbol{U})^{-1}\boldsymbol{U}^{\top}\boldsymbol{V}\boldsymbol{1}_{l\times 1} \end{equation}
最后一个等号假设了\boldsymbol{U}^{\top}\boldsymbol{U}可逆,这通常能满足,如果不满足的话(\boldsymbol{U}^{\top}\boldsymbol{U})^{-1}(\boldsymbol{U}^{\top}\boldsymbol{U})^{\dagger}就行。

现在的问题是直接套用上式的话对原始问题来说计算量太大,因为\boldsymbol{v}_i = \text{vec}(\boldsymbol{x}_i \boldsymbol{y}_i^{\top}),即\boldsymbol{v}_imn维向量,所以\boldsymbol{V}大小为mn\times l\boldsymbol{U}大小为mn\times r,这在m,n较大时比较难受。利用\boldsymbol{v}_i = \text{vec}(\boldsymbol{x}_i \boldsymbol{y}_i^{\top})能帮我们进一步化简上式。不妨设\boldsymbol{u}_i = \text{vec}(\boldsymbol{c}_i \boldsymbol{r}_i^{\top}),那么
\begin{equation}\begin{aligned}(\boldsymbol{U}^{\top}\boldsymbol{V})_{i,j} =&\, \langle \boldsymbol{c}_i \boldsymbol{r}_i^{\top}, \boldsymbol{x}_j \boldsymbol{y}_j^{\top}\rangle_F = \text{Tr}(\boldsymbol{r}_i \boldsymbol{c}_i^{\top}\boldsymbol{x}_j \boldsymbol{y}_j^{\top}) = (\boldsymbol{c}_i^{\top}\boldsymbol{x}_j)(\boldsymbol{r}_i^{\top} \boldsymbol{y}_j) \\[5pt] =&\, [(\boldsymbol{C}^{\top}\boldsymbol{X})\otimes (\boldsymbol{R}\boldsymbol{Y}^{\top})]_{i,j} \end{aligned}\end{equation}
\boldsymbol{U}^{\top}\boldsymbol{V}=(\boldsymbol{C}^{\top}\boldsymbol{X})\otimes (\boldsymbol{R}\boldsymbol{Y}^{\top}),\boldsymbol{U}^{\top}\boldsymbol{U}=(\boldsymbol{C}^{\top}\boldsymbol{C})\otimes (\boldsymbol{R}\boldsymbol{R}^{\top}),这里的\otimesHadamard积,这样恒等变换之后\boldsymbol{U}^{\top}\boldsymbol{V}\boldsymbol{U}^{\top}\boldsymbol{U}的计算量就降低了。

文章小结 #

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

转载到请包括本文地址:https://spaces.ac.cn/archives/10427

更详细的转载事宜请参考:《科学空间FAQ》

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (Oct. 11, 2024). 《低秩近似之路(三):CR 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/10427

@online{kexuefm-10427,
        title={低秩近似之路(三):CR},
        author={苏剑林},
        year={2024},
        month={Oct},
        url={\url{https://spaces.ac.cn/archives/10427}},
}