Epipolar a personal journal

Cholesky 分解(三)

上文我们提到了 Cholesky 分解面临的一个问题,就是对于一些形状的矩阵存在的过度 fill-in 。我们举了两个具体的例子:一个左上箭头形的矩阵和右下箭头形的矩阵。其中,对左上箭头形的矩阵的 Cholesky 分解将得到非常稠密的结果;而对右下箭头形的矩阵的分解则保持了原始的稀疏性。

那么如果一个方程 $Ax=b$ 中的 $A$ 碰巧是左上箭头形状,我们还能使用 Cholesky 分解来加快求解吗?坏消息是:直接 Cholesky 分解的话,得到的将会是稠密的线性系统,计算量会严重增大。好消息是:我们可以通过调整变量的顺序,也就是调整矩阵的行和列,将它变成更加喜闻乐见的(右下箭头)形状。

先来介绍一下更一般的结论:假设 $P$ 是一个置换矩阵(Permutation Matrix),可以通过求解 $PAP^T Px = Pb$ 来得到方程 $Ax = b$ 的解。而通过选择合适的 $P$ ,矩阵 $M = PAP^T$ 的 Cholesky 分解可能比 $A$ 本身的 Cholesky 分解更加稀疏。

为了说明这一点,我们来看看矩阵本身的稀疏模式如何影响 Cholesky 分解的稀疏性。首先,我们忽略矩阵中具体的数值,仅将非零的元素标记出来,得到一个着色的块图,例如下图中的 A 。

根据我们在系列第一篇文章中介绍的 Cholesky 分解算法,我们总是先求出结果的第一行第一列,然后在剩余的部分中减去第一列和第一行的乘积。因此,我们先在矩阵中选择(Pivot)一个目标的行列,将它移动(Permute)到首行首列,然后将它从剩余的部分中去除(Eliminate)。对剩余的部分,我们可以如法炮制继续下去。为了简化讨论,我们还假设计算中不存在任何因为数值舍入产生零的情况。

在上图中,我们展示了两种 Pivot 策略,B1 选择了第 4 行第 4 列,而 B2 则选择了第 5 行第 5 列。随后,将对应的行和列 Permute 到行列首,得到了 C1 和 C2 。

在其后的 Eliminate 中,这两种选择便表现出了差异:对于 B1 的 Pivot 方式,在进行 Eliminate 时,并不会对剩余的部分产生任何的影响。而 B2 的 Pivot 方式导致 C2 中首行首列存在非零元,这使得在 Eliminate 时,出现了引入新的非零元的机会。D2 中,我们用金色格点高亮了这两个新出现的非零元。

B1 相比 B2 的聪明之处是比较容易看出来的:B1 所在的行和列都是尽可能稀疏的,这样就大大降低了在 Eliminate 过程中引入新的非零元的机会。

RefBLAS is Slow

很多的数值应用在核心部分都依赖快速的线性代数计算,BLAS 为这些线性代数计算提供了基本的 Building Block 。然而当前存在着非常多的 BLAS 实现,其中作为参考实现的是 netlib 的 Reference BLAS(我们简称 RefBLAS)。最近发现我们的一些项目中,为了偷懒直接使用了 RefBLAS 。这种做法其实是非常不好的。

RefBLAS 采用了 FORTRAN 语言实现,其作用是明确了 BLAS 的接口定义和行为。但是 FORTRAN 在当今已经不是主流语言了,甚至在 Windows 平台上缺少高效且免费的 FORTRAN 编译器。

RefBLAS 的实现是非常通用的,这就意味着它不会照顾到具体平台的具体特性,这其中就包括了寄存器使用、超标量加速、缓存预存取甚至专用数值指令等等可以带来巨大性能提升的特性。仅仅使用 RefBLAS 的话,等于是将优化寄托于编译器。

随着如今编译器的发展,使用纯 C 语言编写的 BLAS 早就已经赶上甚至超过了 RefBLAS 的性能。也就是说,即使不采用更加针对性的优化,RefBLAS 相较于 C 编写的 BLAS 也没有特别的优势。

进一步的,当各种针对特定计算/访问模式的优化被加入之后,整体效率的提升就更加可观了!BLAS 在各种数值程序中处在非常核心的位置,在计算中占据了相当大的一部分比例,使用 RefBLAS 造成的性能浪费将会是非常严重的。

结论就在标题中了,RefBLAS 是很慢的,绝不能因为它的代码很容易下载到,就使用它。

下面列举一些先进的 BLAS 实现:

  • Intel MKL (商业软件,且仅在 Intel CPU 上有最佳性能)
  • BLIS
  • ulmBLAS
  • OpenBLAS
  • ATLAS

Cholesky 分解(二)

Cholesky 有很多用途,其中之一是求解形如 $LL^T x = Ax = b$ 的线性方程。各类最小二乘算法中,经常会涉及到求解这类问题。

现实中,矩阵 $A$ 通常会根据问题性质具有不同的结构。这些结构有一个共同的特点,就是维度很高的同时,非零元个数很少。也就是说 $A$ 经常是稀疏矩阵。

对于稀疏矩阵的计算,我们可以采用更加高效的算法,这些算法的最原始思路便是跳过所有零元。对于 Cholesky 分解,情况也毫不例外。但是经验告诉我们:不同的稀疏模式, Cholesky 分解后的结果很有可能在稀疏性上有很大的差别。有的稀疏矩阵在分解后会得到稠密的下三角阵。这种现象被称为 fill-in 。

我们举两个例子来说明这个问题:

上面的例子里,$A$ 和 $B$ 具有相同个数的非零元,但在经过 Cholesky 分解后,一个得到的是稠密的下三角矩阵,另一个则保持了与原始矩阵相同的稀疏特性。

ORB 特征提取计时

在我的计算机上对 OpenCV 的 ORB 特征提取速度进行了测试。我首先对测试用数据集的所有图像进行 FAST 特征检测,然后对检测出的特征点进行 ORB 特征提取,测量仅包含了特征提取的时间。结果如下图:

$T(x) = ax+b$ 是对 ORB 时间开销与特征数目关系的线性拟合。 $a=1.55\times10^{-3} \mathrm{ms}$ , $b = 1.67 \mathrm{ms}$ 。

测试环境是 Intel i5-2320 3.00GHz, 8.0GB DDR3, Win 10 x64, VS 2015 U4, OpenCV 3.1.0

同样的测试对 FAST 特征检测就没什么意义了,FAST 检测需要对图像像素进行扫描,因此对图像大小更加敏感。结果的数量对 FAST 的影响比较小(见下图)。

$a = 3.56\times10^{-4} \mathrm{ms}$ ,$b = 0.527 \mathrm{ms}$ 。

从图上可以看出 FAST 检测的时间消耗基本在 5.5 ms 以内。我采用的测试图像大小均为 960x540 。 后面还需要进一步研究 ORB-SLAM 中采用的四叉树特征检测方式对 FAST 特征检测的速度影响。

Cholesky 分解(一)

当 $A$ 是对称正定阵,$L$ 是下三角矩阵,并且 $A = LL^T$,则称 $LL^T$ 是 $A$ 的 Cholesky 分解。此时,如果同时要求 $L$ 的对角元大于零,则这个分解是唯一的。下面我们从 2x2 分块矩阵推导一种 Cholesky 分解算法:

将 $A$ 沿最左和最上一列分块:

若按同样方式分块 $L$ ,即:

将 $LL^T = A$ 计算出来,便有:

于是可知,$l_{11} = \sqrt{a_{11}}$,$L_{21} = \frac{1}{l_{11}}A_{21}$ 。并且 $L_{22}$ 是矩阵 $A_{22}-L_{21}L_{21}^T$ 的 Cholesky 分解。

将这一过程从左上角不断向右下角进行,便可得到 $A$ 的 Cholesky 分解。