简介本资源是一份面向数值计算、科学计算与高性能计算方向学习者与研究者的专业教学讲义聚焦Krylov子空间算法这一求解大规模稀疏线性方程组的核心迭代方法。内容系统覆盖投影原理、Arnoldi与Lanczos过程、FOM/GMRES/CG/MINRES/SYMMLQ/BiCG/QMR等主流算法推导与实现细节并包含Chebyshev多项式理论支撑下的收敛性分析适合研究生、算法工程师及需深入理解迭代法底层逻辑的进阶学习者。资源为单文件PDF文档共1个600KB的高清讲义结构清晰、公式严谨、章节完整含5大主模块与7个子节涵盖从基础投影理论到对称/非对称方程组的全场景算法选型与对比。目前已有422人学习下载可直接用于课堂补充、算法复现参考或工程问题建模前的理论夯实。1. Krylov 子空间算法不是“黑箱迭代器”而是带几何约束的投影求解器你手头正跑着一个含 10⁶ 量级自由度的结构力学仿真直接 LU 分解内存爆掉、时间超限你调用scipy.sparse.linalg.gmres却发现残差曲线在 1e-4 处反复震荡重启几次初始向量也没用你翻开源码想改restart参数却卡在Hessenberg矩阵更新逻辑里——这些都不是配置错误而是你没真正理解Krylov 子空间算法的本质是在低维几何空间中强制施加正交约束的投影过程而非泛泛的“迭代逼近”。它不靠猜解而靠构造一个与残差正交的子空间再把原问题“压扁”到这个子空间里求解。本讲 PDF 不是数学推导汇编它是工程师拆解大规模线性系统时必须掌握的空间映射手册从 Arnoldi 正交化如何避免 Gram-Schmidt 失稳到 GMRES 残差范数为何能实时监控无需显式计算b - Ax再到 CG 在对称正定场景下为何天然满足 A-共轭性——所有算法差异最终都落在“选哪个子空间 K”和“让残差垂直于哪个子空间 L”这两个选择上。适合已写过稀疏矩阵乘法、调试过condest报警、但面对maxiter和tol参数仍凭经验调整的中级以上数值计算实践者。2. 投影算子不是抽象符号而是可显式构造的矩阵映射工具Krylov 方法的根基不在迭代步数而在投影算子 P 的显式构造逻辑。PDF 第 5.1 节明确指出P 是幂等矩阵P² P其像空间 Ran(P) 和零空间 Ker(P) 共同唯一确定该投影。这直接对应工程实现中的两个关键动作基矩阵组装与约束子空间定义。若忽略此点盲目套用scipy接口将无法诊断为何 FOM 在非对称问题中发散或为何 MINRES 对奇异矩阵更鲁棒。2.1 正交投影与斜投影的矩阵表达差异正交投影要求 Ker(P) Ran(P)⊥即 P 必须对称Pᵀ P。此时若 V ∈ ℝⁿˣᵐ 是 Ran(P) 的标准正交基矩阵VᵀV Iₘ则投影矩阵为import numpy as np # V: (n, m) 标准正交基矩阵列满秩 P_orthogonal V V.T # (n, n) 对称幂等矩阵注意此处V V.T计算成本为 O(n²m)仅当 m ≪ n 时可行。实际 Krylov 迭代中我们从不显式构造 P而是通过 Arnoldi 向量v_j隐式操作——这是避免存储大型投影矩阵的核心设计。而斜投影oblique projection不要求对称性其构造依赖两个子空间像空间 V 和零空间 L满足 ℝⁿ V ⊕ L。PDF 式 (5.2) 给出通用形式$$ P V(W^T V)^{-1} W^T $$其中 W 是 L 的基矩阵。此时W^T V必须可逆定理 5.3 给出充分条件A 正定且 L K或 A 非奇异且 L AK。该式揭示了算法选型的底层逻辑FOM/CG取 L K → W V →W^T V I→P V V^T退化为正交投影GMRES取 L AK → W AV →W^T V V^T A^T V非对称需解 Hessenberg 系统BiCG取 L K(Aᵀ, r₀) → W 由 Aᵀ 生成的 Krylov 基构成2.2 Petrov-Galerkin 条件的工程实现路径PDF 式 (5.6) 定义的投影算法find x̃ ∈ K s.t. b - A x̃ ⊥ L在代码中转化为三步硬性操作基矩阵构建设K span{v₁,...,vₘ},L span{w₁,...,wₘ}组装V [v₁...vₘ],W [w₁...wₘ]投影矩阵隐式求解由正交条件Wᵀ(b - A x̃) 0且x̃ x⁽⁰⁾ V y得$$ W^T A V y W^T r_0 \quad \text{where } r_0 b - A x^{(0)} $$小规模系统求解解m × m系统(W^T A V) y W^T r_0再计算x̃ x⁽⁰⁾ V y该流程在scipy.sparse.linalg中被封装但理解其矩阵维度至关重要符号形状工程含义常见陷阱V(n, m)Krylov 基矩阵Arnoldi 向量堆叠m过大导致内存溢出需restartW(n, m)约束子空间基决定算法类型W AV时W^T A V是 Hessenberg 矩阵非对称W^T A V(m, m)投影后的系数矩阵若m1000直接求逆耗时 O(m³)应改用 QR 或 Givens2.3 Arnoldi 过程的数值稳定性实操要点PDF 算法 5.1经典 Gram-Schmidt与算法 5.2修正 Gram-Schmidt, MGS的差异在浮点计算中会导致数量级误差。以下 Python 片段演示 MGS 如何抑制正交性丢失def arnoldi_mgs(A, r0, m): n len(r0) V np.zeros((n, m1)) # v1..v_{m1} H np.zeros((m1, m)) # Hessenberg matrix beta np.linalg.norm(r0) V[:, 0] r0 / beta for j in range(m): w A V[:, j] # Step: A * v_j # MGS: 逐次正交化每次减去当前 v_i 分量 for i in range(j1): H[i, j] w V[:, i] # w, v_i w w - H[i, j] * V[:, i] H[j1, j] np.linalg.norm(w) if H[j1, j] 1e-15: break V[:, j1] w / H[j1, j] return V[:, :j1], H[:j2, :j1] # 使用示例构造 K₃(A, r₀) 的正交基 A_sparse scipy.sparse.csr_matrix(...) # 大型稀疏矩阵 r0 b - A_sparse x0 V, H arnoldi_mgs(lambda x: A_sparse x, r0, m3) print(fKrylov 基维度: {V.shape[1]}, Hessenberg 形状: {H.shape})提示MGS 中w w - H[i,j] * V[:,i]的顺序不可颠倒。经典 GS 先算所有内积再统一减MGS 每算一个内积立即修正w显著提升V列向量间的正交性np.max(np.abs(V.T V - np.eye(V.shape[1])))可降至 1e-13 量级。3. GMRES 与 FOM 的核心分野残差极小化 vs 投影正交化当面对一般非对称矩阵 A 时GMRES 与 FOM 是最常被调用的两个 Krylov 算法但它们的数学目标与工程行为截然不同。PDF 第 5.3.1–5.3.2 节明确区分FOM 强制残差r̃ b - A x̃垂直于 Krylov 子空间 Kₘ而 GMRES 则最小化∥r̃∥₂。这一差异导致二者在收敛性、稳定性及实现细节上存在本质区别。3.1 FOM投影正交化的直接实现与失效场景FOM 的核心是求解Hₘ y β e₁PDF 式 5.13其中Hₘ是m×m上 Hessenberg 矩阵。其残差表达式PDF 定理 5.6为$$ \tilde{r} -h_{m1,m} (e_m^T y) v_{m1} $$这意味着FOM 残差方向始终与下一个 Arnoldi 向量v_{m1}平行。该性质带来两个关键后果收敛不可控若h_{m1,m}极小如 A 有接近不变子空间则∥r̃∥₂可能远大于实际精度需求算法提前终止却未达精度对称矩阵失效当 A 对称时Hₘ为三对角矩阵但 FOM 仍按非对称逻辑求解无法利用 CG 的 A-共轭性优势。以下代码验证 FOM 残差方向特性# 假设已运行 arnoldi_mgs 得到 V, H, beta m H.shape[1] H_m H[:m, :] # (m, m) 上 Hessenberg e1 np.zeros(m); e1[0] 1 y_fom np.linalg.solve(H_m, beta * e1) # 计算残差范数无需显式 x̃ h_next H[m, m-1] # h_{m1,m} r_norm_fom abs(h_next * y_fom[-1]) # PDF 定理 5.6 直接给出 # 对比显式计算残差验证一致性 x_fom x0 V[:, :m] y_fom r_explicit b - A_sparse x_fom print(fFOM 残差范数公式: {r_norm_fom:.2e}) print(fFOM 残差范数显式: {np.linalg.norm(r_explicit):.2e}) # 二者应完全一致浮点误差内3.2 GMRES残差极小化的 QR 分解实现GMRES 的目标是min_{y∈ℝᵐ} ∥β e₁ - H_{m1,m} y∥₂PDF 式 5.15 推导即在m维空间中找y使H_{m1,m} y最接近β e₁。这等价于求解一个(m1)×m超定系统。PDF 第 5.3.3 节强调使用 Givens 旋转进行 QR 分解原因在于H_{m1,m}是上 Hessenberg 矩阵Givens 旋转可将其变为上三角R同时保持Q的隐式结构分解后系统变为R y Q^T (β e₁)R为m×m上三角回代即可得y残差范数∥r̃∥₂直接等于|Q^T (β e₁) 的最后一行元素|无需计算x̃。以下是 GMRES 关键步骤的 NumPy 实现简化版def gmres_residual_norm(H, beta, m): 计算 GMRES 第 m 步残差范数无需显式 QR 分解 H: (m1, m) Hessenberg 矩阵 返回: 残差范数 ∥r̃∥₂ # 初始化 Givens 旋转参数 g np.array([beta] [0] * m) # Q^T * (beta * e1)初始为 [beta, 0, ..., 0] # 对 H 的每一列应用 Givens 旋转同步更新 g for j in range(m): # 计算第 j 列的 Givens 旋转消去 H[j1, j] a, b H[j, j], H[j1, j] if abs(b) 1e-15: continue r np.sqrt(a*a b*b) c, s a/r, b/r # cos, sin # 应用旋转到 H 的第 j 列j 行和 j1 行 H[j, j] r H[j1, j] 0 # 应用旋转到 g 的第 j 和 j1 行 g_j, g_j1 g[j], g[j1] g[j] c * g_j s * g_j1 g[j1] -s * g_j c * g_j1 # 残差范数 |g[m]|最后一行 return abs(g[m]) # 使用示例 r_norm_gmres gmres_residual_norm(H, beta, mH.shape[1]) print(fGMRES 第 {m} 步残差范数: {r_norm_gmres:.2e})关键对比表FOM 与 GMRES 在工程实践中的决策依据特性FOMGMRES数学目标r̃ ⊥ Kₘ残差正交于 Krylov 子空间min ∥r̃∥₂残差 2-范数极小化残差监控∥r̃∥₂ |h_{m1,m} y_m|直接计算∥r̃∥₂ |g_m|Givens 后g向量末元素矩阵需求Hₘ需可逆若奇异则失败H_{m1,m}总可 QR 分解鲁棒性强适用场景A 接近对称正定或需严格正交约束一般非对称矩阵工业级首选scipy 接口无直接对应需手动实现scipy.sparse.linalg.gmres4. 对称正定系统的最优解CG 算法的 A-共轭性本质与收敛边界当系数矩阵 A 满足对称Aᵀ A且正定xᵀAx 0, ∀x ≠ 0时Krylov 方法迎来质变共轭梯度法CG成为理论最优解。PDF 第 5.4.2 节指出CG 的核心并非简单迭代而是在 Krylov 子空间中自动构造 A-共轭方向使得每一步搜索方向d_k满足d_i^T A d_j 0 (i ≠ j)。这一性质直接导致 CG 在至多 n 步内精确收敛n 为矩阵阶数且每步仅需一次矩阵-向量乘法与内积计算。4.1 CG 的 A-共轭性如何从 Lanczos 过程自然导出PDF 第 5.4.1 节的 Lanczos 过程是 Arnoldi 在对称矩阵下的特例因 A 对称Hessenberg 矩阵Hₘ退化为三对角矩阵Tₘ。此时 Arnoldi 关系式A Vₘ Vₘ₊₁ Hₘ₊₁,ₘ变为$$ A V_m V_m T_m t_{m1,m} v_{m1} e_m^T $$其中Tₘ是m×m对称三对角矩阵。Vₘ的列向量v_j构成 A-共轭基v_i^T A v_j 0 (i ≠ j)。CG 的搜索方向d_k正是这些v_j的线性组合因此天然满足 A-共轭性。以下 Python 代码展示 Lanczos 过程如何生成 A-共轭基def lanczos_symmetric(A, r0, m): 对称矩阵专用 Lanczos 过程 n len(r0) V np.zeros((n, m1)) alpha np.zeros(m) # 对角元 beta np.zeros(m) # 次对角元 beta[0] np.linalg.norm(r0) V[:, 0] r0 / beta[0] for j in range(m): w A V[:, j] if j 0: w w - beta[j-1] * V[:, j-1] # 减去前一项 alpha[j] w V[:, j] w w - alpha[j] * V[:, j] if j m-1: beta[j] np.linalg.norm(w) if beta[j] 1e-15: break V[:, j1] w / beta[j] # 构造三对角矩阵 T_m T np.diag(alpha[:m]) if m 1: T np.diag(beta[:m-1], k1) np.diag(beta[:m-1], k-1) return V[:, :len(alpha)], T # 示例验证 A-共轭性 A_sym (A_sparse A_sparse.T) / 2 # 强制对称 V_lan, T lanczos_symmetric(A_sym, r0, m5) # 检查 v_i^T A v_j ≈ 0 (i≠j) A_conjugate_check V_lan.T A_sym V_lan print(A-共轭性检查 (V^T A V):) print(np.round(A_conjugate_check, decimals12)) # 应近似为对角矩阵4.2 CG 收敛速度的显式上界特征值分布决定一切PDF 第 5.5.2 节给出 CG 收敛性的经典估计定理 5.10$$ \frac{|e_k|_A}{|e_0|_A} \leq 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} 1} \right)^k $$其中κ λ_max / λ_min是 A 的条件数∥e∥_A √(eᵀ A e)是 A-范数。该式揭示CG 收敛速度由 A 的特征值分布宽度决定而非矩阵大小。若κ 100则k ≈ 10步即可将误差降低 100 倍若κ 10⁴则需k ≈ 100步。工程实践中可通过预处理preconditioning压缩κ。例如对角预处理M diag(A)的 CG 实现def cg_preconditioned(A, b, x0, M_diag, max_iter100, tol1e-8): x x0.copy() r b - A x z r / M_diag # 预处理解 M z r d z.copy() for k in range(max_iter): Ad A d alpha (r z) / (d Ad) x x alpha * d r r - alpha * Ad z r / M_diag beta (r z) / (r z) if k max_iter-1 else 0 d z beta * d if np.linalg.norm(r) tol * np.linalg.norm(b): print(fCG 收敛于第 {k1} 步) break return x # 使用对角预处理加速 M_diag np.array(A_sparse.diagonal()) x_cg cg_preconditioned(A_sparse, b, x0, M_diag, tol1e-6)提示预处理矩阵M应满足M ≈ A且M⁻¹z易解。对角预处理最简单但块 Jacobi 或不完全 CholeskyIC预处理可进一步压缩κ尤其适用于偏微分方程离散矩阵。5. Krylov 算法实战排错从残差震荡到 Arnoldi 失稳的定位链在真实项目中Krylov 算法失败极少源于理论缺陷而多因子空间构造失稳或约束条件误设。PDF 全文贯穿的 Arnoldi 过程、Hessenberg 矩阵、残差正交性等概念正是排错的黄金线索。以下提供一套可直接复用的诊断流程覆盖从scipy接口报错到自研实现崩溃的全场景。5.1 残差曲线异常的三级归因与验证命令当gmres残差不下降甚至震荡时按优先级执行以下检查一级检查矩阵-向量乘法正确性# 用 scipy 自带的测试函数验证 A x 是否与预期一致 from scipy.sparse.linalg import LinearOperator A_op LinearOperator(shapeA_sparse.shape, matveclambda x: A_sparse x) # 测试几个随机向量 for _ in range(3): x_test np.random.randn(A_sparse.shape[1]) y_scipy A_op x_test y_manual your_matvec_func(x_test) # 替换为你的实现 assert np.allclose(y_scipy, y_manual, rtol1e-10), 矩阵-向量乘法错误二级验证 Arnoldi 向量正交性# 运行 Arnoldi 至 m50检查 V^T V 是否接近单位阵 V, H arnoldi_mgs(lambda x: A_sparse x, r0, m50) orthogonality_error np.max(np.abs(V.T V - np.eye(V.shape[1]))) print(fArnoldi 向量正交性误差: {orthogonality_error:.2e}) # 1e-8 表明 MGS 失效需改用 Householder 或增加精度三级分析 Hessenberg 矩阵病态性# 提取 H_m 并计算条件数 m min(50, H.shape[1]) H_m H[:m, :m] cond_H np.linalg.cond(H_m) print(fH_{m} 条件数: {cond_H:.2e}) # 1e12 表明 Krylov 子空间接近线性相关需重启restart或预处理5.2 “break due to h_{j1,j}0”的深层解读与应对PDF 算法 5.2 第 9–10 行的if h_{j1,j} 0 then break是算法健康的标志而非错误。其物理含义是K_j(A,r₀)已成为 A 的不变子空间PDF 性质 5.9此时r_j 0x_j即为精确解。但若j远小于n如j3而n10000则表明矩阵 A 有极小秩rank(A) ≤ j问题本身欠定初始残差 r₀ 位于低维不变子空间常见于周期性结构或重复单元模型。验证方法# 检查 r_j 是否真为零 j_break 3 # 假设在第 3 步中断 V_j V[:, :j_break] r_j b - A_sparse (x0 V_j y_j) # y_j 为前 j 步解 print(f中断步残差范数: {np.linalg.norm(r_j):.2e}) # 若 ≈ 0则确认为精确解否则检查 Arnoldi 实现5.3 预处理失效的快速检测表预处理是 Krylov 加速的关键但错误预处理会加剧发散。以下表格列出常见预处理类型及其失效信号预处理类型正常表现失效信号验证命令对角预处理M diag(A)残差单调下降残差震荡cond(M⁻¹A)cond(A)np.linalg.cond(np.diag(1/M_diag) A_sparse.toarray())不完全 CholeskyM L Lᵀ∥r_k∥下降速率提升L分解失败nan或infscipy.sparse.linalg.splu(A_sparse, permc_specNATURAL)Jacobi 迭代M D每步M⁻¹r计算稳定M⁻¹r结果含nannp.any(np.isnan(1/np.diag(A_sparse) * r))终极技巧当所有诊断均无异常但算法仍不收敛时强制设置restart10并启用callback函数记录每步y向量def callback(xk): print(fStep {len(callback.history)1}: ∥r∥{np.linalg.norm(b - A_sparse xk):.2e}) callback.history.append(xk) callback.history [] scipy.sparse.linalg.gmres(A_sparse, b, x0x0, restart10, callbackcallback)观察callback.history中xk的变化模式——若xk在固定子空间内循环则证实 Krylov 子空间退化需更换初始向量或引入随机扰动。本文还有配套的精品资源点击获取