« 読了：Nishisato & Clavel (2003) クロス表の対応分析の結果から「行と列の距離」を求めるには | メイン | 読了：水野(2003) 心理学における「信頼」概念レビュー »
2018年5月14日 (月)
仕事の都合で不意に対応分析を使ったりすると、いつも途中で混乱しちゃうので、思い余ってメモをとった。
でもあれなんだよな、メモを取れば覚えるかというと、これがそうでもないんだよな... 思えば、こういうのはじめて勉強したのは20代前半であった。それから時代は流れ、幕府は倒れ、江戸は東京となり、日本は急速な近代化を遂げたが(←ちょっと記憶が混乱してます)、私は全然進歩してないような気がする。
以下、Greenacre(1984) "Theory and Applications of Correspondence Analysis", Chap.4 の前半のメモ。この本だと、固有値分解ではなくて一般化SVDでスーッと説明しちゃうのである。
準備
$I \times J$行列を$\mathbf{N}=\{n_{ij}\}$とする。値は非負、行和と列和は非ゼロとする。そうでなくてもかまわないけれど、$\mathbf{N}$は二元クロス表だと考えるとわかりやすい。
$\mathbf{N}$で総和で割って$\mathbf{P}=\{p_{ij}\}$とし、これを「対応行列」と呼ぶ。$\mathbf{P}$の行和ベクトルと列和ベクトルをそれぞれ$\mathbf{r}=\{r_{i}\}, \mathbf{c}=\{c_{j}\}$とし、これを持たせた対角行列を$\mathbf{D}_r, \mathbf{D}_c$とする。
行プロファイルと列プロファイル
$\mathbf{P}$のある行$\mathbf{r}_i$をその行の和で割ったもの$\tilde{\mathbf{r}}_i$を「行プロファイル」、ある列$\mathbf{c}_j$をその列の和で割ったもの$\tilde{\mathbf{c}}_j$を「列プロファイル」と呼ぶ。プロファイル行列は
$\mathbf{R} \equiv \mathbf{D}^{-1}_r \mathbf{P}$
$\mathbf{C} \equiv \mathbf{D}^{-1}_c \mathbf{P}^T$
行プロファイル空間と列プロファイル空間
行プロファイル$\tilde{\mathbf{r}}_i$を、$J$次元空間におけるその行の座標とみる。その行の行和、つまり行和ベクトル$\mathbf{r}$の$i$番目の要素$r_i$を「質量」と呼ぶことにする。
二点間の距離を、重み$\mathbf{D}^{-1}_c$による加重ユークリッド距離(=カイ二乗距離)とする。つまり、次元$j$の重みは$1/c_j$となるわけである。
そもそも行$i$の行プロファイルの$j$番目の要素は$p_{ij}/r_i$であった。すべての行プロファイルについてその重心を求めると、$\sum_i r_i (p_{ij}/r_i) / \sum_i r_i = \sum_i p_{ij}= c_j$であるから、つまり行プロファイルの重心は列和ベクトル$\mathbf{c}$となる。
以上は行と列をひっくり返しても成立する。つまり$\tilde{\mathbf{c}}_j$を$I$次元空間における座標とみて、二点間の距離を重み$\mathbf{D}^{-1}_r$による加重ユークリッド距離としたとき、その重心は行和ベクトル$\mathbf{r}$である。
全慣性
ちょっと話が逸れるけど... 上の行プロファイル空間なり列プロファイル空間について、各点の重心からの二乗距離の、質量を重みにした加重和を「全慣性」と呼ぶ。
行$i$の重心からの二乗距離は、もし次元に重みがなかりせば
$(\tilde{\mathbf{r}}_i - \mathbf{c})^T (\tilde{\mathbf{r}}_i - \mathbf{c})$
だが、軸に重み$\mathbf{D}^{-1}_c$がつくので
$(\tilde{\mathbf{r}}_i - \mathbf{c})^T \mathbf{D}^{-1}_c (\tilde{\mathbf{r}}_i - c)$
で、この子の質量が$r_i$である。質量で重みづけて合計すると
$in(I) = \sum_i r_i (\tilde{\mathbf{r}}_i - \mathbf{c})^T \mathbf{D}^{-1}_c (\tilde{\mathbf{r}}_i - \mathbf{c})$
同様に
$in(J) = \sum_j c_j (\tilde{\mathbf{c}}_j - \mathbf{r})^T \mathbf{D}^{-1}_r (\tilde{\mathbf{c}}_j - \mathbf{r})$
これが、行プロファイルの空間の全慣性$in(I)$, 列プロファイルの空間の全慣性$in(J)$である。
導出は省略するけど、結局
$in(I)= in(J) = \sum_i \sum_j \{(p_{ij} - r_i c_j)^2/r_i c_j\}$
であることが示せる。この式をよーく見ると、$\mathbf{N}$を二元クロス表とみて、その独立性検定をしたときのカイ二乗値を、総和で割った値になっている。
一般化SVDの登場
さあ、ここからが本番だ！
上の行プロファイル空間・列プロファイル空間から、少数次元の下位空間を取り出したい。その際、オリジナルの空間からの各点の二乗距離の、質量を重みにした加重和を最小化したい。
こういうのは、ご存じ一般化SVD
$\mathbf{A} = \mathbf{N} \mathbf{D}_\mu \mathbf{M}^T$ ただし $\mathbf{N}^T \mathbf{\Omega} \mathbf{N} = \mathbf{M}^T \mathbf{\Phi} \mathbf{M} = \mathbf{I}$
の得意技である。
主はいわれた。あなたがたはこれまで、$\mathbf{A}$の行を点、列を次元とみなしていた。はっきり言っておくが、$\mathbf{N}$の左$K^*$列, $\mathbf{D}_\mu$の左上$K^* \times K^*$, $\mathbf{M}$の左$K^*$列を取ってきて
$\mathbf{A}_{[K^*]} \equiv \mathbf{N}_{[K^*]} \mathbf{D}_{\mu [K^*]} \mathbf{M}^T_{[K^*]}$
とし、$\mathbf{N}_{[K^*]} \mathbf{D}_{\mu [K^*]}$の行を点、列を次元とみなしなさい。なぜなら、ランク$K^*$以下のあらゆる行列$\mathbf{X}$のなかで、二乗距離
$trace\{\mathbf{\Omega}(\mathbf{A}-\mathbf{X}) \mathbf{\Phi}(\mathbf{A}-\mathbf{X})^T\}$
を最も小さくする$\mathbf{X}$こそ、$\mathbf{A}_{[K^*]}$だからである。
行プロファイル行列と列プロファイル行列の一般化SVD
行プロファイル行列$\mathbf{R} \equiv \mathbf{D}^{-1}_r \mathbf{P}$について考えよう。あらかじめ$\mathbf{R}$の各列$j$からその列の重心$c_j$を引き($\mathbf{R} - \mathbf{1} \mathbf{c}^T $)、各列の重心を0にしておく。
おさらいしておくと、点$i$の質量は$r_i$で、次元$j$の重みは$1/c_j$であった。ってことは、主よ、行$i$に重み$r_i$, 列$j$に重み$1/c_j$をつけた一般化SVDで
$\mathbf{D}^{-1}_r \mathbf{P} - 1\mathbf{c}^T = \mathbf{L} \mathbf{D}_\phi \mathbf{M}^T $ ただし $\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{M}^T \mathbf{D}_c^{-1} \mathbf{M} = \mathbf{I}$ [式1]
と分解し、$\mathbf{F} \equiv \mathbf{L} \mathbf{D}_\phi$を座標にすればよろしいのですね。
同様に、列プロファイル行列$\mathbf{C} \equiv \mathbf{D}^{-1}_c \mathbf{P}$については、
$\mathbf{D}^{-1}_c \mathbf{P}^T - 1\mathbf{r}^T = \mathbf{Y} \mathbf{D}_\psi \mathbf{Z}^T $ ただし $\mathbf{Y}^T \mathbf{D}_c \mathbf{Y} = \mathbf{Z}^T \mathbf{D}_r^{-1} \mathbf{Z} = \mathbf{I}$ [式2]
と分解し、$\mathbf{G} \equiv \mathbf{Y} \mathbf{D}_\psi$を座標にすればいいわけだ。
式1と式2は結局同じ行列の一般化SVDである (その1)
さて、式1に左から$\mathbf{D}_r$を掛ける。
$\mathbf{P} - \mathbf{rc}^T = (\mathbf{D}_r \mathbf{L}) \mathbf{D}_\phi \mathbf{M}^T$
左辺の分解は、もはや重み$\mathbf{D}_r, \mathbf{D}^{-1}_c$の一般化SVDではない。ええ、はい、わかってます、確かにそうなんです。しかしですね、$\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{I}$の$\mathbf{L}$に左から$\mathbf{D}_r$を掛けまくるという極悪非道な真似をしますとですね、
$(\mathbf{D}_r \mathbf{L})^T \mathbf{D}_r^{-1} (\mathbf{D}_r \mathbf{L}) = \mathbf{I}$
である。$\mathbf{D}_r$を二回も掛けた分、真ん中を逆数にして取り返していることに注意。また、
$\mathbf{M}^T \mathbf{D}^{-1}_c \mathbf{M} = \mathbf{I}$
であることには変わりない。結局、重み$\mathbf{D}^{-1}_r, \mathbf{D}^{-1}_c$の一般化SVDになっているわけだ。[←この理屈に気が付かなかったせいで、かなり時間を食いました...]
同様に、式2に左から$\mathbf{D}_c$を掛ける。
$\mathbf{P}^T - \mathbf{cr}^T = (\mathbf{D}_c \mathbf{Y}) \mathbf{D}_\psi \mathbf{Z}^T$
$\mathbf{P} - \mathbf{rc}^T = \mathbf{Z} \mathbf{D}_\psi (\mathbf{D}_c \mathbf{Y})^T$
ここで
$\mathbf{Z}^T \mathbf{D}_r^{-1} \mathbf{Z} = \mathbf{I}$
$(\mathbf{D}_c \mathbf{Y})^T \mathbf{D}^{-1}_c (\mathbf{D}_c \mathbf{Y}) = \mathbf{I}$
こちらも重み$\mathbf{D}^{-1}_r, \mathbf{D}^{-1}_c$の一般化SVDになっている。
あらやだわ、二人とも同じ行列$\mathbf{P} - \mathbf{rc}^T$を、同じ重み$\mathbf{D}_r^{-1}, \mathbf{D}_c^{-1}$で一般化SVDしてたのね。おほほほほ。
というわけで、一般化SVDの解が一意に決まるとして、
$\mathbf{D}_\phi = \mathbf{D}_\psi=\mathbf{D}_\mu$
$\mathbf{D}_r \mathbf{L} = \mathbf{Z} = \mathbf{A}$
$\mathbf{M} = \mathbf{D}_c \mathbf{Y} = \mathbf{B}$
と置き、
$\mathbf{P} - \mathbf{rc}^T = \mathbf{A} \mathbf{D}_\mu \mathbf{B}^T$ ただし $\mathbf{A}^T \mathbf{D}_r^{-1} \mathbf{A} = \mathbf{B}^T \mathbf{D}_c^{-1} \mathbf{B} = \mathbf{I}$
厳密にいうと、一般化SVDの解が一意に決まらない場合(特異値のなかに同じ値がある場合)もありうるけど、要するに上のように書ける。
座標は次のようになる。
行プロファイルの座標: $\mathbf{F} = \mathbf{L} \mathbf{D}_\phi = \mathbf{D}_r^{-1} \mathbf{A} \mathbf{D}_\mu$
列プロファイルの座標: $\mathbf{G} = \mathbf{Y} \mathbf{D}_\psi = \mathbf{D}_c^{-1} \mathbf{B} \mathbf{D}_\mu$
式1と式2は結局同じ行列の一般化SVDである (その2)
ここでChap.4 からちょっと離れるけど...
式1に右から$\mathbf{D}^{-1}_c$を掛けたら何が起きるか。
$\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T = \mathbf{L} \mathbf{D}_\phi (\mathbf{M}^T \mathbf{D}^{-1}_c) = \mathbf{L} \mathbf{D}_\phi (\mathbf{D}^{-1}_c \mathbf{M})^T $
ここで
$\mathbf{L}^T \mathbf{D}_r \mathbf{L} = \mathbf{I}$
$(\mathbf{D}^{-1}_c \mathbf{M})^T \mathbf{D}_c (\mathbf{D}^{-1}_c \mathbf{M}) = \mathbf{I}$
というわけで、こんどは重み$\mathbf{D}_r, \mathbf{D}_c$の一般化SVDになる。
同様に、式2に右から$\mathbf{D}^{-1}_r$を掛ける。
$\mathbf{D}^{-1}_c \mathbf{P}^T \mathbf{D}^{-1}_r - 11^T = \mathbf{Y} \mathbf{D}_\psi (\mathbf{Z}^T \mathbf{D}^{-1}_r) = \mathbf{Y} \mathbf{D}_\psi (\mathbf{D}^{-1}_r \mathbf{Z})^T$
$\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T = (\mathbf{D}^{-1}_r \mathbf{Z}) \mathbf{D}_\psi \mathbf{Y}^T$
ここで
$(\mathbf{D}^{-1}_r \mathbf{Z})^T \mathbf{D}_r (\mathbf{D}^{-1}_r \mathbf{Z}) = \mathbf{I}$
$\mathbf{Y}^T \mathbf{D}_c \mathbf{Y} = \mathbf{I}$
というわけで、こんどは重み$\mathbf{D}_r, \mathbf{D}_c$の一般化SVDになる。くどいですね。
あらやだわ二人とも同じ行列$\mathbf{D}^{-1}_r \mathbf{P} \mathbf{D}^{-1}_c - 11^T$を同じ重み$\mathbf{D}_r, \mathbf{D}_c$で一般化SVDしてたのねおほほほほ、というわけで、
$\mathbf{D}_\phi = \mathbf{D}_\psi=\mathbf{D}_\mu$
$\mathbf{L} = \mathbf{D}^{-1}_r \mathbf{Z} = \mathbf{A'}$
$\mathbf{D}^{-1}_c \mathbf{M} = \mathbf{Y} = \mathbf{B'}$
と置き、
$\mathbf{D}^{-1}_c \mathbf{P} \mathbf{D}^{-1}_r - 11^T = \mathbf{A'} \mathbf{D}_\mu \mathbf{B'}^T$ ただし $\mathbf{A'}^T \mathbf{D}_r \mathbf{A'} = \mathbf{B'}^T \mathbf{D}_c \mathbf{B'} = \mathbf{I}$
という一般化SVDで表現できる。
座標は次のようになる。
行プロファイルの座標: $\mathbf{F} = \mathbf{L} \mathbf{D}_\phi = \mathbf{A'} \mathbf{D}_\mu$
列プロファイルの座標: $\mathbf{G} = \mathbf{Y} \mathbf{D}_\psi = \mathbf{B'} \mathbf{D}_\mu$
つまり、一般化SVDで得られた左特異行列と右特異行列に特異値を掛けるだけで、座標が得られるわけである。その点ではこの表現のほうがスマートである。実際、この本のAppendix Aではこういう表現になっている。
Chap.4 に戻り、最後に２つだけメモしておこう。
遷移方程式
導出は省略するけど、
$\mathbf{G} = \mathbf{C} \mathbf{F} \mathbf{D}^{-1}_\mu$
$\mathbf{F} = \mathbf{R} \mathbf{G} \mathbf{D}^{-1}_\mu$
という関係が成り立つ。つまり、列$j$の列プロファイル(=$\mathbf{G}$の$j$行目)は、行プロファイル行列$\mathbf{F}$の各行の重心$\tilde{\mathbf{c}}_j F$を、次元$k$について$1/\mu_k$したものなのだ。うわーすごいわねー[←出来レース]。これを遷移方程式という。
中心化せずに一般化SVDしたら？
さきほどは$\mathbf{P}-\mathbf{rc}^T$を一般化SVDしたが、うっかり$\mathbf{P}$を一般化SVDしちゃったらどうなるか。
$\mathbf{P} = \tilde{\mathbf{A}} \mathbf{D}_\tilde{\mu} \tilde{\mathbf{B}}^T$ ただし$\tilde{\mathbf{A}}^T \mathbf{D}_r^{-1} \tilde{\mathbf{A}} = \tilde{\mathbf{B}}^T \mathbf{D}_c^{-1} \tilde{\mathbf{B}} = I$
めんどくさいので証明は省略するけど、$\tilde{\mathbf{A}}$は「$\mathbf{A}$の左にベクトル$\mathbf{r}$をくっつけた行列」、$\tilde{\mathbf{B}}$は「$\mathbf{B}$の左にベクトル$\mathbf{r}$をくっつけた行列」となり、第一特異値は1となり、元の特異値は一個づつ後ろにずれるのである。ああ、なんかわかるような気がしてきた。中心化してないので「質量因子」みたいのが出てきちゃうわけだ。
雑記：データ解析 - 覚え書き：対応分析を一般化特異値分解で説明する