基于条件独立性的LiNGAM模型剪枝算法

郝志峰　吕宏伟　蔡瑞初　袁　畅

(广东工业大学计算机学院　广东 广州 510006)



基于条件独立性的LiNGAM模型剪枝算法

郝志峰吕宏伟蔡瑞初袁畅

(广东工业大学计算机学院广东 广州 510006)

如何根据观察数据来推断因果网络结构是统计学和机器学习领域的重要问题。近年来学者们取得了许多研究成果，LiNGAM算法是其中一种经典的线性因果推断算法。但LiNGAM算法采用的剪枝策略时间复杂度较高，且在稀疏图上准确率低。为此，提出一种基于条件独立性测试的剪枝算法来解决这个问题。该算法首先将变量根据因果顺序重新排列，再按照该次序采用偏相关系数检验变量之间的条件独立性。大量的实验结果表明，基于条件独立性的剪枝算法在稀疏图上比LiNGAM的剪枝算法获得更高的准确率与执行效率。

线性因果关系偏相关条件独立性剪枝

0　引　言

继贝叶斯网络结构学习之后，因果网络结构推断受到广泛的关注。在因果推断中，因果关系的可识别性是因果推断的基本前提。近些年的研究表明非高斯噪音与非线性关系对于因果关系识别有着重要作用[1,2]。Shimizu等人发现当变量之间为线性连续的关系并且噪音均服从非高斯分布时，变量之间的因果关系是可识别的，并提出线性连续非高斯的因果推断模型LiNGAM模型[1]。后续Hoyer等人[2]将LiNGAM模型扩展为加噪因果模型，并指出若变量之间均是非线性关系或者变量中存在的噪音项均为非高斯分布，那么变量之间的因果关系是可识别的。

与LiNGAM算法和DirectLiNGAM算法所采用的剪枝方法不同，本文提出基于贝叶斯网络结构的推断算法采用的条件独立性原理；将变量按照LiNGAM算法所得出的因果顺序重新排列，采用偏相关系数检验依次测试变量之间的条件独立性的剪枝算法。若将SGS算法[6]或PC算法[7]等贝叶斯结构推断算法直接用作剪枝算法时存在着误剪枝率高的问题。与基于条件独立性原理的SGS算法[6]和PC算法[7]等贝叶斯网络结构推断的经典方法不同，本文算法同时结合了因果顺序与条件独立性原理。大量的实验结果表明，本文算法准确率高，且更为稳定。

1　相关背景

1.1LiNGAM模型

LiNGAM是线性连续、非高斯、有向无环的因果模型。LiNGAM模型中的变量是按照因果顺序依次生成的。因此将变量按照因果顺序进行排列后，位于后面的变量不可能是前面变量的因变量。在实际应用中，观察到的变量的排列是随机的，与因果顺序是不同的。本文根据观察中样本的排列次序(本文记为观察顺序)将变量依次记为{x1,x2,…,xn}，将因果顺序记为k(i),i∈[1,n]。(i)∈[1,n]代表了观察顺序中第i个变量在因果顺序中的位置，那么变量的生成过程可描述为：

(1)

式中，ei为服从非高斯分布的噪音项，噪音项之间两两互相独立；bij不为0表示存在xj指向xi的边。LiNGAM模型用矩阵的形式表述为：

x=Bx+e

(2)

若采用置换矩阵P∈Rn×n来描述因果顺序k(·)，求得B′=PBP，那么根据因果顺序的意义可知，B′为严格的下三角矩阵且对角线的元素均为0。

1.2LiNGAM算法简述

LiNGAM算法是LiNGAM模型中结构推断的经典算法。总的来说，LiNGAM算法分为估计与剪枝两个阶段。

第一阶段是估计阶段，该阶段能够得出因果顺序以及初步估计出整个因果结构(即矩阵系数)。由式(2)可推出：Wx=e,W=(I-B),I为单位矩阵，W称为混合矩阵。LiNGAM算法首先利用ICA算法根据观察数据估计出分离矩阵W′，但是估计得到的W′矩阵存在两个问题：(1)W′中对角线可能为0，这与W′的对角线全为1的约束相矛盾；(2)变量是随机排列的，与因果顺序不符。为了解决这两个问题，LiNGAM算法对W′的行重新排序，使得W′对角线均不为0。然后，LiNGAM算法将W′每一行除以该行对角线的值(归一化)得到W″，并计算出B=I-W″。再根据因果顺序对应的系数矩阵是对角线均为0的下三角矩阵这个约束，LiNGAM算法同时对B的行和列进行重新排列得到B′，使得B′为下三角矩阵，从而得到因果顺序。最后LiNGAM算法根据最小二乘法得出具体矩阵系数。估计阶段所得出的关系矩阵B所对应的有向无环图是全相联的，但是实际中常见的图是稀疏的，因此对B剪枝是必要的步骤。

第二阶段是剪枝阶段，该阶段是对估计阶段得出的关系矩阵B中的每一条不为0的边进行检验。在文献[1]中，Shimizu等人提出了剪枝算法MODELFIT。该算法首先采用Wald检验来检验每条边的显著水平。对于每条未能通过Wald检验的边，MODELFIT均根据卡方检验等方式检验剪掉该边前后的模型差异的显著性。MODELFIT根据剪掉某条边前后模型差异显著水平来决定是否剪枝。该算法有两个缺点：(1) MODELFIT每剪去一条边，都需要检验模型差异的显著性，这项操作需要大量耗时的矩阵求逆和矩阵求根等运算(其复杂度皆为维度的立方)；(2)对比实验表明MODELFIT剪枝算法在稀疏图上的准确率相对较低。

1.3基于偏相关系数的条件独立性测试

本文采用偏相关系数检验作为条件独立性测试的方法。下面引入几个重要的定义。

定义1偏相关系数随机变量x、y在给定受控变量的集合Z时的偏相关系数是Rxz与Ryz之间的相关系数，其中Rxz为变量x与Z线性回归的残差，Ryz为变量y与Z线性回归的残差。

偏相关系数是指两个变量在通过线性回归的方式消除受控变量影响后的相关系数。若所涉及的随机变量均为高斯分布，随机变量x、y在给定受控变量的集合Z时的偏相关系数等于0，等价于变量x与变量y在条件集Z上条件独立；在其他情况下，偏相关系数是条件独立性的一种近似的计算方式。判断偏相关系数是否严格为0常采用假设检验的方式。Fisher Z 变换与T检验是常见的两种方式。

定义2条件独立性设随机变量集合为V={x1,x2,…,xn},X、Y、Z是V的任意非空子集。若有P(X|Y,Z)=P(X|Z),并且(Y,Z)>0，那么称X、Y在给定Z时条件独立。

定义3d-分离准则设X、Y、Z是贝叶斯网络G中任意三个互不相交的节点集合，称Z在图G中d-分离节点集X和Y。若对任意的从X的节点到Y的一个节点的路径P均被Z阻断，也就是路径P上存在一个结点xi满足下列其中一个条件：

(1)xi在P上形成碰撞，即→xi←，且xi及其后代结点都不在Z中；

(2)xi在P上不存在碰撞，即→xi→或←xi←，且xi在Z中。

定义4马尔科夫毯设贝叶斯网络G=，集合V代表节点，集合E代表边。设xi是G中任意节点，xi的马尓科夫毯是由该节点的所有的父亲节点、孩子节点以及配偶节点(xi每个孩子节点的其他父亲节点，称为xi的配偶节点)组成的集合。

贝叶斯网中，节点xi、xj在节点集合Z上条件独立等价于Z阻隔xi、xj之间的所有的路径，即Z集合d-分离了xi、xj。若贝叶斯网中的节点xi、xj之间不存在边，那么必定存在集合Zd-分离了xi、xj；反之，必不存在集合Zd-分离了xi、xj。因此通过逐一测试集合V/{xi,xj}的子集能否d-分离xi、xj，就可断定xi、xj之间关系。该方法称为条件独立性测试。本文中的条件独立性测试是基于偏相关系数，其判断Z集合是否d-分离了节点xi、xj时，根据样本计算节点xi、xj在给定集合Z时的偏相关系数。若偏相关系数显著为0，那么Z集合d-分离了节点xi、xj；反之，Z集合未能d-分离节点xi、xj。节点的马尔科夫毯d分离了该节点与任意不在其马尔科夫毯中的节点。本文将节点xi的马尔科夫毯记为PC(xi)。

2　基于条件独立性的LiNGAM剪枝算法

2.1剪枝问题描述

设LiNGAM模型中变量为{x1,x2,…,xn}，样本矩阵X∈Rn×m，n为样本中变量个数或称为样本维度，m为样本个数。LiNGAM算法的估计阶段的所得到的关系矩阵为B，因果顺序为k(·)。剪枝算法是根据上述已知条件，对B进行剪枝，返回剪枝后的关系矩阵。按因果顺序k(·)排列变量后，式(1)变为：

(3)

其中i、j是变量在因果顺序中的位置。剪枝问题与原结构推断问题的不同之处是因果顺序是已知条件，因此式(3)中的系数矩阵B是对角线均为0的下三角矩阵，并且是稀疏的。

2.2剪枝算法流程

在已知因果顺序基础上，本文提出了基于条件独立性测试的剪枝算法，流程如算法1中所示。

算法1基于条件独立性测试的剪枝算法

输入:样本矩阵X∈Rn×m,因果顺序k(·)与待剪枝系数矩阵B。

输出:剪枝后的矩阵B′(bij=0表示xi与xj之间不存在边)。

预处理:首先计算B′=PBP,X′=PX,P为因果顺序k(·)所对应的置换矩阵。再按照B′中变量的顺序，从变量x1开始剪枝，直到xn结束。记当前处理的变量为xi,i∈[1,n]。

步骤1：根据条件独立性测试，逐一判断当前处理变量xi与xj∈{x1,x2,…,xi-1}的关系。具体做法是根据样本计算每个变量xj,j∈[1,i)与变量xi在给定PC(xi)时的偏相关系数，并判断偏相关系数是否显著为0。若偏相关系数为0，则令bij=0。

步骤2：对集合PC(xi)中的变量进行内部的条件独立性测试。具体做法是通过计算每个变量xt∈PC(xi)与变量xi在给定Z=PC(xi)/xt时的偏相关系数，再次判断xi、xt的关系。若偏相关系数的P值高于显著水平(0.05)，则令bit=0。

步骤3：若已处理完所有变量那么剪枝完毕，返回剪枝后的B′，否则返回步骤一并继续处理下一个变量。

在剪枝之前，本文将B按照因果顺序k(·)重新排列为B′,采用逐步扩展的方式剪枝。

在步骤1中，变量{x1,x2,…,xi-1}之间的关系是已知的，仅有xj∈{x1,x2,…,xi-1}与当前处理变量xi的关系是需要进行判定的。在变量{x1,x2,…,xi}组成的有向无环图中，xi是最后一个变量，它的马尔科夫毯仅由其父变量组成，因此步骤1初步得出了PC(xi)。假设检验采用了Fisher Z变换，显著水平为0.05。

步骤2是进一步测试PC(xi)中的变量是否满足马尔科夫毯的定义。虽然步骤1中已经求出了PC(xi),但是实验表明总体准确率较低，这意味着步骤所求得的PC(xi)中仍存在不属于PC(xi)的变量未能被剪除。从另外一个角度，算法的步骤1是根据PC(xj)去判断xi、xj的关系，步骤2则是利用PC(xi)去判断xi、xj之间的关系。步骤2中采用了T检验来计算P值。

2.3正确性分析

根据条件独立性的定义，在n个变量组成的贝叶斯网络中判定两个变量是否存在边所需的测试次数最多为2n-2次。本文基于引理一将所需测试的变量集合减小为变量的马尔可夫毯，这大幅度降低了所需测试的次数。本文进一步结合剪枝问题中因果顺序的已知条件将判断两个变量关系所需的条件独立性测试次数减少为两次。实验表明有效地降低了误剪枝率。

引理1在贝叶斯网中，对于任意两个不存在边节点xi、xj，一定存在PC(xi)的子集(或PC(xj)的子集)Zd-分离了xi、xj。

根据引理1，判断LiNGAM模型中任意两个变量xi、xj的关系，需要逐一测试的集合是xi或xj的马尔科夫毯集合的所有子集。若存在集合Z是PC(xi)的子集(或PC(xj)的子集)使得xi、xj条件独立，那么xi、xj不存在边；若不存在这样的集合Z，那么xi、xj存在边。采用类似方法的还有BASSUM[9]与HITON[10]等。该方法有效减少了所需的条件独立性测试次数，同时准确率高，然而存在容易误剪枝的问题。这是因为变量的条件独立性一般是根据样本采用某种方法估计得出的，存在一定的误差。同时由于每一条保留的边均需要经过多次条件独立性测试，这使得原本存在的边有较高概率被错误剪掉。若将PC等贝叶斯算法用作剪枝算法也存在类似的误剪枝率高的问题。

由引理1得出的剪枝方法仍存在误剪枝率高的问题，核心问题是如何减少条件独立性测试次数。结合因果顺序，本文提出按照因果顺序依次剪枝，有效地减少了所需条件独立性测试次数。

引理2设变量{x1,x2,…,xn}符合LiNGAM模型，B为系数矩阵，因果顺序k(·)与观察顺序{1,2,…,n}是一致的,那么∀xt,t∈[1,n]，{x1,x2,…,xt}组成了一个完整的LiNGAM模型。

证明：对于∀xt,t∈[1,n]，根据因果顺序的定义可知∀xj,j∈[t+1,n]不会是∀xi,i∈[1,t]的因变量。这使得xt后面的变量xj不会对xt以及xt前面的变量xi产生影响，因此{x1,x2,…,xt}是一个完整的LiNGAM模型。

由引理2可知，判断变量xi、xj之间是否存在关系时，那些位于后面的变量xk,max{i,j}

上文按照因果顺序剪枝，缩减了剪枝问题的规模。这里进一步结合引理1与引理2，提出了仅需两次独立性测试判定两个变量之间关系的方法。实验表明，该方法能够有效剪枝，并且误剪枝率较低。

引理3 在变量{x1,x2,…,xn}组成的LiNGAM模型中，设因果顺序k(·)与观察顺序{1,2,…,n}是一致的，变量{x1,x2,…,xt}是已知的子模型。对于∀xi,i∈[1,t]，变量xi、xt+1存在边近似等价于集合PC(xi)未能d-分离xi、xt+1且集合PC(xt+1)/xi未能d-分离xi、xt+1。

证明：引理3中，判断xi、xt+1的关系时，首先根据已知条件能够判断PC(xi)是否d-分离xi,xt+1。根据马尔科夫毯的性质可知,PC(xi)未能d-分离xi、xt+1表明变量之间有较高的可能性存在边，由此可估计得出xt+1的所有父变量集合，记为Pa(xt+1)。在变量{x1,x2,…,xt,xt+1}组成的模型中，根据马尔科夫毯的定义可知PC(xt+1)=Pa(xt+1)。同时，若PC(xt+1)/xi未能d-分离xi、xt+1，进一步表明xi、xt+1之间可能存在边。反之，若xi、xt+1之间存在边，引理三中所描述的d-分离性质是必须满足的。

依据引理3，判断任意两个变量xi、xj的关系仅需两次测试。该方法优点是仅需两次条件独立性测试即可判断边是否存在；缺点是由于只是一种近似的判断，可能将少数不存在的边误判为存在。然而实验表明，该策略在准确率上并无明显劣势，相反降低了误剪枝率。

根据条件独立性测试解决剪枝问题的难点是如何降低误剪枝率。设当前处理的变量为xt+1，变量{x1,x2,…,xt}是结构已知的LiNGAM模型。根据引理3即可判断出xt+1与∀xi,i∈[1,t]之间的关系，得出变量{x1,x2,…,xt′,xt′+1}的完整结构，从而可以按照因果顺序继续处理下一个变量。本文选取变量的马尔科夫毯作为条件独立性测试的条件集合，提出仅需两次独立性测试来判断两个变量之间的因果关系的方法，较好地解决了误剪枝率高的问题。

3　实　验

实验中本文的剪枝算法与Resampling、OLSBOOT、MODELFIT和Adaptive Lasso四种剪枝方法在模拟数据上做了充分的比较[1]。综合考虑剪枝后关系矩阵的准确率与召回率(召回率高，表明剪枝算法误剪枝率低)，本文采用F1-score作为评价标准。本文算法在实验中命名为PaC与PaCOneStep。PaCOneStep是仅采用步骤1的剪枝算法。Resampling、MODELFIT、OLSBOOT与Adpative Lasso是文献[1,4]中的所采用的剪枝算法。OLSBOOT与Resampling算法均是重采样的方法，其中OLSBOOT是基于BOOTSTRAP的方法。

本文实验在MATLAB 2011b中完成，硬件配置为4 GB内存，酷睿双核2.0 GHz 。

3.1模拟数据的生成以及参数

本文实验中采用了两种方式生成模拟数据。第一种是LiNGAM算法中所提供的模拟数据的生成算法[1]。该算法在高斯分布的基础上采用了随机选取的方式使得干扰变量分布服从超高斯或亚高斯分布。模拟数据在生成后，变量的排列顺序是随机排列的，这使得观察顺序与因果顺序不一致。第二种生成模拟数据的方式是DirectLiNGAM算法[4]所采用的，该算法仅可生成稀疏结构的数据。

模拟数据的几个主要参数为变量的最大入度v(变量最多可以存在的父变量个数)，变量数n，样本数m。

3.2实验效果

首先，图1是在不同稀疏程度下各个剪枝算法对比的实验结果。变量数n=10，样本数目m=1000，变量最大入度v={1,2,3,4,5,FULL}；v=FULL表示采用了全相联关系。将正确的因果顺序作为了已知条件。 由图1中可看出本文算法仅在变量最大入度v=1时(此时的结构退化为直线片段)稍差于Resampling，在其他情况下均优于其他算法。PaCOneStep与Adaptive Lasso算法在变量最大入度v较小的情况下，准确率很低。MODELFIT也存在类似问题。在全相联的结构上各个剪枝算法的表现反应了其误剪枝的情况，由图1可看出本文算法误剪枝率最低。

图1　第一种数据产生方式下各个算法的F1-Score性能比较。

其次，图2-图5是各个剪枝算法作为LiNGAM算法的剪枝算法时在两种不同数据生成方式下对比的实验结果。图2、图3的数量数n=10。图4、图5的数量数n=100,下标1表示给出正确的因果顺序，下标2表示采用LiNGAM估计因果顺序。本文剪枝算法在两种数据生成方式中都得到较优的结果，仅在第二种数据生成方式下略差于Resampling。相比之下，其余剪枝算法存在着不稳定或不适用于变量较多的情况等问题。图2中，MODELFIT、OLSBOOT略差于本文算法；图3中，本文算法仅差于Resampling，OLSBOOT与本文算法的剪枝结果仍然是接近的。Resampling算法在采用第一种数据生成方式时，其剪枝结果是最差的，因此是不稳定的。在变量数目较少时，OLSBOOT与本文算法均能取得较好的剪枝结果，但是OLSBOOT不适于变量多的情况，图4、图5验证了这一点。变量数目为100且采用LiNGAM估计因果顺序时，Resampling与本文算法的剪枝结果是较为接近的，在已知正确的因果顺序后，本文算法具有较大的优势。

图2　采用第一种数据生成方式1

图3　采用第二种数据生成方式1

图4　采用第一种数据生成方式2

图5　采用第二种数据生成方式2

最后，表1是各剪枝算法的耗时情况的比较。本文剪枝算法时间复杂度小于Adaptive Lasso与MODELFIT，高于Resampling与OLSBOOT。本文算法的基本步骤是偏相关系数的计算，其时间复杂度是条件集合中变量个数的三次方。在稀疏图中，变量xi的马尔科夫毯中变量个数远小于n，记为mi,0≤mi≤n-2。本文算法时间复杂度的上界为O(n2×(max{m1,m2,…,mn})3)。

表1　算法所耗时间T的变化(单位：秒，-代表耗时过长)

4　结　语

不同于MODELFIT、Adaptive Lasso等剪枝算法，本文算法依据马尔科夫毯的相关性质，充分利用LiNGAM模型的因果顺序，减少了所需条件独立性测试的次数，较好地解决了将PC等算法作为剪枝算法所存在的误剪枝率高的问题。实验表明，与LiNGAM算法和DirectLiNGAM算法所采用的剪枝算法相比，本文的剪枝算法具有更高的准确率和稳定性。

本文算法在模拟数据集上进行了充分的对比实验，得到了较优的结果。下一步的重要研究内容是本文算法在真实数据上的表现。同时，由于高维度下正确的因果顺序较难获得，高维空间的因果顺序推断算法与不依赖于因果顺序的剪枝算法也是未来的主要研究方向。

[2] Hoyer P O,Janzing D,Mooij J M,et al.Nonlinear causal discovery with additive noise models[C]//Advances in neural information processing systems.2009:689-696.

[5] Shimizu S,Inazumi T,Sogawa Y,et al.DirectLiNGAM:A direct method for learning a linear non-gaussian structural equation model[J].The Journal of Machine Learning Research,2011,12(2):1225-1248.

[6] Glymour C,Spirtes P,Scheines R.Causal inference[J].Erkenntnis,1991,35(1-3):151-189.

[7] Spirtes P,Glymour C N,Scheines R.Causation,prediction,and search [M].MIT press,2000.

[8] Pearl J.Causality:models,reasoning and inference [M].Cambridge:MIT press,2000.

[9] Cai R,Zhang Z,Hao Z.BASSUM:A Bayesian semi-supervised method for classification feature selection[J].Pattern Recognition,2011,44(4):811-820.

[10] Aliferis C F,Tsamardinos I,Statnikov A.HITON:a novel Markov Blanket algorithm for optimal variable selection[C]//AMIA Annu Symp Proc.,2003:21-25.

LINGAM MODEL PRUNING ALGORITHM BASED ON CONDITIONAL INDEPENDENCE

Hao ZhifengLü HongweiCai RuichuYuan Chang

(FacultyofComputerScience,GuangdongUniversityofTechnology,Guangzhou510006，Guangdong,China)

How to conjecture causal network structure according to the observed data is an important problem in the field of statistics and machine learning.Quite a few research achievements have been gained by scholars in recent years,among them LiNGAM algorithm is a classical linear causal inference algorithm.However the pruning policy employed in LiNGAM algorithm requires high runtime complexity and provides poor accuracy on sparse graphs.Therefore in this paper we present a novel pruning method to solve this problem,it is based on conditional independence testing.The algorithm first rearranges the variables according to causal order and then employs partial correlation coefficient to check the conditional independence between variables according to new order.Numerous experimental results indicate that the pruning algorithm based on conditional independence proposed in the paper achieves higher accuracy with better running time on sparse graph than the one of LiNGAM.

Linear causalityPartial correlationConditional independencePruning method

2015-03-12。国家自然科学基金项目(61472089)。郝志峰，教授，主研领域：算法设计与分析，组合优化与算法研究，仿生算法的数学理论，代数学及其应用。吕宏伟，硕士。蔡瑞初，副教授。袁畅，硕士。

TP3

A

10.3969/j.issn.1000-386x.2016.08.056