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

概要

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

$\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}$

以下がそれらについて記載した記事です。 www.iwanttobeacat.com

www.iwanttobeacat.com

実験結果

今回はこれを用いて実際にロジスティック回帰で分類してみたいと思います。

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

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

で、この $\nabla f$ が式(1)、 $\mathbf{H}$ が式(2)です。*1　バーがついているのは、 $\mathbf{x}^{\prime}$ における値という意味です。

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

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

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

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

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

www.iwanttobeacat.com

今回のコードです。

# ロジスティック回帰
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()