機械学習モデルの数学的記述と計算（その１）：バイナリ学習

機械学習は、データからそれらが生成されている法則性を学んで、予測や分類に役立てるものである。
ここでは、スタンフォード大学のAndrew Ng教授の講義ノートをベースに、機械学習の数学的バックグラウンドを整理する。

教師ありバイナリ学習のデータの表現

教師ありバイナリ学習に用いられるデータを、特徴量とターゲットのペアとして考える。これを、 $(x,y)$ により表現する。ここで、 $x \in \mathbb{R}^{n_{x}}$ は特徴量のベクトルを表し、 $y \in \{0,1\}$ はターゲット変数を表す。
m個の学習用データの組は以下で表される

$(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),\cdots,(x^{(m)},y^{(m)})$

特徴量とターゲット変数をそれぞれまとめて以下のように記述する

$x=(x^{(1)},x^{(2)},\cdots,x^{(m)}),\ \ x \in \mathbb{R}^{n_{x} \times m}$

$y=(y^{(1)},y^{(2)},\cdots,y^{(m)}),\ \ y \in \mathbb{R}^{1 \times m}$

ロジスティック回帰モデル
$x \in \mathbb{R}^{n_{x}}$ が与えられたとき、以下の式により $\hat{y}=P(y=1|X),0 \le \hat{y} \le 1$ を得たいとする。
（パラメータとアウトプット）

$w \in \mathbb{R}^{n_{x}},\ \ \ b \in \mathbb{R}$

$\hat{y}=\sigma(w ^{T} x+b)$

ここで、

$\sigma (z)=\frac{1}{1+\exp (-z)}$

ロジスティック回帰モデルの損失関数
上記の定式化で、 $\{(x^{(1)},y^{(1)}),\cdots,(x^{(m)},y^{(m)})\}$ が与えられたとき、 $\hat{y}^{(i)} \approx y^{(i)}$ を得たいとする。
損失（誤差）関数を以下のような「交差エントロピー誤差関数」として定義する。

$J(w, b)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}(\hat{y},y)=-\frac{1}{m}\sum_{i=1}^{m}(y^{(i)}\log \hat{y}^{(i)}+(1-y^{(i)})\log (1-\hat{y}^{(i)}))$

最急降下法
上の式の $J(w,b)$ を最小化する $w,b$ を見つけたいとする。例えば $J$ を $w$ の関数としてみた場合、 $J$ の勾配を計算して、勾配が急な方向に沿って進む。つまり以下を反復計算する

$w:=w-\alpha \frac{dJ(w)}{dw}$

ここで $\alpha$ は学習率(learning rate)である。

計算グラフ
機械学習では以下のような計算グラフを考える。そして、その微分を逆伝搬で考える。これにより計算を局所化し、計算効率が高まる。
ここでは、 $J(a,b,c)=3(a+bc)$ を例として考える。
f:id:nicjps230:20210711130639p:plain

上のロジスティック回帰モデルの計算グラフとその微分（概要のみ）は以下のようになる。
f:id:nicjps230:20210716140912p:plain

m個のデータ例でのロジスティック回帰モデル
損失関数の式を再掲する

$J(w, b)=\frac{1}{m}\sum_{i=1}^{m}\mathcal{L}(a^{(i)},y^{(i)})$

ここで、

$a^{(i)}=\hat{y}^{(i)}=\sigma (w^{T}x^{(i)}+b)$

いま、パラメータが $w_{1},w_{2},b$ のみであるとする。上記の損失関数の、パラメータ $w_{1}$ に関する傾きを求めるためには、以下を計算する。

$\frac{\partial }{\partial w_{1}}J(w,b)=\frac{1}{m}\sum_{i=1}^{m}\frac{\partial }{\partial w_{1}}\mathcal{L}(a^{(i)},y^{(i)})$

上記のように微分を逆伝搬で考える場合、各ノードの変数が1単位変化した時に最終出力変数がどのくらい変化するか、すなわち

$\frac{\partial (final output var)}{\partial (var)}=dvar$

を考えることになる。例えば、 $w_{1}$ についての勾配は、上の図の微分の逆伝搬を伝っていき、

$dw_{1}=\frac{\partial \mathcal{L}}{\partial w_{1}}=\frac{\partial \mathcal{L}}{\partial a}\cdot \frac{\partial a}{\partial z}\cdot x_{1}$

と求められる。いま、

$\mathcal{L}(a,y)=-(y\log a+(1-y)\log (1-a))$

の $a$ についての「勾配」を求めると、

$da=\frac{\partial \mathcal{L}(a,y)}{\partial a}=-\frac{y}{a}+\frac{1-y}{1-a}$

となる。
また、

$a=\sigma (z)=\frac{1}{1+\exp (-z)}$

を $z$ について解くと、

$z=-\log (\frac{1}{a}-1)$

となり、逆関数の微分法より、

$\frac{\partial a}{\partial z}=\frac{1}{\frac{\partial z}{\partial a}}=a(1-a)$

となる。それゆえ、

$dz=\frac{\partial \mathcal{L}}{\partial z}=\frac{\partial \mathcal{L}}{\partial a} \cdot \frac{\partial a}{\partial z}=(-\frac{y}{a}+\frac{1-y}{1-a}) \cdot a(1-a)=a-y$

となる。
さらに、

$z=w_{1}x_{1}+w_{2}x_{2}+b$

より、

$dw_{1}=\frac{\partial \mathcal{L}}{\partial w_{1}}=\frac{\partial \mathcal{L}}{\partial z} \cdot \frac{\partial z}{\partial w_{1}}=dz \cdot x_{1}$

となる。同様に、 $dw_{2}=x_{2} \cdot dz$ , 　 $db=dz$ となる。

アルゴリズムのpythonによる実装
上記のアルゴリズムでロジスティック回帰の係数を求める計算をpythonで実装した。データはKaggle上にあるロジスティック回帰用のデータを用いた。ここでの留意点は、特徴量の線型結合の絶対値が大きくなる（数十程度）と、シグモイド関数が極値に振れてしまい、対数が計算できなくなることである。そのため、ここではあらかじめダミー変数以外の変数を平均0、標準偏差1で基準化した。

# ライブラリのインポート
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib # matplotlib日本語
# データの読み込み
df01 = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Social_Network_Ads.csv')
# 性別に関するダミー変数の作成
df02 = pd.get_dummies(df01, drop_first=True, columns=['Gender'])
# ダミー変数以外の特徴量を標準化する
df02 = df02.assign(salary_std = (df02['EstimatedSalary']-df02['EstimatedSalary'].mean())/df02['EstimatedSalary'].std())
df02 = df02.assign(age_std = (df02['Age']-df02['Age'].mean())/df02['Age'].std())
# 分析に用いるデータをNumpy配列に書き出す
features_list = ['salary_std','Gender_Male','age_std']
X_train = df02[features_list].to_numpy()
y_train = df02['Purchased'].to_numpy()
# バイナリ学習のアルゴリズム
w1 = 1; w2 = 1; w3 = 1; b = 1 # 回帰係数の初期化
alpha = 1 # 学習率
m = len(y_train)
z = np.zeros(m); a=np.zeros(m); dz=np.zeros(m)
for j in range(250): # 指定した回数の反復計算
  J = 0; dw1 = 0; dw2 = 0; dw3 = 0; db = 0 # 損失関数と勾配の初期化
  for i in range(m): # データ行についての処理
    z[i] = w1 * X_train[i,0] + w2 * X_train[i,1] + w3 * X_train[i,2] + b
    a[i] = 1 / (1+np.exp(-z[i]))
    J += -(y_train[i] * np.log(a[i]) + (1-y_train[i]) * np.log(1-a[i]))
    dz[i] = a[i] - y_train[i]
    dw1 += X_train[i,0] * dz[i]
    dw2 += X_train[i,1] * dz[i]
    dw3 += X_train[i,2] * dz[i]
    db += dz[i]
  J /= m; dw1 /= m; dw2 /= m; dw3 /= m; db /= m
  w1 -= alpha * dw1
  w2 -= alpha * dw2
  w3 -= alpha * dw3
  b -= alpha * db
  if j % 20 == 0:
    print('j={}, J={:.3g}, w1={:.3g}, w2={:.3g}, w3={:.3g},b={:.3g}'.format(j,J,w1,w2,w3,b))