循环信念传播理论

1. 概述

信念传播（Belief Propagation, BP） 是概率图模型中最重要的精确/近似推断算法之一。标准BP在树状结构图上可以精确计算边缘概率分布。然而，现实世界的图往往包含环（cycles），这时标准的BP算法不再适用。

循环信念传播（Loopy Belief Propagation, LBP） 将BP的思想推广到带环的图上，虽然失去理论保证，但实践中往往能给出相当准确的结果。

核心挑战：

1. 消息在环中循环可能导致不一致
2. 迭代过程可能发散或振荡
3. 缺乏通用的收敛性保证

研究进展：

• 收敛条件理论分析
• 衰减少数派（Fractional BP）
• 线性规划松弛
• 和-积网络（Sum-Product Networks）

¹²

2. 信念传播基础回顾

2.1 因子图表示

因子图（Factor Graph） 是理解BP算法的关键表示：

$$G=(X,F,E)$$

• $X=\{x_1,\dots ,x_n\}$：变量节点
• $F=\{f_1,\dots ,f_m\}$：因子节点
• $E$：连接变量和因子的边

联合分布分解：

$$p(x_1,\dots ,x_n)=\frac{1}{Z}\prod _{f\in F}f(x_{\text{Ne}(f)})$$

其中 $Z$ 是归一化常数，$\text{Ne}(f)$ 是因子 $f$ 连接的变量集合。

2.2 消息传递规则

和-积算法（Sum-Product Algorithm） 的消息传递：

变量到因子消息：

$${\mu }_{x\to f}(x)=\prod _{h\in \text{Ne}(x)\setminus \{f\}}{\mu }_{h\to x}(x)$$

因子到变量消息：

$${\mu }_{f\to x}(x)=\sum _{x_{\text{Ne}(f)\setminus \{x\}}}\left(f(x_{\text{Ne}(f)})\prod _{y\in \text{Ne}(f)\setminus \{x\}}{\mu }_{y\to f}(y)\right)$$

2.3 边缘分布计算

信念（边缘分布近似）：

$$b_x(x)\propto \prod _{f\in \text{Ne}(x)}{\mu }_{f\to x}(x)$$

归一化：

$$b_x(x)=\frac{{\tilde{b}}_x(x)}{{\sum }_{x^{\prime }}{\tilde{b}}_x(x^{\prime })}$$

3. 循环带来的问题

3.1 消息循环示例

考虑一个简单的4节点环：

x1 --- f1 --- x2
 |              |
f4              f2
 |              |
x4 --- f3 --- x3

消息循环：

1. $x_1$ 发送消息到 $f_1$
2. $f_1$ 发送消息到 $x_2$
3. $x_2$ 发送消息到 $f_2$
4. $f_2$ 发送消息到 $x_3$
5. $x_3$ 发送消息到 $f_3$
6. $f_3$ 发送消息到 $x_4$
7. $x_4$ 发送消息到 $f_4$
8. $f_4$ 发送消息到 $x_1$

此时消息回到起点，形成闭环。后续迭代中，同一消息会被叠加计算。

3.2 不一致性问题

双树图（Two-Tree） 示例：

设 $p(x_1,x_2,x_3)=\frac{1}{Z}{\psi }_1(x_1,x_2){\psi }_2(x_2,x_3)$

树结构：$x_1\to x_2\to x_3$

边缘分布：

$$p(x_2)\propto \sum _{x_1}{\psi }_1(x_1,x_2)\cdot \sum _{x_3}{\psi }_2(x_2,x_3)$$

引入环后（如添加 ${\psi }_3(x_1,x_3)$），上述分解不再成立。

3.3 迭代行为分析

LBP的迭代行为可能表现为：

类型	行为	结果
收敛	消息趋于稳定	得到近似边缘分布
振荡	消息在两个值间跳动	需要阻尼技术
发散	消息值无限增长	需要归一化

4. LBP算法详解

4.1 基本算法

def loopy_belief_propagation(factor_graph, max_iter=100, tol=1e-6, damping=0.0):
    """
    循环信念传播算法
    
    Args:
        factor_graph: 因子图结构
        max_iter: 最大迭代次数
        tol: 收敛容忍度
        damping: 阻尼系数 (0 = 无阻尼, 1 = 完全保持旧值)
    
    Returns:
        beliefs: 边缘分布近似
        messages: 消息历史
    """
    messages = initialize_messages(factor_graph)
    beliefs_history = []
    
    for iteration in range(max_iter):
        old_messages = messages.copy()
        
        # 1. 更新所有因子到变量的消息
        for factor in factor_graph.factors:
            for var in factor.neighbors:
                messages[factor, var] = update_factor_to_var_message(
                    factor, var, messages
                )
        
        # 2. 更新所有变量到因子的消息
        for var in factor_graph.variables:
            for factor in var.neighbors:
                messages[var, factor] = update_var_to_factor_message(
                    var, factor, messages
                )
        
        # 3. 阻尼（可选）
        if damping > 0:
            messages = {
                edge: (1 - damping) * messages[edge] + damping * old_messages[edge]
                for edge in messages
            }
        
        # 4. 检查收敛
        beliefs = compute_beliefs(factor_graph, messages)
        beliefs_history.append(beliefs)
        
        if check_convergence(beliefs_history, tol):
            break
    
    return beliefs, messages
 
 
def update_factor_to_var_message(factor, var, messages):
    """
    更新因子到变量消息
    μ_{f→x}(x) = Σ_{x'~x} [ f(x, x') · Π_{y∈Ne(f)\{x}} μ_{y→f}(y) ]
    """
    neighboring_vars = [v for v in factor.neighbors if v != var]
    
    # 初始化消息域
    message = torch.zeros(var.domain_size)
    
    # 计算消息（边缘化）
    for config in iterate_configs(factor, var):
        # 因子值
        factor_val = factor.evaluate(config)
        
        # 乘以所有其他邻居的消息
        product = factor_val
        for other_var in neighboring_vars:
            msg = messages[other_var, factor][config[other_var]]
            product *= msg
        
        # 累加到对应var状态
        message[config[var]] += product
    
    # 归一化
    if message.sum() > 0:
        message = message / message.sum()
    
    return message
 
 
def update_var_to_factor_message(var, factor, messages):
    """
    更新变量到因子消息
    μ_{x→f}(x) = Π_{h∈Ne(x)\{f}} μ_{h→x}(x)
    """
    message = torch.ones(var.domain_size)
    
    for other_factor in var.neighbors:
        if other_factor != factor:
            message *= messages[other_factor, var]
    
    # 归一化
    if message.sum() > 0:
        message = message / message.sum()
    
    return message

4.2 收敛判定准则

消息变化度量：

$${\Delta }_t=\underset{(i,j)}{\max }\Vert {\mu }_{i\to j}^{(t)}-{\mu }_{i\to j}^{(t-1)}{\Vert }_1$$

收敛当 ${\Delta }_t<ϵ$。

信念变化度量：

$${\Delta }_b=\underset{x}{\max }\Vert b^{(t)}(x)-b^{(t-1)}(x){\Vert }_1$$

4.3 阻尼技术

几何阻尼：

$${\mu }_{i\to j}^{(t)}=(1-\alpha ){\tilde{\mu }}_{i\to j}^{(t)}+\alpha {\mu }_{i\to j}^{(t-1)}$$

其中 $\tilde{\mu }$ 是无阻尼更新，$\alpha \in [0,1]$ 是阻尼系数。

自适应阻尼：

$${\alpha }_t={\alpha }_{\min }+({\alpha }_{\max }-{\alpha }_{\min })\cdot \frac{1}{1+\exp (-\gamma (t-t_0))}$$

5. 收敛性理论

5.1 树状结构的精确性

定理：如果因子图是树状结构（无环），则BP在一次消息传递后收敛，得到精确的边缘分布。

证明：树结构保证每个消息只被计算一次，没有循环依赖。

5.2 环长与收敛性

发现：短环（如三角形、方形）通常导致收敛困难，长环较好。

直觉解释：

• 消息沿环传播的”有效距离”
• 长环中消息衰减更强
• 短环中消息快速叠加导致振荡

量化分析：

• 平均环长 $\bar{L}$
• 收敛概率 $P_{\text{conv}}\approx f(\bar{L})$

5.3 Bethe自由能理论

Bethe自由能：

$$F_{\text{Bethe}}(b,\psi )=\sum _i\sum _{x_i}{b}_i(x_i)\log \frac{b_i(x_i)}{{\psi }_i(x_i)}+\sum _{(i,j)}\sum _{x_{ij}}{b}_{ij}(x_i,x_j)\log \frac{b_{ij}(x_i,x_j)}{{\psi }_{ij}(x_i,x_j)}$$

其中 $b_i$ 是单点信念，$b_{ij}$ 是双点信念，${\psi }_{ij}$ 是势函数。

固定点条件：LBP的收敛点对应Bethe自由能的静止点。

问题：Bethe自由能可能有多个局部最小，每个对应不同的”伪固定点”。

5.4 收敛的充分条件

势函数的性质：

定义：势函数 $\psi$ 是**正则的（regular）**如果对于所有边 $(i,j)$：

$$\sum _{x_i}{\psi }_{ij}(x_i,x_j)=1,\forall x_j$$ $$\sum _{x_j}{\psi }_{ij}(x_i,x_j)=1,\forall x_i$$

定理：如果所有势函数都是正则的，则LBP以指数速度收敛到唯一不动点。

实际意义：归一化势函数（概率势）通常满足正则条件。

6. 变分视角

6.1 变分推断框架

LBP可以理解为最小化Bethe自由能的变分方法：

$$\underset{b\in \mathcal{C}}{\min }{F}_{\text{Bethe}}(b)$$

其中约束集 $\mathcal{C}$ 是信仰传播的一致性约束。

6.2 变分分布族

平均场近似：

$$q(x)=\prod _i{q}_i(x_i)$$

最小化 $\text{KL}(q\Vert p)$ 给出平均场方程。

Bethe近似：

$$q(x)\propto \prod _i{q}_i(x_i)\prod _{(i,j)}\frac{q_{ij}(x_i,x_j)}{q_i(x_i)q_j(x_j)}$$

允许保留双节点相关性。

6.2 能量最小化视角

推土机距离变体（Tree-Reweighted BP）：

$$F_{\text{TRW}}=\underset{b\in {\mathcal{C}}_{\rho }}{\min }\sum _{(i,j)}{\rho }_{ij}{D}_{\text{KL}}(b_{ij}\Vert {\psi }_{ij})+\sum _i{D}_{\text{KL}}(b_i\Vert {\psi }_i)$$

其中 ${\rho }_{ij}$ 是边权重，满足流量约束。

7. 衰减少数派

7.1 Fractional BP

核心思想：将每个环的消息传递”分摊”到多条边上，减少短环的影响。

衰减值 ${\alpha }_{ij}$：

• 每条边分配一个衰减值 ${\alpha }_{ij}\in (0,1)$
• 满足 ${\sum }_{\text{穿过边e的路径}}\prod \alpha =1$

7.2 算法变体

TRBP（Tree-Reweighted BP）：

def trbp(factor_graph, edge_weights, max_iter=100):
    """
    Tree-Reweighted BP
    
    Args:
        edge_weights: 每条边的衰减值 ρ_{ij} ∈ (0,1]
    """
    messages = initialize_messages(factor_graph)
    
    for iteration in range(max_iter):
        for edge in factor_graph.edges:
            var, factor = edge.var, edge.factor
            rho = edge_weights[edge]
            
            # TRBP消息更新
            raw_message = compute_factor_to_var_message(factor, var, messages)
            
            # 衰减
            messages[factor, var] = raw_message ** rho
            
        if check_convergence(messages):
            break
    
    return compute_beliefs(factor_graph, messages)

7.3 收敛保证

定理：TRBP对任意图和任意正的边权重 ${\rho }_{ij}$ 都收敛到唯一不动点。

直觉：权重重参数化将问题转化为”平均树”上的精确推断。

8. 实践应用

8.1 图像去噪

马尔可夫随机场（MRF）：

• 像素作为变量
• 相邻像素间定义势函数
• 观测势函数连接观测和像素

def image_denoising_lbp(observed_image, noise_model, lambda_reg=0.1, max_iter=50):
    """
    使用LBP进行图像去噪
    """
    h, w = observed_image.shape
    
    # 构建MRF
    n_vars = h * w
    variables = [Variable(f'x_{i}', domain_size=256) for _ in range(n_vars)]
    
    # 添加一元势（观测）
    for i in range(n_vars):
        unary_factor = UnaryFactor(variables[i], observed_image.flat[i], noise_model)
        add_factor(unary_factor)
    
    # 添加二元势（平滑先验）
    for i, j in neighbor_pairs(h, w):
        pairwise_factor = PairwiseFactor(
            variables[i], variables[j], 
            lambda_reg * smooth_potential
        )
        add_factor(pairwise_factor)
    
    # 运行LBP
    beliefs, _ = loopy_belief_propagation(fg, max_iter=max_iter, damping=0.5)
    
    # 重建图像
    denoised = torch.stack([b.argmax() for b in beliefs.values()])
    return denoised.reshape(h, w)

8.2 LDPC解码

低密度奇偶校验码天然形成稀疏图，适合LBP。

class LDPCDecoder:
    """LDPC码的置信传播解码器"""
    
    def __init__(self, H, max_iter=100):
        self.H = H  # 校验矩阵
        self.n_vars, self.n_checks = H.shape
        self.max_iter = max_iter
        
    def decode(self, received_llr):
        """
        解码接收到的对数似然比
        
        Args:
            received_llr: 每个比特的LLR (log P(b=0)/P(b=1))
        """
        # 初始化变量到校验节点的消息
        msg_var_to_check = torch.zeros(self.n_vars, self.n_checks)
        
        # 初始化先验消息
        prior = torch.sigmoid(received_llr)  # P(b=0)
        
        for iteration in range(self.max_iter):
            # 1. 变量节点到校验节点
            for v in range(self.n_vars):
                connected_checks = self.H[v].nonzero().flatten()
                
                for c in connected_checks:
                    # LLR域的消息
                    other_checks = [ck for ck in connected_checks if ck != c]
                    if other_checks:
                        # 和积算法（LLR形式）
                        prod = msg_check_to_var[other_checks, v].tanh().prod()
                        msg = 2 * torch.arctanh(prod.clamp(-0.999, 0.999))
                    else:
                        msg = 0.0
                    
                    msg_var_to_check[v, c] = prior[v] + msg
            
            # 2. 校验节点到变量节点
            for c in range(self.n_checks):
                connected_vars = self.H[:, c].nonzero().flatten()
                
                for v in connected_vars:
                    other_vars = [var for var in connected_vars if var != v]
                    
                    # 校验节点消息（硬判决）
                    syndrome = 0
                    for other_v in other_vars:
                        # 变量节点消息的符号
                        sign = 1 if msg_var_to_check[other_v, c] > 0 else -1
                        mag = torch.abs(msg_var_to_check[other_v, c])
                        syndrome += sign * torch.min(mag, torch.tensor(100.0))
                    
                    # 消息 = -sgn(syndrome) * min_magnitude
                    msg_check_to_var[c, v] = -torch.sign(syndrome) * torch.min(
                        torch.abs(syndrome), torch.tensor(100.0)
                    )
            
            # 3. 硬判决
            posterior = prior + msg_var_to_check.sum(dim=1)
            decoded_bits = (posterior < 0).int()
            
            # 4. 校验
            syndrome = (self.H @ decoded_bits.float() % 2).int()
            if syndrome.sum() == 0:
                return decoded_bits.numpy()
        
        return decoded_bits.numpy()

8.3 联合树与LBP比较

选择指南：

场景	推荐方法
树状结构	精确BP
稀疏图（树+少量环）	LBP
密集图	联合树(Junction Tree)
超高树宽	采样方法(MCMC)

9. 与其他推断方法的比较

9.1 方法对比表

方法	精确性	复杂度	适用范围	收敛保证
变量消除	精确	指数（树宽）	小规模	N/A
联合树	精确	指数（树宽）	中等规模	有
LBP	近似	多项式	大规模	通常无
均值场	近似	多项式	大规模	有（单调）
粒子方法	近似	可控	任意	统计保证

9.2 何时使用LBP

适合LBP的场景：

• 图稀疏（每个变量只连接少量因子）
• 树宽较大（联合树不可行）
• 需要快速近似推断
• 势函数是低阶的（便于边缘化）

不适合LBP的场景：

• 高阶势函数（边缘化困难）
• 强耦合的密集图
• 需要精确概率值
• 短环主导的图结构

10. 高级主题

10.1 共识传播（Consensus Propagation）

多层消息传递：

1. 节点聚合邻居信息
2. 因子聚合变量信息
3. 共识迭代直到收敛

收敛性更强，适合更复杂的图结构。

10.2 摊销推断

学习消息函数：

$${\mu }_{f\to x}=g_{\theta }(f,\text{context})$$

使用神经网络参数化消息函数，减少迭代次数。

10.3 随机消息传递

异步更新：

• 不等待所有消息到达
• 随机选择要更新的边
• 类似Gibbs采样

适合大规模分布式系统。

11. 代码实现：完整示例

import torch
import numpy as np
from itertools import product
 
class Factor:
    """因子节点"""
    def __init__(self, variables, potential):
        self.variables = variables  # 连接的变量列表
        self.potential = potential   # 势函数 (tensor或函数)
        
    def evaluate(self, assignment):
        """评估给定赋值下的势函数值"""
        if callable(self.potential):
            return self.potential(assignment)
        else:
            # tensor索引方式
            indices = [assignment[var] for var in self.variables]
            return self.potential[tuple(indices)]
 
 
class Variable:
    """变量节点"""
    def __init__(self, name, domain_size):
        self.name = name
        self.domain_size = domain_size
 
 
class FactorGraph:
    """因子图"""
    def __init__(self):
        self.variables = {}
        self.factors = []
        self.neighbors = {}  # 变量-因子邻接表
        
    def add_variable(self, name, domain_size):
        var = Variable(name, domain_size)
        self.variables[name] = var
        self.neighbors[name] = []
        return var
    
    def add_factor(self, variables, potential):
        factor = Factor(variables, potential)
        self.factors.append(factor)
        
        for var in variables:
            self.neighbors[var.name].append(factor)
        
        return factor
 
 
def run_lbp(fg, max_iter=100, tol=1e-6, damping=0.5):
    """运行LBP"""
    
    # 初始化消息
    messages = {}
    for var_name, var in fg.variables.items():
        for factor in fg.neighbors[var_name]:
            messages[(factor, var)] = torch.ones(var.domain_size) / var.domain_size
            messages[(var, factor)] = torch.ones(var.domain_size) / var.domain_size
    
    # 迭代
    for iteration in range(max_iter):
        old_messages = {k: v.clone() for k, v in messages.items()}
        
        # 更新所有消息
        for var_name, var in fg.variables.items():
            for factor in fg.neighbors[var_name]:
                # 变量 -> 因子
                msg = torch.ones(var.domain_size)
                for other_factor in fg.neighbors[var_name]:
                    if other_factor != factor:
                        msg *= messages[(other_factor, var)]
                
                if damping > 0:
                    msg = (1 - damping) * msg + damping * messages[(var, factor)]
                messages[(var, factor)] = msg / msg.sum()
        
        for factor in fg.factors:
            for var in factor.variables:
                # 因子 -> 变量
                neighbor_vars = [v for v in factor.variables if v != var]
                
                msg = torch.zeros(var.domain_size)
                # 边缘化
                domain_sizes = [v.domain_size for v in factor.variables]
                
                for indices in product(*[range(d) for d in domain_sizes]):
                    assignment = {v: idx for v, idx in zip(factor.variables, indices)}
                    val = factor.evaluate(assignment)
                    
                    # 乘以其他邻居的消息
                    for other_var in neighbor_vars:
                        val *= messages[(other_var, factor)][indices[factor.variables.index(other_var)]]
                    
                    msg[indices[factor.variables.index(var)]] += val
                
                if damping > 0:
                    msg = (1 - damping) * msg + damping * messages[(factor, var)]
                messages[(factor, var)] = msg / (msg.sum() + 1e-10)
        
        # 检查收敛
        max_diff = 0
        for key in messages:
            diff = (messages[key] - old_messages[key]).abs().max().item()
            max_diff = max(max_diff, diff)
        
        if max_diff < tol:
            print(f"Converged at iteration {iteration}")
            break
    
    # 计算信念
    beliefs = {}
    for var_name, var in fg.variables.items():
        belief = torch.ones(var.domain_size)
        for factor in fg.neighbors[var_name]:
            belief *= messages[(factor, var)]
        beliefs[var_name] = belief / belief.sum()
    
    return beliefs, messages
 
 
# 示例：Ising模型
def example_ising_lbp():
    """Ising模型的LBP示例"""
    
    # 2x2网格Ising模型
    fg = FactorGraph()
    
    # 添加变量
    for i in range(4):
        fg.add_variable(f'x{i}', domain_size=2)
    
    # 势函数
    beta = 0.5
    
    # 一元势（外部场）
    def unary_potential(var_assignment):
        h = 0.1
        return torch.exp(beta * h * (1 if var_assignment['x0'] == 1 else -1))
    
    # 二元势（相互作用）
    for (i, j) in [(0, 1), (1, 2), (2, 3), (3, 0)]:  # 环形连接
        def pairwise_potential(var_assignment, i=i, j=j):
            J = 0.3
            val_i = 1 if var_assignment[f'x{i}'] == 1 else -1
            val_j = 1 if var_assignment[f'x{j}'] == 1 else -1
            return torch.exp(beta * J * val_i * val_j)
        
        vars_pair = [fg.variables[f'x{i}'], fg.variables[f'x{j}']]
        fg.add_factor(vars_pair, pairwise_potential)
    
    # 运行LBP
    beliefs, _ = run_lbp(fg, max_iter=100, damping=0.5)
    
    print("\n边缘分布：")
    for var_name, belief in beliefs.items():
        print(f"P({var_name}=0) = {belief[0]:.4f}, P({var_name}=1) = {belief[1]:.4f}")
    
    return beliefs
 
 
if __name__ == '__main__':
    beliefs = example_ising_lbp()

12. 总结

LBP的核心要点：

1. 近似方法：将树状BP推广到带环图
2. 迭代过程：消息不断在环中传播
3. 收敛问题：短环可能导致振荡或发散
4. 阻尼技术：是提高收敛性的实用手段
5. 变分视角：最小化Bethe自由能

实践建议：

1. 优先尝试小阻尼（0.1-0.5）
2. 监控消息变化的范数
3. 限制最大迭代次数避免无限循环
4. 对于关键应用，考虑TRBP或联合树

参考文献

Footnotes

1. Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems. Morgan Kaufmann. ↩

2. Yedidia, J. S., et al. (2005). Constructing Free Energy Approximations and Generalized Belief Propagation Algorithms. IEEE Trans. Info. Theory. ↩