贪婪随机扩展Kaczmarz算法研究

Research on Greedy Random Extended Kaczmarz Algorithm

摘要: 求解大规模线性方程组的问题广泛存在于各个研究领域。求解这一问题的方法有很多，其中，Kaczmarz迭代算法的研究更为广泛，该算法用于超定、欠定系统。本文分析了国内外关于Kaczmarz扩展算法与贪婪算法相结合的分块算法，并提出贪婪随机双块Kaczmarz算法(GRDBK)。GRDBK在REK算法和GRK算法的基础上进行创新，适用于超定不相容的线性系统。

Abstract: The problem of solving large-scale linear systems of equations widely exists in various research fields. There are many methods to solve this problem, among which Kaczmarz iterative algorithm is more widely studied, and it is used for overdetermined and underdetermined systems. In this paper, we analyze the algorithm combining Kaczmarz extension algorithm with greedy algorithm at home and abroad, and propose greedy random double block Kaczmarz algorithm (GRDBK). GRDBK innovates on the basis of REK algorithm and GRK algorithm, and is suitable for overdetermined and inconsistent linear system.

1. 研究背景

在计算机、医学等领域的研究中，往往需要求解多个未知线性方程。问题越复杂，方程中包含的未知数就越多。因此，一种求解线性方程组的有效方法应运而生。目前求解大规模线性方程组的主要方法有代数重建(ART)、同步图像重建(SIRT)、连续超松弛(SOR)、共轭梯度法(CG)、广义最小残差法(GMRES)和双共轭梯度稳定性法(BiCGSTAB)。其中，ART Kaczmarz方法被广泛应用，相关算法的研究也越来越多。本文讨论了一种求解大规模线性方程组的随机迭代法。下面是线性方程组的一般表达式，

$Ax=b$

其中矩阵 $A\in R^{n\times m}$ ，向量 $x\in R^m$ 是未知的，向量 $b\in R^n$ 是给定的。当线性系统是超定( $m>n$ )相容时，满列秩矩阵A对应的线性系统具有最小范数解，当线性系统超定不相容时，满列秩矩阵A对应的线性系统有最小二乘解。对于不同类型的矩阵(密集或稀疏)，方程的近似解不同，选择的方法也不同。经典Kaczmarz方法可以快速求解近似解，且不受方程个数的影响。其主要方法是将初始点正交投影到超平面上，并将迭代点作为下一次迭代的起点。通过循环系数矩阵A的每一行，最终得到无限接近于精确解的近似解。设 $x_0$ 是初始向量， ${\alpha }_k$ 是步长，值为1。根据矩阵A的每行顺序，将当前的近似向量 $x_k$ 正交投影到解的超平面上，即 $\left\{x|\langle A\left(i\right),x\rangle =b_i\right\}$ ，并将得到的向量放入第 $\left(k+1\right)$ 次迭代，具体过程为：

$x_{k+1}:=x_k+{\alpha }_k\frac{b_i-\langle A^{\left(i\right)},x_k\rangle }{{\Vert A^{\left(i\right)}\Vert }_2^2}{A}^{(\; i\; )}$

其中， ${\Vert A^{\left(i\right)}\Vert }_2^2$ 是欧氏范数， $\langle A^i,x_k\rangle$ 为欧氏内积。这个算法一经提出，就衍生出许多改进算法。2009年，Strohmer和Vershynin [1] 提出了指数收敛的随机Kaczmarz方法(RK)。该方法在经典Kaczmarz方法的基础上进行改进，改进后的算法适用于大规模线的超定性系统，并且通过理论证明，该算法以预期的指数速率收敛，RK是第一个不依赖于方程组大小的算法，它甚至不需要知道整个线性系统，只需要知道其中的随机迭代部分即可进行迭代求解。RK算法与经典Kaczmarz算法的不同之处在于，算法不再对矩阵A的所有行进行迭代，而是根据概率随机选择某些工作行，只对选中的下标集合中的行向量进行迭代操作。这样可以极大地节省计算时间，效率也更高。具体的概率运算方法如下：取矩阵A的每一行2-范数与矩阵A的F-范数之比作为选行概率，随机选取一些行，具体迭代过程如下：

$x_{k+1}:=x_k+{\alpha }_k\frac{b_{r\left(i\right)}-\langle A_{r\left(i\right)},x_k\rangle }{{\Vert A_{r\left(i\right)}\Vert }_2^2}{A}_{r(\; i\; )}$

当矩阵A是非数量矩阵时，RK算法明显优于共轭梯度法和经典Kaczmarz法。而当RK方法的系数矩阵A为数量矩阵时，每一行的范数相同，按照上述方法选择活动行的概率也相同，在这种情况下，RK方法成为经典的Kaczmarz方法。针对RK算法这一缺陷，出现了许多改进算法，比如：2015年Needell等人对矩阵A的行列进行分块 [2]、2013年Zouzias和Freris将Kaczmarz迭代过程进行扩展研究 [3]、2018年Bai和Wu将贪婪算法和Kaczmarz相结合 [4] 等。如何随机选取工作行来加快随机Kaczmarz迭代求值仍然是一个值得研究的问题。

2018年Bai和Wu [4] 提出将贪婪算法与Kaczmarz相结合，得到贪心随机Kaczmarz (GRK)算法，该算法适用于求解大规模相容稀疏线性方程组，且近似解收敛到最小范数。GRK在RK的基础上，采用贪婪方法选择工作行。贪婪算法的主要思想是过滤掉残差向量r中较大的分量，最后剩下残差小的那些分量。在这些选定的行向量中，将残差向量r的范数之比作为选行概率，如下：

$\frac{{\Vert {\hat{r}}_k\left(i_k\right)\Vert }^2}{{\Vert {\hat{r}}_k\Vert }_2^2}$ ，

选择工作行中概率较高的行或子矩阵先进入迭代。这样，优先考虑较大分量的方法就是贪婪算法。它的优点是收敛速度快，弥补了RK算法的不足。以此算法为基础，出现了许多改进算法 [5] - [11]，将GRK推广到不相容线性系统中去。

2013年Zouzias和Freris [3] 在超定不相容的线性系统中，基于RK算法提出了扩展的随机Kaczmarz，分析了在不相容的情况下，方程 $Ax=b_{R\left(A\right)}+r$ 中残差向量r的收敛情况，并给出了最小二乘解 $x_k$ 的收敛性定理。这个算法在讨论不相容线性系统的问题中，极具影响力，之后出现了很多基于REK的改进算法，例如2019年Bai和Wu [12] 提出了部分随机扩展Kaczmarz (PREK)，2020年，Du和Si [13] 提出了随机扩展双块Kaczmarz (REABK)。

本文将贪婪算法与REK算法相结合，讨论了REK算法、基于REK算法创新的其他算法和改进算法(GRDBK)在迭代次数和运行时间上的情况。通过数值试验，得出结论，对矩阵A的行列进行分块后使用改进算法进行迭代求解，大大降低了运行时长，提高收敛速度，减少误差。

本文的结构安排如下，第一节分析了算法GRDBK的研究背景，第二节具体讲述了GRDBK算法的创新点和算法步骤，分析了GRDBK的数值实例，第三节总结展望。

2. 贪婪扩展随机块Kaczmarz (GRDBK)

2.1. GRDBK算法

2013年Zouzias和Freris [3] 提出了REK算法，具体算法如下：

算法1. 随机扩展Kaczmarz.

从上面的迭代过程中可以看出，在RK算法基础上，REK算法能够适用于不相容的超定线性系统，原因是它让向量b的那部分“噪声”r趋于零，即让方程组的最小二乘解无限趋近于近似解。利用矩阵A的转置进行选列，让向量 $b_{R\left(A\right)}$ 趋近于 $b-z_k$ 。

2018年，Bai和Wu发表了Kaczmarz贪婪算法，定义了根据残差向量选择活动行下标集的方法。具体算法如下：

算法2. 随机贪婪Kaczmarz.

在上面的算法的基础上，本文做了创新改进。首先对矩阵A的行列分别分块，在扩展算法的第一部分迭代过程中，运用贪婪的思想，对分好的每个行块进行筛选，选择2-范数比值最大的行块作为迭代的活动块。在扩展算法的第二部分迭代中，我们按照对矩阵A列块的分块原始顺序，进行迭代。具体算法如下：

算法3. 贪婪随机双块Kaczmarz.

2.2. 数值实例

本小节比较了REK算法、PREK算法、REABK算法和GRDBK算法，展示了在9000 × 500的线性系统中运行的情况。

从下图1可以分析出，在9000 × 500的超定不相容线性系统中，GRDBK消耗更少的迭代时间，而REK算法和PREK算法仅次其后，REABK算法首先速度最慢，在迭代次数上，GRDBK和REABK相差无几，REK算法和PREK算法收敛步数最多，收敛的很慢。在迭代次数上，GRDBK没有明显优势。

Figure 1. 9000 × 500

图1. 9000 × 500

以下针对5000 × 500、7000 × 500、8000 × 500线性系统，分别作了数值实例。

Figure 2. 5000 × 500

图2. 5000 × 500

从图2~4可以得出结论，GRDBK在迭代时间上更快速，在迭代次数上还有进步的空间。

3. 总结与展望

本文总结了Kaczmarz与贪婪相结合的算法，并在目前研究基础上作出改进。当前有许多关于Kaczmarz在不同系统中求解线性方程的文章，其中大多数基于经典算法REK、RK、RBK [2]，添加贪婪

Figure 3. 7000 × 500

图3. 7000 × 500

Figure 4. 8000 × 500

图4. 8000 × 500

随机、松弛参数、高斯、坐标下降、变换伪逆、添加双块、部分随机、选择工作行等方法，或与其他最小二乘法相结合。贪婪算法在对矩阵下标集进行滤波时非常有利，可以达到加速收敛速度的效果。在贪婪的基础上，分块迭代可以更快。然而，求解大型线性方程的方法并不局限于Kaczmarz方法，还可以结合张量、正则化等领域进行拓展研究。事实上，将求解线性方程的问题推广到三维，并结合Kaczarz方法来解决一些不适定问题。