简介这是一份面向通信、信号处理与概率图模型方向学习者与研究者的MATLAB实践资源聚焦消息传递算法在因子图上的实现与应用。资源完整封装了和积算法Sum-Product Algorithm的核心框架与典型用例可直接支撑LDPC码编解码、隐马尔可夫链推理、卡尔曼滤波建模及数独求解等经典任务适合具备概率论与MATLAB基础的中高级学习者开展原理验证与算法拓展。压缩包共38个文件含31个功能模块化.m脚本如factor_graph_test、LDPC、Hidden_Markov_Chain等、4个预置测试数据.mat文件、2份PDF参考文档含算法原理与接口说明及1个许可说明txt整体仅826KB轻量易部署。目前已有2911人学习下载内容结构清晰从基础节点and_node、or_node、cp_node到复合模型sudoku_9x9、Kalman_Filter_Example再到系统级工具marginal、ls_marginal覆盖建模、消息调度、边缘计算全流程是理解图模型推理机制不可多得的可运行范例。1. 因子图上的消息传递不是黑匣子一份能跑通、能改参数、能看懂每条消息流向的 MATLAB 实战源码你有没有试过在论文里读到“在因子图上执行和积算法进行边缘概率推断”然后打开 MATLAB 写了三行sumprod就卡住不是报错是根本不知道中间那几十次消息更新到底传了什么、从哪来、到哪去——更别说把高斯混合模型或 LDPC 编码器的因子图手动画出来再映射成代码。这份源码就是为这种场景写的它不封装成一个run_inference()黑盒函数而是把因子图构建、变量节点/因子节点定义、消息初始化、迭代调度、消息更新公式、归一化逻辑全部拆成可打断、可打印、可单步调试的 MATLAB 脚本。它面向的是正在啃《Probabilistic Graphical Models》第 9 章、刚学完贝叶斯网络想过渡到因子图、或者正为课程设计实现 LDPC 译码器而抓狂的开发者。所有核心逻辑都在.m文件里没有 mex、不依赖工具箱仅需基础 MATLAB跑起来后你能用disp(msg_v2f{2}{3})直接看到第 2 个变量节点发给第 3 个因子节点的消息向量长什么样——这才是真正能帮你建立直觉的源码。2. 项目结构与技术栈MATLAB 原生实现无外部依赖全模块可读可断点这个项目不是一堆散装.m文件而是一个有明确分层意图的小型工程。它用最朴素的 MATLAB 数据结构cell 数组、struct、double 向量模拟图结构与消息状态避免了面向对象带来的抽象开销也绕开了 Symbolic Math Toolbox 等非标配组件。整个包解压后共 7 个文件核心是 4 个.m脚本 1 个示例配置脚本 2 个辅助函数。我把它按功能划分为三层图建模层定义拓扑、算法层消息更新规则、驱动层调度与可视化。这种分法不是为了炫技而是为了让你在调试时能快速定位——比如发现某次迭代后某个边缘概率异常你不需要翻遍整个文件直接去update_msg_v2f.m查变量节点发消息的逻辑就行。2.1 目录结构与文件职责说明factor_graph_sumprod/ ├── main_demo.m # 入口脚本构建示例图3变量2因子、设置初始势函数、调用主循环、打印结果 ├── build_factor_graph.m # 图建模层输入变量数、因子数、邻接关系输出 graph_struct含 nodes, edges, adj_v2f, adj_f2v ├── init_messages.m # 初始化层根据变量域大小、因子势函数生成初始消息全1向量或均匀分布 ├── update_msg_v2f.m # 核心算法层变量→因子消息更新∏_{k∈N(v)\f} msg_{k→v} × marginal_v ├── update_msg_f2v.m # 核心算法层因子→变量消息更新∑_{x_{-v}} [factor_pot(x) × ∏_{u∈N(f)\v} msg_{u→f}] ├── normalize_messages.m # 归一化层对每个消息向量做 L1 归一化避免数值下溢支持开关控制 └── utils_plot_beliefs.m # 可视化层绘制每次迭代后各变量的后验信念belief支持保存为 gif提示build_factor_graph.m是你二次开发的第一站。它不硬编码图结构而是接受adj_v2f {[1,2], [1,3], [2,3]}这样的 cell 输入——这意味着变量 1 连接因子 1 和 2变量 2 连接因子 1 和 3……这种表示法和《Factor Graphs and the Sum-Product Algorithm》原论文图 1 的编号完全对齐抄论文图时不用脑内转换。2.2 技术选型依据为什么坚持用基础 MATLAB 而非面向对象或 Simulink很多人会问都 2024 年了还手写因子图用 OOP 封装 Node 类不行吗我的血泪经验是在算法理解阶段封装是最大的认知负担。当你连msg_{f→v}(x_v) ∑_{x_{-v}} ψ_f(x) ∏_{u∈N(f)\v} msg_{u→f}(x_u)这个公式里的x_{-v}是怎么枚举的都没搞清时面对f.send_message_to(v)这种方法调用你得先跳进FactorNode类里找send_message_to的实现再确认ψ_f是存在self.potential还是self.table里……这已经偏离了算法本质。而本项目所有消息更新都写成独立函数输入明确msg_in,pot_f,domain_v,neighbors输出明确msg_out你可以把update_msg_f2v.m拖进编辑器加一行fprintf(f%d→v%d: before norm, sum%.6f\n, f_id, v_id, sum(msg_out));就能实时监控数值漂移。至于 Simulink它的图形化建模对理解消息流有帮助但一旦要改一条消息的归一化方式比如换成 max-normalization 做 MAP 推断你就得在图形界面里找半天 block 参数——而本项目改normalize_messages.m里一行msg msg / sum(msg);就够了。2.3 核心数据结构设计用 cell 数组承载异构性用 struct 统一接口因子图里变量节点和因子节点类型不同变量有 domain如{0,1}或1:5因子有 potential table维度随邻居数指数增长。如果强行用矩阵存很快就会内存爆炸或维度错乱。本项目用两个关键结构解耦graph.nodes.varstruct 数组每个元素含domain1×D 向量、prior1×D 向量先验分布graph.nodes.factorstruct 数组每个元素含scope变量 ID 列表如[1,3,4]、tablepotential 表size prod(domain_size(scope))消息则用嵌套 cellmsg_v2f{v_id}{f_id}存变量 v→因子 f 的消息1×|domain_v| 向量msg_f2v{f_id}{v_id}存因子 f→变量 v 的消息1×|domain_v| 向量。这种设计看似啰嗦但带来两个实际好处一是size(msg_v2f{2}{3})永远等于length(graph.nodes.var(2).domain)不会因因子连接多个变量而维度混淆二是你可以用cellfun(length, msg_v2f{2}, UniformOutput, false)一键检查变量 2 发出的所有消息是否长度一致——这是调试时救命的命令。3. 和积算法核心实现逐行解析消息更新公式看清每一步的数学含义与 MATLAB 映射和积算法的本质是让消息在变量节点V和因子节点F之间来回传递最终每个变量的 belief 是其所有相邻消息的乘积。但“乘积”和“求和”的具体操作位置极易混淆。本节带你逐行拆解update_msg_f2v.m和update_msg_v2f.m把公式里的每个符号对应到代码里的每一行、每一个索引、每一次 reshape。3.1 变量→因子消息更新update_msg_v2f.m的完整逻辑链该函数计算msg_{v→f}(x_v) b_v(x_v) / msg_{f→v}(x_v)其中b_v(x_v)是变量 v 的当前 belief即所有入边消息的乘积msg_{f→v}(x_v)是上一轮 f→v 的消息需避免自反馈。代码实现如下function msg_out update_msg_v2f(v_id, f_id, graph, msg_f2v, msg_v2f) % 获取变量 v 的 domain 和当前所有入边消息来自其他因子 domain_v graph.nodes.var(v_id).domain; D length(domain_v); % 构建 belief: ∏_{k ∈ N(v)\f} msg_{k→v}(x_v) % msg_f2v{k}{v_id} 是因子 k → 变量 v 的消息尺寸 1×D neighbors_f graph.adj_v2f{v_id}; % 所有连接变量 v 的因子 ID msg_in ones(1, D); % 初始化为 1乘法单位元 for k 1:length(neighbors_f) f_k neighbors_f(k); if f_k ~ f_id % 排除目标因子 f避免自反馈 if ~isempty(msg_f2v{f_k}{v_id}) msg_in msg_in .* msg_f2v{f_k}{v_id}; else warning(msg_f2v{%d}{%d} is empty, using uniform, f_k, v_id); msg_in msg_in .* (ones(1,D)/D); end end end % 应用变量先验如果定义了 if ~isempty(graph.nodes.var(v_id).prior) msg_in msg_in .* graph.nodes.var(v_id).prior; end % 除以 f→v 的旧消息若存在实现 belief / incoming if ~isempty(msg_f2v{f_id}{v_id}) % 防零除加极小值 msg_denom msg_f2v{f_id}{v_id} eps; msg_out msg_in ./ msg_denom; else msg_out msg_in; % 第一轮无旧消息直接输出 belief end end这段代码的关键在于msg_in是 belief 的近似但它不是真正的 marginal而是所有入边消息的乘积再乘先验。很多初学者误以为这里该调用sum其实变量节点本身不执行求和它只做乘法和除法。eps的加入不是可选项——在多次迭代后某些消息分量可能下溢到 0直接除零会炸掉整个迭代。我一般会把eps改成1e-12并注释“此处非防 NaN是保数值稳定性LDPC 译码中此值影响收敛轮数”。3.2 因子→变量消息更新update_msg_f2v.m的张量收缩实现这才是和积算法的“和”之所在。公式msg_{f→v}(x_v) ∑_{x_{-v}} [ψ_f(x) × ∏_{u∈N(f)\v} msg_{u→f}(x_u)]中x_{-v}表示除x_v外所有变量的联合取值。MATLAB 里没有自动张量收缩必须手动枚举。代码用ndgrid生成全组合再用sub2ind映射到 potential table 索引function msg_out update_msg_f2v(f_id, v_id, graph, msg_v2f, msg_f2v) scope graph.nodes.factor(f_id).scope; % 因子 f 连接的变量列表如 [1,2,4] K length(scope); % 因子邻居数 domain_sizes zeros(1,K); % 获取每个邻居变量的 domain size for k 1:K var_id scope(k); domain_sizes(k) length(graph.nodes.var(var_id).domain); end % 构建所有邻居的消息乘积∏_{u∈N(f)\v} msg_{u→f}(x_u) % msg_v2f{u}{f_id} 是变量 u → 因子 f 的消息1×|domain_u| msg_prod ones(domain_sizes); % K维数组尺寸为 [D1,D2,...,DK] for k 1:K u_id scope(k); if u_id ~ v_id % 排除目标变量 v msg_u msg_v2f{u_id}{f_id}; % 将 msg_u 广播到第 k 维需要 permute repmat dims ones(1,K); dims(k) domain_sizes(k); msg_u_reshaped reshape(msg_u, dims); msg_prod msg_prod .* msg_u_reshaped; end end % 加载因子势函数 ψ_f(x)尺寸应与 msg_prod 一致 pot_f graph.nodes.factor(f_id).table; % 已预处理为 [D1,D2,...,DK] 多维数组 % element-wise multiply: ψ_f(x) × ∏ msg weighted pot_f .* msg_prod; % 对除 v_id 外所有维度求和即 ∑_{x_{-v}} % 找出 v_id 在 scope 中的位置 pos_v find(scope v_id); sum_dims setdiff(1:K, pos_v); % 需要求和的维度 msg_out sum(weighted, sum_dims); % 输出为 (K-1)维不squeeze 后是 1×D_v msg_out squeeze(msg_out); % 降维确保是行向量 msg_out msg_out(:); % 强制为 1×D_v 行向量 end这段代码的玄学点在于sum(weighted, sum_dims)。假设scope[1,2,3]v_id2则pos_v2sum_dims[1,3]sum(..., [1,3])就是对第 1 和第 3 维求和留下第 2 维——这正是x_2的 domain。squeeze和(:)是防御性编程MATLAB 的sum在单维求和时可能返回列向量而后续归一化要求行向量。这个细节不处理normalize_messages.m里sum(msg_out)就会算错。3.3 归一化策略与数值稳定性为什么不能只用msg msg/sum(msg)normalize_messages.m提供两种模式l1默认msg msg / (sum(msg) eps)—— 保证sum(msg)1用于计算概率maxmsg msg / max(msg)—— 保证max(msg)1用于 MAP 推断找最大后验但重点不在模式选择而在何时归一化。本项目在每次update_msg_*返回后立即归一化见main_demo.m循环体内而非在所有消息更新完再统一归一。这是刻意为之如果等全部 20 条消息算完再归一某条消息因数值下溢变成[0,0,1e-200]sum就是0归一化后全NaN。而逐条归一至少保证每条消息自身数值范围可控。我一般还会在normalize_messages.m里加一行if any(isnan(msg) | isinf(msg)) error(NaN/Inf detected in message %d-%d at iter %d, v_id, f_id, iter); end——宁可早报错不传毒数据。4. 运行与调试从零启动、观察消息流、验证收敛性避开新手必踩的 5 个坑别急着改代码先让它跑起来并且跑得明白。main_demo.m自带一个经典示例两个二元变量x1,x2一个二元因子f1定义ψ1(x1,x2)I(x1x2)同或另一个因子f2定义ψ2(x2)I(x21)硬约束 x21。理论后验应为P(x11)1, P(x10)0。我们用这个例子练手再展开排查。4.1 三步启动法确保环境干净、路径正确、参数合理清空工作区与路径缓存在 MATLAB 命令行执行clear; clc; close all; rehash toolboxcache; % 强制刷新函数缓存避免旧版本残留 addpath(genpath(factor_graph_sumprod)); % 添加整个文件夹含子文件运行main_demo.m前先检查关键配置打开main_demo.m找到这几行MAX_ITER 20; CONVERGENCE_TOL 1e-6; NORMALIZE l1;MAX_ITER20足够小例子收敛但如果你接 LDPC建议先设5看前几轮行为CONVERGENCE_TOL不是越小越好1e-8在双精度下可能永远不满足因为sumprod本身有数值误差NORMALIZEl1必须与你的势函数类型匹配——如果ψ_f是 log-potential如高斯核就得先exp(pot_f)再用l1否则得换max。首次运行专注看终端输出main_demo.m末尾有fprintf打印每轮后x1的 belief。正常输出应类似Iter 1: P(x10) 0.3333, P(x11) 0.6667 Iter 2: P(x10) 0.2000, P(x11) 0.8000 ... Iter 6: P(x10) 0.0000, P(x11) 1.0000 → converged!如果第一轮就全是NaN问题出在势函数或 domain 定义如果一直震荡不收敛可能是图结构有环未被识别见避坑章节。4.2 消息流可视化用disp()和plot()看清每条消息的生死不要只信最终 belief。在main_demo.m的迭代循环里插入调试语句if iter 3 fprintf(\n--- Messages at iter 3 ---\n); fprintf(msg_v2f{1}{1} ); disp(msg_v2f{1}{1}); fprintf(msg_f2v{1}{1} ); disp(msg_f2v{1}{1}); fprintf(msg_f2v{1}{2} ); disp(msg_f2v{1}{2}); end你会看到类似msg_v2f{1}{1} [0.5000 0.5000] msg_f2v{1}{1} [0.3333 0.6667] msg_f2v{1}{2} [0.0000 1.0000]这说明变量 1 发给因子 1 的消息还是均匀的因无先验因子 1 回给变量 1 的消息已偏向x11而因子 1 回给变量 2 的消息已是硬约束[0,1]因f2强制x21。这种细粒度观察是理解消息如何“协商”出全局一致性的唯一途径。4.3 收敛性验证不只是看sum(abs(delta)) tolmain_demo.m用max(abs(belief_new - belief_old))做收敛判断这没问题。但更可靠的验证是检查 belief 是否满足局部一致性对每个变量 v其 belief 应等于∏_{f∈N(v)} msg_{f→v}忽略归一化。在收敛后加一段验证% 验证局部一致性 for v_id 1:length(graph.nodes.var) belief_v beliefs{v_id}; msg_product ones(size(belief_v)); for f_id graph.adj_v2f{v_id} msg_product msg_product .* msg_f2v{f_id}{v_id}; end msg_product msg_product / sum(msg_product); % 归一化 fprintf(Var %d consistency error: %.2e\n, v_id, max(abs(belief_v - msg_product))); end如果输出全是1e-15级别说明算法实现正确如果某变量是1e-2说明那个因子的update_msg_f2v有 bug比如sum_dims算错了。5. 避坑指南5 个真实翻车现场与血泪解决方案这些坑我都在深夜三点的实验室里亲手踩过。它们不来自文档错误而来自对因子图语义的微妙误解或 MATLAB 数值特性的忽视。记下它们能帮你省下至少 8 小时 debug 时间。5.1 现象第一次迭代后所有消息都是NaN原因因子势函数ψ_f包含0值且在update_msg_f2v.m中与msg_prod含0相乘导致weighted出现0后续sum正常但若某msg_v2f初始为[]msg_prod初始化为ones(...)0*10没问题但如果ψ_f有Inf或NaN比如log(0)未处理0*InfNaN就爆发了。解决在build_factor_graph.m中加载ψ_f后强制清洗pot_f(isnan(pot_f) | isinf(pot_f)) eps; % 或设为 0视语义而定 pot_f(pot_f 0) eps; % 概率势函数不能为负5.2 现象迭代 100 轮仍不收敛belief 在两个值间震荡原因图中存在短环cycle of length 4而和积算法在有环图上不保证收敛。本项目示例是树状图无环但如果你自己建图时写了adj_v2f {[1,2], [2,3], [3,1]}三角形就必然震荡。解决运行前调用is_tree_graph(graph)可自行实现num_edges num_nodes - 1且连通。若检测到环要么加阻尼msg_new 0.9*msg_new 0.1*msg_old要么换 Junction Tree 算法——但那就超出本项目范围了。5.3 现象msg_f2v{f_id}{v_id}尺寸报错Subscripted assignment dimension mismatch原因scope中变量 domain 长度与ψ_f.table的维度不匹配。例如scope[1,2]domain_11:2长度 2domain_21:3长度 3则ψ_f.table必须是2×3但你传了3×2或6×1。解决在build_factor_graph.m中加入维度校验expected_size 1; for k 1:length(scope) expected_size expected_size * length(graph.nodes.var(scope(k)).domain); end if numel(pot_f) ~ expected_size error(Factor %d table size %d ! expected %d, f_id, numel(pot_f), expected_size); end % 然后用 reshape 强制为多维reshape(pot_f, domain_sizes)5.4 现象update_msg_v2f报错Index exceeds matrix dimensions在msg_f2v{f_k}{v_id}原因msg_f2v是cell数组但f_k或v_id超出graph.nodes.factor或graph.nodes.var的数量。常见于adj_v2f里写了f_id5但graph.nodes.factor只有 4 个。解决在build_factor_graph.m结尾加断言assert(all(cellfun((x) all(x length(graph.nodes.factor)), adj_v2f)), ... Factor ID in adj_v2f exceeds number of defined factors); assert(all(cellfun((x) all(x length(graph.nodes.var)), adj_f2v)), ... Variable ID in adj_f2v exceeds number of defined variables);5.5 现象belief 看似收敛但P(x_v1)是0.999999999而非1.0下游应用报错原因浮点精度累积误差。sumprod多轮乘除后sum(belief)可能是0.999999999999999虽不影响显示但round(belief)会错。解决在输出 belief 前做强制归一并截断belief belief / sum(belief); belief round(belief * 1e12) / 1e12; % 保留 12 位小数或者更鲁棒地belief belief / (sum(belief) eps);已在normalize_messages.m中实现。6. 二次开发实战从 LDPC 译码到高斯混合模型三个改造技巧与验证方法现在你已跑通、看懂、避过坑下一步是让它为你所用。我不会教你“如何写一个新项目”而是给你三个真实场景下的最小改造路径——每个都能在 30 分钟内完成且有明确验证手段。记住所有改造都始于build_factor_graph.m和main_demo.m的两处修改。6.1 场景一接入 LDPC 译码器二元线性码LDPC 的 Tanner 图是因子图的特例变量节点是比特因子节点是校验方程XOR 为 0。改造只需三步定义变量与因子在main_demo.m开头替换图构建% LDPC H 矩阵4×8每行一个校验方程 H [1 1 0 0 1 0 0 0; 0 1 1 0 0 1 0 0; 0 0 1 1 0 0 1 0; 1 0 0 1 0 0 0 1]; [M, N] size(H); % 构建 adj_v2f变量 j 连接所有 H(:,j)1 的因子 i adj_v2f cell(1,N); for j 1:N adj_v2f{j} find(H(:,j)); end graph build_factor_graph(N, M, adj_v2f, []); % N 变量M 因子设置势函数校验因子ψ_f(x) I(⊕_{u∈scope} x_u 0)即 XOR 为 0 时为 1否则为 0。在build_factor_graph.m中对每个因子f_idscope graph.nodes.factor(f_id).scope; % 生成所有 2^K 种取值计算 XOR all_combos dec2bin(0:(2^K-1))-0; % K×2^K 矩阵 xor_vals mod(sum(all_combos, 1), 2); % 1×2^K pot_f double(xor_vals 0); % 1→满足0→不满足 graph.nodes.factor(f_id).table reshape(pot_f, domain_sizes);验证译码正确性不要看最终 belief而要检查校验方程是否满足对每个因子f_id取belief_v最大概率的值mode domain_v(find(belief_vmax(belief_v)))代入H计算H * mode结果应全为 0mod 2。我在某次调试中发现mode是[1,0,1,0,0,0,0,0]但H(1,:)*mode是 1立刻定位到pot_f的 XOR 计算漏了mod。6.2 场景二高斯混合模型GMM参数学习GMM 的 E-step 本质是变量z_i隐变量指示第 i 个数据点属于哪个 component的后验推断。因子图中变量z_i连接两个因子——p(z_i)先验和p(x_i|z_i)似然。改造要点是处理连续变量的离散化近似项目原始离散GMM 改造domain_z[1,2,3]3 个 component[1,2,3]保持离散p(x_iz_ik)查表psi_f(i,k)p(z_i)向量pi向量piEM 中更新在main_demo.m中对每个数据点i构建一个因子f_i其scope[z_i]table arrayfun((k) normpdf(x_i, mu(k), sigma(k)), 1:K)。注意mu和sigma需在每次 EM 迭代中用新 belief 更新所以main_demo.m要包在 EM 外层循环里。验证方法plot(x_i, belief_z_i)应呈现“软分配”且sum(belief_z_i)接近pi_k。6.3 场景三添加 damping阻尼提升收敛鲁棒性当你的图有环或势函数病态时原始和积易震荡。阻尼是简单有效的工程技巧msg_new α * msg_new (1-α) * msg_old。只需修改main_demo.m的更新循环ALPHA 0.7; % 阻尼系数0.5~0.9 for iter 1:MAX_ITER % ... 计算 msg_v2f_new, msg_f2v_new ... % 应用阻尼 for v_id 1:length(msg_v2f) for f_id graph.adj_v2f{v_id} if ~isempty(msg_v2f{v_id}{f_id}) msg_v2f{v_id}{f_id} ALPHA * msg_v2f_new{v_id}{f_id} ... (1-ALPHA) * msg_v2f{v_id}{f_id}; end end end % 同理更新 msg_f2v... end验证阻尼效果画max_abs_deltavsiter曲线。无阻尼时曲线锯齿状震荡加阻尼后应平滑下降。我一般会试ALPHA0.5,0.7,0.9三条线选下降最快且不反弹的。从那以后我每次接手新图结构都强制走一遍build_factor_graph.m的维度校验和init_messages.m的isnan检查——哪怕只是临时 demo。因为消息传递算法的优雅恰恰藏在那些看似琐碎的索引对齐与数值防护里。希望帮到你。本文还有配套的精品资源点击获取
到MRO，掌握多继承与组合设计)
