企业官网建设流程全解析

简介：面向消息传递算法与和积算法学习者的 MATLAB 源码包，围绕因子图建模展开，可用于理解置信传播的基本流程，并延伸到 MPA、LDPC 编解码、卡尔曼滤波与隐马尔可夫链等典型应用。压缩包共 38 个文件、约 826KB，主体为 31 个 .m 脚本，覆盖节点函数、消息更新、边缘概率计算与多个示例入口；另有 4 个 .mat 数据文件、2 份 PDF 参考文档和 1 个说明文本，便于对照公式、运行示例并复现实验。已有 2911 人学习，说明其在相关课程与课题中具有一定参考价值。读者可借助清晰的功能拆分理解因子图节点如何组织、和积消息如何迭代，并从 LDPC、HMM、卡尔曼滤波等案例中提取可复用的建模与调试思路，适合课程实验、算法入门和项目原型验证。

1. 一次乘法和一次加法，凭什么替掉指数级枚举

很多人第一次看到因子图上的消息传递算法，会觉得它像玄学：明明只是每条边上反复做「乘一乘、加一加」，凭什么能替代对整个联合分布的暴力枚举？更反直觉的是，在一棵没有环的因子图上，和积算法只跑两遍就给出精确边际，一遍从叶子往上、一遍从根往下，中间不含任何近似。一份能跑通的 MATLAB 源码，核心函数通常也就两个：一个负责因子节点的「乘完再求和」，一个负责变量节点的「除掉自己再连乘」。理解这套机制的钥匙不在公式本身，而在一个问题上：为什么局部消息的乘积，恰好等于全局边际。

2. 因子图：把联合分布拆成一张二部图

2.1 从联合分布到连乘结构

一切推断问题的起点都是一个联合分布，而几乎所有的概率图模型都在做同一件事——把联合分布写成一堆局部函数的乘积：

p(x1, ..., xn) = (1/Z) * ∏_{s=1..M} f_s(X_s), X_s ⊆ {x1, ..., xn}

这个式子是后面所有机制的底座。f_s是只依赖一小撮变量的局部函数，Z是归一化常数。之所以能这么拆，是因为真实模型里变量之间的耦合几乎总是局部的：一张图像里像素只和邻居强相关，一个序列里当前状态只依赖前一时刻，一个校验方程只约束少数几个比特。全局的复杂性是局部耦合叠出来的，不是凭空长出来的。

这里有个容易被忽略的点：因子分解不唯一。同一个联合分布可以画成粒度很粗的因子图（一个因子吃下十个变量），也可以画成粒度很细的因子图（十个二元因子串起来）。粒度越细，单个因子的表越小，但消息传递的轮数越多；粒度越粗，单步计算越贵，但图更矮。常见的实现思路是尽量用「低元因子」——二元的够用就别写三元的，因为二元因子的局部表是一个矩阵，而矩阵乘法的性能在 MATLAB 里比多维数组的循环收缩高一个数量级。这个取舍在后面第 6 章还会再算一次账。

2.2 两类节点的语义：变量问「我是谁」，因子答「约束是什么」

因子图是一张二部图，节点分两类。变量节点x_i对应一个随机变量，因子节点f_s对应一个局部函数；当且仅当x_i ∈ X_s时，两者之间连一条边。边只存在于异类节点之间，同类节点之间永远不连——这是二部图的定义，也是消息传递能成立的结构保证：信息只能在「变量」和「约束」之间来回倒手，不会出现变量之间直接串话。

从语义上看，两类节点在回答不同的问题。变量节点收集所有邻居因子对它的「意见」，然后把自己的信念广播回去；因子节点则拿住一条硬约束，看看在其他变量已知的情况下，这条约束对本变量的取值有什么说法。这和贝叶斯网络里 CPT 的读法不同：有向图里 head-to-head 结构经常让人误判独立性，而因子图把每个 CPT 直接变成一个因子节点，谁依赖谁一目了然。这也是为什么大量推断代码宁愿多画一层节点，也要统一到因子图上——结构清晰度换来的调试时间是值得的。

提示：因子图不是贝叶斯网络的替代品，而是它的规范化形式。任何有向图模型都能无损转成因子图，反过来不一定。写代码时先统一到因子图，再谈算法，能省掉很多「这个箭头到底指哪边」的争论。

2.3 MATLAB 里的一张图：因子表 + 维度映射表

落到数据结构，一张因子图在 MATLAB 里需要三样东西：每个因子管哪些变量、每个变量有几个取值、每个因子的局部函数值长什么样。前两样是索引，第三样是数值。一份干净的实现通常长这样：

% 变量数 N，因子数 M N = 3; M = 3; % 每个变量的取值个数（二元变量就是 2） card = [2, 2, 2]; % scope{s} 记录因子 s 作用在哪些变量上，顺序即因子表的维度顺序 scope = {[1 2], [2 3], [3]}; % F{s} 是因子 s 的局部函数表，ndims(F{s}) == numel(scope{s}) F = cell(1, M); F{1} = [0.9 0.2; 0.1 0.8]; % 2x2，第1维对应 x1，第2维对应 x2 F{2} = [0.7 0.4; 0.3 0.6]; % 2x2，第1维对应 x2，第2维对应 x3 F{3} = [0.6; 0.4]; % 2x1，一元因子，等价于先验

scope这个元胞是整个实现的地基，它把「第 s 个因子的第 k 维」映射到「第几个变量」。后面所有的消息收缩都靠它来决定对哪一维求和、把消息乘到哪一维上。很多人写消息传递翻车，根子就在这一步：因子表用的是x2在前、x3在后，而消息数组按变量编号排，两边对不上，程序不报错，结果悄悄错。

F{3}这个一元因子值得单独说一句。它的作用等价于给x3一个先验，度数为 1 的因子节点在消息传递里是最简单的情形：它传给变量的消息就是因子表本身，不需要做任何求和。很多实现会偷懒，把先验直接乘到变量节点的初始消息里，省掉一个节点——这在树上没问题，在有环图上会改变消息的调度顺序，收敛行为可能完全不同。

3. 消息更新规则：乘什么、对谁求和

3.1 变量节点：为什么必须「除自己之外」

变量节点的规则一句话能说完：

μ_{x→f}(x) = ∏_{g ∈ N(x) \ {f}} μ_{g→x}(x)

翻译成人话：变量x要把自己的信念告诉因子f，但它不能说「我收到的所有消息的乘积」，而必须把f自己刚才传来的那条消息排除掉。这个「除自己之外」不是数学上的洁癖，而是防止信息回灌。

想象一棵以x为根的树，把邻居因子f剪掉之后，剩下的每个邻居g各自撑起一棵互不相交的子树。每棵子树对x的影响，全部浓缩在μ_{g→x}(x)这一条消息里；子树之间唯一的联系就是x本身，所以在给定x的条件下它们互相独立，联合贡献就是连乘。如果此时把f那条消息也乘进去，等于把f那一侧子树的信息又算了一遍，而f正等着用这条消息去算它那一侧——这就形成了自我引用。在有环图上，这种自我引用正是震荡的发源地。

实现上，「除自己之外」有两种写法。一种是每次算的时候现乘，跳过目标那条边；另一种是维护一个总乘积，算某条边时用总乘积除以那条边的消息。后者看着更快，但在概率域里除法会把数值精度问题放大，一条接近 0 的消息能把整条链带崩。常见的做法是老老实实跳过，度数不大的时候那点循环开销可以忽略。

3.2 因子节点：一次边缘化就是一次张量收缩

因子节点的规则稍微长一点：

μ_{f→x}(x) = ∑_{X_s \ {x}} f_s(X_s) · ∏_{y ∈ N(f) \ {x}} μ_{y→f}(y)

把它拆成三步看。第一步，把因子f_s之外的其他变量y各自传来的消息，按各自的维度乘进因子表——这一步是「把邻居的意见注入约束」。第二步，对所有不属于x的变量求和——这一步是「把不关心的变量边缘化掉」。第三步，剩下的就是一个只依赖x的向量，也就是因子对x的意见。

从张量的角度看，这是一次标准的收缩：f_s是一个d维数组，d是它的度数；消息是若干一维向量；把向量沿对应维度乘进去，再对除目标维以外的所有维度求和。这个操作的本质就是边缘化，只不过被边缘化的分布是「局部约束 × 邻居消息」，而不是全局联合分布。和积算法真正的魔法在这里：它把一次指数级复杂度的全局求和，拆成了M次低维求和，每次的规模只跟单个因子的度数有关。

拿第 2 章那个二元因子举例。f_A(x1, x2)是一个 2×2 矩阵，x1传进来的消息是u1 = [u1(1); u1(2)]。那么f_A发给x2的消息是：

μ_{A→x2}(x2) = ∑_{x1} f_A(x1, x2) · u1(x1)

写成矩阵形式就是μ = F{1}' * u1，一个矩阵乘向量。同理，f_B(x2, x3)发给x2的消息是μ = F{2} * u3（u3是x3传来的消息，维度对得上）。最后x2的边际是两条消息逐元素相乘再归一化。如果你手算一遍∑_{x1,x3} f_A · f_B，会发现结果和消息传递完全一致——这就是「树上精确」最朴素的一个实例。

注意：矩阵化只对二元因子成立。三元因子需要permute加多维求和，或者手工写三重循环。不要为了追求矩阵化硬把三元因子拆成两个二元因子，拆分会引入中间变量，在有环图上可能改变收敛性质。

3.3 归一化与对数域：把下溢这个黑匣子拆开

概率域的消息传递有一个致命弱点：连乘会下溢。一个度数为 300 的变量节点，假设每条消息的量级在 0.1 左右，乘积就是1e-300，刚好卡在双精度浮点数的realmin ≈ 2.2e-308边缘；度数再大一点，直接变成 0。之后所有涉及这个变量的消息全是 0，归一化变成0/0，输出 NaN。这个过程没有任何报错，只有结果莫名其妙地全错，是个典型的黑匣子。

标准的补救办法有两条。第一条是每步归一化：变量节点算完乘积后立刻除以总和，把消息拉回概率单纯形上。对于树结构，归一化只是乘了一个与x无关的常数，不影响最终边际；对于有环图，它会轻微改变消息的动态，但通常利大于弊。第二条是整体转到对数域：乘法变加法，求和变 log-sum-exp。对数域里消息是负实数，量级再大也只是个负数，不会下溢；代价是每次求和都要做一次log(sum(exp(...)))，比乘法贵。

log-sum-exp 的数值稳定实现有个固定套路：先减去当前维度上的最大值，再取指数求和，最后把最大值加回去。

function y = logsumexp_dim(X, d) % 沿第 d 维做 log(sum(exp(X)))，先减最大值防溢出 m = max(X, [], d); % 该维度上的最大值 Y = X - m; % 隐式扩展，需要 R2016b 及以上 y = m + log(sum(exp(Y), d)); % 指数项最大为 1，不会溢出 end

m在这一维上的每个切片里是一个标量，X - m靠 MATLAB 的隐式扩展广播下去。减完之后exp的最大值是exp(0) = 1，求和顶多等于该维长度，稳稳落在浮点范围内。真实代码里这一步通常和归一化合并：log 域归一化就是减去logsumexp沿着整个向量的结果，让最大项归零。

3.4 树上的两遍扫描：collect 与 distribute

树上不需要迭代。因为树没有环，消息不会绕回来污染自己，所以只要保证「算某条消息时，它依赖的消息已经算好」这个拓扑序，一次前向加一次后向就够了。

标准做法是随便选一个变量节点当根，做一次广度优先遍历，记下每个节点的父节点。然后按遍历序的逆序（后序，从叶子往根）算上行消息：每个节点把所有子节点传来的消息合成后发给父节点。再按遍历序（前序，从根往叶子）算下行消息：每个节点把父节点传来的消息分发给各个子节点，分发时同样要排除目标子节点那一条。

order = bfs_order(root); % 广度优先，父节点在前 for k = numel(order):-1:1 % 逆序 = 叶子到根 node = order(k); if parent(node) ~= 0 send_up(node, parent(node)); % 合成所有子消息，发往父节点 end end for k = 1:numel(order) % 正序 = 根到叶子 node = order(k); for each child c of node send_down(node, c); % 排除 c 那一条后下发 end end

两遍跑完，每个变量节点的所有入射消息都齐了，边际就是邻居消息的连乘再归一化。整个过程的消息总数是边上数量的两倍左右，没有任何浪费。这也是为什么树结构的推断可以做到和变量数线性相关——而暴力枚举是指数级。差距不在常数上，在复杂度阶数上。

有环图就没有这个好运了。环意味着消息会绕回自己，拓扑序不存在，只能迭代到（希望）收敛。这是第 6 章要展开的取舍。

4. MATLAB 源码拆解：消息容器、核心函数与调度

4.1 消息怎么存：边编号比 cell 嵌套更抗造

消息存在边上，不在节点上。最直觉的存法是二维元胞msg{i}{j}，表示变量i发给因子j的消息，查找方便，但有两个问题：一是稀疏图里大量空元胞浪费内存，二是遍历「所有边」的时候要写双重循环加判空，代码难看得像一锅粥。

更抗造的写法是边列表。把所有边编号，维护一个ne × 2的矩阵，每行[var_id, fac_id]，再用两个长度为ne的元胞数组分别存两个方向的消息。这样遍历边就是一次for e = 1:ne，没有任何判空。

% edges: ne x 2，每行 [变量编号, 因子编号] % msgV2F{e}: 第 e 条边上，变量 -> 因子的消息（列向量，长度 = 该变量取值数） % msgF2V{e}: 第 e 条边上，因子 -> 变量的消息 edges = [1 1; 2 1; 2 2; 3 2; 3 3]; ne = size(edges, 1); msgV2F = cell(ne, 1); msgF2V = cell(ne, 1); for e = 1:ne msgV2F{e} = ones(card(edges(e,1)), 1); % 初始化为均匀，等价于「无信息」 msgF2V{e} = ones(card(edges(e,1)), 1); end

初始化为均匀分布不是随手写的。均匀消息在概率域里对应「我什么都不知道」，在 log 域里对应全零；它保证了第一轮迭代时，因子节点收到的邻居意见是中性的，输出只由因子表本身决定。如果你用随机数初始化，在有环图上会引入额外的抖动，调试时很难判断不收敛到底是算法问题还是初值问题。

4.2 因子节点核心函数：multiply-and-sum 的矩阵化写法

这是整份源码里最值得逐行读的函数。它要处理任意元数的因子，所以不能写死成矩阵乘法，得用「沿维度乘消息 + 对非目标维求和」的通用形式。

function out = fac2var(F, inMsgs, target) % F : 因子局部函数表，ndims(F) 等于该因子的度数 d % inMsgs : 长度为 d 的元胞，inMsgs{k} 是第 k 个变量传来的消息（列向量） % target : 要输出给第几个变量（1..d） % out : 因子发给 target 变量的消息（列向量） T = F; d = ndims(F); % 第一步：把除 target 外的所有消息，沿各自维度乘进因子表 for k = 1:d if k == target continue; % 排除目标变量自己那条消息 end % 把列向量 reshape 成「第 k 维有效、其余维为 1」的形状，靠隐式扩展广播 shape = ones(1, d); shape(k) = numel(inMsgs{k}); T = T .* reshape(inMsgs{k}, shape); end % 第二步：对除 target 外的所有维度求和（边缘化） sumdims = setdiff(1:d, target); for dd = fliplr(sumdims) % 从大维度往小维度删，索引不会错位 T = sum(T, dd); end out = T(:); % 强制列向量，避免行/列方向混淆 end

几个细节值得停下来看。

第一，reshape(inMsgs{k}, shape)这一手是隐式扩展的关键。消息本来是n×1的列向量，要乘到第k维上，就得把它变成1×...×n×...×1的形状，让 MATLAB 自动广播到其他维度。shape的构造方式保证了只有第k个位置是n，其余全是 1。这个技巧在 R2016b 之后可以直接用，之前的版本得换成bsxfun。

第二，求和循环用fliplr(sumdims)从大维度往小维度删。原因是sum(T, dd)会把第dd维压缩成大小 1 但不消失，如果你从小往大删，删掉第 1 维之后原来的第 2 维变成了第 1 维，循环变量就对不上了。从大往小删，被删维度之前的索引不受影响。这是纯纯的血泪经验，第一次写这个函数的人十个里有八个在这里栽过。

第三，out = T(:)强制列向量。squeeze之后剩下一维的方向取决于原始维度顺序，如果target是第 1 维，结果可能是行向量；如果是第 3 维，可能是列向量。不统一方向，下一步变量节点做逐元素乘法时就会触发隐式扩展，得到一个矩阵而不是向量，错误极其隐蔽。

对于二元因子这个高频情形，可以直接走矩阵乘法，省掉上面所有 reshape 和循环：

% f(x1,x2) 是 2x2 矩阵，u1 是 x1 传来的消息 % 发给 x2 的消息：对 x1 求和 msg_to_x2 = F' * u1; % 等价于 sum_{x1} f(x1,x2) * u1(x1) % 发给 x1 的消息：对 x2 求和 msg_to_x1 = F * u2; % 等价于 sum_{x2} f(x1,x2) * u2(x2)

一行矩阵乘法顶掉十几行通用代码，速度差好几倍。工程上的常见做法是写一个fac2var做通用兜底，再对度数为 2 的情况加一个特判分支，两条路都留着。

4.3 变量节点与边际计算

变量节点比因子节点简单得多，核心就是「排除自己再连乘」，外加一步归一化。

function out = var2fac(inMsgs, target) % inMsgs : 该变量收到的来自各邻居因子的消息（元胞） % target : 要发给第几个邻居（排除它） out = ones(size(inMsgs{1})); for k = 1:numel(inMsgs) if k == target continue; end out = out .* inMsgs{k}; end s = sum(out); if s > 0 out = out / s; % 归一化，防止连乘下溢 end end

if s > 0这个判断不是防御性编程的洁癖。如果上游某条消息因为数值问题退化成全零，out就是全零，0/0直接产生 NaN，然后 NaN 会像传染病一样顺着边扩散到整张图。加一个判断，至少能让错误停留在局部，方便定位。

边际的计算就是不再排除任何一条消息：

function b = belief(inMsgs, card_i) % 变量节点的边际：所有入射消息连乘后归一化 b = ones(card_i, 1); for k = 1:numel(inMsgs) b = b .* inMsgs{k}; end b = b / sum(b); end

拿到边际之后，硬判决就是[~, xhat] = max(b)。如果你要的是最小误符号率，判决每个变量取边际最大的那个值；如果要的是整条序列的最大后验，那得在因子图上另做 max-product，那是另一套算法，和和积算法共享同一张图但消息规则不同——把求和换成取最大值，其他结构一模一样。

4.4 调度、damping 与收敛判据

有环图上没有拓扑序，只能迭代。最朴素的调度叫 flooding：每一轮里，先把所有因子到变量的消息更新一遍，再把所有变量到因子的消息更新一遍，两趟算一轮。

lambda = 0.5; % damping 系数，0 表示完全用旧值，1 表示完全用新值 tol = 1e-6; % 收敛阈值 maxIt = 100; for it = 1:maxIt oldF2V = msgF2V; % 第一趟：所有因子 -> 变量 for e = 1:ne s = edges(e,2); i = edges(e,1); % 收集该因子其他边的变量->因子消息，按 scope 顺序排好 raw = fac2var(F{s}, collect_in(s, e), local_pos(s, i)); msgF2V{e} = (1-lambda)*msgF2V{e} + lambda*raw; % damping end % 第二趟：所有变量 -> 因子 for e = 1:ne i = edges(e,1); s = edges(e,2); msgV2F{e} = var2fac(collect_in(i, e), local_pos(i, s)); end % 收敛判据：所有因子->变量消息的最大变化量 delta = 0; for e = 1:ne delta = max(delta, max(abs(msgF2V{e} - oldF2V{e}))); end if delta < tol break; end end

lambda这一行是 damping，也叫松弛。它的作用是把新算出来的消息和上一轮的消息做加权平均，压制震荡。lambda = 1就是完全不 damping，lambda = 0.5是常见的起点，lambda再小收敛更稳但更慢。有环图上不加 damping 能收敛是运气，加了 damping 收敛是常态——这句话值得记下来。

收敛判据用「消息的最大绝对变化」而不是「边际的变化」，原因是边际可能已经稳定但消息还在缓慢漂移，用消息判据能更早发现异常。阈值1e-6在概率域里是合理的，如果你的消息量级特别小，得改成相对变化delta / max(abs(msg))。

5. 排查清单：和积算法最常见的五种翻车

5.1 现象：消息在两个值之间反复横跳，迭代到上限也不收敛

原因几乎总是有环图加同步更新。flooding 调度下，一个环上的消息会同时更新，两股信息在环上对着跑，撞上就震荡。加上环的长度如果是偶数，震荡会非常规律，看起来像周期解。

解决办法按代价从低到高排：先加 damping，lambda从 0.5 试到 0.2；不行就换串行调度，按边的编号顺序一条一条更新，用刚算出来的新消息去算下一条，收敛速度通常比 flooding 快一倍；再不行就限制迭代次数，取最后若干轮边际的平均值作为输出。第三种做法听着不严谨，但在译码类应用里是常规操作，因为不收敛的轮次往往只影响少数几个变量。

5.2 现象：边际输出全是 0 或 NaN，但图结构检查过没问题

先查下溢。在var2fac里加一句打印，看看连乘前的消息量级。如果出现1e-200以下的数，基本可以确诊。度数为几十的变量节点在概率域里跑几十轮就可能触底。

解决办法是每步归一化，还不行就整体转对数域。转的时候注意因子节点也要改：fac2var里的.乘法变成加法，sum变成logsumexp_dim。只改一半是最常见的错误——变量节点用 log 域、因子节点还用概率域，两边的数值语义直接对不上，结果比不转还离谱。

5.3 现象：明明是一棵树，结果和暴力枚举对不上

树上的和积算法是精确的，对不上一定是实现错了。最高频的错误是消息方向搞反：var2fac里忘了排除target，把目标因子自己传来的消息也乘进去了。这个错误在有环图上表现为收敛慢，在树上表现为结果直接偏掉，而且偏得很微妙——大部分时候数值看着还挺合理。

第二个高频错误是下行消息（根到叶子）算的时候用错了上行消息的版本。上行阶段算出的消息是「子节点合成后发给父节点」，下行阶段父节点发给子节点时，要排除该子节点那一条，用父节点收到的其他消息重新合成。如果图省事直接复用了上行的那条消息，等于把子节点自己的信息又送回给它。

排查手段很土但有效：找一棵只有三个变量的小树，手算一遍，打印每一步的消息和手算结果对比。三变量树的中间结果就那么几个数，对不上立刻能看出来是哪一步。

5.4 现象：换个变量编号顺序，结果就变了

这说明代码里存在对变量顺序的隐式依赖。最常见的是因子表的维度顺序和scope的顺序不一致——F{1}是按(x2, x1)排的，scope{1}却写成[1 2]，于是fac2var把x1的消息乘到了第 2 维上。

这类 bug 的隐蔽性在于它不报错，维度大小恰好相等的时候连形状检查都过不了关。防御手段是强制约定：scope{s}的顺序就是F{s}的维度顺序，没有任何例外；写因子表的时候在注释里显式标出每一维对应哪个变量。更进一步，可以在构图阶段加一个断言，检查size(F{s})和card(scope{s})是否逐维相等。

assert(isequal(size(F{s}), card(scope{s})), ... '因子 %d 的表尺寸与 scope 声明的取值数不一致', s);

十几行断言换几小时的调试时间，这笔账怎么算都划算。

5.5 现象：结果正确，但几百个变量跑了一整夜

性能问题通常出在两个地方。一是因子节点用了通用实现，每个二元因子都在做reshape加sum，而这里本该是一行矩阵乘法。二元因子在图里占大头，这个优化能带来数倍的提速。

二是消息容器用了嵌套元胞，遍历时不断做msg{i}{j}的两级索引，MATLAB 的元胞索引本身就有开销。换成边列表加扁平元胞，配合预分配，能再省一截。如果图更大，可以考虑把同构的因子批量处理——所有 2×2 的因子堆成一个 3 维数组，一次pagefun或者一次reshape后的矩阵乘法全部算完，这是 MATLAB 里做批量小矩阵运算的常用套路。

还有一类隐性开销来自收敛判据。每轮都遍历所有边算最大变化量，边数上万的时候这个循环本身就不便宜。改成每 5 轮检查一次，或者只抽查一部分边，对最终结果几乎没有影响。

6. 设计取舍：树上的精确与环上的近似，差在哪一步

树上精确、环上近似，这个分野的根源在第 3.1 节那句话里：变量节点必须排除目标因子传来的消息，因为那条消息里已经包含了变量自己的信息。在树上，排除之后剩下的子树互不相交，连乘得到的确实是精确的边际贡献；在有环图上，排除一条边并不能切断所有路径，消息绕一圈还是会回来，于是「排除」变成了一个近似动作——你排除的是直接回路，排除不掉间接回路。

这个近似有个名字叫 Bethe 近似，它对应的目标函数是 Bethe 自由能。和积算法的不动点恰好是 Bethe 自由能的驻点。这个结论的解释力在于：它告诉你 loopy BP 不是「凑合用的启发式」，而是在一个明确定义的近似目标上做优化；同时它也告诉你，Bethe 自由能不是凸的，所以驻点可能不唯一，算法收敛到哪个驻点取决于初值和调度。这就解释了为什么同一个图换个初始化，结果会有细微差别。

这张表里最值得盯的是第二行和第三行的差距。很多工程实现之所以能用 loopy BP 解决实际问题，靠的不是理论保证，而是 damping 加串行调度把震荡压住了，再加上问题本身的结构比较「像树」——环长较长、耦合较弱的时候，Bethe 近似的误差很小。低密度校验码的译码就是这个逻辑：校验矩阵对应的因子图环长通常在几十以上，loopy BP 的表现逼近最大似然，这也是它在工程里站住脚的原因。

但反过来，如果图的环很短、耦合很强，loopy BP 就可能给出明显偏差，甚至收敛到一个物理上没意义的解。这时候的补救方向有两个：一是换更高级的近似，比如给短环加上校正项的区域 BP，代价是实现复杂度上升一个量级；二是退回精确算法，用变量消元或者联接树，代价是指数级的树宽。选哪条路取决于你的图有多「不像树」，这个判断没有通用公式，只能拿小规模实例和暴力枚举对齐一次，看误差能不能接受。

我自己的习惯是：不管最终要跑多大的图，先在 8 到 10 个变量的小实例上把消息传递的结果和暴力枚举逐位对齐，误差超过1e-10就不往下走。这一步花不了十分钟，但能挡掉后面几天的怀疑人生。从那以后我每次写消息传递的代码，都强制先跑这遍对齐，再放到有环图上调 damping 和调度。希望帮到你。

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