(R語言)不依賴求導的尋根方法：Nelder-Mead方法

牛頓法需要一階和二階導數，如果導函數不可用，則可以通過僞牛頓(Quasi-Newton)方法近似它們。如果導函數直接不存在怎麼辦？函數中存在不連續則可能發生這種情況。

Nelder-Mead Method

在Nelder-Mead算法下即使函數不連續，Nelder-Mead算法也是十分穩健的，該算法是基于評估n維單純形的頂點處的函數，其中n是函數中輸入變量的數量(在這裡我隻寫計算二維問題的算法)。對于二維問題，n維單純形隻是一個三角形，每個角的一個頂點，通常有n+1個頂點。

下面是算法的具體步驟：

Step1: 下面為三個初始值點 X ( 1 ) = ( x 1 , x 2 ) X(1) = (x_1,x_2) X(1)=(x1​,x2​) X ( 2 ) = ( x 3 , x 4 ) X(2) = (x_3,x_4) X(2)=(x3​,x4​) X ( 3 ) = ( x 5 , x 6 ) X(3) = (x_5,x_6) X(3)=(x5​,x6​)

将上述三個初始值點代入函數 f ( x ) f(x) f(x)中，即 f ( X 1 ) f(X_1) f(X1​)、 f ( X 2 ) f(X_2) f(X2​)、 f ( X 3 ) f(X_3) f(X3​)，再對它們進行排序，假設三個初始值代入的大小排序為： f ( X 1 ) ≤ f ( X 2 ) ≤ f ( X 3 ) f(X_1)\leq f(X_2)\leq f(X_3) f(X1​)≤f(X2​)≤f(X3​)

Step2: 計算這三個點的中心點 X ( 0 ) = ∑ i = 1 n X i n X(0) = \frac {\sum_{i=1}^n X_i} {n} X(0)=n∑i=1n​Xi​​

Step3: 計算對稱點 X ( r ) = X ( 0 ) + α ( X ( 0 ) − X ( 3 ) ) X(r) = X(0) + \alpha(X(0)-X(3)) X(r)=X(0)+α(X(0)−X(3))

通常情況下 α \alpha α取1

Case1: f ( X 1 ) ≤ f ( X r ) ≤ f ( X 3 ) f(X_1)\leq f(X_r) \leq f(X_3) f(X1​)≤f(Xr​)≤f(X3​)此時 X r X_r Xr​不是最優點也不是最差點，用 X r X_r Xr​取代 X 3 X_3 X3​，新的三個構造點變為 X 1 X_1 X1​， X 2 X_2 X2​， X r X_r Xr​，回到Step1，繼續進行疊代運算。

Case2: f ( X r ) &lt; f ( X 1 ) f(X_r) &lt; f(X_1) f(Xr​)<f(X1​)

此時 X r X_r Xr​是最優點，一個好的方向被找出，沿着此方向再找出點 X e X_e Xe​:

X e = X 0 + γ ( X r − X 0 ) X_e = X_0 + \gamma (X_r - X_0) Xe​=X0​+γ(Xr​−X0​)

通常情況下 γ \gamma γ取2

計算 f ( X e ) f(X_e) f(Xe​)，當 f ( X e ) &lt; f ( X r ) f(X_e) &lt; f(X_r) f(Xe​)<f(Xr​)時，擴充點是優于對稱點的，此時用 X 1 X_1 X1​、 X 2 X_2 X2​、 X e X_e Xe​構造；當 f ( X r ) ≤ f ( X e ) f(X_r) \leq f(X_e) f(Xr​)≤f(Xe​)時，對稱點是優于擴充點的，此時用 X 1 X_1 X1​、 X 2 X_2 X2​、 X r X_r Xr​構造，然後回到Step1。

Case3: f ( X r ) ≥ f ( X 3 ) f(X_r) \ge f(X_3) f(Xr​)≥f(X3​)

此時 X r X_r Xr​是最差點，需要構造收縮點 X c X_c Xc​：

X c = X 0 + ρ ( X 3 − X 0 ) X_c = X_0 + \rho (X_3 - X_0) Xc​=X0​+ρ(X3​−X0​)

計算 f ( X c ) f(X_c) f(Xc​)，當 f ( X c ) &lt; f ( X 3 ) f(X_c) &lt; f(X_3) f(Xc​)<f(X3​)時，此時用 X 1 X_1 X1​、 X 2 X_2 X2​、 X c X_c Xc​構造；當 f ( X 3 ) ≤ f ( X c ) f(X_3) \leq f(X_c) f(X3​)≤f(Xc​)時，構造新的 X 2 ′ X_2&#x27; X2′​和 X 3 ′ X_3&#x27; X3′​：

X 2 ′ = X 1 + δ ( X 2 − X 1 ) X_2&#x27; = X_1 + \delta (X_2 - X_1) X2′​=X1​+δ(X2​−X1​) X 3 ′ = X 1 + δ ( X 3 − X 1 ) X_3&#x27; = X_1 + \delta (X_3 - X_1) X3′​=X1​+δ(X3​−X1​)

通常情況下 ρ \rho ρ和 δ \delta δ取0.5

此時用 X 1 X_1 X1​、 X 2 ′ X_2&#x27; X2′​、 X 3 ′ X_3&#x27; X3′​構造，然後回到Step1。

R代碼

此代碼還具備了生成GIF圖的功能，有興趣的朋友可以了解animation這個包，裡面有很多很有趣的函數。

setwd("/Users/Desktop/R programs")
png("NelderMead%03d.png")
library(animation)
NelderMead <- function(X, fun, alpha = 1, gamma = 2, lo = 0.5, 
                       mu = 0.5, eps = 10^(-8), max.loop = 50)
{
  i <- 1
  #build a empty vector to store iterative triangle area 
  area.hist <- vector(mode = "numeric", length = 0)
  area.hist[1] <- 1
  #eps is the absolute convergence criterion for triangle area
  while((i <= max.loop) & (area.hist[i] > eps))
  {
    i <- i + 1
    x1 <- X[1,1]
    x2 <- X[1,2]
    fx.1 <- eval(fun)
    x1 <- X[2,1]
    x2 <- X[2,2]
    fx.2 <- eval(fun)
    x1 <- X[3,1]
    x2 <- X[3,2]
    fx.3 <- eval(fun)
    FX <- c(fx.1, fx.2, fx.3)
    x.x <- X[,1]
    y.y <- X[,2]
    plot(X[,1], X[,2], xlim = c(-1,2), ylim = c(-1,2), xlab = "x", ylab = "y", 
         main = "Nelder-Mead Iteration", pch = 20, cex = 0.5)
    lines(rbind(c(X[1,1],X[1,2]), c(X[2,1],X[2,2]), c(X[3,1],X[3,2]), c(X[1,1],X[1,2])))
    #Get the index of the maximum and minimum、median values in FX
    maximum <- match(max(FX), FX)
    median <- match(median(FX), FX)
    minimum <- match(min(FX), FX)
    #X0: Centriod of the remaining n points
    X0 <- c(sum(X[1,1]+X[2,1]+X[3,1])/3, sum(X[1,2]+X[2,2]+X[3,2])/3)
    #Xr: Reflected point
    Xr <- c(X0[1]+alpha*(X0[1]-X[maximum,1]), X0[2]+alpha*(X0[2]-X[maximum,2]))
    x1 <- Xr[1]
    x2 <- Xr[2]
    fx.xr <- eval(fun)
    
    if((FX[minimum] <= fx.xr) && (fx.xr < FX[maximum]))
    {
      X[maximum,1] <- Xr[1]
      X[maximum,2] <- Xr[2]
      p <- matrix(c(X[1,1],X[1,2],1,X[2,1],X[2,2],1,X[3,1],X[3,2],1), nrow = 3, byrow = TRUE)
      area.hist[i] <- abs(0.5*det(p))
    }
    
    else if(fx.xr < FX[minimum])
    {
      Xe <- c(X0[1]+gamma*(Xr[1]-X0[1]), X0[2]+gamma*(Xr[2]-X0[2]))
      x1 <- Xe[1]
      x2 <- Xe[2]
      if(eval(fun) < fx.xr){
        X[maximum,1] <- Xe[1]
        X[maximum,2] <- Xe[2]
        p <- matrix(c(X[1,1],X[1,2],1,X[2,1],X[2,2],1,X[3,1],X[3,2],1), nrow = 3, byrow = TRUE)
        area.hist[i] <- abs(0.5*det(p))
      }
      else{
        X[maximum,1] <- Xr[1]
        X[maximum,2] <- Xr[2]
        p <- matrix(c(X[1,1],X[1,2],1,X[2,1],X[2,2],1,X[3,1],X[3,2],1), nrow = 3, byrow = TRUE)
        area.hist[i] <- abs(0.5*det(p))
      }
    }
    
    else if(fx.xr >= FX[maximum])
    {
      Xc <- c(X0[1]+lo*(X[maximum,1]-X0[1]), X0[2]+lo*(X[maximum,2]-X0[2]))
      x1 <- Xc[1]
      x2 <- Xc[2]
      if(eval(fun) < FX[maximum]){
        X[maximum,1] <- Xc[1]
        X[maximum,2] <- Xc[2]
        p <- matrix(c(X[1,1],X[1,2],1,X[2,1],X[2,2],1,X[3,1],X[3,2],1), nrow = 3, byrow = TRUE)
        area.hist[i] <- abs(0.5*det(p))
      }
      else{
        X[maximum,1] <- X[minimum,1]+mu*(X[maximum,1]-X[minimum,1])
        X[maximum,2] <- X[minimum,2]+mu*(X[maximum,2]-X[minimum,2])
        X[median,1] <- X[minimum,1]+mu*(X[median,1]-X[minimum,1])
        X[median,2] <- X[minimum,2]+mu*(X[median,2]-X[minimum,2])
        p <- matrix(c(X[1,1],X[1,2],1,X[2,1],X[2,2],1,X[3,1],X[3,2],1), nrow = 3, byrow = TRUE)
        area.hist <- abs(0.5*det(p))
      }
    }
    
  }
  
  return(list(Point_coordinates = X, Area = area.hist, Iteration = i))
}

#by the result, we can find the different initial value lead to different convergence result
X <- matrix(c(1,1,2,1,2,2), nrow = 3)
out <- NelderMead(X, expression(x1^2 + x2^2))
dev.off()
out
bm.files <- sprintf("NelderMead%03d.png", 1:25)
im.convert(files = bm.files, output = "NelderMead.gif")

           

—————————未經允許，謝絕轉載—————————