Nothing
#' @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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.