基于Count Sketch的预处理贪婪Kaczmarz方法

摘要

在贪婪Kaczmarz方法中， 通过对系数矩阵进行正交三角分解引入右预处理子能够提高贪婪Kaczmarz方法的收敛速率。但在系数矩阵的行数远大于列数的情况下， 正交三角分解的成本过高。为降低预处理的成本， 通过引入Count Sketch变换， 提出了基于Count Sketch的预处理贪婪Kaczmarz方法， 并对新方法进行了收敛性分析。理论分析说明了新方法在系数矩阵条件数较大时比已有方法具有更好的收敛速率。数值实验验证了新方法的有效性。

关键词

Kaczmarz方法; 预处理; Count Sketch; 收敛性

考虑求解大规模超定相容线性方程组

$$\mathit{A}\mathit{x}=\mathit{b}$$

（1）

其中， 矩阵$\mathit{A}\in {\mathbf{R}}^{m\times n}$， $m>n$， $\mathit{b}\in {\mathbf{R}}^m$， $\mathit{x}\in {\mathbf{R}}^n$为未知向量。Kaczmarz方法^［1］是求解线性方程组（1）的一种经典的行作用迭代法， 在分布式计算、计算机断层扫描、数字信号处理等领域均有广泛应用^［2-4］。假定初始近似值为$x_{\mathrm{0}}$， 那么经典的Kaczmarz方法可以表示为

$${\mathit{x}}_{k+\mathrm{1}}={\mathit{x}}_k+\frac{{\mathit{b}}_{i_k}-{\mathit{a}}_{i_k}^{\mathbf{T}}{\mathit{x}}_k}{{‖{\mathit{a}}_{i_k}‖}_{\mathrm{2}}^{\mathrm{2}}}{\mathit{a}}_{i_k},\hspace{1em}k=\mathrm{0,1},\mathrm{2}\cdots$$

（2）

式中：$\Vert \cdot {\Vert }_{\mathrm{2}}$表示欧几里得范数； ${\mathit{a}}_{i_k}$为矩阵$\mathit{A}$的第$i_k$行； ${\mathit{b}}_{i_k}$为向量$\mathit{b}$的第$i_k$个元素， $i_k=(k\mathrm{m}\mathrm{o}\mathrm{d}m)+\mathrm{1}$。

2009年， Strohmer和Vershynin^［5］首次证明了当系数矩阵列满秩时， 随机Kaczmarz方法的线性收敛速率， 且该速率只取决于系数矩阵的缩放条件数$\kappa (\mathit{A})=\Vert \mathit{A}{\Vert }_F\Vert {\mathit{A}}^{-\mathrm{1}}{\Vert }_{\mathrm{2}}$， 这再度使得Kaczmarz方法成为科学计算领域的一个热点。后来， Nutini等人^［6］提出基于最大距离规则的贪婪Kaczmarz方法， 它在每次迭代中选取与当前迭代解距离最大的超平面作为工作行。 Bai和Wu^［7］结合贪婪与随机的思想， 提出了贪婪随机Kaczmarz方法， 它在迭代过程中动态跟踪残差向量， 大大提高了随机Kaczmarz方法的效率。更多Kaczmarz 类方法的研究参见文献^［8-11］。

Sketching方法是将Sketch矩阵$\mathit{S}\in {\mathit{R}}^{d\times m}$作用于原始矩阵$\mathit{A}$的过程， 其中$\mathit{S}\in {\mathit{R}}^{d\times m}$是从矩阵的某个分布中指定的一个的随机矩阵， 且$n<d\ll m$ ^［8］ 。该方法使高成本的运算可以在较小的矩阵$\widehat{\mathit{A}}=\mathit{S}\mathit{A}\in {\mathit{R}}^{d\times n}$上执行， 从而加快原始问题的求解。近年来， 有很多sketch矩阵被提出并加以研究。例如： Gaussian Sketch^［9］， 子采样随机Hadamard变换^［10］， Count Sketch^［11-12］以及稀疏Johnson-Lindenstrauss变换^［13］等。最近， 有学者将Sketching方法与Kaczmarz类方法结合起来， 提高Kaczmarz类方法的效率。Katrutsa等^［14］通过行抽样实现sketching过程， 对变换后的矩阵通过正交三角分解得到右预处理子， 再用Kaczmarz方法进行求解， 大大减少预处理时间， 同时加速了Kaczmarz方法的收敛。Needell等^［15］将Gaussian sketch与分块Kaczmarz方法结合， 提出分块Gaussian Kaczmarz方法， 说明了某些特定情况下使用分块Gaussian Kaczmarz方法求解线性方程组的优势。Zhang等^［16］将Count Sketch与最大加权残差Kaczmarz方法结合起来， 提出一种求解高度超定线性方程组的快速随机算法，大大提高了算法的计算效率。

在基于sketching预处理Kaczmarz方法^［14］的基础之上， 结合贪婪距离的思想， 本文提出了一种基于Count Sketch的预处理贪婪Kaczmarz方法。该方法先对超定线性系统进行Count Sketch变换， 然后对变换后的矩阵进行正交三角分解得到右预处理子， 再与贪婪Kaczmarz方法结合求解预处理后的方程组。我们分析了新方法的收敛性， 理论分析说明在系数矩阵条件数较大的情形下， 新方法的收缩因子比已有方法的收缩因子小。数值实验进一步验证了新方法的有效性。

1 基于Count Sketch的预处理贪婪Kaczmarz方法

本文采用$\mathit{R}\left(\mathit{A}\right)$表示由矩阵$\mathit{A}$列向量张成的值空间。其中，$\mathit{R}\left(\mathit{A}\right)=\{\mathit{A}\mathit{x}\left|\mathit{x}\in {\mathit{R}}^n\}\right.$。${\mathit{A}}^{\mathrm{T}}$， ${\mathit{A}}^{†},\Vert \mathit{A}{\Vert }_F$，${\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}(\mathit{A})$和${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}(\mathit{A})$分别表示矩阵$\mathit{A}$的转置， Moore-Penrose伪逆， Frobenius范数， 最小非零奇异值和最大奇异值； ${\mathit{x}}_{\star }={\mathit{A}}^{†}\mathit{b}$是线性代数方程组（1）的最小范数解。

Kaczmarz方法的收敛速率与行指标的选择策略密切相关， 选取与当前迭代点距离最大超平面的贪婪距离准则便是一种有效策略， Nutini等^［6］证明了在某些应用中贪婪选择策略和随机选择策略的计算成本是相当的， 并且在许多情况下贪婪策略比随机策略具有更快的收敛速度。

结合贪婪距离策略得到的贪婪Kaczmarz方法（见算法1）在文献［17］中以最大加权残差Kaczmarz的名称最早出现， 但其给出的收敛估计并不容易计算。后来， 文献［18］在贪婪随机Kaczmarz方法^［7］收敛性证明的基础上更新了贪婪Kaczmarz方法收敛速度的理论估计， 如引理1所示。

算法1  　贪婪Kaczmarz方法^［6］。①置$k:=\mathrm{0}$。②计算$i_k=\underset{i\in \{\mathrm{1},\cdots ,m\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}\frac{{\left|{\mathit{b}}_i-{\mathit{a}}_i^{\mathbf{T}}{\mathit{x}}_k\right|}^{\mathrm{2}}}{{‖{\mathit{a}}_i‖}_{\mathrm{2}}^{\mathrm{2}}}$。③计算${\mathit{x}}_{k+\mathrm{1}}={\mathit{x}}_k+\frac{{\mathit{b}}_{i_k}-{\mathit{a}}_{i_k}^{\mathbf{T}}{\mathit{x}}_k}{{‖{\mathit{a}}_{i_k}‖}_{\mathrm{2}}^{\mathrm{2}}}{\mathit{a}}_{i_k}$。④置$k=k+\mathrm{1}$， 转步骤②。

引理1　  ^［18］    若线性方程组$\mathit{A}\mathit{x}=\mathit{b}$相容， 初始估计${\mathit{x}}_{\mathrm{0}}\in \mathbf{R}({\mathit{A}}^{\mathrm{T}})$。则通过贪婪Kaczmarz方法生成的迭代序列$\{{\mathit{x}}_k{\}}_{k=\mathrm{0}}^{\mathrm{\infty }}$收敛到其最小范数解${\mathit{x}}_{\star }={\mathit{A}}^{†}\mathit{b}$， 且迭代序列$\{{\mathit{x}}_k{\}}_{k=\mathrm{0}}^{\mathrm{\infty }}$的误差满足

$${‖{\mathit{x}}_{\mathrm{1}}-{\mathit{x}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}\left(\mathit{A}\right)}{\Vert \mathit{A}{\Vert }_F^{\mathrm{2}}}\right){‖{\mathit{x}}_{\mathrm{0}}-{\mathit{x}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

对于$k=\mathrm{1,2},\dots$，

$${‖{\mathit{x}}_{k+\mathrm{1}}-{\mathit{x}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le$$

$$\left(\mathrm{1}-\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}\left(\mathit{A}\right)}{\underset{\mathrm{1}\le j\le m}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sum }_{i=\mathrm{1},i\ne j}^m\Vert {\mathit{a}}_i{\Vert }_{\mathrm{2}}^{\mathrm{2}}}\right){‖{\mathit{x}}_k-{\mathit{x}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

贪婪Kaczmarz方法的收敛速度很大程度上依赖于系数矩阵的条件数$k(\mathit{A})$， 当$k(\mathit{A})$较大时， 贪婪Kaczmarz方法的收敛速度会变慢。解决此问题的方法之一是使用预处理矩阵。由于本问题中$m>n$， 因此考虑如下形式的右预处理方法：

$$\mathit{A}\mathit{P}\mathit{y}=\mathit{b},\hspace{1em}\mathit{x}=\mathit{P}\mathit{y}$$

其中， $\mathit{P}\in {\mathbf{R}}^{n\times n}$， 且$k(\mathit{A}\mathit{P})<k(\mathit{A})$。

因此， 最理想的预处理子$\mathit{P}$应满足$k\left(\mathit{A}\mathit{P}\right)=\mathrm{1}$。受到文献［14］的启发， 满足此条件的预处理子可通过矩阵的正交三角分解获取：计算$\mathit{A}$的正交三角分解 $\mathit{A}=\mathit{Q}\mathit{R},\mathit{Q}\in {\mathbf{R}}^{m\times n},\mathit{R}\in {\mathbf{R}}^{n\times n}$， 然后取$\mathit{P}={\mathit{R}}^{†}$。使用这样的$\mathit{P}$进行预处理， 再用贪婪Kaczmarz方法求解， 可得到正交三角分解预处理的贪婪Kaczmarz方法， 该方法在迭代步数方面能够大大加快收敛速率。

注意到对矩阵$\mathit{A}$进行正交三角分解的成本为$O(m{n}^{\mathrm{2}})$， 这对于高度超定的系统（$m\gg n$）来说过高。因此文献［14］使用了基于行抽样的sketching方法以较低的成本获得上述预处理子$\mathit{P}$的近似， 这里的行抽样sketch矩阵$\mathit{S}$作用于系数矩阵$\mathit{A}$后， 相当于随机选择了$\mathit{A}$的某$d$行。但抽样的行数$d$是一个未知的超参数， 且文献［14］没有给出使用此类sketching预处理的Kaczmarz方法的收敛性证明。

Sketching方法中， 一种关键技术是Johnson-Lindenstrauss型嵌入算法^［19］， Count Sketch就是其中一种， 并且在某些情况下它比抽样（sampling）更有优势^［11］。

定义1　^［11-12］定义Count Sketch变换为$\mathit{S}=\mathit{\Phi }\mathit{D}\in {\mathbf{R}}^{d\times m}$。这里， $\mathit{D}$是一个$m\times m$的随机对角矩阵， 每个对角项以相同的概率独立被选为-1或1， $\mathit{\Phi }\in {\{\mathrm{0,1}\}}^{d\times m}$是一个$d\times m$的二元矩阵， 其中${\mathit{\Phi }}_{h(i),i}=\mathrm{1}$， 其余元素为0， 这里$h:[\mathbf{m}]\to [\mathbf{d}]$是一个随机映射， 对于每个$i\in [\mathbf{m}]$， $j\in [\mathbf{d}]$， $h(i)=j$的概率为$\mathrm{1}/d$。

通过定义1可以发现， Count Sketch矩阵由每列有且仅有一个1的矩阵$\mathit{\Phi }$和对角元为1或-1的对角矩阵$\mathit{D}$相乘得到， 它与系数矩阵$\mathit{A}$相乘的成本仅为$O(\mathrm{n}\mathrm{n}\mathrm{z}(\mathit{A}))$， 且$\mathit{\Phi }$和$\mathit{D}$均可以低成本的构造， 其中$\mathrm{n}\mathrm{n}\mathrm{z}(\mathit{A})$为矩阵$\mathit{A}$的非零项个数。

Count sketch变换还满足如下引理：

引理2　  ^［20］    若$\mathit{S}\in {\mathbf{R}}^{d\times m}$是一个Count Sketch变换， $d=O(\left(n^{\mathrm{2}}+n\right)/\left(\delta {\epsilon }^{\mathrm{2}}\right))$， 其中$O<\delta ,\epsilon <\mathrm{1}$， 则在$\mathrm{1}-\delta$的概率下， 有

$$(\mathrm{1}-\epsilon )\Vert \mathit{A}\mathit{x}{\Vert }_{\mathrm{2}}\le \Vert \mathit{S}\mathit{A}\mathit{x}{\Vert }_{\mathrm{2}}\le (\mathrm{1}+\epsilon )\Vert \mathit{A}\mathit{x}{\Vert }_{\mathrm{2}},\mathrm{对所}\mathrm{有}\mathit{x}\in {\mathbf{R}}^n$$

（3）

$$(\mathrm{1}-\epsilon ){\sigma }_i(\mathit{A})\le {\sigma }_i(\mathit{S}\mathit{A})\le (\mathrm{1}+\epsilon ){\sigma }_i(\mathit{A}),\mathrm{对所}\mathrm{有}\mathrm{1}\le i\le n$$

（4）

Count sketch矩阵的这些性质也引起了学者的关注， 例如， 文献［16］和文献［21］分别将Count Sketch与最大加权残差Kaczmarz和斜投影最大加权残差Kaczmarz方法相结合， 在计算时间方面大大提高了原方法的效率。

本文提出基于Count Sketch的预处理贪婪Kaczmarz方法， 如算法2所示。

算法2  　基于Count Sketch的预处理贪婪Kaczmarz方法。① 创建Count Sketch矩阵$\mathit{S}\in {\mathbf{R}}^{d\times m},n<d<m$， 并计算$\widehat{\mathit{A}}=\mathit{S}\mathit{A}\in {\mathbf{R}}^{d\times n}$。②对$\widehat{\mathit{A}}$进行正交三角分解， $\widehat{\mathit{A}}=\mathit{Q}\mathit{R}$， $\mathit{Q}\in {\mathbf{R}}^{d\times n}$， $\mathit{R}\in {\mathbf{R}}^{n\times n}$并将$\mathit{P}={\mathit{R}}^{†}$作为右预处理矩阵。③计算$\tilde{\mathit{A}}=\mathit{A}\mathit{P}$， 并记$\tilde{\mathit{A}}$的第$i$行为${\tilde{\mathit{a}}}_i,i=\mathrm{1,2},\cdots ,m$。④置$k:=\mathrm{0}$。⑤计算$i_k=\underset{i\in \{\mathrm{1},\cdots ,m\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}\frac{{\left|{\mathit{b}}_i-{\tilde{\mathit{a}}}_i^{\mathrm{T}}{\mathit{y}}_k\right|}^{\mathrm{2}}}{{‖{\tilde{\mathit{a}}}_i‖}_{\mathrm{2}}^{\mathrm{2}}}$。⑥计算${\mathit{y}}_{k+\mathrm{1}}={\mathit{y}}_k+\frac{{\mathit{b}}_{i_k}-{\tilde{\mathit{a}}}_{i_k}^{\mathrm{T}}{\mathit{y}}_k}{{‖{{\tilde{\mathit{a}}}_i}_k‖}_{\mathrm{2}}^{\mathrm{2}}}{{\tilde{\mathit{a}}}_i}_k$。⑦置$k=k+\mathrm{1}$， 转步骤⑤。⑧计算$\mathit{x}=\mathit{P}\mathit{y}$。

表1列出基于Count Sketch的预处理贪婪Kaczmarz（PCSGK）方法、贪婪Kaczmarz（GK）方法以及正交三角分解预处理的贪婪Kaczmarz（PGK）方法的计算成本。

表1  PCSGK、GK和PGK方法的计算成本

Tab.1  Computation costs of PCSGK, GK and PGK method

参数	PCSGK	GK	PGK
计算$S$与$A$的乘积	$O(\mathrm{n}\mathrm{n}\mathrm{z}(\mathit{A}))$	--	--
系数矩阵的正交三角分解	$O(d{n}^2)$	--	$O(m{n}^2)$
矩阵$\mathit{R}$求逆	$O(n^3)$	--	$O(n^3)$
对系数矩阵作用右预处理子	$O(m{n}^2)$	--	$O(m{n}^2)$
每步选取行指标	$O(mn)$	$O(mn)$	$O(mn)$
近似解更新	$O(n)$	$O(n)$	$O(n)$
迭代步数	$O(\kappa {(\mathit{A}\mathit{P})}^2$$·\mathrm{l}\mathrm{o}\mathrm{g}(1/\epsilon ))$	$O(\kappa {(\mathit{A})}^2$$·\mathrm{l}\mathrm{o}\mathrm{g}(1/\epsilon ))$	$O(\kappa {(\mathit{A}\mathit{P})}^2$$·\mathrm{l}\mathrm{o}\mathrm{g}(1/\epsilon ))$

从表1可以观察到， 使用sketching方法后， 正交三角分解的复杂度从$O\left(m{n}^{\mathrm{2}}\right)$降低到了$O\left(d{n}^{\mathrm{2}}\right)$， 如果$d=\gamma n$， 则可以重写为$O\left(n^{\mathrm{3}}\right)$， 其中$\gamma \ge \mathrm{1}$是某个给定的常数。这样， 计算预处理子$\mathit{P}$的总复杂度从$O\left(m{n}^{\mathrm{2}}\right)+O\left(n^{\mathrm{3}}\right)$降低为了为$O\left(n^{\mathrm{3}}\right)$， 与矩阵$\mathit{A}$的行数$m$无关。

同时， 注意到当系数矩阵缩放条件数$\kappa (\mathit{A})$较大时， 经典的贪婪Kaczmarz方法可能需要很多次迭代才能收敛， 如果此时用$\mathit{P}$进行预处理， 则有$\kappa (\mathit{A}\mathit{P})~\sqrt[]{n}$， 则算法可以在$O(n)$步内收敛。此时预处理的贪婪Kaczmarz方法在迭代步数上的改进可以远远大于预处理的成本。此外， 在很多问题上， 预处理只需要进行一次， 因此对于条件数较大的系统， 基于Count Sketch的预处理贪婪Kaczmarz方法是一个较为高效的选择。

2 收敛性分析

本节说明基于Count Sketch的预处理贪婪Kaczmarz方法的收敛性。定理1是本文的关键定理：

定理1  　若$\mathit{S}\in {\mathbf{R}}^{d\times m}$是一个Count Sketch变换， $d=O(\left(n^{\mathrm{2}}+n\right)/\left(\delta {\epsilon }^{\mathrm{2}}\right))$， $O<\delta ,\epsilon <\mathrm{1}$， 初始估计${\mathit{y}}_{\mathrm{0}}\in R({\mathit{P}}^{\mathrm{T}}{\mathit{A}}^{\mathrm{T}})$， 记基于Count Sketch的预处理贪婪Kaczmarz方法生成的迭代序列为$\{{\mathit{y}}_k{\}}_{k=\mathrm{0}}^{\mathrm{\infty }}$。 则在$\mathrm{1}-\delta$的概率下， 序列$\{{\mathit{y}}_k{\}}_{k=\mathrm{0}}^{\mathrm{\infty }}$收敛到$\tilde{\mathit{A}}\mathit{y}=\mathit{b}$的最小范数解${\mathit{y}}_{\star }={\tilde{\mathit{A}}}^{†}\mathit{b}$， 且满足

$${‖{\mathit{y}}_{\mathrm{1}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-\frac{(\mathrm{1}-{\epsilon }_{\mathrm{1}}{)}^{\mathrm{2}}}{(\mathrm{1}+{\epsilon }_{\mathrm{2}}{)}^{\mathrm{2}}}\frac{\mathrm{1}}{r}\right){‖{\mathit{y}}_{\mathrm{0}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

对于$k=\mathrm{1,2},\dots$，

$${‖{\mathit{y}}_{k+\mathrm{1}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-(\mathrm{1}-{\epsilon }_{\mathrm{1}}{)}^{\mathrm{2}}\frac{\mathrm{1}}{\underset{\mathrm{1}\le j\le m}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sum }_{i=\mathrm{1},i\ne j}^m\parallel {\tilde{a}}_i{\parallel }_{\mathrm{2}}^{\mathrm{2}}}\right){‖{\mathit{y}}_k-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

（5）

其中， ${\epsilon }_{\mathrm{1}}=\frac{\epsilon }{\mathrm{1}+\epsilon }$， ${\epsilon }_{\mathrm{2}}=\frac{\epsilon }{\mathrm{1}-\epsilon }$， $r$为矩阵$\mathit{A}$的秩， $\tilde{\mathit{A}}=\mathit{A}\mathit{P}$。

证明  　由于基于Count Sketch的预处理贪婪Kaczmarz方法本质上为用贪婪Kaczmarz方法求解线性方程组$\mathit{A}\mathit{P}\mathit{y}=\mathit{b}$， 即$\tilde{\mathit{A}}\mathit{y}=\mathit{b}$， 则由引理1， 有

$${‖{\mathit{y}}_{\mathrm{1}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\tilde{\mathit{A}})}{\Vert \tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}}\right){‖{\mathit{y}}_{\mathrm{0}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

（6）

对于$k=\mathrm{1,2},\dots$，

$${‖{\mathit{y}}_{k+\mathrm{1}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\tilde{\mathit{A}})}{\underset{\mathrm{1}\le j\le m}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sum }_{i=\mathrm{1},i\ne j}^m\parallel {\tilde{a}}_i{\parallel }_{\mathrm{2}}^{\mathrm{2}}}\right){‖{\mathit{y}}_k-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

（7）

接下来对矩阵$\tilde{\mathit{A}}$的性质进行研究。根据引理2中式（4）， $\mathrm{对所}\mathrm{有}\mathrm{1}\le i\le n$， 有

$$(\mathrm{1}-\epsilon ){\sigma }_i(\tilde{\mathit{A}})\le {\sigma }_i(\mathit{S}\tilde{\mathit{A}})\le (\mathrm{1}+\epsilon ){\sigma }_i(\tilde{\mathit{A}})$$

（8）

通过记${\epsilon }_{\mathrm{1}}=\frac{\epsilon }{\mathrm{1}+\epsilon }$， ${\epsilon }_{\mathrm{2}}=\frac{\epsilon }{\mathrm{1}-\epsilon }$， 可将式（8）化为

$$(\mathrm{1}-{\epsilon }_{\mathrm{1}}){\sigma }_i(\mathit{S}\tilde{\mathit{A}})\le {\sigma }_i(\tilde{\mathit{A}})\le (\mathrm{1}+{\epsilon }_{\mathrm{2}}){\sigma }_i(\mathit{S}\tilde{\mathit{A}})$$

（9）

而$\mathit{S}\tilde{\mathit{A}}=\mathit{Q}\mathit{R}{\mathit{R}}^{†}$，

$$n<d<m,\mathit{Q}\in {\mathbf{R}}^{d\times n}$$， $$\mathit{R}\in {\mathbf{R}}^{n\times n}$$，

故对所有$\mathrm{1}\le i\le r$， ${\sigma }_i(\mathit{S}\tilde{\mathit{A}})={\sigma }_i(\mathit{Q})=\mathrm{1}$。

所以

$${\sigma }_i(\tilde{\mathit{A}})\ge (\mathrm{1}-{\epsilon }_{\mathrm{1}}){\sigma }_i(\mathit{S}\tilde{\mathit{A}})=\mathrm{1}-{\epsilon }_{\mathrm{1}}$$

（10）

又因为

$$\Vert \mathit{S}\tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}={\displaystyle \sum _{i=\mathrm{1}}^r}{\sigma }_i^{\mathrm{2}}(\mathit{S}\tilde{\mathit{A}})=r$$

所以

$$\Vert \tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}\le (\mathrm{1}+{\epsilon }_{\mathrm{2}})\Vert \mathit{S}\tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}=(\mathrm{1}+{\epsilon }_{\mathrm{2}})r$$

（11）

将式（10）和式（11）代入式（6）和式（7）即可完成定理证明。

备注 根据定义1， 矩阵$\mathit{S}$的一行或几行元素可能为0， 这意味着$\mathit{S}\tilde{\mathit{A}}$可能存在零行， 则有

$$\underset{\mathrm{1}\le j\le m}{\mathrm{m}\mathrm{a}\mathrm{x}}{\sum }_{i=\mathrm{1},i\ne j}^m\Vert {\tilde{\mathit{a}}}_i{\Vert }_{\mathrm{2}}^{\mathrm{2}}=\Vert \tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}-\underset{\mathrm{1}\le j\le m}{\mathrm{m}\mathrm{i}\mathrm{n}}\Vert {\tilde{\mathit{a}}}_i{\Vert }_{\mathrm{2}}^{\mathrm{2}}$$

$$\le (\mathrm{1}+{\epsilon }_{\mathrm{2}})\Vert \mathit{S}\tilde{\mathit{A}}{\Vert }_F^{\mathrm{2}}+\mathrm{0}=(\mathrm{1}+{\epsilon }_{\mathrm{2}})r$$

因此， 将上式代入式（5）， 在$\mathrm{1}-\delta$的概率下， 有

$${‖{\mathit{y}}_{k+\mathrm{1}}-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}\le \left(\mathrm{1}-\frac{(\mathrm{1}-{\epsilon }_{\mathrm{1}}{)}^{\mathrm{2}}}{(\mathrm{1}+{\epsilon }_{\mathrm{2}}{)}^{\mathrm{2}}}\frac{\mathrm{1}}{r}\right){‖{\mathit{y}}_k-{\mathit{y}}_{\star }‖}_{\mathrm{2}}^{\mathrm{2}}$$

注意到， 对于引理1， 由于

$$\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}{\Vert \mathit{A}{\Vert }_F^{\mathrm{2}}}=\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}{{\displaystyle \sum _{i=\mathrm{1}}^r}{\sigma }_i^{\mathrm{2}}(\mathit{A})}\le \frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}{r{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}=\frac{\mathrm{1}}{r}$$

当系数矩阵病态时， $\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}{\Vert \mathit{A}{\Vert }_F^{\mathrm{2}}}$与$\frac{\mathrm{1}}{r}$相差较大， 若$\frac{\mathrm{1}-{\epsilon }_{\mathrm{1}}}{\mathrm{1}+{\epsilon }_{\mathrm{2}}}\to \mathrm{1}$， 则$\frac{{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}(\mathit{A})}{\Vert \mathit{A}{\Vert }_F^{\mathrm{2}}}$小于$\frac{\mathrm{1}-{\epsilon }_{\mathrm{1}}}{\mathrm{1}+{\epsilon }_{\mathrm{2}}}\frac{\mathrm{1}}{r}$， 基于Count Sketch的预处理贪婪Kaczmarz方法从而具有更快的收敛速度。

3 数值实验

本节通过数值实验比较贪婪Kaczmarz（GK）方法、Count Sketch最大加权残差Kaczmarz（CS-MWRK）方法、正交三角分解预处理的贪婪Kaczmarz（PGK）方法和基于Count Sketch的预处理贪婪Kaczmarz（PCSGK）方法， 数值结果显示新方法在迭代步数（简记为IT）与计算时间（简记为CPU）方面均具有较好的性能。

在实验中， 初始估计${\mathit{x}}_{\mathrm{0}}$取为零向量， 构造一个精确解$\mathit{x}$， 并且右端项$\mathit{b}=\mathit{A}\mathit{x}$。考虑到随机性，本节所有数值结果取20次独立实验的平均值。停机准则为相对残差（简记为RES）

$$\frac{\Vert \mathit{b}-\tilde{\mathit{A}}{\mathit{y}}_k{\Vert }_{\mathrm{2}}}{\Vert \mathit{b}{\Vert }_{\mathrm{2}}}<{\mathrm{10}}^{-\mathrm{3}}$$

或者迭代步数超过10万步。实验表格中， 若在10万步内未达到指定精度， 即用'- -'表示。对于正交三角分解预处理的贪婪Kaczmarz方法和基于Count Sketch的预处理贪婪Kaczmarz方法， 计算时间为预处理时间与迭代时间之和。

例1  　生成一个$m\times n$的随机高斯矩阵并进行奇异值分解， 然后将其奇异值替换为$j^{\alpha }$， $j=\mathrm{1,2},\dots n$， 以获取一个新的系数矩阵， 该矩阵的条件数即为$n^{\alpha }$。通过把$n$取为50， $\alpha$分别取2和2.5将系数矩阵条件数分别控制为2 500.00和17 677.67； 通过把$\gamma$分别取为5， 10和15来改变Sketching的行数$d=\gamma n$。

表2和表3列出了贪婪Kaczmarz方法、Count Sketch最大加权残差Kaczmarz方法、正交三角分解预处理的贪婪Kaczmarz方法以及基于Count Sketch的预处理贪婪Kaczmarz方法的迭代步数与计算时间。

表2  $\alpha =\mathrm{2}$时的数值结果

Tab.2  Numerical results for $\alpha =\mathrm{2}$

矩阵规模		5 000×50			10 000×50			50 000×50
算法		IT	CPU	RES	IT	CPU	RES	IT	CPU	RES
GK		26 075.00	1.927 1	9.99×10-4	14 282.00	3.910 1	9.97×10-4	9 010.00	14.641 7	9.78×10-4
CS-MWRK	$\mathrm{\gamma }=5$	75 153.80	0.584 5	9.90×10-4	42019.45	0.341 5	9.88×10-4	35 374.60	0.347 7	9.87×10-4
	$\mathrm{\gamma }=10$	46 532.65	0.435 0	9.91×10-4	33059.60	0.362 4	9.87×10-4	24 240.15	0.237 7	9.82×10-4
	$\mathrm{\gamma }=15$	43 037.60	0.510 7	9.89×10-4	29658.40	0.393 6	9.90×10-4	22 273.45	0.280 5	9.84×10-4
PGK		50.00	0.009 3	8.94×10-4	44.00	0.025 6	9.16×10-4	37.00	0.141 8	9.86×10-4
PCSGK	$\mathrm{\gamma }=5$	62.60	0.007 9	9.59×10-4	53.80	0.020 7	9.39×10-4	45.00	0.093 3	9.16×10-4
	$\mathrm{\gamma }=10$	55.15	0.008 1	9.41×10-4	49.05	0.018 1	8.97×10-4	40.05	0.080 3	9.08×10-4
	$\mathrm{\gamma }=15$	51.40	0.006 3	9.39×10-4	47.60	0.017 5	9.39×10-4	38.40	0.079 2	9.43×10-4

表3  $\alpha =\mathrm{2.5}$时的数值结果

Tab.3  Numerical results for $\alpha =\mathrm{2.5}$

矩阵规模		5 000×50			10 000×50			50 000×50
算法		IT	CPU	RES	IT	CPU	RES	IT	CPU	RES
GK		16 940.00	1.152 7	9.97×10-4	11 319.00	2.949 0	9.75×10-4	8 396.00	16.146 9	9.72×10-4
CS-MWRK	$\mathrm{\gamma }=5$	40 299.80	0.313 8	9.87×10-4	35 991.45	0.262 3	9.86×10-4	34 652.20	0.324 8	9.87×10-4
	$\mathrm{\gamma }=10$	28 866.85	0.274 8	9.83×10-4	26 393.20	0.250 9	9.89×10-4	26 023.40	0.317 9	9.87×10-4
	$\mathrm{\gamma }=15$	26 065.40	0.332 2	9.73×10-4	23 087.00	0.268 1	9.89×10-4	22 005.00	0.276 7	9.80×10-4
PGK		48.00	0.010 7	8.78×10-4	44.00	0.026 3	9.62×10-4	40.00	0.169 3	7.96×10-4
PCSGK	$\mathrm{\gamma }=5$	61.75	0.007 0	9.52×10-4	58.40	0.019 9	9.13×10-4	44.05	0.096 4	8.52×10-4
	$\mathrm{\gamma }=10$	54.60	0.006 8	9.43×10-4	48.15	0.018 3	8.94×10-4	40.40	0.085 3	9.19×10-4
	$\mathrm{\gamma }=15$	51.60	0.006 0	9.21×10-4	47.00	0.016 3	9.17×10-4	40.80	0.080 7	9.03×10-4

从表2和表3可以看出， PGK方法在迭代步数方面表现的最好， 这是因为PGK方法的预处理子在理论上可以将系统的条件数改善到最优值1， 而PCSGK方法中的Count Sketch过程在加速正交三角分解的过程中牺牲了一定的精度， 所以在迭代步数方面， PCSGK方法相对于PGK方法稍显逊色。

从表2和表3的数据中同样可以观察到， 通过选取适当的sketching行数， PCSGK方法的迭代步数可以十分接近PGK方法的迭代步数， 同时Count Sketch过程有效减少了预处理时间， 这使得PCSGK方法在计算时间方面明显优于PGK方法。综合计算时间和迭代步数两方面的因素， PCSGK在这4种方法中效率是最高的。此外， 在生成的随机矩阵的条件数不变， 即$\alpha$不变时， 随着sketching行数的增加， PCSGK方法在迭代步数的加速效果越明显， 这可以解释为随着$\gamma$的增大， 定理1中$\frac{\mathrm{1}-{\epsilon }_{\mathrm{1}}}{\mathrm{1}+{\epsilon }_{\mathrm{2}}}$这一项越来越接近于1， 从而加快了PCSGK方法的收敛。

图1绘制了矩阵规模为$\mathrm{10}\hspace{0.25em}\mathrm{000}\times \mathrm{50},\alpha =\mathrm{2.5},\gamma =\mathrm{15}$时， 4种方法的log（RES）分别随迭代步数和计算时间变化的曲线。从图1中可看出，4种方法的收敛曲线进一步验证了PCSGK方法的高效性。

图1  矩阵规模为$\mathrm{10}\hspace{0.25em}\mathrm{000}\times \mathrm{50},\mathit{\alpha }=\mathrm{2.5},\mathit{\gamma }=\mathrm{15}$时的收敛曲线

Fig.1  Convergence curves for matrix of $\mathrm{10}\hspace{0.25em}\mathrm{000}\times \mathrm{50}$, with$\mathit{\alpha }=\mathrm{2.5},\mathit{\gamma }=\mathrm{15}$

例2  　选取了佛罗里达稀疏矩阵库中的一些矩阵进行实验， 其中欠定矩阵取其转置。这里包含矩阵lp_fit2d、lp_d6cube、relat5和relat6。

表4列出了矩阵信息， 以及Sketching行数$d$选取一些适当的值时， 4种方法的数值结果。

表4  例2的矩阵信息和数值结果

Tab.4  Matrix information and numerical results for Example 2

矩阵名称		lp_fit2dT	lp_d6cube T	relat5	relat6
$m\times n$		10 524×25	6 184×415	340×35	2 340×157
稠密度/%		49.05	1.47	8.89	2.21
条件数		1.74×103	4.93×1019	Inf	Inf
$d$		300	800	70	300
GK	IT	192.00	7 130.00	196.00	1 359.00
	CPU	0.057 2	0.769 4	0.005 1	0.079 4
	RES	7.76×10-4	9.56×10-4	9.82×10-4	9.82×10-4
CS-MWRK	IT	26 723.15	53 185.70	334.65	3 544.70
	CPU	0.503 4	3.426 3	0.012 1	0.154 9
	RES	9.85×10-4	9.78×10-4	9.53×10-4	9.93×10-4
PGK	IT	31.00	521.00	56.00	283.00
	CPU	0.038 4	12.527 1	0.004 5	0.411 3
	RES	8.09×10-4	9.89×10-4	8.50×10-4	9.91×10-4
PCSGK	IT	35.45	1 631.30	113.20	763.10
	CPU	0.012 7	2.827 1	0.002 1	0.082 1
	RES	9.85×10-4	9.96×10-4	9.68×10-4	9.90×10-4

从表4可以观察到， PCSGK方法相对GK方法和CS-MWRK方法具有明显较少的迭代步数。

图2绘制了对于表4中的矩阵算例relat5， 4种方法的log（RES）分别随迭代步数和计算时间变化的曲线。

图2  矩阵relat5的收敛曲线

Fig.2  Convergence curves for matrix relat5

从图2可以观察到， 对于算例relat5这样的系统， PGK方法所需的预处理时间较长， 甚至出现了PCSGK方法已经收敛但PGK方法还未完成预处理的情形， PCSGK方法虽然牺牲了预处理子的精确度， 在计算时间方面大大提高了算法效率。因此综合迭代步数与计算时间这两方面的因素， 基于Count Sketch的预处理贪婪Kaczmarz方法是一个较好的选择。

4 结语

提出了一种基于Count Sketch的预处理贪婪Kaczmarz方法， 理论分析证明了新方法的收敛性， 且在系数矩阵条件数较大的情况下， 基于Count Sketch的预处理贪婪Kaczmarz方法的收敛因子小于经典的收敛因子， 即新方法具有较快的收敛速度。数值实验的结果表明， 提出的新方法在计算时间与迭代步数方面均有较大的优势。

基于Count Sketch的预处理贪婪Kaczmarz方法能够快速求解超定线性方程组（1）的最小范数解， 主要原因在于采用了正交三角分解进行预处理， 并采用贪婪策略选择工作行。因此如何选用成本更低的预处理以及更高效选取工作行的策略， 仍然是值得研究的问题。

作者贡献声明

叶雨欣：算法设计和算法研究执行，收敛性证明，数值实验和数据分析、论文初稿写作；

殷俊锋：研究构思，实验设计指导，数据分析，论文写作与修改。