在 AI 模型训练的优化器迭代、二阶方法求解或量子模拟中，经常需要计算矩阵函数作用于向量的结果 x = f(A)b，其中 A 是大规模稀疏 Hermite 矩阵，f 可能是逆运算（线性系统求解）、矩阵指数（时间演化）或符号函数（投影算子）。直接计算 f(A) 会导致稠密矩阵存储爆炸，而 Lanczos 算法通过 Krylov 子空间投影将问题降维至 k 维三对角矩阵 Tk，理论上可用 O(k) 次矩阵向量积完成计算。但标准实现需要存储完整的 n×k 基矩阵 Vk：对于 50 万变量、1000 次迭代的问题，仅基矩阵就需要 4GB 内存，这在多任务并行训练或边缘设备部署时成为瓶颈。
本文基于 Luca Lombardo 的工程实践，解析如何通过两遍（two-pass）Lanczos 算法将内存需求从 O(nk) 降至 O(n)，并在 Rust 中利用缓存局部性实现反直觉的性能提升——尽管计算量翻倍，实测运行时间仅增加 20%。
内存墙的根源：基矩阵存储困境
标准 Lanczos 算法通过 Arnoldi 正交化构建 Krylov 子空间 K_k(A, b) = span{b, Ab, A²b, ..., A^(k-1)b} 的标准正交基 V_k = [v₁, v₂, ..., v_k]。当 A 是 Hermite 矩阵时，投影矩阵 H_k 退化为实对称三对角矩阵 T_k，正交化简化为三项递推：
βⱼvⱼ₊₁ = Avⱼ - αⱼvⱼ - βⱼ₋₁vⱼ₋₁
其中 αⱼ = vⱼᴴ(Avⱼ) 是对角元，βⱼ = ||ṽⱼ₊₁||₂ 是次对角元。完成 k 次迭代后，解在降维空间中表示为 x_k = V_k y_k，其中 y_k = f(T_k)e₁||b||₂ 通过求解 k×k 小矩阵获得。
问题在于：y_k 的计算依赖完整的 T_k 矩阵，只能在所有迭代结束后进行。这迫使实现必须存储全部 k 个基向量，导致内存占用 O(nk)。对于典型场景（n=500k, k=1000, 双精度浮点），基矩阵需要 500,000 × 1,000 × 8 字节 ≈ 3.7GB，超过 L3 缓存容量数百倍。
两遍法的核心设计：时间换空间的缓存友好策略
两遍 Lanczos 算法打破时序依赖，将计算分为两个独立阶段：
第一遍：系数收集（Pass One）
运行标准 Lanczos 迭代，但仅保存标量系数 {α₁, β₁, ..., α_k, β_k}，每次迭代后立即丢弃基向量 vⱼ。内存中始终只保留三个 n 维向量：vⱼ₋₁、vⱼ 和工作向量 w。系数数组占用 O(k) 空间（对于 k=1000 仅 16KB），可忽略不计。
关键优化点：
- 零分配热路径：使用 Rust 的
faer 库提供的 MemStack 预分配栈空间，避免迭代中的动态内存分配
- 向量循环复用：通过
core::mem::swap 交换指针（x86-64 上编译为三条 mov 指令），而非拷贝数据
- 提前终止检测：当 βⱼ ≤ breakdown_tolerance（通常为 machine epsilon 的平方根）时停止，表明 Krylov 子空间已不变
第一遍结束后，求解小规模三对角系统 T_k y' = e₁ 得到系数向量 y_k = ||b||₂ · y'。
第二遍：解的重建（Pass Two）
使用第一遍保存的 {αⱼ, βⱼ} 重新运行 Lanczos 递推，但不再计算内积和范数（这些值已知），仅执行正交化和归一化。每生成一个基向量 vⱼ，立即累加其对解的贡献：
x_k += y_k[j] * v_j;
累加完成后丢弃 vⱼ，继续生成下一个向量。整个过程的工作集仍为三个 n 维向量（vⱼ₋₁、vⱼ、x_k），完全驻留在 L1 缓存中。
工程化参数清单：从理论到落地的量化指标
内存布局参数
| 参数 |
标准实现 |
两遍实现 |
说明 |
| 基向量存储 |
n×k×8 字节 |
3n×8 字节 |
两遍法仅保留滑动窗口 |
| 系数数组 |
2k×8 字节 |
2k×8 字节 |
两者相同 |
| 总内存（n=500k, k=1000） |
~3.8 GB |
~11.4 MB |
降低 330 倍 |
计算成本参数
- 矩阵向量积次数：标准法 k 次，两遍法 2k 次
- 内积计算：标准法 k 次（第一遍），两遍法 k 次（仅第一遍）
- 范数计算：标准法 k 次，两遍法 k 次（仅第一遍）
- 解重建：标准法 1 次稠密矩阵向量积（n×k），两遍法 k 次缓存内累加
性能边界阈值
基于 KKT 系统（Karush-Kuhn-Tucker 矩阵）的实测数据：
- 缓存优势显现点：n > 150,000 时，两遍法的缓存局部性开始抵消额外计算
- 最优场景：稀疏矩阵（非零元 < 0.1%）+ 大规模（n > 500k）+ 中等迭代（k = 500-1000）
- 退化场景：稠密矩阵时，矩阵向量积成为瓶颈，两遍法严格慢 2 倍
- 数值稳定性：浮点累加顺序差异导致 O(ε_machine) 级误差，对收敛性无影响
实测数据（n=500k, k=500, 稀疏度 0.01%）：
- 标准实现：3.8 GB 内存，12.3 秒
- 两遍实现：11.4 MB 内存，14.7 秒（仅慢 19.5%）
Rust 实现的关键技术点
1. SIMD 自动向量化
使用 faer 的 zip! 宏融合正交化操作：
zip!(w.rb_mut(), v_prev).for_each(|unzip!(w_i, v_prev_i)| {
*w_i = sub(w_i, &mul(&beta_prev_scaled, v_prev_i));
});
编译器将其优化为 AVX2/AVX-512 指令，单次处理 4-8 个双精度浮点数。
2. 矩阵无关算子抽象
通过 LinOp trait 实现算子多态，支持稀疏矩阵、矩阵自由算子（如有限差分）和 GPU 加速：
pub trait LinOp<T: ComplexField> {
fn apply(&self, y: MatMut<'_, T>, x: MatRef<'_, T>,
par: Par, stack: &mut MemStack);
}
3. 数值稳定性保障
- Breakdown 检测：当 βⱼ < sqrt(ε_machine) 时提前终止，避免除零
- 重正交化触发器：可选的 Paige 准则检测正交性丢失（未在两遍法中实现，因重建保证确定性）
适用场景与限制
推荐使用两遍法的场景：
- 大规模稀疏线性系统求解（共轭梯度法的预处理器）
- 矩阵指数计算（量子模拟、时间演化）
- 内存受限环境（边缘设备、多任务并行训练）
- 迭代次数 k 占问题规模 n 的比例 < 1%
不推荐的场景：
- 稠密矩阵（计算成本严格翻倍）
- 极小规模问题（n < 10k，缓存优势不明显）
- 需要中间基向量的应用（如特征值估计）
资料来源
本文技术细节源自 Luca Lombardo 的工程实现与技术报告：
该实现已在 Rust 生态中集成 faer 线性代数库，可直接用于生产环境的矩阵函数计算任务。对于 AI 训练中的二阶优化器（如 L-BFGS、自然梯度）或物理模拟中的大规模线性系统，两遍 Lanczos 算法提供了内存与性能的实用平衡点。