置信传播（BP）译码原理与Python实现：从因子图到LLR域迭代译码

置信传播这几年在通信、机器学习、图像处理领域出现频率高得吓人。做无线通信的翻LDPC码论文几乎是必见BP做图像分割的也常听说基于马尔可夫随机场的BP求解。可不少人第一次看到BP译码那组变量节点更新和校验节点更新公式时心里多少有点发慌这些tanh和atanh到底是从哪冒出来的为什么一条消息一会儿往左传一会儿往右传先提醒一句这里的BP是Belief Propagation也就是置信传播它和神经网络里那个Backpropagation反向传播虽然都缩写为BP但完全是两码事。这篇笔记想把这个算法讲透先从因子图和最大后验估计出发把变量节点、校验节点、初始化、停止条件一步步推导出来然后给出一份完整的Python实现最后聊聊工程实现里最常见的数值坑和性能权衡。适合刚接触迭代译码的研究生、通信工程师以及想在自研链路里快速搭一个BP译码器验证性能的人。1. 先搞清楚BP译码到底在解决什么问题1.1 通信译码问题的本质任何译码算法的目标都可以写成一句话给定接收向量y估计发送的码字c。最优的意义下有两条路最大后验概率MAP和最大似然ML。在码字等概率发送时MAP退化成ML也就是在所有合法码字里找一个让P(y|c)最大的c。这个定义看起来很美但要真去遍历码字做搜索码长一长就直接爆炸。比如一个长度1000的码哪怕只有一半是信息位候选码字也是2^500这个量级没法枚举。那怎么办通信领域很早就想到一个取巧的思路既然是逐个比特发送的能不能不去找整个码字而是去算每个比特的后验概率P(c_ib|y)最后按每个比特的最大后验概率做硬判决直接算这个边缘概率同样很困难但好消息是通信编码结构本身经常能分解成很多局部约束比如校验方程只涉及少数几个比特。BP正是利用这种“局部约束简单”的结构通过图上消息传递把全局后验概率计算拆成一系列局部计算。1.2 BP属于迭代译码和传统代数译码的思路不一样传统汉明码或BCH码的译码常用伴随式、错位定位这一类代数方法本质上是先定错再进行纠错一次搞定。BP不太一样它属于迭代软译码核心动作是反复交换“软信息”也就是概率或对数似然比LLR。每一轮里变量节点告诉校验节点“我觉得这个比特是什么”校验节点再根据自己连接的多个比特告诉变量节点“按校验约束你更像是啥”。这种信息交换称为消息传递BP就是其中最有代表性的算法。LDPC码、Turbo码、Polar码的译码都能用BP或BP变形实现。比如LDPC码的标准迭代译码就是基于稀疏校验矩阵的BPTurbo码的MAP译码本质上也是在一个更复杂的图结构上做消息传递。所以掌握了BP等于掌握了一大类迭代译码器的公共内核。1.3 为什么BP能做到低复杂度逼近最优BP的复杂度主要取决于图中边的数量以及迭代次数和码字数量没有任何关系。一条消息经过一次迭代只做一次乘法或加法LDPC码的校验矩阵又很稀疏每条边的平均计算量很小。所以一个码长几千的LDPC码BP译码做几十次迭代现代CPU上也就几毫秒的事。另一个更吸引人的性质是在因子图是一棵树的情况下BP计算出来的边缘概率是精确的当图有环的时候BP得到的是一种近似称为loopy BP但工程实践发现它在绝大多数好码上依然能逼近ML性能。LDPC码之所以能逼近香农限一个重要前提就是校验矩阵足够稀疏、环足够长让BP的近似误差被压到很低。2. 公式是从哪来的因子图、和积算法与LLR域推导2.1 把译码问题画成因子图BP的一切都从因子图开始。因子图是一种用于表示多变量函数分解结构的图图上有两类节点圆圈表示变量节点方块表示因子节点。每条边连接一个变量和一个因子表示该因子依赖这个变量。对线性分组码来说假设码字c满足m个线性校验方程每个校验方程是一个因子[ f_j(c) \begin{cases} 1, \bigoplus_{i \in N(j)} c_i 0 \ 0, \text{otherwise} \end{cases} ]其中N(j)是第j个校验方程涉及的比特集合。再把加性白高斯噪声AWGN信道建模成观测因子P(y_i|c_i)整个后验分布就可以写成[ P(\mathbf{c}|\mathbf{y}) \propto \left(\prod_{i1}^{n} P(y_i|c_i)\right)\left(\prod_{j1}^{m} f_j(\mathbf{c}_{N(j)})\right) ]这个式子就是因子图背后的数学来源一个大的全局概率分布被分解成若干个只和局部变量有关的因子相乘。变量节点负责把跟它有关的因子串起来因子节点负责表达局部约束。2.2 消息传递两个节点之间的“对话规则”因子图上要算某个变量的边缘概率利用的是和积算法。基本思想是每个变量节点把从其他邻居那里收到的信息汇总再发给某个因子每个因子节点则把除目标变量以外的所有变量做求和消元再发给目标变量。假设我们想看变量节点i传给因子节点j的消息记为(\mu_{x_i \to f_j}(b))表示“除因子j以外其他所有邻居和信道观测综合起来认为x_ib的可信程度”。更新规则是[ \mu_{x_i \to f_j}(b) P(y_i|x_ib) \prod_{j \in N(i)\setminus{j}} \mu_{f_{j} \to x_i}(b) ]因子节点j传给变量节点i的消息则要对其他所有变量求和[ \mu_{f_j \to x_i}(b) \sum_{\substack{x_{i},i\in N(j)\setminus{i}\ \bigoplus_{i\in N(j)}x_{i}0}} \prod_{i\in N(j)\setminus{i}} \mu_{x_{i} \to f_j}(x_{i}) ]这就是和积算法的两个主角。说人话就是变量节点是“传声筒”它做的只是把别人给它的信念乘起来因子节点是“裁判”它要根据自己手里的约束方程结合其他变量的信念给出对某个变量的修正意见。2.3 变量节点更新LLR直接相加直接用概率做消息传递要频繁做乘法数值上很痛苦。工程上几乎都是用对数似然比LLR表示消息[ L \log \frac{P(x0)}{P(x1)} ]LLR的符号代表倾向的比特绝对值代表可信程度。把所有概率消息换成LLR之后变量节点更新就从“相乘”变成“相加”。因为[ \log \frac{\mu_{x_i\to f_j}(0)}{\mu_{x_i\to f_j}(1)}\log \frac{P(y_i|0)}{P(y_i|1)} \sum_{j\in N(i)\setminus{j}} \log \frac{\mu_{f_{j}\to x_i}(0)}{\mu_{f_{j}\to x_i}(1)} ]所以就有[ L_{v_i \to c_j} L_{\mathrm{ch},i} \sum_{j \in N(i)\setminus{j}} L_{c_{j} \to v_i} ]这就是BP译码里变量节点更新的标准形式。信道LLR在BPSK调制、AWGN信道下是(L_{\mathrm{ch},i}2y_i/\sigma^2)。推导很简单假设比特0映射为1比特1映射为-1则(P(y|x)\propto\exp(-(y-x)^2/(2\sigma^2)))代进LLR定义约掉常数就得到(2y/\sigma^2)。2.4 校验节点更新tanh/atanh是怎么冒出来的校验节点更新比变量节点麻烦不少很多初学者就是卡在这里。拆开看其实逻辑很清晰。假设校验节点j连了k个变量约束是这k个比特异或为0。我们要算的是给定了其他k-1个变量的消息后变量i等于b这件事有多可信。这需要对其他变量做全求和。对二进制变量可以直接利用一个漂亮的小技巧任意独立二进制随机变量Z的奇偶概率满足[ P(\text{异或结果为0}) \frac{1 E[(-1)^Z]}{2} ]而(E[(-1)^X] P(X0) - P(X1))。如果通知变量i固定为b其他每个变量贡献一个(a_{i}q_{i\to j}(0)-q_{i\to j}(1))那么校验约束满足的条件概率就是[ P\left(\text{校验通过} \mid x_ib\right)\frac{1 (-1)^b \prod_{i\neq i} a_{i}}{2} ]这个式子里的a很自然地可以写成LLR形式。因为如果某一方消息的LLR为L那么[ a q(0)-q(1) \frac{e^L - 1}{e^L 1} \tanh\frac{L}{2} ]把条件概率比值代进LLR定义[ L_{c_j \to v_i}\log \frac{1 \prod_{i\neq i} \tanh(L_{v_{i}\to c_j}/2)} {1 - \prod_{i\neq i} \tanh(L_{v_{i}\to c_j}/2)} ]再用反双曲正切的关系(\operatorname{atanh}(x)\frac12\log\frac{1x}{1-x})就得到经典公式[ L_{c_j \to v_i}2\operatorname{atanh}\left( \prod_{i\in N(j)\setminus{i}} \tanh\frac{L_{v_{i}\to c_j}}{2} \right) ]整个过程没有魔法就是“先条件固定一个变量对其他变量做奇偶求和再把概率写成LLR”。记住tanh/atanh不是从天而降的复杂函数它们只是二进制变量在LLR域做校验约束时必然出现的工具。2.5 初始化、判决与迭代停止有了两个更新公式整个BP译码流程还差三个动作。第一是初始化。第一次迭代时变量节点还没有收到任何校验节点的消息所以发给校验节点的消息就是信道LLR本身即(L_{v_i\to c_j}L_{\mathrm{ch},i})。第二是后验LLR计算。每一轮迭代后把某个变量收到的所有校验节点消息加上信道LLR就得到该变量的后验LLR[ L_{\mathrm{post},i}L_{\mathrm{ch},i} \sum_{j \in N(i)} L_{c_j \to v_i} ]如果(L_{\mathrm{post},i}0)判为比特0否则判为比特1。第三是停止条件。BP属于迭代算法不能无限跑下去。常用做法是硬判决后做伴随式计算如果所有校验方程都满足说明已经找到一个合法码字立刻停止如果迭代到最大次数还没满足就按最后一次硬判决结果输出。实际系统中也经常加入CRC辅助进一步确认是否译码正确。3. 工程实现为什么必须用LLR域3.1 概率域为什么容易翻车真正写代码时我强烈建议直接用LLR域理由非常现实数值稳定性。概率域里校验节点要做多个消息的乘积码长一长成百上千个小于1的概率连乘很快下溢成0。比如(0.9^{1000})就是(1.7\times10^{-46})再乘几轮浮点数直接变0后续计算全部失真。LLR域把大量乘法变成了加法变量节点更新尤其干净。校验节点更新本质还是乘法但乘的是双曲正切取值范围被限制在[-1,1]配合裁剪可以控制数值范围。因此工程上的BP译码器十有八九都是LLR域实现。如果你在教科书上看到概率域公式理解思路就好真干活时别直接抄。3.2 LLR域的计算与数值裁剪LLR域最大的隐患藏在(2\operatorname{atanh}(\cdot))这一步。当所有入边的tanh值都接近±1时乘积也会非常接近±1此时atanh的输出去向正负无穷。一旦溢出后一轮变量节点更新就会得到无穷大LLR硬判决虽不一定崩但后续迭代会很不稳。我的建议是在校验节点更新之前把乘积clip到([-110^{-12}, 1-10^{-12}])这个区间。这等于给消息幅度设了一个上限本质上是限幅处理对性能影响很小却能避免绝大多数数值异常。另外要注意在迭代早期消息幅度通常很小tanh值也很小乘积自然在安全区间容易出问题的是迭代后期消息置信度变高的时候。理解了这一点再去看“LLR域 vs 概率域”的争论就很容易看穿概率域适合推导理论LLR域适合写代码和上硬件。3.3 Min-Sum近似的诱人之处工程上还有一个非常常见的变形叫Min-Sum也就是最小和算法。它把校验节点更新里的tanh/atanh直接扔掉改成对消息的符号和幅度分别处理消息符号等于所有入边消息符号的异或消息幅度等于所有入边消息绝对值的最小值。[ L_{c_j \to v_i} \approx \left(\prod_{i\neq i}\operatorname{sign}(L_{v_{i}\to c_j})\right) \cdot \min_{i\neq i}|L_{v_{i}\to c_j}| ]这个近似的直观解释是校验节点的置信度取决于“最不确定的那个邻居”。你把最弱的那个邻居的消息幅度作为自己的消息幅度方向由所有邻居的符号共同决定。Min-Sum完全去掉了atanh在硬件上特别讨喜性能损失通常在0.1-0.3 dB左右后续用Normalized Min-Sum或Offset Min-Sum加上一个补偿系数就能把损失压回来。后面我会专门讲这个优化。4. 代码实战从零写一个BP译码器4.1 选型为什么用Hamming(7,4)做例子BP译码最经典的应用对象是LDPC码但如果我这里直接甩一个码长几千的LDPC矩阵公式还没吃透的读者容易把注意力放在矩阵构造上反而漏掉BP本身。我用一个极小、极容易验证的码Hamming(7,4)。它有7个比特、4个信息位、3个校验位校验矩阵只有3行7列单比特纠错能力用穷举都能验证。更关键的是BP的迭代流程、数值处理、停止条件在这个小码上全部保留没有任何偷工减料。跑通之后再换LDPC码只需要替换矩阵和向量维度代码主体几乎不用改。我给出一份校验矩阵HH [[1 1 1 0 1 0 0], [1 0 1 1 0 1 0], [0 1 1 1 0 0 1]]这个矩阵最后三列是单位阵形状前四列是校验系数。利用系统码关系可以推出生成矩阵G我们代码里直接根据H算出来避免手算出错。4.2 整体代码结构代码分为几个模块生成矩阵构建、编码、BPSK调制、信道LLR计算、BP主循环、单次译码测试和BER仿真。我把完整代码贴在下面逻辑尽量一行行写清楚不做花哨优化方便对照公式理解。import numpy as np def build_h(): # (7,4) Hamming码校验矩阵 H np.array([ [1, 1, 1, 0, 1, 0, 0], [1, 0, 1, 1, 0, 1, 0], [0, 1, 1, 1, 0, 0, 1] ], dtypeint) return H def build_g(H): # 当 H [P | I] 时系统码生成矩阵为 G [I | P^T] m, n H.shape k n - m P H[:, :k] G np.concatenate([np.eye(k, dtypeint), P.T], axis1) return G def encode(u, G): # 信息序列编码为码字 return (u G) % 2 def bpsk_mod(c): # 比特0映射为 1比特1映射为 -1 return 1.0 - 2.0 * c def ch_llr(y, sigma): # AWGN BPSK 的信道LLR return 2.0 * y / (sigma * sigma) def bp_decode(y, H, sigma, max_iter100): m, n H.shape Lch ch_llr(y, sigma) # 邻居索引表 neighbors_of_var [np.where(H[:, i] 1)[0] for i in range(n)] neighbors_of_check [np.where(H[j, :] 1)[0] for j in range(m)] # Lq[j,i]变量节点 i 发给校验节点 j 的消息 # Lr[j,i]校验节点 j 发给变量节点 i 的消息 Lq np.zeros((m, n)) Lr np.zeros((m, n)) # 初始化变量节点消息取信道LLR for i in range(n): for j in neighbors_of_var[i]: Lq[j, i] Lch[i] for it in range(max_iter): # 1. 校验节点更新 for j in range(m): nbrs neighbors_of_check[j] for i in nbrs: prod 1.0 for ii in nbrs: if ii ! i: prod * np.tanh(0.5 * Lq[j, ii]) prod np.clip(prod, -1.0 1e-12, 1.0 - 1e-12) Lr[j, i] 2.0 * np.arctanh(prod) # 2. 变量节点更新 for i in range(n): nbrs neighbors_of_var[i] for j in nbrs: s Lch[i] for jj in nbrs: if jj ! j: s Lr[jj, i] Lq[j, i] s # 3. 后验LLR与硬判决 Lpost np.zeros(n) for i in range(n): s Lch[i] for j in neighbors_of_var[i]: s Lr[j, i] Lpost[i] s xhat (Lpost 0).astype(int) # LLR 0 判为1否则判为0 # 4. 停止条件伴随式全零 if np.all((H xhat) % 2 0): return xhat, it 1 return xhat, max_iter代码里有个细节值得讲一下消息矩阵为什么是(m, n)两维而不是用字典或链表因为H是稀疏的但Hamming(7,4)很小直接用稠密矩阵外加“0值表示无边”最简单直观。换到真正的LDPC码时建议用邻接列表加上两个方向的稀疏数组消息存储会省很多内存。4.3 单次译码测试接下来做一个单次测试随机生成一个消息编码后调制人为把第3个符号翻转再送入BP译码器。if __name__ __main__: H build_h() G build_g(H) u np.array([1, 0, 1, 0]) c encode(u, G) x bpsk_mod(c) # 制造一个接收错误第3个符号翻转 y x.copy() y[2] -y[2] sigma 0.5 xhat, iters bp_decode(y, H, sigma, max_iter50) print(原始消息:, u) print(译码结果:, xhat[:4]) print(迭代次数:, iters)运行正常的话译码结果的前4位应该和原始消息完全一致迭代次数一般在几次到十几次之间。你可以把这个测试当成冒烟测试只要这个用例能过说明BP主循环基本正确如果这个用例都过不了赶紧回头检查消息初始化或变量节点更新公式里的邻居排除逻辑。4.4 性能仿真算一条BER曲线单次测试通过后就可以跑误码率仿真了。下面是完整的BER仿真代码支持输入若干Eb/N0点每个点跑若干个帧统计系统比特误码率。def run_ber_simulation(ebn0_db_list, frame_num2000, max_iter50): H build_h() G build_g(H) k, n G.shape for ebn0_db in ebn0_db_list: snr_lin 10 ** (ebn0_db / 10.0) R k / n sigma np.sqrt(1.0 / (2.0 * R * snr_lin)) bit_err 0 total_bit 0 for _ in range(frame_num): u np.random.randint(0, 2, sizek) c encode(u, G) x bpsk_mod(c) noise np.random.normal(0.0, sigma, sizen) y x noise xhat, _ bp_decode(y, H, sigma, max_iter) bit_err np.sum(xhat[:k] ! u) total_bit k ber bit_err / total_bit print(fEb/N0 {ebn0_db:5.2f} dB, BER {ber:.6f}) if __name__ __main__: run_ber_simulation([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0], frame_num5000)注意sigma的计算公式BPSK信号能量归一化为1码率Rk/n由(E_b/N_0(E_s/N_0)/R)和(E_s/N_01/(2\sigma^2))可以推出(\sigma^21/(2R\cdot E_b/N_0))。这个是仿真里最容易错的点sigma差了根号2信噪比曲线就会整体平移一截。4.5 换用MATLAB需要注意的细节现在很多做通信仿真的同学还是在用MATLAB。Python这段代码移植到MATLAB不难但有三个点要注意。第一MATLAB索引从1开始邻居表生成时要把索引整体加1。第二MATLAB的atanh是atanh函数裁剪乘积时同样用max(min(prod, 1-eps), -1eps)。第三MATLAB的矩阵按列存储循环顺序影响不大但建议先用小码把结果和Python版本逐bit对比确认无误再换大矩阵。我自己经常先写Python调通逻辑再用MATLAB做硬件链路联合仿真两边对照能省很多调试时间。5. 实操中容易踩的坑与排查方法5.1 atanh输出无穷大这是BP译码器最先露出的坑。当你看到atanh报inf或者overflow的warning时不要怀疑编译环境问题基本都出在校验节点更新的乘积上。原因是进入迭代后期消息LLR的绝对值很大tanh值几乎等于±1多个±1相乘浮点数有时候会撞到精确的±1而atanh(±1)正好是无穷大。解决方式就是我前面代码里写的乘积在做atanh之前clip到±(1-1e-12)以内。注意这里的1e-12不是随便选的它对应LLR幅度的上限大概在十几个量级足够满足绝大多数信噪比范围。如果你用定点或更低精度计算可能要用更保守的裁剪值。5.2 迭代不收敛消息来回震荡另一种常见现象是迭代很多次也不满足停止条件硬判决结果在几个码字之间反复横跳。短环是首要嫌疑。当因子图里存在4环时两条消息路径会很快形成正反馈循环消息被不断放大最后震荡。进一步说BP在无环图上是精确的在带环图上只是近似而环越短近似误差越大。这也是LDPC码设计中必须保证“girth至少6”的原因汉明码这种小码天然存在短环所以它的BP性能并不比代数译码强多少但作为教学例子非常合适。如果你的工程码是LDPC建议先用已知无短环的矩阵跑通不要随手生成一个随机稀疏矩阵就当LDPC用。5.3 噪声方差设置不对导致性能崩信道LLR计算里只要σ设错所有消息幅度就全错收敛速度和误码率都会明显恶化。最常见的错误是把σ²当成σ或者上下行链路的信噪比换算搞混。记住公式(L_{ch}2y/\sigma^2)这里的σ是复噪声的实部标准差即噪声功率谱密度(N_0/2)再开方。如果你的系统模型是复数基带还要区分每维噪声功率和总噪声功率别再把系数漏掉一个2。5.4 常见问题速查表现象可能原因处理方式atanh溢出或产生inf乘积为±1对乘积clip到±(1-1e-12)迭代一直不停止硬判决震荡因子图存在短环常见于4环换环长更大的校验矩阵降低迭代上限BER曲线整体偏移信号映射或sigma设置错误检查0/1到±1的映射重推σ与Eb/N0关系单次正常BER却很高初始化错误或变量更新时没排除回传检查邻居排除逻辑确认消息方向高速率下性能差置信度过高改用NMS/OMS或对消息做限幅6. 工程上还能怎么优化6.1 从BP到Min-Sum再到NMS/OMS如果只是做算法验证标准BP已经够用。但真要写进芯片没人愿意为每个校验节点计算还好几个tanh和atanh。Min-Sum用“符号乘积最小值”近似了校验节点更新直接省掉了非线性函数复杂度直线下降。代价是估值偏乐观因为真实校验节点消息幅度比最小值还要小一些。为了补偿这种偏差工程上常见的两种修正Normalized Min-SumNMS在Min-Sum结果上乘一个小于1的系数典型值0.75左右Offset Min-SumOMS则改为减去一个偏移量典型值0.5左右。两个修正系数都需要根据码型和信噪比微调和经验值最接近的结果是最优的。我是建议初学者先用NMS参数少效果稳定调起来最省心。6.2 分层调度与提前终止标准BP的调度方式叫flooding也就是每轮迭代里所有校验节点和变量节点都同时更新。这种方式实现简单但收敛速度一般。分层调度layered scheduling的思想是按顺序逐行更新校验节点每处理完一个校验节点立刻用新消息更新相关变量节点再传给下一个校验节点用。实测下来分层调度大概用一半的迭代次数就能达到flooding的误码率水平而且存储更省很多LDPC码芯片都这么做。提前终止也是工程标配。代码里已经写了伴随式检查只要找到合法码字立即退出。要注意的是如果信道极差BP可能一直找不到合法码字这时需要用一个最大迭代次数兜底避免死循环。6.3 硬件实现视角硬件做BP译码器时消息用定点表示LLR通常量化成6到8比特。量化太狠性能损失大量化太宽面积和功耗又上去了。好在Min-Sum完全没有非线性函数非常适合定点。对于标准BPtanh/atanh会做成长表查询或者在DSP里用多项式逼近代价都不小。所以我见到的商用LDPC译码器绝大多数用的都是Min-Sum或NMS很少直接在硬件里铺标准BP。还有一个容易被忽略的工程问题消息存储。一位带符号定点LLR最少也要6比特一个中等规模的LDPC码有上万条边单是消息矩阵就要好几万比特的存储。所以在硬件上通常不存完整的(m, n)矩阵而是用压缩后的邻居表存储消息配合分层调度把存储需求再砍一半。收尾一点个人经验我自己第一次实现BP时卡得最久的地方不是数学公式而是那个小得不能再小的细节变量节点更新时必须排除回传消息的目标节点。一个变量节点有N个邻居更新对第j个邻居的输出时只能合计其他N-1个邻居的消息绝不能把自己要发往的对象也算进去。把公式和代码逐行对照把维度写清楚这类方向性错误其实非常容易避免。如果让我给初学者一句建议我会说先用Hamming(7,4)这种小码把BP跑通并和暴力ML译码对比再上LDPC。不要一上来就盯着几千块钱的仿真大石头小石头上练好挥锤姿势大石头才能真正敲得动。等你把标准BP的循环结构刻进脑子之后再去看NMS、分层调度、量化定点每一步都有清晰的落点而不会觉得又是一堆黑科技。