随机化Kaczmarz算法在求解欠定线性方程组中的应用 题目描述 考虑一个欠定线性方程组 \( A\mathbf{x} = \mathbf{b} \)，其中 \( A \in \mathbb{R}^{m \times n} \)，\( m < n \)，并且假设方程组是相容的（即存在解）。由于方程个数少于未知数个数，解不唯一，通常希望找到具有最小范数的解（例如最小2-范数解 \( \mathbf{x}^* = A^\dagger \mathbf{b} \)，其中 \( A^\dagger \) 是 \( A \) 的Moore-Penrose伪逆）。随机化Kaczmarz算法是一种迭代方法，适用于求解大规模稀疏线性方程组。传统的Kaczmarz算法按行顺序投影，收敛速度受行顺序影响。随机化Kaczmarz算法通过引入随机行选择策略，可显著提高收敛速度。在欠定情况下，算法需稍作调整，以保证收敛到最小范数解。本题目将详细介绍该算法的原理、实现步骤及收敛性分析。 解题过程 背景与问题形式化 给定欠定系统 \( A\mathbf{x} = \mathbf{b} \)，\( A \in \mathbb{R}^{m \times n} \)（\( m < n \)），秩为 \( m \)（行满秩）。解集是一个仿射空间：\( \{ \mathbf{x} : A\mathbf{x} = \mathbf{b} \} = \mathbf{x}_ 0 + \text{null}(A) \)，其中 \( \mathbf{x}_ 0 \) 是一个特解，\( \text{null}(A) \) 是 \( A \) 的零空间。我们希望找到最小2-范数解 \( \mathbf{x}^* = A^\dagger \mathbf{b} \)，它满足 \( \|\mathbf{x}^* \| 2 = \min {A\mathbf{x} = \mathbf{b}} \|\mathbf{x}\|_ 2 \)。在几何上，这等价于在解空间中寻找最靠近原点的点。 经典Kaczmarz算法回顾 对于超定或方阵系统，经典Kaczmarz算法迭代更新如下：从任意初始猜测 \( \mathbf{x} 0 \) 开始，在第 \( k \) 步，选择一行索引 \( i_ k \)，将当前迭代点 \( \mathbf{x} k \) 投影到该行对应的超平面 \( \langle \mathbf{a} {i_ k}, \mathbf{x} \rangle = b {i_ k} \) 上： \[ \mathbf{x} {k+1} = \mathbf{x} k + \frac{b {i_ k} - \langle \mathbf{a} {i_ k}, \mathbf{x} k \rangle}{\|\mathbf{a} {i_ k}\| 2^2} \mathbf{a} {i_ k} \] 其中 \( \mathbf{a}_ {i_ k}^\top \) 是 \( A \) 的第 \( i_ k \) 行。在超定情况下，若按顺序循环选择行，算法收敛到解（若相容）或最小二乘解（若不相容）。但收敛速度依赖于行之间的角度。 随机化策略引入 随机化Kaczmarz算法以概率分布选取行索引，通常按行范数平方作为权重：设 \( p_ i = \|\mathbf{a}_ i\|_ 2^2 / \|A\|_ F^2 \)，其中 \( \|A\|_ F \) 是Frobenius范数。每一步独立地从分布 \( \{p_ i\} \) 中选取 \( i_ k \)。这种策略已被证明在期望意义下指数收敛，且收敛速率与条件数相关。但对于欠定系统，直接应用上述公式会收敛到解空间中的某个点，但不一定是最小范数解。 欠定系统的调整：对偶视角 考虑对偶问题：最小范数解 \( \mathbf{x}^* = A^\top \mathbf{y}^* \)，其中 \( \mathbf{y}^* \in \mathbb{R}^m \) 满足 \( AA^\top \mathbf{y}^* = \mathbf{b} \)。这是因为 \( A^\dagger = A^\top (AA^\top)^{-1} \) 当 \( A \) 行满秩。因此，可先求解对偶系统 \( AA^\top \mathbf{y} = \mathbf{b} \)，这是一个 \( m \times m \) 的对称正定系统（因为 \( AA^\top \) 正定）。然后通过 \( \mathbf{x} = A^\top \mathbf{y} \) 得到原系统的最小范数解。 随机化Kaczmarz应用于对偶系统 将对偶变量 \( \mathbf{y} \) 初始化为 \( \mathbf{y}_ 0 = \mathbf{0} \)。定义对偶系统的矩阵为 \( C = AA^\top \in \mathbb{R}^{m \times m} \)，其第 \( i \) 行为 \( \mathbf{c}_ i^\top = \mathbf{a}_ i^\top A \)。但更有效的方式是直接在原变量空间操作。注意到对偶系统的第 \( i \) 个方程为 \( \langle \mathbf{a}_ i, A^\top \mathbf{y} \rangle = b_ i \)，即 \( \langle A\mathbf{a}_ i, \mathbf{y} \rangle = b_ i \)？等等，这里需要谨慎推导。实际上： \[ (AA^\top \mathbf{y})_ i = \mathbf{a}_ i^\top (A^\top \mathbf{y}) = \langle \mathbf{a} i, A^\top \mathbf{y} \rangle = b_ i \] 因此，对偶系统的第 \( i \) 个方程是 \( \langle A^\top \mathbf{y}, \mathbf{a} i \rangle = b_ i \)。 将对偶系统的Kaczmarz更新应用于 \( \mathbf{y} \)： \[ \mathbf{y} {k+1} = \mathbf{y} k + \frac{b {i_ k} - \langle A^\top \mathbf{y} k, \mathbf{a} {i_ k} \rangle}{\|A^\top \mathbf{a} {i_ k}\| 2^2} (A^\top \mathbf{a} {i_ k})? \quad \text{不，这是错误的。} \] 正确形式：对偶系统的系数矩阵为 \( AA^\top \)，其第 \( i \) 行是 \( \mathbf{a} i^\top A \)。但计算 \( \| \mathbf{a} i^\top A \| 2^2 \) 代价高。更聪明的方法是保持原变量 \( \mathbf{x} = A^\top \mathbf{y} \)。于是，在第 \( k \) 步，选择行索引 \( i_ k \) 后，对偶变量更新为： \[ \mathbf{y} {k+1} = \mathbf{y} k + \frac{b {i_ k} - \langle \mathbf{a} {i_ k}, \mathbf{x} k \rangle}{\|\mathbf{a} {i_ k}\| 2^2} \mathbf{e} {i_ k} \] 其中 \( \mathbf{e} {i_ k} \) 是第 \( i_ k \) 个标准单位向量。这是因为对偶系统的第 \( i \) 个方程是 \( (AA^\top \mathbf{y})_ i = b_ i \)，而 \( (AA^\top \mathbf{y})_ i = \mathbf{a}_ i^\top (A^\top \mathbf{y}) = \mathbf{a}_ i^\top \mathbf{x} \)。因此残差为 \( b_ i - \mathbf{a}_ i^\top \mathbf{x} \)，且 \( \| \mathbf{a}_ i^\top A \|_ 2^2 \) 实际上等于 \( \mathbf{a}_ i^\top A A^\top \mathbf{a}_ i = \|\mathbf{a}_ i\|_ 2^2 \)？不，注意这里 \( \mathbf{a}_ i \) 是行向量（或列向量），维度混乱。我们重新整理： 令 \( \mathbf{a}_ i^\top \) 表示 \( A \) 的第 \( i \) 行（行向量）。 则对偶系统 \( AA^\top \mathbf{y} = \mathbf{b} \) 的第 \( i \) 个方程为 \( \mathbf{a}_ i^\top A^\top \mathbf{y} = b_ i \)。 但 \( \mathbf{a}_ i^\top A^\top = (A \mathbf{a}_ i)^\top \)，其中 \( A \mathbf{a}_ i \) 是一个 \( m \) 维向量。 该方程可写为 \( \langle A \mathbf{a}_ i, \mathbf{y} \rangle = b_ i \)。 对偶系统的Kaczmarz更新应为： \[ \mathbf{y} {k+1} = \mathbf{y} k + \frac{b {i_ k} - \langle A \mathbf{a} {i_ k}, \mathbf{y} k \rangle}{\|A \mathbf{a} {i_ k}\| 2^2} (A \mathbf{a} {i_ k}) \] 但 \( \|A \mathbf{a} {i_ k}\| 2^2 = \mathbf{a} {i_ k}^\top A^\top A \mathbf{a} {i_ k} \)，计算仍然昂贵。 这里的关键观察是：我们不需要显式维护 \( \mathbf{y} \)。注意 \( \mathbf{x} k = A^\top \mathbf{y} k \)。从 \( \mathbf{y} \) 的更新公式，两边左乘 \( A^\top \)： \[ \mathbf{x} {k+1} = A^\top \mathbf{y} {k+1} = A^\top \mathbf{y} k + \frac{b {i_ k} - \langle A \mathbf{a} {i_ k}, \mathbf{y} k \rangle}{\|A \mathbf{a} {i_ k}\| 2^2} A^\top (A \mathbf{a} {i_ k}) \] 而 \( A^\top \mathbf{y} k = \mathbf{x} k \)，并且 \( \langle A \mathbf{a} {i_ k}, \mathbf{y} k \rangle = \mathbf{a} {i_ k}^\top A^\top \mathbf{y} k = \mathbf{a} {i_ k}^\top \mathbf{x} k \)。同时 \( A^\top (A \mathbf{a} {i_ k}) = (A^\top A) \mathbf{a} {i_ k} \)。因此更新变为： \[ \mathbf{x} {k+1} = \mathbf{x} k + \frac{b {i_ k} - \langle \mathbf{a} {i_ k}, \mathbf{x} k \rangle}{\|A \mathbf{a} {i_ k}\| 2^2} (A^\top A) \mathbf{a} {i_ k} \] 但计算 \( (A^\top A) \mathbf{a} {i_ k} \) 和 \( \|A \mathbf{a}_ {i_ k}\|_ 2^2 \) 仍然代价高。实际上，在欠定情况下，我们通常不直接使用对偶系统的Kaczmarz，而是采用另一种等价形式。 改进的随机化Kaczmarz（RK）用于欠定系统 更常用的方法是：注意到最小范数解位于 \( \text{row}(A) \)（\( A \) 的行张成的空间）。因此，如果初始猜测 \( \mathbf{x} 0 \in \text{row}(A) \)，那么每次Kaczmarz迭代（投影到超平面）会保持迭代点仍在 \( \text{row}(A) \) 中，因为投影方向是行向量 \( \mathbf{a} {i_ k} \in \text{row}(A) \)。因此，从 \( \mathbf{x} 0 = \mathbf{0} \) 开始（它在行空间中），算法将收敛到解空间中在行空间上的分量，而这正是最小范数解。 因此，算法与标准RK相同： \[ \mathbf{x} {k+1} = \mathbf{x} k + \frac{b {i_ k} - \langle \mathbf{a}_ {i_ k}, \mathbf{x} k \rangle}{\|\mathbf{a} {i_ k}\| 2^2} \mathbf{a} {i_ k} \] 其中行索引 \( i_ k \) 以概率 \( p_ i = \|\mathbf{a}_ i\|_ 2^2 / \|A\|_ F^2 \) 随机选取。从 \( \mathbf{x}_ 0 = \mathbf{0} \) 开始，迭代点始终在行空间中，最终收敛到最小范数解 \( \mathbf{x}^* \)。 收敛性分析 可以证明，对于欠定相容系统，上述随机化Kaczmarz算法在期望意义下线性收敛到最小范数解。定义误差向量 \( \mathbf{e} k = \mathbf{x} k - \mathbf{x}^* \)。由于 \( \mathbf{x}^* \) 是解，满足 \( A\mathbf{x}^* = \mathbf{b} \)。在一步更新后： \[ \mathbf{e} {k+1} = \mathbf{e} k - \frac{\langle \mathbf{a} {i_ k}, \mathbf{e} k \rangle}{\|\mathbf{a} {i_ k}\| 2^2} \mathbf{a} {i_ k} \] 两边取范数平方： \[ \|\mathbf{e} {k+1}\|_ 2^2 = \|\mathbf{e}_ k\| 2^2 - \frac{|\langle \mathbf{a} {i_ k}, \mathbf{e} k \rangle|^2}{\|\mathbf{a} {i_ k}\|_ 2^2} \] 对随机索引 \( i_ k \) 取条件期望（给定 \( \mathbf{e} k \)）： \[ \mathbb{E}[ \|\mathbf{e} {k+1}\|_ 2^2 \mid \mathbf{e}_ k] = \|\mathbf{e}_ k\| 2^2 - \sum {i=1}^m p_ i \frac{|\langle \mathbf{a}_ i, \mathbf{e}_ k \rangle|^2}{\|\mathbf{a}_ i\|_ 2^2} \] 代入 \( p_ i = \|\mathbf{a}_ i\|_ 2^2 / \|A\| F^2 \)： \[ \mathbb{E}[ \|\mathbf{e} {k+1}\|_ 2^2 \mid \mathbf{e}_ k] = \|\mathbf{e}_ k\|_ 2^2 - \frac{1}{\|A\| F^2} \sum {i=1}^m |\langle \mathbf{a}_ i, \mathbf{e}_ k \rangle|^2 = \|\mathbf{e}_ k\|_ 2^2 - \frac{\|A \mathbf{e}_ k\|_ 2^2}{\|A\|_ F^2} \] 由于 \( \mathbf{e}_ k \) 在行空间中（因为 \( \mathbf{x}_ 0, \mathbf{x}^* \in \text{row}(A) \)，且更新方向也在行空间），所以 \( A \mathbf{e}_ k \neq \mathbf{0} \) 除非 \( \mathbf{e} k = \mathbf{0} \)。设 \( \sigma {\min}^2(A) \) 是 \( A \) 的最小非零奇异值的平方（即 \( A \) 作为行空间到 \( \mathbb{R}^m \) 的映射的最小奇异值），则有 \( \|A \mathbf{e}_ k\| 2^2 \geq \sigma {\min}^2(A) \|\mathbf{e}_ k\| 2^2 \)。因此： \[ \mathbb{E}[ \|\mathbf{e} {k+1}\|_ 2^2 \mid \mathbf{e} k] \leq \left(1 - \frac{\sigma {\min}^2(A)}{\|A\|_ F^2}\right) \|\mathbf{e}_ k\|_ 2^2 \] 定义 \( \kappa_ F(A) = \|A\| F / \sigma {\min}(A) \)（一种条件数），则收敛因子为 \( 1 - 1/\kappa_ F^2(A) \)。这表明算法线性收敛，且收敛速度依赖于 \( \kappa_ F(A) \)。 算法步骤总结 输入：\( A \in \mathbb{R}^{m \times n} \)（行满秩，\( m < n \)），\( \mathbf{b} \in \mathbb{R}^m \)，最大迭代次数 \( K \) 或容差 \( \epsilon \)。 初始化：\( \mathbf{x}_ 0 = \mathbf{0} \in \mathbb{R}^n \)。 预计算行范数平方 \( \|\mathbf{a}_ i\|_ 2^2 \) 和 \( \|A\| F^2 = \sum {i=1}^m \|\mathbf{a}_ i\|_ 2^2 \)。 对 \( k = 0, 1, \dots, K-1 \)： 以概率 \( p_ i = \|\mathbf{a}_ i\|_ 2^2 / \|A\|_ F^2 \) 随机选取行索引 \( i_ k \)。 计算残差 \( r_ k = b_ {i_ k} - \langle \mathbf{a}_ {i_ k}, \mathbf{x}_ k \rangle \)。 更新：\( \mathbf{x}_ {k+1} = \mathbf{x} k + (r_ k / \|\mathbf{a} {i_ k}\| 2^2) \, \mathbf{a} {i_ k} \)。 输出：\( \mathbf{x}_ K \) 作为最小范数解的近似。 扩展与变体 该算法可推广到不相容系统（即无解）的情况，此时算法收敛到最小二乘解。也可引入松弛参数 \( \omega \)： \[ \mathbf{x} {k+1} = \mathbf{x} k + \omega \frac{b {i_ k} - \langle \mathbf{a} {i_ k}, \mathbf{x} k \rangle}{\|\mathbf{a} {i_ k}\| 2^2} \mathbf{a} {i_ k} \] 适当选择 \( \omega \in (0,2) \) 可加速收敛。此外，对于大规模问题，可使用块随机化Kaczmarz，每次选取一组行进行投影。 数值实验示例（思路） 生成一个随机行满秩矩阵 \( A \in \mathbb{R}^{m \times n} \)（\( m < n \)），计算其最小范数解 \( \mathbf{x}^* = A^\top (AA^\top)^{-1} \mathbf{b} \)。运行上述算法，观察误差 \( \|\mathbf{x}_ k - \mathbf{x}^* \|_ 2 \) 的衰减，验证以速率 \( (1 - 1/\kappa_ F^2(A))^k \) 指数下降。 通过以上步骤，随机化Kaczmarz算法能够高效求解大规模欠定线性方程组，并收敛到最小范数解。该方法特别适用于矩阵 \( A \) 稀疏且行数较少的情况，因为每次迭代仅需访问一行数据，计算成本低，内存需求小。