★ FEATURED ARTICLE 1. 为什么还要手写迭代法从一次内存告急说起很多人第一次接触线性方程组脑子里蹦出来的都是高斯消元——列主元、回代、完事。我当年也这样直到有一次处理一个稀疏矩阵问题矩阵规模大概三千阶非零元占比不到百分之一。用高斯消元跑内存直接顶到几个G机器开始疯狂换页风扇像要起飞。后来换成迭代法内存占用瞬间降到几十兆迭代几百轮就收敛了。那次之后我才真正理解迭代解法和直接解法不是谁替代谁的关系而是各自守着不同的战场。线性方程组的迭代解法说白了就是从一个初始猜测出发反复套用同一个修正公式让近似解一步步逼近真实解。它不追求一次性把矩阵分解干净而是把解方程这件事拆成很多次微调。这个思路在处理大型稀疏矩阵时特别香因为迭代法每次只需要做矩阵和向量的乘法不需要改变矩阵本身的结构稀疏性可以被完整保留下来。这篇文章面向的是已经会写C语言、懂基本矩阵运算但想搞清楚迭代法到底怎么落地的人。我会把雅可比、高斯-赛德尔、超松弛这几类经典方法的C语言实现讲透包括它们各自的收敛条件、代码骨架、实测中的坑以及怎么判断到底收敛了没有。所有代码都是可以直接编译运行的完整片段不是伪代码。你跟着敲一遍基本就能把迭代法这套东西装进自己的工具箱。需要提前说明的是迭代法不是万能的。它能不能收敛取决于系数矩阵的性质。如果矩阵不满足收敛条件你迭代一万次也是白搭甚至数值会越跑越大。所以我会花不少篇幅讲收敛性判断这比代码本身更重要。2. 迭代法的数学骨架把Axb拆成可循环的形式2.1 从矩阵分裂说起所有迭代法的底层逻辑都来自同一个动作——矩阵分裂。把系数矩阵A拆成两部分A M - N其中M要求是可逆的而且最好容易求逆比如对角阵、三角阵。代回原方程Axb得到Mx - Nx b Mx Nx b x M⁻¹Nx M⁻¹b这就自然导出了迭代格式x^(k1) M⁻¹N x^(k) M⁻¹ b令 B M⁻¹Nf M⁻¹b就得到标准形式 x^(k1) Bx^(k) f。这里的B叫迭代矩阵它的谱半径ρ(B)决定了迭代是否收敛。ρ(B) 1 是收敛的充要条件这是整个迭代法理论的基石。不同的M选取方式对应不同的迭代方法。选M为A的对角部分就是雅可比选M为A的下三角部分就是高斯-赛德尔在下三角基础上加个松弛因子就是SOR。理解了这个框架后面几种方法就不再是孤立的公式而是同一个思想的不同变体。2.2 雅可比迭代最直观但收敛最慢雅可比迭代的思路极其朴素解第i个方程的时候把其他所有未知数都用上一轮的旧值代入。对于方程a₁₁x₁ a₁₂x₂ ... a₁ₙxₙ b₁解出x₁x₁ (b₁ - a₁₂x₂ - ... - a₁ₙxₙ) / a₁₁写成迭代格式就是x_i^(k1) (b_i - Σ_{j≠i} a_ij x_j^(k)) / a_ii这个公式的好处是每一轮所有分量的更新互不依赖天然适合并行。坏处也明显它把上一轮所有旧值都用上了信息更新滞后收敛速度通常是几种方法里最慢的。雅可比收敛的一个充分条件是矩阵严格对角占优也就是每一行对角元素的绝对值大于该行其他元素绝对值之和。这个条件在实际问题里经常能满足比如很多离散化后的偏微分方程就天然对角占优。2.3 高斯-赛德尔用上新值快人一步高斯-赛德尔和雅可比的唯一区别在于计算x_i^(k1)时前面已经更新过的x_1到x_{i-1}直接用新值不再用旧值。x_i^(k1) (b_i - Σ_{ji} a_ij x_j^(k1) - Σ_{ji} a_ij x_j^(k)) / a_ii别小看这个改动。它相当于在迭代过程中不断吸收最新信息收敛速度通常能比雅可比快一倍左右。代价是失去了并行性因为每个分量必须等前面的算完才能算。不过在单机串行场景下这个代价完全值得。高斯-赛德尔的收敛条件比雅可比宽松一些如果A是对称正定矩阵或者严格对角占优它一定收敛。注意这里有个反直觉的点——存在雅可比收敛但高斯-赛德尔发散的情况所以不能想当然认为后者一定更好。2.4 SOR给收敛速度加个油门逐次超松弛SOR是在高斯-赛德尔基础上引入松弛因子ωx_i^(k1) (1-ω)x_i^(k) ω * x_i^(GS)其中x_i^(GS)是高斯-赛德尔算出来的值。ω1时退化成高斯-赛德尔0ω1叫低松弛用于某些需要稳定化的场景1ω2叫超松弛用来加速收敛。ω的选取是SOR的灵魂。理论最优值ω_opt 2 / (1 sqrt(1 - ρ(J)²))其中ρ(J)是雅可比迭代矩阵的谱半径。但实际算ρ(J)本身就不便宜所以工程上常用试凑法从1.0开始以0.1为步长往上试观察迭代次数变化找到拐点。我实测下来很多问题在ω1.5到1.8之间效果最好但ω超过某个阈值后会直接发散这个边界必须小心。3. C语言实现从数据结构到收敛判据3.1 矩阵存储二维数组够用但要省着用先明确一点迭代法处理的是稠密小矩阵还是稀疏大矩阵存储策略完全不同。如果是教学或中小规模问题直接用二维数组最省事#define N 100 double A[N][N]; double b[N]; double x[N];但如果矩阵规模上千且稀疏二维数组会浪费大量内存。这时候应该用**压缩稀疏行CSR**格式只存非零元的值、列索引和行偏移。不过CSR的矩阵向量乘法写起来复杂一些本文为了聚焦迭代逻辑本身先用二维数组演示最后会提一下稀疏场景的改造思路。一个完整的迭代求解器需要这些变量double A[N][N]; // 系数矩阵 double b[N]; // 右端向量 double x[N]; // 当前解 double x_new[N]; // 新一轮解 double eps; // 收敛精度 int max_iter; // 最大迭代次数3.2 雅可比迭代的完整实现#include stdio.h #include math.h #define N 4 int jacobi(double A[N][N], double b[N], double x[N], double eps, int max_iter) { double x_new[N]; int iter, i, j; double sum, diff, max_diff; for (iter 0; iter max_iter; iter) { max_diff 0.0; for (i 0; i N; i) { sum b[i]; for (j 0; j N; j) { if (j ! i) { sum - A[i][j] * x[j]; } } x_new[i] sum / A[i][i]; } // 计算最大变化量作为收敛判据 for (i 0; i N; i) { diff fabs(x_new[i] - x[i]); if (diff max_diff) max_diff diff; x[i] x_new[i]; } if (max_diff eps) { return iter 1; // 返回实际迭代次数 } } return -1; // 未收敛 }这段代码有几个细节值得说。收敛判据用的是相邻两轮解的最大分量变化量而不是残差‖Ax-b‖。两者都可以但前者计算量小且能直接反映迭代是否稳定下来。实践中我倾向于两个都算如果变化量小但残差还大说明可能陷入了缓慢漂移需要警惕。另外注意x_new和x的分离。雅可比必须用完整的旧值所以不能边算边覆盖。这个临时数组是必须的别为了省内存直接在x上原地更新那就变成高斯-赛德尔了。3.3 高斯-赛德尔省掉临时数组int gauss_seidel(double A[N][N], double b[N], double x[N], double eps, int max_iter) { int iter, i, j; double sum, diff, max_diff; for (iter 0; iter max_iter; iter) { max_diff 0.0; for (i 0; i N; i) { sum b[i]; for (j 0; j N; j) { if (j ! i) { sum - A[i][j] * x[j]; } } diff fabs(sum / A[i][i] - x[i]); x[i] sum / A[i][i]; if (diff max_diff) max_diff diff; } if (max_diff eps) { return iter 1; } } return -1; }高斯-赛德尔可以直接在x上原地更新因为计算x_i时用到的x_jji已经是新值x_jji还是旧值正好符合公式要求。这个原地更新不仅省内存还让代码更简洁。但要注意收敛判据必须在更新前算好diff否则新旧值混在一起就乱了。3.4 SOR加一个ω参数int sor(double A[N][N], double b[N], double x[N], double omega, double eps, int max_iter) { int iter, i, j; double sum, diff, max_diff, x_gs; for (iter 0; iter max_iter; iter) { max_diff 0.0; for (i 0; i N; i) { sum b[i]; for (j 0; j N; j) { if (j ! i) { sum - A[i][j] * x[j]; } } x_gs sum / A[i][i]; diff fabs(omega * (x_gs - x[i])); x[i] (1.0 - omega) * x[i] omega * x_gs; if (diff max_diff) max_diff diff; } if (max_diff eps) { return iter 1; } } return -1; }SOR的收敛判据我用了omega * (x_gs - x[i])因为实际更新量就是这个。如果直接用x_gs - x[i]当ω很小时会高估收敛速度导致提前退出。这个细节很多教程不讲但实测中确实会带来精度问题。4. 收敛性判断比写代码更关键的环节4.1 对角占优最快的预判手段拿到一个矩阵先别急着写迭代。花三十秒检查一下是否严格对角占优int is_diagonally_dominant(double A[N][N]) { int i, j; double off_sum; for (i 0; i N; i) { off_sum 0.0; for (j 0; j N; j) { if (j ! i) off_sum fabs(A[i][j]); } if (fabs(A[i][i]) off_sum) return 0; } return 1; }如果返回1雅可比和高斯-赛德尔都保证收敛你可以放心用。如果返回0不代表一定发散但需要进一步分析。我遇到过不少矩阵对角占优条件不满足但实际迭代照样收敛只是速度慢一些。4.2 谱半径估计理论上的金标准迭代矩阵B的谱半径ρ(B)小于1才收敛。但直接算特征值代价太高工程上常用幂法估计谱半径反复用B乘以一个随机向量观察向量模长的增长比率。如果比率持续小于1说明收敛。不过说实话在写求解器的时候专门去估谱半径有点本末倒置。更实用的做法是设一个合理的最大迭代次数跑一遍看结果。如果迭代次数内收敛了说明没问题如果没收敛再回头分析矩阵性质。4.3 残差监控运行时的安全网不管理论分析如何运行时都应该监控残差double compute_residual(double A[N][N], double b[N], double x[N]) { int i, j; double r, norm 0.0; for (i 0; i N; i) { r b[i]; for (j 0; j N; j) { r - A[i][j] * x[j]; } norm r * r; } return sqrt(norm); }我的习惯是每迭代50轮打印一次残差。如果残差在下降但很慢可以继续等如果残差开始震荡或上升立刻停说明方法选错了或者矩阵本身有问题。这个监控成本很低但能帮你省下大量调试时间。5. 实测中的坑与调优经验5.1 初始值不是随便给的很多人图省事把初始解全设成0。对于很多问题这没问题但如果b向量里有个别分量特别大从0出发可能导致前几轮迭代数值剧烈跳动甚至溢出。我的做法是用b除以对角元作为初始猜测for (i 0; i N; i) { x[i] b[i] / A[i][i]; }这个初始值虽然粗糙但量级上接近真实解能显著减少迭代次数。实测中好的初始值有时能把迭代次数砍掉三分之一。5.2 收敛精度别设太死eps设成1e-10看起来很严谨但如果你的问题本身数据精度只有1e-6迭代到1e-10纯属浪费。更麻烦的是浮点数在接近机器精度时会陷入假收敛——相邻两轮变化量很小但离真实解还有距离。我一般根据问题规模设eps小规模用1e-8大规模用1e-6配合残差监控一起判断。5.3 最大迭代次数是保险丝max_iter一定要设而且不能太大。我见过有人设成100000结果矩阵不收敛程序跑了半小时才退出。一般设成矩阵阶数的5到10倍就够了。如果这个范围内没收敛基本可以判定方法不适用换直接法或者预处理技术更实际。5.4 对角元为零怎么办如果某行对角元是0雅可比和高斯-赛德尔直接除零崩溃。这时候需要行交换把该行和一个对角元非零的行互换。但行交换会破坏对角占优性质所以交换后要重新检查。如果整个矩阵的对角元都接近零说明这个矩阵根本不适合迭代法考虑换列主元高斯消元。5.5 SOR的ω怎么试前面说了理论最优值但实际试的时候有个技巧先跑高斯-赛德尔记录迭代次数然后从ω1.2开始每次加0.1观察迭代次数变化。迭代次数先降后升最低点附近就是好值。注意ω接近2时非常危险稍微超过就发散所以试到1.9就停别贪心。6. 稀疏矩阵场景的改造思路前面用的二维数组在稠密小矩阵上没问题但真实工程问题往往是稀疏的。这时候需要把存储换成CSR格式typedef struct { double *val; // 非零元值 int *col_idx; // 列索引 int *row_ptr; // 行偏移长度N1 int nnz; // 非零元个数 } CSRMatrix;矩阵向量乘法变成void csr_matvec(CSRMatrix *A, double *x, double *y) { int i, j; for (i 0; i A-n; i) { y[i] 0.0; for (j A-row_ptr[i]; j A-row_ptr[i1]; j) { y[i] A-val[j] * x[A-col_idx[j]]; } } }迭代法的核心循环不用大改只是把内层的A[i][j] * x[j]换成CSR的遍历。这样内存占用能从O(N²)降到O(nnz)对于非零元占比1%的矩阵内存直接省下两个数量级。这是迭代法相比直接法最大的优势所在也是我当年转向迭代法的直接原因。不过CSR的雅可比实现有个小坑找对角元需要额外遍历。可以在构建CSR时单独存一个对角元数组避免每次迭代都去搜索。7. 几种方法的实测对比与选型建议为了让你有个直观感受我用一个4阶严格对角占优矩阵跑了三种方法eps1e-8方法迭代次数特点雅可比42实现最简单可并行收敛最慢高斯-赛德尔23串行最快代码简洁推荐首选SOR(ω1.5)14需要调参调好了最快选型逻辑很直接如果追求实现简单且能并行选雅可比如果单机串行且矩阵对角占优高斯-赛德尔基本够用如果迭代次数成为瓶颈且愿意花时间调ω上SOR。但记住SOR的ω对问题敏感换一个矩阵可能就要重新试通用性不如前两者。还有一个常被忽略的点迭代法的收敛速度还和矩阵的条件数有关。条件数越大收敛越慢。如果发现迭代次数远超预期可以考虑对矩阵做预处理比如对角缩放把每行除以对角元改善条件数。这个改造很简单但效果有时很显著。最后说个我自己的教训。早期写迭代求解器时我只关注代码能不能跑通忽略了收敛性验证。结果有一次处理一个非对角占优矩阵程序收敛了但解完全不对——因为迭代在某个错误值附近震荡变化量很小但残差很大。后来加了残差监控才抓到这个问题。迭代法的输出一定要用残差验证不能只看相邻两轮的变化量这是用多少次调试换来的经验。
阅读完成 · 觉得有帮助？
