AMG代数多重网格详解：V-cycle、粗化与MATLAB实践

简介：经典代数多重网格（AMG）的Matlab演示实现，面向数值计算方向的研究生、工程师及算法爱好者，帮助理解利用AMG求解大型稀疏线性系统的基本流程。代码模块化程度高，包含V循环迭代、强依赖矩阵构造、粗网格生成与插值等核心步骤，便于逐函数阅读。包内附带二维Poisson方程在均匀三角网格上的有限元测试数据，覆盖51×51与101×101两套规模，并给出几何多重网格（GMG）的对比案例，方便通过数值实验观察两类多重网格策略的收敛差异；同时注明了粗网格生成较慢、未包含经典AMG生成粗网格的第二遍处理等边界条件，降低二次开发困惑。整套资源有20个文件，以11个m源码为主，辅以4个mat数据、2张结果图和README说明，压缩包仅309KB，小体积易下载。已有1227人浏览学习，适合想从代码层面入门AMG的读者。

1. 代数多重网格是什么——一个Demo解决什么问题

代数多重网格（Algebraic Multigrid，AMG）是所有做大规模稀疏线性方程组求解的人绕不开的一类方法。和几何多重网格不同，AMG不需要网格的几何信息，只靠系数矩阵A本身来构造一层一层的粗网格算子，这意味着它可以直接当作一个“黑匣子”求解器接到有限元、有限差分或图像处理代码里。Classic_AMG_Demo就是用MATLAB演示经典AMG完整流程的最小工程：从输入矩阵A和右端项b开始，经过粗化、插值、松弛迭代三步，输出一个可用的解x和收敛曲线。适合刚接触多重网格、想在工程里评估AMG是否值得投入的从业者，也适合做预条件子时想理解V-cycle内部行为的同学。

2. 从光滑误差到粗网格校正：拆穿AMG的V-cycle

2.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阶问题，level=6到8是常见区间；对百万级问题，通常设10到12。判断标准是：最粗层矩阵的规模应当降到几个数量级以下，而不是停留在几千个点。

还有一个容易被忽略的点：b的构造必须和A的维度一致。用ones构造右端项没问题，但如果A是7000阶而你写了size(b,1)=10000，MATLAB会在迭代第一次计算r=b-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

我在实际使用中会额外关心一个细节：粗网格的层数是否足够。假如theta=0.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会比直接除法更便宜。参数表如下，方便对照：

4.4 三个参数的联动效应：一套可复用的调参顺序

w、theta、max_levels三个参数不是独立的。theta调大导致V-cycle变短时，需要增加前光滑轮数来补偿；w调小导致单层松弛效率下降时，可以通过增加后光滑轮数来兜底。我常用的调参顺序是：固定theta=0.25和max_levels=6，先调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 稀疏化后性能反降：粗网格密度失控

现象：矩阵明明从满矩阵转成了sparse，AMG整体耗时反而比直接法还高，甚至内存报错。原因：粗化策略选得不对，常见于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再精巧，也只是求解链路里的一个环节，它需要和你的矩阵性质、迭代外层和硬件环境一起工作。把你踩过的每一次残差平台期都记录下来，配好一组参数组合存成配置，这才是把演示代码变成工程工具的正路。希望帮到你。

本文还有配套的精品资源，点击获取