概要
多クラスロジスティック回帰
こちらの記事で多クラスロジスティック回帰の概要と、尤度関数から定義される交差エントロピー誤差の勾配を求めました。
多クラスロジスティック回帰における交差エントロピー誤差のヘッセ行列
こちらの記事でさらに微分してヘッセ行列を求めました。
今回は実際にニュートン法を使って多クラスデータの識別境界を求めてみます。
多クラスロジスティック回帰の実験
まずはおさらい。多クラスロジスティック回帰の勾配は
で、ヘッセ行列は
として、
でした。
多変数のニュートン法の更新式は
で、このが式(1)、
が式(3)です。*1 バーがついているのは、
における値という意味です。
式(4)の更新式を使って解いてみた結果です。
ちょっと微妙。ニュートン法のイテレーションを増やすとソフトマックス関数がオーバーフローしてしまい、その前で止めたので十分な識別になっていません。
思いつきでヘッセ行列の逆行列を、擬似逆行列にしてみたらオーバーフローすることなくきちんと識別できました。
なんで擬似逆行列だとオーバーフローしないのかな。ヘッセ行列の逆行列が正しく計算できなかった、というのが根本の問題なのかな?うーん。
# 多クラスロジスティック回帰 import matplotlib.pyplot as plt import numpy as np from matplotlib import cm # 判別結果 def f(x0, x1): return np.argmax(np.dot(w.T, np.array([x0, x1, 1.0]))) # 交差エントロピー誤差 def Ew(y, T): ew = 0 for n in range(N): for k in range(K): ew -= T[n, k] * np.log(y[n, k]) return ew # 教師データの属するクラスを返す def teaching(x, y): if y > x + 0.5: return 0 elif y > x - 0.5: return 1 return 2 # ソフトマックス関数 def softmax(a, k): return np.exp(a[k])/np.sum(np.exp(a)) # 重み付け行列Rを計算 def R_(y, i, j): I = np.eye(N) R_ij = np.empty(N) for k in range(N): R_ij[k] = y[k, i] * (I[i, j] - y[k, j]) return np.diag(R_ij) # ヘッセ行列を計算 def H_(y): H = np.empty([M * K, M * K]) for i in range(K): for j in range(K): # 重み付け行列R R = R_(y, i, j) H[M * j:M * j + M, M * i:M * i + M] = np.dot(np.dot(X.T, R), X) return H # ランダムシードを固定 np.random.seed(0) # 訓練データ数 N = 90 # 入力次元数(ダミー入力を含めて) M = 3 # クラス数 K = 3 # 教師データ作成 x = np.random.uniform(0, 1, [N, M-1]) t = np.vectorize(teaching)(x[:, 0], x[:, 1]) # 行列X(ダミー入力を加える) X = np.hstack([x, np.ones([N, 1])]) # 1 of K符号 T = np.eye(K)[t] # パラメータw初期化 w = np.zeros([M, K]) # ベクトルy y = np.zeros([N, K]) a = np.dot(X, w) for i in range(N): y[i, :] = np.array([softmax(a[i, :], k) for k in range(K)]) # 勾配ベクトルを計算 dw = np.zeros([M*K]) for k in range(K): dw[k * M:k * M + M] = np.dot(X.T, y[:, k] - T[:, k]) # 最急降下法 # w = w.T.flatten() # for it in range(100): # # 交差エントロピー誤差 # print(Ew(y, T)) # w -= 0.1*dw # a = np.dot(X, w.reshape([K, M]).T) # for i in range(N): # y[i, :] = np.array([softmax(a[i, :], k) for k in range(K)]) # for k in range(K): # dw[k * M:k * M + M] = np.dot(X.T, y[:, k] - T[:, k]) # w = w.reshape([M, K]).T # ヘッセ行列 H = H_(y) # ニュートン法 w = w.T.flatten() for ite in range(10): print(Ew(y, T)) # w -= np.dot(np.linalg.inv(H), dw) w -= np.dot(np.linalg.pinv(H), dw) # yを再計算 a = np.dot(X, w.reshape([K, M]).T) for i in range(N): y[i, :] = np.array([softmax(a[i, :], k) for k in range(K)]) # 勾配を再計算 for k in range(K): dw[k * M:k * M + M] = np.dot(X.T, y[:, k] - T[:, k]) # ヘッセ行列を再計算 H = H_(y) w = w.reshape([M, K]).T # グラフ表示用の判別結果 a, b = np.meshgrid(np.linspace(0, 1, 1000), np.linspace(0, 1, 1000)) vec_f = np.vectorize(f) c = vec_f(a, b) plt.xlim(0, 1) plt.ylim(0, 1) plt.scatter(x[t == 0, 0], x[t == 0, 1], color="blue", alpha=0.5) plt.scatter(x[t == 1, 0], x[t == 1, 1], color="green", alpha=0.5) plt.scatter(x[t == 2, 0], x[t == 2, 1], color="red", alpha=0.5) plt.contourf(a, b, c, alpha=0.2, cmap=cm.coolwarm) plt.show()
*1:参考:ニュートン法(多変数の場合)