linear-solvers

线性求解器

目标

提供通用工作流，用于为数值模拟中的线性系统选择求解器、评估条件数并诊断收敛情况。

要求

• Python 3.10+
• NumPy、SciPy（用于矩阵运算）
• 依赖项请参阅各个脚本

需要收集的输入

输入	描述	示例
矩阵大小	系统维度	n = 1000000
稀疏度	非零元素所占比例	0.01%
对称性	A = Aᵀ 是否成立？	是
定定性	A 是否正定？	是（SPD）
条件数	估计的条件数	10⁶

决策指南

求解器选择流程图

Is matrix dense and small enough to factor in memory (dense float64
storage n²·8 bytes < ~2 GB, i.e. n ≲ 16000)?
├── YES → Use direct solver (Cholesky/LDLᵀ/LU by symmetry)
└── NO → Is matrix symmetric?
    ├── YES → Is it positive definite?
    │   ├── YES → Use CG with AMG/IC preconditioner
    │   └── NO → Use MINRES
    └── NO → Is it nearly symmetric?
        ├── YES → Use BiCGSTAB
        └── NO → Use GMRES with ILU/AMG

快速参考

矩阵类型	求解器	预条件器
SPD，稀疏	CG	AMG、IC
对称不定	MINRES	SPD 预条件器（SSOR、对称块对角，或对 SPD 部分使用 AMG）
非对称	GMRES、BiCGSTAB	ILU、AMG
稠密	LU、Cholesky	无
鞍点	Schur 补、Uzawa	块预条件器

脚本输出（JSON 字段）

脚本	主要输出
scripts/solver_selector.py	recommended、alternatives、notes
scripts/convergence_diagnostics.py	rate、asymptotic_rate、stagnation、recommended_action
scripts/sparsity_stats.py	nnz、density、bandwidth、symmetry
scripts/preconditioner_advisor.py	suggested、notes
scripts/scaling_equilibration.py	row_scale、col_scale、notes
scripts/residual_norms.py	residual_norms、relative_norms、converged

工作流

1. 描述矩阵特征 - 对称性、定定性、稀疏度
2. 分析稀疏性 - 运行 scripts/sparsity_stats.py
3. 选择求解器 - 运行 scripts/solver_selector.py
4. 选择预条件器 - 运行 scripts/preconditioner_advisor.py
5. 应用缩放 - 如果条件数较差，使用 scripts/scaling_equilibration.py
6. 监控收敛情况 - 使用 scripts/convergence_diagnostics.py
7. 诊断问题 - 使用 scripts/residual_norms.py 检查残差历史

对话式工作流示例

用户：我的 GMRES 求解器在迭代 50 次后停滞了。残差降到 1e-3 后就不再改善。

Agent 工作流：

1. 诊断收敛情况：

   python3 scripts/convergence_diagnostics.py --residuals 1,0.1,0.01,0.005,0.003,0.002,0.002,0.002 --json
   

2. 查看预条件建议：

   python3 scripts/preconditioner_advisor.py --matrix-type nonsymmetric --sparse --ill-conditioned --json
   

3. 建议：增大重启参数，尝试使用更高 k 值的 ILU(k)，或切换到 AMG。

求解前检查清单

• 确认矩阵的对称性/定定性
• 根据大小和稀疏度决定使用直接法还是迭代法
• 根据物理量级设置残差容差
• 根据矩阵结构选择合适的预条件器
• 如有需要，应用缩放/均衡化
• 跟踪收敛情况，并在发生停滞时进行调整

CLI 示例

# Analyze sparsity pattern
python3 scripts/sparsity_stats.py --matrix A.npy --json

# Select solver for SPD sparse system
python3 scripts/solver_selector.py --symmetric --positive-definite --sparse --size 1000000 --json

# Get preconditioner recommendation
python3 scripts/preconditioner_advisor.py --matrix-type spd --sparse --json

# Diagnose convergence from residual history
python3 scripts/convergence_diagnostics.py --residuals 1,0.2,0.05,0.01 --json

# Apply scaling
python3 scripts/scaling_equilibration.py --matrix A.npy --symmetric --json

# Compute residual norms
python3 scripts/residual_norms.py --residual 1,0.1,0.01 --rhs 1,0,0 --json

错误处理

错误	原因	解决方法
Matrix file not found	路径无效	检查文件是否存在
Matrix must be square	输入矩阵不是方阵	验证矩阵维度
Residuals must be positive	残差数据无效	检查输入格式

结果解读指南

收敛速率

convergence_diagnostics.py 会报告两种速率：rate（整个历史记录中每次迭代残差比率的均值）和 asymptotic_rate（末尾一个较短窗口内的均值）。stagnation 标志由 asymptotic_rate（> 0.95）决定，因为停滞是尾部特征——前期快速下降可能掩盖尾部的平缓变化。判断下表所示的状态时，请参考 asymptotic_rate：

渐近速率	含义	操作
< 0.1	极佳	当前设置已最优
0.1 - 0.5	良好	对大多数问题而言可以接受
0.5 - 0.95	缓慢	考虑使用更好的预条件器
> 0.95	停滞	更换求解器或预条件器

停滞诊断

模式	可能原因	解决方法
残差趋于平稳	预条件器不佳	改进预条件器
振荡	矩阵接近奇异或不定	检查矩阵，尝试其他求解器
衰减非常缓慢	条件不佳	应用缩放，使用 AMG

验证清单

只有在以下各项都满足且记录了具体数值后，才能信任求解结果；不能只凭“看起来没问题”：

• 记录了 convergence_diagnostics.py 中的 asymptotic_rate，并确认其低于 0.95 的停滞阈值（理想情况下 < 0.5）；仅凭整个历史记录的 rate 较低，并不能排除尾部平稳的情况。
• 将 residual_norms.py 中的相对残差与按物理量纲缩放的 --rel-tol（默认值为 1e-6）进行比较，而不只是查看绝对范数；对于未缩放的 RHS，使用 --require-both，避免过小的 rhs 造成虚假的收敛结果。
• 确认 solver_selector.py 的 recommended 与 sparsity_stats.py 记录的实际矩阵属性（symmetry，以及已知情况下的定性）一致——例如，只有在矩阵对称且正定时才使用 CG；对称不定矩阵使用 MINRES；非对称矩阵使用 GMRES/BiCGSTAB。
• 对标记为条件不佳的系统，运行了 scaling_equilibration.py，并记录 row_scale_max/row_scale_min 和 col_scale_max/col_scale_min；对称矩阵使用 --symmetric（D A D）以保留对称性；非对称矩阵的双边缩放则先应用 row_scale，再应用 col_scale。
• 检查了 scaling_equilibration.py 中来自 sparsity_stats.py 的 notes/zero_rows/zero_cols——任何零行或零列都意味着系统在结构上是奇异的，而将尺度回退为 1 并不能解决问题。
• 确认 preconditioner_advisor.py 给出的预条件器适用于所选的 Krylov 方法——特别是，MINRES 的预条件器必须是 SPD（不定的非完全 LDLᵀ 无效）。

常见陷阱与自我安慰

看似省事的做法	错误原因 / 应采取的做法
“平均 rate 很低，所以已经收敛了。”	rate 是整个历史记录的平均值，容易受早期快速下降的影响；停滞是尾部特征。查看 asymptotic_rate，并确认其低于 0.95。
“绝对残差很小，所以可以结束了。”	如果 RHS 很大或未缩放，较小的绝对范数可能毫无意义。检查 relative_norms / relative_value 是否满足按物理量纲缩放的 --rel-tol。
“矩阵是对称的，所以直接用 CG。”	CG 要求矩阵对称且正定。对称不定矩阵需要使用 MINRES（并配备 SPD 预条件器）；使用 CG 可能会失效或停滞。选择前要确认正定性。
“系统很大，所以直接做分解。”	solver_selector.py 仅在稠密 float64 存储满足条件时才允许使用稠密直接求解器（n²·8 bytes < ~2 GB，n ≈ 16384）；超过该规模后，稠密 Cholesky/LU 不可行，必须改用迭代方法。
“缩放就是将每行除以该行的最大值。”	单边行缩放无法实现均衡缩放。对于非对称矩阵，应从行缩放后的矩阵导出 col_scale 并同时应用两种缩放；对于对称矩阵，应使用对称 D A D 缩放，否则会破坏对称性。
“GMRES 停滞了，那就增加迭代次数。”	尾部平稳意味着问题出在预条件器或重启长度，而不是迭代次数。增强预条件器（提高 ILU 填充级别 / 使用 AMG）、增大重启参数，或更换方法。

安全

输入验证

• 所有数值输入（残差、容差、矩阵元素）都会验证为有限数值
• 以逗号分隔的残差/向量输入最多允许 100,000 个条目
• solver_selector.py 的 --size 参数上限为 10 billion
• --matrix-type 会根据固定的允许列表（spd、symmetric-indefinite、nonsymmetric）进行验证
• 布尔标志（--symmetric、--positive-definite、--sparse、--ill-conditioned）是类型安全的 argparse 标志

文件访问

• sparsity_stats.py 和 scaling_equilibration.py 会读取由 --matrix 指定的单个矩阵文件（.npy 格式）
• 调用 np.load() 时会设置 allow_pickle=False，以防通过构造的 .npy 文件执行任意代码
• 矩阵文件在解析前会被检查；超过 500 MB 的文件将被拒绝
• 矩阵每个维度的大小上限为 100,000，以防内存耗尽
• 所有其他脚本都不会读取外部文件；输入通过 CLI 参数提供

工具限制

• 读取：用于检查脚本源代码、参考资料和矩阵文件
• 写入：用于保存分析结果或求解器建议；写入范围限于用户的工作目录
• Grep/Glob：用于查找相关文件和搜索引用
• 该技能的 allowed-tools 排除了 Bash，以防 Agent 在处理不受信任的矩阵文件或数值输入时执行任意命令

安全措施

• 不使用 eval()、exec() 或动态代码生成
• 所有子进程调用都使用显式参数列表（不设置 shell=True）
• 缩小工具范围（不提供 Bash），将 Agent 的操作限制为读取和写入
• JSON 输出模式会生成结构化、可解析的结果，不包含可由 shell 解释的内容

限制

• 大型稠密矩阵：直接求解器可能会耗尽内存
• 高度不定矩阵：标准预条件器可能会失效
• 鞍点问题：需要专用的分块预条件器

参考资料

• references/solver_decision_tree.md - 选择逻辑
• references/preconditioner_catalog.md - 预条件器选项
• references/convergence_patterns.md - 故障诊断
• references/scaling_guidelines.md - 均衡缩放指南

版本历史

• v1.2.0（2026-06-23）：修复渐近停滞检测、稠密可行性求解器筛选、鞍点/小型稠密矩阵的直接求解器路由、双侧均衡缩放、CG 迭代上限表，以及文档与评估的一致性
• v1.1.0（2024-12-24）：完善文档、决策指南和示例
• v1.0.0：首次发布，包含 6 个求解器分析脚本