分块矩阵的Kaczmarz迭代算法在求解欠定线性方程组中的应用
题目描述：
考虑一个欠定线性方程组 \(A\mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{m \times n}\)，\(m < n\)，并且系数矩阵 \(A\) 是行满秩的（即 \(\text{rank}(A) = m\)）。方程组有无穷多解。我们需要寻找其中范数最小的解（即最小范数解）：\(\min \|\mathbf{x}\|_2\)，满足 \(A\mathbf{x} = \mathbf{b}\)。假设矩阵 \(A\) 是大型稀疏矩阵，或者以分块形式存储（例如按行分块）。请设计一种基于分块矩阵的Kaczmarz迭代算法，逐步逼近这个最小范数解，并解释其原理、收敛性和计算步骤。
解题过程循序渐进讲解：
第一步：理解问题背景与经典Kaczmarz迭代
经典Kaczmarz迭代是用于求解线性方程组 \(A\mathbf{x} = \mathbf{b}\) 的一种逐行投影方法。设 \(a_i^T\) 是 \(A\) 的第 \(i\) 行，\(b_i\) 是 \(\mathbf{b}\) 的第 \(i\) 个分量。迭代格式为：
\[\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \frac{b_i - a_i^T \mathbf{x}^{(k)}}{\|a_i\|_2^2} a_i \]
每次迭代将当前解正交投影到超平面 \(a_i^T \mathbf{x} = b_i\) 上。在超定方程组（\(m > n\)）中，该算法收敛到最小二乘解；在方阵满秩时收敛到唯一解。但对于欠定方程组（\(m < n\)），经典Kaczmarz迭代会收敛到某个解，但不一定是最小范数解。
第二步：欠定方程组的最小范数解
对于欠定方程组 \(A\mathbf{x} = \mathbf{b}\)，行满秩时最小范数解为：
\[\mathbf{x}^* = A^T (A A^T)^{-1} \mathbf{b} \]
这个解位于 \(A\) 的行空间（即 \(\mathbb{R}^n\) 中由 \(A\) 的行张成的子空间）。经典Kaczmarz迭代如果从零向量 \(\mathbf{x}^{(0)} = \mathbf{0}\) 开始，会收敛到最小范数解吗？实际上，如果按固定顺序循环选择行，收敛的解不一定在行空间中。但如果我们每次随机选择行（随机Kaczmarz），在期望意义下会收敛到最小范数解吗？也不是直接的。我们需要修改算法，确保迭代始终保持在行空间中。
第三步：投影到行空间的技巧
关键观察：最小范数解 \(\mathbf{x}^*\) 位于 \(A\) 的行空间，即 \(\mathbf{x}^* = A^T \mathbf{y}\) 对某个 \(\mathbf{y} \in \mathbb{R}^m\) 成立。如果我们定义变量 \(\mathbf{y}\)，则原方程变为 \(A A^T \mathbf{y} = \mathbf{b}\)。这是一个 \(m \times m\) 的对称正定方程组（因为 \(A\) 行满秩，\(A A^T\) 正定）。对这个新方程组应用Kaczmarz迭代，然后通过 \(\mathbf{x} = A^T \mathbf{y}\) 恢复原变量，即可得到最小范数解。具体地：
\[A A^T \mathbf{y} = \mathbf{b} \quad \Rightarrow \quad \text{对第 } i \text{ 行：} a_i^T A^T \mathbf{y} = b_i \]
设 \(\tilde{a}_i^T = a_i^T A^T\)（这是 \(A A^T\) 的第 \(i\) 行），则Kaczmarz迭代更新为：
\[\mathbf{y}^{(k+1)} = \mathbf{y}^{(k)} + \frac{b_i - \tilde{a}_i^T \mathbf{y}^{(k)}}{\|\tilde{a}_i\|_2^2} \tilde{a}_i \]
然后 \(\mathbf{x}^{(k)} = A^T \mathbf{y}^{(k)}\)。但这样需要计算 \(A A^T\) 的行，当 \(m, n\) 很大时计算和存储开销大。
第四步：避免显式形成 \(A A^T\)
注意 \(\tilde{a}_i^T = a_i^T A^T\) 是一个行向量，与向量 \(\mathbf{y}\) 的内积为 \(a_i^T A^T \mathbf{y} = a_i^T \mathbf{x}\)，其中 \(\mathbf{x} = A^T \mathbf{y}\)。并且 \(\|\tilde{a}_i\|_2^2 = a_i^T A^T A a_i\)。因此迭代可以完全用 \(\mathbf{x}\) 表示。从当前 \(\mathbf{x}^{(k)} = A^T \mathbf{y}^{(k)}\) 出发，更新 \(\mathbf{y}\) 后，新解为：
\[\mathbf{x}^{(k+1)} = A^T \mathbf{y}^{(k+1)} = A^T \left( \mathbf{y}^{(k)} + \frac{b_i - a_i^T \mathbf{x}^{(k)}}{a_i^T A^T A a_i} (A a_i) \right) = \mathbf{x}^{(k)} + \frac{b_i - a_i^T \mathbf{x}^{(k)}}{a_i^T A^T A a_i} A^T A a_i \]
这里 \(A^T A a_i\) 是一个向量，可以通过两次矩阵-向量乘计算：先计算 \(w = A a_i\)（\(m\) 维），再计算 \(A^T w\)（\(n\) 维）。这样避免了显式形成 \(A A^T\)。
第五步：分块矩阵情境下的算法扩展
当矩阵 \(A\) 很大时，常按行分块存储，例如 \(A = [A_1; A_2; \dots; A_p]\)，其中 \(A_j \in \mathbb{R}^{m_j \times n}\)，\(\sum_j m_j = m\)。我们可以将分块视为“超行”，对每个块应用块Kaczmarz迭代。对于欠定方程组，类似地，我们希望迭代保持在行空间。将分块 \(A_j\) 视为一个整体，对应的右端项分块为 \(\mathbf{b}_j\)。定义当前解 \(\mathbf{x}^{(k)} = A^T \mathbf{y}^{(k)}\)。块Kaczmarz迭代将当前解投影到子空间 \(\{ \mathbf{x} : A_j \mathbf{x} = \mathbf{b}_j \}\) 上，这需要求解一个最小二乘问题：
\[\mathbf{x}^{(k+1)} = \arg\min_{\mathbf{x}} \|\mathbf{x} - \mathbf{x}^{(k)}\|_2^2 \quad \text{s.t.} \quad A_j \mathbf{x} = \mathbf{b}_j \]
其解析解为：
\[\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + A_j^T (A_j A_j^T)^{-1} (\mathbf{b}_j - A_j \mathbf{x}^{(k)}) \]
注意这里 \(A_j A_j^T\) 是一个 \(m_j \times m_j\) 的矩阵，通常很小，容易求逆。验证：如果从 \(\mathbf{x}^{(k)} = A^T \mathbf{y}^{(k)}\) 开始，更新后：
\[\mathbf{x}^{(k+1)} = A^T \mathbf{y}^{(k)} + A_j^T (A_j A_j^T)^{-1} (\mathbf{b}_j - A_j A^T \mathbf{y}^{(k)}) \]
但 \(A_j A^T\) 是 \(A_j\) 乘以 \(A^T\) 的一个子块，不一定是单位阵。为了保持解在行空间中，我们需要确保更新后仍可写成 \(A^T\) 乘以某个向量。观察更新项：\(A_j^T (A_j A_j^T)^{-1} \mathbf{r}_j\)，其中 \(\mathbf{r}_j = \mathbf{b}_j - A_j \mathbf{x}^{(k)}\)。这个项可以看作 \(A^T\) 乘以一个在 \(A_j\) 行空间中的向量吗？不一定。但我们可以从对偶变量 \(\mathbf{y}\) 的角度推导。
第六步：对偶变量更新公式（分块情形）
设当前 \(\mathbf{x}^{(k)} = A^T \mathbf{y}^{(k)}\)。在块Kaczmarz中，我们选择第 \(j\) 个块，解约束 \(A_j \mathbf{x} = \mathbf{b}_j\) 下的最小范数更新。这等价于求解对偶问题：找增量 \(\Delta \mathbf{y}_j \in \mathbb{R}^{m_j}\)（仅对应第 \(j\) 块非零），使得新解 \(\mathbf{x} = A^T (\mathbf{y}^{(k)} + \Delta \mathbf{y}_j)\) 满足 \(A_j \mathbf{x} = \mathbf{b}_j\)，即：
\[A_j A^T (\mathbf{y}^{(k)} + \Delta \mathbf{y}_j) = \mathbf{b}_j \]
记 \(M = A A^T\)，并分块 \(M_{jj} = A_j A_j^T\)，\(M_{ij} = A_i A_j^T\)。但因为我们只更新 \(\mathbf{y}_j\)（其他块不变），所以方程简化为：
\[A_j A_j^T \Delta \mathbf{y}_j = \mathbf{b}_j - A_j A^T \mathbf{y}^{(k)} = \mathbf{b}_j - A_j \mathbf{x}^{(k)} \]
于是：
\[\Delta \mathbf{y}_j = (A_j A_j^T)^{-1} (\mathbf{b}_j - A_j \mathbf{x}^{(k)}) \]
从而：
\[\mathbf{y}^{(k+1)} = \mathbf{y}^{(k)} + e_j \Delta \mathbf{y}_j \]
其中 \(e_j\) 表示在第 \(j\) 块位置放入增量。相应地：
\[\mathbf{x}^{(k+1)} = A^T \mathbf{y}^{(k+1)} = \mathbf{x}^{(k)} + A_j^T \Delta \mathbf{y}_j = \mathbf{x}^{(k)} + A_j^T (A_j A_j^T)^{-1} (\mathbf{b}_j - A_j \mathbf{x}^{(k)}) \]
这正是上面块投影公式！而且由于 \(\mathbf{x}^{(k+1)} = A^T \mathbf{y}^{(k+1)}\)，它自动保持在行空间中。因此，分块Kaczmarz迭代（块投影更新）在欠定情况下，如果从 \(\mathbf{x}^{(0)} = \mathbf{0}\) 开始，会收敛到最小范数解。
第七步：算法步骤总结
- 初始化：\(\mathbf{x}^{(0)} = \mathbf{0} \in \mathbb{R}^n\)（或任意满足 \(\mathbf{x}^{(0)} = A^T \mathbf{y}^{(0)}\) 的向量，例如 \(\mathbf{y}^{(0)} = \mathbf{0}\)）。
- 对迭代次数 \(k = 0, 1, 2, \dots\)：
a. 选择分块索引 \(j_k\)（可按循环顺序或随机选取）。
b. 计算残差向量：\(\mathbf{r} = \mathbf{b}_{j_k} - A_{j_k} \mathbf{x}^{(k)}\)。
c. 求解小规模线性系统：\((A_{j_k} A_{j_k}^T) \mathbf{z} = \mathbf{r}\)（由于 \(A_{j_k}\) 行满秩，该系数矩阵对称正定，可用Cholesky分解或直接法求解）。
d. 更新解：\(\mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + A_{j_k}^T \mathbf{z}\)。 - 迭代直到残差 \(\|\mathbf{b} - A\mathbf{x}^{(k)}\|_2\) 小于给定容差。
第八步：收敛性分析
- 该算法是线性收敛的，收敛速度依赖于矩阵 \(A A^T\) 的条件数以及分块选取策略。
- 如果每次随机选取分块（随机块Kaczmarz），期望意义下收敛速度有明确界限。
- 由于迭代点始终在行空间中，且更新是到仿射子空间 \(\{ \mathbf{x} : A_j \mathbf{x} = \mathbf{b}_j \}\) 的投影，该算法等价于对正定系统 \((A A^T) \mathbf{y} = \mathbf{b}\) 应用块坐标下降法（或随机坐标下降法）。因此可以利用凸优化理论证明其收敛到最小范数解。
第九步：计算考虑
- 需要预计算每个分块的 \(A_j A_j^T\) 及其Cholesky分解，以加速步骤2c的求解。如果分块很小，这开销不大。
- 适用于大规模稀疏矩阵，因为矩阵-向量乘 \(A_{j_k} \mathbf{x}^{(k)}\) 和 \(A_{j_k}^T \mathbf{z}\) 可以利用稀疏性快速计算。
- 内存友好：只需存储分块 \(A_j\) 及其小矩阵 \(A_j A_j^T\)，无需整个 \(A\) 同时载入内存。
第十步：总结
分块Kaczmarz迭代通过将行投影推广为块投影，并结合对偶变量技巧，使得在欠定情形下能收敛到最小范数解。算法保持了Kaczmarz的简单性和适合大规模数据的特性，同时通过分块减少了迭代次数（因为每步处理多个方程），并确保解在行空间中，从而得到最小范数解。