基于神經(jīng)網(wǎng)絡(luò)的混合數(shù)據(jù)的因果發(fā)現(xiàn)

耿家興，萬亞平,2，李洪飛

(1.南華大學(xué) 計(jì)算機(jī)學(xué)院，湖南 衡陽 421001；2.中核集團(tuán)高可信計(jì)算重點(diǎn)學(xué)科實(shí)驗(yàn)室，湖南 衡陽 421001)

0 引 言

隨著機(jī)器學(xué)習(xí)、深度學(xué)習(xí)等人工智能算法的崛起，人們逐漸意識到不能單純地關(guān)注數(shù)據(jù)的相關(guān)性。2011年圖靈獎得主Judea Pearl[1]認(rèn)為人工智能深陷關(guān)聯(lián)概率的泥潭從而忽視了因果，同時Judea Pearl認(rèn)為研究者應(yīng)該研究因果，因?yàn)檫@可能是實(shí)現(xiàn)真正的人工智能的可能路徑。隨著因果推斷的深入研究，基于因果推斷的應(yīng)用也越來越多，例如在fMRI數(shù)據(jù)上的分析和分子發(fā)現(xiàn)途徑研究；Statnikov等人[2]在2012提出了基因組學(xué)數(shù)據(jù)中因果分子相互作用的鑒定方法；Ma Sisi在2017年的工作中討論了生物醫(yī)學(xué)中觀察數(shù)據(jù)可用于因果發(fā)現(xiàn)[3]，以及在網(wǎng)絡(luò)性能優(yōu)化[4]、網(wǎng)購服裝推薦[5]、社會媒體應(yīng)用等[6]領(lǐng)域都有相關(guān)應(yīng)用。因此因果推理正在成為機(jī)器學(xué)習(xí)領(lǐng)域一個越來越受關(guān)注的研究方向。現(xiàn)階段的因果發(fā)現(xiàn)主要是在不進(jìn)行任何干預(yù)的情況下確定兩個變量之間的因果方向，這是一項(xiàng)具有挑戰(zhàn)性的工作。在近階段的因果發(fā)現(xiàn)研究進(jìn)展中允許在某種假設(shè)下，基于純粹觀察數(shù)據(jù)推斷兩個變量之間的因果方向。然而在現(xiàn)實(shí)世界中觀察到的數(shù)據(jù)往往由多種生成方式產(chǎn)生，這使得傳統(tǒng)的因果推斷方法存在識別率不高和穩(wěn)定性較差的現(xiàn)象。

1 相關(guān)工作

目前在二元變量上的因果推斷方法主要有兩種，第一種是基于時序數(shù)據(jù)的因果方向推斷，第二種是基于非時序數(shù)據(jù)的因果關(guān)系發(fā)現(xiàn)，現(xiàn)階段研究比較多的是基于非時序數(shù)據(jù)的因果推斷。

非時序數(shù)據(jù)上的因果關(guān)系發(fā)現(xiàn)又分為基于統(tǒng)計(jì)分析的因果推斷和基于模型假設(shè)的因果推斷。基于統(tǒng)計(jì)分析：Kano等人[7]于2003年提出路徑分析法用于觀測數(shù)據(jù)的因果結(jié)構(gòu)，該方法通過將隱含協(xié)方差矩陣與樣本協(xié)方差矩陣進(jìn)行比較，對假設(shè)的因果關(guān)系進(jìn)行統(tǒng)計(jì)檢驗(yàn)，從而得出假設(shè)的因果關(guān)系的正確性。2003年Comley[8]提出了一種對稱貝葉斯網(wǎng)絡(luò)通過使用一種簡單的決策樹算法能夠很好地識別連續(xù)型和離散型的聯(lián)合概率分布。Sun等人[9]于2006年提出一種使用馬爾可夫核函數(shù)來選擇多變量之間的所有假設(shè)因果方向中的最可能的一個方向?；谀Ｐ图僭O(shè)：Shimizu[10]在2006年提出了一種線性非高斯非循環(huán)模型LinGAM，該算法假設(shè)數(shù)據(jù)是線性非混淆的，同時干擾變量是非高斯的，該模型使用獨(dú)立成分分析的方法來進(jìn)行因果分析。Zhang等人[11]在2009年通過考慮原因的非線性效應(yīng)，內(nèi)部噪聲效應(yīng)和觀測變量中的測量失真效應(yīng)，提出了后非線性因果模型PNL，并證明了該方法在因果方向的可識別性。在2009年Hoyer提出[12]加性噪聲模型ANM，該模型能夠處理非線性的數(shù)據(jù)，并成為了一種常用的因果發(fā)現(xiàn)方法。Janzing等人[13]提出通過二階指數(shù)模型定義一系列平滑密度和條件密度，即通過最大化受第一和第二統(tǒng)計(jì)矩影響的條件熵來識別包含離散變量和連續(xù)變量的因果結(jié)構(gòu)。2010年Daniusis等人[14]提出一種確定兩個可逆函數(shù)相互關(guān)聯(lián)性來識別確定的因果方向。2011年Jonas[1]提出了利用加性噪聲模型對離散型數(shù)據(jù)的因果推斷算法。2012年Janzing[15]提出了IGCI算法，該算法通過信息空間中的正交性來定義獨(dú)立性的方法描述因果變量值之間的邊緣概率分布和條件概率分布，從而通過因果關(guān)系的不對稱性來確定因果方向。2015年Sgouritsa等人[16]提出了無監(jiān)督的高斯回歸過程來估計(jì)因果變量之間的條件概率CURE算法。2016年Mooij等人[17]總結(jié)了使用觀測數(shù)據(jù)區(qū)分因果關(guān)系的方法和標(biāo)準(zhǔn)。2017年Marx等人[18]采用了一種基于柯爾莫戈羅夫復(fù)雜度的信息論方法，并利用最小值描述長度原理提供了一種實(shí)用的因果關(guān)系發(fā)現(xiàn)方法SLOPE。以上這些方法在進(jìn)行二元變量的混合數(shù)據(jù)的因果推斷中表現(xiàn)較差，2018年Hu[19]提出了一種基于加性噪聲算法的混合數(shù)據(jù)的因果推斷方法ANM-MM(additive noise model-mixture model)模型。

ANM-MM模型可以處理二元變量的混合數(shù)據(jù)的因果推斷，但是在求解分布參數(shù)的過程中沒有考慮到結(jié)果變量對分布參數(shù)的影響。因此文中在ANM-MM模型的基礎(chǔ)上提出了一種使用神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)改進(jìn)的目標(biāo)函數(shù)從而得到分布參數(shù)，原ANM-MM模型的目標(biāo)函數(shù)被看作是該目標(biāo)函數(shù)的一種特例。

2 預(yù)備知識

2.1 模型假設(shè)

文中提出的方法是在ANM-MM模型的假設(shè)下進(jìn)行的，ANM-MM模型的定義是連續(xù)性變量X和Y由有限ANM模型生成，如圖1所示。

圖1 ANM_MM模型

ANM_MM模型中所有的觀測數(shù)據(jù)都是由同一個形式的函數(shù)f不同的參數(shù)θ生成，因?yàn)樵诿總€獨(dú)立的試驗(yàn)中，由于某些難以控制的外部因素的影響，數(shù)據(jù)生成過程略有不同。此外這些因素通常被認(rèn)為是獨(dú)立于觀測變量的。

2.2 HilberSchmidt獨(dú)立性準(zhǔn)則

二元變量的獨(dú)立性測試方法有很多，文中使用的是HilberSchmidt獨(dú)立性準(zhǔn)則。假設(shè)X，Y是兩個一維變量，定義一個非線性映射φ(x)∈F，F(xiàn)屬于再生核希爾伯特空間(reproducing kernel Hilbert space)且x∈X，再定義一個ψ(y)∈G，G同樣屬于再生核希爾伯特空間且y∈Y，則相應(yīng)的核函數(shù)分別記為：

k(x,x')=〈φ(x),φ(x')〉,x,x'∈X

(1)

l(y,y')=〈ψ(y),ψ(y')〉,y,y'∈Y

(2)

對于所有的f∈F,g∈G定義：

〈f,Cxyg〉F=Exy([f(x)-Ex(f(x))][g(y)-

Ey(g(y))])

(3)

交叉協(xié)方差算子Cxy:G→F且：

Cxy=Exy[(φ(x)-μx)?(ψ(y)-μy)]

(4)

其中，μx=Exφ(x)，μy=Eyφ(y)，?為張量積，Cxy可以看作Hilbert-Schmidt算子，然后將Frobenius范數(shù)擴(kuò)展到該算子上，得到Hilbert-Schmidt范數(shù)，將該范數(shù)平方即可得到Hilbert-Schmidt獨(dú)立性準(zhǔn)則：

(5)

其中，Prxy是X，Y的概率分布。

3 模型可識別性和參數(shù)估計(jì)

3.1 模型的可識別性

(6)

圖2 因果推斷示意圖

根據(jù)定理可以通過研究假設(shè)原因與相應(yīng)函數(shù)參數(shù)之間的獨(dú)立性，推斷出模型的因果方向，如果它們在因果方向上是獨(dú)立的，那么它們很可能在反因果方向上是相互依賴的。因此在實(shí)際中兩變量間的因果推斷的方向是變量和中間參數(shù)之間獨(dú)立性更大的方向。

3.2 參數(shù)估計(jì)

在ANM-MM模型中利用目標(biāo)函數(shù)的對數(shù)似然極大化通過梯度下降來尋找潛在的分布參數(shù)的表示形式，然而θ不能直接通過極大似然估計(jì)進(jìn)行求解，所以ANM-MM需要添加X和分布參數(shù)θ之間的HSIC獨(dú)立性到損失函數(shù)得到：

λlogHSICb(X,Θ)]

(7)

由于在求解分布參數(shù)的過程中沒有考慮到結(jié)果變量對分布參數(shù)的影響，文中修改了ANM-MM模型的損失函數(shù)(7)，另外添加Y和分布參數(shù)θ之間的HSIC獨(dú)立性到損失函數(shù)得到：

λlogHSICb(X,Θ)+

α(1-λ)logHSICb(Y,Θ)]

(8)

其中，λ控制原因變量X和分布參數(shù)之間的Hilbert-Schmidt獨(dú)立性，α用來調(diào)節(jié)結(jié)果變量Y和分布參數(shù)的Hilbert-Schmidt獨(dú)立性，當(dāng)λ增大時即增強(qiáng)X和分布參數(shù)的獨(dú)立性，由于結(jié)果變量前又有(1-λ)進(jìn)行約束，則Y和分布參數(shù)獨(dú)立性勢必會減小，因此并未違反因果邏輯性。當(dāng)α=0時，文中使用的目標(biāo)函數(shù)(8)和ANM-MM的目標(biāo)函數(shù)(7)具有一致性。由于該目標(biāo)函數(shù)包含兩個方向上的獨(dú)立性測試，因此文中算法記作Dual-ANMM，另外當(dāng)α=1時目標(biāo)函數(shù)改寫為：

λlogHSICb(X,Θ)+

(1-λ)logHSICb(Y,Θ)]

(9)

其中，目標(biāo)函數(shù)(8)引入α的目的是控制數(shù)據(jù)和參數(shù)之間HilberSchmidt獨(dú)立性的強(qiáng)弱。文中在實(shí)驗(yàn)中使用神經(jīng)網(wǎng)絡(luò)來學(xué)習(xí)，采取梯度下降法最優(yōu)化α=1的目標(biāo)函數(shù)(9)得到混合數(shù)據(jù)的抽象因果分布參數(shù)。

4 算法過程

文中使用了三層網(wǎng)絡(luò)結(jié)構(gòu)來學(xué)習(xí)混合數(shù)據(jù)的分布參數(shù)，第一層是Normlization layer，該層網(wǎng)絡(luò)主要將輸入數(shù)據(jù)正則化、規(guī)范化；第二層是Full-Connect全連接層，文中采用的是10個神經(jīng)元；最后一層是激活層，文中使用的是LRelu激活。最后將通過梯度下降法求得目標(biāo)函數(shù)(3)的最小值從而得到混合數(shù)據(jù)的分布參數(shù)。再通過HilberSchmidt獨(dú)立性測試確定因果方向，如圖3所示。

圖3 網(wǎng)絡(luò)架構(gòu)

當(dāng)HSIC(X,θ1)>HSIC(Y,θ2)時，X→Y；

當(dāng)HSIC(X,θ1)

當(dāng)HSIC(X,θ1)=HSIC(Y,θ2)時無法識別，第三種情況在實(shí)際情況下很少發(fā)生。文中使用的HilberSchmidt獨(dú)立性是Hu Shoubo改進(jìn)的GPPOM_HSIC獨(dú)立性方法。

Dual-ANMM算法過程如下：

Algorithm Dual-ANMM(A kind of causal discovery by neural network)

輸入：D，alpha=1，lamda//其中D是數(shù)據(jù)，alpha和lamda為超參數(shù)

輸出：direct//輸出方向

1.Standardize observation of eachD;

2.Random initialize weight of NN and other argument

3.Optimize target function by NN inX->Yestimate theta1

4.D1=HSIC(theta1,X)

5.Optimize target function by NN inY->Xestimate theta2

6.D2=HSIC(theta2,X)

7.IfD1

8.ElifD1

9.Else direct =0

10.Return direct

5 傳統(tǒng)因果推斷方法對比

LiNGAM模型[10]是一種線性非高斯因果推斷模型，LiNGAM模型假設(shè)數(shù)據(jù)由E=βC+N產(chǎn)生，其中β屬于實(shí)數(shù)域，C⊥N并且N是非高斯的。LiNGAM僅適用于非高斯噪聲的線性關(guān)系，因此如果生成的數(shù)據(jù)違反了該假設(shè)，LiNGAM模型會表現(xiàn)得極差。如果數(shù)據(jù)之間的關(guān)系是線性函數(shù)的關(guān)系，識別率會提高，而且計(jì)算成本相對較低。在實(shí)驗(yàn)中，文中使用的LiNGAM是原始版本中基于獨(dú)立成分分析的算法。

ANM模型[12]是一種加性噪聲模型，ANM模型假設(shè)E=f(C)+N，這里的f是非線性函數(shù)并且C⊥N，該方法能夠處理非線性數(shù)據(jù)。由于現(xiàn)實(shí)世界中的數(shù)據(jù)大多是非線性數(shù)據(jù)，該方法得到了廣泛的應(yīng)用。ANM算法識別因果關(guān)系使用原因和結(jié)果之間的不對稱來檢測，原因和結(jié)果之間的不對稱是通過原因和殘差之間獨(dú)立的假設(shè)實(shí)現(xiàn)。因此該方法需要擬合回歸函數(shù)，并對輸入與的選擇對性能識別率是至關(guān)重要的。在實(shí)驗(yàn)中，文中使用了Mooij等人(2016)提出的一種實(shí)現(xiàn)，該實(shí)現(xiàn)使用高斯過程回歸進(jìn)行預(yù)測，并使用熵估計(jì)進(jìn)行輸入和殘差之間的關(guān)系評估。

PNL模型[11]是一種后線性模型，該模型是對ANM模型的一種泛化。PNL模型假設(shè)數(shù)據(jù)的生成機(jī)制是由E=g(f(C)+N)表示，其中g(shù)是非線性函數(shù)，f是線性函數(shù)，并且C⊥N。為了推斷出因果關(guān)系的方向，在兩個可能的方向上進(jìn)行PNL模型測試，并檢查輸入和擾動之間的獨(dú)立性。然而這里的擾動不同于回歸殘差，擬合PNL模型比擬合ANM要困難得多，準(zhǔn)確率較高于ANM模型。在實(shí)驗(yàn)中，文中使用了基于約束非線性獨(dú)立分量分析來估計(jì)擾動，并使用HilberSchmidt獨(dú)立性準(zhǔn)則進(jìn)行統(tǒng)計(jì)獨(dú)立測試。

6 實(shí)驗(yàn)過程

實(shí)驗(yàn)在CentOS系統(tǒng)上python3.6環(huán)境中進(jìn)行，分別在模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上實(shí)驗(yàn)，其他對比實(shí)驗(yàn)包括LiNGAM，IGCI，PNL，ANM分別按照上面實(shí)現(xiàn)方法在matlab2016a中進(jìn)行。SLOPE實(shí)驗(yàn)是R3.5.2中進(jìn)行。首先比較了λ對實(shí)驗(yàn)準(zhǔn)確率的影響，在實(shí)驗(yàn)中分別取λ為0.5，0.75，1.0，1.25，1.5，其中當(dāng)λ等于1.0的時候該模型等價于ANM-MM模型使用的損失函數(shù)。文中對混合數(shù)據(jù)類型個數(shù)對因果推斷的準(zhǔn)確率的影響進(jìn)行了實(shí)驗(yàn)，然后是數(shù)據(jù)量對ANM-MM，ANM，IGCI，PNL，SLOPE，LiNGAM的影響進(jìn)行了實(shí)驗(yàn)，最后在cause-effect pairs真實(shí)數(shù)據(jù)集上評估了算法的因果推理性能。cause-effect pairs具有用于測試因果檢測算法的不同數(shù)據(jù)，數(shù)據(jù)文件一共包含108個txt文件，每個文件包含兩個變量，一個是原因，另一個是效果。對于每個數(shù)據(jù)文件，都有一個描述性文件，可以在其中找到真實(shí)的因果關(guān)系以及數(shù)據(jù)來源描述，實(shí)驗(yàn)中108組數(shù)據(jù)中有9組被排除在外，因?yàn)樗鼈円从啥嘧兞繑?shù)據(jù)組成，要么由分類數(shù)據(jù)組成(對47、52、53、54，55、70、71、101和105)。每對實(shí)驗(yàn)重復(fù)100次獨(dú)立實(shí)驗(yàn)，記錄不同方法的正確推理率。然后計(jì)算同一數(shù)據(jù)集對的平均百分比作為對應(yīng)數(shù)據(jù)集的精度，實(shí)驗(yàn)結(jié)果如下。

如圖4所示，實(shí)驗(yàn)測試了λ對實(shí)驗(yàn)準(zhǔn)確率的影響，其中當(dāng)λ等于1.0時表示ANM_MM模型對應(yīng)的目標(biāo)函數(shù)，實(shí)驗(yàn)表明α=1.0時，λ=1.25時結(jié)果會達(dá)到最好。

如圖5所示，實(shí)驗(yàn)對比了混合數(shù)據(jù)類型個數(shù)對兩種算法的影響?？梢钥闯鑫闹刑岢龅乃惴ㄊ莾?yōu)于ANM-MM模型的，在實(shí)驗(yàn)中，當(dāng)只有一種因果機(jī)制的時候ANM-MM模型的準(zhǔn)確率低于65%，這是在100次重復(fù)實(shí)驗(yàn)中得到的結(jié)果，而文中提出的Dual-AMMM在參數(shù)α=1.0，λ=1.25時準(zhǔn)確率能達(dá)到77%，即當(dāng)數(shù)據(jù)并不是由多種因果機(jī)制生成，ANM-MM識別率低于Dual-AMMM。在多于一種混合類型的數(shù)據(jù)時，Dual-AMMM算法和ANM-MM模型結(jié)果不相上下。

圖4 λ對實(shí)驗(yàn)準(zhǔn)確率的影響

圖5 混合數(shù)據(jù)類型個數(shù)對準(zhǔn)確率的影響

該實(shí)驗(yàn)對比了幾種不同的傳統(tǒng)的因果推斷算法和ANM-MM模型在不同數(shù)據(jù)量下的模擬數(shù)據(jù)的準(zhǔn)確率，如圖6所示。實(shí)驗(yàn)結(jié)果表明，在不同的數(shù)據(jù)量下，該算法在穩(wěn)定性上優(yōu)于傳統(tǒng)的方法，其中ICGI在數(shù)據(jù)量增加時結(jié)果逼近文中提出的Dual-AMMM算法以及ANM-MM模型，由上圖可以看出Dual-AMMM和ANM-MM原始算法和在穩(wěn)定性和準(zhǔn)確率上都優(yōu)于其它方法。

圖6 數(shù)據(jù)量對模型的影響

該實(shí)驗(yàn)在真實(shí)因果數(shù)據(jù)集cause-effect pairs上進(jìn)行，結(jié)果如圖7所示，柱狀圖的高度表示實(shí)驗(yàn)100次的平均準(zhǔn)確率，其中Y的誤差取上下界誤差的最大。結(jié)果表明文中提出的Dual-AMMM算法優(yōu)于其他傳統(tǒng)的方法，平均準(zhǔn)確率最高約85%左右，ANM-MM算法在82%左右。PNL也表現(xiàn)很好，其次是IGCI算法，但是IGCI和SLOPE的誤差范圍較大。

圖7 在真實(shí)數(shù)據(jù)集上的對比

7 結(jié)束語

在ANM-MM的假設(shè)基礎(chǔ)上提出了一種基于神經(jīng)網(wǎng)絡(luò)的因果推斷方法。該方法使用梯度下降法最優(yōu)化損失函數(shù)得到混合數(shù)據(jù)的抽象因果分布參數(shù)，然后將分布參數(shù)看作是原因變量和結(jié)果變量之間的隱變量，通過比較原因變量和分布參數(shù)之間的HilberSchmidt獨(dú)立性來確定二元變量的因果方向。由于改進(jìn)的ANM-MM求解分布參數(shù)的過程中沒有考慮到結(jié)果變量對分布參數(shù)的影響，在實(shí)驗(yàn)上也表現(xiàn)出較其他傳統(tǒng)的方法具有較好的穩(wěn)定性和準(zhǔn)確率。但是現(xiàn)階段的研究主要是在二元變量之間的因果方向推斷，因此將該方法推廣到多元變量之間的因果圖的識別是下一階段的主要內(nèi)容。