一维搜索技术

本文的插图来自于http://doc.mbalib.com/view/15124720fa0fe598b821b04388b05d37.html

首先，我们要确定搜索区间[a,b]。可采用进退法，算法如下：

有了a,b之后，我们可以采用黄金分割法来寻找最佳步长。

以下为R代码 in 云盘/wangwei/R/richard.R

      df <- function(x)# evaluate the gradient of a function
{


    h= 0.001# step length
    df <- c()#collect the gradient
    n <- length(x) # the number of variables in the x
    for(i in 1:n)
        {
            h1 <- rep(0,n)
            h1[i] <- h   # the step
            xp <- x+h1
            xn <- x-h1
            h2 <- rep(0,n)
            h2[i] <- 2*h   # the double step
            xp2 <- x+h2
            xn2 <- x-h2
            dfnew <- (8*(f(xp)-f(xn))-f(xp2)+f(xn2))/(12*h)
            df <- c(df,dfnew)
        }
    return(df)
}


backandforword <- function(x,h0)# find the linear search interval
    {
        alphak <- rep(0,length(x))
        hk <-rep(h0,length(x))
        t <- 2
        k <- 0
        maxIteration <- 1000
        Iteration <- 0
        phik <- f(x+alphak)
        while(Iteration<maxIteration)
            {


                alphak1 <- alphak+hk
                phik1 <- f(x+alphak1)
                if(phik1<phik)
                    {
                        hk <- t*hk
                        alpha <- alphak
                        alphak <- alphak1
                        phik <- phik1
                        k <- k+1
                        next
                     }
                 else {
                     if(k==0)
                         {
                             hk <- -1*hk
                             alphak <- alphak1
                             phik <-phik1
                             next
                         }
                     else{
                         a <- min(alpha,alphak1)
                         b <- max(alpha,alphak1)
                         break
                        }
                     }
               Iteration <- Iteration+1
            }
        if (abs(a)<10000&abs(b)<10000)
            {
                return(c(a,b))
            }


    }


#interpolation
f <- function(x) #the expression of a function
    {
        value <- x^2+1+3*x
        return(value)
    }


goldmethod(0)
goldmethod <- function(x)
                 {
                     #interval <- backandforword(x,1)
                     interval <- c(-2,2)
                     e <- 0.0001
                     a <- interval[1]
                     b <- interval[2]
                     x1 <- b-0.618*(b-a)
                     x2 <- a+0.618*(b-a)
                     f1 <- f(x+x1*rep(1,length(x)))
                     f2 <- f(x+x2*rep(1,length(x)))
                     while(abs(x1-x2)>e)
                         {
                             if(f1<f2)
                                 {
                                     b <- x2
                                     x2 <- x1
                                     f2 <- f1
                                     x1 <- b-0.618*(b-a)
                                     f1 <- f(x+x1*rep(1,length(x)))
                                 }
                             else
                                 {
                                     a <- x1
                                     x1 <- x2
                                     f1 <- f2
                                     x2 <- a+0.618*(b-a)
                                     f2 <- f(x+x2*rep(1,length(x)))
                                 }


                         }
                     xstar <- 0.5*(x1+x2)
                     print(f(xstar))
                     return(xstar)


                 }