＜はじパタ＞　おまけ：共分散行列のコーディングに苦戦した話

今、はじパタのコーディングをせっせとやってはブログに上げているという作業を続けていますが、いかんせん付け焼刃の知識のために基本的と思われるところで引っかかってしまうことも、多く、今後の備忘録的にも失敗とリカバリの記録を残してみようと思います。

現象：NaNエラーが発生

10.4 「確率モデルによるクラスタリング」をやっていくと、EMモデルの繰り返しによって数値がNaNとなってしまう現象が発生。ひとつNaNがでてくるとその数値を使った値もNaNとなるので真賀田博士もびっくりのすべてがNaNになってしまうのでした。以下がその時のログ。

MU: [[4.2591681450807695, 1.3885253632923946], [3.3363395791482486, 1.0264189358663134], [3.6851362069583895, 1.1830943936011733]]
V:  [[2.8113007747631147, 2.9874538055910858, 3.042873375542265], [0.51266119363098883, 0.56519089148146107, 0.59171106814975671], [1.1252472995735809, 1.2722812742580845, 1.296157580183261]]
N_DET:  -0.860447181357
N_DET:  -0.0603705511147
N_DET:  -0.24839340451
NV: 
nan nan nan
SNV: 
   nan
   nan
   nan
   nan
   nan

原因究明

NaNとなってしまう原因はどこにあるのかを突き止めてみると、ガウス分布の計算内で分散行列の行列式を出すのですが、そこで負の値になってしまっているようでした。その値の平方根を取るものだから当然実数にはならず、NaNとなってしまったのです。なるほど。そもそものこの共分散行列をかくにんしたところ、

{\displaystyle \Sigma } \Sigma は、半正定値行列
分散共分散行列 - Wikipedia

半正定値行列というのは、行列式が０か正の値になる行列のこと。ということはそもそも負の値が出てくること時点でおかしいということです。よし、原因はつかめた。あとは訂正をするだけだ。。。ここからが長かった。

行列式のコーディングのチェック

問題となるのは以下の分散を求める部分のコード。意図的にNumpyの行列を使わず、理解のため分散、共分散をそれぞれ計算してみようということで以下のリスト内包表記のコードとなっています。

        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]

最初のsumと次のsumが2次元の要素それぞれの分散値を計算、最後のsumが共分散を計算しています。なんども確認して問題なさそうだったのですが、その結果が以下のように行列式をとると負の値がでてきてしまう。なぜだ。

#分散1、分散2、共分散、行列式（ad-bc）
     2.81130077476 2.98745380559 3.04287337554 -0.860447181357
     0.512661193631 0.565190891481 0.59171106815 -0.0603705511147
     1.12524729957 1.27228127426 1.29615758018 -0.24839340451

入力値とかも確認してみましたがそもそも共分散行列の行列式が負になること自体が異常なので、やはりこのコードがおかしいんですよね。で、最終的に別の本に書いてあるEMアルゴリズムをコーディングしてみて、その行列要素を出力してみることにしました。すると。。。

V: # Result of ML Book's code 
     2.69263977682 0.640353547689 1.26287084435 0.129398664238
     2.95622008252 0.425917289351 1.08691249439 0.0777264738108
     2.85750222648 0.516449549181 1.19418340666 0.0496817279099

V: # Result of my code
     2.69263977682 2.95622008252 2.85750222648 -0.205283191105
     0.640353547689 0.425917289351 0.516449549181 0.00601751040866
     1.26287084435 1.08691249439 1.19418340666 -0.0534439092193

おおぅ、行列の転置が起きてしまっているようですね。なるほど、リスト内包表記の組み合わせが正しくないようです。基本的な考え方は間違っていないかったですが、こうした順序の認識がまちがっていて、かつ訂正することも難しい（なにが間違っているか理解できない）状況に陥っていました。

訂正

これを踏まえて以下のように訂正をおこなったところ、出力からNaNが出てこなくなりました。

#Corrected
        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)]

#Wrong
        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]