低秩近似之路（一）：伪逆

Author: [苏剑林]

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

最佳排版请看原博客：

低秩近似之路（一）：伪逆 - 科学空间|Scientific Spaces可能很多读者跟笔者一样，对矩阵的低秩近似有种熟悉而又陌生的感觉。熟悉是因为，低秩近似的概念和意义都不难理解，加之目前诸如LoRA等基于低秩近似的微调技术遍地开花，让低秩近似的概念在耳濡目染间就已经深入人心；然而，低秩近似所覆盖的内容非常广，在低秩近似相关的论文中时常能看到一些不熟悉但又让我们叹为观止的新技巧，这就导致了一种似懂非懂的陌生感。

因此，在这个系列文章中，笔者将试图系统梳理一下矩阵低秩近似相关的理论内容，以补全对低秩近似的了解。而在第一篇文章中，我们主要介绍低秩近似系列中相对简单的一个概念——伪逆。

优化视角

伪逆（Pseudo Inverse），也称“广义逆（Generalized Inverse）”，顾名思义就是“广义的逆矩阵”，它实际上是“逆矩阵”的概念对于不可逆矩阵的推广。

我们知道，对于矩阵方程$AB=M$，如果$A$ 方阵且可逆，那么直接得到$B=A^{-1}M$，可如果$A$ 可逆或者干脆不是方阵呢？这种情况下我们很可能找不到$B$ 足$AB=M$，这时候我们如果还想继续下去的话，通常是转为优化问题

$\mathop{\text{argmin} }_B \Vert AB - M\Vert_F^2\label{eq:loss-ab-m} \$

其中$A\in\mathbb{R}^{n\times r}, B\in\mathbb{R}^{r\times m}, M\in\mathbb{R}^{n\times m}$，注意数域是$\mathbb{R}$，表明本系列文章关注的都是实矩阵，而$\Vert\cdot\Vert_F$ 矩阵的$F$ 数（Frobenius Norm），用来衡量矩阵$AB - M$ 全零矩阵的距离，其定义为

$\Vert M\Vert_F = \sqrt{\sum_{i=1}^n\sum_{j=1}^m M_{i,j}^2} \$

说白了，就是从求精确的逆矩阵改为最小化$AB$ $M$ 平方误差。本系列的主题是低秩近似，所以下面假设$r \ll \min(n,m)$，其机器学习意义就是通过低维的、有损的输入矩阵$A$ 线性变换$B$ 重建完整的$M$ 阵。

当$m=n$ $M$ 单位阵$I_n$ ，我们就得到一个只依赖于$A$ 结果，记为

$A^{\dagger} = \mathop{\text{argmin} }_B \Vert AB - I_n\Vert_F^2\label{eq:loss-ab-m-b} \$

它的作用类似于$A$ 逆矩阵，所以称为$A$ “(右)伪逆”。类似地，如果给定的是$B$ 阵，那么我们也可以将优化参数改为$A$，得到$B$ “(左)伪逆”：

$B^{\dagger} = \mathop{\text{argmin} }_A \Vert AB - I_n\Vert_F^2 \$

范数相关

在进一步推导之前，我们先补充一下$F$ 数的相关介绍。向量的范数想必很多读者已经熟知，比较经典的就是“$p$-范数”：对于$x=(x_1,x_2,\cdots,x_m)$，其$p$-范数定义为

$\Vert x\Vert_p = \sqrt[\uproot{10}p]{\sum_{i=1}^m |x_i|^p} \$

而$p$-范数中，最常见的就是$p=2$ 情形，它就是我们经常说的向量模长，也叫“欧几里得范数”，如果省略范数的下标只写$\Vert x\Vert$，那么基本上就是默认$p=2$。

矩阵的范数稍微复杂一些，它至少有两种不同但都常用的范数，其中一种就是上一节已经提到的$F$ 数，它是直接将矩阵展平为向量来计算的范数：

$\Vert M\Vert_F = \Vert \text{vec}(M)\Vert_2 = \sqrt{\sum_{i=1}^n\sum_{j=1}^m M_{i,j}^2} \$

其他矩阵范数我们遇到时再作介绍。由于矩阵范数的多样性，所以$\Vert M\Vert_F$ 下标${}_F$ 常不能省略，以免引起混淆。$F$ 数是将矩阵当成向量然后照搬向量范数的定义而来的，由此启发我们可以尝试把更多的向量运算搬到矩阵中去，比如内积：

$\langle P, Q \rangle_F = \langle \text{vec}(P), \text{vec}(Q) \rangle = \sqrt{\sum_{i=1}^n\sum_{j=1}^m P_{i,j} Q_{i,j} } \$

这称为矩阵$P,Q$ $F$ 积（Frobenius Inner Product），其中$P,Q\in\mathbb{R}^{n\times m}$，它可以用向量的迹运算来表示

$\langle P, Q \rangle_F = \text{Tr}(P^{\top} Q) \$

这可以直接由矩阵乘法和迹的定义来证明（请读者尝试一下）。当$P,Q$ 由多个矩阵连乘而来时，转换为等价的迹运算通常能帮助我们进行化简。比如，利用它我们可以证明正交变换不改变$F$ 数：假设$U$ 一个正交矩阵，利用$\Vert M\Vert_F^2 = \langle M, M \rangle_F$ 及$F$ 积与迹的关系，得到

$\Vert UM\Vert_F^2 = \text{Tr}((UM)^{\top} UM)= \text{Tr}(M^{\top} U^{\top}UM)=\text{Tr}(M^{\top} M) = \Vert M\Vert_F^2 \$

矩阵求导

言归正传，对于一个优化目标，最理想的结果自然是能够通过求导来求出解析解，而$\eqref{eq:loss-ab-m}$ 好能实现这一点！改结论可以简单“目测”出来：$AB-M$ 关于$B$ 线性函数，所以$\Vert AB-M\Vert_F^2$ 关于$A$ $B$ 二次函数，二次函数的最小值有解析解的。

要求$\mathcal{L}=\Vert AB-M\Vert_F^2$ 于$B$ 导数，首先要求$\mathcal{L}$ 于$E=AB-M$ 导数，然后求$E$ 于$B$ 导数，最后通过链式法则组合起来，即

$\frac{\partial \mathcal{L} }{\partial B_{i,j} } = \sum_{k,l}\frac{\partial \mathcal{L} }{\partial E_{k,l} }\frac{\partial E_{k,l} }{\partial B_{i,j} } \$

根据定义$\mathcal{L}=\Vert E\Vert_F^2 = \sum_{i,j} E_{i,j}^2$，很明显在求和的众多平方中只有当$(i,j)=(k,l)$ ，关于$E_{k,l}$ 导数才不为零，所以$\mathcal{L}$ 于$E_{k,l}$ 导数就是$E_{k,l}^2$ 于$E_{k,l}$ 导数，即$\frac{\partial \mathcal{L} }{\partial E_{k,l} }=2E_{k,l}$；接着，根据矩阵乘法的定义有

$E_{k,l} = \left(\sum_{\alpha} A_{k,\alpha}B_{\alpha,l}\right) - M_{k,l} \$

类似地，只有当$(\alpha,l)=(i,j)$ ，上式对$B_{i,j}$ 导数才会产生非零的结果$A_{k,i}$，所以我们可以写出$\frac{\partial E_{k,l} }{\partial B_{i,j} }=A_{k,i}\delta_{l,j}$，这里的$\delta_{l,j}$ Kronecker符号，用来声明$l=j$ 条件。将结果组合起来，我们就得到

$\frac{\partial \mathcal{L} }{\partial B_{i,j} } = \sum_{k,l}E_{k,l}A_{k,i}\delta_{l,j} = \sum_k E_{k,j}A_{k,i} \$

如果我们约定，标量对矩阵的梯度形状跟矩阵本身一致，那么可以写出

$\frac{\partial \mathcal{L} }{\partial B} = 2A^{\top}(AB-M) \$

虽然推导过程破费周折，但好在结果还是很直观的：直觉上$\frac{\partial \mathcal{L} }{\partial B}$ 是$2(AB-M)$ $A$ 乘积（类比标量求导），而我们已经约定$\frac{\partial \mathcal{L} }{\partial B}$ 状跟$B$ 状一致（即$r\times m$），所以就要想办法通过$2(AB-M)\in\mathbb{R}^{n\times m}$ $A\in\mathbb{R}^{n\times r}$ 乘来凑出一个$r\times m$ 结果来，结果就只有上式右端一种方式。根据同样原理，我们可以快速写出

$\frac{\partial \mathcal{L} }{\partial A} = 2(AB-M)B^{\top} \$

基本结果

现在我们已经分别求出了$\frac{\partial \mathcal{L} }{\partial B}$ $\frac{\partial \mathcal{L} }{\partial A}$，让它们等于零便可以解出相应的最优解 $\begin{align} 2A^{\top}(AB-M)=0\quad\Rightarrow\quad (A^{\top} A)^{-1}A^{\top}M = \mathop{\text{argmin} }_B \Vert AB - M\Vert_F^2 \ 2(AB-M)B^{\top}=0\quad\Rightarrow\quad MB^{\top}(B B^{\top})^{-1} = \mathop{\text{argmin} }_A \Vert AB - M\Vert_F^2 \end{align}\$ 代入$M=I_n$，就得到

$\begin{aligned}A^{\dagger} = (A^{\top} A)^{-1}A^{\top} \label{eq:p-inv-a}\ B^{\dagger} = B^{\top}(B B^{\top})^{-1}\label{eq:p-inv-b}\end{aligned}\$

如果$A$ $B$ 可逆方阵，那么容易证明伪逆就等于逆矩阵，即$A^{\dagger}=A^{-1},B^{\dagger}=B^{-1}$。此外，根据上式我们还可以验证：

1、$(A^{\dagger})^{\dagger}=A,(B^{\dagger})^{\dagger}=B$，即伪逆的伪逆等于自身，这意味着伪逆在作为近似逆矩阵的同时，保全了自身的信息；
2、$AA^{\dagger}A=A,B^{\dagger}BB^{\dagger}=B^{\dagger}$，即$AA^{\dagger},B^{\dagger}B$ 然没法成为单位阵$I$，但对$A,B^{\dagger}$ 说它们起到了单位阵的作用。

顺便说一下，矩阵的伪逆实际上是一个很宽泛的概念，它有很多种不同的形式，这里我们介绍的实际上是最常见的“https://en.wikipedia.org/wiki/Moore–Penrose_inverse”，除此之外还有“Drazin逆”、“Bott–Duffin逆”等，但这些笔者也不了解，所以就不作展开，读者可以自行参考维基百科的“https://en.wikipedia.org/wiki/Generalized_inverse”条目。

一般形式

不过，事情还没完。式$\eqref{eq:p-inv-a},\eqref{eq:p-inv-b}$ 立还有个关键前提是相应的$A^{\top} A$ $B B^{\top}$ 逆，如果不可逆呢？

我们以$A^{\dagger}$ 例，假设$A^{\top} A$ 可逆，那么意味着$A$ 秩不足$r$，我们只能从中找到$s < r$ 列向量构成的极大线性无关组构成矩阵$A_s\in\mathbb{R}^{n\times s}$，然后$A$ 以表示成$A_s P$，其中$P\in\mathbb{R}^{s\times r}$ $A_s$ $A$ 矩阵。此时

$\mathop{\text{argmin} }_B \Vert AB - I_n\Vert_F^2 = \mathop{\text{argmin} }_B \Vert A_s PB - I_n\Vert_F^2 \$

如果$B$ 最优解仍然记为$A^{\dagger}$，那么我们只能确定

$PA^{\dagger} = A_s^{\dagger} = (A_s^{\top} A_s)^{-1}A_s^{\top} \$

由于已经假设了$A_s$ 极大线性无关组，所以$A_s^{\top} A_s$ 然可逆，因此上式是良好定义的。然而，从$A^{\dagger}$ $PA^{\dagger}$ 一个降维到过程，这意味着存在多个$A^{\dagger}$ 得$PA^{\dagger} = A_s^{\dagger}$，即此时目标$\eqref{eq:loss-ab-m-b}$ 最优解并不唯一，换言之当$A^{\top} A$ 可逆时我们无法只凭目标$\eqref{eq:loss-ab-m-b}$ 确定唯一的伪逆$A^{\dagger}$。

一个可能的思路是补充$(A^{\dagger})^{\dagger}=A$ $A^{\dagger}AA^{\dagger}=A^{\dagger}$，这样结合$PA^{\dagger} = A_s^{\dagger}$ 可以唯一地确定$A^{\dagger}$。然而，这样打补丁的味道太浓了，实际上我们可以用一个精妙的技巧更优雅和统一地处理这个问题。问题出在当$A^{\top} A$ 可逆时，目标函数$\eqref{eq:loss-ab-m-b}$ 是严格正定的，我们可以加上一个正则项让它变得正定，求出结果后再让正则项的权重趋于零：

$A^{\dagger} = \lim_{\epsilon\to 0},\mathop{\text{argmin} }_B \Vert AB - I_n\Vert_F^2 + \epsilon\Vert B\Vert_F^2 \$

这里$\epsilon > 0$，$\epsilon\to 0$ 指正向趋于零。由上式可以解得

$A^{\dagger} = \lim_{\epsilon\to 0},(A^{\top} A + \epsilon I_r)^{-1}A^{\top}\label{eq:a-pinv-lim} \$

当$\epsilon > 0$ ，$A^{\top} A + \epsilon I_r$ 然是可逆的（请读者证明一下），因此上式是良好定义的。由于$\epsilon\to 0$ ，正则项可以忽略不计，所以上述极限必然是存在的。注意，这里说的是整体极限的存在性，当$A^{\top} A$ 可逆时，极限$\lim\limits_{\epsilon\to 0},(A^{\top} A + \epsilon I_r)^{-1}$ 不存在的（结果会出现无穷大），只有乘上$A^{\top}$ 整体再取极限才是正常的。

式$\eqref{eq:a-pinv-lim}$ 为伪逆的一般推广有什么优点呢？首先，我们已经有了$A^{\top} A$ 逆时的$A^{\dagger}$ 达式$\eqref{eq:p-inv-a}$，式$\eqref{eq:a-pinv-lim}$ 为它的推广，具有直观且形式一致的理论优雅性；其次，形式上的一致也使得$A^{\top} A$ 逆时$A^{\dagger}$ 性质如$(A^{\dagger})^{\dagger}$ 够得以保留，从而让我们在讨论$A^{\dagger}$ 几乎可以完全不考虑$A^{\top} A$ 可逆性。

数值计算

当然，目前的式$\eqref{eq:a-pinv-lim}$ 是一个形式化的定义，如果直接利用它来数值计算的话，就必须取一个足够小的$\epsilon$ 后把$(A^{\top} A + \epsilon I_r)^{-1}$ 出来，这样必然会面临严重的数值不稳定性。为了得到一个稳定的计算方式，我们利用实对称矩阵总可以正交对角化这一特点（https://en.wikipedia.org/wiki/Spectral_theorem），对$A^{\top} A$ 如下分解：

$A^{\top} A = U\Lambda U^{\top} \$

其中$U$ 正交矩阵，$\Lambda=\text{diag}(\lambda_1,\lambda_2,\cdots,\lambda_r)$ 特征值组成的对角矩阵，由于$A^{\top} A$ 半正定性，它的特征值总是非负的。利用这个分解，我们有

$\begin{aligned} (A^{\top} A + \epsilon I_r)^{-1}A^{\top} =&, (U\Lambda U^{\top} + \epsilon I_r)^{-1} A^{\top} \ =&, [U(\Lambda + \epsilon I_r) U^{\top}]^{-1}A^{\top} \ =&, U(\Lambda + \epsilon I_r)^{-1} U^{\top}A^{\top} \end{aligned}\$

对于$(\Lambda + \epsilon I_r)^{-1}$ 们有

$(\Lambda + \epsilon I_r)^{-1} = \text{diag}\Big((\lambda_1 + \epsilon)^{-1},(\lambda_2 + \epsilon)^{-1},\cdots,(\lambda_r + \epsilon)^{-1}\Big) \$

如果$\lambda_i > 0$，那么$\lim\limits_{\epsilon\to 0},(\lambda_i + \epsilon)^{-1}=\lambda_i^{-1}$，这是有限的结果，不妨碍计算，问题出现在$\lambda_i = 0$ $\lim\limits_{\epsilon\to 0},(\lambda_i + \epsilon)^{-1}=\lim\limits_{\epsilon\to 0},\epsilon^{-1}\to\infty$ 。然而，我们知道$\epsilon\to 0$ 正则项的影响就会消失，所以断定极限$\eqref{eq:a-pinv-lim}$ 然不会出现无穷大值，因此如果存在$\lambda_i=0$，那么右端所乘的$U^{\top}A^{\top}$ 然有办法抵消$\lim\limits_{\epsilon\to 0}, \epsilon^{-1}$ 来的无穷大。而能抵消这种无穷大的，唯有“乘以0”，即$\lim\limits_{\epsilon\to 0}, \epsilon^{-1}\times 0 = 0$。

换句话说，如果$\lambda_i=0$，那么$U^{\top}A^{\top}$ $(\lambda_i+\epsilon)^{-1}$ 乘的因子必然是$0$。既然如此，由于“0乘任何数都得0”，所以其实$\lambda_i=0$ $(\lambda_i+\epsilon)^{-1}$ 取值反而不重要，我们可以简单地让它等于0。这样一来，我们就得到了一种通用的计算$A^{\dagger}$ 简单方法：

$A^{\dagger} = U\Lambda^{\dagger}U^{\top}A^{\top}, \quad A^{\top} A = U\Lambda U^{\top} \$

其中$\Lambda^{\dagger}$ 示对角线上的元素如果等于零则不变，非零则取倒数。

可能有读者疑问，既然“0乘任何数都得0”，那么为什么等于零的$\lambda_i$ 不变呢？随意取一个别的值可以吗？其实这里随便取一个值也不会影响结果的，但由于我们用了$\Lambda^{\dagger}$ 个记号，那么就要保持它跟式$\eqref{eq:a-pinv-lim}$ 一致性，即它跟将对角阵$\Lambda$ 入式$\eqref{eq:a-pinv-lim}$ 直接计算结果要一致

$\Lambda^{\dagger} = \lim_{\epsilon\to 0},(\Lambda^{\top} \Lambda + \epsilon I_r)^{-1}\Lambda^{\top} = \text{diag}(\kappa_1,\kappa_2,\cdots,\kappa_n),\quad \kappa_i = \left{\begin{aligned}\lambda_i^{-1}, ,,\lambda_i\neq 0 \ 0, ,,\lambda_i=0 \end{aligned} \right. \$