読者です 読者をやめる 読者になる 読者になる

共役勾配法を実装してみた

最適化理論 R

本当は去年みたいなplotをしたいんだけど、面倒(ry。

工学基礎 最適化とその応用 (新・工科系の数学)の4.6章で遊んでいる。二次関数の簡単なのをとりあえず。

f <- function(x) {
  x1 <- x[1]; x2 <- x[2]
  3 / 2 * x1^2 + x1 * x2 + x2^2 - 6 * x1 - 7 * x2
}

nabla_f <- function(x) {
  x1 <- x[1]; x2 <- x[2]
  c(3 * x1 + x2 -6, x1 + 2 * x2 -7)
}

alpha_k <- function(x, d) {
  A <- matrix(c(3,1,1,2), 2, 2) # 正定値対称行列
  - sum(d * nabla_f(x)) / (t(d) %*% A %*% d)
}

x0 <- c(2, 1)
d <- - nabla_f(x0)
x <- x0
k <- 1

while(TRUE){
  old_nabla_f <- sum(nabla_f(x)^2)
  old_x <- x
  x <- x + alpha_k(x, d) * d
  if(sum((x - old_x)^2) < 0.001) break 
  beta <- sum(nabla_f(x)^2) / old_nabla_f
  d <- - nabla_f(x) + beta * d
  k <- k + 1
}

ちゃんと最適解にいっていて、反復の回数も正しいことを確認。

> x
[1] 1 3
> k
[1] 3

なんとなくoptimの共役勾配法についても試してみる。もちろん同じ。

> optim(x0, fn=f, method="CG")
$par
[1] 1.000001 3.000000

$value
[1] -13.5

$counts
function gradient 
      59       25 

$convergence
[1] 0

$message
NULL

非線形共役勾配法(直線探索付き)

今回の問題では二次関数と分かっているからステップ幅が陽に決まるし簡単なんだけど、一般の場合はそうもいかない。そこで

  • アルミホの条件でステップ幅を決定
  • Bletcher-Reeves公式を使ってパラメータ\beta_kを決定

するようにした。基本的にコードに変更はなくて、こんな感じ。

armijo <- function(x, d, xi, tau) {
  beta <- 1
  while(f(x + beta * d) > f(x) + xi * beta * sum(nabla_f(x) * d)) {
    beta <- tau * beta
  }
  return(beta)
}

x0 <- c(2, 1)
d <- - nabla_f(x0)
x <- x0
k <- 1

xi <- 0.9
tau <- 0.1

while(TRUE){
  old_nabla_f <- sum(nabla_f(x)^2)
  old_x <- x
  x <- x + armijo(x, d, xi, tau) * d
  if(sum((x - old_x)^2) < 0.0001) break 
  beta <- sum(nabla_f(x)^2) / old_nabla_f
  d <- - nabla_f(x) + beta * d
  k <- k + 1
}

最適解はもちろん一緒だけど、反復回数に大分違いがあることが分かる。

> x
[1] 1.035180 2.939214
> k
[1] 32

工学基礎 最適化とその応用 (新・工科系の数学)

工学基礎 最適化とその応用 (新・工科系の数学)