Как оптимизировать для целочисленных параметров (и другого разрывного пространства параметров) в R?

Ответ 1

Вот несколько идей.

1. Палированная оптимизация. Вы можете объединить аргументы целевой функции и добавьте штраф за нецелые числа. Но это создает множество локальных экстремумов, поэтому вы можете предпочесть более надежную процедуру оптимизации, например, дифференциальная эволюция или оптимизация рой частиц.

fr <- function(x) {
  x1 <- round( x[1] )
  x2 <- round( x[2] )
  value <- 100 * (x2 - x1 * x1)^2 + (1 - x1)^2
  penalty <- (x1 - x[1])^2 + (x2 - x[2])^2
  value + 1e3 * penalty
}

# Plot the function
x <- seq(-3,3,length=200)
z <- outer(x,x, Vectorize( function(u,v) fr(c(u,v)) ))
persp(x,x,z,
  theta = 30, phi = 30, expand = 0.5, col = "lightblue", border=NA,
  ltheta = 120, shade = 0.75, ticktype = "detailed")

library(RColorBrewer)
image(x,x,z, 
  las=1, useRaster=TRUE,
  col=brewer.pal(11,"RdYlBu"),
  xlab="x", ylab="y"
)

# Minimize
library(DEoptim)
library(NMOF)
library(pso)
DEoptim(fr, c(-3,-3), c(3,3))$optim$bestmem
psoptim(c(-2,1), fr, lower=c(-3,-3), upper=c(3,3))
DEopt(fr, list(min=c(-3,-3), max=c(3,3)))$xbest
PSopt(fr, list(min=c(-3,-3), max=c(3,3)))$xbest

2. Исчерпывающий поиск. Если пространство поиска мало, вы также можете использовать поиск по сетке.

library(NMOF)
gridSearch(fr, list(seq(-3,3), seq(-3,3)))$minlevels

3. Локальный поиск с указанными пользователем районами. Не настраивая целевую функцию, вы можете использовать некую форму локального поиска, в котором вы можете указать, какие точки следует проверять. Это должно быть намного быстрее, но чрезвычайно чувствительно к выбору функции окрестности.

# Unmodified function
f <- function(x) 
  100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2

# Neighbour function
# Beware: in this example, with a smaller neighbourhood, it does not converge.
neighbour <- function(x,...)
  x + sample(seq(-3,3), length(x), replace=TRUE)

# Local search (will get stuck in local extrema)
library(NMOF)
LSopt(f, list(x0=c(-2,1), neighbour=neighbour))$xbest
# Threshold Accepting
TAopt(f, list(x0=c(-2,1), neighbour=neighbour))$xbest

4. Поиск в Tabu. Чтобы избежать повторного использования одних и тех же пунктов, вы можете использовать tabu search, то есть помните последние k пунктов, чтобы не посещать их снова.

get_neighbour_function <- function(memory_size = 100, df=4, scale=1){
  # Static variables
  already_visited <- NULL
  i <- 1
  # Define the neighbourhood
  values <- seq(-10,10)
  probabilities <- dt(values/scale, df=df)
  probabilities <- probabilities / sum(probabilities)
  # The function itself
  function(x,...) {
    if( is.null(already_visited) ) {
      already_visited <<- matrix( x, nr=length(x), nc=memory_size )
    }
    # Do not reuse the function for problems of a different size
    stopifnot( nrow(already_visited) == length(x) )
    candidate <- x
    for(k in seq_len(memory_size)) {
      candidate <- x + sample( values, p=probabilities, length(x), replace=TRUE )
      if( ! any(apply(already_visited == candidate, 2, all)) )
        break
    }
    if( k == memory_size ) {
      cat("Are you sure the neighbourhood is large enough?\n")
    } 
    if( k > 1 ) {
      cat("Rejected", k - 1, "candidates\n")
    }
    if( k != memory_size ) {
      already_visited[,i] <<- candidate
      i <<- (i %% memory_size) + 1
    }
    candidate
  }
}

В следующем примере это не работает: мы переходим только к ближайшему локальному минимуму. И в более высоких измерениях все становится еще хуже: окрестности настолько велики, что мы никогда не попадаем в кеш уже посещенных точек.

f <- function(x) {
  result <- prod( 2 + ((x-10)/1000)^2 - cos( (x-10) / 2 ) )  
  cat(result, " (", paste(x,collapse=","), ")\n", sep="")
  result
}
plot( seq(0,1e3), Vectorize(f)( seq(0,1e3) ) )

LSopt(f, list(x0=c(0,0), neighbour=get_neighbour_function()))$xbest
TAopt(f, list(x0=c(0,0), neighbour=get_neighbour_function()))$xbest
optim(c(0,0), f, gr=get_neighbour_function(), method="SANN")$par

Дифференциальная эволюция работает лучше: мы получаем только локальный минимум, но он лучше, чем ближайший.

g <- function(x) 
  f(x) + 1000 * sum( (x-round(x))^2 )
DEoptim(g, c(0,0), c(1000,1000))$optim$bestmem

Поиск табу часто используется для чисто комбинаторных задач (например, когда пространство поиска представляет собой набор деревьев или графов) и не кажется отличной идеей для целочисленных проблем.

Ответ 2

Целочисленное программирование (IP) имеет свои собственные правила и алгоритмы. Использование непрерывного решателя не имеет большого смысла. R не имеет специализированного решателя целых программ, но вы можете попробовать:

• Если ваша функция линейна, используйте одно из смешанного целочисленного программирования такие как lp_solve как "lpSolve" в R или GLPK как "Rglpk" в R.

• В противном случае вы можете попробовать оптимизировать метод "SANN", имитируемый отжиг подход, о котором говорится в документации:

"It uses only function values but is relatively slow... If a function to generate a new candidate point is given, method 'SANN' can also be used to solve combinatorial optimization problems... Note that the 'SANN' method depends critically on the settings of the control parameters."

Вот пример с переведенной сферической функцией в [-10,10]x[-10,10]:

fun <- function(x) sum((x-c(3.2, 6.7))^2)
nextfun <- function(x) sample(-10:10, 2, replace=TRUE)

optim(fn=fun, par=c(-10,-10), gr=nextfun, method="SANN", 
      control=list(maxit=1000,fnscale=1,trace=10))

# sann objective function values
# initial       value 458.000000
# iter      999 value 0.000000
# final         value 0.000000
# sann stopped after 999 iterations
# $par
# [1] 3 7
# $value
# [1] 0.13

Но вы должны применять более интеллектуальный "градиент", который использует случайную выборку, или полный поиск через ваш домен целых чисел, если ничего не помогает. Конечно, в более высоких измерениях потребуется специализированный подход.

Ответ 3

В R доступны новые пакеты, которые позволяют использовать прерывистые входные параметры (например, integer) в программах оптимизации. Один из них rgenoud

Используя опцию "data.type.int = TRUE" и задав правильные границы, функция будет использовать только целые числа, чтобы минимизировать или максимизировать заданную функцию.

Под rgenoud использует статистику:: optim() Для оптимизации. Поэтому пользователь может передавать любые параметры rgenoud, которые обычно передаются в optim()