機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

ロジスティック回帰を最急降下法とニュートン法で解く

まずはおさらい。2クラスの分類問題に確率的生成モデルを適用すれば、出力はシグモイド関数で表すことができ、それはデータがクラスC_kに属する確率として解釈することができました。そして最尤推定によりパラメータ\mathbf{w}を求めるため、交差エントロピー誤差E(\mathbf{w})を定義し、この最適解を求めることを考えます。最急降下法、ニュートン法を用いるために交差エントロピー誤差の勾配、ヘッセ行列を求めると、以下のようになることを確認しました。(参考: ロジスティック回帰における交差エントロピー誤差のヘッセ行列ロジスティック回帰確率的生成モデル(1)確率的生成モデル(2)


\displaystyle \frac{\partial E(\mathbf{w})}{\partial \mathbf{w}} = \sum_{n=1}^{N}(y_n - t_n)\mathbf{x}_n \tag{1}


\mathbf{H} =  \mathbf{X}^{T} \mathbf{R} \mathbf{X} \tag{2}

今回はこれを用いて実際に解いてみたいと思います。

多変数のニュートン法の更新式は


\mathbf{x} = \mathbf{x}^{\prime} - \bar{\mathbf{H}}^{-1} \nabla \bar{f} \tag{3}

で、この\nabla fが式(1)、\mathbf{H}が式(2)です。(参考:ニュートン法(多変数の場合))バーがついているのは、 \mathbf{x}^{\prime}における値という意味です。

早速、適当なデータを作ってロジスティック回帰を試してみたのが下のグラフ。 f:id:opabinia2:20180816130935p:plain きちんと分類できていますね。

次に、最小二乗法による線形識別では正しく分類できなかったパターンを試してみます。(参考:線形識別を最小二乗法で解くf:id:opabinia2:20180816130943p:plain こちらも正しく分類できています。最小二乗法のように外れ値に引っ張られたりすることがなく、偏ったデータでもきちんと分類できるみたいです。

最後に、ニュートン法と最急降下法の最適化のスピードの違いを見てみます。

f:id:opabinia2:20180816130950p:plain

最急降下法は途中でサチってしまっているように見えますが、反復回数を増やせば下がっていきます。ちなみに10万回反復したら4.0E-3くらいまでになりました。学習率を変えると多少スピードが変化しますが、どんなパラメータにしたところでニュートン法のスピードには遠く及びませんでした。ニュートン法すげー。

入力を非線形変換したロジスティック回帰も試してみました→入力を非線形変換したロジスティック回帰

今回のコードです。

# ロジスティック回帰
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


# 判別結果
def f(x0, x1):
    if sigmoid(np.dot(w, np.array([x0, x1, 1.0]))) < 0.5:
        return 0
    else:
        return 1


# 交差エントロピー誤差
def Ew(y, t):
    ew = 0
    for i in range(N):
        ew -= np.log(t[i]*y[i] + (1-t[i])*(1-y[i]))
    return ew


# 教師データの属するクラスを返す
def teaching(x, y):
    if y > x+0.5:
        return 0
    else:
        return 1


# シグモイド関数
def sigmoid(a):
    return 1.0/(1 + np.exp(-a))


# 重み付け行列Rを計算
def R_(y):
    R_ij = np.empty(N)
    for i in range(N):
        R_ij[i] = y[i] * (1 - y[i])
    return np.diag(R_ij)


# ランダムシードを固定
np.random.seed(0)

# 訓練データ数
N = 50

# 入力次元数(ダミー入力を含めて)
M = 3

# クラス数
K = 2

# 最急降下法の学習率
ALPHA = 1.0

# 教師データ作成
T = np.zeros([N, K])
t = np.empty(N)
xt = np.random.uniform(0, 1, [N, M-1])
for i in range(N):
    t[i] = teaching(xt[i, 0], xt[i, 1])
    T[i, int(t[i])] = 1

# 行列X(ダミー入力を加える)
X = np.hstack([xt,  np.ones([N, 1])])

# パラメータw初期化
w = np.random.uniform(-1, 1, M)

# 出力y
y = np.vectorize(sigmoid)(np.dot(w, X.T))

# 重み付け行列R
R = R_(y)

# ヘッセ行列
H = np.dot(np.dot(X.T, R), X)

# ニュートン法
for i in range(20):
    w -= np.dot(np.dot(np.linalg.inv(H), X.T), (y -t))
    # yを再計算
    y = np.vectorize(sigmoid)(np.dot(w, X.T))
    # R,Hを再計算
    R = R_(y)
    H = np.dot(np.dot(X.T, R), X)
    print(Ew(y, t))

# 最急降下法
# for i in range(20):
#  w -= ALPHA*np.dot(y - t, X)
#  # yを再計算
#  y = np.vectorize(sigmoid)(np.dot(w, X.T))
#  print(Ew(y, 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(xt[t == 0, 0], xt[t == 0, 1], color="blue", alpha=0.5)
plt.scatter(xt[t == 1, 0], xt[t == 1, 1], color="red", alpha=0.5)
plt.contourf(a, b, c, alpha=0.2, cmap=cm.coolwarm)
plt.show()