因子图

作者: Jacob Schreiber 联系方式: jmschreiber91@gmail.com

因子图是一种强大但在我看来未被充分利用的概率模型。其中一个潜在原因是它们在概念上不如贝叶斯网络或混合模型直观。在这里，我们将消除所有困惑。

因子图与贝叶斯网络类似，都由一组概率分布和连接它们的图结构组成。但不同于贝叶斯网络的是，因子图是二分图结构，其概率分布集合不仅包含数据集中的变量，而是由各变量的边缘分布与因子分布（通常为多元分布）共同组成。在图中，边缘分布位于一侧，因子分布位于另一侧，无向边仅连接因子分布与其构成变量的边缘分布。

一种理解方式是先从贝叶斯网络入手。如果将所有的条件概率分布转换为联合概率分布，并保持单变量分布不变，就得到了因子分布集。接着，为网络中的每个变量添加一个边缘分布，这就构成了因子图中的节点集。最后，在每个因子分布与组成该联合分布的变量对应的边缘分布之间添加边。当因子是单变量分布时，只需在该因子与其边缘分布之间添加一条边。

一旦因子图实例化后，推理过程涉及一个迭代的消息传递算法。第一步，边缘分布会发出消息，这些消息是它们自身的副本。第二步，因子分布接收每个边缘分布的估计值，将它们与联合概率参数结合，并向每个边缘分布发回对各变量的新估计值。最后在第三步，边缘分布接收这些估计值，将它们进行平均，并向各变量发回新的估计值。第二和第三步会重复执行，直到收敛为止。

注意：虽然因子图已实现参数拟合，但目前还不支持结构学习。

%pylab inline
import seaborn; seaborn.set_style('whitegrid')
import torch

%load_ext watermark
%watermark -m -n -p torch,pomegranate

Populating the interactive namespace from numpy and matplotlib
torch      : 1.13.0
pomegranate: 1.0.0

Compiler    : GCC 11.2.0
OS          : Linux
Release     : 4.15.0-208-generic
Machine     : x86_64
Processor   : x86_64
CPU cores   : 8
Architecture: 64bit

初始化与拟合

因子图由两组分布组成：因子分布（即联合概率与单变量概率的混合）和边缘分布（数据中每个变量对应一个单变量概率分布）。边缘分布通常设为均匀分布，而因子分布则编码了底层数据的统计信息。这两组分布通过一个无向无权的二分图相互连接，使得因子可以与边缘分布相连但不能与其他因子相连。数据中的每个变量都必须有一个对应的边缘分布，并至少出现在因子侧的某个分布中。每个变量可以出现在多个因子分布中。

让我们从实现最简单的因子图开始：两个变量通过联合分类分布连接在一起。

from pomegranate.distributions import Categorical
from pomegranate.distributions import JointCategorical
from pomegranate.factor_graph import FactorGraph

m1 = Categorical([[0.5, 0.5]])
m2 = Categorical([[0.5, 0.5]])

f1 = JointCategorical([[0.1, 0.3], [0.2, 0.4]])

model = FactorGraph([f1], [m1, m2], [(m1, f1), (m2, f1)])

请注意，向模型添加边的顺序表示变量在JointCategorical分布中的排序方式，即第一条包含f1的边（本例中就是第一条边）表明该边的父节点是因子张量中的第一个维度。

类似地，将边缘分布添加到模型中的顺序对应于数据中变量的顺序。具体来说，m1对应第一列，m2对应第二列，而f1则由第一列和第二列的数据组成。因子中变量的顺序不需要排序。如果边是以[(m2, f1), (m1, f1)]的形式添加的，将会生成一个有效的因子图，唯一的区别是f1张量的第一个维度将对应第二列数据。

如果您希望以更编程化的方式构建因子图，可能会倾向于使用add_edge、add_marginal和add_factor方法来逐步构建图。这些方法在初始化时传入这些值时会在内部使用，但也可以直接调用它们。

model = FactorGraph()
model.add_factor(f1)

model.add_marginal(m1)
model.add_marginal(m2)

model.add_edge(m1, f1)
model.add_edge(m2, f1)

如果您有数据，就可以像调整其他参数一样训练因子参数。

X = torch.randint(2, size=(11, 2))
X

tensor([[0, 0],
        [1, 0],
        [1, 0],
        [0, 0],
        [1, 1],
        [0, 1],
        [1, 0],
        [1, 1],
        [0, 0],
        [1, 1],
        [0, 1]])

model.fit(X)
model.factors[0].probs

Parameter containing:
tensor([[0.2727, 0.1818],
        [0.2727, 0.2727]])

重要的是，这些更新不会改变边缘分布值。这是因为因子值编码了你可以从先前数据中学到的所有信息（似然函数），而边缘分布表示你对该变量符号的先验概率。当你在知道某个变量的符号时进行推断，你会将边缘分布设置为该值的100%。

model.marginals[0].probs

Parameter containing:
tensor([[0.5000, 0.5000]])

概率与对数概率

与其他方法类似，可以通过给定模型计算每个样本的概率。但与其他模型不同的是，这个概率会按照因子和边缘分布进行分解，当数据包含f个因子和d个维度时，公式为\(P(X) = \prod\limits_{i=0}^{f} P(X_{p_i} | F_i) \prod\limits_{i=0}^{d} P(X_i | M_i)\)。本质上，您需要评估数据在模型似然下的概率（因子部分的乘积），再乘以数据在边缘分布下的概率（边缘部分的乘积）。在没有任何先验信息的情况下，边缘概率是常数可以忽略不计。

model.probability(X[:1])

tensor([0.0682])

请注意，上述值包含了与边缘分布的乘积运算，因此可能不会直接与因子中的任何值匹配。

model.log_probability(X)

tensor([-2.6856, -2.6856, -2.6856, -2.6856, -2.6856, -3.0910, -2.6856, -2.6856,
        -2.6856, -2.6856, -3.0910])

预测

与贝叶斯网络类似，因子图可以对数据集中的缺失值进行预测。实际上，贝叶斯网络和马尔可夫网络在后端经常构建因子图来进行实际的推理。这些方法使用和积算法（也称为循环信念传播）。该算法的基本工作原理如下：

• 初始化从每个边缘分布到每个因子的消息，这些消息是边缘分布的副本

• 对于每个因子，遍历其连接的边缘分布，并使用所有其他传入消息（忽略来自目标边缘分布的消息）计算该因子的边缘分布，然后将其对该分布应为何值的置信度发送回目标边缘分布

• 对于每个边缘分布，将所有传入的消息相乘，以估算新的边缘值应该是多少。

消息收敛后，将返回新的边缘分布，表示每个分布的概率

当已知某些变量的预期取值时，例如处理不完整矩阵的情况，可以设置初始边缘概率分布（不是消息传递，而是实际分布）与之保持一致（即在分类情况下，观测到2意味着将该类别的概率设为100%，而非使用均匀概率分布）。

通过使用torch.masked.MaskedTensor来指定不完整的矩阵。

X_torch = torch.tensor(X[:4])
mask = torch.tensor([[True, False],
                     [False, True],
                     [True, True],
                     [False, False]])

X_masked = torch.masked.MaskedTensor(X_torch, mask=mask)

/tmp/ipykernel_89475/916229149.py:1: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
  X_torch = torch.tensor(X[:4])
/home/jmschr/anaconda3/lib/python3.9/site-packages/torch/masked/maskedtensor/core.py:156: UserWarning: The PyTorch API of MaskedTensors is in prototype stage and will change in the near future. Please open a Github issue for features requests and see our documentation on the torch.masked module for further information about the project.
  warnings.warn(("The PyTorch API of MaskedTensors is in prototype stage "

掩码通过True值表示哪些索引是已知的。当掩码为False时，所取的数据值无关紧要。

model.predict(X_masked)

tensor([[0, 0],
        [0, 0],
        [1, 0],
        [1, 0]])

虽然计算了完整的概率分布，但这只是对sum-product算法结果的离散化展示。如果想获取完整的预测概率，可以这样做：

model.predict_proba(X_masked)

[tensor([[1.0000, 0.0000],
         [0.5000, 0.5000],
         [0.0000, 1.0000],
         [0.4545, 0.5455]]),
 tensor([[0.6000, 0.4000],
         [1.0000, 0.0000],
         [1.0000, 0.0000],
         [0.5455, 0.4545]])]

请注意，这里的输出是一个二维列表。列表中的每个张量对应基础数据的一个维度，当使用分类数据类型时，这些张量的形状为(n_examples, n_categories)。基本上，第一个张量中的第一行包含元素X_masked[0, 0]取值为0和1的概率。由于X_masked[0, 0]是一个观测值（即掩码值为True），这意味着在观测值处的概率被固定为1。

如果我们只想要对数概率，也可以计算这些值。

model.predict_log_proba(X_masked)

[tensor([[ 0.0000,    -inf],
         [-0.6931, -0.6931],
         [   -inf,  0.0000],
         [-0.7885, -0.6061]]),
 tensor([[-0.5108, -0.9163],
         [ 0.0000,    -inf],
         [ 0.0000,    -inf],
         [-0.6061, -0.7885]])]

总结

与其他模型类似，因子图也可以利用摘要API进行训练。基本上，您可以使用summarize方法汇总一个或多个数据批次，然后调用from_summaries方法来更新分布的参数。这将仅更新因子的参数，而保持边缘分布不变。

model.summarize(X[:3])
model.from_summaries()

然后，我们可以检查参数。

model.factors[0].probs

Parameter containing:
tensor([[0.3333, 0.0000],
        [0.6667, 0.0000]])