前三篇文章我们分别介绍了伪逆、SVD、CR 近似,它们都可以视为寻找特定结构的低秩近似:其中 。当不再添加其他约束时,最优解由 SVD 给出;当约定 并且 之一已经给出求另一半的最优解时,最优解可以通过伪逆来给出;如果约定且,那么就是 CR 近似关心的问题。CR 近似通过选择原本矩阵的行/列来构建低秩近似,这使得近似结果更具解释性,同时也适用于一些非线性场景,但 CR 近似的前提是矩阵 本身是由两个矩阵相乘而来,它的初衷是降低矩阵计算量,而对于直接给出矩阵 的场景,类似的低秩近似则由 ID 给出。具体来说,在 ID 中有,其中是 的若干列, 是任意的,即以 若干列为骨架来逼近它自己:根据《低秩近似之路(一):伪逆》的结果,如果 已经确定,那么 的最优解就是 ,所以 ID 的实际难度只有 S 的优化,即列的选取,这是一个组合优化问题,精确求解是 NP-Hard 的,所以主要目前是寻找效率和精度都适当的近似算法。
几何意义
在试图求解之前,我们先来进一步了解一下 ID 的几何意义,这有助于我们更好地理解它的应用场景和求解思路。我们先将 表示为列向量形式 ,那么对于任意列向量 ,我们有所以 的几何意义是 的列向量的线性组合。注意 选自 ,所以 ID 就是说选择若干列作为(近似的)基向量,将剩下的列都表示为这些基的线性组合,这就是 ID 中的 “I(Interpolative,插值)”的含义。我们知道,“Interpolative” 更准确的含义是“内插”,为了更好地突出“内插”这一特性,有些文献会给 ID 的定义加上 的条件( 是矩阵 的任意元素)。当然这个条件实际上也比较苛刻,保证它严格成立的难度可能也是 NP-Hard 的,所以很多文献会将其放宽为 ,大多数近似算法的实际表现都能让这个界成立,如果没有其他需求,只考虑逼近误差的最优,那么也可以去掉这个限制。 QR分解
ID 的求解算法分为确定性算法和随机算法两大类,其中确定性算法的计算量更大但近似程度往往更优,反之随机算法计算效率更高但精度稍次。注意它们都只是实际表现尚可的近似算法,并且都不排除有完全失效的极端例子的可能性。
第一个被视为标准的近似算法是基于 QR 分解的,更准确地说是 Column-Pivoting 的 QR 分解,中文常译作“列主元 QR 分解”(不过笔者感觉还不如干脆意译为“列驱 QR 分解”),它是一个确定性算法。为什么 ID 会跟 QR 分解联系在一起呢?我们可以从 的求法出发来理解。前面提到,如果 已经给定,那么 的最优解就是 ,这个答案当然是对的,但不够直观。不失一般性,假设 线性无关,那么从几何的角度看,求 形式的最优近似实际上就是将 的每个列向量都投影到 构成的 维子空间中,而为了求出这个投影结果,我们可以先将 执行 Gram-Schmidt 正交化,使其变为一组标准正交基,然后在标准正交基上投影就简单多了,而正交化的过程,自然就对应 QR 分解了。Gram-Schmidt 正交化是递归执行如下步骤:其结果将 表示成:有了 ,那么矩阵 的第 列 在 上的最优逼近和误差就分别是 列驱 QR当然,上述结果是在已知 的前提下得到的,那怎么从 中挑出比较优的 列构成 呢?列驱 QR 分解给出了一个参考答案。一般来说,如果我们要对 做 Gram-Schmidt 正交化的话,都是按照顺序来的,即从 出发,接下来是 ,而列驱 QR 分解则根据模长来修改了正交化顺序,写成公式是说白了,列驱 QR 分解就是每一步都选择剩下的误差最大的列来执行正交归一化。除了执行顺序有所变化外,列驱 QR 分解跟普通 QR 分解的计算并无其他不同之处,所以列驱 QR 分解的最终形式可以表示为其中 是列置换矩阵。根据每一步都选择误差(模长)最大列的操作,我们可以得到对于任意 ,子矩阵 的第一列模长是最大的,它不小于剩余任意一列的模长,即由此还可以推得。这些性质让我们相信,如果想要 的一个 秩近似,只保留 的前 行应该是一个不错的选择,即注意我们之前约定过切片的优先级高于求逆,所以这里的含义是 。不难发现实际上就是矩阵 的 列,所以上式实际上给出了一个 ID 近似:以上就是基于列驱 QR 分解的 ID 求解算法,也是 SciPy 内置的求解算法(rand=False)。注意该算法是无法保证 或者 的,但根据很多文献的反馈,在实践中它几乎不会出现 ,所以这算是一个比较良好的求解算法。此外,SciPy 也内置了列驱 QR 分解,在 scipy.linalg.qr 中设置 pivoting=True 即可打开。 随机求解列驱 QR 分解每一步正交化操作,需要遍历剩余的所有向量取误差最大者,这在 很大时往往难以接受,另一方面如果 很大,那么模长、内积的计算量也会很高,于是随机算法便应运而生,它设法减少 或 的值来降低计算复杂度。首先我们来看降低 n 的思路,即降低 的每个列向量的维度,常用的方法是随机投影,跟《让人惊叹的 Johnson-Lindenstrauss 引理:理论篇》介绍的 “JL 引理”如出一辙。具体来说,假设 是某个随机投影矩阵(其中 ),它的元素是从某个分布如 独立重复采样出来的,那么我们考虑在小矩阵 上执行列驱 QR 分解来确定被选中的 列的位置。更详细的介绍可以参考《Randomized algorithms for pivoting and for computing interpolatory and CUR factorizations》[3]。根据笔者有限的调研显示,SciPy 求解 ID 的随机算法也是用的是类似思路,只是把随机采样的矩阵换成了更加结构化的 “Subsampled Randomized Fourier Transform(SRFT)”,使得 这一步的计算量可以从 降到 。不过 SRFT 以及 SciPy 的实现细节笔者也不了解,有兴趣的读者可以参考《Enabling very large-scale matrix computations via randomization》[4]、《A brief introduction to Randomized Linear Algebra》[5] 等资料进一步深究。没有深究 SRFT 等复杂随机投影方法的另一个原因,是论文《Efficient Algorithms for Constructing an Interpolative Decomposition》[6] 发现更简单的列采样往往能得到更优的结果,而且还特别好理解,就是从 中随机采样 列,然后用列驱 QR 分解从这 列中选出 列作为 ,最后再来根据 求解 ,这样就将列驱 QR 分解的矩阵大小从 降低到 。实验显示这样的简单思路顶多在个别任务上增大了 的风险,但在误差上有明显优势:
从上述表格可以留意到一个也许会让人觉得意外的结果:随机列采样的 Optim-RID,在误差方面不仅优于同样是随机算法的 SciPy-RID,在个别任务上甚至还优于确定性算法 SciPy-ID 和 Optim-ID(它们数学上是等价的,都是基于完整的列驱 QR 分解,只是实现上的效率有所不同)。
这个看似反直觉的现象,实则说明了一个事实:列驱 QR 分解虽然可以作为 ID 的一个较好的 baseline,但它选择基的能力可能跟随机选择差不了多少,列驱 QR 分解的主要作用只是保证大概率成立 。其实这也不难理解,我们以 为例,这时候列驱 QR 分解就是返回模长最大的那一列,可是模长最大的那列一定是好的(能使得重构误差最小的)基底吗?显然不是,好的基底应该是多数向量共同指向的方向,模长最大不能体现这一点。对于 ID 来说,列驱 QR 分解本质上是一种贪心算法,它将选 列贪心地转化为多个选 1 列的递归,而当 时,在 不算太大或者要求高精度的场景,通过枚举来精确求解的复杂度是可以接受的即遍历所有的 ,将剩下的所有列都投影到 上来计算总误差,选择总误差最小的 ,其复杂度跟 成正比。如果将列驱 QR 分解每一步选择模长最大的操作改为选择总误差最小的上式,那么就能找出更好的基底,从而实现更低的重构误差(代价自然是复杂度更高了,而且更加无法保证 了)。
总的来说,由于精确求解的 NP-Hard 性,所以 ID 有非常多的求解思路,上面列举的只是非常有限的几种,有兴趣的读者可以可以围绕着 Randomized Linear Algebra、Column Subset Selection 等关键词深入搜索。
特别要指出的是,Randomized Linear Algebra,旨在通过随机方法来加速各种矩阵运算,本身已经成为了一个内容丰富的学科,本文的随机 ID 和上一篇基于采样的 CR 近似,都是这个学科的经典例子。 文章小结本文介绍了 ID(Interpolative Decomposition,插值分解),它通过从原矩阵中选择若干列来作为“骨架”来逼近原矩阵,是一种具有特定结构的低秩分解,几何意义相对来说更加直观,其核心难度是列的选择,本质上是一个 NP-Hard 的离散优化问题。