LiNGAMモデルの識別可能性について

統計的因果探索の重要なモデルであるLiNGAMモデルの識別可能性について簡単にまとめます。
主に機械学習プロフェッショナルシリーズ「統計的因果探索」のP95-105をまとめたものになります。
(当方は統計的因果探索初学者ですので、間違った解釈してるかもです。その時はご指摘いただけると幸いです。)

## 統計的因果探索の数理モデル

$n$個の観測変数 $ \boldsymbol{x} = ${ $x_{1}, x_{2}, \cdots , x_{n}$ } があったとき、どの変数がどの変数にどれくらい影響を与えているか（因果関係）を推定するのが統計的因果探索かと思います。図で描くと以下のような感じになります。（このような図を「因果グラフ」といいます。）
![image.png](https://qiita-image-store.s3.amazonaws.com/0/240999/adc30a9e-ddbd-09c4-ea5a-7ac0394f8ac5.png)

これを数理的に表現すると以下のようになります。

```math
\boldsymbol{x} = f(\boldsymbol{x}, \boldsymbol{e})
```
ただし、$\boldsymbol{e}=${$e_{1}, e_{2},\cdots,e_{n}$} は未観測変数（外生変数）といって、観測変数に影響を与える観測されなかった変数を表します。ノイズと考えて良さそうです。$f$ は各変数の関係性を表しています。このような方程式を構造方程式と呼びます。

## LiNGAMモデルとは

LiNGAM（linear non-Gaussian acyclic model）モデルは上記に以下の仮定を与えたモデルです。
$f$ : 線型
$e_{1},e_{2},\cdots,e_{n}$　は互いに独立で非ガウス連続分布に従う
各 $e_{i}$ に独立性を仮定することは、各観測変数に対して、未観測共通原因がないことを意味します。（$x_{i}$ に影響を与える未観測変数 $e_{i}$ は他の観測変数 $x_{j}, j \neq i$　に影響を与えない。上の図で言うと、$e_{1}, e_{2}, e_{3}$ を結ぶ矢印がない。）
さらに各観測変数の因果関係に非巡回性を仮定します。これは因果グラフを考えた時、どの変数から出発して関係性をたどっても、元の変数に戻ってこれないことを意味します。例えば因果グラフは以下のような感じです。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/240999/32436a53-2f37-7b29-0b92-a32772d87b3f.png)


$f$ に線型性を仮定したので、上記のモデルを行列を用いて以下のように記述することができます。

```math
\boldsymbol{x} = B\boldsymbol{x} + \boldsymbol{e}
```
ただし、 $B$ は $n\times n$-行列です。

この $B$ を推定するのが統計的因果探索の目的になるかと思います。例えば

```math
B = \left (
\begin{matrix}
0 & 0 & 0 \\
b_{21} & 0 & 0 \\
b_{31} &  b_{32} & 0
\end{matrix}
\right )
```
であった場合、各観測変数は、

```math
\begin{align}
x_{1} &= e_{1} \\
x_{2} &= b_{21}x_{1} + e_{2} \\
x_{3} &= b_{31}x_{1} + b_{32}x_{2} + e_{3}
\end{align}
```
という関係性になることがわかります。上記式からわかるように行列のゼロ・非ゼロパターンが因果グラフの矢印の向きを表し、非ゼロの成分は各観測変数が左辺の観測変数にどれくらい影響を与えるかの尺度を表します。


## LiNGAMモデルの識別可能性
LiNGAMモデル $\boldsymbol{x} = B\boldsymbol{x} + \boldsymbol{e}$ の行列 $B$ を推定できれば、上記の例のように各観測変数の因果関係を特定することができます。LiNGAMモデルはこの行列 $B$ を識別することができます。

### 解説
まず、非巡回性の仮定から行列 $B$ を厳密な下三角行列と仮定しても一般性を失わないことが知られています。(K.A.Bollen, Structural equations with latent variables. 1989)（厳密な下三角行列とは、対角成分が全て $0$ であるような下三角行列です。）　なので、以下では行列 $B$ は厳密な下三角行列と仮定します。LiNGAMモデル $\boldsymbol{x} = B\boldsymbol{x} + \boldsymbol{e}$ の右辺の $B\boldsymbol{x}$ を左辺に移行して、$$(I-B)\boldsymbol{x} = \boldsymbol{e}$$とします。行列 $I-B$ の逆行列を両辺左から掛けると、

```math
\begin{align}
(I-B)^{-1}(I-B) \boldsymbol{x} &= (I-B)^{-1}\boldsymbol{e} \\
\boldsymbol{x} &= (I-B)^{-1} \boldsymbol{e} \\
\boldsymbol{x} &= A\boldsymbol{e}　　　（A:= (I-B)^{-1}）
\end{align}
```
のように変形できます。ここで、LiNGAMモデルの誤差変数ベクトル $\boldsymbol{e}$ の非ガウスかつ独立性の仮定から、独立成分分析（independent component analysis, ICA）の考え方によると、行列 $W := A^{-1} (= I-B)$ が行の順序と尺度を除いて識別可能であることが知られています。つまり、ある置換行列 $P$、対角行列 $D$ が存在して、以下の行列 $W_{ICA}$ が識別可能です。

```math
W_{ICA} = PDW\ (=PD(I-B))
```

よって $P$ と $D$ が推定できれば、これらの逆行列を左から掛けることで、

```math
\begin{align}
D^{-1}P^{-1}W_{ICA} &= D^{-1}P^{-1}PDW \\
 &= W \\
 &= I-B
\end{align}
```
となり、 $B = I - D^{-1}P^{-1}W_{ICA}$ と求めることができます。

#### P の求め方

非巡回性から行列 $B$ は厳密な下三角行列と仮定してよいので、行列 $W = I-B$ は、対角成分が全て $1$ である下三角行列です。なので、行列 $DW$ は対角成分が全て $0$ でない下三角行列であることがわかります。よって、$W_{ICA}$ は対角成分が $0$ でない下三角行列になにかしら置換が施された行列であることがわかります。

推定済みの行列 $W_{ICA}$ の対角成分に $0$ が現れなくなるような置換行列 $\bar{P}$ を探します。すると、 $\bar{P}W_{ICA} = \bar{P}PDW$ となるわけですが、LiNGAM本家の論文
[S.Shimizu, P.O.Hoyer, A.Hyvarinen and A.Kerminen. A linear non-Gaussian acyclic model for causal discovery (2006)](https://www.cs.helsinki.fi/u/ahyvarin/papers/JMLR06.pdf)
のAppendix Lemma1により対角成分に $0$ がない下三角行列に単位行列でない置換行列（どこかしら順番を入れ替える置換行列）を掛けると、必ず対角成分のどこかに $0$ がくることが知られています。
つまり、 $\bar{P}P \neq I$ だったとすると、対角成分が $0$ でない下三角行列 $DW$ に何かしら置換が施されたことになり、$\bar{P}W_{ICA}$ の対角成分が全て $0$ でないことに矛盾します。つまり $\bar{P}P = I （\bar{P} = P^{-1}）$ であることがわかります。
よって、求めたい $P$ は推定済み行列 $W_{ICA}$ の対角成分に $0$ が現れなくなるような置換行列 $\bar{P}$ の逆行列として求まりました。

#### D の求め方
行列 $W$ の対角成分が全て $1$ であることを思い出すと、 $diag(D) = diag(DW)$ であることがわかります。( $diag(\cdot)$ は対角成分だけ抜き出した行列を意味します。)
つまり、

```math
D = diag(D) = diag(DW) = diag(\bar{P}W_{ICA})
```
となるので、 $D$ が $diag(\bar{P}W_{ICA})$ として求まりました。

## まとめ
LiNGAMモデルの誤差変数の非ガウス独立性から独立成分分析の結果が使えて、非巡回性から係数行列が厳密な下三角行列と見なすことができる、ということが重要だと思いました。LiNGAMモデルの美しさがちょっとだけ体感できた気がします。
今度は実際にプログラミングで実データを扱って見たいと思います。

## 参考文献
どう考えても本家本元の清水先生のスライドを見たほうがわかりやすいです。
[構造方程式モデルによる因果推論: 因果構造探索に関する最近の発展](https://www.slideshare.net/sshimizu2006/bsj2012-tutorial-finalweb)
[因果探索: 基本から最近の発展までを概説](https://www.slideshare.net/sshimizu2006/ss-60015497)
[統計的因果推論への招待 -因果構造探索を中心に-](https://www.slideshare.net/sshimizu2006/ss-77876940)