求解大型线性最小二乘问题的贪婪随机坐标下降法

A Greedy Randomized Coordinate Descent Method for Solving Large Linear Least Squares Problem

摘要: 贪婪随机坐标下降法(GRCD)是求解大型线性最小二乘问题的有效迭代方法之一。本文在GRCD算法中引入松弛因子，构造了一种含参数的贪婪随机坐标下降法。并证明了当线性最小二乘问题的系数矩阵为列满秩时该方法依期望的收敛性。数值实验表明，当选取适当的松弛因子时，该算法在迭代步数和计算时间比GRCD方法更有效。

Abstract: The greedy randomized coordinate descent method (GRCD) is one of the effective iterative methods to solve large linear least squares problem. A greedy randomized coordinate descent method with parameters was constructed by introducing a relaxation parameter in the GRCD algorithm. It is also proved that the method has the expected convergence when the coefficient matrix of the linear least squares problem is of full column rank. Numerical experiments show that the proposed algorithm is more effective than the GRCD method in terms of iterative steps and calculation time when the appropriate relaxation parameter is selected.

1. 引言

对于大型线性最小二乘问题

$f\left(x\right)=\underset{x\in R}{\min }{\Vert b-Ax\Vert }_2^2$ (1)

其中系数矩阵$A\in R^{m\times n}$ 列满秩，$b\in R^m$ 。$x_{\ast }=A^{+}b$ ，为其唯一最小范数最小二乘解，其中$A^{+}={\left(A^{T}A\right)}^{-1}{A}^T$ 表示$A$ 的Moore-Penrose广义逆。众所周知，线性最小二乘问题在数学上等同于正规方程$A^{T}Ax=A^{T}b$ ，也等价于无约束二次优化问题

$\underset{x\in R^n}{\min }f\left(x\right):=\frac{1}{2}{x}^T{A}^{T}Ax-b^{T}Ax$ (2)

对于矩阵$G=\left(g_{ij}\right)\in R^{m\times n}$ ，$G_{\left(j\right)}$ 、${\Vert G\Vert }_2$ 和${\Vert G\Vert }_F$ 分别代表矩阵的第j列、2范数以及F范数，$e_i$ 表示第i个坐标向量，${(·)}^{\left(j_k\right)}$ 表示第k次迭代向量的第j个元素，用${\lambda }_{\min }\left(G^{T}G\right)$ 表示$G^{T}G$ 最小的正特征值。

最小二乘问题在岭回归[1]-[3]，计算机断层扫描[4]-[6]，机器学习[7]，生物特征选择[8]及其他领域都有广泛的应用[9]-[12]。对于最小二乘问题(1)，坐标下降法是求解这一类问题的重要方法之一，它将经典的Gauss-Seidel迭代方法直接应用于求解(2) [13]。在2009年Strohmer和Vershynin [14]最先通过随机选择系数矩阵A的一行作为工作行提出随机Kaczmarz算法。2010年Leventhal和Lewis [15]受此启发，通过随机选择系数矩阵A的一列作为工作列提出随机坐标下降(RCD)方法。之后，随机坐标下降法被进一步推广或拓展[16]-[18]。2019年Bai和Wu [19]用了一种更有效的概率准则选取系数矩阵A的列，提出了贪婪随机坐标下降法(GRCD)。本文中在GRCD的算法分析中引入松弛因子，构造求解系统(1)的含参数的GRCD算法，并给出了该算法的收敛性。

2. 参数的贪婪随机坐标下降算法

Bai和Wu提出了求解大型线性最小二乘问题(1)的GRCD方法，该算法如下。

算法1 GRCD算法

输入：$A,b,l$ 和$x_0$

输出：$x_l$

1：对$k=0,1,\cdots ,l-1$ 做

2：计算

${\delta }_k=\frac{1}{2}\left(\frac{1}{{\Vert A^T{r}_k\Vert }_2^2}\right)\underset{1\le j\le n}{\max }\left(\frac{{\left|A_{\left(j\right)}^T{r}_k\right|}^2}{{\Vert A_{\left(j\right)}\Vert }_2^2}+\frac{1}{{\Vert A\Vert }_F^2}\right)$

3：构造集合

${\nu }_k=\left\{j|{\left|A_{\left(j\right)}^T{r}_k\right|}^2\ge {\delta }_k{\Vert A^T{r}_k\Vert }_2^2{\Vert A_{\left(j\right)}\Vert }_2^2\right\}$

4：令$s_k=A^T{r}_k$ ，根据如下公式计算${\tilde{s}}_k^{\left(j_k\right)}$

${\tilde{s}}_k^{\left(j_k\right)}=\{\begin{array}{l}{s}_k^{\left(j_k\right)},j_k\in {\nu }_k,\\ 0,j_k\notin {\nu }_k,\end{array}$

5：依概率$\mathrm{Pr}\left(column=j_k\right)=\frac{{\left|{\tilde{s}}_k^{\left(j_k\right)}\right|}^2}{{\Vert {\tilde{s}}_k\Vert }_2^2}$ 选择$j_k\in {\nu }_k$

6：计算$x_{k+1}=x_k+\frac{s_k^{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}$

7：结束

本文中含参数的GRCD算法(GRCD(ω))是在GRCD算法迭代公式中添加一个参数$\omega \in \left(0,2\right)$ ，见如下算法2。

算法2 含参数的贪婪随机坐标下降法

输入：$A,b,l,\omega$ 和$x_0$

输出：$x_l$

1：对$k=0,1,\cdots ,l-1$ 做

2：计算

${\delta }_k=\frac{1}{2}\left(\frac{1}{{\Vert A^T{r}_k\Vert }_2^2}\right)\underset{1\le j\le n}{\max }\left(\frac{{\left|A_{\left(j\right)}^T{r}_k\right|}^2}{{\Vert A_{\left(j\right)}\Vert }_2^2}\text{+}\frac{1}{{\Vert A\Vert }_F^2}\right)$

3：构造集合

${\nu }_k=\left\{j|{\left|A_{\left(j\right)}^T{r}_k\right|}^2\ge {\delta }_k{\Vert A^T{r}_k\Vert }_2^2{\Vert A_{\left(j\right)}\Vert }_2^2\right\}$

4：令$s_k=A^T{r}_k$ ，根据如下公式计算${\tilde{s}}_k^{\left(j_k\right)}$

${\tilde{s}}_k^{\left(j_k\right)}=\{\begin{array}{l}{s}_k^{\left(j_k\right)},j_k\in {\nu }_k,\\ 0,j_k\notin {\nu }_k,\end{array}$

5：依概率$\mathrm{Pr}\left(column=j_k\right)=\frac{{\left|{\tilde{s}}_k^{\left(j_k\right)}\right|}^2}{{\Vert {\tilde{s}}_k\Vert }_2^2}$ 选择$j_k\in {\nu }_k$

6：计算$x_{k+1}=\left(1-\omega \right)x_k+\omega \left(x_k+\frac{s_k^{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}\right)$

7：结束

注1：上述算法中依概率选取工作列时采用Matlab的randsrc命令。

注2：当$\omega =1$ 时，GRCD(ω)算法简化为GRCD算法。

3. GRCD(ω)算法收敛性分析

关于算法2有如下收敛定理。

定理1对于线性最小二乘问题(1)，设$\omega \in \left(0,2\right)$ 。那么由算法2迭代生成的序列${\left\{x_k\right\}}_{k=0}^{\infty }$ 依期望指数收敛到$x_{\ast }$ ，且解误差依期望满足

$E_k{\Vert x_{k+1}-x_{\ast }\Vert }_{A^{T}A}^2\le \left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right){\Vert x_k-x_{\ast }\Vert }_{A^{T}A}^2$ (3)

以及

$E{\Vert x_k-x_{\ast }\Vert }_{A^{T}A}^2\le {\left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)}^k{\Vert x_0-x_{\ast }\Vert }_{A^{T}A}^2$ (4)

这里用$E_k$ 表示前k次迭代的期望值，即$E_k[\cdot ]=E_k\left[\cdot \left|j_0,j_1,\cdots ,j_{k-1}\right|\right]$ ，$j_l\left(l=0,1,\cdots ,k-1\right)$ 表示算法在第l次迭代时选择的第$j_l$ 列，同时也易知$E\left[E_k[\cdot ]\right]=E[\cdot ]$ 。

证明：令$P_{j_k}=\frac{A_{\left(j_k\right)}{A}_{\left(j_k\right)}^T}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}$ ，则

$\begin{array}{c}A\left(x_{k+1}-x_{\ast }\right)=A\left(x_k-x_{\ast }\right)+\omega A\frac{s_k^{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}\\ =A\left(x_k-x_{\ast }\right)+\omega \frac{A{A}_{\left(j_k\right)}^T{r}_k}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{e}_{j_k}=A\left(x_k-x_{\ast }\right)+\omega \frac{A_{\left(j_k\right)}^T{r}_k{A}_{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\\ =A\left(x_k-x_{\ast }\right)+\omega \frac{A_{\left(j_k\right)}^T\left(A{x}_{\ast }-A{x}_k\right)A_{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\\ =A\left(x_k-x_{\ast }\right)-\omega \frac{A_{\left(j_k\right)}{A}_{\left(j_k\right)}^T}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}A\left(x_k-x_{\ast }\right)\\ =A\left(x_k-x_{\ast }\right)-\omega P_{j_k}A\left(x_k-x_{\ast }\right)\end{array}$

由$P_{j_k}$ 的定义，可知$P_{j_k}^T{P}_{j_k}=P_{j_k}$ ，因此

$\begin{array}{l}{\Vert A\left(x_{k+1}-x_{\ast }\right)\Vert }_2^2={\Vert A\left(x_k-x_{\ast }\right)-\omega P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2\\ ={\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-2\omega {\left[A\left(x_k-x_{\ast }\right)\right]}^T{P}_{j_k}A\left(x_k-x_{\ast }\right)+{\omega }^2{\left[P_{j_k}A\left(x_k-x_{\ast }\right)\right]}^T{P}_{j_k}A\left(x_k-x_{\ast }\right)\\ ={\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-2\omega {\left[A\left(x_k-x_{\ast }\right)\right]}^T{P}_{j_k}^T{P}_{j_k}A\left(x_k-x_{\ast }\right)+{\omega }^2{\left[P_{j_k}A\left(x_k-x_{\ast }\right)\right]}^T{P}_{j_k}A\left(x_k-x_{\ast }\right)\\ ={\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-2\omega {\Vert P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2+{\omega }^2{\Vert P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2\\ ={\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-\omega \left(2-\omega \right){\Vert P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2\end{array}$ (5)

又

$\begin{array}{c}{E}_k{\Vert P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2=E_k{\Vert \frac{A_{\left(j_k\right)}{A}_{\left(j_k\right)}^T}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}A\left(x_k-x_{\ast }\right)\Vert }_2^2\\ =E_k{\Vert \frac{A_{\left(j_k\right)}}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|\Vert }_2^2\\ =E_k\left({\left(\frac{\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{A}_{\left(j_k\right)}\right)}^T\left(\frac{\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{A}_{\left(j_k\right)}\right)\right)\\ =E_k\left(\frac{1}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|}^2\right)\end{array}$

$\begin{array}{c}={\displaystyle \sum _{j_k=1}^n\frac{{\left|{\tilde{s}}_k^{\left(j_k\right)}\right|}^2}{{\displaystyle \sum _{j\in {\nu }_k}{\left|{\tilde{s}}_k^{\left(j\right)}\right|}^2}}}\left(\frac{1}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|}^2\right)\\ ={\displaystyle \sum _{j_k\in {\nu }_k}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\displaystyle \sum _{j\in {\nu }_k}{\left|s_k^{\left(j\right)}\right|}^2}}}\left(\frac{1}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{\left|A_{\left(j_k\right)}^{T}A\left(x_k-x_{*}\right)\right|}^2\right)\\ ={\displaystyle \sum _{j_k\in {\nu }_k}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\displaystyle \sum _{j\in {\nu }_k}{\left|s_k^{\left(j\right)}\right|}^2}}}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\end{array}$

由${\nu }_k$ 的定义可知，若$j_k\in {\nu }_k$ ，则${\left|s_k^{\left(j_k\right)}\right|}^2\ge {\delta }_k{\Vert s_k\Vert }_2^2{\Vert A_{\left(j_k\right)}\Vert }_2^2$ ，即$\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\ge {\delta }_k{\Vert s_k\Vert }_2^2$

因此${\displaystyle \sum _{j_k\in {\nu }_k}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\displaystyle \sum _{j\in {\nu }_k}{\left|s_k^{\left(j\right)}\right|}^2}}}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}\ge {\displaystyle \sum _{j_k\in {\nu }_k}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\displaystyle \sum _{j\in {\nu }_k}{\left|s_k^{\left(j\right)}\right|}^2}}{\delta }_k{\Vert s_k\Vert }_2^2}={\delta }_k{\Vert s_k\Vert }_2^2$

即$E_k{\Vert P_{j_k}A\left(x_k-x_{*}\right)\Vert }_2^2\ge {\delta }_k{\Vert s_k\Vert }_2^2={\delta }_k{\Vert A^{T}A\left(x_k-x_{*}\right)\Vert }_2^2$

对(5)式左右两边求期望可得

$\begin{array}{l}{E}_k{\Vert A\left(x_{k+1}-x_{\ast }\right)\Vert }_2^2\\ ={\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-\omega \left(2-\omega \right)E_k{\Vert P_{j_k}A\left(x_k-x_{\ast }\right)\Vert }_2^2\\ \le {\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2-\omega \left(2-\omega \right){\delta }_k{\Vert A^{T}A\left(x_k-x_{*}\right)\Vert }_2^2\\ \le \left(1-\omega \left(2-\omega \right){\delta }_k{\lambda }_{\min }\left(A^{T}A\right)\right){\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2\end{array}$ (6)

上式最后一步是由下面的不等式得到

${\Vert A^{\text{T}}u\Vert }_2^2\ge {\lambda }_{\min }\left(A^{T}A\right){\Vert u\Vert }_2^2,u\in R\left(A^T\right)$

又由于

$\begin{array}{c}{\delta }_k{\Vert A\Vert }_F^2=\frac{1}{2}\left(\frac{1}{{\Vert A^T{r}_k\Vert }_2^2}\underset{1\le j_k\le n}{\max }\frac{{\left|A_{\left(j_k\right)}^T{r}_k\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}+\frac{1}{{\Vert A\Vert }_F^2}\right){\Vert A\Vert }_F^2\\ =\frac{\underset{1\le j_k\le n}{\max }\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}}{2{\displaystyle \sum _{j_k=1}^n\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A\Vert }_F^2}}}+\frac{1}{2}=\frac{\underset{1\le j_k\le n}{\max }\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}}{2{\displaystyle \sum _{j_k=1}^n\frac{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{{\Vert A\Vert }_F^2}\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j_k\right)}\Vert }_2^2}}}+\frac{1}{2}\\ \ge \frac{\underset{1\le j_k\le n}{\max }\frac{{\left|s_k^{\left(j_k\right)}\right|}^2}{{\Vert A_{\left(j\right)}\Vert }_2^2}}{2{\displaystyle \sum _{j_k=1}^n\frac{{\Vert A_{\left(j_k\right)}\Vert }_2^2}{{\Vert A\Vert }_F^2}\underset{1\le i\le n}{\max }\frac{{\left|s_k^{\left(i\right)}\right|}^2}{{\Vert A_{\left(i\right)}\Vert }_2^2}}}+\frac{1}{2}=\frac{1}{2}+\frac{1}{2}=1\end{array}$

因此${\delta }_k\ge \frac{1}{{\Vert A\Vert }_F^2}$ ，$k=0,1,\cdots$ ，

对(6)式进一步估计，可得到

$\begin{array}{c}{E}_k{\Vert A\left(x_{k+1}-x_{\ast }\right)\Vert }_2^2\le \left(1-\omega \left(2-\omega \right){\delta }_k{\lambda }_{\min }\left(A^{T}A\right)\right){\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2\\ \le \left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }(A^{T}A)}{{\Vert A\Vert }_F^2}\right){\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2\end{array}$

因此(3)式得证。

对上式两边求全期望可得

$E{\Vert A\left(x_{k+1}-x_{\ast }\right)\Vert }_2^2\le \left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)E{\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2$

又由上式可得

$\begin{array}{c}E{\Vert A\left(x_k-x_{\ast }\right)\Vert }_2^2\le \left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)E{\Vert A\left(x_{k-1}-x_{\ast }\right)\Vert }_2^2\\ \le {\left(1-\omega (2-\omega )\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)}^{2}E{\Vert A\left(x_{k-2}-x_{\ast }\right)\Vert }_2^2\\ \le \cdots \le {\left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)}^{k}E{\Vert A\left(x_0-x_{\ast }\right)\Vert }_2^2\end{array}$

即$E{\Vert x_k-x_{*}\Vert }_{A^{T}A}^2\le {\left(1-\omega \left(2-\omega \right)\frac{{\lambda }_{\min }\left(A^{T}A\right)}{{\Vert A\Vert }_F^2}\right)}^{k}E{\Vert x_0-x_{*}\Vert }_{A^{T}A}^2$

最后(4)式得证。

4. 数值实验

在本节中，我们通过数值实验来验证GRCD(ω)方法的有效性，所有的数值结果均是使用软件MATLAB(版本R2022b)完成的，所使用的计算机搭载Windows10操作系统和16 G内存以及3.00 GHz中央处理器(Intel I CoITM) i5-12500 CPU。

用GRCD方法和GRCD(ω)方法求解线性最小二乘系统(1)的算法算例是由Matlab函数randn服从标准正态分布n(0, 1)随机生成的矩阵，以及从文献中选取的5个具有应用背景的稀疏矩阵[20]。分别是Worldcities、divorce、cage5、Cities、abtaha1。

GRCD(ω)和GRCD方法的迭代步骤表示为IT，运行时间以秒为单位表示为CPU。表格中IT和CPU是重复运行50次所需迭代步骤数和CPU时间的中位数，为了直观地展示GRCD(ω)方法的优势，我们将GRCD(ω)方法相对于GRCD方法的加速定义为：$\text{Speed-up}=\frac{\text{CPU}\text{of}\text{GRCD}}{\text{CPU}\text{of}\text{GRCD}\left(\omega \right)}$ ，以此表示GRCD相对于GRCD(ω)在计算时间这一指标上的加速程度，其数值越大表明GRCD(ω)对同一测试矩阵加速得越快。对于选取的稀疏矩阵，表7也给出了测试矩阵的一些基本信息，例如矩阵的条件数cond(A)，稠密性$\text{density}=\frac{m\times n矩阵A中非零元素的个数}{m\times n}$ 。在数值实验中，解向量$x_{\ast }$ 是利用MATLAB函数randn(n, 1)随机生成的。当线性系统相容时，令$b=A{x}_{\ast }$ 。当线性系统不相容时取$b=A{x}_{\ast }-r_0$ ，其中$r_0$ 是属于$A^{\text{T}}$ 的零空间中的非零向量，null ($A^{\text{T}}$ )是利用MATLAB函数null生成的。在迭代运算中选取的初始向量为$x_0=0$ ，当相对误差$\text{RSE}=\frac{{\Vert x_k-x_{\ast }\Vert }_2^2}{{\Vert x_{\ast }\Vert }_2^2}\le {10}^{-6}$ 或者迭代步数超过10万步迭代停止。

对于随机生成的矩阵，我们将线性系统相容时GRCD和GRCD(ω)方法的迭代步数和计算时间列于表1~3，线性系统不相容时列于表4~6。从这些表中，线性系统是相容的还是不相容的，GRCD(ω)方法在迭代步数和CPU时间方面几乎是优于GRCD方法。我们观察到，当线性系统相容时，加速最大达到1.28；当线性系统不相容时，加速最大为1.26。

对于Davis和Hu的稀疏全秩矩阵，在表7和表8中分别列出了线性系统相容时和线性系统不相容时GRCD和GRCD(ω)方法的迭代步数和计算时间。数值实验结果表明GRCD(ω)方法在迭代步数和CPU时间方面都明显优于GRCD方法。在表7中，加速至少为2.55，最大达到5.34；在表8中，加速度至少为1.98，最大达到3.91。

Table 1. The matrix A is randomly generated, when n = 50, m is different, the numerical result of the compatible system is solved

表1. 随机生成矩阵A当n = 50，m不同时求解相容系统的数值结果

$m\times n$		1000 × 50	2000 × 50	3000 × 50	4000 × 50	5000 × 50
ω		1.04	1.03	1.01	1.01	1.01
GRCD	IT CPU	130.5 0.002578	114 0.002676	100 0.002805	98 0.003731	103 0.005030
GRCD(ω)	IT CPU	120 0.002013	108 0.002365	99 0.002564	97 0.003694	100 0.004784
Speed-up		1.28	1.13	1.09	1.03	1.05

Table 2. The matrix A is randomly generated, when n = 100, m is different, the numerical result of the compatible system is solved

表2. 随机生成矩阵A当n = 100，m不同时求解相容系统的数值结果

$m\times n$		1000 × 100	2000 × 100	3000 × 100	4000 × 100	5000 × 100
ω		1.09	1.03	1.03	1.01	1.02
GRCD	IT CPU	348.5 0.006807	246 0.006185	239 0.008076	213 0.010942	216 0.015263
GRCD(ω)	IT CPU	277.5 0.005704	236 0.005671	226 0.007313	210 0.010639	210 0.014583
Speed-up		1.20	1.09	1.10	1.03	1.05

Table 3. The matrix A is randomly generated, and when n = 150, m is different, the numerical result of the compatible system is solved

表3. 随机生成矩阵A当n = 150，m不同时求解相容系统的数值结果

$m\times n$		1000 × 150	2000 × 150	3000 × 150	4000 × 150	5000 × 150
ω		1.15	1.07	1.04	1.03	1.02
GRCD	IT CPU	612.5 0.013390	423 0.013721	371 0.018007	365.5 0.029131	338 0.032732
GRCD(ω)	IT CPU	476 0.012004	389 0.012360	354 0.015663	353 0.027459	325 0.031827
Speed-up		1.12	1.11	1.15	1.06	1.03

Table 4. The matrix A is randomly generated, and when n = 50, m is different, the numerical result of the incompatible system is solved

表4. 随机生成矩阵A当n = 50，m不同时求解不相容系统的数值结果

$m\times n$		1000 × 50	2000 × 50	3000 × 50	4000 × 50	5000 × 50
ω		1.03	1.01	1.01	1.01	1.02
GRCD	IT CPU	122 0.001985	106 0.002142	105 0.002708	99 0.003863	100 0.004740
GRCD(ω)	IT CPU	118 0.001877	105 0.002094	103 0.002662	97 0.003615	96 0.004498
Speed-up		1.06	1.02	1.02	1.07	1.05

Table 5. The matrix A is randomly generated, and when n = 100, m is different, the numerical result of the incompatible system is solved

表5. 随机生成矩阵A当n = 100，m不同时求解不相容系统的数值结果

$m\times n$		1000 × 100	2000 × 100	3000 × 100	4000 × 100	5000 × 100
ω		1.05	1.03	1.03	1.01	1.01
GRCD	IT CPU	325 0.005651	268 0.006633	242 0.008443	205 0.010361	210 0.014445
GRCD(ω)	IT CPU	283 0.005175	241 0.006070	226 0.007721	201 0.010205	207 0.013776
Speed-up		1.09	1.09	1.09	1.02	1.05

Table 6. The matrix A is randomly generated, and when n = 150, m is different, the numerical result of the incompatible system is solved

表6. 随机生成矩阵 × 当n = 150，m不同时求解不相容系统的数值结果

$m\times n$		1000 × 150	2000 × 150	3000 × 150	4000 × 150	5000 × 150
ω		1.1	1.08	1.04	1.03	1.02
GRCD	IT CPU	602 0.012388	457 0.013413	386 0.016221	358 0.024975	333 0.322180
GRCD(ω)	IT CPU	475 0.009764	400 0.011613	353 0.015425	340 0.023833	322 0.030974
Speed-up		1.26	1.16	1.05	1.05	1.04

Table 7. The GRCD(ω) and GRCD methods solve the numerical results of the compatible system in the specific example matrix

表7. GRCD(ω)和GRCD两种方法在具体算例矩阵时求解相容系统的数值结果

名称		Worldcities	divorce	cage5	Cities	abtaha1
$m\times n$		315 × 100	50 × 99	37 × 37	55 × 46	14596 × 209
density		23.87%	50.00%	17.02%	50.34%	1.68%
Cond(A)		66.00	19.39	15.42	207.15	12.23
ω		1.8	1.8	1.6	1.8	1.7
GRCD	IT CPU	3834.5 0.067084	594 0.005716	2235 0.022035	60,142 0.669865	10,342 0.554940
GRCD(ω)	IT CPU	1185.5 0.020597	162.5 0.001580	760 0.007437	11,260 0.123660	4063.5 0.554940
Speed-up		3.26	3.66	2.94	5.34	2.55

Table 8. The GRCD(ω) and GRCD methods solve the numerical results of the incompatible system in the specific example matrix

表8. GRCD(ω)和GRCD两种方法在具体算例矩阵时求解不相容系统的数值结果

名称		Worldcities	divorce	Cities	abtaha1
$m\times n$		315 × 100	50 × 99	55 × 46	14,596 × 209
density		23.87%	50.00%	50.34%	1.68%
Cond(A)		66.00	19.39	207.15	12.23
ω		1.8	1.8	1.8	1.7
GRCD	IT CPU	4568 0.077293	601 0.005446	46,131 0.495157	13816.5 0.710776
GRCD(ω)	IT CPU	1267 0.021641	173.5 0.002749	11,973 0.130001	4283 0.228498
Speed-up		3.57	1.98	3.91	3.11

5. 总结

本文在GRCD算法的基础上提出了GRCD(ω)算法，分析了算法的收敛性。数值实验结果表明，在选择适当的参数ω后，GRCD(ω)算法在迭代步数和计算时间方面明显优于GRCD算法。