因子图与消息传递算法

1 引言

因子图（Factor Graph）是概率图模型中一种优雅而强大的表示形式，它将联合分布分解为多个局部函数的乘积，并通过消息传递算法（Message Passing）进行高效推断。这种表示方法最早由 Kschischang、Frey 和 Loeliger 在 2001 年的开创性论文中系统阐述¹，现已成为现代概率推断的核心框架。

1.1 为什么要用因子图？

在概率图模型中，我们经常需要计算边缘分布、边缘概率或进行最大后验推断。对于 $n$ 个变量的联合分布，直接计算的时间复杂度是指数级的 $O (d^{n})$ （ $d$ 为每个变量的取值数）。因子图通过以下两个关键思想将复杂度降低到多项式级别：

稀疏分解：利用条件独立性，将联合分布分解为多个局部函数的乘积
消息传递：利用动态规划/递归分解，将全局计算分解为局部消息的迭代

1.2 消息传递的统一视角

消息传递是一种通用计算范式，其思想渗透于现代深度学习的多个领域：

领域	消息传递形式
GNN 中的消息传递	邻居节点特征聚合
Transformer 中的注意力	Query-Key-Value 交互
变分自编码器	近似后验的消息传递
贝叶斯神经网络	参数推断的消息传递

理解因子图的消息传递机制，是理解这些现代方法本质的基础。

2 因子图基础

2.1 形式化定义

定义 2.1（因子图）：因子图是一个二部图 $G = (X, F, E)$ ，由以下三部分组成：

变量节点集合 $X = {x_{1}, x_{2}, \dots, x_{n}}$ ，每个 $x_{i}$ 取值于离散集合 $X_{i}$ 或连续空间
因子节点集合 $F = {f_{1}, f_{2}, \dots, f_{m}}$ ，每个因子是一个非负函数 $f_{a} : \prod_{i \in N (a)} X_{i} \to R_{\geq 0}$
边集合 $E \subseteq X \times F$ ， $(x_{i}, f_{a}) \in E$ 当且仅当 $x_{i}$ 是因子 $f_{a}$ 的参数

记号约定： $N (a) = {i : (x_{i}, f_{a}) \in E}$ 表示与因子 $f_{a}$ 相连的变量索引集合， $N (i) = {a : (x_{i}, f_{a}) \in E}$ 表示与变量 $x_{i}$ 相连的因子索引集合。

定义 2.2（全局函数分解）：给定因子图 $G$ ，全局函数（通常是未归一化的联合分布）可以分解为：

g (X) = a = 1 \prod m f_{a} (X_{a})

其中 $X_{a} = {x_{i} : i \in N (a)}$ 是因子 $f_{a}$ 的参数变量集合。

归一化分布：如果 $g (X)$ 是一个概率分布，则需要归一化：

p (X) = \frac{1}{Z} a = 1 \prod m f_{a} (X_{a}), Z = X^{'} \sum a = 1 \prod m f_{a} (X_{a}^{'})

其中 $Z$ 称为配分函数（Partition Function）。

2.2 树结构因子图

树结构因子图是一类重要的特殊情形，其变量图（忽略因子节点）是树状的。在树结构上，消息传递算法可以精确计算所有边缘分布。

定义 2.3（树结构判定）：一个因子图是树结构的，当且仅当：

忽略因子节点后，变量节点形成的图是连通的且无环
每个变量节点最多连接两个因子节点（可选，此条件针对链式结构）

引理 2.1：树结构因子图上的消息传递可以精确计算所有边缘分布，且算法终止后不产生歧义。

2.3 从因子图到因子节点图的映射

因子图 $G = (X, F, E)$ 可以转换为因子节点图 $G^{'} = (F, E^{'})$ ，其中边 $(f_{a}, f_{b}) \in E^{'}$ 当且仅当存在变量 $x_{i}$ 同时连接 $f_{a}$ 和 $f_{b}$ 。这种转换对于理解消息传递中的计算顺序非常有用。

# 伪代码：因子图到因子节点图的转换
def factor_graph_to_factor_graph(G):
    """
    将因子图转换为因子节点图
    G: (variables, factors, edges)
    返回: 因子节点之间的连接关系
    """
    variable_to_factors = defaultdict(list)
    for (var, factor) in G.edges:
        variable_to_factors[var].append(factor)
    
    factor_edges = set()
    for var, factors in variable_to_factors.items():
        # 对于同时连接同一个变量的因子对，建立边
        for i in range(len(factors)):
            for j in range(i + 1, len(factors)):
                factor_edges.add((factors[i], factors[j]))
    
    return factor_edges

2.4 实例：图像去噪模型

考虑一个经典的图像去噪问题，这是因子图在实际应用中的典型案例。

问题设定：设观测图像为 $Y = {y_{1}, y_{2}, \dots, y_{n}}$ ，真实图像为 $X = {x_{1}, x_{2}, \dots, x_{n}}$ ，其中每个像素 $x_{i} \in {0, 1}$ 为二值变量。

马尔可夫随机场（MRF）分解：根据 Hammersley-Clifford 定理，正规化 MRF 可以分解为：

p (X ∣ Y) \propto 一元势函数 i \prod ϕ_{i} (x_{i}) \cdot 二元势函数 (i, j) \in E \prod ψ_{ij} (x_{i}, x_{j})

势函数设计：

数据项（一元势函数）： $ϕ_{i} (x_{i}) = p (y_{i} ∣ x_{i})$ ，衡量像素值与观测的一致性
平滑项（二元势函数）： $ψ_{ij} (x_{i}, x_{j}) = exp (- β \cdot 1 [x_{i} \neq = x_{j}])$ ，鼓励相邻像素取相同值

因子图表示：

     y_1        y_2        y_3
      |         |         |
   [φ_1]      [φ_2]      [φ_3]
      \         |         /
       \        |        /
      (x_1)---ψ_12---(x_2)---ψ_23---(x_3)

每个二元势函数 $ψ_{ij}$ 对应一个因子节点，连接变量 $x_{i}$ 和 $x_{j}$ 。

3 和-积算法（Sum-Product Algorithm）

3.1 问题表述

给定因子图 $G = (X, F, E)$ ，和-积算法的目标是高效计算边缘分布：

p (x_{i}) = X ∖ {x_{i}} \sum p (X) = \frac{1}{Z} X ∖ {x_{i}} \sum a = 1 \prod m f_{a} (X_{a})

以及最大后验推断：

\hat{X} = ar g X max a = 1 \prod m f_{a} (X_{a})

3.2 消息传递规则的推导

3.2.1 从变量到因子的消息

设变量 $x_{i}$ 连接到因子 $f_{a}$ ，我们需要计算消息 $μ_{x_{i} \to f_{a}} (x_{i})$ ，表示 $x_{i}$ 从其他因子获得的信息汇总。

推导：变量 $x_{i}$ 收到的所有消息来自其邻居因子 $N (i) = {b : (x_{i}, f_{b}) \in E}$ 。当 $x_{i}$ 向 $f_{a}$ 发送消息时，需要聚合来自除 $f_{a}$ 以外所有邻居的消息：

μ_{x_{i} \to f_{a}} (x_{i}) = b \in N (i) ∖ {a} \prod μ_{f_{b} \to x_{i}} (x_{i})

直觉理解：这是一个乘积-聚合（product-aggregation）操作。每个邻居因子提供一条关于 $x_{i}$ 可能取值的信息，各条信息通过乘积（相当于概率的乘法/证据的累积）进行组合。

3.2.2 从因子到变量的消息

设因子 $f_{a}$ 连接到变量 $x_{i}$ ，我们需要计算消息 $μ_{f_{a} \to x_{i}} (x_{i})$ 。

推导：因子 $f_{a}$ 的局部函数为 $f_{a} (X_{a})$ ，其中 $X_{a} = {x_{j} : j \in N (a)}$ 。为了得到关于 $x_{i}$ 的消息，我们需要对其他所有变量进行边缘化：

μ_{f_{a} \to x_{i}} (x_{i}) = X_{a} ∖ {x_{i}} \sum f_{a} (X_{a}) \cdot j \in N (a) ∖ {i} \prod μ_{x_{j} \to f_{a}} (x_{j})

关键观察：

输入消息 $μ_{x_{j} \to f_{a}}$ 已经聚合了除 $f_{a}$ 以外所有其他因子对 $x_{j}$ 的信息
乘积 $f_{a} (X_{a}) \cdot \prod_{j \neq = i} μ_{x_{j} \to f_{a}} (x_{j})$ 可以看作一个未归一化的联合分布
对该联合分布边缘化 $x_{j} (j \neq = i)$ 得到关于 $x_{i}$ 的消息

3.2.3 消息计算的终止条件

当所有消息传递完成后，每个变量 $x_{i}$ 的边缘分布（未归一化）为：

\tilde{p} (x_{i}) = a \in N (i) \prod μ_{f_{a} \to x_{i}} (x_{i})

归一化后：

p (x_{i}) = \frac{p ~ ( x _{i} )}{\sum _{x_{i}^{'}} p ~ ( x _{i}^{'} )}

3.3 算法伪代码

def sum_product(factor_graph):
    """
    和-积算法主循环
    factor_graph: (variables, factors, edges)
    返回: 所有变量的边缘分布
    """
    # 初始化：从因子节点到变量节点的消息初始化为均匀分布
    messages = initialize_messages(factor_graph)
    
    # 迭代消息传递直到收敛（或达到最大迭代次数）
    for iteration in range(max_iterations):
        # 1. 变量 -> 因子 消息更新
        for variable in factor_graph.variables:
            for neighbor_factor in variable.neighbors:
                messages[variable -> neighbor_factor] = \
                    product_of_other_messages(variable, neighbor_factor)
        
        # 2. 因子 -> 变量 消息更新
        for factor in factor_graph.factors:
            for neighbor_variable in factor.neighbors:
                messages[factor -> neighbor_variable] = \
                    sum_over_other_variables(factor, neighbor_variable)
        
        # 3. 检查收敛（可选）
        if has_converged(messages):
            break
    
    # 计算边缘分布
    marginals = {}
    for variable in factor_graph.variables:
        marginals[variable] = product_of_incoming_messages(variable)
    
    return marginals

3.4 树结构上的消息传递顺序

对于树结构因子图，消息传递存在拓扑序（Topological Order）：

引理 3.1：在树结构因子图中，从叶子节点开始的消息传递可以在有限步骤内完成。

算法：

找到所有叶子变量节点和叶子因子节点
从叶子节点向内部节点发送消息
重复直到所有消息传递完成

复杂度分析：对于有 $n$ 个变量和 $m$ 个因子的树结构因子图，每个消息的计算复杂度为 $O (d^{k})$ ，其中 $k$ 是因子的最大度数， $d$ 是变量的取值数。总复杂度为 $O (n \cdot d^{m a x_{a} ∣ X_{a} ∣})$ ，相对于指数级的暴力计算有显著降低。

3.5 连续变量的扩展

当变量取值于连续空间时，求和替换为积分：

μ_{f_{a} \to x_{i}} (x_{i}) = \int_{X_{a} ∖ {x_{i}}} f_{a} (X_{a}) j \in N (a) ∖ {i} \prod μ_{x_{j} \to f_{a}} (x_{j}) d X_{a}^{∖ i}

高斯因子的特殊性质：如果所有因子都是高斯分布，则消息传递可以解析求解，边缘分布仍是高斯分布。这在因子图的高斯消息传递（Gaussian Message Passing）中非常重要，可用于线性动态系统和卡尔曼滤波的推导。

3.6 环状因子图与近似方法

当因子图包含环（cycle）时，精确消息传递不再可行，因为消息传递规则可能产生不一致的结果或无限循环。

loopy belief propagation (LBP)：将消息传递规则直接应用于带环图，虽然不保证收敛，但通常能获得不错的近似结果。

收敛条件：对于二元势函数满足亚模性（submodularity）条件的问题，LBP 收敛到全局最优解。

4 Max-积算法（Max-Product Algorithm）

4.1 问题与和-积算法的联系

Max-积算法与和-积算法几乎完全相同，唯一的区别是将求和操作替换为取最大值操作：

μ_{f_{a} \to x_{i}}^{m a x} (x_{i}) = X_{a} ∖ {x_{i}} max f_{a} (X_{a}) \cdot j \in N (a) ∖ {i} \prod μ_{x_{j} \to f_{a}}^{m a x} (x_{j})

4.2 最大后验推断

Max-积算法的输出是最大后验估计：

\hat{X} = ar g X max p (X) = ar g X max a \prod f_{a} (X_{a})

在树结构图上，Max-积算法可以同时给出 MAP 估计和各变量的边缘最大概率。

4.3 追踪最优解：反向指针

为了从消息中恢复最优解，需要记录产生最大值的中间变量：

def max_product_with_backpointers(factor_graph):
    """
    带反向指针的Max-积算法
    """
    messages = {}
    backpointers = {}  # 记录反向指针
    
    # 前向消息传递
    for variable in topological_order(factor_graph):
        for neighbor_factor in variable.neighbors:
            messages[variable -> neighbor_factor] = \
                max_product_of_other_messages(variable, neighbor_factor)
    
    for factor in topological_order(factor_graph):
        for neighbor_variable in factor.neighbors:
            msg, bp = max_sum_over_variables(factor, neighbor_variable)
            messages[factor -> neighbor_variable] = msg
            backpointers[factor -> neighbor_variable] = bp
    
    # 从根节点开始反向追踪
    solution = {}
    for root in find_root_variables(factor_graph):
        solution[root] = argmax(messages[root])
        trace_back(solution, root, backpointers)
    
    return solution

4.4 数值稳定性

直接使用 max 操作可能导致数值下溢。对于概率分布，通常使用对数空间变换：

lo g μ_{f_{a} \to x_{i}}^{m a x} (x_{i}) = X_{a} ∖ {x_{i}} max lo g f_{a} (X_{a}) + j \in N (a) ∖ {i} \sum lo g μ_{x_{j} \to f_{a}}^{m a x} (x_{j})

乘积变为求和，避免了数值问题。

5 变分消息传递（Variational Message Passing）

5.1 变分推断框架

当精确推断困难或不可行时，变分消息传递提供了一种近似推断的框架。

目标：近似后验分布 $p (X ∣ Y)$ ，寻找一个易于处理的变分分布 $q (X)$ 使其最小化 KL 散度：

q^{*} (X) = ar g q min KL (q (X) ∥ p (X ∣ Y))

等价于最大化证据下界（ELBO）：

L (q) = E_{q} [lo g p (X, Y)] - E_{q} [lo g q (X)]

5.2 平均场近似

平均场假设：变分分布可以被完全分解：

q (X) = i = 1 \prod n q_{i} (x_{i})

这假设所有变量在给定其他变量时是条件独立的，是最常用的变分近似形式。

另一种分解形式：基于因子图的分解：

q (X) = a = 1 \prod m q_{a} (X_{a})

其中 $q_{a}$ 是与因子 $f_{a}$ 相关的势函数。

5.3 坐标上升变分推断（CAVI）

坐标上升变分推断（Coordinate Ascent Variational Inference）通过交替优化每个变分因子来最大化 ELBO。

更新规则推导：考虑优化 $q_{i} (x_{i})$ ，固定其他变分因子 $q_{j \neq = i} (x_{j})$ ：

lo g q_{i}^{*} (x_{i}) = E_{q_{- i}} [lo g p (X, Y)] + const

其中 $E_{q_{- i}} [\cdot]$ 表示对 $X ∖ {x_{i}}$ 的期望。

因子分解下的更新：如果变分分布按因子分解，则：

lo g q_{a}^{*} (X_{a}) \propto E_{q_{- a}} [lo g f_{a} (X_{a})] + i \in N (a) \sum E_{q_{i}} [lo g q_{i} (x_{i})] + const

5.4 期望传播（Expectation Propagation）

期望传播（Expectation Propagation, EP）是另一种变分消息传递方法，它在局部近似和全局一致性之间取得平衡。

核心思想：将每个因子的精确影响近似为一个指数族分布，然后迭代更新直到收敛。

EP 更新：

删除：从全局近似中移除因子 $f_{a}$ 的影响
求积：将移除后的分布与精确因子 $f_{a}$ 结合
近似：将结果投影回指数族分布

def expectation_propagation(factor_graph, num_iterations=100):
    """
    期望传播算法
    """
    # 初始化：每个因子的变分近似为均匀分布
    approximations = {f: ExponentialFamily() for f in factor_graph.factors}
    
    for iteration in range(num_iterations):
        for factor in factor_graph.factors:
            # 1. 删除步骤：计算去掉该因子后的全局近似
            cavity = compute_cavity_distributions(approximations, factor)
            
            # 2. 乘积步骤：将精确因子与空腔分布结合
            tilted = cavity * factor
            
            # 3. 投影步骤：将倾斜分布投影回指数族
            approximations[factor] = project_to_exponential_family(tilted)
        
        if has_converged(approximations):
            break
    
    return approximations

5.5 变分消息传递与和-积算法的联系

定理 5.1：当变分分布族包含精确后验分布时，平均场变分推断的消息传递规则与和-积算法在不动点处等价。

直觉理解：

变分消息传递寻找 ELBO 的驻点
和-积算法寻找信念传播的不动点
在适当的假设下，这两个不动点是相同的

这揭示了变分推断和置信传播之间的深刻联系。

6 与神经网络的联系

6.1 消息传递作为神经网络层

消息传递规则可以重新解释为可学习的神经网络操作，这为图神经网络（Graph Neural Networks, GNN）提供了理论基础。

消息传递范式（来自 GNN 框架）：

消息计算（Message Computation）：
$m_{ij}^{(l)} = MSG^{(l)} (h_{i}^{(l)}, h_{j}^{(l)}, e_{ij})$
消息聚合（Message Aggregation）：
$a_{i}^{(l)} = AGG^{(l)} ({m_{ij}^{(l)} : j \in N (i)})$
节点更新（Node Update）：
$h_{i}^{(l + 1)} = UPDATE^{(l)} (h_{i}^{(l)}, a_{i}^{(l)})$

与因子图消息传递的对应关系：

因子图	图神经网络
因子节点 $f_{a}$	消息函数 $MSG$
变量节点 $x_{i}$	图节点 $v_{i}$
消息 $μ_{x_{j} \to f_{a}}$	邻居表示 $h_{j}$
乘积-聚合	求和/注意力聚合

6.2 门控消息传递

门控循环单元（Gated Recurrent Unit, GRU）和门控图神经网络使用门控机制控制消息流动：

h_{i}^{(l + 1)} = GRU h_{i}^{(l)}, j \in N (i) \sum m_{ij}^{(l)}

这可以理解为软版本的消息过滤，类似于因子图中根据消息强度自适应调整贡献。

6.3 注意力机制

Transformer 中的自注意力机制是一种软消息传递：

Attention (Q, K, V) = softmax (\frac{Q K ^{T}}{d _{k}}) V

与因子图的联系：

Query $Q$ 相当于因子节点”请求”信息
Key $K$ 相当于变量节点”提供”信息的能力
注意力权重 $softmax (Q K^{T} / d_{k})$ 相当于消息传递中的归一化权重

class FactorGraphAttention(nn.Module):
    """因子图风格的注意力机制"""
    
    def __init__(self, dim, num_heads=8):
        super().__init__()
        self.num_heads = num_heads
        self.head_dim = dim // num_heads
        
        # 消息函数：Query, Key, Value 投影
        self.W_q = nn.Linear(dim, dim)
        self.W_k = nn.Linear(dim, dim)
        self.W_v = nn.Linear(dim, dim)
        
        # 聚合后的更新函数
        self.update_net = nn.Sequential(
            nn.Linear(dim, dim),
            nn.ReLU(),
            nn.Linear(dim, dim)
        )
    
    def forward(self, node_features, edge_index=None):
        """
        node_features: [num_nodes, dim]
        edge_index: [2, num_edges]，可选的稀疏连接
        """
        batch_size, num_nodes, dim = node_features.shape
        
        # 1. 计算 Query, Key, Value
        Q = self.W_q(node_features)  # [B, N, D]
        K = self.W_k(node_features)
        V = self.W_v(node_features)
        
        # 2. 计算注意力权重（消息强度）
        # QK^T / sqrt(d_k) 相当于计算因子节点到变量节点的消息
        attention_scores = torch.matmul(Q, K.transpose(-2, -1))
        attention_scores = attention_scores / (self.head_dim ** 0.5)
        
        # 3. 归一化得到注意力权重
        attention_weights = F.softmax(attention_scores, dim=-1)
        
        # 4. 消息传递：聚合邻居信息
        # 相当于 sum_j (message_j * attention_weight_j)
        messages = torch.matmul(attention_weights, V)
        
        # 5. 更新节点表示
        updated_features = self.update_net(messages)
        
        return updated_features, attention_weights

6.4 可微分消息传递

可微分消息传递是概率编程和神经概率推断的核心思想。通过将消息函数参数化，消息传递规则可以通过反向传播进行端到端学习。

class DifferentiableMessagePassing(nn.Module):
    """
    可微分消息传递模块
    结合因子图结构和神经网络的可学习性
    """
    
    def __init__(self, node_dim, edge_dim, hidden_dim):
        super().__init__()
        
        # 消息函数：从因子到变量的消息网络
        self.factor_to_node_msg = nn.Sequential(
            nn.Linear(node_dim + edge_dim, hidden_dim),
            nn.LayerNorm(hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, node_dim)
        )
        
        # 消息函数：从变量到因子的消息网络
        self.node_to_factor_msg = nn.Sequential(
            nn.Linear(node_dim, hidden_dim),
            nn.LayerNorm(hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, node_dim)
        )
        
        # 聚合函数：可以是求和、平均或注意力
        self.aggregation = AttentionAggregator(hidden_dim)
        
        # 更新函数：GRU风格的更新
        self.update_fn = nn.GRUCell(hidden_dim, node_dim)
    
    def forward(self, node_features, factor_features, edge_index):
        """
        node_features: [num_nodes, node_dim]
        factor_features: [num_factors, factor_dim]
        edge_index: [2, num_edges]，变量-因子连接
        """
        # 假设 edge_index[0] 是变量索引，edge_index[1] 是因子索引
        var_idx, fac_idx = edge_index[0], edge_index[1]
        
        # 1. 变量 -> 因子消息
        var_msgs = self.node_to_factor_msg(node_features[var_idx])
        # 按因子聚合
        fac_incoming = scatter_add(var_msgs, fac_idx, dim=0, dim_size=factor_features.shape[0])
        
        # 2. 因子 -> 变量消息
        combined = torch.cat([fac_incoming[fac_idx], factor_features[fac_idx]], dim=-1)
        node_incoming = self.factor_to_node_msg(combined)
        
        # 3. 按变量聚合消息
        aggregated = self.aggregation(node_incoming, var_idx, num_nodes=node_features.shape[0])
        
        # 4. 更新变量表示
        updated_features = self.update_fn(aggregated, node_features)
        
        return updated_features

6.5 神经消息传递与变分推断的统一

近年来，研究者提出了神经变分推断（Neural Variational Inference）框架，将变分推断中的消息传递与神经网络的消息传递统一起来。

关键思想：

用神经网络参数化变分分布 $q_{ϕ} (X ∣ Y)$
用消息传递网络实现变分参数的更新
通过反向传播学习消息函数

7 代码实现

7.1 完整实现：离散变量的和-积算法

#include <bits/stdc++.h>
using namespace std;
 
/**
 * 离散变量和-积算法实现
 * 支持任意结构的因子图（树结构和带环图）
 */
 
struct FactorGraph {
    // 变量信息
    struct Variable {
        string name;
        int domain_size;                    // 取值范围 {0, 1, ..., domain_size-1}
        vector<int> neighbor_factors;       // 连接的因子索引
        
        Variable(string n, int d) : name(n), domain_size(d) {}
    };
    
    // 因子信息
    struct Factor {
        string name;
        vector<int> neighbor_vars;          // 连接的变量索引
        vector<double> potential;            // 势函数值，索引对应变量取值组合
        
        Factor(string n) : name(n) {}
    };
    
    vector<Variable> variables;
    vector<Factor> factors;
    
    // 添加变量
    int add_variable(const string& name, int domain_size) {
        variables.emplace_back(name, domain_size);
        return variables.size() - 1;
    }
    
    // 添加因子
    int add_factor(const string& name, const vector<int>& var_indices, 
                   const vector<double>& potential) {
        int idx = factors.size();
        factors.emplace_back(name);
        factors.back().neighbor_vars = var_indices;
        factors.back().potential = potential;
        
        // 更新变量的邻居信息
        for (int vi : var_indices) {
            variables[vi].neighbor_factors.push_back(idx);
        }
        
        return idx;
    }
    
    // 计算势函数的索引
    int potential_index(const Factor& f, const vector<int>& values) const {
        int idx = 0, multiplier = 1;
        for (int i = 0; i < f.neighbor_vars.size(); i++) {
            int vi = f.neighbor_vars[i];
            idx += values[vi] * multiplier;
            multiplier *= variables[vi].domain_size;
        }
        return idx;
    }
};
 
/**
 * 和-积算法
 */
class SumProductAlgorithm {
private:
    FactorGraph& graph;
    vector<vector<double>> var_to_factor_msgs;  // [var_idx][factor_idx]
    vector<vector<double>> factor_to_var_msgs;   // [factor_idx][var_idx_local]
    int num_iterations;
    double convergence_threshold;
    
    // 初始化消息
    void initialize_messages() {
        int V = graph.variables.size();
        int F = graph.factors.size();
        
        var_to_factor_msgs.assign(V, vector<double>(F, 1.0));
        factor_to_var_msgs.assign(F, vector<double>());
        
        for (int fi = 0; fi < F; fi++) {
            int num_neighbors = graph.factors[fi].neighbor_vars.size();
            factor_to_var_msgs[fi].assign(num_neighbors, 1.0);
        }
    }
    
    // 计算变量到因子的消息
    vector<double> compute_var_to_factor_msg(int vi, int fi) {
        const auto& var = graph.variables[vi];
        vector<double> msg(var.domain_size, 1.0);
        
        for (int neighbor_fi : var.neighbor_factors) {
            if (neighbor_fi == fi) continue;
            
            // 找到该因子在变量消息中的位置
            const auto& factor = graph.factors[neighbor_fi];
            auto it = find(factor.neighbor_vars.begin(), factor.neighbor_vars.end(), vi);
            int local_idx = it - factor.neighbor_vars.begin();
            
            // 消息乘积
            for (int val = 0; val < var.domain_size; val++) {
                msg[val] *= factor_to_var_msgs[neighbor_fi][local_idx];
            }
        }
        
        return msg;
    }
    
    // 计算因子到变量的消息
    vector<double> compute_factor_to_var_msg(int fi, int vi_local) {
        const auto& factor = graph.factors[fi];
        int num_vars = factor.neighbor_vars.size();
        
        // 构建完整的变量取值组合
        vector<int> domain_sizes(num_vars);
        vector<int> current_values(num_vars, 0);
        
        for (int i = 0; i < num_vars; i++) {
            int vi = factor.neighbor_vars[i];
            domain_sizes[i] = graph.variables[vi].domain_size;
        }
        
        // 边缘化
        vector<double> result(domain_sizes[vi_local], 0.0);
        
        // 遍历所有变量取值组合
        auto iterate = [&](auto&& self, int var_idx) -> void {
            if (var_idx == num_vars) {
                // 计算势函数值和输入消息的乘积
                double product = factor.potential[graph.potential_index(factor, current_values)];
                
                for (int i = 0; i < num_vars; i++) {
                    if (i == vi_local) continue;
                    product *= factor_to_var_msgs[fi][i];
                }
                
                result[current_values[vi_local]] += product;
                return;
            }
            
            for (int v = 0; v < domain_sizes[var_idx]; v++) {
                current_values[var_idx] = v;
                self(self, var_idx + 1);
            }
        };
        
        iterate(iterate, 0);
        return result;
    }
    
public:
    SumProductAlgorithm(FactorGraph& g, int max_iter = 100, double threshold = 1e-6)
        : graph(g), num_iterations(max_iter), convergence_threshold(threshold) {
        initialize_messages();
    }
    
    // 运行消息传递
    vector<vector<double>> run() {
        for (int iter = 0; iter < num_iterations; iter++) {
            double max_change = 0.0;
            
            // 1. 因子 -> 变量消息更新
            for (int fi = 0; fi < graph.factors.size(); fi++) {
                for (int li = 0; li < graph.factors[fi].neighbor_vars.size(); li++) {
                    int vi = graph.factors[fi].neighbor_vars[li];
                    auto new_msg = compute_factor_to_var_msg(fi, li);
                    
                    double change = 0.0;
                    for (int v = 0; v < new_msg.size(); v++) {
                        change = max(change, abs(new_msg[v] - factor_to_var_msgs[fi][li]));
                    }
                    max_change = max(max_change, change);
                    
                    factor_to_var_msgs[fi][li] = new_msg;
                }
            }
            
            // 2. 变量 -> 因子消息更新
            for (int vi = 0; vi < graph.variables.size(); vi++) {
                for (int fi : graph.variables[vi].neighbor_factors) {
                    auto new_msg = compute_var_to_factor_msg(vi, fi);
                    var_to_factor_msgs[vi][fi] = new_msg;
                }
            }
            
            // 检查收敛
            if (max_change < convergence_threshold) {
                cout << "Converged at iteration " << iter + 1 << endl;
                break;
            }
        }
        
        return get_marginals();
    }
    
    // 获取边缘分布
    vector<vector<double>> get_marginals() {
        vector<vector<double>> marginals;
        
        for (int vi = 0; vi < graph.variables.size(); vi++) {
            const auto& var = graph.variables[vi];
            vector<double> unnorm(var.domain_size, 1.0);
            
            // 乘积所有输入消息
            for (int fi : var.neighbor_factors) {
                auto it = find(graph.factors[fi].neighbor_vars.begin(),
                               graph.factors[fi].neighbor_vars.end(), vi);
                int local_idx = it - graph.factors[fi].neighbor_vars.begin();
                
                for (int v = 0; v < var.domain_size; v++) {
                    unnorm[v] *= factor_to_var_msgs[fi][local_idx];
                }
            }
            
            // 归一化
            double Z = accumulate(unnorm.begin(), unnorm.end(), 0.0);
            for (double& p : unnorm) p /= Z;
            
            marginals.push_back(unnorm);
        }
        
        return marginals;
    }
};
 
// 示例：图像去噪
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    
    FactorGraph graph;
    
    // 添加变量：3个像素
    int x0 = graph.add_variable("x0", 2);
    int x1 = graph.add_variable("x1", 2);
    int x2 = graph.add_variable("x2", 2);
    
    // 添加一元势函数（数据项）
    // phi(x_i) = p(y_i | x_i)
    vector<double> phi0 = {0.9, 0.1};  // P(x=0)高, P(x=1)低
    vector<double> phi1 = {0.5, 0.5};  // 不确定
    vector<double> phi2 = {0.2, 0.8};  // P(x=0)低, P(x=1)高
    
    graph.add_factor("phi0", {x0}, phi0);
    graph.add_factor("phi1", {x1}, phi1);
    graph.add_factor("phi2", {x2}, phi2);
    
    // 添加二元势函数（平滑项）
    // psi(x_i, x_j) = exp(beta * 1[x_i == x_j])
    vector<double> psi01 = {2.0, 0.5, 0.5, 2.0};  // [x0=0,x1=0], [x0=0,x1=1], [x0=1,x1=0], [x0=1,x1=1]
    vector<double> psi12 = {2.0, 0.5, 0.5, 2.0};
    
    graph.add_factor("psi01", {x0, x1}, psi01);
    graph.add_factor("psi12", {x1, x2}, psi12);
    
    // 运行和-积算法
    SumProductAlgorithm solver(graph, 100, 1e-10);
    auto marginals = solver.run();
    
    // 输出结果
    cout << "Marginal distributions:" << endl;
    for (int i = 0; i < marginals.size(); i++) {
        cout << "P(x" << i << "=0) = " << marginals[i][0] << ", P(x" << i << "=1) = " << marginals[i][1] << endl;
    }
    
    return 0;
}

7.2 PyTorch 实现：可学习消息传递

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.utils import scatter
 
 
class LearnableMessagePassing(nn.Module):
    """
    可学习的消息传递层
    结合因子图结构和神经网络的可学习性
    """
    
    def __init__(
        self,
        node_dim: int,
        edge_dim: int = 0,
        message_dim: int = 64,
        num_heads: int = 4,
        aggregation: str = 'sum'  # 'sum', 'mean', 'attention'
    ):
        super().__init__()
        
        self.node_dim = node_dim
        self.message_dim = message_dim
        self.num_heads = num_heads
        self.aggregation = aggregation
        
        # 消息函数网络
        self.message_net = nn.Sequential(
            nn.Linear(node_dim * 2 + edge_dim, message_dim),
            nn.LayerNorm(message_dim),
            nn.ReLU(),
            nn.Linear(message_dim, message_dim)
        )
        
        # 注意力机制（可选）
        if aggregation == 'attention':
            self.attention_net = nn.Sequential(
                nn.Linear(message_dim, 1)
            )
        
        # 更新函数
        self.update_net = nn.Sequential(
            nn.Linear(node_dim + message_dim, node_dim),
            nn.LayerNorm(node_dim)
        )
        
        # 可选的跳跃连接
        self.skip_weight = nn.Parameter(torch.tensor(0.5))
    
    def forward(
        self,
        node_features: torch.Tensor,        # [num_nodes, node_dim]
        edge_index: torch.Tensor,          # [2, num_edges]
        edge_features: torch.Tensor = None # [num_edges, edge_dim]
    ) -> torch.Tensor:
        """
        消息传递前向传播
        
        Args:
            node_features: 节点特征矩阵
            edge_index: 边索引 [2, num_edges]
            edge_features: 边特征（可选）
            
        Returns:
            更新后的节点特征
        """
        num_nodes = node_features.shape[0]
        src_idx, tgt_idx = edge_index[0], edge_index[1]
        
        # 1. 计算消息
        # 消息从源节点发送到目标节点
        if edge_features is not None:
            msg_input = torch.cat([
                node_features[src_idx],
                node_features[tgt_idx],
                edge_features
            ], dim=-1)
        else:
            msg_input = torch.cat([
                node_features[src_idx],
                node_features[tgt_idx]
            ], dim=-1)
        
        messages = self.message_net(msg_input)  # [num_edges, message_dim]
        
        # 2. 聚合消息
        if self.aggregation == 'sum':
            aggregated = scatter(
                messages, tgt_idx, dim=0, dim_size=num_nodes, reduce='sum'
            )
        elif self.aggregation == 'mean':
            aggregated = scatter(
                messages, tgt_idx, dim=0, dim_size=num_nodes, reduce='mean'
            )
        elif self.aggregation == 'attention':
            # 计算注意力权重
            attention_scores = self.attention_net(messages)  # [num_edges, 1]
            attention_weights = F.softmax(
                scatter(attention_scores, tgt_idx, dim=0, dim_size=num_nodes, reduce='max'),
                dim=0
            )
            # 重新索引注意力权重
            attention_weights = attention_weights[tgt_idx]
            # 加权聚合
            weighted_messages = messages * attention_weights
            aggregated = scatter(
                weighted_messages, tgt_idx, dim=0, dim_size=num_nodes, reduce='sum'
            )
        
        # 3. 更新节点表示
        update_input = torch.cat([node_features, aggregated], dim=-1)
        updated_features = self.update_net(update_input)
        
        # 4. 跳跃连接
        updated_features = (1 - self.skip_weight) * node_features + self.skip_weight * updated_features
        
        return updated_features
 
 
class MultiLayerMessagePassing(nn.Module):
    """
    多层消息传递网络（类似GraphSAGE/GIN）
    """
    
    def __init__(
        self,
        node_dim: int,
        hidden_dim: int = 64,
        num_layers: int = 3,
        dropout: float = 0.1
    ):
        super().__init__()
        
        self.layers = nn.ModuleList()
        current_dim = node_dim
        
        for i in range(num_layers):
            next_dim = hidden_dim if i < num_layers - 1 else node_dim
            self.layers.append(
                LearnableMessagePassing(
                    node_dim=current_dim,
                    message_dim=hidden_dim,
                    aggregation='sum'
                )
            )
            current_dim = next_dim
        
        self.dropout = nn.Dropout(dropout)
        self.layer_norms = nn.ModuleList([
            nn.LayerNorm(hidden_dim) for _ in range(num_layers - 1)
        ])
    
    def forward(self, node_features, edge_index):
        initial_features = node_features
        
        for i, layer in enumerate(self.layers):
            node_features = layer(node_features, edge_index)
            
            if i < len(self.layers) - 1:
                node_features = self.layer_norms[i](node_features)
                node_features = F.relu(node_features)
                node_features = self.dropout(node_features)
        
        # 残差连接
        return node_features + initial_features
 
 
# 示例用法
if __name__ == "__main__":
    # 模拟数据
    num_nodes = 10
    node_dim = 32
    
    node_features = torch.randn(num_nodes, node_dim)
    edge_index = torch.randint(0, num_nodes, (2, 30))  # 随机生成30条边
    
    # 创建模型
    model = MultiLayerMessagePassing(
        node_dim=node_dim,
        hidden_dim=64,
        num_layers=3
    )
    
    # 前向传播
    output = model(node_features, edge_index)
    print(f"Input shape: {node_features.shape}")
    print(f"Output shape: {output.shape}")

8 现代应用场景

8.1 计算机视觉

图像分割与语义标注：因子图用于建模像素之间的关系，一元势函数基于 CNN 特征，二元势函数编码平滑性约束。CRF（条件随机场）后处理是典型应用。

目标检测中的后处理：非极大值抑制（NMS）可以理解为基于置信度的消息传递，soft-NMS 是其软化版本。

8.2 自然语言处理

句法分析与依存关系：句子的句法结构可以用因子图表示，因子包括一元标签势函数（POS标签）和二元依存关系势函数。

序列标注：BiLSTM-CRF 模型使用 CRF 层进行序列级别的推断，本质上是链式因子图的解码。

8.3 推荐系统

协同过滤：矩阵分解可以解释为用户-物品因子图上的消息传递，用户和物品作为变量，评分作为因子。

图神经网络推荐：PIN-SAGE 等方法使用随机游走和消息传递学习用户/物品嵌入。

8.4 概率编程

PyMC3/JAXProb：自动消息传递框架，将概率模型编译为高效的消息传递算法。

变分自编码器：重参数化技巧使变分推断可以端到端学习。

9 总结与展望

因子图和消息传递算法是概率推断的核心框架，其思想贯穿于现代人工智能的多个领域：

理论基础：提供了从联合分布到边缘分布的高效计算方法
算法统一：和-积、max-积、变分消息传递是同一范式的不同变体
深度学习联系：GNN、注意力机制都可以理解为消息传递的实例
工程实践：在计算机视觉、NLP、推荐系统等领域有广泛应用

未来方向：

消息传递与 Transformer 的进一步统一
可扩展的近似消息传递算法
消息传递与神经符号计算的结合

参考文献

Kschischang, F. R., Frey, B. J., & Loeliger, H. A. (2001). Factor graphs and the sum-product algorithm. IEEE Transactions on Information Theory, 47(2), 498-519. ↩

Metaphor

探索

因子图与消息传递算法

因子图与消息传递算法

1 引言

1.1 为什么要用因子图？

1.2 消息传递的统一视角

2 因子图基础

2.1 形式化定义

2.2 树结构因子图

2.3 从因子图到因子节点图的映射

2.4 实例：图像去噪模型

3 和-积算法（Sum-Product Algorithm）

3.1 问题表述

3.2 消息传递规则的推导

3.2.1 从变量到因子的消息

3.2.2 从因子到变量的消息

3.2.3 消息计算的终止条件

3.3 算法伪代码

3.4 树结构上的消息传递顺序

3.5 连续变量的扩展

3.6 环状因子图与近似方法

4 Max-积算法（Max-Product Algorithm）

4.1 问题与和-积算法的联系

4.2 最大后验推断

4.3 追踪最优解：反向指针

4.4 数值稳定性

5 变分消息传递（Variational Message Passing）

5.1 变分推断框架

5.2 平均场近似

5.3 坐标上升变分推断（CAVI）

5.4 期望传播（Expectation Propagation）

5.5 变分消息传递与和-积算法的联系

6 与神经网络的联系

6.1 消息传递作为神经网络层

6.2 门控消息传递

6.3 注意力机制

6.4 可微分消息传递

6.5 神经消息传递与变分推断的统一

7 代码实现

7.1 完整实现：离散变量的和-积算法

7.2 PyTorch 实现：可学习消息传递

8 现代应用场景

8.1 计算机视觉

8.2 自然语言处理

8.3 推荐系统

8.4 概率编程

9 总结与展望

参考文献

Footnotes

关系图谱

目录

反向链接