初歩からの機械学習：最急降下法による重回帰モデル～PythonとRでスクラッチから～

機械学習の教師あり学習の中でも、重回帰モデルはとても有名です。統計学でも有名なこのモデルですが、機械学習では、最急降下法というもっとも基本的かつ、重要なアルゴリズムを使ってパラメーターを求めることができます。

今回は、最急降下法を使って重回帰モデルのパラメーターを推定し、回帰問題を解いていきます。コードはPythonとRの両方を使って、スクラッチからコードを書いていきますので、機械学習アルゴリズムを0から理解したいという方は参考にしてみてください。（Rコードは最後にまとめて記述してあります。）

最急降下法とは？
最急降下法を使った単回帰モデル
本題の重回帰モデルのその前に
最急降下法による重回帰モデル
まとめ
参考書籍
Rコード

最急降下法とは？

最急降下法とは、ある関数を最小にするようなパラメーターを求めるためのアルゴリズムで、複雑な関数に対しても対応できます。とはいっても、いきなり複雑な問題では理解しにくいので、まずは、 $\,f(\theta) = (\theta-3)^{2}+5\,$ を最小化する $\,\theta\,$ を求めてみます。こんな感じの関数になります。

f:id:h-wadsworth02:20181229120721p:plain

計算するまでもなく、この関数の値が最小になる $\theta$ は3ですが、これを最急降下法を使って求めていきましょう。まずは事前準備として、微分を使って $\theta$ を求めてみます。最急降下法には微分の考えがとても重要なのです。

$\dfrac{d}{d\theta}f(\theta) = 2\theta-6$

傾きが０となるところが関数を最小にする値なので、

$2\theta- 6 = 0\\ \theta = 3$

となってこの関数を最小とする $\theta$ は3であると求めることができました。ここまではただの数学です。ここで、適当に $\theta$ の値決めて、それぞれの傾きを求めてみます。

最急降下法

みてみると $\theta$ の値が3から離れるほど傾きの絶対値は大きくなっています。この傾き利用して、少しずつ $\theta$ を更新しながら関数を最小とする値を求めるのが最急降下法のメカニズムです。傾きがもっとも急なところに降りていくという意味からこの名前がついています。例えば、マイナス側に大きく値がずれているときは、傾きは大きなマイナスになるため、 $\theta$ は大きくプラス側に更新されますし、プラス側に小さくずれている場合は、傾きは小さなプラスの値になるため、小さくマイナス方向に更新されます。

イメージとしてはボールを転がして、最終的には凹みの最下部にくる感じですね。最急降下法を式で表すとこんな感じになります。

$\theta = \theta-\alpha \dfrac{d}{d\theta}f(\theta)=\theta-\alpha(2\theta- 6)$

ここで $\alpha$ は学習率と呼ばれ、求めた傾きをどれぐらいパラメーターの更新に反映させるかを表し、1より小さい値にします（1以上だと発散してしまうため）。実際に値を入れて確かめてみましょう。 $\alpha$ の値は0.01とし、 $\theta$ の初期値は6とします。

$6-0.01 \times 6 = 5.94$

少しだけ $\theta$ の値がマイナス側に更新されて小さくなりました。もう一度更新してみましょう。

$6-0.01 \times (2\times5.94 -6) = 5.8812$

さらに小さくなりましたね。これを繰り返していきます。Pythonで実装してみましょう。

import numpy as np
import matplotlib.pyplot as plt
% matplotlib inline


alpha = 0.01　#学習率
itera = 1000　#繰り返し回数(パラメータ更新回数)
theta_lis = np.zeros(itera)　#更新ごとにパラメータを格納するためのもの

theta = 6 #thetaの初期値

#繰り返し処理
for i in range(itera):
    theta_lis[i] = theta
    theta = theta - alpha*(2*theta-6)

plt.plot(theta_lis)
plt.ylabel("theta")
plt.xlabel("iteration")
plt.show()

最急降下法

1000回値を更新しました。大体200回ぐらい繰り返したところで、 $\theta$ の値は3に近づいているようです。値を確認してみると

theta
>3.0000000050489004

ほぼ3になっていますね。

最急降下法を使った単回帰モデル

さて、お次は視覚的にも理解しやすい単回帰モデルを例に最急降下法を使ってパラメーターを求めてみましょう。単回帰モデルのパラメーター（切片と傾き）は最小二乗法で求められます。

www.medi-08-data-06.work

最小二乗法とは、それぞれの点からの距離がもっとも近くなるところを通るように回帰直線のパラメーターを求める方法なのですが、直線からそれぞれの点の距離を二乗して全て足した値の平均値（平均二乗誤差:Mean Square Error又の名をコスト）を最小にするように切片と傾きを求めます。つまり、

$cost= \dfrac{1}{m} \sum (y - (\theta_{0} + \theta_{1}x))^{2}$

$mはサンプルサイズ$

を最小にするパラメーター（ここでは切片： $\theta_{0}$ と傾き： $\theta_{1}$ ）を最急降下法で求めることになります。これを機械学習の世界ではコスト関数を最小にするパラメーターを求めると呼ぶこともあります。どうするかというと、上記のコスト関数をそれぞれのパラメーターで偏微分して、それぞれの傾きを使って値を更新していきます。まずは以下のようにコスト関数を定義します。

$cost = f(\theta_{0},\theta_{1})= \dfrac{1}{2m} \sum ( (\theta_{0} + \theta_{1}x)-y)^{2}$

コスト関数は微分の便宜上、1/mではなく、1/2mを使うことが一般的です。これをそれぞれのパラメーターで偏微分してみます。偏微分は難しそうに聞こえて、普通の微分となんら変わりません。偏微分する変数以外を全て定数として扱うだけのなので、合成関数の微分を使って（合成関数の微分も難しくないので是非調べてみてください。）

$\theta_{0}$ で偏微分

$\dfrac{\partial}{\partial \theta_{0}}f(\theta_{0},\theta_{1}) = \dfrac{1}{m} \sum ((\theta_{0} + \theta_{1}x)-y) \times 1$

$\theta_{1}$ で偏微分

$\dfrac{\partial}{\partial \theta_{1}}f(\theta_{1},\theta_{1}) = \dfrac{1}{m} \sum ((\theta_{0} + \theta_{1}x)-y) \times x$

となります。通常であれば、それぞれを傾き０として求めるのですが、ここに最急降下法を導入して、

$\theta_{0} = \theta_{0} - \alpha \dfrac{\partial}{\partial \theta_{0}}f(\theta_{0},\theta_{1}) = \theta_{0} - \alpha\dfrac{1}{m} \sum ((\theta_{0} + \theta_{1}x)-y)$

$\theta_{1} = \theta_{1} - \alpha \dfrac{\partial}{\partial \theta_{1}}f(\theta_{0},\theta_{1}) = \theta_{1} -\alpha \dfrac{1}{m} \sum ( (\theta_{0} + \theta_{1}x)-y)x$

とすれば、先ほど同様にコスト関数を最小とするパラメーターを求めることができますね。Pythonでこれを実装してみましょう。データは有名なボストンデータセットを使って、部屋の広さと家賃の関係を単回帰モデルで実装します。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
% matplotlib inline

#ボストンデータセットの読込
boston = load_boston()
data = pd.DataFrame(boston.data,columns=boston.feature_names)

とりあえず、家賃と部屋の広さのグラフにしてみます。横軸が部屋の広さ、縦軸が家賃です。

plt.grid(which='major',color='gray',linestyle='-')
plt.scatter(data["RM"],boston.target)
plt.xlabel("Room")
plt.ylabel("Medv")
plt.show()

最急降下法

みたところ回帰直線が引けそうです。いよいよ最急降下法を実装していきます。

X = data["RM"]
y = boston.target

alpha = 0.04 
itera = 10000
cost = np.zeros(itera)
theta0 = 1　#theta0の初期値
theta1 = 1　#theta1の初期値
m = len(y)

for i in range(itera):
    cost[i] = (1/(2*m))*np.sum(np.square(theta0+theta1*X-y))
    temp0 =theta0 - alpha*(1/(m))*np.sum(theta0+theta1*X-y)
    temp1 =theta1 - alpha*(1/(m))*np.sum((theta0+theta1*X-y)*X)
    
    theta0 = temp0
    theta1 = temp1

plt.plot(cost)
plt.xlabel("Cost")
plt.ylabel("Iteration")
plt.show()

最急降下法

今回の学習率は0.04にしました。ここを色々変えて結果がどうなるかも試してみてください。また、それれぞれの繰り返しごとにコストも格納していきます。thetaをいったんtempという変数に格納するのは、theta0もtheta1も更新される前の値を使うためです。

さて、今回は10000回パラメータを更新しました。コストは前半に大きく減少し、徐々に減っているように見えます。パラメーターを確認してみます。

print(theta0)
print(theta1)

>-34.377040288363766
>9.055956640093719

切片が-34、傾きが9の回帰直線パラメータが求まりました。グラフに加えてみます。

plt.grid(which='major',color='gray',linestyle='-')
plt.scatter(data["RM"],boston.target)
plt.plot(data["RM"],theta0+theta1*data["RM"],color="red")
plt.xlabel("Room")
plt.ylabel("Medv")
plt.show()

最急降下法

おお！もっともらしい回帰直線が引けました。これが最急降下法を使った単回帰モデルパラメーターの求め方です。

本題の重回帰モデルのその前に

さて、本題の重回帰を最急降下法を使って求めていく前に、ベクトル化の概念をご説明します。単回帰モデルでは、Xもy も一次元でしたが、重回帰モデルはXの変数がたくさんあるため、 $\theta_{0} + \theta_{1}x1+\theta_{2}x2\dots$ と書くのも面倒ですし、パラメーターを一つ一つ更新するのも手間です。そこで、変数がいくつあっても一つの文字で表してしまおうというのがベクトル化です。一つ一つパラメーターを扱う必要がないため計算もとても速くなります。まず、通常の重回帰モデルは以下のようになります。

$y=\theta_{0} + \theta_{1}x1+\theta_{2}x2$

これをベクトル化するとどうなるかというと

f:id:h-wadsworth02:20190614212457p:plain

と表すことが出来ます。ここでXは1行3列の行列、 $\theta$ は3行1列の行列です。さらにこれをサンプル（m）だけならべればよいので、

f:id:h-wadsworth02:20190614212807p:plain

となって、同様に同じ式で表すことができます。上矢印は縦に並ぶベクトルであることを表しています。yはm行1列の行列で、Xはm行3列の行列になります。ベクトル化は慣れが必要かと思いますが、よくわからんという方は以下の本がおすすめです。

統計学が最強の学問である[数学編]――データ分析と機械学習のための新しい教科書

作者:西内啓
出版社/メーカー: ダイヤモンド社
発売日: 2017/12/21
メディア: 単行本（ソフトカバー）

このベクトル化を使って、最急降下法の式も書き直してみましょう。

最急降下法

$\boldsymbol{X}^{T}$ は転置を表し、m行3列の行列を3行m列に変換しています。こうすることで、最終的に $\boldsymbol{\theta}$ が3行1列の行列となるため、一気にパラメーターを更新することが出来るのです。

最急降下法による重回帰モデル

いよいよ今までの知識を総動員して、最急降下法による重回帰モデルのパラメーターを求めてみましょう。先ほど同様にボストンデータセットを使い、部屋の広さに加えて、低所得者の割合を表す変数の2変数による重回帰モデルを作ります。まずは、Xとｙを定義しておきましょう。

import statsmodels.api as sm

X = data[["RM",'LSTAT']]
X = sm.add_constant(X) #最初の列に1を加える
y = boston.target

Xはm行2列の行列ですが、切片の項を加える必要があるため最初の列に1の列を加えます。次に $\boldsymbol\theta$ の初期値を定義します。 $\boldsymbol\theta$ は切片と変数を合わせて3行1列の行列になります。

init_theta = np.zeros(3).reshape(3,)　#reshapeで3行1列にする。

パラメーターの更新を行うコードを書いてみます。

alpha = 0.001
m = len(y)
init_theta - alpha*(1/m)*np.dot(X.T,(np.dot(X,init_theta)-y))

>array([0.02253281, 0.14609502, 0.23675723])

numpyのdotは行列の掛け算を行う関数です。これを実行すると更新後の3行1列の $\theta$ が返ってきます。

さて、これらの処理をまとめて、コストの格納とパラメーター更新を行う関数を定義しておきます。

def grad_desc(theta,alpha,itera,X,y):
    m = len(y)
    cost = np.zeros(itera)
    
    for i in range(itera):
        cost[i] = (1/(2*m))*np.sum(np.square(np.dot(X,theta)-y))
        theta = theta - alpha*(1/m)*np.dot(X.T,(np.dot(X,theta)-y))   
    return [cost,theta]

これで準備が整いました！いよいよパラメーターを求めていきます。今回は1000回の更新を行います。

alpha = 0.001
itera = 1000
cost , theta = grad_desc(init_theta,alpha,itera,X,y)

plt.plot(cost)
plt.xlabel("Cost")
plt.ylabel("Iteration")
plt.show()

最急降下法

最終的なcostとthetaの値を確認してみると

print(theta)
print(cost[itera-1])

>[ 0.64799681  4.81706945 -0.66204958]
>15.268362109644587

$家賃=0.64 + 4.81\times 部屋の広さ-0.66\times 低所得者の割合$

という関係性にあるようです。試しにsklearnを使って重回帰を行った結果と比較してみましょう。

from sklearn import linear_model

X = data.loc[:,["RM","LSTAT"]]
lreg = linear_model.LinearRegression()
lreg.fit(X,y)

print(lreg.intercept_)
print(lreg.coef_)

>-1.3582728118745315
>[ 5.09478798 -0.64235833]

$家賃=-1.35 + 5.09\times 部屋の広さ-0.64\times 低所得者の割合$

若干値は違いますが、大枠はあっているようです。ついでにコストも比較してみましょう。

#結果をthetaに代入
theta = np.array([lreg.intercept_,lreg.coef_[0],lreg.coef_[1]]).reshape(3,)
X = data[["RM",'LSTAT']]
X = sm.add_constant(X)

#コスト関数で求める
(1/(2*m))*np.sum(np.square(np.dot(X,theta)-y))
>15.256234388649734

こちらの方が若干性能が良いようです。最急降下法では学習率や更新回数を調整することで、性能を上げることができるので、いろいろやってみてください。

補足：最急降下法の使い道

最急降下法を使わなくても重回帰モデルのパラメーターを求めることが出来るのですが、なぜ最急降下法のアルゴリズムが重要なのでしょうか？実は簡単な線形モデルであれば、正規方程式という行列演算を用いてそれぞれのパラメーターを求めています。すべての偏微分方程式をイコール0として求めているイメージです。今回のような簡単なモデルであれば、正規方程式によるパラメーター推定の方が効率がよいです。しかし、正規方程式は変数の数が膨大になると、解をもとめることが困難になります。そのため、ビックデータを前提とする機械学習モデルでは、最急降下法を基本とする様々な勾配アルゴリズムと呼ばれるものを駆使して、パラメーターを求めるのです。

まとめ

今回は機械学習モデルの中でも、もっとも基本的な重回帰モデルを、勾配アルゴリズムの中でも最も基本的な最急降下法を使って実装しました。今回の内容はシンプルではありますが、教師あり学習の様々なモデルのベースとなるところですので、理解を深めておきたいところです。

※本記事は筆者が個人的に学んだことをまとめた記事なります。所属する組織の意見・見解とは無関係です。また、数学の記法や詳細な理論、用語等で誤りがあった際はご指摘頂けると幸いです。

参考書籍

ベクトルや行列の基礎から、最急降下法などのアルゴリズム、ニューラルネットワークまで機械学習を学ぶための前提知識が網羅されています。本書の立ち位置は機械学習を学ぶための知識を学ぶための本とされ、まさにその通りの内容です。

統計学が最強の学問である[数学編]――データ分析と機械学習のための新しい教科書

作者:西内啓
出版社/メーカー: ダイヤモンド社
発売日: 2017/12/21
メディア: 単行本（ソフトカバー）

www.medi-08-data-06.work

Rコード

library(tidyverse)
library(MASS)

#最急降下法による二次関数の最小化
alpha = 0.01　#学習率
itera = 1000　#繰り返し回数(パラメータ更新回数)
theta_lis = numeric(itera)　#更新ごとにパラメータを格納するためのもの

theta = 6
for (i in 1:itera){
  theta_lis[i] = theta
  theta = theta - alpha*(2*theta-6)
}
    
plot(theta_lis,xlab = "Iteration",ylab="Theta")


#最急降下法による単回帰
X <- Boston$rm
y <- Boston$medv

alpha <- 0.04 
itera <- 10000
cost <- numeric(itera)
theta0 <- 1　#theta0の初期値
theta1 <- 1　#theta1の初期値
m <- length(y)

for (i in 1:itera){
  cost[i] <- (1/(2*m))*sum((theta0+theta1*X-y)^2)
  temp0 =theta0 - alpha*(1/(m))*sum(theta0+theta1*X-y)
  temp1 =theta1 - alpha*(1/(m))*sum((theta0+theta1*X-y)*X)
  
  theta0 = temp0
  theta1 = temp1
  
}
    

plot(cost,xlab="Cost",ylab="Iteration")
#結果の可視化
plot(Boston$rm,Boston$medv)
abline(a=theta0,b=theta1,col="red")


#最急降下法による重回帰

X <- matrix(c(rep(1,nrow(Boston)),Boston$rm,Boston$lstat),ncol=3)
y <- Boston$medv
init_theta <- rep(0,3)

#関数の定義
grad_desc <- function(theta,alpha,itera,X,y){
  m <- length(y)
  cost <- numeric(iter)
  
  for(i in 1:itera){
    cost[i] <-  (1/(2*m))*sum((y-X %*% theta)^2)
    theta <- theta - alph*(1/m)*(t(X)%*%(X %*% theta-y))
  }

#結果はリストで返す
  return(list(cost,theta))
}


alph <- 0.001
itera <- 1000
res <- grad_desc(init_theta,alph,itera,X,y)

#結果の格納
cost <- res[[1]]
theta <- res[[2]]

#glmによる重回帰モデル
fit <- glm(data=Boston,medv~rm+lstat)

#コストの計算
(1/(2*length(y)))*sum((fit$residuals)^2)