Python
numpy
機械学習

行と列を意識しながら多クラスロジスティック回帰を実装する

More than 1 year has passed since last update.

最初に

numpyでイチから書いてみると理解が捗るぞ〜ということではじパタ片手に実装してみました。
ただ、ベクトルで定式化されていることで中々手こずりました。
数式をコードに落とし込む場合、しっかりと行×列を意識することが重要だと実感。。。
今回は行と列を中心に見ながら多クラスロジスティックについてまとめていきたいと思います。

やりたいこと

例のごとくIrisを分類したい!
2クラスではなく3クラスを分類したい!したい!したい!したい!
(線形分離できそうなpatal lengthとpatal widthを使います)

iris.png

数式の確認

まずは、細かい導出は除いてベクトル化された状態でざっと見ていきます。

仮説関数

h_k = x\theta_k^T

この場合、3クラスなので$k=1,2,3$です。

ソフトマックス関数

P(C_k|x) = \pi_k(x) = \frac{e^{h_k}}{\sum_{j=1}^Ke^{h_j}}

多クラスなので2クラス分類で使用していたシグモイド関数ではなくソフトマックス関数を使用して事後確率を求めます。各クラスで事後確率を計算し、最大事後確率を返すクラスへ分類します。

コスト関数(負の対数尤度関数)

E(\theta_1, ・・・, \theta_K) = - \frac{1}{M} \sum_{i=1}^N \sum_{k=i}^K t_{ik} \log\pi_{ik}

ここで、$t_{ik}$はoneHotEncodingされた正解ラベルになります。
Irisデータセットでは初期のクラスラベルが下図のように0, 1, 2で表現されています。

スクリーンショット 2017-12-13 18.53.23.png

これをこんな感じにしてあげる必要があります。

スクリーンショット 2017-12-13 18.54.32.png

このようなoneHotEncodingの実装についてはこちらでまとめてあります。
ありがたいことにたくさんの方からスタイリッシュなコメントをいただいています。
https://qiita.com/hatunina/items/78656e8736b93a980acb

重みパラメータの更新

\theta_k = \theta_k - \alpha\frac{\partial E}{\partial \theta_k} \\
= \theta_k - \alpha\frac{1}{M}\sum_{i=1}^N (\pi_{ik} - t_{ik})x_i

これを最急降下法の中で使用して$\theta$を求めていきます。

数式を行列で見る

ここからはfor文を使わずに行列でカッコよく計算するために、ベクトル化されていた数式を行×列を意識しながら見ていきます。

ここで、各パラメータを今回の場合に当てはめると下記のようになります。
クラスは3つ、説明変数はpatal lengthとpatal widthに定数項を足して3つです。
トレーニングセットは各クラスから30個ずつデータを抽出し合計90個としました。

K = クラス数 = 3
N = 説明変数 = 3(定数項含む)
M = トレーニングセット数 = 90

仮説関数

h_k = x\theta_k^T \\
ベクトルから行列へ \\
↓  \\
h = x\theta^T \\
右辺を行×列で整える\\
↓  \\
 M×N・(K×N)^T \\
↓  \\
 M×N・N×K \\
↓  \\
 M×K

$x$はデータ数×説明変数、$\theta$はクラス数×説明変数です。
出力はドット積を用いることでデータ数×クラス数となり、各データ、各クラスの仮説関数作ることができます。

def hypothesis(X, theta):
    h = np.dot(X, theta.T)
    return h

ソフトマックス関数

P(C_k|x) = \pi_k(x) = \frac{e^{h_k}}{\sum_{j=1}^Ke^{h_j}}\\
ベクトルから行列へ\\
↓  \\
P(C|x) = \pi(x) = \frac{e^{h}}{\sum_{i=1}^Me^{h_i}}\\
右辺を行×列で整える\\
↓  \\
\frac{M×K}{M×1} \\
↓  \\
 M×K

分母が少しやっかいですね。$e^h$は$M×K$行列なので行ごとに列方向に合計を取ります。すると、$M×1$の行列になるので、これで$e^h$の各要素を割ります。出力は$M×K$行列になります。データごとに各クラスの事後確率を持っている形になります。

下記コードの np.sum へ axis=1 の引数を持たせることで列方向の合計が取れます。

def softmax(h):
    return np.exp(h) / np.sum(np.exp(h), axis=1).reshape([h.shape[0], 1])

コスト関数(負の対数尤度関数)

E(\theta_1, ・・・, \theta_K) = - \frac{1}{M} \sum_{i=1}^N \sum_{k=i}^K t_{ik} \log\pi_{ik} \\
ベクトルから行列へ\\
↓  \\
E(\theta) = - \frac{1}{M} \sum t \log\pi \\
右辺を行×列で整える\\
↓  \\
\sum (M×K) × (M×K) \\
↓  \\
\sum (M×K) \\
↓  \\
1つの値

コスト関数は全データの予測と正解の誤差を出力するため、最終的にその時点の$\theta$で計算した1つの数値を返します。ここでは、$\pi$の各要素に対応した$t$の積を求め、それの全ての和を取るためドット積ではなく、対応した要素ごとの積を使用しています。

図解するとこんな感じ。

スクリーンショット 2017-12-16 18.55.45.png

# trainYはonehot
def costFunction(trainX, trainY, theta, M):
    # M*K行列
    h = hypothesis(trainX, theta)
    # M*K行列
    pi = softmax(h)
    # M*K と M*K の要素ごとの"積"を求め合計を出す
    cost = (-1/M) * np.sum(trainY * np.log(pi))
    return cost

重みパラメータの更新

\theta_k = \theta_k - \alpha\frac{1}{M}\sum_{i=1}^N (\pi_{ik} - t_{ik})x_i \\
ベクトルから行列へ\\
↓ \\
\theta = \theta - \alpha\frac{1}{M} (\pi - t)x \\
右辺を行×列で整える\\
↓  \\
K×N - (M×K - M×K)^T・M×N\\
↓  \\
K×N - (M×K)^T・M×N\\
↓  \\
K×N - K×M・M×N\\
↓  \\
K×N - K×N\\
↓  \\
K×N

シグマはドット積で代用することができます。
また、ここでは、$\theta$が$K×N$行列のため、$K×N$を作るという意識で整えてもいいかもしれません。

def updateTheta(trainX, trainY, theta, alpha, M):
    gradList = []
    # M*K
    h = hypothesis(trainX, theta)
    # M*K
    pi = softmax(h)
    # M*K.T ・ M*N → K*N 
    grad = np.dot((pi-trainY).T, trainX)
    # theta更新 K*N 
    theta = theta - (alpha/M) * grad

    return theta

最急降下法

必要な数式は上で見てきたので、あとは最急降下法を書いてあげるだけです。さすがにここのfor文は外せません。
...外せないよね??

def gradientDescent(trainX, trainY, theta, iter_num, alpha, M):
    costList = []

    for i in range(iter_num):
        cost = costFunction(trainX, trainY, theta, M)
        costList.append(cost)
        theta = updateTheta(trainX, trainY, theta, alpha, M)

    return theta, costList
# ハイパーパラメータ設定と最急降下法の実行
K = 3
N = 3
theta = np.ones([N, K])
iter_num = 10000
alpha = 0.1
M = len(trainX)

theta, costList = gradientDescent(trainX, trainY, theta, iter_num, alpha, M)

結果

先ほどのcostListを描画してみます。

スクリーンショット 2017-12-16 20.00.05.png

いい感じに下がってくれているので、予測と正解率も計算してみます。

スクリーンショット 2017-12-16 20.05.42.png

正解率96%となりました。
テストセットは各クラスを20個ずつ抽出しており、そのまま順番通りになっているので、どうやらクラス2をクラス1と誤分類している箇所が2つあるっぽいですね。

ということで、最後に識別境界を描画してみます。

スクリーンショット 2017-12-16 20.13.11.png

やはり緑のvirginicaが二つ分類できていませんね!
まあ線形なのでしょうがないです。

最後に

数式をコードに落とし込む際に自分がハマった観点からまとめてみました。
次元を確認する上で紙とペンで数式を解いていくということも大事なのかなと思います。
ここちょっとおかしくない?とかもっとスタイリッシュに書けるでしょ!という意見がありましたら、ぜひコメントいただけると嬉しいです!