R/lnsrch.R

In fda: Functional Data Analysis

lnsrch <- function (xold, fold, g, p, func, dataList, stpmax, 
                    itermax=20, TOLX=1e-10, dbglev=0) {
  # Arguments:
  # xold     ... The current parameter vector value.
  # fold     ... The current function value.
  # g        ... The current gradient vector.
  # p        ... The current search direction vector.
  # func     ... The name of the function being optimized.
  # dataList ... A list object containing objects specifying the 
  #              function to be minimized.
  # stpmax   ... The maximum step size.
  # itermax  ... The maximum number of iterations. Default 20.
  # TOLX     ... Convergence criterion.
  # dbglev   ... The level of output.  0: no output; 1: function and slope
  #   at each iteration.  2 and above: also results within iteration.
  
  #  Last modified 28 March 2024
  
  n     <- length(xold)
  check <- FALSE
  f2    <- 0
  alam2 <- 0
  ALF   <- 1e-4
  psum  <- sqrt(sum(p^2))
  #  scale if attempted step is too big
  if (psum > stpmax) {
    p <- p*(stpmax/psum)
  }
  slope <- sum(g*p)
  if (dbglev > 1) {
    cat("\n")
    cat("      ")
    cat(0)
    cat("      ")
    cat(round(slope,5))
    cat("      ")
    cat(round(fold,5))
  } 
  if (slope >= 0) {
    stop('Initial slope not negative.')
  }
  # compute lambdamin
  test <- 0
  for (i in 1:n) {
    temp <- abs(p[i])/max(abs(xold[i]),1)
    if (temp > test) {
      test <- temp
    }
  }
  alamin <- TOLX/test
  #  always try full Newton step first
  alam   <- 1
  #  start of iteration loop
  iter <- 0
  while (iter <= itermax) {
    iter <- iter + 1
    x <- xold + alam*p
    #  -------------  function evaluation  -----------
    result <- func(x, dataList)
    f <- result$PENSSE
    g <- result$DPENSSE
    if (dbglev > 1) {
      cat("\n")
      cat("      ")
      cat(iter)
      cat("      ")
      cat(round(slope,5))
      cat("      ")
      cat(round(fold,5))
    }
    #  -----------------------------------------------
    #  convergence on x change.
    if (alam < alamin) {
      x     <- xold
      check <- TRUE
      print("alam < alamin ")
      return(list(x=x, check=check))
    } else {
      #  sufficient function decrease
      if (f <= fold + ALF*alam*slope) {
        return(list(x=x, check=check))
      }
      #  backtrack
      if (alam == 1) {
        #  first time
        tmplam <- -slope/(2*(f-fold-slope))
      } else {
        #  subsequent backtracks
        rhs1 <- f  - fold - alam *slope
        rhs2 <- f2 - fold - alam2*slope
        a <- (rhs1/alam^2 - rhs2/alam2^2)/(alam-alam2)
        b <- (-alam2*rhs1/alam^2 + alam*rhs2/(alam*alam2))/(alam-alam2)
        if (a == 0) {
          tmplam <- -slope/(2*b)
        } else {
          disc <- b^2 - 3*a*slope
          if (disc < 0) {
            tmplam <- 0.5*alam
          } else {
            if (b <= 0) {
              tmplam <- (-b+sqrt(disc))/(3*a)
            } else {
              tmplam <- -slope/(b+sqrt(disc))
            }
          }
          if (tmplam > 0.5*alam) {
            tmplam <- 0.5*alam
          }
        }
      }
      alam2 <- alam
      f2    <- f
      #  lambda > 0.1 lambda1
      alam <- max(tmplam, 0.1*alam)
    }
    #  try again
  }
  
  return(list(x=x, check=check))
  
}