Optimización univariada: Método de fase de acotamiento en R

El método de fase de acotamiento es usado para "encerrar" el mínimo de una función. Este método garantiza encerrar el mínimo de una función unimodal. El algoritmo comienza con una apuesta inicial y partir de ahí encuentra una dirección de búsqueda al efectuar dos evaluaciones más en la vecindad del punto inicial. Posteriormente, una estrategia de búsqueda exponencial es adoptada to alcanzar el óptimo. Asumiendo que buscamos el mínimo, a continuación se presenta el algoritmo:

Algoritmo 

Paso 1: Elige un punto inicial x⁽⁰⁾ y un incremento ∆. Asigna k=0.
Paso 2: Si f(x⁽⁰⁾ − |∆|) ≥ f(x⁽⁰⁾ ) ≥ f(x⁽⁰⁾ + |∆|), luego ∆ es positivo;
si no, si f(x⁽⁰⁾ − |∆|) ≤ f(x⁽⁰⁾ ) ≤ f(x⁽⁰⁾ + |∆|), luego ∆ es negativo;

si no, ir al paso 1.

Paso 3: Asigna x^(k+1) = x^(k) + 2^k ∆.
Paso 4: Si f(x^(k+1) ) < f(x^(k)), asigna k = k + 1 y ve al paso 3;
si no, el mínimo está en el intervalo (x^(k−1) , x^(k+1)) y termina.

Si el valor de ∆ es grande, la precisión del intervalo al punto mínimo es baja pero es encontrado rápidamente. Y por el contrario, si el valor de ∆ es pequeño, la precisión del intervalo es mejor, pero más evaluaciones de la función pueden ser necesarias para encerrar al mínimo. Se presenta un ejemplo para encontrar el intervalo que encierra al mínimo de la función:

f(x) = x² + 54/x

 Ahora presentamos la implementación del algoritmo en R.
fase_acotamiento <- function(x0 = 0.6, delta = 0.5, evalua = function(x){x^2 + 54/x}, trace = FALSE){
  solucion <- NULL
  x <- x0
 
  #  Paso 1:
  k = 0

  # Paso 2:
  fx_minus_delta = evalua(x0 - abs(delta))
  fx0 = evalua(x0)
  fx_plus_delta = evalua(x0 + abs(delta))
  if (fx_minus_delta  >= fx0  & fx0 >= fx_plus_delta){
    delta = abs(delta)
  }else if(fx_minus_delta  <= fx0  & fx0 <= fx_plus_delta){
    delta = abs(delta) * -1
  } else{
    stop("Elige un nuevo punto inicial")
  }

  while(is.null(solucion)){
    # Paso 3:
    assign( paste0("x", k+1), get(paste0("x",k)) + 2^k*delta)

    # Paso 4:
    if (evalua( get(paste0("x",k+1)) ) < evalua( get(paste0("x",k))) ){
      k = k + 1
      if (trace) cat("Iteracion:", k, "\tx: ", get(paste0("x",k)) , "\tf(x)", evalua(get(paste0("x",k))) ,"\n")
    }else{
      solucion <- c( get(paste0("x", k-1)),   get(paste0("x", k+1))  )
    }
  }
  if (trace)  cat("\nEl mínimo se encuentra en el intervalo:\n")
  sort(solucion)
}

En el siguiente ejemplo, aplicamos nuestro algoritmo sobre la función presentada, con un valor de incremento ∆ relativamente pequeño, así como un punto de inicio x⁽⁰⁾ en 0.6.

> # Ejemplo
> fase_acotamiento(x0=0.6, delta = 0.000001, trace = TRUE)
Iteracion: 1     x:  0.600001     f(x) 90.35985
Iteracion: 2     x:  0.600003     f(x) 90.35955
Iteracion: 3     x:  0.600007     f(x) 90.35896
Iteracion: 4     x:  0.600015     f(x) 90.35777
Iteracion: 5     x:  0.600031     f(x) 90.35539
Iteracion: 6     x:  0.600063     f(x) 90.35063
Iteracion: 7     x:  0.600127     f(x) 90.34111
Iteracion: 8     x:  0.600255     f(x) 90.32207
Iteracion: 9     x:  0.600511     f(x) 90.28403
Iteracion: 10     x:  0.601023     f(x) 90.20804
Iteracion: 11     x:  0.602047     f(x) 90.05645
Iteracion: 12     x:  0.604095     f(x) 89.75484
Iteracion: 13     x:  0.608191     f(x) 89.15779
Iteracion: 14     x:  0.616383     f(x) 87.9878
Iteracion: 15     x:  0.632767     f(x) 85.73986
Iteracion: 16     x:  0.665535     f(x) 81.58067
Iteracion: 17     x:  0.731071     f(x) 74.3987
Iteracion: 18     x:  0.862143     f(x) 63.37791
Iteracion: 19     x:  1.124287     f(x) 49.29446
Iteracion: 20     x:  1.648575     f(x) 35.47336
Iteracion: 21     x:  2.697151     f(x) 27.29575

El mínimo se encuentra en el intervalo:
[1] 1.648575 4.794303