基于偏相關性測試的遞歸式因果推斷算法

陳銘杰，張浩，彭昱忠，謝峰，龐悅

（1.東莞理工學院 計算機科學與技術學院，廣東 東莞 523808；2.廣東石油化工學院 計算機學院，廣東 茂名 525099；3.復旦大學 計算機科學技術學院，上海 200433；4.南寧師范大學 計算機與信息工程學院，南寧 530001；5.北京大學 數學科學學院，北京 100871；6.中國銀聯博士后科研工作站，上海 201201）

0 概述

在因果關系推斷問題中，研究者通常使用統計獨立性測試［1］和條件獨立性測試［2］對數據集間蘊含的因果關系進行推斷，因此，對于數據中因果信息的獲取尤為重要［3］。以PC 算法［4］為例，作為一個經典的基于約束的因果推斷算法，其算法原理是根據數據集間變量的聯合概率分布與條件概率分布特點，通過條件獨立性（Conditional Independence，CI）測試判斷變量之間是否滿足獨立性或條件獨立性，在進行一系列基于方向規則的判斷后構造出一個大致的因果網絡圖。此圖通常是一個部分有向無環圖（Partial Directed Acyclic Graph，PDAG），其中包含一組馬爾可夫（Markov）等價類，而不是一個完全準確的因果網絡圖。

因為CI 測試的條件集會較大程度影響到測試結果的準確率，所以CI 測試比統計獨立性測試的難度要大得多［5］。對于CI 測試，一般有三種使用場景：一種是數據離散分布的情況，可以根據對應的條件概率P(X，Y|Z)=P(X|Z)P(Y|Z)將條件集合Z變為單個條件變量，然后使用常用的條件互信息或卡方檢驗技術判斷條件獨立性關系是否成立［5］；另外兩種分別針對數據連續分布和混合連續離散分布的情況，會采用離散化的方法對數據進行離散化處理，然后通過傳統的CI 測試方法即第一種場景中的條件互信息和卡方檢驗等方法進行檢驗，但是存在的問題是在離散化的過程中容易丟失一些數據信息，導致最終推斷結果不夠準確。雖然可以通過計算條件密度估計值PX|YZ=PX|Z間的距離來檢驗CI 測試是否準確［6］，但是當數據量較大時，離散的區間比較難設定，過大容易導致數據丟失，過小則需要很多的樣本量，在Z較大的情況下，離散化數據和估算條件密度都是需要解決的重要問題。

近年來，研究者提出一系列基于核函數的CI 測試方法。這些方法可以利用高階矩信息將觀測變量映射到再生核希爾伯特空間（RKHSs），此時的變量具有幫助推斷高維變量分布的特點，如獨立性和同質性［7］。文獻［8］提出一種利用希爾伯特-施密特范數的條件交叉協方差算子，它是可以對應到RKHSs函數下x和y條件協方差的一種度量。該方法證實了如果與之對應的RKHSs 是特征核［9-11］，則算子范數為0且x⊥y|Z。相關文獻以及深入研究表明，基于核方法的CI 測試可以從給定的變量中獲得比基于離散數據集的CI 測試更多且更完整的信息，同時也表明基于核方法的CI 測試在因果推斷場景中能夠發現更準確的因果關系。

基于核函數的方法，通過將變量集的條件概率分布空間映射到基于平方可積函數空間的高維的線性空間，找出一種單向或雙向映射關系，如KCIT［12］。這種方法準確率較高，但運算復雜度也非常高，一般來說，當變量個數超過10 或樣本量大于2 000 時，就很難在一般的計算環境中運行。因此，該方法難以在高維數據場景中高效地解決因果推斷問題。由此可見，因果推斷的一個關鍵難題就是設計一個合適的CI 測試方法，使之能夠在高維數據場景下也能高效運行。

本文提出一種基于偏相關性測試的因果推斷算法CIPCT，從遞歸分治的角度應對數據集和條件集過大的情況。采用變量/特征分割方法將變量集遞歸式地分割成多個子數據集，使用因果推斷算法對每個子數據集進行因果推斷，得到對應的子結果。在此基礎上，通過合并策略對所有的子結果進行合并整理，得到與原始數據集相對應的因果網絡。此外，在“分-治-合”過程中采用較為高效的偏相關性測試，避免計算復雜度高的核密度估算，進一步提高因果推斷的效率。

1 相關研究

本節將介紹一些關于提升因果推斷速度的工作。在研究因果關系相關問題時，函數模型性能也是其中一個決定性因素。本文采用加性噪聲模型（Additive Noise Model，ANM）作為因果函數模型，以此作為函數模型的優勢在于：一是有獨特的模型特性；二是數據都有一個特殊的加噪聲項，可以借此獲取額外的信息［13］。由于基于最小殘差平方的回歸方法估算可得到關于Z的非線性平方可積函數，因此文獻［14］對（xi，Z）和（xj，Z）分別做非線性回歸得到函數f和g，通過檢驗f和g是否存在將CI 測試簡化為兩到三個統計獨立性測試［15］。當加性噪聲模型的特性被加以利用，即可避免估算條件概率密度，減少條件集Z所引起的高維難題，從而獲取完整的因果信息，加速CI 測試并打破Markov 等價類限制。

在打破Markov等價類限制的前提下，文獻［16-17］提出一種利用殘差獨立性發現條件獨立性的方法，其對于Z維度具有更好的魯棒性。Markov 等價類包含三個結構：x1→x2→x3，x1←x2→x3，x1←x2←x3，這三個結構有著相同的條件獨立性和統計獨立性，在傳統的因果推斷算法下判斷異常艱難，但是這三個結構所具有的獨特的外加噪聲項是完全不同的。文獻［18］通過最小二乘回歸計算出殘差，利用外加噪聲項驗證條件獨立性等價于殘差獨立性，使得CI 測試進一步簡化成一個獨立性測試，從而進一步提升速度。

2 預備知識

2.1 因果網絡

因果網絡是一種概率推理模型，由一組隨機變量X=｛x1，x2，…，xn｝以及對應的頂點集V=｛v1，v2，…，vn｝和頂點間邊的集合E=｛e1，e2，…，em｝組成，記作G（X，V，E）。可以看出，因果網絡是由一組隨機變量X與其相對應的一個有向無環圖G（V，E）構成的，其要求X的聯合概率分布P（X）與G（V，E）是對應的，其中，E=｛e（xi，xj）|xi，xj∈X｝表示有向無環圖G中兩個變量xi和xj之間的邊集合，e（xi，xj）表示變量xi和xj之間的關系，可以是xi-xj、xi→xj或xi←xj。需要注意的是，在因果推斷領域中，為了方便或提高閱讀性，在不影響理解的情況下，通常用X同時表示變量與節點，不分開論述。

2.2 馬爾可夫等價

當有兩個有向無環圖是馬爾可夫等價類時，即兩個有向無環圖滿足包含同樣節點、邊（不考慮方向）和V結構的條件時，它們具有完全一樣的獨立性與條件獨立性。以圖1、圖2 為例，這兩個網絡結構是屬于馬爾可夫等價的，他們的獨立性和條件獨立性一致：1）獨立性：X不獨立Y，X不獨立Z和Y不獨立Z；2）條件獨立性：X不獨立Y|Z，X⊥Z|Y，Y不獨立Z|X。因此，對于馬爾可夫等價類，由于它們都具有相同的獨立性與條件獨立性，因此CI 測試無法正確判斷出真實結構。

圖1 順連結構Fig.1 Consequent structure

圖2 分連結構Fig.2 Disjuncition structure

2.3 忠實性假設

已知一個因果網絡G，其中有三個變量（集）X、Y、Z，他們的聯合概率分布為P（XYZ），如果滿足從X⊥Y|Z可以推導出X、Y被Z集合D分離的條件，那么就稱這個聯合概率分布P（XYZ）滿足關于因果網絡G的因果忠實性假設。

3 偏相關性測試

本節將討論如何在因果推斷中應用偏相關性測試檢測條件獨立性是否成立，這里先給出DAUDIN的關于CI 的一個經典理論［19-20］，其中描述了如何在平方可積函數空間中把相關性與獨立性相關聯，將以該理論支撐本文后續的理論分析。

3.1 條件獨立性特性

條件獨立性特性（Characterization of Conditional Independence，CCI）的相關定義如下：

令X、Y、Z為3 個隨機變量或隨機變量集，定義：

CCI 描述的是一種在平方可積函數空間中相關性與獨立性的等價關系。后續將給出一些相關的理論分析與結果，即分別在高斯和非高斯情況下偏相關性與條件（不）獨立性的關系。

3.2 偏相關性與條件獨立性的關系

定理1給定m+2個隨機變量：x、y和Z=｛z1，z2，…，zm｝，它們都是基于獨立的隨機變量si（i=1，2，…，l）的線性組合，如果所有的si都符合聯合高斯分布，那么x⊥y|Z成立，當且僅當σxy.Z=0。

證明考慮到給定變量集Z后x和y的偏相關性系數為ρxy.Z=，其中σ**.Z（給定Z的協方差）可以認為是投影到Z上的x與y殘差的協方差。因此，由獨立與相關關系可得出σxy.Z=cov(x-E(x|Z)，y-E(y|Z))=0。在聯合高斯分布下，σxy.Z=0 ?x⊥y|Z。

另一方面，如果x⊥y|Z，那么根據式（7）進行推斷，因為E(x-E(x|Z)|Z)=0和E(y-E(y|Z)|Z)=0，所以有x-E(x|Z)∈E1和y-E(y|Z)∈E2。此外，考慮到給定Z關于x和y的偏相關性，由于給定Z后的偏協方差σxy.Z等價于x和y在Z形成的線性空間上投影的協方差，因此有：

定理2偏相關性不成立與獨立成立在聯合高斯分布條件下是等價的，即獨立性測試可以由偏相關性測試替換。

推論1給定m+2 個隨機變量：x、y以及Z=｛z1，z2，…，zm｝，它們都是基于獨立的隨機變量si（i=1，2，…，l）的線性組合，如果所有的si符合聯合非高斯分布，那么σxy.Z≠0 ?x不獨立y|Z。

推論1 可以由定理1 直接推導出，表明在非高斯場景下，偏相關性不成立是條件不獨立的充分條件，即可以通過偏相關性測試判定獨立性不成立。所以，在非常苛刻的條件——非高斯且偏相關性成立的情況下，無法判斷條件獨立的情況。因此，在絕大多數情況下，都可以用偏相關性對獨立性進行判斷，從而大幅縮短檢測獨立性是否成立的時間。

4 因果推斷算法

本節提出基于偏相關性測試的遞歸式因果推斷（CIPCT）算法，對算法基本框架、因果分割、因果方向學習與子圖合并4 個部分分別進行介紹。

4.1 算法基本框架

CIPCT 算法是一種遞推式算法，其采用傳統的“分-治-合”方法［21-22］進行因果學習，并在具體的算法步驟中融合偏相關性測試。CIPCT 算法流程如圖3所示，具體步驟如下：

圖3 CIPCT 算法流程Fig.3 Procedure of CIPCT algorithm

1）采用分治策略將龐大的原始數據集進行遞歸，分解成更小的子數據集，使用高效的偏相關性測試提高效率。

2）在分出若干子數據集的基礎上，對每一個子數據集進行因果方向學習，形成子結果。

3）以比較顯著性值的方式作為合并策略，對所有的子結果進行合并整理后，得到原始數據集對應的完整因果網絡圖。

4.2 因果分割

由于傳統的因果推斷算法在面對高維數據時運行時間較長，因此本文在進行因果分割時采用分治策略，并融合偏相關性測試提升效率。因果分割階段的目的就是盡可能地將數據集分割成子數據集，子數據集會通過與父數據集相同或更高階的偏相關性測試繼續進行因果分割，從而減少偏相關性測試的次數，縮短算法在高維環境下的耗時。在基于CI 測試的因果推斷中，一般將CI 測試的執行次數作為算法時間復雜度，因為與CI測試相比，其他操作在執行時間上可以忽略不計。假設原始變量集V=｛v1，v2，…，vn｝被遞歸分割成m個子集｛V1，V2，…，Vm｝，其中，對于所有m，有|Vm|≤n。當假設CIPCT 使用PC 算法進行因果推斷時，則求解子問題的算法時間復雜度為，其中kmax=max（|V1|，|V2|，…，|Vm|）。另一方面，需要在每次迭代中都構造一個CI 測試表，在最壞的情況下，需要針對原始變量集V計算一個σ階CI 測試表（σ是預先設定的參數）。因此，遞歸式推斷算法CIPCT 的時間復雜度為如果只使用PC 算法解決因果推斷問題，則算法復雜度為O（n22n-2）。由于kmax通常要比n小得多，因此CIPCT 算法的算法時間復雜度較傳統基于CI 測試的算法大幅降低。

因果分割的具體實現如算法1 所示，主要步驟如下：

1）構造一個0 階偏相關性測試表M，其通過對數據集V中v1～vn進行偏相關性測試后得到。Mij=1 表示在0 階條件下vi和vj相關系數為0，即corr（vi，vj）=0；Mij=0 表示在0 階條件下vi和vj不相關。

2）從Mij=1 中取vi∈A，vj∈B，把數據集V分為子數據集ABC，其中C中斷了AB之間所有聯系，要求C必須是滿足條件的最小集，但不一定為是A或者B的D 分離集。令D=V，從D中去除掉所有與C不相關的變量，令V1=A∪D，V2=B∪D。

3）若在上一步中無法對V進行因果分割，則構造更高一階的偏相關性測試表，再執行步驟2。通過以上操作，最后可得一個因果分割V=｛V1，V2｝。

算法1因果分割

4.3 因果方向學習

經過因果分割后將會得到若干個子數據集，此處將使用V結構性質進行因果方向學習。首先從局部結構x-z-y入手，如果在給定某個條件集Z之后，能令x、y滿足偏相關性，即pcorr（x，y|Z）=0，根據D 分離準則，此時若z不屬于Z，就可以判斷出x-z-y結構為x→z←y；然后進行一致性傳播［3］，推導出剩下邊的方向。在因果分割過程中會將這幾個相似的局部結構保存，以此來完成下一步子圖合并。

4.4 子圖合并

本文在文獻［21］研究的基礎上，將子圖合并中的CI 測試替換為偏相關測試，根據邊與邊之間的顯著性值sv 的大小決定連接方向。顯著性值sv 表示了方向被接受的概率，給出相關的理論分析如下：

定理3對于任意一個無向圖或局部無向圖，假設給定一個局部結構v1-v2-v3，在對該結構進行基于V結構的方向判斷時，接受兩個方向v1→v2和v3→v2的概率相同，均為svvstructure：=P(v1⊥v3|S)×P(v1⊥v3|(S∪v2))，其中S為滿足CI 的條件集。

證明根據V結構定義，v1與v3在給定某一個條件集S滿足CI，而條件集增加中間節點集后不滿足CI 這兩個特征即可直接推斷得到。首先假設給定一個局部V結構v1-v2-v3，當v1→v2和v3→v2時，它們的顯著性值分別定義為：

其中：S是對應的D 分離集。然后對比顯著性值，取顯著性值更大的方向進行方向連接。最后通過合并所有的子因果網絡，得到關于原始數據集的完整因果網絡結構圖。在子圖合并過程中，顯著性值的計算是在偏相關測試的基礎上進行的，通過檢測上一步因果分割中的偏相關測試得到顯著性值再保存結果，因此，不會增加額外的計算時間。子圖合并的具體實現如算法2 所示。

算法2子圖合并

5 實驗結果與分析

本文所有的實驗均在MATLAB R2016a 中進行，運行環境是一臺CPU 為i5-7200U，主頻為2.50 GHz、內存12 GB 的64 位操作系統的筆記本，將本文提出的CIPCT 算法和目前具有代表性的遞歸式因果推斷算法CAPA［21］進行對比。

實驗所用數據為10 組不同的真實因果網絡結構［22-23］，涵蓋了多個領域，表1 中列出了這些數據的網絡結構信息，其中，節點數代表維度。

表1 10 個因果網絡結構的統計特性Table 1 Statistical characteristics of ten causal network structures

實驗分為兩種方式進行：首先在不同數據場景下固定相同樣本量，將樣本量控制在2|V|=2n，這是一種經典的評估樣本量、節點個數與準確率關系的方法［14］，并同時對比兩種推斷算法的耗時；然后在相同數據環境下，分別比較不同樣本量下的實驗結果，評估樣本量對準確率、耗時的影響。本文實驗使用召回率R、精確率P和F1值F1對算法的綜合性能進行評估。其中：召回率反映算法發現因果關系的查全率；精確率反映算法發現因果關系的查準率；F1 值反映算法的綜合的準確率。這三個指標的計算公式如下：

5.1 不同數據相同樣本量的實驗

在不同數據場景下固定相同樣本量進行對比實驗，選取全部10 個因果網絡結構并將樣本量控制在2|V|=2n，將提出的CIPCT 算法與CAPA 算法［21］進行對比。對比兩種推斷算法的準確率，并列出兩種算法的運行時間進行分析比較，實驗結果如表2 所示，其中加粗數據表示最優值。

表2 CIPIT 和CAPA 算法在不同數據集中的性能對比Table 2 Performance comparison of CIPIT algorithm and CAPA algorithm on different datasets

分析表2 結果可以發現，CIPCT 算法的耗時相對較少，但是精確率相對略低，這是因為偏相關性測試對于高維數據時存在錯誤添加節點的現象。但是由于其采用了計算偏相關性系數的方式，因此在維度增加時可以篩選出更合適的節點，從而加快構造網絡結構的速度，進而加快對高維數據的處理速度。

5.2 相同數據不同樣本量的實驗

設計在相同數據場景下不同樣本量的對比實驗，挑選Barley（兩者準確率相近）和Andes（CAPA 優于CIPIT）分別在不同樣本量（n、2n、3n、5n、7n）的情況下進行對比實驗，實驗結果如表3 所示，其中加粗數據表示最優值。

表3 CIPIT 和CAPA 算法在不同樣本量下的性能對比Table 3 Performance comparison of CIPIT algorithm and CAPA algorithm under different sample sizes

從表3 可以看出，兩種算法在Barley 數據集上準確率幾乎持平，而在Andes 數據集上CAPA 優于CIPIT，因為偏相關性測試是一種近似的CI 測試，在特定的網絡結構下，如Andes，這個網絡結構要求嚴格的CI，則偏相關性測試效果會一定程度下降。此外，可以看到當維度為n時，在Andes 結構中CI 測試的速度較快，原因是CAPA 采用的是啟發式因果分割，在低樣本量情況下，加上CI 測試的耗時也可能要比CIPIT 快。

通過對比表3 中的F1 值和運行時間可以看出，雖然本文算法在Andes 網絡結構中的F1值略低于CAPA 算法，但對比兩種算法的處理時間可知，CIPIT可以在高維網絡結構中進行更快速的有效構建，并且在維度越高的數據環境下這種時間差距越明顯。由于目前其他的因果推斷算法都難以在高維數據環境下進行快速有效的處理，使用偏相關性測試后耗時減少明顯，說明本文算法在運算時間上具有顯著優勢。

6 結束語

本文提出的CIPCT 算法是一種基于偏相關性測試的快速因果推斷算法，其在傳統的“分-治-合”框架中融入偏相關性測試，從而保證準確率同時提高運算效率。在對因果結構進行降維處理“分”的過程中，將整體的因果網絡快速拆分成若干個子因果網絡，對每一個子因果網絡都進行因果推斷，使得因果推斷能夠在高維數據環境下運行，避免傳統CI 測試運算復雜度和時間復雜度高的缺點。實驗結果表明，在高維數據環境下，本文算法與目前具有代表性的遞歸式因果推斷算法CAPA 在準確率幾乎持平的情況下，速度可提升2～10 倍。后續將改進推斷過程及算法框架，將本文算法與神經網絡相結合，進一步提升準確率和降低算法復雜度。