R語言最優化(一維)

　　最優化問題是普遍存在的，以前上運籌學課的時候也接觸過最優化相關的問題，當時主要是理論課，並且關注的重點是單純形法、運輸問題以及圖論等，這里指的最優化是指函數的最優化，即函數的極值，由於尋找一個局部最優比尋找全局最優要簡單得多，所以這里的最優解也是指的局部最優解。

• 牛頓最優化方法

　　　　僅給出代碼，公式什么的。。。我不知道博客園怎么插入公式，，，，

newton <- function(f3, x0, tol=1e-9, n.max = 100){
  x <- x0
  f3.x <- f3(x)
  n <- 0
  while((abs(f3.x[2]) > tol) & (n < n.max)){
    x <- x - f3.x[2] / f3.x[3]
    f3.x <- f3(x)
    n <- n + 1
  }
  if(n == n.max){
    cat("newton failed to converge\n")
  }
  else{
    return(x)
  }
}

gamma.2.3 <- function(x){
  if(x < 0) return(c(0,0,0))
  if(x == 0) return(c(0,0,NaN))
  y <- exp(-2*x)
  return(c(4*x^2*y, 8*x*(1-x)*y, 8*(1-2*x^2)*y))
}

x0 <- seq(0,10,0.01)
x <- c()
for(i in 1:length(x0)){
  x[i] <- newton(gamma.2.3,x0[i])
}

　　代碼在各種初始值下取到的極值點，在某些點可能會出現明顯錯誤的極值點，因此使用的時候應該謹慎，多試試幾組值。

　　

• 黃金分割法

　　　　黃金分割法只能在一維情況下使用，但是不必知道函數的導數。其核心思想與夾逼求根的方法是一樣的，即不斷縮小區間范圍。

　　　　為了加速運算，取了含1+ρ的項。

gsection <- function(ftn, x.l, x.r, x.m, tol=1e-9){
  ###黃金分割率
  gr1 <- 1 + (1+sqrt(5))/2
  f.l <- ftn(x.l)
  f.r <- ftn(x.r)
  f.m <- ftn(x.m)
  while((x.r - x.l) > tol){
    if((x.r - x.m) > (x.m - x.l)){
      y <- x.m + (x.r - x.m)/grl
      f.y <- ftn(y)
      if(f.y >= f.m){
        x.l <- x.m
        f.l <- f.m
        x.m <- y
        f.m <- f.y
      }
      else{
        x.r <- y
        f.r <- f.y
      }
    }else{
      y <- x.m - (x.m - x.l)/grl
      f.y <- ftn(y)
      if(f.y >= f.m){
        x.r <- x.m
        f.r <- f.m
        x.m <- y
        f.m <- f.y
      }
      else{
        x.l <- y
        f.l <- f.y
      }
    }
  }
  return(x.m)
}
f.4 <- function(x){return((x-1.7)^2+1)}
gsection(f.4, 1.5, 2, 1.8)