【ロジスティック回帰の思わぬ罠！】学習結果の発散とその回避方法について

ロジスティック回帰 f:id:wakaba-mafin:20181118205651j:plain:w150:right

気軽に使えるクラス分類の基本的なツール。

しかし、油断しているとうっかりはまってしまう落とし穴があります。うまく学習率を調整しても、学習結果が発散する場合があるという罠です。

今回の記事ではこの罠とその回避方法についてご紹介します！

Logistic回帰とは
Logistic回帰の罠
原因の解説
- 線形分離可能性
- 係数が発散する場合のロス関数
回避方法
まとめ

Logistic回帰とは

Logistic回帰については前回の記事

wakaba-mafin.hatenablog.com

にて紹介いたしました。簡単に言ってしまうと、Logistic回帰とはモデル

$p(C|x) = \sigma(w^{T} x) := \frac{1}{1+\exp(-w^{T} x)}$

に基づいて、特徴 $x$ からクラス $C$ を識別する手法です。

特にこの識別では特徴空間に超平面を張ることによって線形にクラスを識別します。

計算の方法も容易で多数のパッケージが提供されています！

Logistic回帰の罠

Logistic回帰は非常に扱いやすい識別器でありますが、実は油断しているととんでもない罠にはまります。ここではその一例をPythonを使ってご紹介します。

今回使うライブラリはこちら：

# library
import numpy as np
import matplotlib.pyplot as plt

とりあえず、numpyとmatplotlibにとどめておきましょう。

で、ここで問題が生じるデータは次のようにして作ります！

num_of_data = 10000
# 特徴ベクトル
X = np.random.randn(num_of_data,2)
# 教師ラベル
C = (X@np.ones(2)>0).astype(dtype=np.int32)
# 特徴ベクトルにバイアスに関する項を加える
Xb = np.hstack([X, np.ones((num_of_data,1))])

これはこんな感じのデータ。ちょうど $X_1 + X_2 = 0$ を境界とするようなデータですね！

plt.scatter(X[C==1, 0], X[C==1, 1], c="r", label="label 1", alpha=0.25)
plt.scatter(X[C==0, 0], X[C==0, 1], c="b", label="label 0", alpha=0.25)
plt.legend()
plt.xlabel("X1")
plt.ylabel("X2")
plt.show()

f:id:wakaba-mafin:20181219234951p:plain

そして、Logistic回帰に関する前回の記事

wakaba-mafin.hatenablog.com

に従って、このデータに対してLogistic回帰を適用してみると...

def sigmoid(a):
    return 1.0/(1.0 + np.exp(-a))

eta = 0.01
num_of_learning = 10000
w = np.zeros((num_of_learning, 3))
w[0,:] = np.random.randn(3)
for n in range(1, num_of_learning):
    # Gradient Descent Method
    w[n,:] = w[n-1,:] - eta * Xb.T @ (sigmoid(Xb@w[n-1,:])-C)
    
print("result : ", w[-1,:])

result :  [90.91425211 91.22758037  0.10616866]

なるほど。少し係数が大きいような気もしないでもないですが... ま、目的の係数 $(1,1,0)$ の定数倍程度なのでよしとしましょう！

$w$ の学習回数ごとの推移を見てみましょう。Parameter 0,1,2がそれぞれw[0], w[1], w[2]に対応しています。

for i in range(w.shape[1]):
    plt.plot(w[:,i], "-", label="parameter {0}".format(i), lw=9-4*i)
plt.xlabel("iteration")
plt.ylabel("parameter")
plt.legend()
plt.show()

f:id:wakaba-mafin:20181219235037p:plain

おやおや、まだ係数は収束してなさそうですね...

もう少し学習を回してみましょう。

eta = 0.01
num_of_learning = 100000
w = np.zeros((num_of_learning, 3))
w[0,:] = np.random.randn(3)
for n in range(1, num_of_learning):
    # Gradient Descent Method
    w[n,:] = w[n-1,:] - eta * Xb.T @ (sigmoid(Xb@w[n-1,:])-C)
    
print("result : ", w[-1,:])

result :  [1.92573126e+02 1.93011457e+02 1.04457667e-01]

for i in range(w.shape[1]):
    plt.plot(w[:,i], "-", label="parameter {0}".format(i), lw=9-4*i)
plt.xlabel("iteration")
plt.ylabel("parameter")
plt.legend()
plt.show()

f:id:wakaba-mafin:20181219235125p:plain

なるほど。まだまだ増加して発散しそうですね...

ところで、ロスはどうなっているんでしょう？

def crossentropy(x,y):
    return -x[y>0]@np.log(y[y>0])-(1-x[y<1])@np.log(1-y[y<1])
    
plt.plot(1+10*np.arange(w[::10].shape[0]), np.array([crossentropy(C,sigmoid(Xb@v)) for v in w[::10]])/num_of_data, "ko-")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("iteration")
plt.ylabel("loss(cross-entropy)")
plt.show()

f:id:wakaba-mafin:20181219235232p:plain

学習を重ねても減少していて収束してなさそうですね...学習データに過剰にフィッティングしていることを疑って、別のデータを新たに生成しても...

num_of_new_data = 100
# 特徴ベクトル
newX = np.random.randn(num_of_data,2)
# 教師ラベル
newC = (X@np.ones(2)>0).astype(dtype=np.int32)
# 特徴ベクトルにバイアスに関する項を加える
newXb = np.hstack([X, np.ones((num_of_data,1))])

plt.plot(1+10*np.arange(w[::10].shape[0]), np.array([crossentropy(newC,sigmoid(newXb@v)) for v in w[::10]])/num_of_data, "ko-")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("iteration")
plt.ylabel("loss(cross-entropy)")
plt.show()

f:id:wakaba-mafin:20181219235308p:plain

学習データでも検証データでも、学習回数ごとにロスは正しく減少しています。このままでは学習結果の $w$ が発散してしまいます...💦

一体何が起こっているのでしょうか？

原因の解説

実はこの問題はデータが持っている性質「線形分離可能性」によって生じています。

いやいや、線形分離可能性って何やねんって話ですね笑

ということでまずは線形分離可能性についてお話します。

線形分離可能性

線形分離可能性とは、その名が示す通りデータを線形に分離できることを指します。つまり、特徴空間内の２つのクラスのデータをある超平面によって分離することができるということです。

例えば、二次元平面上でクラス1のデータは全て第一象限に、クラス0のデータは全て第三象限にあったとしましょう（境界線は含みません）。

f:id:wakaba-mafin:20181219235340p:plain

このとき、これらのデータは

$x + y = 0$

という直線（超平面に対応）によって分離することができます。このような場合をデータが線形分離可能といいます。

一般的な次元 $n$ において数式を以て線形分離可能性を表現をするのであれば、ある $a=(a_1, \cdots, a_n) \in \mathbb{R}^{n}$ と定数 $b\in\mathbb{R}$ が存在して、

$a_1 x_1 + \cdots + a_n x_n \gt b \Rightarrow x\mbox{はクラス1}$

$a_1 x_1 + \cdots + a_n x_n \lt b \Rightarrow x\mbox{はクラス0}$

が成り立つときに、 $x$ とクラス $C$ の組み合わせである学習データは線形分離可能であると言います。

係数が発散する場合のロス関数

さて、線形分離可能な場合に係数が発散することを簡単な例で確認しましょう！（※危ない議論もありますが...）

今、“無限個”のデータ $\{(x_i, C_i)\}_{i=1}^{\infty}$ （ $i$ に関して独立）があるとして、これが線形分離可能としましょう。特に簡単のため、特徴 $x$ は一次元としておき、クラスは決定論的に

$C = 1 \;\; {\rm if} \;\; x \gt 0\\ C = 0 \;\; {\rm if} \;\; x \lt 0$

によって決まっていて、 $x\sim\mathcal{N}(0,1)$ を満たしているとしましょう。このとき、ロジスティック回帰でモデル

$p(C|x) = \sigma(wx)$

の $w$ を求めてみましょう！（ここでも簡単のためバイアスに相当する項は無視しています）

ロジスティック回帰におけるロス関数は対数尤度、つまりクロスエントロピーで定義されます。特に今、“無限個”のデータを考えているので、ロス関数は

$l(w|(x_1, C_1), \cdots) = \mathbb{E}_X \left[ -C \log \sigma(wX) -(1-C)\log(1-\sigma(wX)) \right]$

となります。もう少し詳しく書くと

$l(w|(x_1, C_1), \cdots) = -\int_{-\infty}^{0} \log(1-\sigma(wx)) n(x)dx -\int_{0}^{\infty} \log(\sigma(wx)) n(x)dx$

となります（ここで $n(x)$ は標準正規分布の密度関数です）。

このロス関数はグラフを書くと...

f:id:wakaba-mafin:20181219235605p:plain

となります。つまり、 $w$ が大きくなればなるほど、ロス関数 $l$ がどんどん小さくなっていきます！

実際に $l$ を $w$ で微分してみると（ちょっと怪しいですが...💦）

$\frac{dl}{dw} = \int_{-\infty}^{0} x\sigma(wx) n(x)dx -\int_{0}^{\infty} x(1-\sigma(wx))n(x)dx \lt 0$

となって、常に単調減少であることがわかります。したがって、この簡単な例ではロス関数を最小化しようとして $w$ を求めたとしても、それよりも大きい $w$ の方がさらにロス関数が小さいので、永遠に $w$ を大きくしていってしまうという問題が生じます。

ま、もっと簡単に説明するのであれば、この例の場合、実は真のモデルが存在して

$p(C=1|x) = H(x)$

となります。ここで $H(x)$ はHeavisideのステップ関数で、ロジスティックシグモイド関数 $\sigma(wx)$ でそれを表現しようとすると $w\rightarrow \infty$ としなければならないということが簡単な説明ですかね？

f:id:wakaba-mafin:20181219235747p:plain

以上のようなことが、他の線形分離可能な場合にも生じて、学習結果の係数が発散するという問題が起こっています。

回避方法

このような問題の回避策の１つとしては事前分布を用意することが挙げられます。これはロジスティック回帰のロス関数（対数尤度）にパラメータ $w$ の正則化項を導入することに相当します。

正則化項としては $l_1$ 正則化項や $l_2$ 正則化項など色々あります。

$\tilde{l}_1(w|\{(x_i,C_i) \}_{i=1}^{N}) = -\frac{1}{N}\sum_{i=1}^{N}\left\{ C_i \log \sigma(wx_i) +(1-C_i)\log \left( 1-\sigma(wx_i)\right)\right\} + \lambda \sum_{j} |w_j|$

$\tilde{l}_2(w|\{(x_i,C_i) \}_{i=1}^{N}) = -\frac{1}{N}\sum_{i=1}^{N}\left\{ C_i \log \sigma(wx_i) +(1-C_i)\log \left( 1-\sigma(wx_i)\right)\right\} + \lambda w^{T}w$

このようにすれば正則化項が $w$ があまり大きくならないように働き、 $w$ の発散を抑えることができるようになります！

例えば、正則化係数 $\lambda=0.1$ の $l_2$ 正則化項を一番最初に示した例に適用するのであれば、次のようになります：

num_of_data = 10000
# 特徴ベクトル
X = np.random.randn(num_of_data,2)
# 教師ラベル
C = (X@np.ones(2)>0).astype(dtype=np.int32)
# 特徴ベクトルにバイアスに関する項を加える
Xb = np.hstack([X, np.ones((num_of_data,1))])

# Logistic回帰の学習
eta = 0.01               # 学習率
lam = 0.1                # 正則化係数
num_of_learning = 10000  # 学習回数
w = np.zeros((num_of_learning, 3))
w[0,:] = np.random.randn(3)
for n in range(1, num_of_learning):
    # Gradient Descent Method
    w[n,:] = w[n-1,:] - eta * Xb.T @ (sigmoid(Xb@w[n-1,:])-C) - 2*eta*lam*w[n-1,:]
    
print("result : ", w[-1,:])

for i in range(w.shape[1]):
    plt.plot(w[:,i], "-", label="parameter {0}".format(i), lw=9-4*i)
plt.xlabel("iteration")
plt.ylabel("parameter")
plt.legend()
plt.show()

result :  [22.59205216 22.51610759  0.15557752]

f:id:wakaba-mafin:20181219235948p:plain

今度はちゃんと学習した係数が収束しましたね！

これで問題解決です！

まとめ

というわけで今回の記事では、ロジスティック回帰でうっかりしてるとはまってしまう罠についてご紹介いたしました！

学習データ自体が線形分離可能であると、ロジスティック回帰では無限によい識別器を作ることができるので学習結果が発散してしまうという問題です！この回避方法の１つとしては正則化項を加えると割とうまく行くということをご紹介しました。

ロジスティック回帰で発散して何か変だなって思った際には、ぜひ思い出していただければうれしいです！

なお、Pythonの機械学習ツールscikit-learnのロジスティック回帰ではデフォルトで正則化項がついているので安心です。

また、今回の記事の内容は機械学習のバイブル「Pattern Recognition and Machine Learning」でも触れられていますので、ぜひ興味のある方は一度読んでみるといいと思います！↓

Amazon CAPTCHA

というわけで、とやかく書きましたが、この辺りで執筆を終わりたいと思います。

ご質問やご意見などございましたら、ぜひコメントのほどよろしくお願いします！

若葉の技術メモ

若葉の技術メモ

コンピュータやプログラミング・数理に関して調べたり、取り組んだりしたことをまとめる若葉のノート。

コンピュータ・プログラミング・数理が初めてって方も一緒に勉強していきましょう！

初心は大事！いつでも若葉☘のような意気込みで！

【ロジスティック回帰の思わぬ罠！】学習結果の発散とその回避方法