求解大型线性系统的K-均值最大残差块方法

K-Means Maximum Residual Block Method for Solving Large Linear Systems

摘要: 求解大型线性系统，带K-均值聚类的贪婪随机块Kaczmarz方法是近几年被广受关注的一类方法。本文在该方法基础上做了进一步的研究即在每一次迭代中优先消除残差向量中的最大块，构建了最大残差块Kaczmarz方法及其加速版本并进行了收敛性分析。数值实验证实了本文算法的有效性。

Abstract: The greedy random block Kaczmarz method with K-means clustering for solving large linear systems has been widely studied in recent years. This article conducted further research on this method by prioritizing the elimination of the largest block in the residual vector in each iteration, constructing the Kaczmarz method for the maximum residual block and its accelerated version, and conducting convergence analysis. Numerical experiments have confirmed the effectiveness of the algorithm proposed in this paper.

1. 引言

随着数据科学的发展和现实世界问题规模的不断扩大，研究求解大型线性系统的有效算法成为一个重要的课题。具体地说，我们考虑如下形式的线性系统

$Ax=b$ (1)

其中$A\in R^{m\times n}$ 是一个$m\times n$ 矩阵，是已知向量且系统(1)是相容的。研究者通常使用选代方法求解此类问题。其中一种方法是Kaczmarz方法[1]，该方法的选代形式如下：

$x_{k+1}=x_k+\frac{b^{\left(i_k\right)}-A^{\left(i_k\right)}{x}_k}{{\Vert A^{\left(i_k\right)}\Vert }_2^2}{\left(A^{\left(i_k\right)}\right)}^{\text{T}}$

其中$A^{\left(i\right)}$ 表示A的第i行，$b^{\left(i\right)}$ 表示b的第i个元素。

对于任何向量$y\in {ℂ}^n$ ，${\left(y\right)}^{\left(j_k\right)}$ 表示第k次选代的第j个元素，用$e_i$ 表示除了第i个位置的元素为1其余为零的列向量，给定任意矩阵B，对于矩阵$B=\left(b_{ij}\right)\in {ℝ}^{m\times n}$ ，$B_{\left(j\right)}$ 表示B的第j列，如果B对称的，用${\sigma }_{\min }\left(B\right)$ 表示B的最小奇异值，用${\sigma }_{\max }\left(B\right)$ 表示B的最大奇异值。此外，我们定义[m]为$1,2,\cdots ,m$ ，其中m是任意的正整数。我们考虑集合$\mathcal{V}=\left\{{\mathcal{V}}_1,{\mathcal{V}}_2,\cdots ,{\mathcal{V}}_k\right\}$ 作为[m]的一个划分，如果索引集$\mathcal{V}$ ，其中$i=1,2,\cdots ,k$ ，满足条件${\mathcal{V}}_i\cap {\mathcal{V}}_j=\varnothing$ ，其中$\varnothing$ 是空集。对于$i\ne j$ ，以及$U_{i=1}^k{V}_i=\left[m\right]$ 。进一步地，给定一个行索引集$\mathcal{V}$ ，我们使用$A_{{\mathcal{V}}_i}$ 表示由$\mathcal{V}$ 索引的矩阵A的子矩阵，使用$b{\nu }_i$ 表示向量b的子向量，$A_{{\mathcal{V}}_i}$ 和$b_{{\mathcal{V}}_i}$ 的大小由K-means聚类后的结果决定。

为了便于解释，我们简要说明一下K-means方法。K-means聚类是一种迭代算法，他将输入的数据X分为k组，每个项分配给具有最近质心的聚类。这些组$C=\left\{C_1,\cdots ,C_k\right\}$ 满足

${\displaystyle \underset{v=1}{\overset{k}{\cup }}{C}_v}=X,C_v\cap C_j=\varnothing$

其中$1\le v\ne j\le k$ ，$\varnothing$ 是一个空集。

在K-means均值聚类中，聚类的数量由K决定，是我们指定的正整数。最初，k个聚类是从给定数据中随机选择的然后将一个项分配给质心最近的聚类所谓“最接近”就是指对它们的距离或者相似度最小。对于距离我们定义为：

${\displaystyle \sum _{r=1}^k{\displaystyle \sum _{x_i\in C_r}{\Vert x_i-{\rho }_r\Vert }_2^2}}$

其中$X=\left\{x_i\right\}$ 是给定的数据，是聚类的质心，$i=1,2,\cdots ,m$ ；$v=1,2,\cdots ,k$ 。对于相似系数，我们使用以下标准来衡量项目对的接近度，该标准为

${\displaystyle \sum _{v=1}^k{\displaystyle \sum _{x_1\in C_v}\left|1-\frac{x_i^{\text{T}}{\rho }_v}{{\Vert x_i\Vert }_2{\Vert {\rho }_v\Vert }_2}\right|}}$

上述过程相继进行，直到不发生分配结束。Jiang和Zhanng在[2]中给出了完整过程。

一个列划分函数：${\mathcal{V}}_1,{\mathcal{V}}_2,\cdots ,{\mathcal{V}}_k$ 是由算法开始的初始步所定得到。$X=\left[A,b\right]\in C^{m\times \left(n+1\right)}$ 通过K-means方法我们得到关于A和b的行指标$\left\{1,\cdots ,m\right\}$ 的划分$T=\left\{v_1,\cdots ,v_k\right\}$ ，也就是：$A=\left[A_{v_1};\cdots ;A_{v_k}\right]$ 和$b=\left[b_{v_1};\cdots ;b_{v_k}\right]$ ；${\overline{A}}_{v_i}$ 和${\overline{b}}_{v_i}$ 分别是块$A_{v_i}$ 和$b_{v_i}$ 相对应的簇中心，$i=1,2,\cdots ,k$ 。我们可以写成紧凑形式即：$\overline{A}=\left[\begin{array}{c}{\overline{A}}_{v_1};\cdots ;A_{v_k}\end{array}\right]$ 和$\overline{b}=\left[\begin{array}{c}{\overline{b}}_{v_1};\cdots ;{\overline{b}}_{v_k}\end{array}\right]$ 。

求解线性和非线性问题在各个领域的应用都比较广泛，例如在自动控制[3] [4]，医学[5] [6]，数据测绘[7]-[10]等其他领域[11]-[15]。对于大型线性系统的求解方法。随机Kaczmarz算法是一种比较有效的解决办法。Strohmer和Vershynin [1]提出了随机Kaczmarz算法。遵循RK的方法，Bai和Wu提出了贪婪随机化Kaczmarz方法(RGRK) [16]，在[2]中作者使用K-均值聚类来划分系数矩阵，同时采用了GRK中使用的贪心技术，构造了求解线性系统(1)的K-均值聚类确定的随机块Kaczmarz方法(RBK(k))。本文在(RBK(k))的基础上进行了改进。在每一次的选代中根据$i_k=\arg {\max }_{1\le i\le k}{\Vert b{\nu }_i-A{\nu }_i{x}_k\Vert }_2^2$ 选择工作行块${\mathcal{V}}_{i_k}$ 确保首先消除最大的残差。在本文中提出了两种算法并给出了收敛性的证明，并在最后的数值实验部分与RBK(k)进行了比较，实验结果表明了算法的优越性。

2. 带K-均值最大残差块方法

算法1是文献[2]提出的随机块Kaczmarz方法，该算法如下。

算法1：带K-均值随机块Kaczmarz方法(RBK(k))

输入：$A,b,x_0,k,l$ 和$\theta$

输出：${\text{x}}_{\text{l}}$

1：对A和b使用K-均值方法，得到行索引$1,\cdots ,m$ 的划分$\left\{v_1,\cdots ,v_k\right\}$ ，得到k组$A_{v_i}$ 和$b_{v_i}$ ，$i=1,2,\cdots ,k$ 。$A_{v_i}$ 和$b_{v_i}$ 再以紧凑形式分配即：$\overline{A}=\left[{\overline{A}}_{v_1};\cdots ;{\overline{A}}_{v_k}\right]$ 和$\overline{b}=\left[\begin{array}{c}{\overline{b}}_{v_1};\cdots ;{\overline{b}}_{v_k}\end{array}\right]$

2：对$j=1,2,\cdots ,l$ 做

3：计算：${ϵ}_j=\frac{\theta }{{\Vert \overline{b}-\overline{A}{x}_j\Vert }_2^2}\underset{1\le v_j\le k}{\max }\left\{\frac{{\left|{\overline{b}}_{v_j}-{\overline{A}}_{v_j}{x}_j\right|}^2}{{\Vert {\overline{A}}_{v_j}\Vert }_2^2}\right\}+\frac{1-\theta }{{\Vert A\Vert }_F^2},\theta \in \left(0,1\right)$

4：令$U_j=\left\{v_j|{\left|{\overline{b}}_{v_j}-{\overline{A}}_{v_j}{x}_j\right|}^2\ge {ϵ}_j{\Vert \overline{b}-\overline{A}{x}_j\Vert }_2^2{\Vert {\overline{A}}_{v_j}\Vert }_2^2\right\}$

5：根据下面规则计算${\tilde{r}}_j^{\left(v\right)}$

${\tilde{r}}_j^{\left(v\right)}=\{\begin{array}{ll}{\overline{b}}_v-{\overline{A}}_v{x}_j,\hfill & 如果v\in U_j\hfill \\ 0,\hfill & 其他\hfill \end{array}$

6：依概率$\mathrm{Pr}\left(v=v_j\right)=\frac{{\left|{\tilde{r}}_j^{\left(v\right)}\right|}^2}{{\Vert {\tilde{r}}_j\Vert }_2^2}$ 选择$v_j\in U_j$

7：计算$x_{j+1}=x_j+A_{v_j}^{†}\left(b_{v_j}-A_{v_j}{x}_j\right)$

8：结束

算法2是在算法1的基础上，直接选择了具有最大残差的行块。

算法2：K-均值最大残差块方法

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

输出：$x_l$

1：使用K-均值方法对A和b处理，得到行索引$1,\cdots ,m$ 的划分$\left\{v_1,\cdots ,v_k\right\}$ ，就有k组$A_{v_i}$ 和$b_{v_i}$ ，$i=1,2,\cdots ,k$ 。

2：对$j=1,2,\cdots ,l$ 做

3：选择$i_j=\arg {\max }_{1\le i\le k}{\Vert b{\nu }_i-A{\nu }_i{x}_j\Vert }_2^2$

4：计算$x_{j+1}=x_j+A_{{\mathcal{V}}_{i_j}}^{†}\left(b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\right)$

5：结束

算法3是在算法2的基础上，避免了求伪逆的计算，大大节省了计算成本。

算法3：K-均值最大残差块的加速方法

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

输出：$x_l$

1：使用K-均值方法对A和b处理，得到行索引$1,\cdots ,m$ 的划分$\left\{v_1,\cdots ,v_k\right\}$ ，就有k组$A_{v_i}$ 和$b_{v_i}$ ，$i=1,2,\cdots ,k$ 。

2：对$j=1,2,\cdots ,l$ 做

3：选择$i_j=\arg {\max }_{1\le i\le k}{\Vert b{\nu }_i-A{\nu }_i{x}_j\Vert }_2^2$

4：计算${\alpha }_j=\omega \frac{{\Vert b_{{\nu }_{i_j}}-A{\nu }_{i_j}{x}_j\Vert }_2^2{\Vert A{\nu }_{i_j}\Vert }_F^2}{{\Vert A_{{\nu }_{i_j}}^{\text{T}}\left(b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\right)\Vert }_2^2}$

5：计算$x_{j+1}=x_j+{\alpha }_j\frac{A_{{\nu }_{i_j}}^{\text{T}}\left(b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\right)}{{\Vert A{\nu }_{i_j}\Vert }_F^2}$

6：结束

3. 算法2和算法3的收敛性分析

引理1 [17]：设向量$z\in ℛ\left(A^{\text{T}}\right)$ ，则有${\Vert Az\Vert }_2^2\ge {\sigma }_{\min }^2\left(A\right){\Vert z\Vert }_2^2$ 。

定理1：对于相容线性系统(1)，由算法2选代生成的序列${\left\{x_j\right\}}_{i=0}^{\infty }$ 收敛到$x_{\star }=A^{†}b$ ，对于任意的$j\ge 0$ ，满足下列关系：

${\Vert x_1-x_{\star }\Vert }_2^2\le \left(1-\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)k}\right){\Vert x_0-x_{\star }\Vert }_2^2$

以及

${\Vert x_{j+1}-x_{\star }\Vert }_2^2\le {\left(1-\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\right)}^j\left(1-\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)k}\right){\Vert x_0-x_{\star }\Vert }_2^2$

证明：

$x_{j+1}-x_{\star }=x_j-x_{\star }+A_{{\mathcal{V}}_{i_j}}^{†}\left(b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\right)$

$A_{{\mathcal{V}}_{i_j}}\left(x_{j+1}-x_{\star }\right)=A_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)-A_{{\mathcal{V}}_{i_j}}{A}_{{\mathcal{V}}_{i_j}}^{†}{A}_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)=0$

所以

$A_{{\mathcal{V}}_{i_j}}^{†}{A}_{{\mathcal{V}}_{i_j}}\left(x_{j+1}-x_{\star }\right)=0$

由于

$x_{j+1}-x_j=A_{{\nu }_{i_j}}^{†}{A}_{{\nu }_{i_j}}\left(x_{\star }-x_j\right)$

所以

$x_{j+1}-x_j$ 与$x_{j+1}-x_{\star }$ 正交

${\Vert x_{j+1}-x_{\star }\Vert }_2^2={\Vert x_j-x_{\star }\Vert }_2^2-{\Vert A_{{\nu }_{i_j}}^{†}{A}_{{\nu }_{i_j}}\left(x_j-x_{\star }\right)\Vert }_2^2$

下面第一个等式是根据引理1得到的

$\begin{array}{c}{\Vert A_{{\nu }_{i_j}}^{†}A{\nu }_{i_j}\left(x_j-x_{\star }\right)\Vert }_2^2\ge {\sigma }_{\min }^2\left(A_{{\mathcal{V}}_{i_j}}^{†}\right){\Vert A_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)\Vert }_2^2\\ =\frac{1}{{\sigma }_{\max }^2\left(A_{{\mathcal{V}}_{i_j}}\right)}{\Vert A_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)\Vert }_2^2\\ \ge \frac{1}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}{\Vert A_{{\nu }_{i_j}}\left(x_j-x_{\star }\right)\Vert }_2^2\\ =\frac{1}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}{\Vert b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\Vert }_2^2\\ =\frac{1}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}\underset{1\le i\le k}{\max }{\Vert b_{{\nu }_i}-A_{{\nu }_i}{x}_j\Vert }_2^2\end{array}$

此外

$\begin{array}{c}b{\nu }_{i_j}-A_{{\nu }_{i_j}}{x}_{j+1}=b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}\left(x_j+A_{{\mathcal{V}}_{i_j}}^{†}\left(b_{{\mathcal{V}}_i}-A_{{\mathcal{V}}_i}{x}_j\right)\right)\\ =b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j-A_{{\mathcal{V}}_{i_j}}{A}_{{\mathcal{V}}_{i_j}}^{†}\left(b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\right)\\ =b_{{\mathcal{V}}_{ij}}-A_{{\mathcal{V}}_{ij}}{x}_j-A_{{\mathcal{V}}_{ij}}{A}_{{\mathcal{V}}_{ij}}^{†}{b}_{{\mathcal{V}}_{ij}}+A_{{\mathcal{V}}_{ij}}{A}_{{\mathcal{V}}_{ij}}^{†}{A}_{{\mathcal{V}}_{ij}}{x}_j\\ =A_{V_{i_j}}{x}_{\star }-A_{V_{i_j}}{A}_{V_{i_j}}^{†}{A}_{V_{i_j}}{x}_{\star }\\ =A_{{\mathcal{V}}_{i_j}}{x}_{\star }-A_{{\mathcal{V}}_{i_j}}{x}_{\star }\\ =0\end{array}$

因此我们可以得到对于$j=1,2,\cdots$ 有

${\Vert b-A{x}_j\Vert }_2^2={\displaystyle \sum }_{{\mathcal{V}}_{i_j}\in \mathcal{V}-{\mathcal{V}}_{i_{j-1}}}{\Vert b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\Vert }_2^2\le \left(k-1\right)\underset{1\le i\le k}{\max }{\Vert b_{{\mathcal{V}}_i}-A_{{\mathcal{V}}_i}{x}_j\Vert }_2^2$

当$j=0$ 我们还可以得到

${\Vert b-A{x}_0\Vert }_2^2\le k\underset{1\le i\le k}{\max }{\Vert b_{{\mathcal{V}}_i}-A_{{\mathcal{V}}_i}{x}_0\Vert }_2^2$

因此有

$\underset{1\le i\le k}{\max }{\Vert b_{{\mathcal{V}}_i}-A_{{\mathcal{V}}_i}{x}_j\Vert }_2^2\ge \frac{1}{k-1}{\Vert b-A{x}_j\Vert }_2^2$

综上所述，我们最终就可以得到

$\begin{array}{c}{\Vert x_{j+1}-x_{\star }\Vert }_2^2\le {\Vert x_j-x_{\star }\Vert }_2^2-\frac{1}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}\frac{{\Vert b-A{x}_j\Vert }_2^2}{k-1}\\ ={\Vert x_j-x_{\star }\Vert }_2^2-\frac{1}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}\frac{{\Vert A\left(x_j-x_{\star }\right)\Vert }_2}{k-1}\\ \le {\Vert x_j-x_{\star }\Vert }_2^2-\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}{\Vert x_j-x_{\star }\Vert }_2^2\\ =\left(1-\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\right){\Vert x_j-x_{\star }\Vert }_2^2\end{array}$

定理2：对于对于相容线性系统(1)，设$\omega \in \left(0,2\right)$ 由算法3选代生成的序列${\left\{x_j\right\}}_{j=0}^{\infty }$ 收敛到$x_{\star }=A^{†}b$ ，对于任意的$j\ge 0$ ，满足下列关系：

${\Vert x_1-x_{\star }\Vert }_2^2\le \left(1-\left(2\omega -{\omega }^2\right)\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)k}\right){\Vert x_0-x_{\star }\Vert }_2^2$

以及

${\Vert x_{j+1}-x_{\star }\Vert }_2^2\le {\left(1-\left(2\omega -{\omega }^2\right)\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\right)}^j\left(1-\left(2\omega -{\omega }^2\right)\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)k}\right){\Vert x_0-x_{\star }\Vert }_2^2$

证明：从算法3的第5步和$A_{{\mathcal{V}}_{i_j}}{x}_j=b_{{\mathcal{V}}_{i_j}}$ ，我们可以得到

$\begin{array}{c}{x}_{j+1}-x_{\star }=x_j-x_{\star }-{\alpha }_j\frac{A_{{\nu }_{i_j}}^{\text{T}}\left(b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\right)}{{\Vert A{\nu }_{i_j}\Vert }_F^2}\\ =x_j-x_{\star }-{\alpha }_j\frac{A_{{\nu }_{i_j}}^{\text{T}}{A}_{{\nu }_{i_j}}\left(x_j-x_{\star }\right)}{{\Vert A_{{\nu }_{i_j}}\Vert }_F^2}\\ =\left(I-{\alpha }_j\frac{A_{{\mathcal{V}}_{i_j}}^{\text{T}}A{v}_{i_j}}{{\Vert A_{{\mathcal{V}}_{i_j}}\Vert }_F^2}\right)\left(x_j-x_{\star }\right)\end{array}$

对上述等式取欧几里得范数的平方，可以得到

$\begin{array}{c}{\Vert \begin{array}{c}{x}_{j+1}-x_{\star }\end{array}\Vert }_2^2={\begin{array}{c}\Vert \left(I-{\alpha }_j\frac{A_{{\mathcal{V}}_{i_j}}^{\text{T}}A{v}_{i_j}}{{\begin{array}{c}\Vert A_{{\mathcal{V}}_{i_j}}\Vert \end{array}}_F^2}\right)\left(x_j-x_{\star }\right)\Vert \end{array}}_2^2\\ ={\Vert \begin{array}{c}{x}_j-x_{\star }\end{array}\Vert }_2^2-2{\alpha }_j\frac{A_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)}{{\Vert A_{{\mathcal{V}}_{i_j}}\Vert }_F^2}+{\alpha }_j^2\frac{A_{{\mathcal{V}}_{i_j}}^{\text{T}}{A}_{{\mathcal{V}}_{i_j}}\left(x_j-x_{\star }\right)}{{\Vert \begin{array}{c}{A}_{{\mathcal{V}}_{i_j}}\end{array}\Vert }_F^2}\end{array}$

将${\alpha }_j$ 代入上式可以得到

$\begin{array}{c}{\Vert x_{j+1}-x_{\star }\Vert }_2^2={\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{{\Vert b_{{\nu }_{i_j}}-A{\nu }_{i_j}{x}_j\Vert }_2^4}{{\Vert A_{{\nu }_{i_j}}^{\text{T}}\left(b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\right)\Vert }_2^2}\\ \le {\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{{\Vert b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\Vert }_2^4}{{\sigma }_{\max }^2\left(A{\nu }_{i_j}\right){\Vert b{\nu }_{i_j}-A{\nu }_{i_j}{x}_j\Vert }_2^2}\\ \le {\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{{\Vert b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_j}{x}_j\Vert }_2^2}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}\\ ={\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{\underset{1\le i\le k}{\max }{\Vert b_{{\mathcal{V}}_i}-A_{{\mathcal{V}}_i}{x}_j\Vert }_2^2}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)}\end{array}$

第一个不等式成立因为$2\omega -{\omega }^2>0$ ，$\omega \in \left(0,2\right)$ 以及根据引理1可以得到的下面不等式

${\Vert A_{{\mathcal{V}}_{i_j}}^{\text{T}}\left(b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\right)\Vert }_2^2\le {\sigma }_{\max }^2\left(A_{{\mathcal{V}}_{i_j}}\right){\Vert b_{{\mathcal{V}}_{i_j}}-A_{{\mathcal{V}}_{i_j}}{x}_j\Vert }_2^2$

从定理1的证明可以看出，对于$i=1,2,\cdots$ 有下面的式子成立0

$\begin{array}{c}{\Vert x_{j+1}-x_{\star }\Vert }_2^2\le {\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{{\Vert b-A{x}_j\Vert }_2^2}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\\ \le {\Vert x_j-x_{\star }\Vert }_2^2-\left(2\omega -{\omega }^2\right)\frac{{\sigma }_{\min }^2\left(A\right){\Vert x_j-x_{\star }\Vert }_2^2}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\\ =\left(1-\left(2\omega -{\omega }^2\right)\frac{{\sigma }_{\min }^2\left(A\right)}{{\sigma }_{\max }\left(A_v{A}_v^{\text{T}}\right)\left(k-1\right)}\right){\Vert x_j-x_{\star }\Vert }_2^2\end{array}$

4. 数值实验

用算法1~3方法求解线性系统(1)的例子是由Matlab函数randn服从标准正态分布随机生成的矩阵，以及从文献中选取的5个具有应用背景的稀疏矩阵[18]。分别是Franz1、mk11 b2、Trefethen-300、gen4、p6000。

用IT表示选代次数，CPU表示运行时间(单位：秒)，表格中的IT和CPU是循环运行5次取的平均值，为了更加直观地显示算法2和算法3较算法-1的有效性。我们定义$\text{speed-up}1=\frac{\text{CPU}\text{of}\text{Alogrithm}2}{\text{CPU}\text{of}\text{Alogrithm}1}$ ，$\text{speed-up}2=\frac{\text{CPU}\text{of}\text{Alogrithm}3}{\text{CPU}\text{of}\text{Alogrithm}1}$ 。对于选取的具有应用背景的矩阵，我们也给出了这些矩阵的基本信息，其中cond(A)和density(A)分别表示矩阵的条件数和稠密性。从表1可以知道该信息。我们定义$\text{density}=\frac{\text{number}\text{of}\text{nonzeros}\text{of}\text{an}m\times n\text{matrix}}{\text{mn}}$ 。在数值实验中，解向量是由$x_{\star }$ 是由Matlab函数randn随机生成的，其元素是遵循独立的标准正态分布。根据右侧向量b，我们取$b=A{x}_{\star }$ ，所有计算的初始向量从$x_0=0$ 开始的。我们取RES表示选代值和真实值的相对误差设置停机准则为$\text{RSE}=\frac{{\Vert x_j-x_{\star }\Vert }_2^2}{{\Vert x_{\star }\Vert }_2^2}\le {10}^{-6}$ 或超过规定的最大选代次数$\text{IT}=200000$ 。

在下列表中我们分别列出了算法1~3方法的选代步数(IT)和计算时间(CPU)。对于随机生成的矩阵，从表2，表3我们可以看出算法2和算法3的方法在选代步数(IT)和计算时间(CPU)都是优于算法1的方法。来自Davis和Hu工作的稀疏列满秩矩阵，我们可以从表4可以看出，算法2和算法3在比较的时候是有竞争力的。

Table 1. Specific matrix related information

表1. 具体矩阵的相关信息

矩阵	Franz1	mk11-b2	Trefethen-300	gen4	p6000
m × n	2240 × 768	6390 × 990	300 × 300	1537 × 4298	2095 × 7967
density	0.003	0.003	0.052	0.0162	0.0012
cond(A)	7.13E+15	5.04E+15	1772.7	39.19	4710

Table 2. Randomly generated matrix A, numerical results for solving linear systems when m = 8000 and n is different

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

m × n		8000 × 100	8000 × 200	8000 × 500	8000 × 1000	8000 × 2000
RBK(100)	IT	4	18.2	65	148.6	867.9
	CPU	0.3073	0.4655	0.7049	1.2617	3.444
MRBK(100)	IT	2	2	2	2	867.9
	CPU	0.2594	0.406	0.7016	1.6365	3.444
MARBK(100)	IT	13.2	13	10	14.2	1711.6
	CPU	0.1973	0.3211	0.5318	1.0977	2.4463
RBK(200)	IT	20.6	52	150.2	589.1	310
	CPU	0.3884	0.5219	0.8537	1.142	1.8901
MRBK(200)	IT	2	2	2	1136.7	2
	CPU	0.2792	0.457	0.7249	1.1415	1.8686
MARBK(200)	IT	17.2	15	12.6	1136.7	16.2
	CPU	0.2111	0.3506	0.5740	1.1415	1.3085
RBK(400)	IT	68.2	140.6	289.2	1136.7	604
	CPU	0.4616	0.705	1.3092	1.1415	2.9269
MRBK(400)	IT	2	2	2	1136.7	2
	CPU	0.3494	0.4697	0.8575	1.1415	1.797
MARBK(400)	IT	19.2	19	14.2	1136.7	19
	CPU	0.3	0.3436	0.5987	1.1415	1.3888

Table 3. Randomly generated matrix A, numerical results for solving linear systems when n = 8000 and m is different

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

m × n		100 × 8000	200 × 8000	500 × 8000	1000 × 8000	2000 × 8000
RBK(2)	IT	5	6	7	9	13
	CPU	0.0655	0.1124	0.4639	1.1651	2.976
MRBK(2)	IT	5	6	7	9	14
	CPU	0.0503	0.1108	0.3963	0.8916	2.4732
MARBK(2)	IT	8	9	12	17	27
	CPU	0.043	0.0892	0.3412	0.8569	2.4311
RBK(4)	IT	10.4	12.8	18	22	34.6
	CPU	0.0669	0.1505	0.5937	1.3685	3.3253
MRBK(4)	IT	8	9.6	12	15.2	23.2
	CPU	0.0491	0.0992	0.4762	0.9208	2.539
MARBK(4)	IT	11.2	12.2	16.8	22	36.8
	CPU	0.0440	0.0993	0.4071	0.8747	2.8612
RBK(6)	IT	16.4	20.8	28	32.2	57
	CPU	0.0519	0.1121	0.5102	1.3914	3.649
MRBK(6)	IT	12	13.2	15.8	22.2	36.4
	CPU	0.0517	0.1165	0.3808	1.0045	2.757
MARBK(6)	IT	13.4	15.8	20.4	28.4	48.2
	CPU	0.0419	0.0993	0.3554	0.932	2.8281

Table 4. Numerical results of solving linear systems using specific matrices

表4. 具体矩阵求解线性系统的数值结果

矩阵		Franz1	mk11-b2	Trefethen-300	gen4	p6000
RBK(k)	IT	15.4	5	52	521.8	30
	CPU	2.2341	0.8508	0.1003	1.1908	1.7317
MRBK(k)	IT	4	3.8	52	399.4	29.4
	CPU	0.2671	0.8474	0.0117	1.081	1.71
MARBK(k)	IT	34.4	10	74	931.4	49.6
	CPU	0.0075	0.0362	0.0063	0.8079	0.0244
k	/	4	4	20	4	4
speed-up1	/	8.36428304	1.00401227	8.57264957	1.10157262	1.01269006
speed-up2	/	297.88	23.5027624	15.9206349	1.4739448	70.9713115

5. 总结

本文在RBK(k)算法基础上进行改进构建了最大残差块Kaczmarz算法及其加速版本，分析了算法的收敛性。数值实验证实了该算法的有效性。最优k值的选择仍是未来值得研究的方向。