随机化块Kaczmarz算法在稀疏线性方程组中的概率分布与收敛加速策略
我们先看算法背景。
随机化Kaczmarz算法是一种求解大规模线性方程组 \(A x = b\) 的迭代方法，尤其适用于超定或欠定稀疏矩阵。
当矩阵是分块结构时，采用随机化块Kaczmarz能显著加速收敛。今天介绍其概率分布设计与收敛加速策略。
1. 问题描述
设线性方程组
\[A x = b, \]
其中 \(A \in \mathbb{R}^{m \times n}\)，\(b \in \mathbb{R}^m\)，且 \(A\) 是稀疏矩阵。
我们将矩阵的行分成 \(t\) 个块（子矩阵）：
\[A = \begin{bmatrix} A_1 \\ A_2 \\ \vdots \\ A_t \end{bmatrix}, \quad b = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_t \end{bmatrix}, \]
每个块 \(A_i \in \mathbb{R}^{m_i \times n}\)，\(\sum_{i=1}^t m_i = m\)。
我们的目标是设计一种随机选取块迭代的规则，使得收敛比传统Kaczmarz更快。
2. 经典Kaczmarz与随机化Kaczmarz回顾
Kaczmarz迭代逐行更新：
\[x_{k+1} = x_k + \frac{b_i - a_i^T x_k}{\|a_i\|^2} a_i, \]
其中 \(a_i\) 是 \(A\) 的第 \(i\) 行。
随机化版本（Strohmer & Vershynin, 2009）按概率
\[p_i = \frac{\|a_i\|^2}{\|A\|_F^2} \]
选取行，其中 \(\|A\|_F\) 是Frobenius范数，这样能保证指数收敛期望。
3. 分块Kaczmarz算法
若选取块 \(i\)，则更新公式变为：
\[x_{k+1} = x_k + A_i^\dagger (b_i - A_i x_k), \]
其中 \(A_i^\dagger\) 是 \(A_i\) 的伪逆（若 \(A_i\) 行满秩，则 \(A_i^\dagger = A_i^T (A_i A_i^T)^{-1}\)）。
这相当于将当前残差在块 \(i\) 对应的行空间上做最小二乘投影。
4. 随机化块选择的概率设计
为了加速收敛，需要考虑块行范数与条件数两个因素。
方法1：依块Frobenius范数采样
定义
\[p_i = \frac{\|A_i\|_F^2}{\sum_{j=1}^t \|A_j\|_F^2}. \]
这是直接推广，但未考虑不同块的“困难程度”。
方法2：依块条件数加权的概率
记 \(\sigma_{\min}(A_i)\) 为 \(A_i\) 的最小奇异值，若行满秩，则收敛速度与 \(\sigma_{\min}^2(A_i)\) 有关。
定义
\[p_i \propto \|A_i\|_F^2 \cdot \frac{1}{\kappa_i^2}, \]
其中 \(\kappa_i = \frac{\sigma_{\max}(A_i)}{\sigma_{\min}(A_i)}\) 是 \(A_i\) 的条件数。
这样更倾向于选择条件数较小的块，避免迭代在病态块上进展缓慢。
5. 自适应权重调整策略
我们可以在迭代过程中动态调整概率：
- 开始时所有块等概率。
- 记录每次选择块 \(i\) 后残差下降量 \(\Delta r_i = \|b_i - A_i x_k\|^2\)。
- 若某块残差下降显著，增大其概率；否则降低。
具体：设 \(w_i^{(k)}\) 为块 \(i\) 在第 \(k\) 轮权重，归一化得概率
\[p_i^{(k)} = \frac{w_i^{(k)}}{\sum_j w_j^{(k)}}. \]
权重更新：
\[w_i^{(k+1)} = w_i^{(k)} \cdot \exp\left( \eta \cdot \frac{\Delta r_i}{\text{avg}(\Delta r)} \right), \]
\(\eta\) 是小正数，\(\text{avg}(\Delta r)\) 是最近几次下降量的滑动平均。
这样算法会自适应地更频繁选择当前残差大的块，类似“重点攻关”。
6. 收敛加速策略
除了自适应概率，还可结合以下技巧：
(1) 块内使用随机化行迭代
若块 \(A_i\) 较大，计算 \(A_i^\dagger\) 开销大，可以改为在块内做若干步随机化Kaczmarz行迭代，近似实现块投影。
(2) 过度松弛参数
更新公式引入松弛因子 \(\omega \in (0, 2)\)：
\[x_{k+1} = x_k + \omega A_i^\dagger (b_i - A_i x_k). \]
通常 \(\omega \approx 1.5\) 可加速，但需在实验中调节。
(3) 动量加速
引入动量项（类似随机梯度下降的加速）：
\[d_k = A_{i_k}^\dagger (b_{i_k} - A_{i_k} x_k), \]
\[ v_{k+1} = \beta v_k + d_k, \]
\[ x_{k+1} = x_k + \omega v_{k+1}. \]
参数 \(\beta\) 可取 0.9 左右。
7. 理论保证（简述）
随机化块Kaczmarz的收敛速度依赖于块延伸条件数（block extension of scaled condition number）：
\[\kappa_{\text{block}}^2 = \|A\|_F^2 \cdot \max_i \sigma_{\min}^{-2}(A_i), \]
收敛速率与 \(\kappa_{\text{block}}^2\) 成反比。
自适应概率可降低有效条件数，从而提升速度。
8. 算法伪代码（自适应权重版本）
输入：分块矩阵 A_1,...,A_t，右端项 b_1,...,b_t，初始解 x0，迭代次数 K，初始权重 w_i=1
输出：近似解 x
x = x0
for k = 1 to K do
计算概率 p_i = w_i / sum(w)
依概率 p_i 选取块索引 i
计算残差 r_i = b_i - A_i x
更新方向 d = A_i^\dagger r_i （若块大，可用若干步随机行迭代近似）
x = x + ω * d
记录 Δr = ||r_i||^2
更新平均残差下降量 avg_Δr（滑动平均）
更新 w_i = w_i * exp(η * Δr / avg_Δr)
end for
9. 数值示例说明
假设 \(A\) 是 1000×200 稀疏矩阵，分 20 块，每块 50 行。
若某块的条件数是其他块的 10 倍，等概率采样会使其拖慢收敛。
自适应权重后，算法会减少选取该块的频率，而更多选择“好”的块，使平均每步残差下降更大，从而加快整体收敛。
10. 实际应用场景
这种改进的随机化块Kaczmarz可用于：
- 大规模CT图像重建（投影矩阵分块对应不同角度）
- 分布式求解线性系统（每块在不同计算节点）
- 在线学习中的随机坐标下降法（数据点分块）
通过结合概率分布的自适应调整与动量加速，随机化块Kaczmarz在稀疏线性方程组中的收敛速度可以得到明显提升，同时保持计算每步的低成本。