这一步大大减少了内积运算的数量，从 $O(jn)$ 降至 $O(kn)$。
c. **规范化**：计算 $h_{j+1, j} = \|\mathbf{w}\|_2$，如果 $h_{j+1, j}$ 接近于0，则提前终止。否则，令 $\mathbf{v}_{j+1} = \mathbf{w} / h_{j+1, j}$。

经过这个过程，我们得到了一组近似正交的基向量 $V_m = [\mathbf{v}_1, \mathbf{v}_2, ..., \mathbf{v}_m] \in \mathbb{R}^{n \times m}$，以及一个上Hessenberg矩阵 $H_m \in \mathbb{R}^{(m+1)\times m}$，其元素 $h_{i,j}$ 记录了正交化过程中的系数。它们近似满足关系：

等式之所以是近似的，是因为我们采用了不完全正交化。

其中 $U_m, W_m \in \mathbb{R}^{m \times m}$ 是正交矩阵，$\Sigma_m = \text{diag}(\sigma_1, \sigma_2, ..., \sigma_m)$，且 $\sigma_1 \ge \sigma_2 \ge ... \ge \sigma_m \ge 0$。

或者直接指定 $r$。然后，我们保留前 $r$ 个奇异值和对应的左右奇异向量：

其中 $U_m^{(r)}, W_m^{(r)} \in \mathbb{R}^{m \times r}$ 由 $U_m, W_m$ 的前 $r$ 列组成，$\Sigma_m^{(r)} = \text{diag}(\sigma_1, ..., \sigma_r) \in \mathbb{R}^{r \times r}$。

这一步捕捉了 $T_m$（即 $A$ 在近似Krylov子空间上投影）最主要的谱信息，并过滤掉了可能由近似正交化引入的数值噪声或不重要信息，增强了数值稳定性。

即

由于 $U_m^{(r)}$ 是列正交的，我们可以用最小二乘法求解这个方程。令 $\mathbf{z} = (W_m^{(r)})^T \mathbf{y}_m \in \mathbb{R}^r$，则方程变为：

其最小二乘解为：

这里 $(\Sigma_m^{(r)})^{-1} = \text{diag}(1/\sigma_1, ..., 1/\sigma_r)$。由于我们进行了截断，忽略的小奇异值对应的分量被舍弃，这起到了正则化的效果，对病态问题有益。

基于Krylov子空间的FABLAS算法解大规模稀疏线性方程组 题目描述 FABLAS（Fast Approximate Block Lanczos with Approximate SVD）算法是一种结合了Krylov子空间方法与低秩近似技巧的迭代算法，专门用于求解大规模稀疏线性方程组 \(A\mathbf{x} = \mathbf{b}\)。与标准Krylov子空间方法（如GMRES、CG）不同，FABLAS算法在构建Krylov子空间的过程中，采用了近似技术来加速正交化过程，并利用近似的奇异值分解（SVD）来有效捕捉系数矩阵 \(A\) 的主导谱信息，从而在保证一定精度的前提下，显著降低计算复杂度和内存开销，特别适用于那些矩阵维数巨大、但有效数值秩相对较低的问题。 解题过程 求解问题为 \(A\mathbf{x} = \mathbf{b}\)，其中 \(A \in \mathbb{R}^{n \times n}\) 是大规模稀疏矩阵。FABLAS算法的核心思想是：在一个低维的Krylov子空间中，通过构建一个近似于 \(A\) 在该子空间上投影的、规模小得多的稠密矩阵，并求解这个约化后的小规模问题来得到原问题的近似解。其核心加速技巧在于对Krylov基向量的正交化过程进行近似，以及对投影矩阵进行低秩近似分解。 步骤一：生成近似Krylov子空间 标准的Krylov子空间方法（如Arnoldi过程）通过完全的Gram-Schmidt正交化来构建一组标准正交基，这需要 \(O(m^2n)\) 的运算量（\(m\) 是子空间维数）。FABLAS算法采用“近似正交化”来加速。 初始向量 ：选择初始猜测解 \(\mathbf{x}_ 0\)（常取零向量或根据问题设定），计算初始残差 \(\mathbf{r}_ 0 = \mathbf{b} - A\mathbf{x}_ 0\)。令第一个基向量 \(\mathbf{v}_ 1 = \mathbf{r}_ 0 / \|\mathbf{r}_ 0\|_ 2\)。 递推生成近似正交基 ： 对于 \(j = 1, 2, ..., m\)（\(m \ll n\) 是设定的子空间最大维数）： a. 矩阵向量乘 ：计算 \(\mathbf{w} = A\mathbf{v} j\)。 b. 不完全正交化 ：不进行完全的Gram-Schmidt正交化（即对所有前 \(j\) 个基向量做正交化）。相反，FABLAS算法采用一种“部分”或“选择性”的正交化策略。一种常见做法是只对最近的 \(k\) 个基向量（\(k\) 是一个小的常数，如5或10）执行Gram-Schmidt正交化，即： \[ \mathbf{w} := \mathbf{w} - \sum {i=\max(1, j-k+1)}^{j} (\mathbf{w}^T\mathbf{v} i) \mathbf{v} i \] 这一步大大减少了内积运算的数量，从 \(O(jn)\) 降至 \(O(kn)\)。 c. 规范化 ：计算 \(h {j+1, j} = \|\mathbf{w}\| 2\)，如果 \(h {j+1, j}\) 接近于0，则提前终止。否则，令 \(\mathbf{v} {j+1} = \mathbf{w} / h_ {j+1, j}\)。 经过这个过程，我们得到了一组近似正交的基向量 \(V_ m = [ \mathbf{v}_ 1, \mathbf{v} 2, ..., \mathbf{v} m] \in \mathbb{R}^{n \times m}\)，以及一个上Hessenberg矩阵 \(H_ m \in \mathbb{R}^{(m+1)\times m}\)，其元素 \(h {i,j}\) 记录了正交化过程中的系数。它们近似满足关系： \[ A V_ m \approx V {m+1} H_ m \] 等式之所以是近似的，是因为我们采用了不完全正交化。 步骤二：构建近似投影并利用低秩近似 在标准Krylov方法中，我们寻找近似解 \(\mathbf{x}_ m = \mathbf{x}_ 0 + V_ m \mathbf{y}_ m\)，其中 \(\mathbf{y}_ m \in \mathbb{R}^m\) 通过求解投影方程 \((V_ m^T A V_ m) \mathbf{y}_ m = V_ m^T \mathbf{r}_ 0\) 得到。由于基向量是近似正交的，投影矩阵 \(V_ m^T A V_ m\) 并非精确的 \(H_ m\) 的前 \(m\) 行。FABLAS算法的关键创新在于巧妙地处理这个投影。 构造稠密投影矩阵 ：虽然基向量近似正交，但我们可以直接计算一个更可靠的稠密投影矩阵 \(T_ m = V_ m^T A V_ m \in \mathbb{R}^{m \times m}\)。注意，由于 \(V_ m\) 的列是近似正交的，且 \(A\) 是稀疏的，计算 \(T_ m = V_ m^T (A V_ m)\) 的代价是 \(O(m^2n)\)，但由于 \(m\) 很小，这个代价是可接受的。 对 \(T_ m\) 进行低秩SVD近似 ： FABLAS算法的另一个核心思想是假设 \(A\) 在其主导的Krylov子空间上的作用可以用一个低秩算子很好地近似。因此，我们对小稠密矩阵 \(T_ m\) 进行经济型奇异值分解（SVD）： \[ T_ m = U_ m \Sigma_ m W_ m^T \] 其中 \(U_ m, W_ m \in \mathbb{R}^{m \times m}\) 是正交矩阵，\(\Sigma_ m = \text{diag}(\sigma_ 1, \sigma_ 2, ..., \sigma_ m)\)，且 \(\sigma_ 1 \ge \sigma_ 2 \ge ... \ge \sigma_ m \ge 0\)。 截断SVD（低秩近似） ： 根据一个设定的容差 \(tol\)（例如 \(10^{-6}\)）或保留前 \(r\) 个最大的奇异值（\(r < m\)），我们对 \(T_ m\) 进行低秩近似。找到最小的 \(r\) 使得： \[ \sum_ {i=1}^{r} \sigma_ i^2 \ge (1 - tol^2) \sum_ {i=1}^{m} \sigma_ i^2 \] 或者直接指定 \(r\)。然后，我们保留前 \(r\) 个奇异值和对应的左右奇异向量： \[ T_ m \approx U_ m^{(r)} \Sigma_ m^{(r)} (W_ m^{(r)})^T \] 其中 \(U_ m^{(r)}, W_ m^{(r)} \in \mathbb{R}^{m \times r}\) 由 \(U_ m, W_ m\) 的前 \(r\) 列组成，\(\Sigma_ m^{(r)} = \text{diag}(\sigma_ 1, ..., \sigma_ r) \in \mathbb{R}^{r \times r}\)。 这一步捕捉了 \(T_ m\)（即 \(A\) 在近似Krylov子空间上投影）最主要的谱信息，并过滤掉了可能由近似正交化引入的数值噪声或不重要信息，增强了数值稳定性。 步骤三：求解约化的小规模最小二乘问题 现在，我们在近似的Krylov子空间中寻找解。我们希望最小化残差范数： \[ \| \mathbf{b} - A(\mathbf{x}_ 0 + V_ m \mathbf{y}_ m) \|_ 2 = \| \mathbf{r}_ 0 - A V_ m \mathbf{y} m \| 2 \] 由于 \(A V_ m \approx V {m+1} H_ m\) 且基向量近似正交，最小化问题可以近似为： \[ \min {\mathbf{y}_ m} \| \beta \mathbf{e}_ 1 - H_ m \mathbf{y}_ m \|_ 2, \quad \text{其中 } \beta = \|\mathbf{r}_ 0\|_ 2, \mathbf{e}_ 1 = [ 1, 0, ..., 0 ]^T \in \mathbb{R}^{m+1} \] 然而，FABLAS算法利用我们构造的稠密矩阵 \(T_ m\) 和其低秩近似，构建了一个更稳健的模型。 构建投影方程 ：利用关系 \(A V_ m \approx V_ {m+1} H_ m\) 和基向量的近似正交性，一个更精确的逼近是考虑： \[ V_ m^T A V_ m \mathbf{y}_ m = V_ m^T \mathbf{r}_ 0 \] 即 \[ T_ m \mathbf{y}_ m = \mathbf{g}_ m，\quad 其中 \mathbf{g}_ m = V_ m^T \mathbf{r}_ 0 = \beta V_ m^T \mathbf{v}_ 1 = \beta \mathbf{e}_ 1^{(m)}，\mathbf{e}_ 1^{(m)} = [ 1, 0,...,0 ]^T \in \mathbb{R}^m \] 利用低秩SVD近似求解 ：将 \(T_ m\) 的低秩近似代入方程： \[ U_ m^{(r)} \Sigma_ m^{(r)} (W_ m^{(r)})^T \mathbf{y}_ m \approx \mathbf{g}_ m \] 由于 \(U_ m^{(r)}\) 是列正交的，我们可以用最小二乘法求解这个方程。令 \(\mathbf{z} = (W_ m^{(r)})^T \mathbf{y}_ m \in \mathbb{R}^r\)，则方程变为： \[ U_ m^{(r)} (\Sigma_ m^{(r)} \mathbf{z}) \approx \mathbf{g}_ m \] 其最小二乘解为： \[ \mathbf{z} = (\Sigma_ m^{(r)})^{-1} (U_ m^{(r)})^T \mathbf{g}_ m \] 这里 \((\Sigma_ m^{(r)})^{-1} = \text{diag}(1/\sigma_ 1, ..., 1/\sigma_ r)\)。由于我们进行了截断，忽略的小奇异值对应的分量被舍弃，这起到了正则化的效果，对病态问题有益。 恢复 \(\mathbf{y}_ m\) ： \[ \mathbf{y}_ m = W_ m^{(r)} \mathbf{z} = W_ m^{(r)} (\Sigma_ m^{(r)})^{-1} (U_ m^{(r)})^T \mathbf{g}_ m \] 步骤四：更新近似解并检查收敛 计算近似解 ： \[ \mathbf{x}_ m = \mathbf{x}_ 0 + V_ m \mathbf{y}_ m \] 计算真实残差 （或通过递归关系估计）：计算 \(\mathbf{r}_ m = \mathbf{b} - A\mathbf{x}_ m\)。 检查收敛 ：如果 \(\|\mathbf{r}_ m\|_ 2 / \|\mathbf{b}\|_ 2 < \epsilon\)（预设容差），则停止迭代，输出 \(\mathbf{x}_ m\) 作为近似解。 若不收敛，重新启动（Restart） ：如果未达到收敛且子空间维数达到最大值 \(m\)，则进行重启。令新的初始猜测 \(\mathbf{x}_ 0 := \mathbf{x}_ m\)，然后返回步骤一，基于新的残差重新开始构建Krylov子空间。重启策略是控制内存使用的关键，因为基向量矩阵 \(V_ m\) 的存储与 \(m\) 成正比。 总结 FABLAS算法的核心优势在于： 近似正交化 ：通过只对最近少数基向量进行正交化，大幅降低了构建Krylov子空间的计算成本。 低秩SVD滤波 ：对稠密投影矩阵 \(T_ m\) 进行低秩SVD近似，不仅有效提取了矩阵在子空间上作用的主导信息，还抑制了由近似正交化可能引入的误差，并隐式地提供了正则化，增强了对病态问题的鲁棒性。 灵活的重启策略 ：允许在固定的小内存开销下处理大规模问题。 该算法适用于那些系数矩阵 \(A\) 虽然是稀疏的，但其在相关子空间上的行为可以用一个低秩算子有效近似的问题，常见于某些偏微分方程离散化、数据科学和机器学习中的大规模线性系统。