f <- function(x) { return((x[1]-1)^2 + 100*(x[1]^2-x[2])^2) } # Rosenbrock関数

g <- function(x) { # gradient df/dx
  dy <- numeric(2)
  dy[1] <- 400*x[1]^3 + (2-400*x[2])*x[1] - 2
  dy[2] <- -200*x[1]^2 + 200*x[2]
  return(dy)
}

H <- function(x) { # Hesse matrix
  ddy <- matrix(numeric(2*2), ncol=2)
  ddy[1,1] <- 1200*x[1]^2 -400*x[2] + 2
  ddy[1,2] <- -400*x[1]
  ddy[2,1] <- -400*x[1]
  ddy[2,2] <- 200
  return(ddy)
}

maxg <- function(b) { return((t(b)%*%b)[1,1]) }


Newton <- function(x0, f, g, H, n) {
# x0: 解xの初期値，f: 目的関数，g: fの勾配ベクトル，H: fのヘッセ行列，n: 繰り返しの上限
  x <- x0 # 初期値
  for (i in 1:n) {
    y <- f(x)
    dy <- g(x)
    if (maxg(dy) < 1.0e-3) { # g(x)が零ベクトルに近ければ収束したとみなす
      break;
    }
    x1 <- x - solve(H(x)) %*% g(x) # Hの逆行列 × g
    x <- x1
  }
  return( list(y,x,i) ) # y: fの最小値，x: 最適解，i: 収束までの繰り返し数
}


# example

x0 <- c(10, 10)
res <- Newton(x0, f, g, H, 100)