简介：围绕因子图与和积算法的 Matlab 源码包，面向通信、机器学习与贝叶斯推断方向的学生和研究人员，适合课程设计、科研复现与算法自学。资源基于 SumProductLab 设计，覆盖消息传递算法（MPA）、LDPC 编解码、卡尔曼滤波、隐马尔可夫链（HMC）等典型应用，示例包括数独求解、布尔因子图推理、噪声模型、动态系统滤波等；通过阅读和运行源码，可观察变量节点与因子节点之间消息的生成、传递与更新过程，加深对概率图模型和近似推断的理解，且脚本中预置了多种基本因子节点（如与、或、等式约束、噪声等），并封装了边际化、模板等工具函数，便于组合出自定义模型。压缩包共 38 个文件，包含 31 个 .m 脚本、4 个 .mat 数据文件、2 个 PDF 参考文档及 1 个 license 说明，整体仅 826KB；脚本按节点定义、模型构建、示例演示等层次组织，MAT 数据文件提供预置模型与状态，PDF 文档对算法推导和使用方式做了系统说明，便于对照源码逐行学习。已有 2911 人学习下载，对于需要动手验证 MPA、实现 LDPC 译码或搭建自定义因子图模型的读者，这份源码提供了可直接修改的完整框架，可显著降低入门门槛。
1. 把指数级联合概率求和降成局部消息迭代：和积算法与因子图源码的落地玩法
做概率模型推理时，最朴素的思路是把所有变量的联合分布展开，再对目标变量做边缘化求和。但变量一过十五六个，全概率展开的计算量就指数爆炸，8GB内存的机器直接卡死。消息传递算法（Belief Propagation）及其核心形态和积算法（Sum-Product Algorithm），正是用来对付这种爆炸的：它把全局计算拆成因子图上的局部消息传递，每个节点只算相邻边上的信息，复杂度从指数级降到多项式级。这份MATLAB源码把因子图建模、双向消息更新、边缘分布计算和收敛判定全部封装成可独立调用的函数，适合做LDPC解码、马尔可夫随机场推断或贝叶斯网络推理的从业者直接改来用，也适合刚接触概率图模型的开发者照着代码理解BP到底在算什么。
2. 因子图与和积算法拆解：局部消息凭什么能还原全局边缘分布
2.1 因子图是怎样把联合分布拆成可计算单元的
因子图是一种二分图，节点分两类：变量节点（用圆圈表示）和因子节点（用方框表示）。每条边连接一个变量节点和一个因子节点，表示该变量出现在该因子对应的函数里。联合分布写成所有因子的乘积：
P(x1, x2, ..., xn) = (1/Z) * ∏ fa(xa)
这里的fa是因子函数，Z是归一化常数。每个因子只涉及一部分变量，变量之间的耦合关系通过共享变量节点体现。以某个简单贝叶斯网络P(A)P(B|A)P(C|B)为例，映射到因子图就是三个因子节点分别挂到A、AB、BC这三组变量上。这种拆法的核心价值在于：原本要对所有变量做一次性联合求和，现在可以沿图的边逐节点展开，每个局部只需要处理与自己相邻的因子。
在MATLAB里建模时，我习惯用结构体数组来组织这两类节点。变量节点存维度、取值域、邻居因子编号；因子节点存涉及变量列表和因子值表。这种组织方式比用稀疏矩阵存邻接关系更直观，调试时打印结构体就能看清楚图的拓扑。
2.2 消息更新规则：和积算法的两个方向
消息传递的核心是两类消息：变量节点发给因子节点的消息，表示变量对因子提供的"信念"；因子节点发给变量节点的消息，表示因子对变量取值倾向的"证据更新"。和积算法的更新规则只有两条：
变量节点到因子节点的消息，等于除目标因子外，所有邻居因子传来的消息逐元素相乘。因为变量节点不改变概率值，它只是把各方证据汇总。
因子节点到变量节点的消息，等于因子函数乘以所有相邻变量节点传来的消息，然后对除目标变量外的所有变量做求和（边缘化）。这里"和"是对离散变量求和，"积"是消息相乘，算法因此得名。
这两条规则交替迭代，直到消息不再变化。收敛后，任意变量的边缘分布等于它所有邻居因子传来的消息逐元素相乘，再归一化。
2.3 复杂度对比：从指数爆炸到线性增长
为了说清楚这个算法省了多少计算量，我用一个链式因子图举例：N个二元变量依次相连，每个因子只涉及相邻两个变量。全概率展开需要枚举2^N种状态组合然后逐个求和，N=20时已经是百万级；N=30就到了十亿级，MATLAB跑完得按小时计。而和积算法在链状图上的复杂度是O(N·d)，d是最大因子阶数，这里d=2，N=30只需要几十次乘加运算，毫秒级出结果。
在树状因子图（无环图）上，和积算法是精确推理，收敛后的边缘分布和暴力枚举完全一致。在带环图上，它退化为近似推理算法，消息可能不收敛或收敛到近似解，这一点放到后面的避坑章节细说。
3. MATLAB源码拆解：建图、消息迭代与边缘概率计算
3.1 数据结构设计：用结构体数组构建因子图
拿到这份源码，先看它怎么定义图结构。源码里用两个结构体数组分别管理变量节点和因子节点，这是最直接、也最符合MATLAB习惯的做法。
% 变量节点结构体 % id: 节点唯一标识 % dim: 变量取值个数（比如二元变量dim=2） % neighbors: 关联的因子节点id数组 % msg_to_factor: 发给各因子节点的消息（cell数组，每个元素是dim*1向量） var_node = struct('id', [], 'dim', [], 'neighbors', [], 'msg_to_factor', {}); % 因子节点结构体 % id: 节点唯一标识 % var_ids: 涉及的变量节点id数组（按顺序） % table: 因子值表，维度为[dim1, dim2, ...]，比如二元变量两两组合就是2x2矩阵 % msg_to_var: 发给各变量节点的消息（cell数组） fac_node = struct('id', [], 'var_ids', [], 'table', [], 'msg_to_var', {});这里有个关键设计：因子值表table的维度顺序必须和var_ids一一对应。比如var_ids = [3, 7]，table(1,2) 表示变量3取第一个值、变量7取第二个值时的因子函数值。这个对应关系在后续边缘化时极其重要，搞反了消息就全错。源码里用注释明确标了这一条，实际使用时我建议把var_ids和table的维度顺序打印出来核对一遍再跑迭代。
3.2 消息初始化与迭代更新：核心循环的实现
消息初始化的规则很简单：所有变量节点到因子节点的消息初始化为全1向量（长度等于变量维度）。因子节点到变量节点的消息初始化为全1向量（长度等于目标变量维度）。全1的意思是"还没有任何证据，所有取值同等可能"。
% 初始化所有消息为全1向量 % 变量节点侧 for i = 1:numel(var_nodes) k = numel(var_nodes(i).neighbors); var_nodes(i).msg_to_factor = cell(k, 1); for j = 1:k var_nodes(i).msg_to_factor{j} = ones(var_nodes(i).dim, 1); end end % 因子节点侧同理，msg_to_var初始化为全1接下来是迭代更新。和积算法的迭代框架是：依次更新每个因子节点到变量节点的消息，然后依次更新每个变量节点到因子节点的消息，反复循环直至收敛。源码把两步拆成两个函数，便于单独调试：
function msg = update_fac_to_var(fac, var_nodes) % fac: 因子节点结构体 % var_nodes: 变量节点结构体数组 % 返回值msg: 一个cell数组，每个元素是对应变量收到的因子消息 num_vars = numel(fac.var_ids); msg = cell(num_vars, 1); for v = 1:num_vars % 获取当前目标变量在因子表中的维度和索引 target_idx = v; other_idx = setdiff(1:num_vars, target_idx); % 构造乘积消息：因子表逐元素乘以其他变量节点传来的消息 % 先复制因子表，然后按其他变量的维度逐维扩展 prod_table = fac.table; for oi = other_idx var_id = fac.var_ids(oi); % 找到该变量在var_nodes中的位置 var_pos = find([var_nodes.id] == var_id); msg_vec = var_nodes(var_pos).msg_to_factor{fac_to_msg_idx(fac, var_id)}; % 按目标维度扩展并逐元素相乘 % 这里用reshape循环扩展，MATLAB的bsxfun或隐式扩展也可 expand_dim = ones(1, num_vars); expand_dim(oi) = numel(msg_vec); msg_reshaped = reshape(msg_vec, expand_dim); prod_table = prod_table .* msg_reshaped; end % 对除目标变量外的所有变量求和(边缘化) sum_dims = other_idx; tmp = sum(prod_table, sum_dims); % 对指定维度求和 % 去掉多余的单一维度 msg{v} = squeeze(tmp); % 归一化：消息是概率信念，归一化后数值更稳定 msg{v} = msg{v} / sum(msg{v}); end end这段代码的逻辑说明：prod_table逐元素乘上其他变量传来的消息，这一步对应"积"；sum对指定维度求和，对应"和"。squeeze去掉长度为1的维度，得到一个长度等于目标变量维度的向量。最后归一化，这行很关键——虽然归一化不改变边缘分布的相对比例，但能防止数值在下一次迭代中指数级膨胀或下溢。
变量节点到因子节点的更新更简单，就是把除了目标因子之外的所有邻居消息逐元素相乘：
function msg = update_var_to_fac(var, fac_nodes, factor_idx) % var: 变量节点结构体 % fac_nodes: 因子节点结构体数组 % factor_idx: 当前要发送消息的目标因子在var.neighbors中的下标 % 返回值msg: 长度为var.dim的列向量 msg = ones(var.dim, 1); nbrs = var.neighbors; for k = 1:numel(nbrs) if k == factor_idx continue; % 跳过目标因子 end % 找到该因子发来的消息 fac_pos = find([fac_nodes.id] == nbrs(k)); fac_msg = fac_nodes(fac_pos).msg_to_var{...}; % 细节见源码 msg = msg .* fac_msg; end msg = msg / sum(msg); % 同样做归一化 end逻辑说明：变量节点不改变信息，它只是把来自不同因子的证据汇总。跳过目标因子是为了避免消息循环自证——这是消息传递算法收敛的一个前提条件。
3.3 边缘分布计算：收敛后怎么提取结果
迭代收敛后，每个变量的边缘分布等于它所有邻居因子消息的逐元素乘积，再归一化：
function bel = compute_belief(var, fac_nodes) % var: 变量节点结构体 % 返回值bel: 归一化后的边缘概率向量 bel = ones(var.dim, 1); nbrs = var.neighbors; for k = 1:numel(nbrs) fac_pos = find([fac_nodes.id] == nbrs(k)); fac_msg = fac_nodes(fac_pos).msg_to_var{k}; bel = bel .* fac_msg; end bel = bel / sum(bel); end这一步本质上是把因子节点发来的所有证据汇总，得到该变量的后验信念。在树状图上，这个结果是精确的；在有环图上，它是环近似的结果，但通常足够好。
4. 参数配置：迭代次数、收敛阈值、消息调度怎么设才稳
4.1 调度策略：同步更新与异步更新的取舍
和积算法有两种消息调度方式。同步更新指每一轮迭代中，所有消息同时基于上一轮的消息计算，全部更新完后进入下一轮。异步更新（又称洪水序）指一个一个更新节点，每更新一个就立即用新值参与后续计算。
我做实验下来，同步更新在MATLAB里实现简单但收敛慢，往往需要多迭代50%到100%的轮数；异步更新收敛快，但结果依赖节点遍历顺序。源码默认用异步更新，遍历顺序按节点id从小到大。如果你发现结果不稳定，可以改成随机打乱遍历顺序，多做几次取平均，或者改用同步更新对比一下。对于链状和树状图，两种方式收敛结果一致；对于有环图，异步更新的收敛性通常更好，但可能收敛到不同的近似值。
4.2 迭代次数与收敛阈值：怎么判定消息"不再变化"
源码提供两种停止条件：固定迭代次数和阈值判停。固定次数适合结构稳定的场景，比如LDPC解码常常直接迭代50次或100次。阈值判停则需要计算相邻两轮迭代消息的差异量。
% 收敛判定：计算所有消息在两轮迭代间的最大绝对变化量 max_change = 0; for i = 1:numel(var_nodes) for j = 1:numel(var_nodes(i).msg_to_factor) diff = abs(var_nodes(i).msg_to_factor{j} - prev_msg{i}{j}); max_change = max(max_change, max(diff)); end end if max_change < tol % tol通常取1e-6或1e-8 converged = true; end阈值取多少合适？我一般这样设：如果只是做边缘分布估计，1e-4足够；如果要拿去做后续精确计算（比如解码比特判决），我会压到1e-6以下。注意，阈值压得太低会导致迭代次数显著增加，而收益往往微乎其微，因为消息收敛是渐进的，最后几百轮可能只改变了1e-7级别的量。
4.3 阻尼系数：消息震荡时的后悔药
在带环图上跑和积算法，最常见的情况是消息在两个取值之间反复横跳，残差一直降不下来。这时候不要盲目加大迭代次数，先检查是不是缺了阻尼。阻尼的数学形式很简单：
% 更新后的消息与旧消息做加权平均 new_msg = damping * new_msg + (1 - damping) * old_msg;damping的取值范围是(0,1)，通常取0.3到0.7之间。取0.5是通用起步值——让新旧消息各占一半。阻尼的作用是给消息更新加上惯性，防止大步修正造成震荡。实际经验是：树状图不需要阻尼，加阻尼反而会拖慢收敛；有环图通常需要阻尼，特别是当因子表的数值起伏较大时。代码里把damping设成了全局变量，调整时直接在配置段改，不用动核心循环。
提示：阻尼系数不是越大越好。超过0.7时，收敛速度会变得极慢，因为每次消息更新只挪动30%的步长。遇到震荡先试试0.5，不行再往上加。
5. 避坑指南：和积算法在MATLAB里最常见的五个翻车现场
5.1 消息矩阵维度对不上，乘不出来
现象：运行到消息更新时报错"Matrix dimensions must agree"，或者squeeze之后消息维度变成0×1。
原因：因子值表的维度顺序与var_ids不一致，或者边缘化后没有正确去掉单一维度。比如因子表是2×3（变量A取2个值、B取3个值），但var_ids写成了[B, A]，消息相乘时维度就对不上。
解决：在初始化因子节点后，立即写一段自检代码，打印每个因子的table尺寸和var_ids对应的维度，确认一致再进迭代。我在源码里加了一个verify_factor_table函数，就是干这个的。血的教训：建图时图省事少写了一对括号，排错排了两小时。
5.2 迭代不收敛，消息在两个值之间反复横跳
现象：残差曲线呈锯齿状，迭代50轮后max_change还在1e-2量级，没有下降趋势。
原因：带环图上的消息更新容易产生振荡，尤其是因子表的条件概率分布尖锐（接近0或1）时。
解决：先加阻尼系数0.5，多数情况能压住震荡。如果还不行，检查图的环结构：把因子图中的环找出来，看看是不是存在一条消息要绕环一周才能传回起点。此时可以改用max-product算法（最大积算法）做MAP推断，或者用环上BP的变种（如树重加权BP）。源码里提供了max-product版本的切换开关，切换后不用改数据结构。
5.3 消息被乘到错误的边上
现象：结果明显错误，但调试时每个函数单独跑都对，联合起来就错。
原因：因子节点发消息给变量时，用了数组下标而非节点id取消息，导致消息错位。特别是多个变量维度相同的时候，MATLAB不会报错，只会静默算错。
解决：所有跨节点取消息统一用id查找，不用位置索引。源码里已经把所有用下标的地方改成了id映射，这种改法多花点时间但值得。另一个习惯是给每个消息打上"来源节点id+目标节点id"的标签，出错时一眼定位。
5.4 数值下溢，概率全变成0
现象：迭代几轮后，消息向量里出现NaN或全部为0，归一化时报"Divide by zero"。
原因：概率值连乘多次之后下溢到双精度浮点的最小正数以下。比如20个0.1连乘的结果是1e-20，已经逼近MATLAB的realmin。
解决：两个方案。第一，每次消息更新后立即归一化，这是最省事的做法，但极端情况下仍需警惕；第二，改用log域消息传递：把乘法变成加法，求和时用log-sum-exp技巧。源码默认用方案一，我一般建议在因子个数超过15个时切换到log域。log域实现的核心代码如下：
% log域因子消息更新：避免乘法下溢 % log_prod = log(fac_table) + 其他变量log消息之和 % 边缘化后: log_msg = logsumexp(log_prod, sum_dims)5.5 因子图存在环，直接当成树状图跑
现象：BP结果和暴力枚举的正确答案对不上，差得不多但不是精确解。
原因：和积算法在树状图上精确，在有环图上只是近似。很多初用者拿到一张马尔可夫随机场图，看结构像树就直接跑，没检查是否有环。
解决：用MATLAB的graph对象构建邻接图，调用graph.dfs或直接检测边数是否等于节点数-1。如果边数大于节点数-1，必然有环。有环时根据场景选策略：边缘概率估计用BP图近似可接受；要求精确解时改用联结树算法（Junction Tree）或变量消元法。源码附了一个check_acyclic脚本，返回环上有哪些节点，至少能让你知道近似程度的来源。
6. 验证与进阶：三个技巧让源码真正变成你的工具
6.1 用暴力枚举做小规模对照验证
拿到源码后第一件事，不是直接跑大数据，而是先做一个极小规模的交叉验证：构造一个3到4个变量的因子图，用源码跑和积算法得到边缘分布，再写一个三重循环暴力枚举所有状态组合，直接从联合分布里边缘化出精确答案，两者对比。
% 暴力枚举对照示例: 3个二元变量, 因子表已知 % 期望输出: bel_brutal = [p(x1=0), p(x1=1)] bel_bp = run_bp_solver(example_graph); bel_brutal = brute_force_marginal(example_graph, 1); disp(max(abs(bel_bp - bel_brutal))); % 期望输出接近1e-15如果两者差超过1e-10，说明实现有bug，按第5章逐条排查。这个习惯能帮你把"看起来能跑"变成"确实算对"。我每次拿到新的消息传递实现，都会先跑这一遍，十分钟能省一天排错时间。
6.2 把因子换成校验约束，代码秒变LDPC解码雏形
源码的因子节点并不限定于概率表，它可以是任意非负函数。把因子换成校验方程的指示函数（满足校验取1，不满足取0），变量消息当信道的先验概率，和积算法就变成了置信传播解码器。这是这份源码最值得玩的一个扩展。改动点只有两处：一是因子值表从概率表换成校验约束的0/1表；二是因子归一化得谨慎处理，因为约束因子不满足时消息可能变成全0。我常用它跑一个7,4汉明码的软判决解码，效果很直观。
6.3 容易被忽略的细节：消息归一化的时机
消息归一化可以在两个位置做：一是在消息更新函数的末尾归一化，二是在计算边缘分布时统一归一化。很多人图省事只在最后归一化，中间消息逐轮膨胀或缩小。正确的做法是每次消息更新后立即归一化，原因有两个：防止中间量溢出或下溢；加快收敛速度。归一化改变了消息的绝对尺度但不改变相对比例，等价于每一步都重新缩放，实际收敛性能比不归一化好得多——在带环图上可以快两到三倍。
从那以后，我每次搭因子图模型，都强制走一遍先暴力枚举验证、再跑小规模对照、最后才上真实数据的流程。这套三层保险，已经帮我挡住了好几轮逻辑错位和维度搞反的翻车事故。希望这份源码和这些经验对你有用。
本文还有配套的精品资源，点击获取
