(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")

           

—————————未经允许，谢绝转载—————————