好的,让我们来探讨一个在之前的长列表中尚未出现的重要算法。
随机化块坐标下降法(Randomized Block Coordinate Descent, RBCD)求解线性方程组
题目描述
我们考虑求解一个大型的线性方程组:
\[A\mathbf{x} = \mathbf{b} \]
其中 \(A \in \mathbb{R}^{n \times n}\) 是一个对称正定 矩阵,\(\mathbf{b} \in \mathbb{R}^n\) 是已知向量,\(\mathbf{x} \in \mathbb{R}^n\) 是未知向量。
当矩阵 \(A\) 的维度 \(n\) 非常大(例如数百万甚至更高),并且 \(A\) 可能是稠密或具有某种特殊结构时,传统的直接解法(如Cholesky分解)会因为计算复杂度和存储开销过高而变得不可行。
随机化块坐标下降法 是一种迭代优化算法。它将求解线性方程组 \(A\mathbf{x} = \mathbf{b}\) 转化为一个等价的无约束二次凸优化问题,然后通过随机选择一个坐标块(变量子集),在该块上精确最小化目标函数来更新解。这种方法特别适合并行化,并且在某些问题结构下,其收敛速度优于经典的确定性坐标下降法。
问题转化:从线性方程组到优化问题
对于对称正定矩阵 \(A\),求解 \(A\mathbf{x} = \mathbf{b}\) 等价于最小化如下二次凸函数:
\[f(\mathbf{x}) = \frac{1}{2} \mathbf{x}^T A \mathbf{x} - \mathbf{b}^T \mathbf{x} \]
为什么等价?
因为 \(f(\mathbf{x})\) 的梯度是 \(\nabla f(\mathbf{x}) = A\mathbf{x} - \mathbf{b}\)。梯度为零的点 \(\nabla f(\mathbf{x}) = 0\) 正好是 \(A\mathbf{x} = \mathbf{b}\) 的解。又由于 \(A\) 正定,\(f(\mathbf{x})\) 是强凸函数,其驻点是唯一的全局最小值点。
因此,我们的目标就是找到 \(\mathbf{x}^*\),使得:
\[\mathbf{x}^* = \arg\min_{\mathbf{x} \in \mathbb{R}^n} f(\mathbf{x}) \]
算法核心思想与步骤详解
算法的核心是 “分而治之” 和 “随机抽样”。我们不是一次更新所有的变量,而是每次只更新一个精心挑选的变量子集(一个“块”)。
步骤1:问题分区(Block Partitioning)
我们将变量索引集合 \(\{1, 2, ..., n\}\) 划分为 \(M\) 个互不相交的块(blocks):
\[\{1, 2, ..., n\} = \mathcal{B}_1 \cup \mathcal{B}_2 \cup ... \cup \mathcal{B}_M \]
每个块 \(\mathcal{B}_i\) 包含一组索引,其大小为 \(n_i = |\mathcal{B}_i|\)。例如,如果我们有12个变量,可以均匀地划分为3个块,每块4个变量:\(\mathcal{B}_1 = \{1,2,3,4\}, \mathcal{B}_2 = \{5,6,7,8\}, \mathcal{B}_3 = \{9,10,11,12\}\)。
分块可以是任意的,但通常根据问题的自然结构(如网格中的一片区域、图中的一个簇)或计算效率(使子问题易于求解)来划分。
步骤2:随机块选择与子问题定义
在第 \(k\) 次迭代,我们随机均匀地 从 \(\{1, 2, ..., M\}\) 中选取一个块索引 \(i_k\)。设当前解估计为 \(\mathbf{x}_k\)。
我们只允许属于块 \(\mathcal{B}_{i_k}\) 的变量改变,而固定其他所有块中的变量。定义坐标方向子空间:
我们只考虑更新向量 \(\mathbf{d}\),其仅在块 \(\mathcal{B}_{i_k}\) 对应的位置有非零值。记 \(U_{i_k}\) 为一个 \(n \times n_{i_k}\) 的矩阵,其列是单位矩阵中对应块 \(\mathcal{B}_{i_k}\) 索引的那些列。那么任何这样的更新 \(\mathbf{d}\) 都可以写成 \(\mathbf{d} = U_{i_k} \mathbf{y}\),其中 \(\mathbf{y} \in \mathbb{R}^{n_{i_k}}\) 是块内的新值与原值的差。
于是,子优化问题 是在这个低维子空间上最小化 \(f\):
\[\min_{\mathbf{y} \in \mathbb{R}^{n_{i_k}}} f(\mathbf{x}_k + U_{i_k} \mathbf{y}) \]
步骤3:精确求解子问题(块内最小化)
将 \(f(\mathbf{x}_k + U_{i_k} \mathbf{y})\) 展开:
\[f(\mathbf{x}_k + U_{i_k} \mathbf{y}) = \frac{1}{2} (\mathbf{x}_k + U_{i_k}\mathbf{y})^T A (\mathbf{x}_k + U_{i_k}\mathbf{y}) - \mathbf{b}^T (\mathbf{x}_k + U_{i_k}\mathbf{y}) \]
这是关于 \(\mathbf{y}\) 的二次函数。为了找到最小值,我们对其求梯度并令其为零:
\[\frac{\partial f}{\partial \mathbf{y}} = U_{i_k}^T (A(\mathbf{x}_k + U_{i_k}\mathbf{y}) - \mathbf{b}) = 0 \]
这等价于:
\[(U_{i_k}^T A U_{i_k}) \mathbf{y} = U_{i_k}^T (\mathbf{b} - A\mathbf{x}_k) \]
关键观察:
- \(U_{i_k}^T A U_{i_k}\) 是矩阵 \(A\) 的一个 \(n_{i_k} \times n_{i_k}\) 主子矩阵,即由 \(A\) 中行和列都属于块 \(\mathcal{B}_{i_k}\) 的元素构成的子矩阵。记作 \(A_{[i_k]}\)。由于 \(A\) 正定,它的任意主子矩阵 \(A_{[i_k]}\) 也是正定的,因此可逆。
- \(\mathbf{r}_k = \mathbf{b} - A\mathbf{x}_k\) 是当前迭代的残差向量。
- \(U_{i_k}^T \mathbf{r}_k\) 是残差向量在块 \(\mathcal{B}_{i_k}\) 对应位置上的分量,记作 \(\mathbf{r}_k^{[i_k]}\)。
因此,子问题的解为:
\[\mathbf{y}_k^* = (A_{[i_k]})^{-1} \mathbf{r}_k^{[i_k]} \]
步骤4:更新当前解
得到最优的块内更新 \(\mathbf{y}_k^*\) 后,我们更新解向量:
\[\mathbf{x}_{k+1} = \mathbf{x}_k + U_{i_k} \mathbf{y}_k^* \]
即,只更新块 \(\mathcal{B}_{i_k}\) 中的变量:对于索引 \(j \in \mathcal{B}_{i_k}\),\(x_{k+1}^{(j)} = x_k^{(j)} + y_k^{*(j)}\);对于其他索引,值保持不变。
步骤5:高效更新残差(可选但关键)
为了在下一步中不用重新计算整个 \(A\mathbf{x}_{k+1}\) 来得到新残差,我们可以利用更新是局部的这一事实:
\[\mathbf{r}_{k+1} = \mathbf{b} - A\mathbf{x}_{k+1} = \mathbf{b} - A(\mathbf{x}_k + U_{i_k}\mathbf{y}_k^*) = \mathbf{r}_k - A U_{i_k} \mathbf{y}_k^* \]
而 \(A U_{i_k} \mathbf{y}_k^*\) 只需要计算矩阵 \(A\) 中与块 \(\mathcal{B}_{i_k}\) 相关的那些列(即列索引在 \(\mathcal{B}_{i_k}\) 中的列)与向量 \(\mathbf{y}_k^*\) 的乘积。这比全矩阵-向量乘法便宜得多。
算法伪代码总结
输入:对称正定矩阵 A,向量 b,分块方案 {B_1, ..., B_M},初始解 x_0,最大迭代次数 K
输出:近似解 x_K
r_0 = b - A * x_0 // 初始化残差
for k = 0 to K-1:
// 1. 随机均匀选择一个块
i = 随机从 {1, 2, ..., M} 中选取
// 2. 提取主子矩阵和块残差
A_block = A[B_i, B_i] // 索引为 B_i 的行和列构成的子矩阵
r_block = r_k[B_i] // 残差在 B_i 索引上的分量
// 3. 精确求解块子问题
y_star = 求解线性方程组 A_block * y = r_block
// (通常使用针对小块矩阵的Cholesky分解或直接求逆,因为块很小)
// 4. 更新解
x_{k+1}[B_i] = x_k[B_i] + y_star
// x_{k+1} 在非 B_i 索引上的值等于 x_k 的值
// 5. 高效更新残差
r_{k+1} = r_k - A[:, B_i] * y_star // 仅需计算 A 中相关列与 y_star 的乘积
end for
核心优势与特性
- 低每次迭代成本:每次迭代只需求解一个维度为 \(n_i\)(块大小)的小型线性方程组,并更新残差。如果 \(A\) 是稀疏的,\(A[:, B_i]\) 与向量的乘法成本也很低。
- 随机性带来的收敛加速:与循环遍历所有块的确定性坐标下降法相比,随机均匀选择块具有理论上的优势。其期望收敛速度通常为 \(O(1/k)\)(对于强凸光滑函数),并且对问题条件数的依赖可能更小。
- 并行与分布式潜力:不同块的更新在计算上是独立的。虽然基本算法是串行的,但其框架易于扩展到异步并行或分布式设置,多个工作节点可以同时处理不同的块。
- 易于处理特殊结构:如果矩阵 \(A\) 具有分块对角、分块带状或某种可分离结构,那么主子矩阵 \(A_{[i]}\) 可能具有更简单的形式(甚至是对角阵),使得步骤3的求解极其快速。
典型应用场景
- 大规模机器学习:如线性支持向量机(SVM)、逻辑回归的优化问题,其目标函数常包含一个大规模数据矩阵的二次项。
- 求解由有限元法或有限差分法离散化偏微分方程产生的大型稀疏线性系统,特别是当系统矩阵具有天然的分块结构(如来自区域分解)时。
- 大规模最小二乘问题。
与已讲题目的区别
您提供的列表中虽然出现了“随机化块坐标下降法在凸优化问题求解中的应用”,但本讲解专注于其在线性方程组 \(A\mathbf{x}=\mathbf{b}\)(即二次函数最小化)这一特定、基础且重要的应用背景下的详细步骤、原理推导和计算技巧,阐明了如何从方程转化为优化问题,以及如何利用矩阵的分块结构实现高效迭代。这是理解更一般凸优化应用的基础。