简介经典代数多重网格AMG的Matlab演示实现面向数值计算方向的研究生、工程师及算法爱好者帮助理解利用AMG求解大型稀疏线性系统的基本流程。代码模块化程度高包含V循环迭代、强依赖矩阵构造、粗网格生成与插值等核心步骤便于逐函数阅读。包内附带二维Poisson方程在均匀三角网格上的有限元测试数据覆盖51×51与101×101两套规模并给出几何多重网格GMG的对比案例方便通过数值实验观察两类多重网格策略的收敛差异同时注明了粗网格生成较慢、未包含经典AMG生成粗网格的第二遍处理等边界条件降低二次开发困惑。整套资源有20个文件以11个m源码为主辅以4个mat数据、2张结果图和README说明压缩包仅309KB小体积易下载。已有1227人浏览学习适合想从代码层面入门AMG的读者。1. 代数多重网格是什么——一个Demo解决什么问题代数多重网格Algebraic MultigridAMG是所有做大规模稀疏线性方程组求解的人绕不开的一类方法。和几何多重网格不同AMG不需要网格的几何信息只靠系数矩阵A本身来构造一层一层的粗网格算子这意味着它可以直接当作一个“黑匣子”求解器接到有限元、有限差分或图像处理代码里。Classic_AMG_Demo就是用MATLAB演示经典AMG完整流程的最小工程从输入矩阵A和右端项b开始经过粗化、插值、松弛迭代三步输出一个可用的解x和收敛曲线。适合刚接触多重网格、想在工程里评估AMG是否值得投入的从业者也适合做预条件子时想理解V-cycle内部行为的同学。2. 从光滑误差到粗网格校正拆穿AMG的V-cycle2.1 为什么单层迭代在工程矩阵上掉速低频误差是主犯经典的逐点雅可比或Gauss-Seidel迭代表面上看起来每轮都在推进但收敛曲线到后期几乎是一条水平的直线。原因在误差的谱分解上迭代法对高频误差分量的衰减速度很快对低频误差分量基本无感。所谓“光滑误差”指的不是误差本身变得光滑而是经过几轮松弛之后残差里剩下的全是低频分量——这些分量波长长局部修正几乎碰不动它们。AMG的基本策略是三步先用松弛迭代把误差里的高频分量压掉然后把低频误差放到一个更粗的网格上去看在粗网格上原本低频的分量相对变得高频可以用同样的松弛方法继续压最后把修正量插值回细网格。整个过程形成的递归路径叫V-cycle。Classic_AMG_Demo演示的正是这条路径在MATLAB里如何一步步落到代码上而不是停留在公式层面。看一个实际场景求解泊松方程离散出的五点差分矩阵规模从1000阶涨到100万阶。Jacobi迭代的迭代次数随网格细化成倍增加而AMG因为拥有从细到粗整个层级链收敛率对规模几乎不敏感。这个性质在工程上的价值非常直接——矩阵规模翻10倍求解时间可能只翻有限几倍。这也是AMG能在CFD、结构力学等领域被反复使用的根本原因。2.2 RS粗化的C/F划分把未知量分成“留在这层”和“下放到粗层”经典AMG的粗化coarsening通常用Ruge-Stüben策略简称RS粗化。这一步不关心几何网格只依赖矩阵的强连接关系。强连接的判定是粗化最重要的参数之一。对第i行如果某个非对角元a_ij满足-a_ij theta * max(-a_ik, 0) (k不等于i)就认为变量j强影响着变量i。常见阈值为0.25通常不好用的参数区间为0.5以上或0.1以下我后面专门讲这个参数。C/F划分的粗化过程是先根据每个点强连接的邻居数量建立一个影响度估计从影响度最高的点开始把该点选为C点同时把受它强影响的所有点标记为F点。继续选直到所有点都标完。MATLAB里RS粗化的实现通常包含一个点数较少的排序循环示意代码如下% 计算强连接矩阵S(i,j)1 表示 第i行第j列强影响 S abs(A) theta * repmat(max(abs(A), [], 2), 1, n); S S - eye(n); % 去掉对角线自环 % 每个点的影响度 它影响多少个其他点 它被多少其他点影响 lambda full(sum(S, 1)) full(sum(S, 2)); % 标记初始全部为F点按影响度从大到小依次选C isC false(n, 1); for iter 1:n [~, idx] max(lambda .* (1 - isC)); % 挑当前影响度最大的未标记点 isC(idx) true; S_idx find(S(idx, :)); % 受该点强影响的未标记F点 % 更新剩余点的影响度粗化循环的核心 lambda(S_idx) lambda(S_idx) 1; end这段代码是RS粗化的极简示意不是Demo里的原样函数但逻辑一致。核心要知道三点强连接阈值太大会让C点几乎挑不出来层级只有两三层阈值太小会把几乎所有点都放进C点粗网格矩阵和细网格一样大失去意义。另外这段代码每轮用max扫描来选择下一个C点在百万级矩阵上会非常慢实际工程代码要用排序加队列或优先队列来改写但作为演示版已经足够说明原理。2.3 插值与限制算子Galerkin粗化怎么把矩阵“投影”到粗层选完C/F点之后需要构造插值算子P作用是把粗网格上的误差修正量映射回细网格。经典RS插值的规则是对每个F点用其强连接的C点和F点邻居加权组合。构造完成后限制算子R通常直接取P的转置这被称为对称配置保证粗网格系统对称正定。粗网格矩阵通过Galerkin公式构造% Galerkin粗化R取P的转置A_c完全由A_f、P决定 R P; A_c R * A_f * P;这行代码几乎是所有AMG实现的中心。它意味着粗网格矩阵完全从细网格矩阵投影而来不涉及任何几何坐标。也正因如此AMG可以把一个来源不明的稀疏矩阵当作输入而不需要知道它来自有限差分还是无网格法。我在做某跨平台模拟时遇到过一种常见误用有人直接把A_f去掉对角线当成粗化结果或者手工采样一行行拼A_c导致收敛率崩坏。正确做法是保持R P并严格用A_c RAP。这个恒等式一旦破坏后面的V-cycle几乎必然发散。Classic_AMG_Demo这类演示包之所以收敛曲线漂亮也是因为严格保持了这个统一框架。V-cycle的递归结构可以简洁表示为function x vcycle(A, b, level, max_level) % A_c和P由外部粗化阶段预先构造 if level max_level x A \ b; % 最粗层直接求解 else x relax(A, b, x); % 前光滑 r b - A * x; % 计算残差 e_c vcycle(A_c, R * r, level 1, max_level); % 粗网格校正 x x P * e_c; % 插值回加 x relax(A, b, x); % 后光滑 end end前光滑负责压掉高频误差粗网格递归负责修正低频误差后光滑再清理一遍插值过程引入的高频误差。三者缺一不可。把这个递归展开到第4层就形成一条从细网格一路下沉再一路上升的路径看上去像一个V字V-cycle因此得名。这里有一个经常被忽略的细节前光滑和后光滑的轮数设置不必一样。经典设置是前光滑1次、后光滑1次。对对角占优强的矩阵可以只做后光滑对刚性较强的结构问题往往前1后2的效果更好。参数要跟着矩阵性质走而不是一套配置打天下。3. 跑通Classic_AMG_Demo环境检查与最小操作路径3.1 运行Demo前要做的三件事每个标题里带“Demo”的MATLAB工程第一件事就是确认你的MATLAB版本能跑脚本还是必须跑函数。Classic_AMG_Demo这类经典演示通常由三个文件组成一个主脚本负责生成测试矩阵A和右端项b一个粗化函数实现RS粗化或类似策略一个V-cycle函数。三者的调用关系是主脚本调用粗化函数得到层级结构然后进入迭代循环反复调用V-cycle。我一般会把Demo目录放在一个不含中文和空格的路径下例如D:\work\amg_demo然后用addpath命令把目录加进MATLAB搜索路径。避免中文路径不是玄学老版本的MATLAB在中文路径下处理稀疏矩阵的IO偶尔会翻车。这个坑跟AMG没直接关系但足够让人卡在第一步。第二步是确认矩阵的存储类型。AMG的操作对象必然是稀疏矩阵。如果输入一个满矩阵粗化计算会发现所谓“强连接”关系几乎满篇都是生成的粗网格矩阵毫无稀疏性可言程序会慢到像死机。主脚本里通常写死了sparse的构造方式但如果读者想接入自己的矩阵必须先用sparse()显式转换。第三步是做一个最朴素的测试把AMG的迭代次数设为1直接看输出的x和MATLAB内置的A\b结果相差多少。这一步的意义是验证配置正确而不在乎收敛速度。如果连这一步都不一致后面的参数调优全是无用功。3.2 最小运行流程把一个Demo从下载变成收敛曲线以最常见的Demo启动方式为例操作步骤如下clear; close all; clc; % 1. 把Demo目录加入搜索路径 addpath(D:\work\amg_demo); % 2. 生成测试问题典型五点差分矩阵 n 100; A gallery(poisson, n); b ones(n*n, 1); % 3. 调用主求解入口返回解和残差历史 [x, hist] amg_solve(A, b, ... max_levels, 6, ... relax_factor, 0.8, ... pre_smooth, 1, ... post_smooth, 1, ... tol, 1e-8); % 4. 画残差收敛曲线 semilogy(hist.residual); xlabel(迭代步数); ylabel(相对残差); grid on;这段流程有两个值得解释的点。第一gallery(poisson, n)生成的是标准的n×n五点差分矩阵规模是10000阶。这个矩阵是AMG的理想测试样例对称正定对角占优又有足够的低频分量能直观看到V-cycle对低频的修正效果。第二amg_solve是Demo的主入口。如果原始代码里没有这样一个统一入口而是把主脚本的循环部分直接摊开你可以把它封装成函数。封装的好处是可以反复换矩阵做实验而不用改主脚本里的数据生成代码。几乎所有人第一次都会在max_levels这个参数上踩坑。max_levels设得太小比如3对10000阶矩阵还达不到最粗层V-cycle根本没走完收敛率跟单层迭代差不多。设得太大比如20粗网格矩阵小到直接求解毫无意义还白白多做几次矩阵乘法。对10000阶问题level6到8是常见区间对百万级问题通常设10到12。判断标准是最粗层矩阵的规模应当降到几个数量级以下而不是停留在几千个点。还有一个容易被忽略的点b的构造必须和A的维度一致。用ones构造右端项没问题但如果A是7000阶而你写了size(b,1)10000MATLAB会在迭代第一次计算rb-A*x时报维度错误。我见过因为这种低级错误排查半小时的案例报错信息极其难看因为错误已经发生在多层递归的某一层里。3.3 怎么确认“解对了”残差、真实误差和收敛率Demo跑完之后光看屏幕输出“converged”是不足以说明问题的。这里讲三个验证方法。第一看相对残差曲线是否呈指数下降。AMG的收敛曲线应当是接近直线的指数衰减如果曲线出现明显的平台期说明某个层次的松弛或插值配置有问题下面排错章节会详细讲。第二把AMG的解和直接法结果做一次绝对比较x_direct A \ b; % MATLAB内置直接法 err_inf norm(x - x_direct, inf); fprintf(最大误差: %.3e\n, err_inf);这里的误差会受松弛因子和迭代终止条件的双重影响但量级应当和tol处于同一水平附近。如果差了好几个量级问题几乎一定出在粗化或插值上而不是松弛参数。第三计算每步迭代实际达到的收敛率。对于对称正定矩阵AMG的收敛率稳定在0.05到0.3之间是健康的。如果收敛率超过0.9说明V-cycle基本没有在做事得回头查粗化。我这里的收敛率指相邻两步迭代残差范数的比值可以直接在hist.residual数组上计算diff(log(residual))观察每个周期的衰减斜率。4. 把AMG参数调到肉眼可见的快松弛因子、粗化阈值与层级数4.1 松弛因子0.8起步0.9以上开始不稳定V-cycle内部用得最多的松弛迭代是加权Jacobi或SOR。Classic_AMG_Demo通常默认前光滑1次、后光滑1次内部权重取0.8。这个配置对大多数模型问题已经能稳定收敛。松弛因子w的作用原理是w越接近1单次松弛对高频误差的压制越强但迭代格式从正定逐渐走向不稳定边界w太小比如0.5每轮光滑之后残差下降太少需要增加光滑轮数才能保证整体收敛率。工程上常见做法是在0.7到0.9之间先打底观察残差曲线。如果前几步衰减很明显但后期出现锯齿状波动说明松弛偏激进把w降到0.75再试。这里有个经验规律w每降0.1收敛率大约会恶化20%到50%但换来的是鲁棒性。对于对角占优较弱、条件数很大的矩阵我倾向于0.75起步而不是0.85起步省得后期排错。4.2 粗化阈值theta粗网格大小由它决定theta是RS粗化中最难调的参数因为它对矩阵性质非常敏感。0.25是经典默认值但这个默认值是从五点差分矩阵上调出来的不代表所有矩阵都适用。theta偏大比如0.5以上强连接关系的判定变严格能被选入粗层的C点数量减少粗网格矩阵偏小V-cycle变短每轮代价下降但插值精度下降整体收敛率恶化。theta偏小比如0.1几乎所有非对角元都被判为强连接C点数量膨胀粗网格矩阵和细网格矩阵规模差距不明显V-cycle变得臃肿等于没粗化。判断theta是否合适的工具不是理论公式而是在固定其他参数下观察两个指标C点占比和收敛率。C点占比通常在20%到40%是健康区间。如果超过50%说明粗化太温和低于10%则粗化太激进插值矩阵会非常稀疏低频误差修正能力不足。可以用下面这段代码扫描阈值% 扫描粗化阈值看每层C点占比和迭代步数 thetas [0.1 0.25 0.4 0.5]; for k 1:numel(thetas) [~, stats] amg_solve(A, b, theta, thetas(k), ... max_levels, 6, relax_factor, 0.8); fprintf(theta%.2f C占比%.2f%% 迭代步数%d\n, ... thetas(k), stats.C_ratio*100, stats.iter); end我在实际使用中会额外关心一个细节粗网格的层数是否足够。假如theta0.5的时候6层V-cycle实际只走完了4层就到最粗说明粗化速度太快。宁可theta从0.25向0.2方向微调也不要一步跳到0.1因为C点占比变化在0.2以下会发生突变收敛率曲线经常出现断崖。4.3 最后再动max_levels让V-cycle走到真正的底层max_levels的设置在某种程度上是theta的下游指标。合理的流程是先固定theta再根据粗层的规模来确定需要多少层。一行命令可以打印所有层的规模% 查看每层矩阵规模确认V-cycle走到底 for L 1:max_levels fprintf(第%d层: 阶数%d, 非零元%d\n, ... L, sizes(L), nnzs(L)); end如果最后一层的阶数还在500以上说明层级不够。如果到了第5层之后阶数已经小于50直接求解的成本很低再往下递归意义不大。对对称正定矩阵最粗层甚至可以设为1阶直接做一次除法对非对称矩阵我一般保留较多直接求解空间因为粗网格矩阵用不完全LU会比直接除法更便宜。参数表如下方便对照参数建议区间对收敛率影响常见误用relax_factor0.75 ~ 0.9每次松弛压高频的强度设到0.95以上易发散theta0.2 ~ 0.4决定C点占比与层级厚度0.1以下粗化失效max_levels6 ~ 12决定最粗层是否到底小矩阵设过大空耗pre/post_smooth1/1 ~ 1/2控制高低频修正平衡两边都设0必发散4.4 三个参数的联动效应一套可复用的调参顺序w、theta、max_levels三个参数不是独立的。theta调大导致V-cycle变短时需要增加前光滑轮数来补偿w调小导致单层松弛效率下降时可以通过增加后光滑轮数来兜底。我常用的调参顺序是固定theta0.25和max_levels6先调w等残差曲线稳定平滑再回去调theta看C点占比最后根据层级规模表决定max_levels。这套流程下来基本能覆盖大多数对称正定问题。遇到非对称问题比如对流占优的离散矩阵AMG经典配置很可能收敛不了那就要考虑换成AMG作为预条件子的GMRES而不是在松弛因子上继续耗。Classic_AMG_Demo能演示的边界就在这里——它教的是AMG的骨架不是所有矩阵的银弹。5. AMG常见坑与排查收敛翻车、索引错位与性能误判5.1 残差先快后蔫低频误差没被粗网格接住现象AMG前两三步残差下降很快然后曲线变成一条接近水平的线怎么加大迭代次数都没用。原因前光滑已经把高频分量压到底但低频分量需要粗网格校正而粗网格校正没有生效。这通常是两个原因之一theta设得过大导致C点不足或者max_levels太小导致V-cycle根本没递归到足够粗的层。解决先打印每层规模表确认最粗层阶数已经降到很小再看C点占比如果低于15%把theta往0.2方向回调。5.2 索引越界在递归内层爆出MATLAB从1开始参考伪码从0开始现象运行时报“Index exceeds matrix dimensions”而且报错位置出现在多层递归的内部不做深挖根本找不到是哪一行。原因AMG参考书里的伪码和Fortran实现大量使用从0开始的索引而MATLAB的数组索引从1开始。尤其是C/F划分里记录强连接邻居的下标数组照搬参考代码时忘了对所有下标加1导致在插值构造阶段访问了不存在的行。解决在所有强连接索引构造的地方统一加1再在粗网格构造结束后立即用size(A_c,1)检查确认A_c的行数等于粗层未知量总数。这个检查只需要一行能省下大量排查时间。5.3 残差收敛但解不对限制算子与插值算子失配现象残差范数明显下降最后也提示收敛但把解和A\b的结果对比发现误差大得离谱。原因限制算子R和插值算子P不匹配常见于手工改了插值规则但忘了同步R或者直接在外部对P做了截断却还在用旧R做粗网格投影。Galerkin粗化的基石是R P任何一边单独改动都会破坏粗网格矩阵的定义。解决每次修改插值逻辑后强制重新赋值R P再做一次A_c R * A * P。另外打印粗网格残差r_c R * r的范数观察它是否随层数下降如果某一层突然变大基本就是R和P脱节。5.4 稀疏化后性能反降粗网格密度失控现象矩阵明明从满矩阵转成了sparseAMG整体耗时反而比直接法还高甚至内存报错。原因粗化策略选得不对常见于theta取得太小C点占比膨胀每一层的非零元数量相比上一层下降不到一半。粗网格矩阵密度过大导致V-cycle单轮代价飙升整体迭代步数虽少但每一轮都慢。解决打印每一层的非零元数量如果相邻两层规模比低于60%说明粗化效率太差回退theta重跑。还要确认A本身在进入AMG之前不是被sparse(A)转成了满矩阵的稀疏存储却带着大量显式零元素这种情况要先做A sparse(A)再调A spones(A)清结构。5.5 做PCG预条件子时外层不收敛V-cycle和预条件目标脱节现象AMG单独跑迭代收敛正常PCG外层迭代次数却没有下降和不用预条件子几乎一样。原因AMG作为预条件子的核心评估标准是PCG外层迭代次数而不是AMG内部V-cycle的收敛残差。如果只保证AMG内部收敛而外层PCG无感说明AMG没有捕捉到矩阵的低频空间粗网格处理不够彻底。解决把PCG与外层AMG解耦验证。单独跑一次AMG求解记录收敛率是否在健康区间再单独跑无预条件的PCG确认问题本身的外层收敛率。两者都正常但合起来异常就在预条件初始化阶段打印预条件残差检查是不是AMG的初始猜测没有正确传递。常见问题出在把AMG的迭代解当作预条件算子的一部分而不是只把一次V-cycle的结果作为预条件作用。6. 把Demo改造成自己的求解器从函数封装到小矩阵回归测试我一般会用演示包做的事不是拿它跑一两次测试就结束而是把核心逻辑抽成一个带参数透传的接口然后塞进自己的求解流水线。第一步是把主脚本末尾的求解循环抽成函数保持默认参数和外部参数分离示意代码function [x, info] my_amg_solve(A, b, opts) % 默认参数 if nargin 3, opts struct(); end theta getfield_default(opts, theta, 0.25); max_levels getfield_default(opts, max_levels, 6); relax_factor getfield_default(opts, relax_factor, 0.8); tol getfield_default(opts, tol, 1e-8); % 构造层级、进入V-cycle迭代 levels setup_amg(A, theta, max_levels); [x, info.residual] amg_iterate(A, b, levels, ... relax_factor, tol); info.iter numel(info.residual); end这个封装最大的价值是让同一套AMG配置能在不同矩阵上做A/B测试而不是每次改一行参数就改一遍主脚本。我在本地保留着一个名为matrix_gallery的小测试集里面收集了从五点差分到对流扩散的七八种矩阵。每次改AMG代码我第一件事不是跑大型算例而是把这套小矩阵列表全部过一遍看收敛率有没有回退。这种回归式习惯帮我避开了不知多少次“改了一处插值弄坏了另一类矩阵”的事故。对于想把这个Demo用到真实项目里的读者我的建议是先不要急着调参追求最快收敛先把你自己的业务矩阵导入跑通稳定收敛再去优化参数和层级结构。AMG再精巧也只是求解链路里的一个环节它需要和你的矩阵性质、迭代外层和硬件环境一起工作。把你踩过的每一次残差平台期都记录下来配好一组参数组合存成配置这才是把演示代码变成工程工具的正路。希望帮到你。本文还有配套的精品资源点击获取
