基于線性非高斯無環(huán)模型的魯棒因果發(fā)現(xiàn)算法

曾 艷，郝志峰，2，蔡瑞初，謝 峰

(1.廣東工業(yè)大學(xué)計算機(jī)學(xué)院， 廣東 廣州 510006；2.佛山科學(xué)技術(shù)學(xué)院數(shù)學(xué)與大數(shù)據(jù)學(xué)院，廣東 佛山 528000；3.北京大學(xué)數(shù)學(xué)科學(xué)學(xué)院，北京 100080)

1 引言

探索數(shù)據(jù)背后的因果機(jī)制已經(jīng)成為當(dāng)前的熱點問題，同時也是計算機(jī)仿真的重要研究內(nèi)容，被廣泛應(yīng)用于機(jī)器學(xué)習(xí)等領(lǐng)域中[1]-[4]。傳統(tǒng)的因果發(fā)現(xiàn)算法大致分為兩類：基于約束的算法，如PC[5]和基于評分的算法，如GES[6]。雖然上述方法[5]，[6]已經(jīng)被用于許多重要領(lǐng)域，但是其不能夠找到一個完整的有向無環(huán)圖(Directed Acyclic Graph，DAG)。Shimizu等人[7]，[8]提出的線性非高斯無環(huán)模型(Linear Non-Gaussian Acyclic Model，LiNGAM)很好地解決了上述問題，該模型能夠唯一地識別出完整的DAG并成為最重要的因果模型之一?，F(xiàn)有估計LiNGAM的算法包括了基于獨立成分分析的ICA-LiNGAM[7]、基于獨立性分析的DirectLiNGAM[9]算法以及Pairwise-LiNGAM算法[10]、GPL算法[11]等。最近Xie等人[12]利用信息熵提出了ETPIA算法。

然而，上述因果發(fā)現(xiàn)方法[7]-[12]由于搜索方式不恰當(dāng)存在算法性能下降的問題。也即是說， 它們采取全局的搜索方式，從整個變量集合出發(fā)，在整個變量集合中逐層地選擇外生變量(請見定義1)，確定因果次序，進(jìn)而獲得整個因果網(wǎng)絡(luò)。而且，在逐層選擇外生變量的過程中，這些方法均需要將每一個變量與其余所有變量成對地進(jìn)行對比。因此，當(dāng)數(shù)據(jù)的維度增高時，對比的變量個數(shù)將增大，累積誤差也會相應(yīng)增加，使得在逐層選擇外生變量的步驟上易出錯。同理，當(dāng)樣本量不足時，額外的估計誤差也會引入。錯誤地選擇外生變量會產(chǎn)生級聯(lián)效應(yīng)，使得整個網(wǎng)絡(luò)的估計誤差會隨著層數(shù)的增大而變得越來越大，導(dǎo)致算法性能大幅度降低。

因此，本文基于LiNGAM，提出了一種魯棒的因果發(fā)現(xiàn)算法。不同于上述方法[7]-[12]，文中的算法采用d分離準(zhǔn)則確定因果網(wǎng)絡(luò)骨架，根據(jù)網(wǎng)絡(luò)骨架提供的相鄰節(jié)點信息(請見定義2)，采取局部搜索的方式依次尋找外生變量。算法的主要特點主要體現(xiàn)在搜索空間的減少，即，當(dāng)數(shù)據(jù)的變量個數(shù)增加時，算法在選擇外生變量時的搜索空間不是集中在整個變量集合，而是在每個節(jié)點的相鄰節(jié)點集合，因此減小了搜索空間，從而使得算法能在數(shù)據(jù)的維度增大或樣本量不足時展示較好的魯棒性。此外，本文提供了算法的理論證明。且不同于以往方法的理論[9]-[12]，證明了在引入相鄰節(jié)點的信息后外生變量的最強(qiáng)獨立性以及算法的一致性。同時，該有效性也通過仿真得到驗證。

2 線性非高斯無環(huán)模型

線性非高斯無環(huán)模型(Linear Non-Gaussian Acyclic Model， LiNGAM)是因果關(guān)系發(fā)現(xiàn)中的重要模型[]，其基本形式為

(1)

其中i={1，…，n}，xi表示觀測變量；bij表示有向無環(huán)圖DAG中變量xj到xi的連接強(qiáng)度；ei表示噪聲變量，服從非零方差的非高斯分布，且彼此間獨立；K(i)表示變量xi的因果次序。

LiNGAM需滿足以下三個基本假設(shè)：1)數(shù)據(jù)產(chǎn)生過程滿足式(1)，即線性假設(shè)；2)因果充分性假設(shè)，即不存在隱藏變量；3)噪聲變量ei服從非零方差的非高斯分布。

將(1)式寫成矩陣的形式如下：

X=BX+e

(2)

其中X=[x1，…，xn]T，e=[e1，…，en]T，連接矩陣B存儲了變量xi與變量xj之間的連接強(qiáng)度bij。值得注意的是，由于DAG的性質(zhì)，連接矩陣可以通過對應(yīng)的行列變換，轉(zhuǎn)換成一個嚴(yán)格的下三角矩陣[7](嚴(yán)格的下三角矩陣指的是一個對角線上的元素均為0的下三角矩陣)。

本文的目標(biāo)是如何有效地估計LiNGAM模型。即在已知觀察數(shù)據(jù)X的情況下，求解因果次序K和連接矩陣B，進(jìn)而得到完整的因果網(wǎng)絡(luò)結(jié)構(gòu)。

3 S2-LiNGAM因果發(fā)現(xiàn)算法

3.1 基本定義

在闡述本算法的基本思想與理論分析之前，為了便于描述，首先給出外生變量、相鄰節(jié)點和d分離準(zhǔn)則的基本定義。

定義1[9](外生變量)：外生變量是指在因果關(guān)系網(wǎng)絡(luò)中只起解釋作用的變量，即在因果模型內(nèi)部沒有直接的原因節(jié)點。

外生變量有時也叫外生變量。由定義1可知，在具體的因果關(guān)系網(wǎng)絡(luò)圖中，外生變量只指向其它節(jié)點，而不存在其它節(jié)點指向外生變量。

定義2[5](相鄰節(jié)點)：在因果關(guān)系網(wǎng)絡(luò)中，如果節(jié)點xi與節(jié)點xj由一條有向邊直接相連，則稱節(jié)點xi是節(jié)點xj的相鄰節(jié)點，或節(jié)點xj是節(jié)點xi的相鄰節(jié)點。

由定義2可知，相鄰節(jié)點包括了一個節(jié)點所有子節(jié)點與父節(jié)點。

定義3[5](d分離準(zhǔn)則)：在因果關(guān)系網(wǎng)絡(luò)中，一條路徑p被節(jié)點集Zd分離當(dāng)且僅當(dāng)

(a)p包含一條鏈i→m→j或者碰撞點i→m←j使得m∈Z，或者

(b)p包含碰撞點i→m←j使得m?Z并且不存在節(jié)點m的后代節(jié)點屬于Z。

換句話說，節(jié)點集Zd分離X和Y當(dāng)且僅當(dāng)Z阻隔了從X到Y(jié)的所有路徑。由定義3可知，d分離準(zhǔn)則能被用于發(fā)現(xiàn)因果網(wǎng)絡(luò)骨架，從而獲得每個節(jié)點的相鄰節(jié)點信息。

3.2 算法的理論分析

在這一小節(jié)中，將給出外生變量的判定定理(請見定理2)，再給出算法的一致性證明(請見定理3)，保證算法的有效性。在描述定理2和3前，先提供定理1，D-S定理，因為它會被用在后續(xù)的定理證明中。

定理1[9](D-S定理)：定義兩個隨機(jī)變量x1和x2是由獨立的隨機(jī)變量εi，i=1，2，…，q線性組合而成，即

(3)

那么，如果x1和x2統(tǒng)計獨立，則對于αiβi≠0，所有的εi是高斯分布。換句話說，定理1說明如果存在服從非高斯分布的εi滿足αiβi≠0，那么x1和x2統(tǒng)計依賴。

(4)

證明：從1)充分性與2)必要性證明。

xj=ej

xi=bijej+ei

其中，每一個節(jié)點xk都(包括xi)可以表示為噪音變量e的線性組合(噪音變量不包括ej)

(5)

(6)

定理1描述了給定相鄰節(jié)點信息下的外生變量性質(zhì)，它幫助本文對因果網(wǎng)絡(luò)結(jié)構(gòu)中的外生變量進(jìn)行判定，同時它為本文逐層選擇外生變量的正確性提供了初步的理論基礎(chǔ)。此外，在選擇第一個外生變量后，本文需要更新原始觀察數(shù)據(jù)，再繼續(xù)選擇下一個外生變量，進(jìn)而最終確定因果次序。定理3為下一個外生變量的正確識別，即因果次序的正確識別，提供了完整的理論支撐。

3.3 算法流程

基于上述的理論分析，本文提出了基于LiNGAM的魯棒因果關(guān)系發(fā)現(xiàn)算法(A RobuSt CauSal Discovery Algorithm Based on LiNGAM)簡稱為：S2-LiNGAM。算法的基本步驟如下：

輸入：觀察數(shù)據(jù)X={x1，x2，…，xn}。

輸出：因果次序K和因果網(wǎng)絡(luò)B。

1) 初始化因果次序集合K=?與集合U={1，…，n}；

2) 利用d分離準(zhǔn)則，尋找每個節(jié)點的相鄰節(jié)點，并存儲在相鄰矩陣D={d1，…，dn}中；

3) Repeat

6) 根據(jù)式(4)，更新數(shù)據(jù)X，為每個節(jié)點除去外生變量對其產(chǎn)生的影響，同時，U=U/{j}；

7) Until U的節(jié)點個數(shù)為1；

8) 將最后一個節(jié)點的下標(biāo)加入至K中；

9) 根據(jù)得到的因果次序K和相鄰矩陣D，使用剪枝算法估計連接矩陣B。為了更清晰地說明流程，提供例子如下。

圖1 一個簡單LiNGAM模型的因果網(wǎng)絡(luò)圖

x3⊥r(3)1,x3⊥r(3)4,x2⊥r(2)1,x2⊥

r(2)4,x1⊥

r(1)2,x1⊥

4 仿真研究

為了說明算法的有效性，本文從貝葉斯網(wǎng)絡(luò)存儲庫中選擇了6種不同領(lǐng)域下的真實網(wǎng)絡(luò)結(jié)構(gòu)來進(jìn)行仿真，如表1所示(1)真實網(wǎng)絡(luò)的詳細(xì)信息請見：http：∥www.cs.huji.ac.il/～galel/Repository/。首先，本文根據(jù)真實網(wǎng)絡(luò)結(jié)構(gòu)產(chǎn)生服從LiNGAM的數(shù)據(jù)，其中噪聲變量服從非高斯的拉普拉斯分布，連接矩陣B從[-1，-0.5]∪[0.5，1]的均勻分布隨機(jī)產(chǎn)生。選取常用的ICA-LiNGAM算法[7]， Pairwise-LiNGAM算法[10]以及ETPIA算法[12]作為本實驗的對比方法。此外，針對衡量標(biāo)準(zhǔn)，本文選擇了貝葉斯因果網(wǎng)絡(luò)的常用指標(biāo)召回率(Recall)、準(zhǔn)確率(Precision)和F1得分(F1 score)。F1得分反映了因果網(wǎng)絡(luò)結(jié)構(gòu)的總體估計優(yōu)劣，具體為

Recall=TP/(TP+FN)，

Precision=TP/(TP+FP)，

其中，TP表示正確估計的邊的個數(shù)；FP表示未估計到的邊的個數(shù)；FN表示錯誤估計的邊的個數(shù)。結(jié)果取20次的平均值。

表1 真實網(wǎng)絡(luò)的介紹

本文分別為6種不同大小的真實網(wǎng)絡(luò)設(shè)計了樣本量在50、100、200、500下的對比實驗，結(jié)果如圖2至圖7所示。從總體來看，隨著樣本量的逐漸增加，本算法S2-LiNGAM在不同網(wǎng)絡(luò)下的召回率、準(zhǔn)確率與F1得分都逐漸上升，說明了算法理論的正確性。特別地，當(dāng)樣本量小(=50)或者節(jié)點數(shù)量高(WIN95PTS的維度為76)時，算法的F1得分都明顯優(yōu)于其余三個算法，驗證了算法良好的魯棒性，同時也說明了本文采取的局部搜索選擇外生變量的方式的有效性。雖然S2-LiNGAM的召回率在Barley等部分網(wǎng)絡(luò)結(jié)構(gòu)上不及對比的算法，但其準(zhǔn)確率和F1得分都大大地超過于對比算法，這說明了算法的剪枝操作學(xué)習(xí)到的邊盡管少，但大多數(shù)是準(zhǔn)確的。因此，在現(xiàn)實應(yīng)用中，可以根據(jù)實際情況選擇合適的剪枝操作，從而得到更有效的因果網(wǎng)絡(luò)結(jié)構(gòu)。Pairwise-LiNGAM和ETPIA算法是基于Direct-LiNGAM框架的，它們采取全局搜索的方式選擇外生變量和學(xué)習(xí)因果網(wǎng)絡(luò)結(jié)構(gòu)，因此在樣本量不足或維度增高時算法性能降低，F(xiàn)1得分較低。

綜上分析，相較于其余三個算法，本方法都能夠取得較好的效果，驗證了S2-LiNGAM在因果關(guān)系發(fā)現(xiàn)中具有良好的魯棒性。

圖2 四種算法在ALARM真實結(jié)構(gòu)上的實驗效果

圖3 四種算法在MILDEW真實結(jié)構(gòu)上的實驗效果

圖4 四種算法在BARLEY真實結(jié)構(gòu)上的實驗效果

圖5 四種算法在HEPAR2真實結(jié)構(gòu)上的實驗效果

圖6 四種算法在WIN95PTS真實結(jié)構(gòu)上的實驗效果

圖7 四種算法在HAILFINDER真實結(jié)構(gòu)上的實驗效果

5 結(jié)束語

本文基于LiNGAM模型，提出了一種魯棒的因果關(guān)系發(fā)現(xiàn)算法。具體的創(chuàng)新包括：

1)算法充分融合了d分離準(zhǔn)則和局部搜索方法的優(yōu)勢，能高效處理高維或樣本量不足的數(shù)據(jù)，展現(xiàn)了較好的魯棒性；

2)本文分析了外生變量最強(qiáng)獨立性的性質(zhì)，并提供了算法的一致性證明，從理論上保證了算法的有效性；

3)基于真實因果結(jié)構(gòu)的實驗驗證了算法的正確性與魯棒性。特別地，當(dāng)樣本量低至50或維度高至76時，算法的性能優(yōu)勢明顯。