基于條件獨立性的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