正则化模型在复杂网络结构中的研究

Research on Regularization Model in Complex Network Structure

摘要: 研究基于高斯图模型下的复杂网络结构，考虑惩罚项为The seamless-L0惩罚，提出了新的正则化模型。利用Minorize-Maximization算法，构造所提正则化模型的Minorization函数，再对该函数利用Maxdet算法求解模型。该模型同时进行模型选择和参数估计，并且估计结果具有稀疏性。数值模拟研究和基因调控网络的实例说明该模型有良好的参数估计和模型选择效果。

Abstract: Based on the complex network structure under the Guassian model, the penalty term is considered as the seamless-L0 penalty, and a new regularization model is proposed. The Minorize-Maximization algorithm was used to construct the Minorization function of the proposed regularization model, and then the Maxdet algorithm was used to solve the model. The model simultaneously performs model selection and parameter estimation. And the result is sparse. Numerical simulation studies and examples of gene regulatory networks show that the model has good parameter estimation and model selection effects.

1. 引言

复杂网络是由数据庞大的节点及节点之间错综复杂的关系所构成的网络结构。复杂网络普遍存在于现实生活中，如社会关系网络 [1] ，交通运输网络，基因网络等。因此，对复杂网络的研究成为近年来的热点问题 [2] 。

本文基于高斯图 [3] 模型(Gaussian Graphical Model, GGM)对复杂网络结构进行研究。其数学模型如下，假设 $X=\left(X^{\left(1\right)},\cdots ,X^{\left(p\right)}\right)$ 服从p维正态分布 $N_p\left(\mu ,\Sigma \right)$ ， $\mu$ 为均值， $\Sigma$ 为协方差矩阵。利用 $X$ 的样本 $X_1,X_2,\cdots ,X_n$ 估计 $A={\Sigma }^{-1}$ [3] 。当 $a_{ij}=0$ 时，表示 $X^{\left(i\right)}$ 与 $X^{\left(j\right)}$ 条件独立。高斯图模型中，随机向量用无向图表示 $G=VE$ [4] 。其中，V代表节点，每个节点分别对应随机向量的一个分量。 $E$ 代表边，两个节点之间无边，当且仅当这两个节点对应的随机变量条件独立，即 $a_{ij}=0$ 。因此，判断协方差矩阵的逆 $A$ 中的每个元素是否为0即是判断网络的结构。

在高斯图模型中，标准的模型选择方法是Greedy stepwise。该方法下参数估计结果是基于被选择的模型 [5] 。在每一步中，边的选择或删除通常通过在某个水平 $\alpha$ 的假设检验来完成。这种方法的缺点是没有正确地考虑到多重比较。另一个缺点是该方法在计算上具有复杂性 [5] 。为了解决以上问题，Drton和Perlman [6] 提出产生保守同步的 $1-\alpha$ 置信区间，并在每一步中利用这些置信区间进行模型选择。该方法的提出是基于对渐近性的考虑。Meinshausen和Buhlman [7] 为协方差选择提出具有较强计算能力的方法，可用于更庞大的高斯图。该方法为图中的每个节点进行领域选择，并结合选择结果来研究高斯浓度图模型的结构。Meinshausen和Buhlman说明了他们的方法对于稀疏高维图模型具有一致性。上述的这些方法，在模型选择和参数估计上都是分开的。并且，参数的估计是基于被选择的模型。Breiman [8] 也证明了这些方法在估计上具有不稳定性，即数据发生小的变动时，会导致估计结果有很大的不同。

本文提出的正则化模型在高斯图中可同时进行模型选择和参数估计，并使得估计结果具有稀疏性。正则子为The seamless-L₀ (SEL0) [9] 。该正则化模型对所要估计矩阵非对角线上的元素进行SEL0惩罚。该思想类似于SEL0惩罚在线性回归模型中的应用。本文使用MM算法和Maxdet算法求解该模型。最后，通过数值模拟和实例分析说明本文提出的正则化模型有良好的模型选择和参数估计效果。

2. 正则化模型

设 $X=\left(X^{\left(1\right)},\cdots ,X^{\left(p\right)}\right)$ 服从p维正态分布 $N_p\left(\mu ,\Sigma \right)$ ， $\mu$ 为均值， $$ 为协方差矩阵。设 $A={\left(a_{ij}\right)}_{1\le i,j\le p}={\Sigma }^{-1}$ 为协方差矩阵的逆。基于随机变量 $X$ 样本观测值 $X_1,X_2,\cdots ,X_n$ 的对数似然函数是

$\frac{n}{2}\log \left|A\right|-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\mu \right)}^{\prime }A\left(X_i-\mu \right)}$ (1)

当不依赖于常数A和 $\mu$ 时，设 $\left(\mu ,\Sigma \right)$ 的最大似然估计为 $\left(\bar{X},\bar{C}\right)$ 。其中

$C=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\bar{X}\right)}^{\prime }\left(X_i-\bar{X}\right)}$

通常利用样本协方差矩阵 $S=n\bar{C}/\left(n-1\right)$ 进行协方差矩阵的估计。因此， $A={\Sigma }^{-1}$ 可利用 $S^{-1}$ 或者 ${\bar{C}}^{-1}$ 进行估计。但是，当变量的维数很大时，矩阵 $S^{-1}$ 稳定性差。并且，当 $n\ge p$ 时，矩阵 $S^{-1}$ 是正定矩阵，矩阵中的元素不包含零项。所以，以此方法进行的估计结果是不具有稀疏性的图模型。

因此，通过引入L₁正则子进行模型选择和参数估计。本文选择的L₁正则子为SEL0惩罚，改惩罚形式是对L₀惩罚的形式的模仿 [10] 。

但当变量维数很大时，该方法在计算上具有复杂性，并且具有不连续性和复杂性 [9] 。为解决上述问题，Dicker L.和Huang B.等提出SEL0惩罚。该惩罚形式为

$p_{\text{SEL0}}\left({\beta }_j\right)=\frac{\lambda }{\log \left(2\right)}\left(\frac{\left|{\beta }_j\right|}{\left|{\beta }_j\right|+\tau }+1\right)$ (2)

其中， $\tau >0$ 是调节参数，当 $\tau$ 很小时， $p_{\text{SEL0}}\left({\beta }_j\right)\approx \lambda I\left\{{\beta }_j\ne 0\right\}$ ，并且，当 $\tau =0.01$ 时，该方法呈现更好的效果。

因此，为了得到具备稀疏性的图模型结构，并更加准确的估计 $A={\Sigma }^{-1}$ 。本文结合The seamless-L₀惩罚(SEL0)和 $-\log \left|A\right|+\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\mu \right)}^{\prime }A\left(X_i-\mu \right)}$ 寻找 $\left(\hat{\mu },\hat{A}\right)$ 的最大值。即

$\log \left|A\right|-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_i-\mu \right)}^{\prime }A\left(X_i-\mu \right)}$ st ${\displaystyle \underset{i\ne j}{\sum }\frac{1}{\log \left(2\right)}}\log \left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)\le \lambda$ (3)

其中， $\lambda \ge 0$ 是调节参数。

因为对观测值作中心化处理。所以 $\hat{\mu }=0$ 。因此，(1)式可写为

$\log \left|A\right|-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{X^{\prime }}_{i}A{X}_i}$ st ${\displaystyle \underset{i\ne j}{\sum }\frac{1}{\log \left(2\right)}}\log \left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)\le \lambda$ (4)

进一步地，我们可以将(2)式写为

$\log \left|A\right|-tr\left(A\bar{C}\right)$ st ${\displaystyle \underset{i\ne j}{\sum }\frac{1}{\log \left(2\right)}}\log \left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)\le \lambda$ (5)

由于(3)式中的目标函数和可行区域都是凸的，因此，可以利用拉赫朗日形式将(3)式等价得写为

$L\left(A\right)=\log \left|A\right|-tr\left(A\bar{C}\right)-\lambda {\displaystyle \underset{i\ne j}{\sum }\frac{1}{\log \left(2\right)}}\log \left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)$ (6)

因此，最大化(4)式中的 $L\left(A\right)$ ，即可得到图模型的协方差矩阵的逆的估计。

3. Maxdet算法与参数选择

本文利用Maxdet算法 [11] 求解目标函数(6)。但由于本文中目标函数 $L\left(A\right)$ 的SEL0惩罚项是 $\log$ 和 $L_1$ 惩罚函数的复合形式，因此，目标函数较难直接进行优化，需要利用Minorize-Maximization (MM) [12] 算法，构造 $L\left(A\right)$ 的Minorization函数 $H\left(a|a^n\right)$ ，其中， $a^n$ 是MM算法的第n部估计，也是参数a的一个定值。该函数逼近于原目标函数 $L\left(A\right)$ 且易于优化。因此，本节先构造出函数 $H\left(a|a^n\right)$ ，再通过对函数 $H\left(a|a^n\right)$ 利用Maxdet算法，得到逼近于 $L\left(A\right)$ 的估计值。最后使用BIC方法进行参数估计。

首先构造 $L\left(A\right)$ 的Minorization函数 $H\left(a|a^n\right)$ 。

若满足以下条件，则目标函数 $L\left(A\right)$ 可用构造函数 $H\left(a|a^n\right)$ 代替。

$c_1.L\left(A\right)\ge H\left(a|a^n\right)$

$c_2.L\left(A^n\right)=H\left(a|a^n\right)$

由于当 $x>0$ ， $-\log \left(1+x\right)\ge -x$ 故

$\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\log \left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)}-\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\log \left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)}\ge \frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\left(\frac{\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1}{\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1}\right)}$

不等式两边同时加上

$\log \left|A\right|-tr\left(A\bar{C}\right)-\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\log \left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)}$

得

$L\left(A\right)\ge \log \left|A\right|-tr\left(A\bar{C}\right)-\frac{\lambda }{\log \left(2\right)}\left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)-\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\left(\frac{\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1}{\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1}\right)}$

令不等式右边为 $H\left(a|a^n\right)$ ，即

$H\left(a|a^n\right)=\log \left|A\right|-tr\left(A\bar{C}\right)-\frac{\lambda }{\log \left(2\right)}\left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)-\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\left(\frac{\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1}{\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1}\right)}$ (7)

故 $L\left(A\right)\ge H\left(a|a^n\right)$ ，条件 $c_1$ 满足。

又因

$H\left(a^n|a^n\right)=\log \left|A\right|-tr\left(A\bar{C}\right)-\frac{\lambda }{\log \left(2\right)}\left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)-\frac{\lambda }{\log \left(2\right)}{\displaystyle \underset{i\ne j}{\sum }\left(\frac{\frac{\left|a_{ij}{}^n\right|}{\left|a_{ij}{}^n\right|+\gamma }+1}{\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1}\right)}$

即

$H\left(a^n|a^n\right)=\log \left|A\right|-tr\left(A\bar{C}\right)-\frac{\lambda }{\log \left(2\right)}\left(\frac{\left|a_{ij}^n\right|}{\left|a_{ij}^n\right|+\gamma }+1\right)$

故 $H\left(a^n|a^n\right)=L\left(A\right)$ ，条件 $c_2$ 满足。

因此， $L\left(A\right)\ge H\left(a|a^n\right)$ ，且 $H\left(a^n|a^n\right)=L\left(A\right)$ 故满足 $c_1$ ， $c_2$ 两个条件。所以 $L\left(A\right)$ 的Minorization函数为 $H\left(a|a^n\right)$ 。

下面利用Maxdet算法求解式(7)。Maxdet问题的形式如下：

${\min }_{x\in R^m}{b}^{\prime }x-\log \left|G\left(x\right)\right|$

故对(7)式取负号，将其转化为求最小值问题。假设 $a_{ij}$ 的符号是已知的，设为 $s_{ij}$ 。根据Maxdet问题的形式，可将转化后的(5)式写成如下形式：

$\min 2{\displaystyle \underset{i<j}{\sum }{c}_{ij}{a}_{ij}}+{\displaystyle \underset{i}{\sum }{c}_{ii}{a}_{ii}}-\log \left|{\displaystyle \underset{i}{\sum }{a}_{ii}{I}^{\left(i\right)}}+{\displaystyle \underset{i<j}{\sum }{c}_{ij}{I}^{\left(ij\right)}}\right|st$

${\lambda }^{\prime }-2{\displaystyle \underset{i<j}{\sum }\left(\frac{\left|a_{ij}\right|}{\left|a_{ij}\right|+\gamma }+1\right)s_{ij}}\ge 0，s_{ij}{a}_{ij}\ge 0$ (8)

其中， $I^{\left(ij\right)}$ 为第 $\left(i,j\right)$ 项为1，其余项为0的 $n\times n$ 的矩阵， $I^{\left(i\right)}$ 为第 $\left(i,i\right)$ 项为1，其余项为0的 $n\times n$ 的矩阵。

关于参数 $\lambda$ 的选择，常用的方法是多重交叉验证法，但计算上更有效的替代方案是用于模型选择和估计的BIC准则。因此，本文选择如下BIC准则选择参数 $\lambda$ [13] ：

$\text{BIC}\left(\lambda \right)=-\log \left|A\left(\lambda \right)\right|+tr\left\{A\left(\lambda \right)\bar{C}\right\}+\frac{\log n}{n}{\displaystyle \underset{i\ne j}{\sum }{\hat{e}}_{ij}\left(\lambda \right)}$ 其中，如果 $a_{ij}=0$ ，则 $e_{ij}=0$ ；其余情 ${\hat{e}}_{ij}=1$ 。

下面式算法的具体步骤：

第一步：令 ${\hat{A}}_{old}={\bar{C}}^{-1}$ 并且对任意的 $i\ne j$ ， ${\hat{A}}_{old}={\bar{C}}^{-1}i\ne j{s}_{ij}=sign\left\{{\left({\hat{A}}_{old}\right)}_{ij}\right\}$ 。

第二步：令 ${\hat{A}}_{new}$ 在正定矩阵集上解决(6)式。

第三步：如果 ${\hat{A}}_{new}={\hat{A}}_{old}$ ，停止迭代并令 $\hat{A}={\hat{A}}_{new}$ 。否则，对所有的 $\left(i,j\right)$ 设 ${\hat{A}}_{new}={\hat{A}}_{old}$ ，并且 $s_{ij}=-s_{ij}$ ，所以 ${\hat{c}}_{ij}=0$ 。然后返回第二步。

第四步：使用BIC准则选择 $\lambda$ ，并输出相应的A的估计值。

4. 数值模拟

本文中的网络图模型由R语言huge包中的huge.generator函数生成，所生成的图形模式为随机模型。并且在此函数的输出结果中可得到该图形模式下随机生成的样本数据及其协方差矩阵。本文设置样本量为n = 200，维数p分别为15，30，50，100。分别比较本文提出的方法和Meinshausen & Buhlmann (mb)提出的方法在不同维数下参数估计和模型选择的情况。用KL损失量以及FPR (false positive rate)，TPR (true positive rate)作为评价正则化模型在参数估计和模型选择上优良性的标准。其中，KL损失量(KL loss)为 [14] ：

$\text{KL}=-\log \left|\hat{A}\right|+tr\left(\hat{A}\hat{\Sigma }\right)-\left(-\log \left|{\Sigma }^{-1}\right|+p\right)$

FPR和TPR为：

$\text{FPR}=\text{FP}/\left(\text{FP}+\text{TN}\right)$ ； $\text{TPR}=\text{TP}/\left(\text{TP}+\text{FN}\right)$ 。其中，TPR越大，FPR越小说明模拟效果越好。

模拟结果如下(表1~表4)。

Table 1. p = 15

表1. p = 15

Table 2. p = 30

表2. p = 30

Table 3. p = 50

表3. p = 50

Table 4. p = 100

表4. p = 100

模型选择的目标是在高维数据中挑选出更具代表性的变量。表数据结果表明，在不同维数下，SEL0方法相比于mb方法在模型选择和参数估计方面均是最优的。并且在p = 15时，SEL0方法下所得到的KL loss最小，FPR达到最小，TPR达到最大。随着维数的增加，虽然KL loss值和FPR值逐渐增大，TPR值逐渐减小，但总体来说，gSEL0方法仍有较好的模型选择和参数估计能力。故综合以上结果，本文提出的gSEL0方法有良好的模型选择和参数估计能力。

5. 实例分析

本文将提出的正则化模型应用到基因调控网络的数据中。数据来源于Faith等人的论文数据 [15] 。该数据包含两个对象。第一个对象是40行，153列的矩阵，该矩阵的行和列代表了Escherichaia coil (E. coil)细菌在40种不同实验下，153种基因的基因表达水平，基因表达的测量结果指示了153种基因在每个特定的实验干扰下的活动水平。该153种基因是转录因子，在基因调控关系中有重要作用。第二个对象是一个邻接矩阵，概括了生物学文献中发现的E. coil细菌的153种基因之间已经实际存在的调控关系。结果如图1所示。

下面通过统计软件，首先根据邻接矩阵统计出生物学文献中已经发现的153种基因之间实际存在的调控关系的个数；然后利用本文提出的正则化模型，统计出在该模型下得到的生物学文献中已发现的实际调控关系及生物学文献中未被发现的调控关系的总量。步骤和结果如下：

1) 对邻接矩阵中基因之间的实际调控关系的个数进行统计，所得结果如下：

## IGRAPH UN-153 209

## attr: name (v/c)

结果表明：生物学文献中已经发现的基因之间的实际调控关系共209对。

2) 利用第一个数据对象中40行，153列的矩阵，根据本文提出的正则化模型，统计出在该模型下得到的生物学文献中已发现的实际调控关系及生物学文献中未被发现的调控关系的总量。所得结果如下：

Figure 1. Actual regulatory relationships among 153 genes

图1. 153个基因之间实际存在的调控关系

IGRAPH 68d9b0b U-153 554--

结果表明：根据本文提出的正则化模型共找到554对调控关系。

3) 最后统计在步骤2中利用本文提出的正则化模型得到的生物学文献中已发现的实际调控关系的个数。所得结果如下：

## IGRAPH UN-153 192--

结果表明：通过本文提出的正则化模型在已发现的209对实际调控关系中共得到了192对调控关系。这一结果进一步说明了该正则化模型在实际应用中的优良性。

6. 结论

本文基于高斯图模型进行复杂网络结构的研究，根据SEL0惩罚在线性回归模型中应用的思想，将其应用到图模型的模型选择和参数估计中，提出了新的正则化模型。本文使用Minorize-Maximization (MM)算法和Madex算法求解该模型。数值模拟实验结果显示，本文提出的正则化模型相比于Meinshausen & Buhlmann (mb)方法，KL loss值和FPR值更小，TPR值更大。说明本文提出的正则化模型在参数估计和模型选择方面都有良好的效果。最后，将该模型应用到E. coil细菌基因调控关系的实际数据当中，结果表明，本文提出的正则化模型对基因之间实际存在的调控关系有较好的估计效果，说明该方法在应用中具有良好的适应性。