17.3 连续变量无向图模型

本节将讨论所有变量都是连续变量的马尔科夫网络。在这样的图模型中，一般总会使用高斯分布的假设，这是因为它比较方便的分析性质。假设观测样本服从一个均值为 $\mu$ 协方差矩阵为 $\mathbf{\Sigma}$ 的多元高斯分布。由于高斯分布中最高只有二阶的交互关系，所以它自然地符合一个两两高斯分布图。图 17.1 中的图模型就是一个高斯图模型的例子。

高斯分布有一个性质：所有的条件分布仍然是高斯分布。协方差矩阵的逆矩阵 $\mathbf{\Sigma}^{-1}$ 包含了关于变量之间 偏协方差（partial covariance） 的信息，即在所有其他随机变量的条件下，一对变量 $i$ 和 $j$ 之间的协方差。特别地，如果 $\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$ 的第 $ij$ 元素为零，那么变量 $i$ 和 $j$ 在给定其他变量的条件下是独立的（练习 17.3）。

考察给定其他变量后一个变量的条件分布有助于理解这个模型，在这里能够明显看出 $\mathbf{\Theta}$ 的作用。假设我们将随机向量划分为 $X=(Z,Y)$，其中的 $Z=(X_1,\dots,X_{p-1})$ 包含了前 $p-1$ 个变量，$Y=X_p$ 是最后一个变量。那么给定 $Z$ 后 $Y$ 的条件概率分布是（可参考 Mardia et al., 1979）：

其中我们将 $\mathbf{\Sigma}$ 划分成了：

式 17.6 中的条件期望与 $Y$ 在 $Z$ 上的多元线性回归的总体（population）表达式完全一致，回归的系数为 $\beta=\mathbf{\Sigma}_{ZZ}^{-1}\sigma_{ZY}$（参考 19 页的式 2.16）。如果将 $\mathbf{\Theta}$ 以同样的方式划分，并且因为 $\mathbf{\Sigma}\mathbf{\Theta}$，根据分块矩阵的逆矩阵公式可得：

并且 $1/\theta_{YY}=\sigma_{YY}-\sigma_{ZY}^T\mathbf{\Sigma}_{ZZ}^{-1}\sigma_{ZY}>0$。因此：

从这里可知两点：

• 式 17.6 中 $Y$ 对 $Z$ 的依赖（$Z=z$）只体现在均值项中。从上述推导中可看出，$\beta$ 中为零的元素也就是 $\theta_{ZY}$ 中为零的元素，其在 $Z$ 中对应的元素在给定其他变量的条件下与 $Y$ 独立。
• 条件分布中的依赖结构可以通过多元线性回归体现出来。

在描述每个节点在给定其他变量条件下的条件分布时，$\mathbf{\Theta}$ 包含了所需的所有结构性和数值的二阶（相关性）信息，因此它也被称为高斯图模型的“自然”参数¹。

协方差图（covariance graph） 或 相关性网络（relevance network） 是另一种不同类型的图模型，其中如果两个变量之间的协方差（而不是偏协方差）不为零，则对应的顶点之间由二向边连接。这类模型在基因组学（genomics）中比较常见，可参考 Butte et al. (2000)。这些模型的负对数几率函数是非凸的，导致它们的计算更加困难（Chaudhuri et al., 2007）。

17.3.1 已知图结构时的参数估计

给定 $X$ 的一些样本取值，我们想要估计出近似它们的联合概率分布的无向图中的参数。首先假设这个图是完全图（所有点之间都有连接）。并且假设 $x_i,i=1,\dots,N$ 为一个总体均值 $\mu$ 协方差 $\mathbf{\Sigma}$ 的多元高斯分布的 $N$ 个随机样本值。

式 17.10 是经验协方差矩阵，其中 $\bar{x}$ 是样本均值向量。忽略常数项，数据样本的对数几率函数可写为：

在式 17.11 中我们已经针对均值向量 $\mu$ 对函数进行了部分最大化。数值 $-\ell(\mathbf{\Theta})$ 是对 $\mathbf{\Theta}$ 的一个凸函数。易证明 $\mathbf{\Sigma}$ 的最大似然估计就是 $\mathbf{S}$。

为了使图模型（特别是在高维问题中）更实用，现假设一些顶点之间没有边；例如图 17.1 中 PIP3 和 Erk 之间，以及其他。如之前所述，对于高斯分布来说这意味着 $\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$ 中对应位置的元素为零。所以我们可以在对式 17.11 最大化时，添加某个预设子集中的参数为零的假设。这是一个等式约束的凸优化问题，已有一些求解它的方法了，比如迭代比例拟合过程（iterative proportional fitting procedure，IPFP，Speed and Kiiveri, 1986）。Whittaker (1990) 和 Lauritzen (1996) 概述了这个方法和其他一些方法。这些方法都是利用了前一节所介绍的将图分解成极大团从而得到的简化。我们在这里介绍另外一个简单的方法，它以另一个方式利用了模型的稀疏性。在后续对图结构的估计问题的讨论中，我们可以更清楚地看到这个方法的优势。

这个方法的思路是基于线性回归，它受式 17.6 和 17.9 所启发。具体来说，假设我们想要估计的是连接到某个顶点 $i$ 的顶点对应边的参数 $\theta_{ij}$，并且约束那些无连接对应的参数为零。那么顶点 $i$ 的值对其他相关顶点的值的线性回归看起来可能会得到一个合理的估计。但是这样就忽略了这个回归中自变量之间的相关依赖结构。我们发现如果在进行这个回归时使用了当前对自变量交叉积矩阵（基于模型）的估计，就可以得到正确的解并且正确地求解了带约束的最大似然问题。以下为具体细节。

在式 17.11 的对数几率函数中，添加所有无连接边的拉格朗日约束项：

式 17.12 最大化的梯度方程可写为：

这里用到了 $\log\det\mathbf{\Theta}$ 的导数等于 $\mathbf{\Theta}^{-1}$ 的结论（例如 Boyd and Vandenberghe, 2004，641 页)。$\mathbf{\Gamma}$ 是拉格朗日参数的矩阵，所有缺失边的点对对应的参数值都非零。

我们将说明如何利用回归方法来逐行逐列地求解 $\mathbf{\Theta}$ 和它的逆矩阵 $\mathbf{W}=\mathbf{\Theta}^{-1}$。简单起见，我们先看最后一行和最后一列。则式 17.13 中的右上分块可写为：

与式 17.7 中类似，这里我们将矩阵划分成了两个部分：第一部分是前 $p-1$ 个行和列，第二部分是第 $p$ 行和列。将 $\mathbf{W}$ 和它的逆矩阵 $\mathbf{\Theta}$ 都按这种方式划分，可得：

可解出：

其中的 $\beta=-\theta_{12}/\theta_{22}$，与式 17.9 中一样。然后将式 17.17 带入到 17.14 中，可得：

上式可被看成是 $X_p$ 对其他自变量带约束回归的 $p-1$ 个估计方程式，只不过将其中的交叉积矩阵的观测样本均值 $\mathbf{S}_{11}$ 替换成了从模型中得出的当前协方差矩阵估计 $\mathbf{W}_{11}$。

然后可以用简单的子集回归求解式 17.18。假设 $\gamma_{12}$ 中有 $p-q$ 个非零元素，即有 $p-q$ 个边被约束为零。这 $p-q$ 个行没有携带信息，可以被移除。我们可进一步地从 $\beta$ 中移除对应的 $p-q$ 个零值元素，将其缩短为 $\beta^*$，这样就得到了简化后的 $q\times q$ 方程组：

它的解为 $\hat{\beta}^*={\mathbf{W}_{11}^*}^{-1}s_{12}^*$。再填充回 $p-q$ 个零后，就得到了 $\hat{\beta}$。

尽管从式 17.16 中可见我们只是推导出了元素 $\theta_{12}$ 与 $\theta_{22}$ 的比例，但通过分区逆矩阵公式可容易推导出：

并且由于式 17.13 中的 $\mathbf{\Gamma}$ 的对角线元素为零，所以 $w_{22}=s_{22}$。

这样就得出了算法 17.1 概括出的简单迭代计算过程，来计算缺失边约束下的 $\hat{\mathbf{W}}$ 和逆矩阵 $\hat{\mathbf{\Theta}}$ 的估计。

算法 17.1：估计已知图结构的无向高斯图模型的修改回归算法

1. 初始化 $\mathbf{W}=\mathbf{S}$。
2. 对 $j=1,2,\dots,p,1,\dots$ 重复下列步骤，直到结果收敛。
   1. 将矩阵 $\mathbf{W}$ 划分成两部分：第一部分是除了 $j$ 行和列的所有其他元素，第二部分是第 $j$ 行和列。
   2. 利用式 17.19 中的简化方程组，从 $\mathbf{W}_{11}^*\beta^*-s_{12}^*=0$ 解出无约束的边参数 $\beta^*$。在 $\hat{\beta}^*$ 相应的位置中填充零元素，获得了 $\hat{\beta}$。
   3. 更新 $w_{12}=\mathbf{W}_{11}\hat{\beta}$。
3. 在（每个 $j$ 的）最后一次迭代中，求解 $\hat{\theta}_{12}=-\hat{\beta}\cdot\hat{\theta}_{22}$，其中的 $1/\hat{\theta}_{22}=s_{22}-w_{12}^T\hat{\beta}$。

值得注意的是这个算法从概念上是合理的。图的估计问题并不是 $p$ 个独立的回归问题，而是 $p$ 个耦合的问题。在步骤 2.2 中使用的不是观测样本的交叉积矩阵，而是当前步骤共同的 $\mathbf{W}$ 矩阵，这样就以一个恰当的方式将所有问题耦合在了一起。出乎我们的意料，我们在文献中没有找到这种方法。不过它与 Dempster (1972) 中的协方差选择过程（covariance selection procedure）有关联，并且与 Chaudhuri et al. (2007) 提出的迭代条件拟合过程（iterative conditional fitting procedure，ICFP）有相似之处。

这里给出一个来自 Whittaker (1990) 的例子。假设模型如图 17.4 所描绘，以及图中给出了它的经验协方差矩阵 $\mathbf{S}$。在这个问题上应用算法 17.1；举个例子，在步骤 2.2 中对第一个变量的修改回归中，变量三需要被排除（与变量一之间没有边）。迭代过程很快地收敛到下式的解：

需要注意 $\hat{\mathbf{\Sigma}}^{-1}$ 中的零元素，对应着缺失的边 $(1,3)$ 和 $(2,4)$。另外也注意到 $\hat{\mathbf{\Sigma}}$ 中这两个边对应的元素就是它与 $\mathbf{S}$ 全部的不同之处。矩阵估计 $\hat{\mathbf{\Sigma}}$ 是有时被称为 $\mathbf{S}$ 的正定还原（positive definite completion）的一个例子。

17.3.2 图结构的估计

在大多数的场景中我们并不知道图模型中哪些边可以被去掉，所以我们希望可以从数据本身发掘这个信息。近年一些研究者提出了使用 $L_1$（套索）正则化来达到这个目的。

Meinshausen and Bühlmann (2006) 提出了一个简单的方法：不再完整地估计 $\mathbf{\Sigma}$ 或 $\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$，而是只估计哪些 $\theta_{ij}$ 元素是非零的。他们将每个变量作为输出变量、其他变量作为自变量来拟合一个套索回归。如果变量 $i$ 对变量 $j$ 的估计系数非零，或者变量 $j$ 对变量 $i$ 的估计系数非零，那么 $\theta_{ij}$ 就被估计为非零（或者说他们使用的是“并集”规则）。他们证明了这个计算过程在渐进条件下对 $\mathbf{\Theta}$ 中非零元素的估计是一致的。

继续上一节中的推导，我们可以更系统化地在其中加入套索惩罚项。对带惩罚项的对数几率函数求最大化：

其中的 $\|\mathbf{\Theta}\|_1$ 是 $L_1$ 范数，即矩阵 $\mathbf{\Sigma}^{-1}$ 中所有元素绝对值的和。上式中省略掉了常数项。这个带惩罚项的似然函数的负数是 $\mathbf{\Theta}$ 的凸函数。

结果是通过套索回归可以得到这个带惩罚对数几率函数的精确最大值解。具体来说，只需要简单地将算法 17.1 步骤 2.2 中的修改回归变更为一个修改套索回归。以下介绍具体细节。

梯度方程组 17.13 的类比方程组：

这里我们使用了 次梯度（sub-gradient） 的概念。当 $\theta_{jk}\neq0$，则 $\operatorname{Sign}(\theta_{jk})=\operatorname{sign}(\theta_{jk})$；当 $\theta_{jk}=0$，则 $\operatorname{Sign}(\theta_{jk})\in[-1,1]$。继续与上述类似的推导，则可得到式 17.18 的类比：

（注意 $\beta$ 和 $\theta_{12}$ 的符号相反。）我们可说明这个方程组完全等价于一个套索回归的估计方程。

考虑一个常规的套索回归问题，输出变量为 $\mathbf{y}$、输入变量矩阵为 $\mathbf{Z}$。则套索回归的最小化准则函数是：

（见 68 页的式 3.52；这里为了方便起见添加了一个 $\frac{1}{2}$ 因子。）这个表达式的梯度为：

所以忽略一个 $1/N$ 因子后，$\mathbf{Z}^T\mathbf{y}$ 是 $s_{12}$ 的类比；而且我们用当前模型的交叉积矩阵估计 $\mathbf{W}_{11}$ 代替了 $\mathbf{Z}^T\mathbf{Z}$。

这就是 Friedman et al. (2008b) 基于 Banerjee et al. (2008) 之上提出的 图套索（graphical lasso） 方法。算法 17.2 概括了这个计算过程。

算法 17.2：图套索（graphical lasso）

1. 初始化 $\mathbf{W}=\mathbf{S}+\lambda\mathbf{I}$。在后续的步骤中，$\mathbf{W}$ 的对角线元素保持不变。
2. 对 $j=1,2,\dots,p,1,\dots$ 重复下列步骤，直到结果收敛。
   1. 将矩阵 $\mathbf{W}$ 划分成两部分：第一部分是除了 $j$ 行和列的所有其他元素，第二部分是第 $j$ 行和列。
   2. 利用修改套索回归的循环坐标下降（cyclical coordinate descent）算法（式 17.26），求解估计方程组 $\mathbf{W}_{11}\beta-s_{12}+\lambda\cdot\operatorname{Sign}(\beta)=0$
   3. 更新 $w_{12}=\mathbf{W}_{11}\hat{\beta}$。
3. 在（每个 $j$ 的）最后一次迭代中，求解 $\hat{\theta}_{12}=-\hat{\beta}\cdot\hat{\theta}_{22}$，其中的 $1/\hat{\theta}_{22}=s_{22}-w_{12}^T\hat{\beta}$。

Friedman et al. (2008b) 在每一个阶段求解修改套索问题时，使用了逐路径坐标下降方法（第 3.8.6 节）。以下介绍以下图套索算法中逐路径梯度下降的细节。令 $\mathbf{V}=\mathbf{W}_{11}$，则更新的形式为：

在 $j=1,2,\dots,p-1,1,2,\dots,p-1,\dots,$ 的每次迭代中，其中的 $S$ 是软阈值算子：

对所有自变量循环这个过程，直到结果收敛。

易证明矩阵 $\mathbf{W}$ 的解中的对角线元素 $w_{jj}$ 就等于 $s_{jj}+\lambda$，而且这在算法 17.2 的第一步中就已经确定了²。

图套索算法的计算非常快，可以在一分钟之内求解出 1000 个节点、稀疏性适中的模型。这个算法可以很容易地被修改成对边进行惩罚的参数 $\lambda_{jk}$ 的形式；由于 $\lambda_{jk}$ 将会意味着 $\hat{\theta_{jk}}$ 变为零，所以这个算法也包含了算法 17.1。将对稀疏协方差求逆矩阵的问题转换成一系列的回归问题，我们可以快速地求解并分析作为惩罚参数 $\lambda$ 的函数的解路径。更多细节可参考 Friedman et al. (2008b)。

图 17.1 展示了在流式细胞计量（flow-cytometry）数据集上应用图套索的结果。其中套索惩罚参数 $\lambda$ 设置为 14。在实践中，随着 $\lambda$ 的变化可以得到不同的图模型，从这些结果中可能看出一些信息。图 17.5 展示了四个不同的解。随着惩罚参数的增加，图模型变得更稀疏了。

最后需要注意在一个图模型中可能无法观测到某些节点的值；它们可能是缺失或是隐藏（hidden）变量。如果只是在一个节点上存在一些缺失值，那么可以使用 EM 算法来插补（impute）缺失值（练习 17.9）。但是有时整个节点会是隐藏（hidden）或潜伏（latent）的。在高斯分布模型中，如果一个节点所有数据都缺失，那么根据线性性质，可以对缺失的节点进行简单的平均，从而在观测到的节点上生成另一个高斯分布模型。所以引入隐藏节点不能为观测节点所产生的模型添加更多信息；实际上这样反而对其协方差矩阵添加了额外的结构约束。然而在（下一节）离散模型中，其内在的非线性使得隐藏节点称为扩展模型的一个强大的工具。

本节练习

练习 17.3

令 $\mathbf{\Sigma}$ 为 $p$ 个随机变量 $X$ 的协方差矩阵。将 $X$ 分成两个子集：$X_a=(X_1,X_2)$ 包含了前两个随机变量，$X_b$ 则为其他的随机变量。两个子集之间的偏协方差矩阵为 $\mathbf{\Sigma}_{a.b}=\mathbf{\Sigma}_{aa}-\mathbf{\Sigma}_{ab}\mathbf{\Sigma}_{bb}^{-1}\mathbf{\Sigma}_{ba}$。这等于 $X_a$ 中两个随机变量在对其他变量进行了线性调整后的协方差矩阵。如果假设高斯分布，那么这等于条件概率分布 $X_a|X_b$ 的协方差矩阵。 The partial correlation coefficient $\rho_{jk|\text{rest}}$ between the pair $X_a$ conditional on the rest $X_b$, is simply computed from this partial covariance. 记 $\mathbf{\Theta}=\mathbf{\Sigma}^{-1}$。

1. 证明 $\mathbf{\Sigma}_{a.b}=\mathbf{\Theta}_{aa}^{-1}$。
2. 证明如果 $\mathbf{\Theta}$ 的任意非对角线元素都为零，那么 the partial correlation coefficient between the corresponding variables is zero.
3. Show that if we treat Θ as if it were a covariance matrix, and compute the corresponding “correlation” matrix $$\mathbf{R} = \operatorname{diag}(\mathbf{\Theta})^{-1/2} \cdot \mathbf{\Theta} \cdot \operatorname{diag}(\mathbf{\Theta})^{-1/2} \tag{17.40}$$ 则 $r_{jk}=-\rho_{jk|\text{rest}}$。

练习 17.5

Consider the setup in Section 17.3.1 with no missing edges. Show that

are the estimating equations for the multiple regression coefficients of the last variable on the rest.

练习 17.6

Recovery of $\hat{\mathbf{\Theta}}=\hat{\mathbf{\Sigma}}^{-1}$ from Algorithm 17.1. Use expression (17.16) to derive the standard partitioned inverse expressions

Since $\hat{\beta}=\mathbf{W}_{11}^{-1}w_{12}$ , show that $\hat{\theta}_{22}=1/(w_{22}-w_{12}^T\hat{\beta})$ and $\hat{\theta}_{12}=-\hat{\beta}\hat{\theta}_{22}$ Thus $\hat{\theta}_{12}$ is a simply rescaling of $\hat{\beta}$ by $-\hat{\theta}_{22}$.

练习 17.7

Write a program to implement the modified regression procedure in Algorithm 17.1 for fitting the Gaussian graphical model with pre-specified edges missing. Test it on the flow cytometry data from the book website, using the graph of Figure 17.1.

练习 17.8

1. Write a program to fit the lasso using the coordinate descent procedure (17.26). Compare its results to those from the lars program or some other convex optimizer, to check that it is working correctly.
2. Using the program from (a), write code to implement the graphical lasso (Algorithm 17.2). Apply it to the flow cytometry data from the book website. Vary the regularization parameter and examine the resulting networks.

练习 17.9

Suppose that we have a Gaussian graphical model in which some or all of the data at some vertices are missing.

1. Consider the EM algorithm for a dataset of N i.i.d. multivariate ob- servations xi ∈ IRp with mean µ and covariance matrix Σ. For each sample i, let oi and mi index the predictors that are observed and missing, respectively. Show that in the E step, the observations are imputed from the current estimates of µ and Σ: $$\begin{align} \hat{x}_{i,m_i} &= \operatorname{E}(x_{i,m} | x_{i,o_i}, \theta) \\ &= \hat{\mu}_{m_i} + \hat{\Sigma}_{m_i,\sigma_i}\hat{\Sigma}_{o_i,o_i}^{-1} (x_{i,o_i} - \hat{\mu}_{o_i}) \tag{17.43} \end{align}$$ while in the M step, µ and Σ are re-estimated from the empirical mean and (modified) covariance of the imputed data: $$\begin{align} \hat{\mu}_j &= \sum_{i=1}^N \hat{x}_{ij} / N \\ \hat{\Sigma}_{jj'} &= \sum_{i=1}^N [(\hat{x}_{ij} - \hat{\mu}_j)(\hat{x}_{ij'} - \hat{\mu}_{j'}) + c_{i,jj'}] / N \tag{17.44}\end{align}$$ where ci,jj ′ = Σ̂jj ′ if j, j ′ ∈ mi and zero otherwise. Explain the reason for the correction term ci,jj ′ (Little and Rubin, 2002).
2. Implement the EM algorithm for the Gaussian graphical model using the modified regression procedure from Exercise 17.7 for the M-step.
3. For the flow cytometry data on the book website, set the data for the last protein Jnk in the first 1000 observations to missing, fit the model of Figure 17.1, and compare the predicted values to the actual values for Jnk. Compare the results to those obtained from a regression of Jnk on the other vertices with edges to Jnk in Figure 17.1, using only the non-missing data.