@2007__JMLR__Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm

# Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm > [!abstract] 概要 > 本論文は、ガウス分布に対応する非常に高次元の有向非巡回グラフ(DAG)のスケルトンと同値類を推定するための PC アルゴリズム(Spirtes et al., 2000)を扱う。PC アルゴリズムは、ノード(変数)が多い疎な問題に対して計算上実行可能であり、しばしば非常に高速である。また、真の DAG の疎さの関数として高い計算効率を自動的に達成するという好ましい性質を持つ。ノード数が標本数 n とともに、任意の 0 < a < ∞ について O(n^a) と同程度に急速に増加してもよい、非常に高次元で疎な DAG に対して、このアルゴリズムの一様一致性を証明する。疎さの仮定はかなり最小限であり、DAG における近傍が標本数 n より低いオーダーであることだけを要求する。シミュレーションデータに対する PC アルゴリズムの実証も行う。 ## 論文情報 - 著者: [[Markus Kalisch]]、[[Peter Bühlmann]](ともに [[ETH Zürich]] Seminar für Statistik) - 掲載: Journal of Machine Learning Research 8 (2007) 613-636。編集は David Maxwell Chickering。24 ページ。 - 実装: R パッケージ pcalg(著者らが公開)。 ## 概要 [[PCアルゴリズム]]を高次元の統計的設定(変数数 p が標本数 n より大きくなりうる)へ拡張して解析した論文である。完全無向グラフから条件付き独立性の判定で辺を再帰的に削除してスケルトンを得て、そこから有向化ルールで CPDAG(同値類の代表)へ進む。論文の貢献は、疎性の仮定のもとで一様一致性を証明した点と、計算量が疎さに応じて自動的に下がることを示した点にある。著者らは、高次元 DAG に対して計算可能かつ証明可能に正しい方法は、この時点で PC アルゴリズムだけだと位置づける。 ## 問題設定 - DAG の推定は、DAG の数がノード数に対して超指数的に増えるため難しい。探索とスコアに基づく既存手法(木に制限する MWST、貪欲探索、同値類上を探索する GES)は、ノード数が小〜中程度の問題に向き、ベイズ的手法は計算負荷が高い。 - スケルトン(DAG の向きを無視した無向グラフ)は、一般には条件付き独立グラフ(CIG)と一致せず、CIG の推定法をそのまま借りられない。スケルトンの辺は、他の変数を考慮しても説明し切れない強い依存を表す。 - 忠実性だけを仮定すると一様一致性は得られず、各点での一致性しか得られない(Spirtes et al., 2000; Robins et al., 2003)。本論文はこれを、追加の仮定で一様一致性に拡張し、p・近傍サイズ・偏相関の下限が n に依存して変化する設定へ広げる。 - 同値類が求まればマルコフブランケットが読み取れるため、特徴選択とも重なる。 ## 提案手法新しいアルゴリズムの提案ではなく、既存の PC アルゴリズムの高次元での定式化と解析である。 ### 準備 - DAG の分布 P が **忠実**であるとは、条件付き独立性が d 分離と一対一に対応することをいう。ガウス族では非忠実な分布はルベーグ零集合になる。 - 2 つの DAG が同値であることと、スケルトンと v 構造が同じであることは同値である(Verma と Pearl)。同値類は CPDAG で一意に表せる。 - CPDAG の推定は、(1) スケルトンの推定と (2) 辺の部分的な向き付けの 2 段階に分かれる。統計的推論はすべて (1) で行い、(2) は (1) の結果への決定論的な規則の適用にすぎない。したがって (1) が正しければ (2) は失敗しない。 ### スケルトンの推定(アルゴリズム 1) 完全無向グラフから始め、条件付けの集合の大きさ ℓ = 0, 1, 2, ... の順に、隣接する各ノード対 (i, j) について、i の隣接集合から j を除いた集合の大きさ ℓ の部分集合 k を選び、i と j が k で条件付き独立なら辺を削除して分離集合 S(i, j) に k を記録する。ℓ が、どのノードの隣接集合の大きさよりも大きくなったら止める。到達した最大の ℓ を m_reach と書く。母集団版では、真のスケルトンが得られ、m_reach ∈ {q − 1, q}(q は最大近傍サイズ)となる(命題 1)。標本版では、ガウス分布のもとで偏相関 ρ_{i,j|k} = 0 が条件付き独立と同値であることを使い、Fisher の z 変換 $Z(i,j|k)=\tfrac12\log\frac{1+\hat\rho_{i,j|k}}{1-\hat\rho_{i,j|k}}$ に対し、√(n − |k| − 3)·|Z| ≤ Φ⁻¹(1 − α/2) なら独立と判定する。偏相関は再帰式で計算する。調整パラメータは有意水準 α だけである。 ### スケルトンから CPDAG へ(アルゴリズム 2) 分離集合 S を使って、隣接しないノード対 (i, j) の共通近傍 k が S(i, j) に含まれなければ i − k − j を i → k ← j(v 構造)に置き換える。その後、規則 R1〜R4 を繰り返し適用してできるだけ多くの無向辺に向きを与える。出力が CPDAG であることは Meek(1995)が示している。 ### 一致性の仮定と定理 - (A1) P_n は多変量ガウスで、G_n に忠実である。 - (A2) p_n = O(n^a)(0 ≤ a < ∞)。 - (A3) 最大近傍サイズ q_n = O(n^{1−b})(0 < b ≤ 1)。これが疎性の仮定である。 - (A4) 非零の偏相関の絶対値は c_n 以上(c_n⁻¹ = O(n^d)、0 < d < b/2)で、上限は M < 1 である。定理 1 は、α_n = 2(1 − Φ(n^{1/2} c_n / 2)) と選べば、推定スケルトンが真のスケルトンと一致する確率が 1 − O(exp(−C n^{1−2d})) で 1 に収束することを述べる。定理 2 は CPDAG についても同じ収束を述べる。証明の要点は、標本偏相関が一様に一致すること(Hotelling の結果を用いた指数型の裾確率の評価)と、検定の回数が高々 O(p^{m_reach}) であることによる。α_n は未知の c_n に依存するため、実用上は構成的でない。 ## 新規性 - 忠実性のみでは各点一致性しか得られないという既知の限界に対し、疎性と偏相関の下限という条件を加えて一様一致性を示した。 - p と近傍サイズが標本数とともに増え、偏相関の下限が n とともに減る漸近設定で一致性を示した。著者らは、高次元で計算可能かつ証明可能に正しい DAG 推定法として初めての解析だと述べる。 - 高次元では CPDAG より局所的な解釈ができるスケルトンを、より現実的な推定目標として推奨する。 ## 実験設定 - ノードの順序を固定し、隣接行列の下三角の各成分を成功確率 s のベルヌーイ乱数で 0/1 にし、1 の成分を Uniform[0.1, 1] の重みに置き換えて DAG を生成する。期待近傍サイズは E[N] = s(p − 1) である。データは X(i) = Σ_k A_ik X(k) + ε(i)、ε ~ N(0, 1) で生成する。 - 指標は真陽性率(TPR)、偽陽性率(FPR)、構造ハミング距離(SHD。推定 CPDAG を正しい CPDAG に変えるための辺の挿入・削除・反転の回数)。 - α の選択実験: α は 8 通り(0.00005〜0.1)、p ∈ {7, 15, 40, 70, 100}、n ∈ {30〜30000} の 7 通り、E[N] ∈ {2, 5} で、各 40 反復。 - 高次元実験: p を 9 から 2187 まで 3 倍ずつ、n を 50 から 300 まで線形、E[N] = 0.2√n で増やし、α = 0.05、20 反復。 - 計算時間: n = 1000 固定、p = 10〜1000、E[N] ∈ {2, 8}、α = 0.01、各 10 反復。Athlon 64 X2 2.6 GHz、4 GB、R 2.4.1。 ## 実験結果 ### 有意水準の選択 ![[wiki/sources/_attachments/2007__Estimating-High-Dimensional-Directed-Acyclic-Graphs-with-the-PC-Algorithm/fig01-alpha-vs-ave-shd.png]] *Figure 1: α と平均 SHD(95% 信頼区間つき)。平均 SHD は α = 0.005〜0.01 付近で最小になる。* 70 通りの設定で平均した SHD は α が 0.005 と 0.01 で最小になり、大きな α への悪化のほうが顕著である。Wilcoxon 検定と Bonferroni 補正では、この 2 つは他の値より有意に小さい。 ### パラメータ別の性能 ![[wiki/sources/_attachments/2007__Estimating-High-Dimensional-Directed-Acyclic-Graphs-with-the-PC-Algorithm/fig02-performance-by-parameter.png]] *Figure 2: p = 7, 40, 100 での TPR・FPR・SHD の平均(三角は E[N] = 5、丸は E[N] = 2)。* α = 0.01 では、密なグラフ(E[N] = 5)のほうが疎なグラフ(E[N] = 2)より当てはまりが悪い。TPR と SHD は n の増加で明確に改善するが、FPR の挙動は明確でない。全 n で同じ α を使ったためだと説明されている。 ### 高次元での挙動 | p | n | E[N] | TPR | FPR | |---|---|---|---|---| | 9 | 50 | 1.4 | 0.61 (0.03) | 0.023 (0.005) | | 27 | 100 | 2.0 | 0.70 (0.02) | 0.011 (0.001) | | 81 | 150 | 2.4 | 0.753 (0.007) | 0.0065 (0.0003) | | 243 | 200 | 2.8 | 0.774 (0.004) | 0.0040 (0.0001) | | 729 | 250 | 3.2 | 0.794 (0.004) | 0.0022 (0.00004) | | 2187 | 300 | 3.5 | 0.805 (0.002) | 0.0012 (0.00002) | *Table 1: p を指数的、n を線形、E[N] を劣線形に増やす設定(α = 0.05、20 反復、括弧内は標準偏差)。* ![[wiki/sources/_attachments/2007__Estimating-High-Dimensional-Directed-Acyclic-Graphs-with-the-PC-Algorithm/fig03-high-dim-tpr-fpr.png]] *Figure 3: 表1の設定での TPR と FPR の箱ひげ図。TPR は上がり FPR は下がる。* p が急増しても TPR は増え FPR は減る。これは理論と整合する。一方、近傍サイズは n とともに増えるため、辺の総数に占める真の辺の割合は下がる。疎さの意味(近傍サイズか辺の割合か)によって見え方が変わる。 ### 計算量計算量は最悪で O(p^{m̂_reach})、高い確率で O(p^q) に抑えられる。この上界は多くの分布でかなり緩く、最大近傍サイズ 30 のノードを含む比較的密なグラフでも計算できた。 | p | E[N] | スケルトン(秒) | CPDAG(秒) | |---|---|---|---| | 10 | 2 | 0.037 (0.004) | 0.072 (0.005) | | 10 | 8 | 0.093 (0.005) | 0.124 (0.006) | | 30 | 2 | 0.15 (0.02) | 0.23 (0.02) | | 30 | 8 | 0.84 (0.05) | 0.93 (0.05) | | 50 | 2 | 0.33 (0.01) | 0.48 (0.02) | | 50 | 8 | 2.2 (0.06) | 2.4 (0.06) | | 100 | 2 | 1.03 (0.05) | 1.49 (0.05) | | 100 | 8 | 8.9 (0.3) | 9.4 (0.27) | | 300 | 2 | 8.3 (0.1) | 13.8 (0.13) | | 300 | 8 | 89 (3) | 95 (3) | | 1000 | 2 | 116 (0.5) | 262 (0.8) | | 1000 | 8 | 1300 (60) | 1445 (59) | *Table 2: 平均プロセッサ時間(n = 1000、α = 0.01、括弧内は標準誤差)。* ![[wiki/sources/_attachments/2007__Estimating-High-Dimensional-Directed-Acyclic-Graphs-with-the-PC-Algorithm/fig04-processor-time.png]] *Figure 4: p に対する平均プロセッサ時間(両対数。三角は E[N] = 8、丸は E[N] = 2)。* p = 100 までは約 1 秒(疎な場合)で推定でき、p = 1000 かつ E[N] = 8 でも約 25 分である。両対数の傾きはおおよそ 2(二乗オーダー)である。疎な場合に曲線が上に凸なのは、p とともに最大近傍サイズが(制御されずに)増えるためだと推測されている。スケルトンから CPDAG への追加時間は、スケルトンの数%〜ほぼ 100% で、p が大きいほど割合は下がる。 ### R パッケージ pcalg による例 ![[wiki/sources/_attachments/2007__Estimating-High-Dimensional-Directed-Acyclic-Graphs-with-the-PC-Algorithm/fig05-pcalg-example.png]] *Figure 5: (a) 真の DAG、(b) 推定スケルトン(α = 0.05、n = 10000)、(c) 推定 CPDAG。辺の太さは z 値(信頼度)を表す。* p = 10、n = 10000 で pcAlgo によりスケルトンを、udag2cpdag により CPDAG を推定する。CPDAG の両矢印は向きを決められない無向辺である。 ## 考察 - 疎さ(最大近傍サイズ)が、統計的一致性と計算可能性の両方を決める鍵になる。 - 推定の全推論はスケルトンの段階に集約される。高次元ではスケルトンのほうが CPDAG より回復しやすく、解釈も局所的で信頼しやすいため、より単純だが現実的な目標として使える。 - α の理論値は未知の下限 c_n に依存するため、実用では n が大きいほど小さくする調整が必要になる(実験では 0.005〜0.01 が良好)。 - Lasso による無向グラフの推定は同様の計算効率をもつが、次元が固定でも非一致になりうる(Zhao と Yu)と述べられている。 ## 強み / 弱点・課題 - 強み: 高次元・疎という現実的な設定で一様一致性を示した。疎さに応じて計算量が自動的に下がる。調整パラメータが α の 1 つで、公開実装がある。 - 弱点・課題: ガウス性と忠実性(A1)が必要で、実データでは成り立たない場合がある。一致性の速度は緩い評価で、定数が明示されない。α_n の構成的な選び方が与えられない。真の DAG が疎でない場合(近傍が n と同程度)の保証はない。図の実験は人工データだけで、実データ上の検証はない。検定の誤りが v 構造の向きに伝播しやすい点は、CPDAG の頑健性が低い理由として述べられるのみで定量化されない。