PRML4章の解説と実装


PRML学習記

この度「パターン認識と機械学習」4章の輪講発表担当になったので、勉強したことやちょっとした解説などを書いていきたいと思う。自分もこの本に苦戦した人の一人なので、今後似たような境遇の人がいたときに参考になればとても嬉しい。もし数理的な誤り等を見つけたり、もっとこうした方がいいといった指摘があれば遠慮なくして頂けると助かります。

フィッシャーの線形判別

2クラス

識別関数の項は最小二乗から始まっているがそもそも最小二乗は「うまく使えないのは当たり前」という結論なので割愛。ということで2クラスのフィッシャーから。ここでは線形識別を次元削減の視点から見る。

入力としてD次元ベクトルを得て、以下の式で1次元に射影

y = \boldsymbol{w}^T\boldsymbol{x}

$y$に閾値を設定して$y \ge -w_0 $ のときクラス$C_1$に分類しそうでない時は$C_2$に分類する。次元を落とした分情報の損失が発生するから$\boldsymbol{w}$を調整してクラスの分離を最大にしていきたい。

ここで、クラス$C_1$の点が$N_1$個、$C_2$の点が$N_2$個あるとすると、それぞれのクラスの平均ベクトルは

\boldsymbol{m}_1 = \frac{1}{N_1}\sum_{n \in C_1}\boldsymbol{x}_n, \quad
\boldsymbol{m}_2 = \frac{1}{N_2}\sum_{n \in C_2}\boldsymbol{x}_n

この時、「クラスの平均同士がもっとも離れるところに射影しよう」という考えに基いて、以下の式を最大にする$\boldsymbol{w}$を選択

m_2 - m_1 = \boldsymbol{w}^T(\boldsymbol{m}_2 - \boldsymbol{m}_1)

ここで、$m_k$は$C_k$から射影されたデータの平均を表す。$\boldsymbol{w}$をいくらでも大きくできると意味がないのでノルム=1という制約を加える。いわゆるラグランジュの未定乗数法の出番。ベクトルの微分の基礎が分かっていれば何の問題もなし。

L = \boldsymbol{w}^T(\boldsymbol{m}_2 - \boldsymbol{m}_1) + \lambda(\boldsymbol{w}^T\boldsymbol{w}-1)\\
\\
\nabla L=(\boldsymbol{m}_2 - \boldsymbol{m}_1)+2\lambda\boldsymbol{w}\\
\\
\boldsymbol{w}=-\frac{1}{2\lambda}(\boldsymbol{m}_2 - \boldsymbol{m}_1)\propto(\boldsymbol{m}_2 - \boldsymbol{m}_1)

ただ実際これではまだクラス同士が重なり合ってしまう場合がある。なので、「射影後に同じクラスはまとまっていて、かつクラス同士は離れている」ような方法をとりたい。そこでフィッシャーの判別基準を導入。各クラスのクラス内分散は

s_k^2 = \sum_{n \in C_k}(y_k - m_k)^2

よって判別基準は以下

J(\boldsymbol{w}) = \frac{(m_2-m_1)^2}{s_1^2 + s_2^2}

分母は総クラス内分散で、各クラスの分散の和で定義。分子はクラス間分散。本節ではこれを次のように書き直している。

J(\boldsymbol{w}) = \frac{\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w}}{\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w}}

ここで

\boldsymbol{S}_\boldsymbol{B} = (\boldsymbol{m}_2 - \boldsymbol{m}_1)(\boldsymbol{m}_2 - \boldsymbol{m}_1)^T\\
\\
\boldsymbol{S}_\boldsymbol{W} =\sum_{k}\sum_{n\in C_k}(\boldsymbol{x}_n-m_k)(\boldsymbol{x}_n-m_k)
^T

前者はクラス間共分散行列、後者は総クラス内共分散行列と呼ばれる。自分には少し取っ付き難い見た目をしていて戸惑ったが、分母も分子も$y=\boldsymbol{w}^T\boldsymbol{x}$であることを利用して展開して見ればもとの式と同じになることがわかる。

よって、J(w)をwに関して微分してゼロとおくことで、Jが最大になるようなwを求められる。


\frac{\partial J}{\partial w}=\frac{(2(\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w})(\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w})-2(\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w})(\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w}))}{(\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w})^2}=0\\
\\\\
(\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w})(\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w}) = (\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w})(\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w})

$\boldsymbol{w}^T\boldsymbol{S}_\boldsymbol{W}\boldsymbol{w}$がスカラーであることと、二次形式を微分する際に共分散行列が対称行列であることを利用しているのがポイントです。これについては今度別の記事に書きます。

先ほどと同じく、今回も重要なのは$\boldsymbol{w}$の向きであって大きさではないので、$\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w}$が

\boldsymbol{S}_\boldsymbol{B}\boldsymbol{w} = (\boldsymbol{m}_2 - \boldsymbol{m}_1)(\boldsymbol{m}_2 - \boldsymbol{m}_1)^T\boldsymbol{w}

より$(\boldsymbol{m}_2 - \boldsymbol{m}_1)$と同じ方向のベクトルであることを利用して

\boldsymbol{w} \propto \boldsymbol{S}_\boldsymbol{W}^-1(\boldsymbol{m}_2 - \boldsymbol{m}_1)

これでwの方向が定まったのでおしまい!

番外編: コードにしてみました

fisher_2d.py
# Class 1
mu1 = [5, 5]
sigma = np.eye(2, 2)
c_1 = np.random.multivariate_normal(mu1, sigma, 100).T

# Class 2
mu2 = [0, 0]
c_2 = np.random.multivariate_normal(mu2, sigma, 100).T

# Average vectors
m_1 = np.sum(c_1, axis=1, keepdims=True) / 100.
m_2 = np.sum(c_2, axis=1, keepdims=True) / 100.

# within-class covariance matrix 
S_W = np.dot((c_1 - m_1), (c_1 - m_1).T) + np.dot((c_2 - m_2), (c_2 - m_2).T)

w = np.dot(np.linalg.inv(S_W), (m_2 - m_1))
w = w/np.linalg.norm(w)

plt.quiver(4, 2, w[1, 0], -w[0, 0], angles="xy", units="xy", color="black", scale=0.5)
plt.scatter(c_1[0, :], c_1[1, :])
plt.scatter(c_2[0, :], c_2[1, :])

quiverを使って求めたベクトルをプロットした結果がこちら

方向はいい感じですね。ということで次回は多クラス版を記事にしようと思います。