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

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

0 引 言

隨著機(jī)器學(xué)習(xí)、深度學(xué)習(xí)等人工智能算法的崛起，人們逐漸意識(shí)到不能單純地關(guān)注數(shù)據(jù)的相關(guān)性。2011年圖靈獎(jiǎng)得主Judea Pearl[1]認(rèn)為人工智能深陷關(guān)聯(lián)概率的泥潭從而忽視了因果，同時(shí)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]、社會(huì)媒體應(yīng)用等[6]領(lǐng)域都有相關(guān)應(yīng)用。因此因果推理正在成為機(jī)器學(xué)習(xí)領(lǐng)域一個(gè)越來越受關(guān)注的研究方向?，F(xiàn)階段的因果發(fā)現(xiàn)主要是在不進(jìn)行任何干預(yù)的情況下確定兩個(gè)變量之間的因果方向，這是一項(xiàng)具有挑戰(zhàn)性的工作。在近階段的因果發(fā)現(xiàn)研究進(jìn)展中允許在某種假設(shè)下，基于純粹觀察數(shù)據(jù)推斷兩個(gè)變量之間的因果方向。然而在現(xiàn)實(shí)世界中觀察到的數(shù)據(jù)往往由多種生成方式產(chǎn)生，這使得傳統(tǒng)的因果推斷方法存在識(shí)別率不高和穩(wěn)定性較差的現(xiàn)象。

1 相關(guān)工作

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

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

2.1 模型假設(shè)

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

圖1 ANM_MM模型

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

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

二元變量的獨(dú)立性測(cè)試方法有很多，文中使用的是HilberSchmidt獨(dú)立性準(zhǔn)則。假設(shè)X，Y是兩個(gè)一維變量，定義一個(gè)非線性映射φ(x)∈F，F(xiàn)屬于再生核希爾伯特空間(reproducing kernel Hilbert space)且x∈X，再定義一個(gè)ψ(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)

對(duì)于所有的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í)別性和參數(shù)估計(jì)

3.1 模型的可識(shí)別性

(6)

圖2 因果推斷示意圖

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

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

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

λlogHSICb(X,Θ)]

(7)

由于在求解分布參數(shù)的過程中沒有考慮到結(jié)果變量對(duì)分布參數(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)λ增大時(shí)即增強(qiáng)X和分布參數(shù)的獨(dú)立性，由于結(jié)果變量前又有(1-λ)進(jìn)行約束，則Y和分布參數(shù)獨(dú)立性勢(shì)必會(huì)減小，因此并未違反因果邏輯性。當(dāng)α=0時(shí)，文中使用的目標(biāo)函數(shù)(8)和ANM-MM的目標(biāo)函數(shù)(7)具有一致性。由于該目標(biāo)函數(shù)包含兩個(gè)方向上的獨(dú)立性測(cè)試，因此文中算法記作Dual-ANMM，另外當(dāng)α=1時(shí)目標(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個(gè)神經(jīng)元；最后一層是激活層，文中使用的是LRelu激活。最后將通過梯度下降法求得目標(biāo)函數(shù)(3)的最小值從而得到混合數(shù)據(jù)的分布參數(shù)。再通過HilberSchmidt獨(dú)立性測(cè)試確定因果方向，如圖3所示。

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

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

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

當(dāng)HSIC(X,θ1)=HSIC(Y,θ2)時(shí)無法識(shí)別，第三種情況在實(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)因果推斷方法對(duì)比

LiNGAM模型[10]是一種線性非高斯因果推斷模型，LiNGAM模型假設(shè)數(shù)據(jù)由E=βC+N產(chǎn)生，其中β屬于實(shí)數(shù)域，C⊥N并且N是非高斯的。LiNGAM僅適用于非高斯噪聲的線性關(guān)系，因此如果生成的數(shù)據(jù)違反了該假設(shè)，LiNGAM模型會(huì)表現(xiàn)得極差。如果數(shù)據(jù)之間的關(guān)系是線性函數(shù)的關(guān)系，識(shí)別率會(huì)提高，而且計(jì)算成本相對(duì)較低。在實(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算法識(shí)別因果關(guān)系使用原因和結(jié)果之間的不對(duì)稱來檢測(cè)，原因和結(jié)果之間的不對(duì)稱是通過原因和殘差之間獨(dú)立的假設(shè)實(shí)現(xiàn)。因此該方法需要擬合回歸函數(shù)，并對(duì)輸入與的選擇對(duì)性能識(shí)別率是至關(guān)重要的。在實(shí)驗(yàn)中，文中使用了Mooij等人(2016)提出的一種實(shí)現(xiàn)，該實(shí)現(xiàn)使用高斯過程回歸進(jìn)行預(yù)測(cè)，并使用熵估計(jì)進(jìn)行輸入和殘差之間的關(guān)系評(píng)估。

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

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

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

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

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

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

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

該實(shí)驗(yàn)對(duì)比了幾種不同的傳統(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ù)量增加時(shí)結(jié)果逼近文中提出的Dual-AMMM算法以及ANM-MM模型，由上圖可以看出Dual-AMMM和ANM-MM原始算法和在穩(wěn)定性和準(zhǔn)確率上都優(yōu)于其它方法。

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

該實(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ù)集上的對(duì)比

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é)果變量對(duì)分布參數(shù)的影響，在實(shí)驗(yàn)上也表現(xiàn)出較其他傳統(tǒng)的方法具有較好的穩(wěn)定性和準(zhǔn)確率。但是現(xiàn)階段的研究主要是在二元變量之間的因果方向推斷，因此將該方法推廣到多元變量之間的因果圖的識(shí)別是下一階段的主要內(nèi)容。