简介：面向学习消息传递算法、和积算法与因子图建模的读者，这份MATLAB源码包提供从理论到实现的完整练习材料，适合通信、信号处理、概率图模型及机器学习方向的学生与工程人员参考。压缩包共38个文件，约826KB，以31个.m脚本和函数为主体，配合4个.mat数据文件、2份PDF参考文档及1个许可说明，兼顾代码运行、样例数据与算法查阅。已有2911人学习，说明其在相关课程与自学场景中具有一定参考价值。代码围绕因子图上的和积消息传递展开，包含节点函数、边缘概率计算与图结构组织等基础模块，并给出LDPC编解码、卡尔曼滤波、隐马尔可夫链、数独求解以及与或、噪声或、证据推理等示例，便于读者对照运行、修改参数并观察消息更新过程。两份参考文档可辅助理解算法推导与实现细节，适合作为课程实验、算法复现和进一步扩展到置信传播应用的起点。
1. 从指数枚举到传小纸条：消息传递算法凭什么算得动
把一个看起来必须枚举指数多个组合才能算出来的边缘概率，拆成一张二部图上反复传递的几条小向量，就能一分不差地算出来——这是消息传递算法与和积算法最反直觉的地方。它的主线是在因子图上做纯局部运算：变量节点只做乘法，因子节点只做「乘完再求和」，迭代若干轮后，每个变量上攒出来的那一条小向量就是它的边缘分布。更反直觉的是，真正让它在工程里站住脚的恰恰是有环的图，理论保证全部失效，实测却好得离谱。理解了这套机制，再去看那份 matlab 源码，你看到的就不再是一堆嵌套循环，而是一张图在内存里的投影，每条边对应一个数组槽位，每次迭代就是把这些槽位按拓扑顺序读一遍、算一遍、写回去。
2. 因子图：为什么把联合分布拆成二部图就能算得动
2.1 联合分布里那个不肯消失的指数
先算一笔账。假设有 n 个离散变量，每个变量取 q 种值，联合分布 p(x₁,…,xₙ) 就要用 qⁿ 个数才能描述完整。想求其中一个变量的边缘分布 p(x₁)，必须把其余 n−1 个变量全部边缘化掉，求和次数是 qⁿ⁻¹ 量级。取 q=2、n=30，那就是五亿多次加法，而且前提是你已经把那张 2³⁰ 项的联合表塞进了内存——用双精度存，光这张表就要 8 GB 上下。
现实里的概率模型从来不是任意分布。图像里相邻像素才相关，卷积码里只有滑动窗口内的比特互相约束，隐马尔可夫模型里当前状态只依赖前一状态。这种局部性让联合分布可以写成若干个只依赖少量变量的函数的乘积：
p(x) = (1/Z) · ∏ fₐ(xₐ)
其中每个 xₐ 是全体变量 x 的一个很小的子集，fₐ 就是所谓的因子，Z 是归一化常数。这个式子本身就是因子图的定义式：把每个变量画成圆形节点，把每个因子画成方形节点，某个变量出现在某个因子里就连一条边——二部图就这么来了。
关键转变在于，一旦分布写成乘积形式，就再没有什么「联合表」需要被存储，模型的大小从 qⁿ 掉到了「因子个数 × 每个因子的表大小」。这跟把一张稠密矩阵换成稀疏矩阵是同一类操作，只不过稀疏的是概率依赖关系本身。
2.2 分配律才是真正的主角
乘积形式本身并不自动带来计算上的好处，真正的杠杆是分配律：求和号可以推进因子内部。拿最小的例子看，三个二值变量，联合分布
p(x₁,x₂,x₃) ∝ fₐ(x₁,x₂) · f_b(x₂,x₃)
求 p(x₁)，直接做法是对 x₂ 和 x₃ 做双重求和，共 2² = 4 项相乘再相加：
p(x₁) = (1/Z) ∑_{x₂} ∑_{x₃} fₐ(x₁,x₂)·f_b(x₂,x₃)
而 f_b 里根本没有 x₁，可以把它从内层求和里挪出来，先单独对 x₃ 求和：
μ_{b→x₂}(x₂) = ∑_{x₃} f_b(x₂,x₃) p(x₁) = (1/Z) ∑_{x₂} fₐ(x₁,x₂) · μ_{b→x₂}(x₂)
这就是一步消息传递。μ_{b→x₂} 是因子 b 发给变量 x₂ 的消息，它是对 x₂ 的一个函数——x₂ 有两种取值，所以它是两个数。然后 fₐ 收到这条消息，乘进去再对 x₂ 求和，得到 x₁ 上的边缘。
这几个变量看不出效率差异，但把它排成一条链，结论就很刺眼了：n 个变量组成的链，直接求边缘需要 qⁿ⁻¹ 次运算，按消息传递的做法只需要大约 n·q² 次。n=30、q=2 时，前者是五亿，后者是一百多次。七位数量级的差距，全部来自「把对多个变量同时求和拆成逐个变量求和」。能拆的前提只有一个：被边缘化的那个变量，只出现在相邻的少数几个因子里。图的结构决定了哪些求和能拆、按什么顺序拆。
顺手说一句，这就是变量消除算法的复杂度来源。因子图上的消息传递并不是另一套东西，它把变量消除过程中产生的中间量缓存下来复用——每求一个边缘都从头消一遍太浪费，消一次、把中间结果沿着边传下去，所有边缘就都拿到了。
2.3 变量节点、因子节点和一条边的三重身份
很多资料一上来就讲和积算法的两条递推式，读者容易卡在一个地方：为什么消息要分两个方向、为什么节点要分两种形状。答案藏在「一条边到底代表什么」这件事里，它有三种身份，理解了这三种身份，后面所有公式都是自然推论。
对概率模型来说，因子是一条局部约束或者一个局部似然，变量是待推断的量，边表示「这个变量是这个因子的第几个参数」。一个变量可以同时是多个因子的参数，这条边就是它参与某个局部函数的凭证。
对算法来说，节点是运算单元，边是消息通道。变量节点负责聚合，因子节点负责边缘化，消息只沿边走，绝不跨越节点。
对 matlab 源码来说，边是一个数组下标。每条边上有两个方向的向量，变量节点持有一个邻居边的下标列表，因子节点也持有同样的列表，整个图在内存里就是两张邻接表和两个消息矩阵。这个映射关系是读源码时最重要的一根线。
顺便把三种常见的图区分清楚，它们经常被混着用，但在实现层面差别很大。
| 图类型 | 节点 | 边 | 局部函数的归属 | 实现上的特点 |
|---|---|---|---|---|
| 贝叶斯网络 | 全是变量 | 有向 | 条件概率表挂在「子节点 + 父节点」上 | 一个变量对应一张大表，父子关系耦合在一起 |
| 马尔可夫随机场 | 全是变量 | 无向 | 势函数挂在团上 | 团越大表越膨胀，高阶团难处理 |
| 因子图 | 变量 + 因子两类 | 无向 | 因子节点显式承载 | 每个因子是独立函数对象，变量只剩标号 |
表格最后一行是因子图真正的价值。在贝叶斯网络里，一个变量如果有好几个父节点，关于它的条件概率表是一个整体，「它如何参与各个局部函数」这件事被压在一个节点里没法拆；因子图把每个局部函数拆成一个独立节点，变量退化成纯粹的参数位。落到代码上，因子可以用函数句柄或者数组实现，加一个新的因子类型不需要动任何已有代码，这跟面向对象里的策略模式是一个思路。
3. 和积算法：一条消息为什么是「乘完再求和」的函数
3.1 两条递推式，一条做乘，一条做边缘化
消息传递的规则只有两条，区别只在于「谁在发」。
变量节点 x 发给因子 f 的消息：
μ_{x→f}(x) = ∏_{g ∈ N(x){f}} μ_{g→x}(x)
因子节点 f 发给变量 x 的消息：
μ_{f→x}(x) = ∑_{x_N(f){x}} f(x_f) · ∏_{y ∈ N(f){x}} μ_{y→f}(y)
第一条说的是：变量节点把自己收到的所有消息乘起来，但是要把「正要发往的那个因子」传过来的那条排除掉。第二条说的是：因子节点把因子函数和其余所有入边消息相乘，然后对除目标变量以外的所有变量求和。
排除目标边这条规则容易被当成技术细节略过，实际上它是整套算法能工作的前提，业界管它叫外信息原则。设想不排除会发生什么：因子 f 上一轮发给 x 的消息，这一轮又被 x 原样乘回去，f 再把它当成「外部证据」收进来，同一条信息在环路上来回反射、自我放大。指数级增长的虚假置信度就是这么来的。排除自己刚收到的部分，才能保证每条消息里携带的都是「来自这一侧子树的新信息」。
至于因子节点那边的求和，本质就是一次边缘化——把因子定义的联合函数投影到单个变量上。所以和积算法的名字非常直白：求和（sum）与乘积（product）。变量节点做乘积，因子节点做乘积再求和。
3.2 消息不是一个数，是一张表
初学者最容易在实现阶段卡住的地方，是把消息当成标量。不是的。消息是对某个变量的函数，变量取 q 种值，消息就是一个长度 q 的向量。
一个二值变量上，消息长这样：(0.7, 0.3)。如果建模的是 QPSK 符号或者一个四值的状态变量，消息就是四个数。所以「传一条消息」在内存里是搬一个 q 维列向量，整张图跑一轮的开销是「边数 × q」，这两个数决定了源码能不能跑在大规模问题上。
理解这一点之后，两条递推式的代码形态就非常清楚了。变量节点那一步，是若干个 q 维向量逐元素相乘；因子节点那一步，是把因子表（一个 d 维数组）和第 j 条以外的入消息做外积式相乘，再沿除第 j 维以外的所有维度求和，最后得到一个 q 维向量。
真正会踩坑的地方是维度顺序。因子表 f 是 d 维数组，它的第 j 维必须与第 j 条边上变量的取值顺序严格对应，一旦顺序错位，代码照样跑、不报错、结果全错。这种错误在 matlab 里尤其隐蔽，因为数组索引错误往往不会导致越界，只是悄悄算出另一个函数。
3.3 树上为什么精确：消息只来自一侧
和积算法在树结构上是精确的，这不是近似。理解这一点不需要复杂的证明，只要抓住一条图论事实：树上任意一条边，把它切断，图正好分成两半，任意从出发点 x 到其他变量的路径都必须经过 f，绝无第二条路。
于是可以归纳。基础情形是叶子：叶子因子只连一个变量，它发出的消息就是自己对那个变量的求和，显然等于「这半张图里所有变量边缘化之后，对 x 的贡献」。归纳步骤是：假设离 x 更远的那些边都已经发出正确消息，那么因子 f 收到的每一条入消息就精确等于对应子树对 y 的贡献，f 把它们和因子函数乘起来、对其余变量求和，就精确得到「不含 x 那半张图对 x 的贡献」。
因为路径唯一，信息不会被重复计算两次；因为每半张图只从一条边流进来，最后把所有入消息乘起来，得到的就是整个模型在 x 上的边缘，只差一个归一化常数。
这也直接解释了树上的收敛性质：从叶子开始往内推一轮、再往外推一轮，两遍扫描之后每条边上的消息都到位了，再多迭代也不会变。所以如果你在树上跑和积算法，发现迭代了五十轮还没收敛，那一定是代码有问题，不是算法有问题。这个性质是最好的单元测试，后面会用到。
3.4 把求和换成取最大：只差一个算子
和积算法求的是边缘概率，如果目标变成最大后验（MAP），只需把因子节点的求和换成取最大：
μ_{f→x}(x) = max_{x_N(f){x}} f(x_f) · ∏ μ_{y→f}(y)
其他一切不变，这叫最大积算法。再取对数，乘法变成加法，最大值不变。如果进一步把 log-sum-exp 用「只取最小项」近似，就得到通信领域用得极广的最小和算法——它牺牲了少量精度，换来的是彻底摆脱超越函数。
这几种变体之间的差别，在代码上就是几行的事，但带来的性能差异是量级级别的。下面这段展示了同一个校验因子在两种算子下的实现，能直观看出精度换速度的取舍：
% inLLR: 入边的对数似然比向量，长度为 d % k : 本次要计算发往第 k 条边的消息 % 和积算法的双曲正切形式：理论上精确 function out = spCheck(inLLR, k) others = inLLR([1:k-1, k+1:end]); out = 2 * atanh(prod(tanh(others / 2))); % tanh 在大值时饱和，需裁剪 end % 最小和算法：把 log-sum-exp 换成取最小绝对值 function out = msCheck(inLLR, k) others = inLLR([1:k-1, k+1:end]); sgn = prod(sign(others)); % 符号取所有入边的乘积 out = sgn * min(abs(others)); % 幅度取绝对值最小的一条 endspCheck 里的 atanh(tanh(x/2)) 是二值变量下和积公式的经典化简，它成立的前提是变量只有两个取值、因子是偶校验。msCheck 则完全省掉了 tanh 和 atanh 两个超越函数，每次消息更新只剩比较和乘法。常见的经验值是精度损失在零点几个 dB 量级，而译码吞吐能翻倍——这就是为什么几乎所有硬件译码器都选最小和而不是和积。
4. matlab 源码：消息在内存里长什么样、迭代怎么转
4.1 数据结构：给每条边开一个槽位
读这类源码，第一眼不要看迭代循环，先看数据结构。因为循环里所有的下标运算都建立在数据结构之上，结构看不懂，循环就是一团乱麻。
最直观但最慢的写法是用三维数组按 (变量, 因子) 索引，问题在于绝大多数位置是空的，而且每次取消息都要判断「这两个节点之间到底有没有边」。工程上通行的做法是先建一张边表，把图的拓扑一次性固定下来，之后所有消息都以边编号为下标。
% ---- 边表：把图的拓扑一次性固定下来 ---- % adjV{i} : 第 i 个变量节点相连的因子节点编号 % e2v / e2f : 每条边两端的变量节点与因子节点 nE = sum(cellfun(@numel, adjV)); % 总边数 e2v = zeros(nE, 1); e2f = zeros(nE, 1); v2e = cell(nVar, 1); f2e = cell(nF, 1); cnt = 0; for i = 1:nVar for a = adjV{i}(:)' cnt = cnt + 1; e2v(cnt) = i; e2f(cnt) = a; v2e{i}(end+1) = cnt; % 变量节点持有的邻边列表 f2e{a}(end+1) = cnt; % 因子节点持有的邻边列表 end end % ---- 消息容器：每条边、每个方向各一个 q 维列向量 ---- q = 2; % 变量取值个数，二值场景 msgs.v2f = ones(q, nE) / q; % 第 e 列 = 边 e 上 变量→因子 的消息 msgs.f2v = ones(q, nE) / q; % 第 e 列 = 边 e 上 因子→变量 的消息几个参数值得说清楚。nVar 是变量节点数，nF 是因子节点数，nE 是边数，三者满足握手定理：nE 同时等于所有变量节点度数之和、也等于所有因子节点度数之和。q 是变量取值个数，二值变量取 2，这条不变量会直接影响后面所有数组的维度。
把消息按「边」而不是按「(变量, 因子) 对」索引，好处是初始化一次、之后零查找。代价是需要多维护 v2e 和 f2e 两张邻接表，但这两张表只在建图时算一次，之后的迭代里全是纯下标访问。这就是用一点预处理换掉每次迭代里的线性搜索，规模一大，差距非常明显。
4.2 因子节点的通用更新：一次外积乘、一次降维求和
因子节点的更新是整套代码里最需要小心的地方，因为它同时涉及数组维度和变量顺序。下面这个函数是通用的，任何因子只要给出取值表和维度信息就能用，不需要为每种因子单独写一套。
function out = spFactorMsg(fTab, dims, inMsgs, k) % fTab : 因子的取值表，元素个数为 prod(dims) % dims : 各邻接变量的取值个数，长度为 d，顺序必须与因子内部变量顺序一致 % inMsgs : 1×d 元胞，每个元素是该变量→本因子的消息（q 维列向量） % k : 本次计算发往第 k 个变量的消息 prod_ = reshape(fTab, dims); % 因子表摊成 d 维数组 for j = 1:numel(dims) if j ~= k % 把第 j 条入消息 reshape 成只在该维非单的 d 维数组，再逐元素乘 prod_ = prod_ .* reshape(inMsgs{j}, ... [ones(1, j-1), dims(j), ones(1, numel(dims)-j)]); end end % 对除第 k 维以外的所有维求和，得到关于第 k 个变量的函数 out = squeeze(sum(prod_, setdiff(1:numel(dims), k))); out = out(:); % 统一成列向量，防止维度塌陷 end这段代码对应第二条递推式，一步一步对得上。reshape(fTab, dims) 把因子表变成 d 维张量；循环里跳过了 j==k，这就是「排除目标变量」那条规则；reshape 后面的 [ones(1,j-1), dims(j), ones(1,...)] 是在构造一个只在第 j 维非单的形状，这样广播乘法正好等价于把消息沿该维度展开成外积；最后的 sum 带 setdiff 参数，等价于对除第 k 维以外的所有维做边缘化。
有一点必须强调：dims 的顺序不能随便给，它必须和因子内部定义的变量顺序严格一致。比如校验因子 f(x₁,x₃,x₅) 写成表的时候是按 (x₁,x₃,x₅) 的字典序展开的，那传进来的 dims 就必须是 [q₁,q₃,q₅]，传成 [q₃,q₁,q₅] 不报错，只是算出来的东西没有意义。这是这类代码里最常见、也最难查的一类错误。
还有 squeeze 之后的维度塌陷问题。当 d=2 时，squeeze 会得到一个向量，没问题；但当某条消息的维度恰好是 1（退化情况）时，squeeze 可能把本该保留的维度也压掉。最后那句 out = out(:) 就是兜底，把所有情况统一成列向量，避免下游按列索引时错位。
4.3 调度：flooding 与串行，为什么后者往往收敛更快
消息传递有两种更新顺序。第一种是 flooding：把这一轮所有变量节点的消息全算完、写下去，再统一算所有因子节点的消息，一轮结束所有边同步更新。第二种是串行或者分层：按节点或按边逐个更新，算出一条消息立刻用上，不等整轮结束。
代码结构上差别很小，都是那两层循环换一下位置：
for it = 1:maxIter % ---- 第一段：变量节点聚合并发消息 ---- for i = 1:nVar es = v2e{i}; for k = 1:numel(es) others = es([1:k-1, k+1:end]); m = ones(q, 1); for e = others m = m .* msgs.f2v(:, e); % 除目标边以外的入消息相乘 end msgs.v2f(:, es(k)) = m / sum(m); % 每步归一化，压住数值漂移 end % 顺手算这个变量当前的边缘，用于判决和提前终止 b = ones(q, 1); for e = v2e{i} b = b .* msgs.f2v(:, e); end belief{i} = b / sum(b); end % ---- 第二段：因子节点做边缘化并发消息 ---- for a = 1:nF es = f2e{a}; for k = 1:numel(es) in = cell(1, numel(es)); for j = 1:numel(es) if j ~= k in{j} = msgs.v2f(:, es(j)); end end m = spFactorMsg(fTab{a}, dims{a}, in, k); msgs.f2v(:, es(k)) = m / sum(m); end end % ---- 终止判据：信念稳定或约束全部满足 ---- if converged(belief, it), break; end end为什么串行往往更快？因为在 flooding 里，一条消息从图的一端传到另一端，至少要经过「与图的直径成正比」轮数，每一轮都要等所有节点同步。串行更新时，同一轮扫描内部后半段节点就已经用上了前半段的新消息，信息的实际传播距离接近翻倍。工程上常见的观察是，达到同样误码性能所需迭代数能减少三成到五成，代价是并行度下降——硬件上这两种调度的取舍完全是另一套账，但纯软件实现几乎没有理由选 flooding。
belief 那一段值得单独留意。边缘信念是所有入边消息的乘积再归一化，它和发出去的消息用的是同一批数据，顺手算一下几乎不增加开销，却给提前终止和调试提供了抓手。判断算法有没有在正常工作的最直接办法，就是盯住 belief 里最大分量的走势：如果它在两个值之间规律地跳，说明图上有环在振荡。
4.4 数值：归一化、对数域与提前终止
概率连乘的宿命是下溢。一张有几百条边的图，每条消息乘 0.5，几百次之后就是双精度也扛不住，直接变 0；一旦变 0，后面所有消息全部归零，边缘变成 0/0，输出 NaN。源码里那两行除以 sum 的归一化就是在处理这个问题：因为因子图只把分布定义到常数倍，每条消息乘任意正常数都不影响最终结果，所以可以放心地每步归一化到和为 1。
归一化解决的是漂移，解决不了动态范围的问题。更彻底的办法是转对数域，把消息存成对数似然比：
function s = logsumexp2(a) % 二值变量的 log(exp(a1) + exp(a2))，先减最大值防止 exp 溢出 m = max(a); s = m + log(sum(exp(a - m))); end减最大值这一步不是可选优化，是必需的。直接算 exp(a) 时，只要 LLR 超过 700 左右就会溢出成 Inf，而实际迭代中 LLR 轻松破千。转对数域之后，乘法变加法、归一化变减去均值、求和变 logsumexp，数值稳定性会好很多，代价是每次因子更新多两三个超越函数调用——这个代价在 matlab 里不小，最小和算法之所以流行，很大程度上就是在省这笔开销。
提前终止一般有两个判据。一是信念变化量：本轮 belief 和上轮 belief 的最大差异小于阈值就停，简单但对环上的缓慢振荡不敏感。二是约束满足：把每个变量的硬判决代进所有因子，看是否全部满足，这个判据在编码类问题里特别有效，因为满足全部校验方程就意味着找到了一个合法码字，再迭代也不会有改进。实践里两个一起用，先看约束满足，再看信念稳定性。
5. 避坑：五个让消息传递结果全错的常见问题
5.1 现象：环上反复振荡，永远不收敛
现象很典型：跑了几十轮，某个变量的 belief 在两个值之间规律跳动，振幅不衰减，错误位置每轮都在变。把图打印出来一看，存在长度为 4 或 6 的短环。
原因在于树假设被破坏。3.3 节那个归纳证明的关键前提是「任意两点之间路径唯一」，有环之后，一条信息绕一圈回到原点，被当成新证据再一次乘进去，形成正反馈。环越短、耦合越强，振荡越明显。
解决办法有几种，按代价从低到高：限制最大迭代次数并取历史最优解，这是最省事的兜底；加阻尼，把新消息和旧消息做凸组合 m ← (1−λ)m_old + λ m_new，λ 取 0.3 到 0.7 之间，实测对抑制振荡非常有效；改调度，把 flooding 换成串行，缩短反馈路径；如果结构允许，做因子图变换，把短环消掉。我的经验是先加阻尼，如果还不行，那就是图结构本身不适合直接跑 BP，该换算法了。
5.2 现象：消息全部变成 0，最后输出 NaN
现象是跑到某一轮之后，所有消息突然归零，belief 计算出现 0/0，结果全是 NaN。加打印一看，前几轮还正常，某一轮之后集体塌陷。
原因是概率连乘下溢。双精度的最小正规数在 1e−308 量级，一条长链上的消息每步乘一个小于 1 的数，几百步就到底了。
解决办法就是 4.4 节说的两条：每一步之后归一化到和为 1，把常数因子提走；或者干脆整体转到对数域，用 LLR 存储消息。要注意的是，即使做了归一化，如果某条消息的某个分量本身就是极小值，归一化也救不回来——这时候要看是不是因子表里有零项导致某条路径被彻底排除，那属于模型问题不是数值问题。
5.3 现象：边缘概率之和不为 1，判决阈值失效
现象是最后拿到的 belief 加起来不等于 1，或者做硬判决时 0.5 这个阈值不再好使，需要手动找一个偏移量才能对齐。
原因是因子图只把联合分布定义到常数倍。所有递推式里的等号其实都是正比于，缺的那个 Z 从头到尾都没被算过。如果中间某一步忘了归一化，误差会一路累积下去。
解决方式很直接：每个变量节点的 belief 在使用前必须归一化，别指望它自己收敛到和为 1。这件事在做硬判决时尤其重要，因为二值情况下判决规则就是比两个分量的大小，尺度无所谓，但一旦你要拿 belief 当软信息喂给下一级（比如级联的译码器），尺度就必须统一。
5.4 现象：图结构完全正确，结果就是全错
这种最气人：用暴力枚举算出来的精确边缘做对比，发现差得离谱，但检查图结构、因子表、迭代公式，全都没问题。
原因通常是消息槽位没有对齐。具体说，v2e{i} 的第 k 个元素和 f2e{a} 的第 j 个元素，必须指向同一条边，而且这条边在两端各自的消息数组中要落在同一列。如果建边表的时候变量端和因子端各数各的，编号对不上，消息就会从一条边发出去、从另一条边被读走，整套流程看起来完全正常。
解决办法只有一个：显式建边表，并且加一条断言检查——对每条边 e，e2v(e) 必须出现在 v2e{e2v(e)} 里、e2f(e) 必须出现在 f2e{e2f(e)} 里，两边指向的边编号必须一致。这几行检查代码在调试阶段的价值，比后面所有的打印都高。
5.5 现象：和积算法跑出 Inf，最小和算法却正常
现象是同一张图、同一组数据，和积版本迭代几轮之后消息变成 Inf 或 −Inf，换最小和版本一切正常。
原因是 tanh 饱和。3.4 节那段 spCheck 里，当入边 LLR 的绝对值超过 20 左右，tanh(x/2) 在双精度下就已经等于 1，几个 1 相乘还是 1，atanh(1) 直接返回 Inf。这不是算法错，是浮点数的表达能力到顶了。
解决办法是在调用 atanh 之前把参数裁剪到一个安全范围，比如夹到 [−1+1e−12, 1−1e−12]；或者改用 log 域下的等价形式，用 log((1+t)/(1−t)) 展开避免直接算 atanh。更常见的做法是干脆切到最小和或者它的归一化变体，因为那条路径上根本没有超越函数，不存在饱和问题。性能上还更快，算是把坑绕过去顺便捡了个便宜。
6. 设计取舍：什么时候该信 loopy BP，怎么验证它没算错
把前面几章的东西收一收，会看到一个很别扭的局面：理论上精确的只有树，而工程里真正好用的全是有环的图。这个矛盾没法用「近似算法」四个字糊弄过去，得具体看取舍。
如果有环的图，BP 一般不会收敛到真实边缘，它收敛到的是所谓 Bethe 近似下的驻点。这个近似什么时候好？经验规律是：环长越长、耦合越弱、图越稀疏，效果越好。稀疏长环的图（典型的 LDPC 校验矩阵就是这样构造的）上，BP 的表现接近理论极限；而短环密集、耦合强的图上，它可能连不动点都找不到。这不是调参能解决的，是模型和算法匹配的问题。
| 场景 | 推荐算子 | 消息表示 | 主要理由 |
|---|---|---|---|
| 树结构图模型 | 和积 | 线性域或对数域均可 | 一次遍历即精确，无需迭代 |
| 稀疏长环（编码类） | 最小和 / 归一化最小和 | 对数域 LLR | 省掉超越函数，吞吐优先，精度损失可接受 |
| 短环密集（视觉/统计物理） | 和积 + 阻尼 + 串行调度 | 对数域 | 抑制振荡，用调度换收敛性 |
| 需要 MAP 而非边缘 | 最大积 | 对数域 | 求和换取最大，其余框架复用 |
验证手段比调参更重要，我一般按下面三层来查。第一层是树上的自检：随机生成一棵树，跑 BP，结果必须在一次双向扫描后就和暴力枚举完全一致，误差小于 1e−12。这一步能抓出绝大部分索引错位和维度顺序错误。第二层是小规模暴力对比：变量数不超过 15 时直接枚举所有 qⁿ 种组合，算出精确边缘，和 BP 结果逐项比较，有环的情况下允许有偏差，但偏差应该在合理范围内、且随环长增加而减小。第三层是不变量检查：每条边的消息归一化后和为 1、每个 belief 归一化后和为 1、每次迭代后残差单调不增（有环时允许波动但不应发散），这三条只要有一条破了，就说明实现有问题。
还有一个容易被忽略的验证角度：跑同一份数据、只改变变量编号，结果应该完全不变。因子图是对称的，任何与编号相关的差异都是 bug 的信号。这一条几乎不花成本，却能抓出一类非常隐蔽的问题。
说个自己的教训。早期写这类代码，我最喜欢盯着误码率曲线看，曲线一旦正常就认为对了，结果有一次在一张中等规模的图上，曲线看起来完全合理，但换一组随机种子结果就崩了，查了两天才发现是边表在两个端点上的编号没对齐——错误被图的对称性掩盖了。从那以后我每次改完消息容器的结构，都强制先跑一遍树上的精确性自检和编号置换不变性检查，两分钟的事，省掉的是两天。希望帮到你。
这份 matlab 源码把这些机制都落成了可运行的代码：从边表构建、通用因子更新、两种调度方式，到最小和近似和收敛判据，每一段都对应上面某个具体的设计决定。想把它跑起来、或者想改成自己的图模型，直接拿去对着改比从头写要快得多。
本文还有配套的精品资源，点击获取
