2019-06-29

NumPyでロジスティック回帰を実装する

Python 機械学習

ディープラーニング以前 (〜2010年) の機械学習について、はじパタを使って整理しています。

今回は「第6章線形識別関数」を参考にしながら、ロジスティック回帰をNumPyで実装してみます。

はじめてのパターン認識

作者: 平井有三
出版社/メーカー: 森北出版
発売日: 2012/07/31
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 7回
この商品を含むブログ (5件) を見る

ロジスティック回帰

実装に入る前に、ロジスティック回帰の理論を整理します。ここでは、2値 $\{C_1, C_2\}$ に分類するタスクを考えます。

シグモイド関数

まずクラス $C_1$ に属する確率を計算するための定式化を行います。

2値分類ですので、 $C_1$ の確率と $C_2$ の確率の合計は1です

入力 $\vec{x}$ (列ベクトル) が与えられたときにクラス $C_1$ に分類される事後確率は以下となります。

$P(C_1 | \vec{x}) = \frac{P(\vec{x} | C_1)P(C_1)}{P(\vec{x} | C_1)P(C_1) + P(\vec{x} | C_2)P(C_2)}$

このときaを導入すると、

$a := ln \frac{P(\vec{x} | C_1)P(C_1)}{P(\vec{x} | C_2)P(C_2)}$

事後確率 $P(C_1 | \vec{x})$ は以下のシンプルな式になります。

$P(C_1 | \vec{x}) = \frac{1}{1 + exp(-a)} = \sigma(a)$

$\sigma(a)$ をシグモイド関数といいます。シグモイド関数の値は [0, 1] となりますので、確率を表現するために好都合な関数です。

このシグモイド関数の逆関数をロジット関数といいます。

$a = ln \frac{\sigma(a)}{1 - \sigma(a)} = ln \frac{P(C_1 | \vec{x})}{P(C_2 | \vec{x})}$

ロジスティック回帰モデル

次にaの計算方法を定義します。

ここで重み $\vec{w}$ (列ベクトル) を導入し $a=\vec{w}^T\vec{x}$ で表現します。

$\sigma(a) = \frac{1}{1 + exp(-\vec{w}^T\vec{x})} = \frac{exp(\vec{w}^T\vec{x})}{1 + exp(\vec{w}^T\vec{x})}$

ロジスティック回帰モデルは、二項分布をリンク関数とした一般化線形モデルに分類されます。こちらの投稿も参考にしてみてください。緑本とstatsmodelsを使って、他の一般化線形モデルと比較しています。

Pythonで実装しながら緑本を学ぶ (第6章 GLMの応用範囲を広げる -ロジスティック回帰など-) - け日記

交差エントロピー誤差関数

最後に、重み $\vec{w}$ を教師データから獲得するために、最適化すべき尤度関数を定義します。

あるクラスに分類される確率で表現できる2値分類のため、尤度関数として交差エントロピー誤差関数を導入します。

$N$ は教師データ数
$t_i$ は教師データ
$y_i$ は入力 $\vec{x_i}$ に対して予測した確率

$L(\{y_1, ..., y_N\}) = \prod_{i=1}^{N} y_{i}^{t_i}(1 - y_{i})^{(1-t_{i})}$

この尤度関数を対数化して、負の符号をつけます。

$lL(\{y_1, ..., y_N\})=-\sum_{i=1}^{N}(t_{i} ln y_{i} + (1-t_i) ln (1-y_i))$

$y_i$ に代入します。

$lL(\vec{w})=-\sum_{i=1}^{N}(t_{i} \vec{w}^T \vec{x_i} - ln(1 + exp(\vec{w}^T\vec{x_i})))$

これを $\vec{w}$ で微分します。この $\frac{\partial lL}{\partial \vec{w}}$ が0となるように $\vec{w}$ を学習していきます。

$\frac{\partial lL}{\partial \vec{w}} = -\sum_{i=1}^N (t_i \vec{x_i} - \frac{\vec{x_i}exp(\vec{w}^T \vec{x_i})}{1 + exp(\vec{w}^T\vec{x_i})}) = \sum_{i=1}^N \vec{x_i}(y_i - t_i)$

最急降下法

学習方法として、今回は最急降下法で実装します。 $\alpha$ は更新量をコントロールするためのハイパパラメータです。

初期化: ランダムな重み $\vec{w_0}$ を設定
$w$ で勾配 $\frac{\partial dL}{\partial \vec{w}}$ を計算
$w$ <- $\vec{w} - \alpha \frac{\partial dL}{\partial \vec{w}}$ で重みを更新
収束するまで2.3.を繰り返す

NumPyで実装

scikit-learn付属のirisデータセットを使って、上を実装してみます。

入力は3次元 (萼片の長さ・幅とバイアス)
- 重み $\vec{w}$ も3次元
setosaとversicolorのデータのみを使う
- versicolorの場合に $t_i = 1$ とする

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

# irisデータセットのロード
iris = load_iris()

# setosaとversicolorの萼片の長さ・幅のみを使う
X = iris.data[iris.target != 2][:, :2].T
X = np.concatenate([X, np.ones((1, X.shape[1]))])  # バイアス項
t = iris.target[iris.target != 2]

# 分布をチェック
plt.scatter(X[0][t == 0], X[1][t == 0], label=iris.target_names[0]) 
plt.scatter(X[0][t == 1], X[1][t == 1], label=iris.target_names[1])
plt.legend()
plt.show()

f:id:ohke:20190626111425p:plain

学習用のコードは以下です。上の式をナイーブに実装しています。

# wを乱数で初期化
def initialize_w():
    return np.random.normal(0, size=3).reshape((3, 1))

# versicolorの確率を予測
def predict_probability(w, X):
    return (np.exp(np.dot(w.T, X)) / (1.0 + np.exp(np.dot(w.T, X)))).reshape(-1)

# 損失関数
def loss_function(w, X, t):
    loss = 0

    for i in range(t.shape[0]):
        loss += t[i]*np.dot(w.T, X[:, i]) - np.log(1 + np.exp(np.dot(w.T, X[:, i])))

    return -loss

# 損失関数の勾配
def gradient_loss(w, X, t):
    y = predict_probability(w, X)
    n = t.shape[0]
    loss = np.zeros_like(w)
  
    for i in range(n):
        loss += ((y[i] - t[i]) * X[:, i])[:, np.newaxis]

    return loss

# 最急降下法で学習
def gradient_descent(X, t, alpha=0.001):
    w_t = initialize_w()
    loss_t = loss_function(w_t, X, t)

    w_history, loss_history = [w_t], [loss_t]

    for i in range(1000000):  # lossに変化がなくなるまでループする
        w_t -= gradient_loss(w_t, X, t) * alpha
        loss_t = loss_function(w_t, X, t)

        # 収束判定
        if loss_t >= loss_history[-1]:
            break

        w_history.append(w_t)
        loss_history.append(loss_t)
        
    return w_history, loss_history


# 学習
w_history, loss_history = gradient_descent(X, t)

損失関数の値の推移はこちらで、最初にぐっと下がって、あとはダラダラと少しずつ下げていってることがわかります。

f:id:ohke:20190626115100p:plain

得られた重みを使って分類させてみましたところ、精度100%で完全に分離できていることを確認できました。

# 重みの確認
w_hat = w_history[-1]
print(w_hat)
# [[ 12.30194381]
#  [-12.92596168]
#  [-26.27142723]]

# 得られた重みで確率を計算
y_probability = predict_probability(w_hat, X)

# 精度
y = np.where(y_probability > 0.5, 1, 0)
sum(y == t) / t.shape[0]  # 1.0

まとめ

今回はNumPyでナイーブなロジスティック回帰を実装しました。

2019-06-22

ResNetでCIFAR-10を分類する

CV Python Keras

KerasでResNetを作ってCIFAR-10を分類し、通常のCNNモデルと比較します。

ResNet

ResNetはCNNのモデルの1つです。
Microsoft ResearchのKaiming Heらが2015年に提案¹し、その年のILSVRCではResNetで学習したモデルが優勝しました。

VGGやGoogLeNetにて、畳み込み層を重ねることでより良い感じの特徴抽出ができることが明らかになっていましたが、多層になればなるほど前段の層で勾配消失問題が顕著になります。

ResNetは勾配消失問題を解消するために、複数の畳み込み層をスキップするshortcut connectionを導入しました。このresidual blockの図を以下に示します。

今までの畳み込み層は、入力を各レイヤに通すdeep path (青) のみ
ResNetでは、deep pathの出力と、入力をそのまま伝搬するshortcut connection (赤) を加算する

学習過程で畳み込み層F(x)は、畳み込み層の出力H(x)をそのまま学習するのではなく、入力xとの残差H(x)-xで重みを学習するようになります (ResNetのResidualはここから来てます) 。

deep pathを通過することで勾配が小さくなっていってましたが、shortcut connectionで畳み込み層をスキップできるので、前段の層にも勾配を伝搬可能となった、ということです。

論文中では、deep path内でサイズを変えて合流させるBottleneckアーキテクチャも触れられていましたが、今回はそれを行わないPlainアーキテクチャを組みます

また、以下の特徴も性能へ好影響を与えていることが指摘されてます。

表現力が高すぎて役に立たない畳み込み層を残差0でスキップできる
複数のパスを通しているためアンサンブル学習に近い

CIFAR-10

32x32のカラー画像 (3チャネル) を10クラスに分類するタスクで、Kerasのデータセット付属のデータセットです。

学習データ50,000サンプル、テストサンプル10,000サンプルからなります
各クラスのサンプル数は均等

Kerasで実装

以下の3パターンを実装します。

thinモデル: 畳み込み層x3 (Shortcut connection無し, 出力層8x8x32)
plainモデル: 畳み込み層x13 (Shortcut connection無し, 出力層4x4x64)
residualモデル: 畳み込み層x13 (Shortcut connection有り, 出力層4x4x64)

import matplotlib.pyplot as plt
import keras
from keras import backend as K
from keras.models import Sequential, Model
from keras.layers import Input, Conv2D, MaxPool2D, BatchNormalization, ReLU, Flatten, Dense, Add, Dropout
from keras.datasets import cifar10

# データのロード
(X_train, y_train), (X_test, y_test) = cifar10.load_data()
print(X_train.shape, y_train.shape, X_test.shape, y_test.shape)
# (50000, 32, 32, 3) (50000, 1) (10000, 32, 32, 3) (10000, 1)

# 0〜1へ正規化
X_train = X_train.astype('float32') / 255.0
X_test = X_test.astype('float32') / 255.0

# one-hotエンコーディング
Y_train = keras.utils.to_categorical(y_train)
Y_test = keras.utils.to_categorical(y_test)

input_shape = X_train.shape[1:]
train_samples = X_train.shape[0]
test_samples = X_test.shape[0]

# ハイパパラメータ
epochs = 200
batch_size = 50

畳み込み層以外の実装

最初に畳み込み層以外の実装を見ていきます。共通して使うモデルため、生成関数generate_modelを定義してます。

最初の畳み込み層で (32, 32, 3) -> (16, 16, 32) にします
畳み込み層はblock_sets*blocks*block_layers層で構成します
- block_fでブロックの実装を関数で渡します (後述)
- 1ブロックセット終了時に、2x2でプーリングし、次のブロックセットのフィルタ数を倍に拡張
全結合層で畳み込み層の出力 -> 100 -> 10 で絞っていきます
畳み込み後にBatchNormalization, ReLU, Dropoutを行う

# モデルの生成
def generate_model(input_shape, block_f, blocks, block_sets, block_layers=2, first_filters=32, kernel_size=(3,3)):
  inputs = Input(shape=input_shape)
  
  # 入力層
  x = Conv2D(filters=first_filters, kernel_size=kernel_size, padding='same')(inputs)
  x = BatchNormalization()(x)
  x = ReLU()(x)
  x = MaxPool2D((2, 2))(x)
  x = Dropout(0.2)(x)
  
  # 畳み込み層
  for s in range(block_sets):
    filters =  first_filters * (2**s)
    
    for b in range(blocks):
      x = block_f(x, kernel_size, filters, block_layers)
      
    x = MaxPool2D((2, 2))(x)
    x = Dropout(0.2)(x)
  
  # 出力層
  x = Flatten()(x)
  x = Dropout(0.4)(x)
  x = Dense(100)(x)
  x = ReLU()(x)
  outputs = Dense(10, activation='softmax')(x)
  
  model = Model(input=inputs, output=outputs)
  
  return model

shortcut connection無しのブロック

layersの数だけ畳み込み層を追加する実装です。畳み込み後にBatchNormalizationとReLUを付加してます。

# shortcut connection無しのブロック
def plain_block(x, kernel_size, filters, layers):
  for l in range(layers):
    x = Conv2D(filters, kernel_size, padding='same')(x)
    x = BatchNormalization()(x)
    x = ReLU()(x)
    
  return x

shortcut connection有りのブロック (residual block)

shortcut connection有りのブロックも同様にlayers層の畳み込み層を構成しますが、入力信号をshortcut_xに保持しておき、最後のReLU関数の前にAddでdeep pathの出力信xと加算してます。

加算ではフィルタ数が一致する必要があるため、ズレている場合は1x1の畳み込み層を噛ませてフィルタ数を揃えてます

# shortcut path有りのブロック (residual block)
def residual_block(x, kernel_size, filters, layers=2):
  shortcut_x = x
  
  for l in range(layers):
    x = Conv2D(filters, kernel_size, padding='same')(x)
    x = BatchNormalization()(x)
    
    if l == layers-1:
      if K.int_shape(x) != K.int_shape(shortcut_x):
        shortcut_x = Conv2D(filters, (1, 1), padding='same')(shortcut_x)  # 1x1フィルタ
      
      x = Add()([x, shortcut_x])
      
    x = ReLU()(x)
    
  return x

検証

3つのモデルを生成・学習し、train lossを比較します。

引数block_fをthin/plainとresidualで切り替えていることに着目してください

# ハイパパラメータ
epochs = 200
batch_size = 50

# thinモデル
thin_model  = generate_model(input_shape, plain_block, blocks=1, block_sets=1)
thin_model.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.Adam(), metrics=['accuracy'])
thin_history = thin_model.fit(X_train, Y_train, epochs=epochs, batch_size=batch_size, validation_data=(X_test, Y_test))

# plainモデル
plain_model  = generate_model(input_shape, plain_block, blocks=3, block_sets=2)
plain_model.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.Adam(), metrics=['accuracy'])
plain_history = plain_model.fit(X_train, Y_train, epochs=epochs, batch_size=batch_size, validation_data=(X_test, Y_test))

# residualモデル
residual_model  = generate_model(input_shape, residual_block, blocks=3, block_sets=2)
residual_model.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.Adam(), metrics=['accuracy'])
residual_history = residual_model.fit(X_train, Y_train, epochs=epochs, batch_size=batch_size, validation_data=(X_test, Y_test))

CIFAR-10で実験

それではtrain lossを示します。今回の実験では、畳み込み層が浅く、精度・収束速度の両方においてResNetの持ち味は発揮できていないようです。

thinとplain/residualを比較すると、一段小さくなっており、畳み込み層を厚くすることで特徴抽出が改善できていることを確認できます
plainとresidualを比較すると、わずかですが一貫してplainの方が小さい値となってます

f:id:ohke:20190621163549p:plain

畳み込み層をさらに増やして実験 (13層 -> 25層)

上の結果を受けて、今度はplainとresidualの畳み込み層を13層から25層まで増やして学習させてみます。

1個のブロックセットが持つ畳み込み層を6層から12層に増やしています
- 1+12*2で、合計25層となります
それ以外は変更を加えていません

# plainモデル
plain_model  = generate_model(input_shape, plain_block, blocks=6, block_sets=2, first_filters=32)

# residualモデル
residual_model  = generate_model(input_shape, residual_block, blocks=6, block_sets=2, first_filters=32)

同じくtrain lossを比較すると、常にplainモデルよりもresidualモデルが下回っており、また収束も速くなってます。
一方でplainモデルは、12層バージョンよりもlossが大きい傾向にあり、層が増えると学習が難しくなることも伺えます。

f:id:ohke:20190621201147p:plain

検証データでのaccuracyで比較しても、plainモデルより3〜4%程度改善されています

13層のplainモデルと比較しても、約2%改善しています

f:id:ohke:20190621201655p:plain

まとめ

ResNetをKerasを使って実装し、CIFAR-10タスクでプレーンなCNNと比較しました。畳み込み層を深くすると、精度・収束速度の両方でResNetの方が良くなったことを確認しました。

Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun, Deep Residual Learning for Image Recognition (2015) (arXiv)↩

2019-06-15

numpyで無相関化・白色化する

Python 機械学習

最近通っているディープラーニングの講習会にて、BatchNormalizationの文脈でデータの白色化なるものについて触れましたので、「はじめてのパターン認識」を読みながらnumpyで実装してみます。

はじめてのパターン認識

作者: 平井有三
出版社/メーカー: 森北出版
発売日: 2012/07/31
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 7回
この商品を含むブログ (5件) を見る

白色化とは

白色化 (whitening)は $\vec{x}$ (d次元) の各特徴量を無相関化し、かつ、平均0・標準偏差1にすることです。

無相関化

最初に行う無相関化では、各特徴量の相関係数が0になるように線形変換します。

サンプル数n・次元数dの入力行列 ${\bf X}$ (n行×d列) から、平均ベクトル $\vec{\mu}$ (1×d) と共分散行列 ${\bf \Sigma}$ (d×d) を計算します。

${\bf \Sigma}$ の各要素 $\sigma_{i,j}$ は、i=jなら分散、i≠jなら共分散です

$$ \vec{\mu} = (\mu_1, ..., \mu_d) = (E[\vec{X_1}], ..., E[\vec{X_d}]) $$

$$ {\bf \Sigma} = E[(\vec{x} - \vec{\mu})(\vec{x} - \vec{\mu})^T] $$

共分散行列を固有値 $\lambda_i$ ・固有ベクトル $\vec{s}_i$ に分解し、固有ベクトルを列ベクトルとして結合した行列 ${\bf S}$ を得ます。

${\bf \Sigma} \vec{s_i} = \lambda_i \vec{s_i}$ の関係となります
固有値分解はnumpy.linalg.eigで計算できます

あとは ${\bf S}$ を使って変換後の $\vec{x'}$ を計算します。

${\bf S}$ は元の座標系から小昨夜靴方向へ回転させる回転行列となってます
固有値は対応する固有ベクトル方向の分散となります

$$ \vec{x'} = {\bf S}^T \vec{x} $$

これで各特徴量間の相関係数は0となります。

白色化

無相関化後に、平均0・標準偏差1にします。

固有値 $\lambda_i$ を対角化した行列 ${\bf \Lambda}$ を定義し、以下の計算を行います。

$$ \vec{x'} = {\bf \Lambda}^{-\frac{1}{2}} {\bf S}^{T} (\vec{x} - \vec{\mu}) $$

numpyでの実装

irisデータセットを使った白色化の実装は以下です。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

d = 2
n = 50

# データセットのロード
iris = load_iris()

# setosaの萼片の長さ・幅のみを使う
X = iris.data[iris.target == 0][:, :d].T  # (2, 50)
plt.scatter(X[0], X[1])  # 元の散布図

# 平均ベクトル
mu = np.mean(X, axis=1).reshape((d, 1))  # (2, 1)

# 共分散行列
Sigma = np.dot((X - mu), (X - mu).T) / n  # (2, 2)

# 無相関化
lam, S = np.linalg.eig(Sigma)  # 固有値・固有ベクトルに分解

X_dash = np.dot(S.T, X)  # (50, 2)
plt.scatter(X_dash[0], X_dash[1])  # 無相関化後の散布図

# 白色化
Lambda = np.diag(lam)  # 固有値の対角化

X_dash = np.dot(
    np.linalg.inv(np.sqrt(Lambda)),
    np.dot(S.T, (X - mu))
)
plt.scatter(X_dash[0], X_dash[1])  # 白色化後の散布図