R/cmaes.R

In jakobbossek/cmaesr: Covariance Matrix Adaptation Evolution Strategy

#' @title Covariance-Matrix-Adaptation
#'
#' @description
#' Performs non-linear, non-convex optimization by means of the Covariance
#' Matrix Adaptation - Evolution Strategy (CMA-ES).
#'
#' @details
#' This is a pure R implementation of the popular CMA-ES optimizer for continuous
#' black box optimization [2, 3]. It features a flexible system of stopping conditions
#' and enables restarts [1], which can be triggered by arbitrary stopping conditions
#' and can lead to superior performance on multimodal problems.
#'
#' You may pass additional parameters to the CMA-ES via the \code{control} argument.
#' This argument must be a named list. The following control elements will be considered
#' by the CMA-ES implementation:
#' \describe{
#'   \item{lambda [\code{integer(1)}]}{Number of offspring generated in each generation.}
#'   \item{mu [\code{integer(1)}]}{Number of individuals in each population. Defaults to \eqn{\lfloor \lambda / 2\rfloor}.}
#'   \item{weights [\code{numeric}]}{Numeric vector of positive weights.}
#'   \item{sigma [\code{numeric(1)}]}{Initial step-size. Default is 0.5.}
#'   \item{restart.triggers [\code{character}]}{List of stopping condition codes / short names (see
#'   \code{\link{makeStoppingCondition}}). All stopping conditions which are placed in this vector do trigger a restart
#'   instead of leaving the main loop. Default is the empty character vector, i.e., restart is not triggered.}
#'   \item{max.restarts [\code{integer(1)}]}{Maximal number of restarts. Default is 0. If set
#'   to >= 1, the CMA-ES is restarted with a higher population size if at least one of the
#'   stoppping conditions is defined as a restart trigger \code{restart.triggers}.}
#'   \item{restart.multiplier [\code{numeric(1)}]}{Factor which is used to increase the population size after restart.
#'   Default is 2.}
#'   \item{stop.ons [\code{list}]}{List of stopping conditions. The default is to stop after 10 iterations or after a
#'   kind of a stagnation (see \code{\link{getDefaultStoppingConditions}}).}
#'   \item{log.population [\code{logical(1L)}]}{Should each population be stored? Default is \code{FALSE}.}
#' }
#'
#' @note Internally a check for an indefinite covariance matrix is always performed, i.e.,
#' this stopping condition is always prepended internally to the list of stopping conditions.
#'
#' @references
#' [1] Auger and Hansen (2005). A Restart CMA Evolution Strategy With Increasing
#' Population Size. In IEEE Congress on Evolutionary Computation, CEC 2005, Proceedings,
#' pp. 1769-1776.
#' [2] N. Hansen (2006). The CMA Evolution Strategy: A Comparing Review. In J.A. Lozano,
#' P. Larranaga, I. Inza and E. Bengoetxea (Eds.). Towards a new evolutionary computation.
#' Advances in estimation of distribution algorithms. Springer, pp. 75-102.
#' [3] Hansen and Ostermeier (1996). Adapting arbitrary normal mutation distributions in evolution
#' strategies: The covariance matrix adaptation. In Proceedings of the 1996 IEEE
#' International Conference on Evolutionary Computation, pp. 312-317.
#'
#' @keywords optimize
#'
#' @param objective.fun [\code{smoof_function}]\cr
#'   Continuous objective function of type \code{smoof_function}. The function
#'   must expect a vector of numerical values and return a scaler numerical value.
#' @param start.point [\code{numeric}]\cr
#'   Initial solution vector. If \code{NULL}, one is generated randomly within the
#'   box constraints offered by the paramter set of the objective function.
#'   Default is \code{NULL}.
#' @param monitor [\code{cma_monitor}]\cr
#'   Monitoring object.
#'   Default is \code{\link{makeSimpleMonitor}}, which produces a console output.
#' @param control [\code{list}]\cr
#'   Futher paramters for the CMA-ES. See the details section for more in-depth
#'   information. Stopping conditions are also defined here.
#'   By default only some stopping conditions are passed. See \code{\link{getDefaultStoppingConditions}}.
#' @return [\code{cma_result}] Result object. Internally a list with the following
#'   components:
#'   \describe{
#'     \item{par.set [\code{\link[ParamHelpers]{ParamSet}}]}{Parameter set of the objective function.}
#'     \item{best.param [\code{numeric}]}{Final best parameter setting.}
#'     \item{best.fitness [\code{numeric(1L)}]}{Fitness value of the \code{best.param}}.
#'     \item{n.evals [\code{integer(1L)}]}{Number of function evaluations performed.}
#'     \item{past.time [\code{integer(1L)}]}{Running time of the optimization in seconds.}
#'     \item{n.restarts [\code{integer(1L)}]}{Number of restarts.}
#'     \item{population.trace [\code{list}]}{Trace of population.}
#'     \item{message [\code{character(1L)}]}{Message generated by stopping condition.}
#'   }
#'
#' @examples
#' # generate objective function from smoof package
#' fn = makeRosenbrockFunction(dimensions = 2L)
#' res = cmaes(
#'   fn,
#'   monitor = NULL,
#'   control = list(
#'     sigma = 1.5,
#'     lambda = 40,
#'     stop.ons = c(list(stopOnMaxIters(100L)), getDefaultStoppingConditions())
#'   )
#' )
#' print(res)
#'
#' @export
cmaes = function(
  objective.fun,
  start.point = NULL,
	monitor = makeSimpleMonitor(),
  control = list(
    stop.ons = c(
      getDefaultStoppingConditions()
    )
  )) {
	assertClass(objective.fun, "smoof_function")

	# extract relevant data
	par.set = ParamHelpers::getParamSet(objective.fun)
  lb = getLower(par.set); ub = getUpper(par.set)
	n = getNumberOfParameters(objective.fun)

	# sanity checks
	if (isNoisy(objective.fun)) {
		stopf("Noisy optimization is not supported at the moment.")
	}

	if (!isNumeric(par.set, include.int = FALSE)) {
		stopf("CMA-ES only works for objective functions with numeric parameters.")
	}

	if (isMultiobjective(objective.fun)) {
		stopf("CMA-ES can only handle single-objective functions.")
	}

	if (!is.null(start.point)) {
		assertNumeric(start.point, len = n, any.missing = FALSE)
	} else {
		if (!hasFiniteBoxConstraints(par.set)) {
			stopf("No start point provided. Cannot generate one, because parameter set cannot sample with Inf bounds!")
		}
		start.point = unlist(sampleValue(par.set))
	}

	if (!is.null(monitor)) {
		assertClass(monitor, "cma_monitor")
	}

  # get stopping conditions
  stop.ons = getCMAESParameter(control, "stop.ons", NULL)
  if (is.null(stop.ons)) {
    stopf("There must be at least one stopping condition!")
  }
  assertList(stop.ons, min.len = 1L, types = "cma_stopping_condition")

  # alwas check for indefinit covariance matrix first
  stop.ons = c(list(stopOnIndefCovMat()), stop.ons)

  # restart mechanism (IPOP-CMA-ES)
  restart.triggers = getCMAESParameter(control, "restart.triggers", character(0L))
  stop.ons.names = sapply(stop.ons, function(stop.on) stop.on$code)
  if (!isSubset(restart.triggers, stop.ons.names)) {
    stopf("Only codes / short names of active stopping conditions allowed as restart trigger, but '%s' are no stopping conditions.", collapse(setdiff(restart.triggers, stop.ons.names), sep = ", "))
  }
  restart.multiplier = getCMAESParameter(control, "restart.multiplier", 2)
  assertNumber(restart.multiplier, lower = 1, na.ok = FALSE, finite = TRUE)
  max.restarts = getCMAESParameter(control, "max.restarts", 0L)
  assertInt(max.restarts)

  #FIXME: default value should be derived from bounds
  sigma = getCMAESParameter(control, "sigma", 0.5)
  assertNumber(sigma, lower = 0L, finite = TRUE)

  # Precompute E||N(0,I)||
	chi.n = sqrt(n) * (1 - 1 / (4 * n) + 1 / (21 * n^2))

	# bookkeep best individual
	best.param = rep(NA, n)
	best.fitness = Inf

  # set initial distribution mean
	m = start.point

  # logs
  log.population = getCMAESParameter(control, "log.population", FALSE)
  assertFlag(log.population, na.ok = FALSE)
  population.trace = list()

  # init some termination criteria stuff
	iter = 0L
  n.evals = 0L
	start.time = Sys.time()

	callMonitor(monitor, "before")

  # somehow dirty trick to "really quit" if stopping condition is met and
  # now more restart should be triggered.
  do.terminate = FALSE

  for (run in 0:max.restarts) {
    # population and offspring size
    if (run == 0) {
      lambda = getCMAESParameter(control, "lambda", 4L + floor(3 * log(n)))
      assertInt(lambda, lower = 4)
      mu = getCMAESParameter(control, "mu", floor(lambda / 2))
      assertInt(mu)
    } else {
      lambda = getCMAESParameter(control, "lambda", 4L + floor(3 * log(n)))
      # increase population size (IPOP-CMA-ES)
      lambda = ceiling(restart.multiplier^run * lambda)
      mu = floor(lambda / 2)
    }

    # path for covariance matrix C and stepsize sigma
    pc = rep(0, n)
    ps = rep(0, n)

    # initialize recombination weights
    weights = getCMAESParameter(control, "weights", log(mu + 0.5) - log(1:mu))
    if (any(weights < 0)) {
      stopf("All weights need to be positive, but there are %i negative ones.", sum(which(weights < 0)))
    }
    if (length(weights) != mu) {
      stopf("You need to pass %i 'weights', but you passed %i.", mu, length(weights))
    }

    # normalize weights
    weights = weights / sum(weights)

    # variance-effectiveness / variance effective selection mass of sum w_i x_i
    mu.eff = sum(weights)^2 / sum(weights^2) # chosen such that mu.eff ~ lambda/4

    # step-size control
    cs = (mu.eff + 2) / (n + mu.eff + 5)
    ds = 1 + 2 * max(0, sqrt((mu.eff - 1) / (n + 1)) - 1) + cs # damping factor

    # covariance matrix Adaptation parameters
    cc = (4 + mu.eff / n) / (n + 4 + 2 * mu.eff / n)
    c1 = 2 / ((n + 1.3)^2 + mu.eff)
    alpha.mu = 2L
    cmu = min(1 - c1, alpha.mu * (mu.eff - 2 + 1/mu.eff) / ((n + 2)^2 + mu.eff))

    # covariance matrix
    sigma = getCMAESParameter(control, "sigma", 0.5)
    B = diag(n)
    D = diag(n)
    BD = B %*% D
    C = BD %*% t(BD) # C = B D^2 B^T = B B^T, since D equals I_n
    Cinvsqrt = B %*% diag(1 / sqrt(diag(D))) %*% t(B)

    # no restart trigger fired until now
    restarting = FALSE

    # break inner loop if terminating stopping condition active or
    # restart triggered
  	while (!restarting) {
      iter = iter + 1L

      # create new population of search points
  		arz = matrix(rnorm(n * lambda), ncol = lambda) # ~ N(0, I)
      ary = BD %*% arz # ~ N(0, C)
      arx = m + sigma * ary # ~ N(m, sigma^2 C)

      # Here we apply a penalization of violated bounds
      arx.repaired = ifelse(arx < lb, lb, ifelse(arx > ub, ub, arx))

      # Prepare penalization based on distance to repaired points (see Eq. 51)
      penalty.alpha = 1L
      penalties = penalty.alpha * colSums((arx - arx.repaired)^2)
      penalties[is.infinite(penalties)] = .Machine$double.max / 2

      # compute fitness values of repaired points
      fitn.repaired = if (isVectorized(objective.fun)) {
        objective.fun(arx.repaired)
      } else {
        apply(arx.repaired, 2L, function(x) objective.fun(x))
      }

      # apply penalization (see Eq. 51)
      fitn = fitn.repaired + penalties

      # update evaluation
  		n.evals = n.evals + lambda

      # order fitness values
      fitn.ordered.idx = order(fitn, decreasing = FALSE)
      fitn.ordered = fitn[fitn.ordered.idx]

      # update best solution so far
      valid = (penalties == 0)
      if (any(valid)) {
        #stopifnot(all(fitn[valid] == fitn.repaired[valid]))
        #stopifnot(all(arx[, valid, drop = FALSE] == arx.repaired[, valid, drop = FALSE]))
        min.valid.idx = which.min(fitn.repaired[valid])
        if (fitn.repaired[valid][min.valid.idx] < best.fitness) {
          best.fitness = fitn.repaired[valid][min.valid.idx]
          best.param = arx[, valid, drop = FALSE][, min.valid.idx]
        }
      }

      # update mean value / center of mass
      new.pop.idx = fitn.ordered.idx[1:mu]
      x.best = arx[, new.pop.idx, drop = FALSE]
      m.old = m
      m = drop(x.best %*% weights)

      y.best = ary[, new.pop.idx, drop = FALSE]
      y.w = drop(y.best %*% weights)
      z.best = arz[, new.pop.idx, drop = FALSE]
      z.w = drop(z.best %*% weights)

      # log population
      if (log.population) {
        population.trace[[iter]] = x.best
      }

  		# Update evolution path with cumulative step-size Adaptation (CSA) / path length control
      # For an explanation of the last factor see appendix A in https://www.lri.fr/~hansen/cmatutorial.pdf
      ps = (1 - cs) * ps + sqrt(cs * (2 - cs) * mu.eff) * (Cinvsqrt %*% y.w)
  		h.sigma = as.integer(norm2(ps) / sqrt(1 - (1 - cs)^(2 * (iter + 1))) < chi.n * (1.4 + 2 / (n + 1)))

  		# Update covariance matrix
      pc = (1 - cc) * pc + h.sigma * sqrt(cc * (2 - cc) * mu.eff) * y.w
      y = BD %*% z.best
      delta.h.sigma = as.numeric((1 - h.sigma) * cc * (2 - cc) <= 1)
  		C = (1 - c1 - cmu) * C + c1 * (pc %*% t(pc) + delta.h.sigma * C) + cmu * y %*% diag(weights) %*% t(y)

      # Update step-size sigma
      sigma = sigma * exp(cs / ds * ((norm2(ps) / chi.n) - 1))

      # Finally do decomposition C = B D^2 B^T
      e = eigen(C, symmetric = TRUE)
      B = e$vectors
      D = diag(sqrt(e$values))
      BD = B %*% D
      Cinvsqrt = B %*% diag(1 / diag(D)) %*% t(B) # update C^-1/2

      callMonitor(monitor, "step")

      # escape flat fitness values
      if (fitn.ordered[1L] == fitn.ordered[ceiling(0.7 * lambda)]) {
        sigma = sigma * exp(0.2 + cs / ds)
        if (!is.null(monitor)) {
          warningf("Flat fitness values; increasing mutation step-size. Consider reformulating the objective!")
        }
      }

      # CHECK STOPPING CONDITIONS
      # =========================
      stop.obj = checkStoppingConditions(stop.ons)

      n.stop.codes = length(stop.obj$codes)
      if (max.restarts > 0L && any(stop.obj$codes %in% restart.triggers)) {
        if (!is.null(monitor)) {
          messagef("Restart trigger fired! Restarting!!!")
        }
        n.stop.codes = sum(!(stop.obj$codes %in% restart.triggers))
        restarting = TRUE
      }

      # check if CMA-ES should really quit, i.e., is there a stopping condition,
      # that is active and does not trigger a restart?
      if (!restarting && (n.stop.codes > 0L)) {
        do.terminate = TRUE
        break
      }
  	}

    # really quit without more restarts
    if (do.terminate) {
      break
    }
  }

  callMonitor(monitor, "after")

	makeS3Obj(
		par.set = par.set,
		best.param = best.param,
		best.fitness = best.fitness,
		n.evals = n.evals,
		past.time = as.integer(difftime(Sys.time(), start.time, units = "secs")),
		n.iters = iter - 1L,
    n.restarts = run,
    population.trace = population.trace,
		message = stop.obj$stop.msgs,
		classes = "cma_result"
	)
}

#' @export
print.cma_result = function(x, ...) {
	best.param = list(x$best.param)
	names(best.param) = getParamIds(x$par.set)
	catf("Best parameter      : %s", paramValueToString(x$par.set, best.param))
	catf("Best fitness value  : %.6g", x$best.fitness)
	catf("Termination         : %s", x$message)
	catf("  #Iterations       : %i", x$n.iters)
	catf("  #Evaluations      : %i", x$n.evals)
	catf("  Time              : %i (secs)", x$past.time)
}