Rで勾配降下(上昇)法を行う

UPDATE: 2022-08-30 13:37:26

はじめに

このノートではRを使った勾配降下法(最急降下法)の方法についてまとめておく。

勾配降下法とは

勾配降下法(Gradient Descent)は関数の傾き(1階微分)を利用して、目的関数の最小値を探索するアルゴリズム。下記の更新式を用いて行われる。\(\eta\)は学習率であり、係数をかけることで学習速度を早めることもできる。

\[ x_{n+1} = x_{n} - \eta f^{\prime}(x_{n}) \]

1次関数で考えると、この式は、ある点で微分した際に、点の傾きがマイナスであれば右に動かし、傾きがプラスなら左に動かせるようにしている。例えば、\(\eta=1\)で\(f(x)=x^{2}\)を\(x=-2\)で微分すると、\(f^{\prime}(-2)=2x=2(-2)=-4\)となるため、更新式では\(x_{n} - \eta(-4)=-2 + 4 = +2\)となり、右に動くことになる。反対に、\(x=2\)で微分すると、\(f^{\prime}(2)=2x=2(2)=4\)となるため、更新式では\(x_{n} - \eta(4)=2 - 4 = -2\)となって、左に動く。

以前ノートにまとめたニュートン法が機械学習であまり使われないのかは、下記が詳しい。

• なぜニュートン法は機械学習であまり使われないのですか？

勾配降下法をRで実装していく。今回は、可視化用の更新履歴データを保存している。

f <- function(x) 1.2 * (x-2)^2 + 3.2
iter <- 200
eta <- 0.01
h <- 0.01
x <- 10

xtrace <- ftrace <- vector('numeric', length = iter)
for (i in seq_along(1:iter)) {
  dx <- ((f(x + h) - f(x)) / h)
  x <- x - eta * dx
  xtrace[i] <- x
  ftrace[i] <- f(x)
  if (i %% 10 == 0) print(sprintf("%d times: (x=%2.3f, f(x)=%2.3f)", i, x, f(x)))
}

## [1] "10 times: (x=8.274, f(x)=50.429)"
## [1] "20 times: (x=6.919, f(x)=32.241)"
## [1] "30 times: (x=5.857, f(x)=21.055)"
## [1] "40 times: (x=5.024, f(x)=14.176)"
## [1] "50 times: (x=4.371, f(x)=9.946)"
## [1] "60 times: (x=3.859, f(x)=7.345)"
## [1] "70 times: (x=3.457, f(x)=5.746)"
## [1] "80 times: (x=3.141, f(x)=4.763)"
## [1] "90 times: (x=2.894, f(x)=4.159)"
## [1] "100 times: (x=2.700, f(x)=3.788)"
## [1] "110 times: (x=2.548, f(x)=3.561)"
## [1] "120 times: (x=2.429, f(x)=3.421)"
## [1] "130 times: (x=2.335, f(x)=3.335)"
## [1] "140 times: (x=2.262, f(x)=3.282)"
## [1] "150 times: (x=2.204, f(x)=3.250)"
## [1] "160 times: (x=2.159, f(x)=3.230)"
## [1] "170 times: (x=2.124, f(x)=3.218)"
## [1] "180 times: (x=2.096, f(x)=3.211)"
## [1] "190 times: (x=2.074, f(x)=3.207)"
## [1] "200 times: (x=2.057, f(x)=3.204)"

この結果を可視化すると、うまく関数の最小値を探索できていそうです。

df <- data.frame(
  i  = 1:iter,
  x  = xtrace,
  fx = ftrace
)

x_base <- seq(-3,15,0.1)
y_base <- f(x_base)
df_base <- data.frame(x_base, y_base)

library(ggplot2)
ggplot() + 
  geom_line(data = df_base, aes(x_base, y_base)) + 
  geom_point(data = df, aes(x, fx), col = 'red', size = 2) +
  geom_hline(yintercept = df[max(df$i),'fx']) +
  scale_x_continuous(breaks = seq(-3, 15, 1)) +
  scale_y_continuous(breaks = seq(0, 100, 10), limits = c(0,100)) +
  labs(title = 'f(x) = 1.2 * (x-2)^2 + 3.2', x = 'x', y = 'y') + 
  theme_classic()

線形回帰と勾配降下法

勾配降下法を用いて、線形回帰のパラメタを探索してみる。線形回帰では2乗誤差を目的関数として、この関数を最小化できれば良いので、下記の関数を最小化する。

\[ s = \sum(y_{i} - (\alpha + \beta x_{i}))^{2} = \sum(y_{i} - \alpha - \beta x_{i})^{2} \]

サンプルデータを擬似的に生成し、得られるべきパラメタを確認しておく。

set.seed(1989)
x <- rnorm(100,0,1)
y <- 5 + 2.5*x + rnorm(100, 10, 1)
lm(y ~ x)

## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      14.930        2.493

\(\alpha, \beta\)ともに初期値は1として、50回の繰り返しを行う。ここでは関数\(s\)を偏微分して得られた式を予め使うのではなく、目的関数をそのまま利用し、微分しながら更新する。

# 目的関数
f <- function(alpha, beta) sum((y - alpha-beta*x)^2)

# 初期設定
beta <- 1
alpha <- 1
iter <- 50
eta <- 0.001
h <- 0.01

for(i in seq_along(1:iter)){
  # 更新前のalpha, beta
  params <- c(alpha, beta)
  
  # alphaを更新
  dalpha <- (f(params[1] + h, params[2]) - f(params[1], params[2])) / h
  alpha <- alpha - eta * dalpha
  # betaを更新
  dbeta <- (f(params[1], params[2] + h) - f(params[1], params[2])) / h
  beta  <-  beta - eta * dbeta
  
  if (i %% 5 == 0) print(sprintf("%d times: (alpha=%2.3f, beta=%2.3f)", i, alpha, beta))
}

## [1] "5 times: (alpha=10.288, beta=1.558)"
## [1] "10 times: (alpha=13.366, beta=2.048)"
## [1] "15 times: (alpha=14.396, beta=2.302)"
## [1] "20 times: (alpha=14.744, beta=2.414)"
## [1] "25 times: (alpha=14.862, beta=2.459)"
## [1] "30 times: (alpha=14.903, beta=2.477)"
## [1] "35 times: (alpha=14.917, beta=2.484)"
## [1] "40 times: (alpha=14.922, beta=2.486)"
## [1] "45 times: (alpha=14.924, beta=2.487)"
## [1] "50 times: (alpha=14.924, beta=2.487)"

最小二乗法で正規方程式を利用するlm()の結果と似たような結果が得られている。

ポアソン分布と勾配上昇法

関数の形によっては、勾配を上昇する必要がある。その例としてポアソン分布の尤度を取り上げる。ポアソン分布の尤度関数\(L\)は、

\[ L = \prod {\frac{\lambda^{x_{i}}}{x_{i}!}e^{-\mu}} \]

であり、対数尤度関数\(log L\)に直すと、

\[ log L = \sum x_{i} log(\lambda) - n\lambda - \sum log(x_{i}!) \]

となる。これを勾配上昇法で最大となる点を探索する。ポアソン分布の場合、最尤推定量はmean()で計算できるので、探索結果として得られだろう値を確認しておく。

# Data From 
set.seed(1989)
x <- rpois(n = 100, lambda = 3.5)
mean(x)

## [1] 3.4

さきほどの対数尤度関数をそのまま書いて、勾配上昇法を実行すると、最尤推定値とほとんど同じ値が得られている。

loglik <- function(lam) sum(x)*log(lam) - length(x)*lam - sum(log(factorial(x)))
lambda <- 1
iter <- 500
eta <- 0.001
h <- 0.01

for(i in seq_along(1:iter)){
  lambda <- lambda + eta * (loglik(lambda + h) - loglik(lambda)) / h

  if (i %% 50 == 0) print(sprintf("%d times: (loglik=%2.3f, lambda=%2.3f)", i, loglik(lambda), lambda))
}

## [1] "50 times: (loglik=-204.702, lambda=3.114)"
## [1] "100 times: (loglik=-203.487, lambda=3.336)"
## [1] "150 times: (loglik=-203.430, lambda=3.382)"
## [1] "200 times: (loglik=-203.426, lambda=3.392)"
## [1] "250 times: (loglik=-203.426, lambda=3.394)"
## [1] "300 times: (loglik=-203.426, lambda=3.395)"
## [1] "350 times: (loglik=-203.426, lambda=3.395)"
## [1] "400 times: (loglik=-203.426, lambda=3.395)"
## [1] "450 times: (loglik=-203.426, lambda=3.395)"
## [1] "500 times: (loglik=-203.426, lambda=3.395)"

関数の最大値、最小値を探索できるoptimize()でも同じような結果が得られている。

optimize(loglik, c(0, 10), maximum = TRUE)

## $maximum
## [1] 3.399999
## 
## $objective
## [1] -203.4255

正規分布と勾配上昇法

次の例として正規分布を取り上げる。正規分布は、

\[ f(x; \mu, \sigma^2)=\frac{1}{\sqrt{2 \pi \sigma^2}}exp \left[-\frac{(x-\mu)^2}{2 \sigma^2}\right] \]

であり、正規分布の対数尤度関数および偏微分の値は、

\[ \begin{eqnarray} log L &=& - \frac{n}{2}log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum (x_{i}-\mu)^2 \\ \frac{ \partial L(\mu, \sigma^2) }{ \partial \mu } &=& - \frac{1}{\sigma^2} \sum (x_{i}-\mu) = 0 \\ \frac{ \partial L(\mu, \sigma^2) }{ \partial \sigma^2 } &=& - \frac{n}{2\sigma^2} + \frac{1}{2\sigma^4}\sum (x_{i}-\mu)^2 = 0 \\ \end{eqnarray} \]

である。これを勾配上昇法で探索することで、パラメタを求められる。Rで実装していく前に、サンプルデータを擬似的に生成し、得られるべきパラメタを確認しておく。

# mu = 2, sigma2 = 9
set.seed(1989)
x <- rnorm(100, mean = 2, sd = sqrt(9))
list(
  mean = mean(x),
  sigma2 = var(x)
  )

## $mean
## [1] 1.735709
## 
## $sigma2
## [1] 9.703747

まずは、偏微分の結果を利用して更新していく。下記では\(\sigma\)を更新しているので、分散\(\sigma^2\)は\(\sigma\)の結果を2乗する必要がある

iter <- 500
h <- 0.01
eta <- 0.01
mu <- 1
sigma <- 1
n <- length(x)

# 下記ではσを更新しているので、分散σ^2はsigmaの結果を2乗する必要がある
for (i in seq_along(1:iter)) {
  mu <- mu + eta * sum(x - mu)/(sigma^2)
  sigma <- sigma + eta * (-n/(2*sigma^2)+sum((x - mu)^2)/(2*sigma^4))
  
  if (i %% 50 == 0) print(sprintf("%d times: (mu=%2.3f, sigma2=%2.3f, sigma=%2.3f)", i, mu, sigma^2, sigma))

}

## [1] "50 times: (mu=1.736, sigma2=22.110, sigma=4.702)"
## [1] "100 times: (mu=1.736, sigma2=16.454, sigma=4.056)"
## [1] "150 times: (mu=1.736, sigma2=12.247, sigma=3.500)"
## [1] "200 times: (mu=1.736, sigma2=10.274, sigma=3.205)"
## [1] "250 times: (mu=1.736, sigma2=9.739, sigma=3.121)"
## [1] "300 times: (mu=1.736, sigma2=9.631, sigma=3.103)"
## [1] "350 times: (mu=1.736, sigma2=9.611, sigma=3.100)"
## [1] "400 times: (mu=1.736, sigma2=9.608, sigma=3.100)"
## [1] "450 times: (mu=1.736, sigma2=9.607, sigma=3.099)"
## [1] "500 times: (mu=1.736, sigma2=9.607, sigma=3.099)"

対数尤度関数をそのまま使う方法もメモしておく。この式は少し分解すると、下記の通り分解でき、第1項はパラメタが含まれていないので、最大化する際はあってもなくてもよい。

\[ \begin{eqnarray} log L &=& - \frac{n}{2}log(2 \pi \sigma^2) - \frac{1}{2 \sigma^2} \sum (x_{i}-\mu)^2 \\ &=& - \frac{n}{2}log(2 \pi) - \frac{n}{2}log(\sigma^2)- \frac{1}{2 \sigma^2} \sum (x_{i}-\mu)^2 \\ \end{eqnarray} \]

そのため、第1項を除いた下記の式をここでは利用する。

\[ \begin{eqnarray} log L &=& - \frac{n}{2}log(\sigma^2)- \frac{1}{2 \sigma^2} \sum (x_{i}-\mu)^2 \\ \end{eqnarray} \]

さきほどと同様の結果が得られている。

mu <- 1
sigma <- 1

loglik <- function(mu, sig){
  - (n/2)*log(sig^2) - (1/(2*sig^2))*sum((x - mu)^2)
}

for(i in seq_along(1:iter)){
  mu <- mu + eta * (loglik(mu + h, sigma) - loglik(mu, sigma)) / h
  sigma <- sigma + eta * (loglik(mu, sigma + h) - loglik(mu, sigma)) / h

  if (i %% 50 == 0) print(sprintf("%d times: (mu=%2.3f, sigma2=%2.3f, sigma=%2.3f)", i, mu, sigma^2, sigma))
}

## [1] "50 times: (mu=1.731, sigma2=16.012, sigma=4.001)"
## [1] "100 times: (mu=1.731, sigma2=9.576, sigma=3.095)"
## [1] "150 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "200 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "250 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "300 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "350 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "400 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "450 times: (mu=1.731, sigma2=9.576, sigma=3.094)"
## [1] "500 times: (mu=1.731, sigma2=9.576, sigma=3.094)"

optim()でも同様の計算結果が得られている。

loglik2 <- function(x){
  return(function(par){
    mu <- par[1]
    sig <- par[2]
    - n / 2*log(sig^2) - 1/2 * (sum((x - mu)^2) / sig^2)
    # 分散sigma2を直接計算したい場合
    # sigma2 <- par[2]
    # - n / 2*log(sigma2) - 1/2 * (sum((x - mu)^2) / sigma2)
})
}

res <- optim(
  par = c(1, 10), 
  fn = loglik2(x),
  # 最大化するための引数
  control = list(fnscale = -1)
)

list(
  mu = res$par[1],
  sigma2 = (res$par[2])^2,
  sigma = res$par[2]
)

## $mu
## [1] 1.735846
## 
## $sigma2
## [1] 9.605287
## 
## $sigma
## [1] 3.09924

分散\(\sigma^{2}\)を直接計算したい場合は下記の通り書き直せば良い。

loglik3 <- function(x){
  return(function(par){
    mu <- par[1]
    # 分散sigma2を直接計算したい場合
    sigma2 <- par[2]
    - n / 2*log(sigma2) - 1/2 * (sum((x - mu)^2) / sigma2)
})
}

res <- optim(
  par = c(1, 10), 
  fn = loglik3(x),
  # 最大化するための引数
  control = list(fnscale = -1)
)

list(
  mu = res$par[1],
  sigma2 = res$par[2],
  sigma = sqrt(res$par[2])
)

## $mu
## [1] 1.736
## 
## $sigma2
## [1] 9.609064
## 
## $sigma
## [1] 3.099849