简介：消息传递算法与和积算法的 Matlab 源码包，面向通信、信号处理与概率图模型学习者。通过构造因子图实现 Sum-Product Algorithm，可直接用于消息传递算法（MPA）、LDPC 编解码、卡尔曼滤波、隐马尔可夫链等典型场景，适合算法入门与二次开发。资源共 38 个文件，以 31 个 .m 脚本为核心，涵盖 LDPC、Kalman Filter、Hidden Markov Chain、数独求解等完整示例，并配有 4 个 .mat 数据文件、2 份 PDF 参考文档和 1 个 license 说明，压缩包仅 826KB，轻量易下载。代码按模块组织，多种因子节点实现与模板可帮助理解因子图构造、消息传递和边缘概率计算逻辑，PDF 参考文档又补充了算法原理与节点定义；示例附有运行数据，便于直接复现结果。目前已有 2911 人学习下载，适合需要从公式到代码落地验证的学习者参考，尤其适合通信、概率图模型方向的初学者与工程师。

1. 消息传递算法实战：和积算法在因子图上到底怎么算边缘概率

拿到一个概率推理需求，如果只知道先验和一堆因子表，怎么求某个变量的边缘概率？暴力枚举在 10 多个二元变量上就跑不动了。这份 Matlab 源码包实现的是消息传递算法里的和积算法，核心思路是把联合概率拆成一张因子图，让变量节点和因子节点沿着边互相“发消息”，迭代若干轮后，边缘概率直接从某个节点的消息乘积里读出来。我最早是在调一个 LDPC 译码简化模拟时碰到它，后来发现贝叶斯网络、隐马尔可夫模型里的推理也能套进同一个框架。适合想彻底看懂因子图推理、并且希望直接在 Matlab 里改参数跑实验的开发者。

2. 核心机制拆解：因子图建模与和积算法的两条消息更新规则

2.1 因子图建模：从联合概率到局部因子

概率图模型里最常见的诉求是：已知 p(x1, x2, ..., xn)，求某个变量的边缘分布。直接对联合概率求和，计算量随变量数量指数增长。因子图的想法是把联合概率写成若干个局部因子的乘积：

p(x) = (1/Z) * ∏ f_a(x_a)

其中每个 x_a 是变量集合的一个子集，Z 是归一化常数。比如一个三变量链式图，联合概率可以写成 p(x1,x2,x3) ∝ fA(x1,x2) * fB(x2,x3)。因子图里用圆形节点表示变量，方形节点表示因子，变量与因子之间画边表示“这个因子依赖哪些变量”。信源编码里的校验关系、贝叶斯网络里的条件概率表、马尔可夫随机场里的势函数，都能化成这种形式。

这份源码里的build_factor_graph就是干这个事的。它接收两个主要参数：vars_dim是每个变量的取值个数，tables是因子结构的 cell 数组，每个元素包含两个字段：vars表示该因子涉及哪些变量编号，table是多维概率表。我从实际调试的角度提醒一句：因子表每个维度的顺序必须和你填入的vars顺序一致，否则消息计算出来的边缘概率会错得很隐蔽。

2.2 和积算法的两条消息更新规则

假设图里有一个变量节点 v 和一个因子节点 f。变量发给因子的消息记为 μ_{v→f}(x_v)，因子发给变量的消息记为 μ_{f→v}(x_v)。和积算法的全部内容，就是两类消息的交替更新。

第一类，变量节点发给因子节点的消息：它等于该变量所有其他邻居因子发来的消息逐项相乘。直观理解为“我听了其他所有因子的意见，综合出来一个对 x_v 各取值的置信度，现在告诉你因子 f”。注意这里会把因子 f 自己的消息排除掉，防止信息在同一轮迭代里被重复使用。

第二类，因子节点发给变量节点的消息：它把因子表里除目标变量外的所有变量求和（边缘化），同时乘上那些变量发给该因子的消息。这一步是“因子对环境做校验”：如果其他变量传给因子的置信度与因子表本身严重矛盾，乘出来的消息会把对应取值压得很低。

2.3 边缘概率计算与证据因子

迭代收敛后，任意变量 x_i 的边缘概率正比于所有相邻因子发给它的消息乘积：

p(x_i) ∝ ∏_{f∈N(i)} μ_{f→i}(x_i)

得到未归一化的向量后，除以总和即为边缘概率。如果你有观测证据，比如已知 x3 = 0，不需要改主循环——把它表示成一个单变量因子 table = [1, 0]，塞进因子列表就行。观测证据在因子图里就是普通因子，这是这套框架很顺手的地方。

我在跑这个源码包自带的三元链例子时，初始消息全部设为全 1，跑了 30 多轮后输出 x2 的边缘概率稳定在 [0.57, 0.43] 附近。与 brute_force_marginal.m 暴力枚举结果完全一致。树形图上和积算法一定收敛到精确解，这也是为什么它常被当作概率图推理的基准实现。

3. 源码结构解析：Matlab 实现里的数据结构、消息更新与迭代调度

3.1 文件目录与函数职责

这份源码包的目录组织得很直接，没有绕弯子。关键文件如下表：

sum_product.m是核心调度器，update_v2f.m和update_f2v.m是真正干活的函数。我把这套结构拎出来的原因很简单：如果你想改成 max-product 做 MAP 推断，只需要改动update_f2v.m；如果你想换调度策略，只需要改sum_product.m。职责分得清楚，二次开发成本就低。

3.2 核心源码：变量到因子消息的消息更新

先看变量到因子的实现。源码包里的消息统一约定为msgs_f2v{f}{v}，表示因子 f 发给变量 v 的消息，这样的数据结构在取“某个变量的全部入边消息”时非常方便。

function msgs_v2f = update_v2f(msgs_f2v, nbrs_v) % 约定: msgs_f2v{f}{v} 是因子 f 发给变量 v 的消息 % nbrs_v{v} 是变量 v 的相邻因子编号数组 for v = 1:numel(nbrs_v) neigh = nbrs_v{v}; if isempty(neigh) continue; end for k = 1:numel(neigh) f = neigh(k); acc = ones(size(msgs_f2v{neigh(1)}{v})); for j = 1:numel(neigh) if j ~= k acc = acc .* msgs_f2v{neigh(j)}{v}; end end msgs_v2f{v}{f} = acc / sum(acc(:)); % 归一化, 防止数值下溢 end end end

这里的逻辑就是 2.2 节讲的乘法规则：对变量 v 来说，除了目标因子 f 外，把所有其他邻居因子发来的消息逐元素相乘，得到一个新的消息向量。acc初始化为全 1 是关键——乘法单位元是 1，不是 0，我最初改成 0 导致过整条消息被清零的翻车现场。

归一化那行是血泪经验。消息的绝对大小没有意义，只要各取值之间的比例正确就行。每轮归一化能避免迭代几十轮后消息值小到接近 0 甚至下溢成 NaN。常见做法是除以sum(acc(:))，也就是让消息的和为 1，相当于把消息当作一个概率分布处理。

3.3 核心源码：因子到变量消息的消息更新

因子到变量的消息是计算量最大的部分，也最容易写错。教学版实现先看二维因子表的情况，多维推广由marginalize_table.m完成。

function msg = factor_to_var_msg_2d(tab, other_msg, dim) % tab 是 k1 x k2 的二维因子表 % other_msg 是除目标变量外另一个变量的消息, 列向量 % dim 指定 other_msg 对应 tab 的第几个维度(1 或 2) if dim == 1 msg = tab' * other_msg; else msg = tab * other_msg; end msg = msg(:) / sum(msg(:)); end

这段代码把“乘上邻居消息再求和”这两步合并成了一次矩阵乘法。以dim == 2为例，因子表 tab 的第 1 维是目标变量，第 2 维是邻居变量。tab * other_msg表示对邻居维求和，结果向量里第 i 个元素就是目标变量取第 i 个值时，因子表和其他消息共同给出的未归一化置信度。用矩阵乘法代替显式求和，是这里最值得抄的实现技巧。

源码包里的marginalize_table.m把同样逻辑推广到了任意维：先通过 reshape 让其他变量的消息沿着对应维度广播，再对非目标维度循环求和，最后 squeeze 去掉长度为 1 的维度。如果你要处理十多个变量的因子表，直接看marginalize_table.m比看二维版更接近实际需求。

3.4 迭代调度与收敛判断

sum_product.m里的主循环默认使用并行调度，也叫 flooding 调度：每轮先更新所有变量到因子的消息，再更新所有因子到变量的消息。另一种常见做法是串行调度，也就是每算出一条新消息立刻用于后续计算，收敛速度通常更快，但实现时要小心消息的更新顺序不能乱。

收敛判断我习惯用消息变化量：主循环里可以新增如下判断逻辑。

% 伪代码: 每一轮结束后比较新旧消息 delta = 0; for f = 1:numel(msgs_f2v_new) for v = 1:numel(msgs_f2v_new{f}) delta = max(delta, abs(msgs_f2v_new{f}{v} - msgs_f2v_old{f}{v})); end end if delta < tol break; end

这里比较的是消息本身而不是边缘概率，原因是消息是所有后续计算的基础，消息收敛了边缘概率自然收敛，反过来则不一定成立。tol默认 1e-8 对大多数离散概率模型足够，如果你在跑带连续变量离散化的模型，建议先设成 1e-6 观察几轮再收紧。

4. 运行与复现：从初始化环境到换用自定义因子表

4.1 环境要求与初始化

这份源码是纯 Matlab 实现，基础函数只用到了矩阵乘法、reshape、sum 这类标准运算，不需要额外工具箱。我最初直接用命令行验证过，Matlab R2018a 以上版本都能跑，更老的版本只要支持 cell 数组和结构体也没问题。没有 Matlab 环境的开发者可以用 GNU Octave 跑，大部分语法兼容，但要注意 Octave 对某些隐式展开的告警会比 Matlab 严格。

拿到源码包后第一步是加入路径：

cd /path/to/sum_product_source addpath(genpath(pwd))

然后在命令行直接执行main_demo就能看到跑通结果。genpath(pwd)会把当前目录下所有子目录都加进来，避免 models 目录里的示例模型文件加载失败。

4.2 跑通内置 demo 并对照结果

main_demo.m里构造的是一个三变量链式因子图，每个变量取值为 0 或 1，两个因子分别约束 x1-x2 和 x2-x3。核心调用代码如下：

% main_demo.m 核心片段 vars_dim = [2 2 2]; tables{1}.vars = [1 2]; tables{1}.table = [0.9 0.1; 0.4 0.6]; tables{2}.vars = [2 3]; tables{2}.table = [0.8 0.2; 0.3 0.7]; fg = build_factor_graph(vars_dim, tables); res = sum_product(fg, 'max_iter', 100, 'tol', 1e-8); for i = 1:3 fprintf('x%d 边缘概率: %.4f %.4f\n', i, res.marginals{i}); end

tables{1}.table的第一维对应变量 1，第二维对应变量 2。第 (i,j) 个元素表示 x1=i 且 x2=j 时的因子值。跑完会看到三个变量的边缘概率，总和为 1。为了验证结果不是瞎写的，源码包里的brute_force_marginal.m会把所有 2^3 种变量取值逐一枚举，算出精确边缘概率。我在第 2.3 节说的 [0.57, 0.43] 就是通过这个对拍流程确认的。

这里有个实用的验证习惯：每次改完模型，先用暴力枚举对照小规模结果，再放大变量规模。暴力枚举脚本保留在源码包里不是鸡肋，它是你改因子表时最重要的安全网。

4.3 把内置模型换成自己的因子表

换模型只需要改vars_dim和tables两个变量，主循环完全不动。步骤大致如下：

第一步，确定每个变量的取值个数。比如你想表达一个三分类变量加两个二分类变量，就设vars_dim = [3 2 2]。第二步，按因子涉及的范围填tables。每个因子的table维度必须等于对应vars里各变量的取值个数，顺序一一对应。第三步，关注归一化。你的因子表可以是不归一化的（比如从统计频次转过来的），和积算法会在迭代过程中自动处理比例关系。第四步，跑完检查每个边缘概率是否求和为 1，以及和暴力枚举结果是否一致。

需要特别留意的是因子表中 0 值的处理。比如某个因子表里存在为 0 的项，消息计算时会直接把对应取值概率清零。这在概率上没问题，但如果你后续要做平滑，最好在填表之前就把 0 替换成一个小正数，否则一些节点可能在迭代中失去所有消息导致 NaN。

5. 避坑与排查：消息震荡、归一化错误与有环图近似

5.1 消息震荡不收敛

现象：迭代 20 多轮后边缘概率一直在 [0.38, 0.62] 和 [0.62, 0.38] 之间来回跳，tol设多小都停不下来。我用带强耦合因子的模型第一次跑出这种现象时，第一反应是代码写错了。

原因：二元变量加并行调度天然容易产生振荡。两个因子互相强烈约束对方时，消息会在两个极端之间摆荡，这是 loopy BP 的经典问题，不是 bug。解决方式两条路：一是换成串行调度，也就是每计算一条消息立即用于同一轮后续计算；二是给消息加阻尼，把新的消息按比例和历史消息做加权平均。阻尼实现的常见做法是new_msg = alpha * new_msg + (1 - alpha) * old_msg，alpha 取 0.5 左右能很好地压住振荡。

5.2 边缘概率总和不为 1

现象：输出的三个边缘概率各自求和不是 1，比如 [0.8, 0.1] 加起来只有 0.9。这种问题通常伴随其他变量边缘概率偏大，整体比例是乱的。

原因：最常发生在计算边缘概率时漏掉最后一步归一化，或者是因子表维度顺序与vars声明不一致。比如你心里想着第 1 维是 x2，但源代码里第 1 维是 x1，消息算出来自然张冠李戴。排查时先打印size(tables{i}.table)，确认每一维长度和vars_dim(tables{i}.vars)完全对应，再检查marginalize_table的 squeeze 是否把长度为 1 的维度错误剔除。最后在所有消息计算完成后补上msg = msg / sum(msg)，问题基本就消失了。

5.3 有环图结果和暴力枚举对不上

现象：把模型改成三角形环后，和积算法输出的边缘概率与暴力枚举结果有肉眼可见的偏差，比如某个变量概率相差 0.1 左右，而且增大迭代次数也不能完全消除。

原因：和积算法在树形因子图上保证收敛到精确边缘概率，但在有环图上只保证收敛到一个近似解，这就是 loopy BP 的近似特性。它在很多实际问题中效果不错，但不要指望精确。如果对精度有硬要求，路线有两条：一是用团树算法做一次树分解，把环变成树再跑；二是直接改用 max-product 做 MAP 推断，求最可能的联合取值而不是边缘概率。我在这个源码包上跑三角形环模型时就踩过这个坑，差点把正确的update_f2v当 bug 重写了一遍。

5.4 连续变量离散化导致计算量爆炸

现象：把连续变量离散成 50 个取值后，两个变量的因子表就有 2500 个元素，三个变量直接变成 125000 个，循环迭代跑起来明显卡顿，内存占用也一路飙升。

原因：用离散概率表近似连续分布，维度一高就指数爆炸。解决思路分两层：模型层面，把连续变量按高斯因子处理，消息只传递均值和精度两个参数，因子到变量的消息变成高斯消息更新公式；工程层面，用粒子消息近似，也叫非参数信念传播，通过一组带权样本代表消息。如果只是做教学验证，建议先把离散 bin 压到 10 以下跑通逻辑，再考虑连续因子扩展。

6. 二次开发：把消息引擎改成 max-product 并接入观测证据

6.1 从和积到最大积：一行改动做 MAP 推断

边缘概率解决的是“每个变量取值分别是什么概率”，但实际工程里经常需要的是“最可能的联合取值组合是什么”，也就是 MAP 推断。实现上，把因子到变量消息更新里的求和操作换成求最大值即可，其他逻辑完全不动。

function msg = factor_to_var_msg_2d_max(tab, other_msg, dim) % tab 是 k1 x k2 二维因子表 % other_msg 是邻居变量的消息, 列向量 if dim == 1 msg = max(tab .* other_msg, [], 1)'; else msg = max(tab .* other_msg', [], 2); end msg = msg(:); end

和 3.3 节的和积版本对比，核心区别就是tab * other_msg变成了max(tab .* other_msg, [], dim)。矩阵乘法求和变成逐元素相乘后取最大值，语义从“边缘化求和”变成了“保留最可能路径”。注意这个版本不需要归一化，因为我们要的是最大值的取值位置而不是概率比例，归一化不会改变 argmax。拿到消息后，从图的一端回溯 argmax 就能得到完整的 MAP 推断结果。

6.2 快速验证与改造顺序

我每次改完这类消息函数都会强制走一遍三个步骤：先在刚才的三元链上对拍和积与最大积的结果，确认语义切换正确；再在一个只有一个环的模型上跑一遍看是否收敛；最后才丢到真实规模的数据上。其中第一步最省时间也最值得做，暴力枚举脚本就在手边，跑一次不到一毫秒。从那以后我每次碰消息传递类算法的改动，都会先在最小模型上做一次对拍，再放量测试，这个习惯帮我避开了至少三次隐蔽的维度顺序错误。希望帮到你。

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