R语言)不依赖求导的寻根方法：Nelder-Mead方法

(R语言)不依赖求导的寻根方法：Nelder-Mead方法

牛顿法需要一阶和二阶导数，如果导函数不可用，则可以通过伪牛顿(Quasi-Newton)方法近似它们。如果导函数直接不存在怎么办？函数中存在不连续则可能发生这种情况。

Nelder-Mead Method

在Nelder-Mead算法下即使函数不连续，Nelder-Mead算法也是十分稳健的，该算法是基于评估n维单纯形的顶点处的函数，其中n是函数中输入变量的数量(在这里我只写计算二维问题的算法)。对于二维问题，n维单纯形只是一个三角形，每个角的一个顶点，通常有n+1个顶点。

下面是算法的具体步骤：
Step1: 下面为三个初始值点X(1)=(x1,x2)X(1) = (x_1,x_2)X(1)=(x1,x2) X(2)=(x3,x4)X(2) = (x_3,x_4)X(2)=(x3,x4) X(3)=(x5,x6)X(3) = (x_5,x_6)X(3)=(x5,x6)
将上述三个初始值点代入函数f(x)f(x)f(x)中，即f(X1)f(X_1)f(X1)、f(X2)f(X_2)f(X2)、f(X3)f(X_3)f(X3)，再对它们进行排序，假设三个初始值代入的大小排序为：f(X1)≤f(X2)≤f(X3)f(X_1)\leq f(X_2)\leq f(X_3)f(X1)≤f(X2)≤f(X3)
Step2: 计算这三个点的中心点 X(0)=∑i=1nXinX(0) = \frac {\sum_{i=1}^n X_i} {n}X(0)=n∑i=1nXi
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(X1)≤f(Xr)≤f(X3)f(X_1)\leq f(X_r) \leq f(X_3)f(X1)≤f(Xr)≤f(X3)此时XrX_rXr不是最优点也不是最差点，用XrX_rXr取代X3X_3X3，新的三个构造点变为X1X_1X1，X2X_2X2，XrX_rXr，回到Step1，继续进行迭代运算。

Case2: f(Xr)<f(X1)f(X_r) < f(X_1)f(Xr)<f(X1)
此时XrX_rXr是最优点，一个好的方向被找出，沿着此方向再找出点XeX_eXe:
Xe=X0+γ(Xr−X0)X_e = X_0 + \gamma (X_r - X_0)Xe=X0+γ(Xr−X0)
通常情况下γ\gammaγ取2
计算f(Xe)f(X_e)f(Xe)，当f(Xe)<f(Xr)f(X_e) < f(X_r)f(Xe)<f(Xr)时，扩展点是优于对称点的，此时用X1X_1X1、X2X_2X2、XeX_eXe构造；当f(Xr)≤f(Xe)f(X_r) \leq f(X_e)f(Xr)≤f(Xe)时，对称点是优于扩展点的，此时用X1X_1X1、X2X_2X2、XrX_rXr构造，然后回到Step1。

Case3: f(Xr)≥f(X3)f(X_r) \ge f(X_3)f(Xr)≥f(X3)
此时XrX_rXr是最差点，需要构造收缩点XcX_cXc：
Xc=X0+ρ(X3−X0)X_c = X_0 + \rho (X_3 - X_0)Xc=X0+ρ(X3−X0)
计算f(Xc)f(X_c)f(Xc)，当f(Xc)<f(X3)f(X_c) < f(X_3)f(Xc)<f(X3)时，此时用X1X_1X1、X2X_2X2、XcX_cXc构造；当f(X3)≤f(Xc)f(X_3) \leq f(X_c)f(X3)≤f(Xc)时，构造新的X2′X_2'X2′和X3′X_3'X3′：
X2′=X1+δ(X2−X1)X_2' = X_1 + \delta (X_2 - X_1)X2′=X1+δ(X2−X1) X3′=X1+δ(X3−X1)X_3' = X_1 + \delta (X_3 - X_1)X3′=X1+δ(X3−X1)
通常情况下ρ\rhoρ和δ\deltaδ取0.5
此时用X1X_1X1、X2′X_2'X2′、X3′X_3'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 areawhile((i <= max.loop) & (area.hist[i] > eps)){i <- i + 1x1 <- 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 FXmaximum <- match(max(FX), FX)median <- match(median(FX), FX)minimum <- match(min(FX), FX)#X0: Centriod of the remaining n pointsX0 <- 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 pointXr <- 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")

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