所有变量都是离散的无向马尔科夫网络的使用很普遍，其中最常用的是二值变量的两两马尔科夫网络（pairwise Markove network）。在统计力学领域有时被称为 伊辛模型（Ising model），在机器学习领域则被称为 玻尔兹曼机（Boltzmann machine），其中二值的顶点被称作“节点（node）”或“单元（unit）”。
另外，每个节点的取值可能是观测到的（“可见的”）或没有观测到的（“隐藏的”）。通常这些节点会被组织为层（layer）的结构，这与神经网络类似。玻尔兹曼机在监督学习和无监督学习中都有很多应用，特别是对图片这样结构化的输入数据上。但是它们在计算上的困难限制了它们的应用。
图 17.6 展示了一个受限玻尔兹曼机（稍后会介绍），其中一部分变量是隐藏的，并且只有一部分节点之间是连接的。我们先从一个简单的情况入手，其中所有 $p$ 个节点都是可见的，$E$ 是边 $(j,k)$ 的集合。
将节点 $j$ 处的二值变量记为 $X_j$，那么联合分布的伊辛模型可写为：
$$p(X, \mathbf{\Theta}) = \exp \left[ \sum_{(j,k)\in E} \theta_{jk} X_j X_k - \Phi(\mathbf{\Theta}) \right] \text{ for } X \in \mathcal{X} \tag{17.28}$$其中的 $\mathcal{X}=\{0,1\}^p$。与上一节中的高斯模型类似，这里只考虑变量两两（pairwise）之间的交互。伊辛模型是在统计力学中发展来的，现在被更广泛地用作两两交互的联合分布的模型。$\Phi(\mathbf{\Theta})$ 是配分函数（partition function）的对数，定义为：
$$\Phi(\mathbf{\Theta}) = \log \sum_{x\in\mathcal{X}} \left[ \exp \left( \sum_{(j,k) \in E} \theta_{jk} x_j x_k \right) \right] \tag{17.29}$$配分函数保证了样本空间上的概率和为一。其中的 $\theta_{jk}X_jX_k$ 项代笔了（对数）势函数（potential function，式 17.5）的一个特定的参数化。由于技术上的原因，模型中需要包括一个“常数”节点 $X_0\equiv1$（练习 17.10），它与其他所有节点之间都有边。在统计学中，这个模型等价于多维计数表的一个一阶交互泊松对数线性模型（Bishop et al., 1975；McCullagh and Nelder, 1989；Agresti, 2002）。
伊辛模型假设了给定其他节点后每个节点的条件概率是一个对数几率的形式：
$$\operatorname{Pr}(X_j=1|X_{-j}=x_{-j}) = \frac{1} {1 + \exp(-\theta_{j0} - \sum_{(j,k) \in E} \theta_{jk} x_k)}$$ $$\tag{17.30}$$其中的 $X_{-j}$ 代表着除 $j$ 外的所有节点。所以参数 $\theta_{jk}$ 衡量了给定其他节点后 $X_j$ 对 $X_k$ 的依赖程度。
17.4.1 已知图结构时的参数估计
给定这个模型的一个数据集后，如何估计它的参数？假设观测样本点 $x_i=(x_{i1},x_{i2},\dots,x_{ip})\in\{0,1\}^p$，$i=1,\dots,N$。则对数似然函数是：
$$\begin{align} \ell(\mathbf{\Theta}) &= \sum_{i=1}^N \log \operatorname{Pr}\_\mathbf{\Theta} (X_i=x_i) \\ &= \sum_{i=1}^N \left[ \sum_{(j,k) \in E} \theta_{jk} x_{ij} x_{ik} - \Phi(\mathbf{\Theta}) \right] \tag{17.31}\end{align}$$对数似然函数的梯度函数是：
$$\frac{\partial \ell(\mathbf{\Theta})}{\partial \theta_{jk}} = \sum_{i=1}^N x_{ij} x_{ik} - N \frac{\partial \Phi(\mathbf{\Theta})}{\partial \theta_{jk}} \tag{17.32}$$并有：
$$\begin{align} \frac{\partial \Phi(\mathbf{\Theta})}{\partial \theta_{jk}} &= \sum_{x \in \mathcal{X}} x_j x_k \cdot p(x, \mathbf{\Theta}) \\ &= \operatorname{E}_{\mathbf{\Theta}} (X_j X_k) \tag{17.33}\end{align}$$将梯度设为零，可得：
$$\hat{\operatorname{E}} (X_j X_k) - \operatorname{E}_{\mathbf{\Theta}} (X_j X_k) = 0 \tag{17.34}$$其中的第一项定义为：
$$\hat{\operatorname{E}} (X_j X_k) = \frac{1}{N}\sum_{i=1}^N x_{ij} x_{ik} \tag{17.35}$$其中的期望是取在数据的经验概率分布上的。从式 17.34 中可见在最大似然估计中，节点之间（期望）内积的估计就是它们在观测样本上的（平均）内积。这是指数分布族的模型的得分（score）或梯度（gradient）方程式的标准形式，在指数分布模型中充分统计量（sufficient statistics）就等于统计量在模型假设下的期望。
我们可以使用梯度搜索或牛顿方法来求解最大似然的估计。然而对 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$ 的计算需要将 $p(X,\mathbf{\Theta})$ 在 $X$ 的 $|\mathcal{X}|=2^p$ 个可能取值中的 $2^{p-2}$ 个取值上遍历计算一遍，一般来说当 $p$ 比较大时（例如大于 30）这个计算不现实。在 $p$ 不大时，已有一些标注的统计计算方法：
- 泊松对数线性模型（Poisson log-linear modeling）
将问题当成很多的回归问题来解决（练习 17.12）。输出向量 $\mathbf{y}$ 是长度为 $2^p$ 的向量，即数据多维表格的每个单元格内的值1。自变量矩阵 $\mathbf{Z}$ 有 $2^p$ 行，并且最多有 $1+p+p^2$ 个列，描述了每一个单元格；不过具体的数字还取决于图模型的稀疏程度。它的计算量基本等同于一个同规模回归问题的计算量，等于 $O(p^42^p)$，当 $p<20$ 时是可以接受的。 一般会通过迭代重加权最小二乘（iterative reweighted least square，IRLS）来计算牛顿更新步骤，通常迭代的步骤数量不大于 10。更多细节可参考 Agresti (2002) 以及 McCullagh and Nelder (1989)。可以使用通用软件（例如 R 扩展包glm）来拟合这个模型 - 梯度下降（Gradient descent）
这个方法中梯度的计算最多需要 $O(p^22^{p-2})$ 的计算量，但可能会比二阶牛顿方法需要更多的梯度步骤。它仍然可以处理稍微高维度的问题，$p\leq30$。利用稀疏图中特殊的团结构，以及联合树（junction tree）算法，可以进一步缩减这些计算量。在这里不做更多介绍。 - 迭代比例拟合（Iterative proportional fitting，IPF）
对梯度方程式 17.34 使用了循环坐标下降（cyclical coordinat descent）方法。在每个步骤中，只更新一个参数并使得梯度方程式等于零。循环地进行这个过程，直到所有的梯度方程式都等于零。一次完全的循环的计算量消耗与一次梯度计算相同，但有使其更高效的可能性。Jirouśek and Přeučil (1995) 利用联合树模型提出了一个更高效版本的 IPF。
当 $p$ 比较大时（大于 30），有一些其他的方法来近似梯度。
- 根据平均场近似（mean field approximation），用 $\operatorname{E}_\mathbf{\Theta}(X_j)\operatorname{E}_\mathbf{\Theta}(X_k)$ 来估计 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$（Peterson and Anderson, 1987），并用输入变量的均值来替代这些输入变量，这样就得到了参数 $\theta_{jk}$ 的一个非线性方程组。
- 使用吉布斯采样可以获得近精确的解（第 8.6 节），通过从估计的模型概率 $\operatorname{Pr}_\mathbf{\Theta}(X_j|X_{-j})$ 中的连续采样来近似出 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$（参考 Ripley (1996)）。
这里没有介绍 可分解模型（decomposable model），它们可以得出最大似然估计的闭式解，而不需要任何的迭代过程。树就是这种模型的一个例子，即带有树结构拓扑的特殊的图。在关系到计算可行性的问题中，树结构代表了一类有效的模型，它们可以绕过这些场景中出现的计算量障碍。更多细节可以参考 Whittaker (1990) 的第十二章。
17.4.2 隐藏节点 😱
可以通过加入潜在或隐藏节点来增加一个离散马尔科夫网络的复杂程度。假设输入变量中的一个子集 $X_\mathcal{H}$ 是没有观测到的，或“隐藏的”，其他的变量 $X_\mathcal{V}$ 是观测到的，或“可见的”。那么观测数据的对数似然函数为：
$$\begin{align} \ell(\mathbf{\Theta}) &= \sum_{i=1}^N \log [\operatorname{Pr}_\mathbf{\Theta} (X_{\mathcal{V}}=x_{i\mathcal{V}})] \\ &= \sum_{i=1}^N \left[ \log \sum_{x_\mathcal{H} \in \mathcal{X}_\mathcal{H}} \exp \sum_{(j,k) \in E} (\theta_{jk} x_{ij} x_{ik} - \Phi(\mathbf{\Theta})) \right] \tag{17.36}\end{align}$$对 $x_\mathcal{H}$ 的求和是在对隐藏单元的所有可能的 $\{0,1\}$ 取值上进行求和。梯度函数为：
$$\frac{d \ell(\mathbf{\Theta})}{d \theta_{jk}} = \hat{\operatorname{E}}_\mathcal{V}\operatorname{E}_\mathbf{\Theta} (X_j X_k | X_\mathcal{V}) - \operatorname{E}_\mathbf{\Theta} (X_j X_k) \tag{17.37}$$如果 $X_j$ 和 $X_k$ 都是可见的，那么第一项就是 $X_jX_k$ 的经验分布平均；如果其中之一或者两者都是隐藏的，则首先要用可见数据对隐藏变量进行插补（imputed），然后在隐藏变量上进行平均。第二项是 $X_jX_k$ 的无条件期望。
第一项中内部的期望可以使用条件期望的基本规则以及伯努利（Bernoulli）随机变量的性质来计算取值。具体来说，对观测样本 $i$：
$$\begin{align} & \operatorname{E}_\mathbf{\Theta} (X_j X_k | X_\mathcal{V}=x_{i\mathcal{V}}) \\ =& \begin{cases} x_{ij} x_{ik} & \text{如果 } j,k \in \mathcal{V} \\ x_{ij} \operatorname{Pr}_\mathbf{\Theta} (X_k=1 | X_\mathcal{V} = x_{i\mathcal{V}}) & \text{如果 } j \in \mathcal{V}, k \in \mathcal{H} \\ \operatorname{Pr}_\mathbf{\Theta} (X_j=1, X_k=1 | X_\mathcal{V} = x_{i\mathcal{V}}) & \text{如果 } j,k \in \mathcal{H} \end{cases} \tag{17.38}\end{align}$$这样就需要分别进行两次吉布斯采样。第一次采样通过从上面模型的采样来估计 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$，第二次采样来估计 $\operatorname{E}_\mathbf{\Theta}(X_jX_k|X_\mathcal{V}=x_{i\mathcal{V}})$。在后者的第二次采样中，可见单元被固定在了它们的观测值上，只对隐藏变量进行采样。在梯度搜寻的每个步骤中，需要对训练集中的每个观测样本点分别进行吉布斯采样。结果就导致了这个计算过程即使是在一个中等大小的模型上也可能很慢。第 17.4.4 节 中会给模型添加更多约束条件，从而使这个计算过程更可行。
17.4.3 图结构的估计
Lee et al. (2007) 和 Wainwright et al. (2007) 已经提出过在二值两两马尔科夫网络中对套索惩罚项的使用。第一个文章的作者研究了通过一个共轭梯度（conjugate gradient）过程对一个带惩罚项的对数似然函数的精确最大化求解。它的瓶颈是对梯度中 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$ 的计算。对于稀疏的图来说，通过联合树算法的精确计算还是可行的，但是对于密集的图就变得不再可行了。
第二个文章的作者提出了一个近似解，它与 Meinshausen and Bühlmann (2006) 求解高斯图模型的方法相类似。他们为每个节点对其他节点进行了一个 $L_1$ 惩罚的对数几率回归，然后用某些方式将边参数进行对称化处理。例如，假设 $\tilde{\theta}_{jk}$ 是以节点 $j$ 为输出变量的对数几率模型中得出的对 $(j,k)$ 边的参数的估计，那么“最小值”对称化规则就令 $\hat{\theta}_{jk}$ 等于 $\tilde{\theta}_{jk}$ 和 $\tilde{\theta}_{kj}$ 两者中绝对值最小的那个。相似地，也可以定义出“最大值”对称化规则。他们证明了在某些条件下、随着样本量趋于无穷大，任一近似方法都可以正确地估计出非零的边。Hoefling and Tibshirani (2008) 将图套索（graphical lasso）扩展到了离散马尔科夫网络中，得出了一个从某些程度上比共轭梯度方法更快的计算过程，但仍然需要面对 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$ 的计算。他们也通过大量的模拟研究对比了精确解和近似解；并发现在估计非零边以及估计边参数的值时，“最小值”或“最大值”近似的准确性只会比精确计算过程的稍差一点，但却要快得多。此外，这种近似方法可以用在密集的图上，因为它不需要计算 $\operatorname{E}_\mathbf{\Theta}(X_jX_k)$。
最后，我们说明一个高斯模型和二值模型之间的一个重要的区别。在高斯模型中，通常需要同时估计出 $\mathbf{\Sigma}$ 和它的逆矩阵，图套索计算过程可以得出这两者的估计。不过 Meinshausen and Bühlmann (2006) 在高斯图模型中的近似，与 Wainwright et al. (2007) 在二值模型中的近似相类似，都只得出了 $\mathbf{\Sigma}^{-1}$ 的估计。与之相反，在二值数据的马尔科夫模型中，只需要估计出 $\mathbf{\Theta}$，而不需要估计它的逆矩阵。Wainwright et al. (2007) 中的近似方法可以高效地估计出 $\mathbf{\Theta}$，因此是二值问题中一个有吸引力的求解方法。
17.4.4 受限玻尔兹曼机
本节介绍一个受神经网络所启发的一个图模型的特别的结构，其中的单元被组织到层（layer）的结构中。一个 受限玻尔兹曼机（restricted Boltzmann machine，RBM） 包含了一层可见的单元和一层隐藏的单元，在每层之内都没有连接。假设隐藏单元之间没有连接，那么条件期望（类似式 17.37 和式 17.38）的计算就变得简单得多了2。
图 17.6 展示了一个例子。可见层被分成了输入变量 $\mathcal{V}_1$ 和输出变量 $\mathcal{V}_2$ 两个部分，并有一个隐藏层 $\mathcal{H}$。这个模型可写为：
$$\mathcal{V}_1 \leftrightarrow \mathcal{H} \leftrightarrow \mathcal{V}_2 \tag{17.39}$$例如，$\mathcal{V}_1$ 可能是一个手写数字图片的二值像素点，而 $\mathcal{V}_2$ 可能是十个单元，每个单元对应着观测到的 0~9 的类别标签。
由于在给定了其他层的变量后，每个层内部的变量彼此独立，所以这个模型的受限形式简化了用于估计式 17.37 中的期望的吉布斯采样。可以使用表达式 17.30 中的条件概率对他们同时采样。
这样（受限）的模型没有玻尔兹曼机的通用性强，但仍然很实用；例如它可以从图片中提取出一些有趣的特征变量。在图 17.6 所示的受限玻尔兹曼机（RBM）中，交替地对每层中变量进行采样，就可以生成来自联合概率密度模型中的样本。如果在交替采样中 $\mathcal{V}_1$ 部分的可见层被固定在了一个特定的特征向量取值，就可以从给定了 $\mathcal{V}_1$ 后类别标签的条件分布进行采样。或者，也可以通过观测样本中每个类别标签和特征变量的非标准化的联合分布密度之间的比较，来得出对测试样本的分类。我们不需要计算配分函数，因为它对所有这些组合来说都是一样的。
可见受限玻尔兹曼机在一般形态上与单隐藏层的神经网络（第 11.3 节）一样。神经网络中的边是由方向的，隐藏单元通常是实数取值的，并且拟合准则函数也不同。神经网络是对给定输入特征变量条件下目标变量与模型对它们的预测之间的误差（交叉熵）进行最小化。与之相反，受限玻尔兹曼机是对所有可见单元（即特征变量和目标变量）的联合分布的对数似然函数进行最大化。它可以从输入特征变量中提取出对分类标签预测有效的信息，但与监督学习方法不同的是，它也可以利用其中一些隐藏单元来为那些看似与标签预测不直接相关的特征向量来建模。不过当与从其他隐藏层中衍生出的特征结合起来之后，这些特征变量可能会是有效的。
不过可惜的是，受限玻尔兹曼机里面的吉布斯采样可能会非常慢，它可能需要花费很长的时间才能达到稳定状态。随着网络的权重增大，链条混合得更慢，为了得到无条件估计则需要花费更多的步骤。Hinton (2002) 在实验中发现，在估计式 17.37 中的第二项期望时，如果以样本数据作为马尔可夫链的起始点然后只进行少数的几步（而不是到收敛），那么训练的结果仍然很好。他将这个过程称为 对比发散（contrastive divergence）：先给定 $\mathcal{V}_1$ 和 $\mathcal{V}_2$ 后对 $\mathcal{H}$ 采样，然后给定 $\mathcal{H}$ 后对 $\mathcal{V}_1$ 和 $\mathcal{V}_2$ 采样，最后再次给定 $\mathcal{V}_1$ 和 $\mathcal{V}_2$ 后对 $\mathcal{H}$ 采样。它的想法是当参数与解的距离比较远时，将吉布斯采样迭代到稳态可能是一个浪费，然而只需要一次迭代就可以得到改进估计结果的一个不错的方向。
我们用一个例子来演示受限玻尔兹曼机（RBM）。利用对比发散，可以训练 RBM 模型来识别 MNIST 数据集（LeCun et al., 1998）中的手写数字。当使用了 2000 个隐藏单元、784 个代表二值的像素亮度的单元、和 1 个代表类别标签的 10 维多项单元，RBM 在测试集上可以达到 1.9% 的错误率。这比支持向量机的 1.4% 略高一些，与反向传播训练出的神经网络的错误率类似。不过如果用 500 个从图片中提取而没有用到任何标签信息的特征变量替换原来的 784 个像素亮度特征，RBM 的错误率可以被降低到 1.25%。首先，在图片数据集上通过对比发散方法训练一个有 784 个可见单元和 500 个隐藏单元的 RBM。然后将这第一个 RBM 中的隐藏单元作为输入数据，再训练第二个有 500 个可见单元和 500 个隐藏单元的 RBM。最后将第二个 RBM 中的隐藏单元作为特征变量，训练一个有 2000 个隐藏单元的 RBM 作为最后的联合分布密度模型。Hinton et al. (2006) 介绍了这个贪婪的、逐层的特征训练方式的细节和背后的原理。图 17.7 展示了以这种方式训练的合成模型的一个示意图，以及它所能够应对的各种类型的图片变形的一些示例。
本节练习
练习 17.10
Using a simple binary graphical model with just two variables, show why it is essential to include a constant node X0 ≡ 1 in the model.
练习 17.11
Show that the Ising model (17.28) for the joint probabilities in a discrete graphical model implies that the conditional distributions have the logistic form (17.30).
练习 17.12
Consider a Poisson regression problem with p binary variables xij , j = 1, . . . , p and response variable yi which measures the number of observations with predictor xi ∈ {0, 1}p . The design is balanced, in that all n = 2p possible combinations are measured. We assume a log-linear model for the Poisson mean in each cell
$$\log \mu(X) = \theta_{00} + \sum_{(j,k) \in E} x_{ij} x_{ik} \theta_{jk} \tag{17.45}$$using the same notation as in Section 17.4.1 (including the constant variable xi0 = 1∀i). We assume the response is distributed as
$$\operatorname{Pr}(Y=y | X=x) = \frac{e^{-\mu(x)} \mu(x)^y}{y!} \tag{17.46}$$Write down the conditional log-likelihood for the observed responses yi , and compute the gradient.
- Show that the gradient equation for θ00 computes the partition func- tion (17.29).
- Show that the gradient equations for the remainder of the parameters are equivalent to the gradient (17.34).