简介面向消息传递算法与和积算法学习者的 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) * ∏_{s1..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} * u3u3是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不会溢出 endm在这一维上的每个切片里是一个标量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 endif 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 endlambda这一行是 damping也叫松弛。它的作用是把新算出来的消息和上一轮的消息做加权平均压制震荡。lambda 1就是完全不 dampinglambda 0.5是常见的起点lambda再小收敛更稳但更慢。有环图上不加 damping 能收敛是运气加了 damping 收敛是常态——这句话值得记下来。收敛判据用「消息的最大绝对变化」而不是「边际的变化」原因是边际可能已经稳定但消息还在缓慢漂移用消息判据能更早发现异常。阈值1e-6在概率域里是合理的如果你的消息量级特别小得改成相对变化delta / max(abs(msg))。5. 排查清单和积算法最常见的五种翻车5.1 现象消息在两个值之间反复横跳迭代到上限也不收敛原因几乎总是有环图加同步更新。flooding 调度下一个环上的消息会同时更新两股信息在环上对着跑撞上就震荡。加上环的长度如果是偶数震荡会非常规律看起来像周期解。解决办法按代价从低到高排先加 dampinglambda从 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 自由能不是凸的所以驻点可能不唯一算法收敛到哪个驻点取决于初值和调度。这就解释了为什么同一个图换个初始化结果会有细微差别。场景消息规则收敛性结果精度典型代价树结构和积两遍扫描必然收敛有限步结束精确边际O(边数 × 单因子规模)有环图 flooding和积迭代无保证可能震荡近似通常不错O(迭代数 × 边数 × 因子规模)有环图 damping 串行和积迭代大幅改善仍无理论保证近似通常更稳同上常数更大暴力枚举无不涉及精确指数级变量数上 20 就动不了这张表里最值得盯的是第二行和第三行的差距。很多工程实现之所以能用 loopy BP 解决实际问题靠的不是理论保证而是 damping 加串行调度把震荡压住了再加上问题本身的结构比较「像树」——环长较长、耦合较弱的时候Bethe 近似的误差很小。低密度校验码的译码就是这个逻辑校验矩阵对应的因子图环长通常在几十以上loopy BP 的表现逼近最大似然这也是它在工程里站住脚的原因。但反过来如果图的环很短、耦合很强loopy BP 就可能给出明显偏差甚至收敛到一个物理上没意义的解。这时候的补救方向有两个一是换更高级的近似比如给短环加上校正项的区域 BP代价是实现复杂度上升一个量级二是退回精确算法用变量消元或者联接树代价是指数级的树宽。选哪条路取决于你的图有多「不像树」这个判断没有通用公式只能拿小规模实例和暴力枚举对齐一次看误差能不能接受。我自己的习惯是不管最终要跑多大的图先在 8 到 10 个变量的小实例上把消息传递的结果和暴力枚举逐位对齐误差超过1e-10就不往下走。这一步花不了十分钟但能挡掉后面几天的怀疑人生。从那以后我每次写消息传递的代码都强制先跑这遍对齐再放到有环图上调 damping 和调度。希望帮到你。本文还有配套的精品资源点击获取
