ロジスティック回帰の自然勾配降下法

Introduction

機械学習の学習アルゴリズムではさまざまな数理最適化が活躍しています。今回はロジスティック回帰を例に、
・標準的な学習アルゴリズムである勾配降下法
・勾配降下法のconsを解決し得る自然勾配降下法
を紹介し、その比較対照実験を行いましょう。

なお勾配降下法や自然勾配降下法の実装はR言語によって行っています。環境はMac OS version 10.14, R version 3.5.0です。また、このノートは数値計算 Advent Calendar 2018, 14日目の記事になります。作成者の@ceptreeさんに感謝申し上げます。なお、急用により投稿が遅くなってしまったことをお詫び申し上げます。

デモデータの生成

ロジスティック回帰の学習に用いるデータを乱数で生成しておきましょう。

・ラベルy = 1のデータポイントを母平均(1, 1), 母分散Eで生成する。
・ラベルy = 0のデータポイントを母平均(-1, -1), 母分散Eで生成する。

として、それぞれ100個ずつサンプルしてきます。

# デモデータの生成
library(MASS)
each_size <- 100
X_1 <- mvrnorm(n = each_size, mu = c(1, 1), Sigma = diag(c(1, 1)))
X_0 <- mvrnorm(n = each_size, mu = c(-1, -1), Sigma = diag(c(1, 1)))
X <- rbind(X_1, X_0)    # 2次元の入力行列X
y <- rep(c(1, 0), each = each_size)    # 出力ベクトルt
plot(X, col = factor(y))    # 可視化

なお、あとで学習結果を比較するためにtrain test splitしておきましょう。

# train test split : train_size = 0.75
train_index <- sample(x = 1:(2*each_size), size = 0.75*2*each_size)
test_index <- setdiff(x = 1:(2*each_size), y = train_index)
X_train <- X[train_index, ]; X_test <- X[test_index, ]
y_train <- y[train_index]; y_test <- y[test_index]

ロジスティック回帰と勾配降下法による学習

ロジスティック回帰モデルは、以下のような確率の推定器のパラメータwについて、訓練データで交差エントロピー損失の平均値が最小になるような解を求めることで学習しています。

ここでσ(z)はsigmoid関数です。交差エントロピー損失の最小化問題は解析的に解くことができないため、実際には勾配降下法などの最適化アルゴリズムで解くのが一般的でした。

L(w)を交差エントロピー損失の訓練データ上での平均値とするとき、勾配降下法は次のように実装されます。

1. パラメータwの初期値を決める。
2. 更新式w <- w - eta * grad(L(w))を繰り返す。

勾配降下法では、事前にetaと損失の勾配grad(L(w))を準備しておく必要があります。etaは、事前に正の数でスケジュールするハイパーパラメータで、学習率と呼ばれています。交差エントロピー損失の勾配grad(L(w))の計算の実装では、合成関数の微分公式を用いるのが便利です。ここで、sigmoid関数の微分公式が役に立ちます。

σ'(z) = σ(z)(1-σ(z))

では、実際に実装を見てみましょう。

# 勾配降下法
# 事前に必要な関数を準備しておく。
sigmoid <- function(z){1/(1+exp(-z))}
d_sigmoid <- function(z){sigmoid(z)*(1-sigmoid(z))}
cross_entropy <- function(y, p){mean(-y*log(p)-(1-y)*log(1-p))}

# 各更新時の損失のhistoryを記録するオブジェクトを準備する。
loss_GD <- c()

# ここからが学習アルゴリズム
# Step 1 : 初期値と学習率ηを決める。
w <- c(0, 0)   # 偏回帰係数の初期値
eta <- 0.1    # 学習率
size <- length(y_train)    # 訓練データのサンプルサイズを取得しておく。
# Step 2 : 更新のiteration, 今回は10回更新する。
for(i in 1:10){
    # forward計算
    z <- X_train %*% w
    p <- sigmoid(z)
    loss_GD[i] <- cross_entropy(y_train, p)
    # 勾配の計算
    dp <- (p - y_train) / (size * p * (1-p))
    dz <- d_sigmoid(z)
    dw <- t(X_train) %*% (dz * dp)
    # wの更新
    w <- w - eta * dw
}
# Step 3 : 学習結果の確認
plot(x = 1:length(loss_GD), y = loss_GD, type = "l")
cross_entropy(y_test, sigmoid(X_test %*% w))    # テストデータでの損失

自然勾配降下法

勾配降下法は本来、目的関数の最小値問題を解くための数理最適化アルゴリズムです。一方で、ロジスティック回帰モデルには母集団モデルがあります。二項分布を用いて書くとBinom(1, σ(w^Tx))です。勾配降下法に、更新中の現状の解が真のモデルからどれだけ異なっていそうかという情報を反映させることができればより解に早く収束するのではないでしょうか。

これを実現したのが自然勾配降下法です。確率分布の離れ具合はKullback-Leibler divergenceによって測る方法があります。Kullback-Leibler divergenceの2次近似がFisher情報行列Fによる2次形式によって与えられることに注意すると、勾配降下法を以下のように修正すれば良いことがわかります。

1. パラメータwの初期値を決める。
2. 更新式w <- w - eta * F^{-1} grad(L(w))を繰り返す。

実装例は以下のようになります。

# 自然勾配降下法
# 事前に必要な関数を準備しておく。
sigmoid <- function(z){1/(1+exp(-z))}
d_sigmoid <- function(z){sigmoid(z)*(1-sigmoid(z))}
cross_entropy <- function(y, p){mean(-y*log(p)-(1-y)*log(1-p))}

# 各更新時の損失のhistoryを記録するオブジェクトを準備する。
loss_NGD <- c()

# ここからが学習アルゴリズム
# Step 1 : 初期値と学習率ηを決める。
w <- c(0, 0)    # 偏回帰係数の初期値
eta <- 0.1    # 学習率
size <- length(y_train)    # 訓練データのサンプルサイズを取得しておく。
# Step 2 : 更新のiteration, 今回は10回更新する。
for(i in 1:10){
    # forward計算
    z <- X_train %*% w
    p <- sigmoid(z)
    loss_NGD[i] <- cross_entropy(y_train, p)
    # 勾配の計算
    dp <- (p - y_train) / (size * p * (1-p))
    dz <- d_sigmoid(z)
    dw <- t(X_train) %*% (dz * dp)
    # Fisher情報行列の計算
    grad_logLL_z <- (y_train - p) / (p * (1-p)) * dz
    grad_logLL_w <- matrix(rep(grad_logLL_z, each = ncol(X_train)), nrow(X_train), ncol(X_train), byrow = TRUE) * X_train 
    F <- cov(grad_logLL_w)
    # wの更新
    w <- w - eta * solve(F) %*% dw
}
# Step 3 : 学習結果の確認
plot(x = 1:length(loss_NGD), y = loss_NGD, type = "l")
cross_entropy(y_test, sigmoid(X_test %*% w))    # テストデータでの損失

結果の比較

これは訓練データでの損失最小化の様子を更新ごとに記録したグラフです。

# 結果の可視化
library(ggplot2)
res <- data.frame(method = rep(c("Gradient Descent", "Natural Gradient Descent"), each = 10), loss = c(loss_GD, loss_NGD))
res$iteration <- c(1:10, 1:10)
ggplot(res, aes(x = iteration, y = loss, color = method)) + geom_line()

同じ初期値から更新を開始したのですが、自然勾配降下法は勾配降下法に比べ圧倒的なスピードで収束しています。テストデータでの損失も

・勾配降下法：0.390
・自然勾配降下法：0.083

とたった10回のiterationで桁違いの損失最小化を達成できていました。もはや雲泥の差です。確率分布の近さを測ることで勾配降下法の収束を改善させることが数値実験的に見て取れる結果になりました。やったね！

注意点

さて、自然勾配降下法のすばらしさを見てきましたが、ここでアルゴリズムにFisher乗法行列の逆行列が登場することには注意が必要です。まず空間計算量はパラメータのサイズdに対して2乗でかかってきます。また、一般に逆行列の計算には、例えばガウスの消去法などを用いるとしてO(d³)の時間計算量が要求されます。要は、パラメータサイズdが大きくなると、このままの実装では立ち行かなくなるということには注意を払う必要があります。（そのあたりの話はまた後日書きます。）

おわりに

深層学習の学習アルゴリズムにAdamというoptimizerがよく使われていることをご存知かもしれません。機械学習では学習アルゴリズムに限らず様々なところで数理最適化が活躍しています。今回の自然勾配降下法のように、統計的な要素が絡むことで独自の発展をしている部分もまた、機械学習にまつわる数理最適化のとても面白いところだと思っています。

機械学習と数理最適化の交点には本当に面白いトピックがたくさんあります。この記事をきっかけに、少しでも多くの数理最適化好きの皆さんに、機械学習をより楽しいと思っていただければ幸いです。