巫文婷

（北京理工大学数学与统计学院，北京 100081）

对于系数矩阵为A∈Cm×n且右端向量为b∈Cm的大规模相容线性代数方程组

的求解，Kaczmarz方法［1］是经典的行处理迭代方法，在信号与图像处理领域有着广泛的应用。其每步迭代只需按照给定的循环顺序选取系数矩阵的某一行，并将当前迭代向量正交投影至由该行所形成的超平面上。Strohmer等［2］提出按照与系数矩阵每一行的欧氏范数平方成比例的概率准则随机选取系数矩阵的行，得到了收敛更为快速的随机Kaczmarz方法。若用(·)*表示相应矩阵或向量的共轭转置，则当初始迭代向量在A*的列空间中时，随机Kaczmarz方法期望线性收敛［2-5］到线性代数方程组（1）的最小欧氏范数解x⋆=A†b，其中A†表示系数矩阵A的Moore-Penrose伪逆。当线性代数方程组（1）的右端向量发生扰动时，Needell［6］给出了随机Kaczmarz方法的期望解误差的上界，并说明了随着迭代步数的增长，随机Kaczmarz方法的期望解误差会以线性速率下降至一个误差阈值。之后，Zouzias等［7］对Needell所提出的期望解误差的上界进行了改进。

影响随机Kaczmarz方法收敛速率的关键在于其中所蕴含的用于选取每步迭代所需调用的系数矩阵行的概率准则。为了提高随机Kaczmarz方法的收敛速率，Bai等［8-9］提出了一个可以获取每步迭代中残向量的模较大分量的概率准则，并基于该概率准则构造出了贪婪随机Kaczmarz方法。当初始迭代向量在A*的列空间时，贪婪随机Kaczmarz方法所产生的迭代序列收敛到线性代数方程组（1）的最小范数解x⋆=A†b且具有期望线性收敛速率

当线性代数方程组（1）的右端向量b被加上一个不为零的扰动向量r∈Cm时，实际求解的线性代数方程组问题将变为

其中，y=b+r。若对任意矩阵G∈Cm×n和任意向量u∈Cm，用G(i)和u(i)分别表示矩阵G的第i行和向量u的第i个分量，由于贪婪随机Kaczmarz方法的第k步迭代将当前迭代向量x k投影至超平面=上，而线性代数方程组（1）的最小范数解x⋆=A†b满足Ax⋆=y-r，故其并不能收敛到x⋆。在这种情况下，本文对贪婪随机Kaczmarz方法的期望解误差进行分析，得到了其上界，证明了随着迭代步数的增长，由贪婪随机Kaczmarz方法所产生的迭代解与原线性代数方程组（1）的最小范数解x⋆之间的误差以线性速率下降至一个给定阈值。数值实验表明本文所给出的阈值能够很好地估计贪婪随机Kaczmarz方法的迭代解误差所能达到的最小值。

1 贪婪随机Kaczmarz方法

不失一般性，在本文的讨论中，总是假设系数矩阵A不存在零行，即A(i)≠0，i=1，2，…，m。当线性代数方程组（1）的右端向量发生扰动时，求解线性代数方程组（2）的贪婪随机Kaczmarz方法［8］如下：

输入：A，y，ℓ与x0

输出：xℓ

1：fork=0，1，…，ℓ-1 do

2：计算

3： 确定正整数指标集

4： 计算向量r͂k∈Cm的第i个分量

7：endfor

令R(A)和R(A)⊥分别为系数矩阵A的像空间和其像空间的正交补子空间，则有r=r R(A)+r R(A)⊥，其中r R(A)和r R(A)⊥分别表示向量r在R(A)和R(A)⊥上的正交投影。若线性代数方程组（1）右端向量的扰动r在系数矩阵A的像空间R(A)中，则线性代数方程组（2）等价于线性代数方程组

易知线性代数方程组（3）是相容的且其最小范数解为=A†(b+r R(A))。此时，若初始迭代向量在A*的列空间中，则求解线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代序列期望线性收敛到x͂⋆。

2 误差分析

与求解原线性代数方程组（1）的贪婪随机Kaczmarz方法的分析［8］类似，根据求解线性代数方程组（2）的贪婪随机Kaczmarz方法中εk的定义可知，对于k=1，2，…，由于

有

而对于k=0，有

因此，可得引理1。

引理1求解系数矩阵为A∈Cm×n且右端向量为y∈Cm的线性代数方程组（2）的贪婪随机Kaczmarz方法的概率准则中的量εk，k=0，1，2，…，满足

和

当线性代数方程组（1）右端向量的扰动r在R(A)中时，若初始迭代向量在A*的列空间中，则求解线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代向量期望线性收敛到线性代数方程组（3）的最小范数解x͂⋆。对于一般的扰动情形，求解线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代向量与x͂⋆之间的期望误差的上界可以由引理2给出。

引理2如果初始迭代向量x0∈Cn在A*的列空间中，则贪婪随机Kaczmarz方法求解扰动后的线性代数方程组（2）所产生的迭代序列与线性代数方程组（3）的最小范数解之间的期望误差满足

其中

证明：由求解线性代数方程组（2）的贪婪随机Kaczmarz方法的定义和Ax͂⋆=b+r R(A)可知

其中I表示具有合适阶数的单位矩阵，则有

用Ek表示固定前k步迭代的条件期望，即Ek[·]=E[·|i0，i1，…，ik-1]，其中il(l=0，1，…，k-1)为第l步迭代所选取的行指标，则对等式（5）两边取条件期望可得

由于

可知

则根据正整数指标集Uk的定义可得

对于在矩阵A*的像空间R(A*)中的任意向量u，有不等式

成立。利用该不等式，通过引理1、向量b+r R(A)-Ax k=A(-x k)∈R(A)和向量r R(A)⊥的正交性、x k-∈R(A*)以 及可 知，对 于k=1，2，…，有

而对于k=0，有

对不等式两边取期望可得对于k=1，2，…，有

而对于k=0，有

因此，由于0＜α＜1，β≥0，对于k=1，2，…，有

由Jensen不等式可知

基于引理2，关于求解扰动后的线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代近似解与原线性代数方程组（1）的最小范数解x⋆之间的期望误差的上界，可以给出如下定理。

定理1如果初始迭代向量x0∈Cn在A*的列空间中且为线性代数方程组（3）的最小范数解，则贪婪随机Kaczmarz方法求解扰动后的线性代数方程组（2）所产生的迭代序列与原线性代数方程组（1）的最小范数解x⋆=A†b之间的期望误差满足

其中，α、α0和β如（4）式中定义。

证明：因为-x⋆∈R(A*)，则通过不等式（6）可知

又因为有

成立，则由引理2可得

定理1说明了当迭代步数k趋于无穷时，贪婪随机Kaczmarz方法求解扰动后的线性代数方程组（2）所产生的迭代序列的期望相对解误差以的线性速率下降至某个阈值，并给出了该阈值的估计

若r R(A)⊥为零，则扰动后的线性代数方程组（2）即为相容的线性代数方程组（3）。因此，求解线性代数方程组（2）的贪婪随机Kaczmarz方法的迭代序列收敛到线性代数方程组（3）的最小范数解，且其与原线性代数方程组（1）的最小范数解x⋆的差趋于-x⋆。此时，如果初始迭代向量在A*的列空间中，则求解线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代序列满足

若r R(A)为零，则线性代数方程组（3）即为线性代数方程组（1），且线性代数方程组（3）的最小范数解即为x⋆=A†b。此时，如果初始迭代向量在A*的列空间中，则求解线性代数方程组（2）的贪婪随机Kaczmarz方法所产生的迭代序列满足

3 数值实验

通过数值实验测试贪婪随机Kaczmarz方法求解带右端项扰动的线性代数方程组（2）时的数值表现，并将所估计的阈值与贪婪随机Kaczmarz方法所产生的真实迭代解误差进行比较。

所测试的系数矩阵A分为两类：一类为利用MATLAB函数randn随机产生的元素服从标准正态分布的矩阵，矩阵维数分别为200×100 000和100000×200；另一类为取自稀疏矩阵库［10］的稀疏矩阵bibd_16_8和ash958，矩阵维数分别为120×12870和958×292。

线性代数方程组（1）的右端向量b取为Ax*，其中x*∈Rn是利用MATLAB函数randn随机产生的。线性代数方程组（1）的右端项的扰动向量r分为三类：一类为随机扰动，一类为在系数矩阵的像空间中的扰动，另一类为在系数矩阵的像空间的正交补空间中的扰动。这三类扰动均由MATLAB函数randn生成，并且其欧氏范数为原始右端向量b的欧氏范数的0.0005倍。由于200×100000的随机矩阵和稀疏矩阵bibd_16_8为行满秩矩阵，故其所对应的右端项扰动r一定在其像空间中。所有计算均从初始向量x0=0开始。

图1描绘了当线性代数方程组（1）的系数矩阵为随机产生的矩阵且右端项扰动r为随机扰动时，贪婪随机Kaczmarz方法重复运行50次所产生的相对解误差ERS的中值和阈值τ分别取以10为底的对数之后相对于迭代步数的曲线，其中

从图1可以看出，贪婪随机Kaczmarz方法所产生的相对解误差下降至10-3左右之后不再减小，且阈值τ能够很好地给出该最小值的估计。特别地，当系数矩阵为200×100 000的随机矩阵时，τ与真实相对解误差所能达到的最小值非常接近。

图1 当m=200，n=100 000和m=100 000，n=200时，lg E RS和lgτ相对于迭代步数的曲线Fig.1 Pictures of lg E RS and lgτversus the iteration step when m=200,n=100000 and m=100000,n=200

图2描绘了对于100 000×200的随机产生的系数矩阵，当线性代数方程组（1）的右端项扰动r分别在系数矩阵的像空间和系数矩阵的像空间的正交补空间中时，贪婪随机Kaczmarz方法重复运行50次所产生的相对解误差的中值和阈值τ分别取以10为底的对数之后相对于迭代步数的曲线。从图2可以看出，贪婪随机Kaczmarz方法所产生的相对解误差下降至10-3左右之后不再减小，且阈值τ能够很好地给出该最小值的估计，特别当扰动向量r在系数矩阵的像空间中时。

图2 当m=100 000，n=200且r∈R(A)和r∈R(A)⊥时，lg E RS和lgτ相对于迭代步数的曲线Fig.2 Pictures of lg E RS and lgτversus the iteration step when m=100 000，n=200,and r∈R(A)and r∈R(A)⊥

当线性代数方程组（1）的系数矩阵为稀疏矩阵bibd_16_8和ash958时，图3描绘了右端项扰动r为随机扰动时，贪婪随机Kaczmarz方法重复运行50次所产生的相对解误差的中值和阈值τ分别取以10为底的对数之后相对于迭代步数的曲线；图4描绘了右端项扰动r分别在系数矩阵的像空间和系数矩阵的像空间的正交补空间中时的相应曲线。类似地，从图3和图4可以看出，贪婪随机Kaczmarz方法所产生的相对解误差下降至10-3左右之后不再减小，且阈值τ能够很好地给出该最小值的估计，特别当系数矩阵为bibd_16_8和当系数矩阵为ash958且扰动向量r在系数矩阵的像空间中时。

图3 当矩阵为bibd_16_8和ash958时，lg E RS和lgτ相对于迭代步数的曲线Fig.3 Pictures of lg E RS and lgτversus the iteration step when the matrices are bibd_16_8 and ash958

图4 当矩阵为ash958且r∈R(A)和r∈R(A)⊥时，lg E RS与lgτ相对于迭代步数的曲线Fig.4 Pictures of lg E RS and lgτversus the iteration step when the matrix is ash958,and r∈R(A)and r∈R(A)⊥

4 结论

当线性代数方程组的右端向量发生扰动时，证明了贪婪随机Kaczmarz方法的期望解误差以线性速率下降至一个给定阈值。数值实验表明该阈值能够很好地估计贪婪随机Kaczmarz方法的迭代解误差所能达到的最小值。可以发现当扰动不是很大且对解的精度要求不是很高时，贪婪随机Kaczmarz方法仍然能够很好地给出原线性代数方程组最小范数解的近似。