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

本当は去年みたいな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
}