CMF, Machine Learning 3, Nelder-Mead simplex algorithm

nelder-mead.R

# Encoding UTF-8

# ЦМФ МГУ, математические финансы, весна 2015
# Machine Learning 3, Nelder-Mead simplex algorithm
# Орлов Евгений, 19.03.2015

f <- function(x) {
  x[1]^2 + x[2]^2
}

# функция для определения начального симплекса
init.simplex <- function(x0) {
  # матрица nrow == length(x0) + 1
  # по строкам - вершины симплекса
  
  # константы
  tau_nonzero <- 5 * 10^-2
  tau_zero <- 2.5 * 10^-4
  len <- length(x0)
  # инициализируем матрицу
  init_mat <- matrix(nrow = len + 1, ncol = len)
  # заполняем строки матрицы вершинами симплекса
  init_mat[1, ] <- x0
  for (i in 1:len) {
    tau <- ifelse(x0[i], tau_nonzero, tau_zero)
    e_vec <- rep(0, len)
    e_vec[i] <- 1
    init_mat[i + 1, ] <- x0 + tau * e_vec
  }
  # возвращаем матрицу
  init_mat
}

# сортировка вершин по значению целевой функции
sort.vertex <- function(simplex, fn) {
  # список из двух элементов: vertex и vertex.value
  # vertex - матрица вершин симплекса, отсортированных по
  # возрастанию значений функции на них
  # vertex.value - вектор значений функций на вершинах, отсортированный
  # по возрастанию
  
  # Вычисляем значения функции на вершинах симплекса
  fn_values <- rep(NA, nrow(simplex))
  for (i in 1:nrow(simplex)) {
      fn_values[i] <- fn(simplex[i, ]) 
  }
  # Возвращаем отсортированные данные
  rows_order <- order(fn_values)
  list(vertex = simplex[rows_order, ], vertex.value = sort(fn_values))
}

# начальная точка и начальный симплекс
x0 <- c(2.5, 2.5)
z <- sort.vertex(init.simplex(x0), f)
s <- z$vertex
vert.value <- z$vertex.value

# параметры оптимизации
d <- length(x0)
nm.par <- c(1, (d + 2) / d, (3 * d - 2) / (4 * d), (d - 1) / d)
names(nm.par) <- c("alpha","beta","gamma","delta")
val <- mean(vert.value)
delta.val <- Inf
delta <- 10^-12
values <- val

while (delta.val >= delta) {
  tmp <- mean(vert.value)  # сохраняем предыдущее значение
  # центроид наилучших вершин
  x.bar <- apply(s[1:d, ], 2, mean)
  # отражённая точка
  x.r <- x.bar + nm.par["alpha"] * (x.bar - s[d + 1, ])
  f.r <- f(x.r)
  shrinked <- FALSE  # индикатор, было ли сокращение симплекса
  # проверка условий 2. и 3.
  if ((f.r >= vert.value[1]) && (f.r < vert.value[d])) {
    s[d + 1, ] <- x.r
    f.new.vert <- f.r
  } else { 
    if (f.r < vert.value[1]) {
      # продвинутая точка
      x.e <- x.bar + nm.par["beta"] * (x.r - x.bar)
      f.e <- f(x.e)
      if (f.e < vert.value[1]) {
        s[d + 1, ] <- x.e
        f.new.vert <- f.e
      } else {
        s[d + 1, ] <- x.r
        f.new.vert <- f.r
      }
    }
  }
  # проверка условия 4.
  if ((f.r >= vert.value[d]) && (f.r < vert.value[d + 1])) {
    # внешнее сжатие
    x.oc <- x.bar + nm.par["gamma"] * (x.r - x.bar)
    f.oc <- f(x.oc)
    if (f.oc <= f.r) {
      s[d + 1, ] <- x.oc
      f.new.vert <- f.oc
    } else {
      # сокращение симплекса
      for (i in 2:(d + 1))
        s[i, ] <- s[1, ] + nm.par["delta"]*(s[i, ] - s[1, ])
      shrinked <- TRUE
    }
  }
  # проверка условия 5.
  if (f.r >= vert.value[d + 1]) {
    # внутреннее сжатие
    x.ic <- x.bar - nm.par["gamma"] * (x.r - x.bar)
    f.ic <- f(x.ic)
    if (f.ic < vert.value[d + 1]) {
      s[d + 1, ] <- x.ic
      f.new.vert <- f.ic
    } else {
      # сокращение симплекса
      for (i in 2:(d+1))
        s[i, ] <- s[1, ] + nm.par["delta"] * (s[i, ] - s[1, ])
      shrinked <- TRUE
    }
  }
  # переоценка значений функции на вершинах, пересортировка симплекса
  if (shrinked) {
    z <- sort.vertex(s, f)
    s <- z$vertex
    vert.value <- z$vertex.value
  } else {
    vert.value[d + 1] <- f.new.vert
    z <- order(vert.value)
    s <- s[z, ]
    vert.value <- vert.value[z]
  }
  delta.val <- abs(tmp - mean(vert.value))
  values <- c(values, mean(vert.value))
}

opt <- s[1, ]
cat("Точка минимума: ", opt, "\n")
opt.value <- vert.value[1]
cat("Минимум:", opt.value, "\n")
plot(values, main = "Изменение среднего значения на вершинах", type = 'l')