TAopt: Optimisation with Threshold Accepting

Description

The function implements the Threshold Accepting algorithm.

Usage

TAopt(OF, algo = list(), ...)

Value

TAopt returns a list with four components:

xbest: the solution
OFvalue: objective function value of the solution, ie, OF(xbest, ...)
Fmat: if algo$storeF is TRUE, a matrix with one row for each iteration (excluding the initial algo$nD steps) and two columns. The first column contains the objective function values of the neighbour solution at a given iteration; the second column contains the value of the current solution. Since TA can walk away from locally-optimal solutions, the best solution can be monitored through cummin(Fmat[ ,2L]).
xlist: if algo$storeSolutions is TRUE, a list; else NA. Contains the neighbour solutions at a given iteration (xn) and the current solutions (xc). Example: Fmat[i, 2L] is the objective function value associated with xlist[[c(2L, i)]].
initial.state: the value of .Random.seed when the function was called.

If algo$classify was set to TRUE, the resulting list will have a class attribute TAopt.

Arguments

OF: The objective function, to be minimised. Its first argument needs to be a solution x; it will be called as OF(x, ...).
algo: A list of settings for the algorithm. See Details.
...: other variables passed to OF and algo$neighbour. See Details.

Author

Enrico Schumann

Details

Threshold Accepting (TA) changes an initial solution iteratively; the algorithm stops after a fixed number of iterations. Conceptually, TA consists of a loop than runs for a number of iterations. In each iteration, a current solution xc is changed through a function algo$neighbour. If this new (or neighbour) solution xn is not worse than xc, ie, if OF(xn,...) <= OF(xc,...), then xn replaces xc. If xn is worse, it still replaces xc as long as the difference in ‘quality’ between the two solutions is less than a threshold tau; more precisely, as long as OF(xn,...) - tau <= OF(xc,...). Thus, we also accept a new solution that is worse than its predecessor; just not too much worse. The threshold is typically decreased over the course of the optimisation. For zero thresholds TA becomes a stochastic local search.

The thresholds can be passed through the list algo (see below). Otherwise, they are automatically computed through the procedure described in Gilli et al. (2006). When the thresholds are created automatically, the final threshold is always zero.

The list algo contains the following items.

nS: The number of steps per threshold. The default is 1000; but this setting depends very much on the problem.
nT: The number of thresholds. Default is 10; ignored if algo$vT is specified.
nI
nD: The number of random steps to compute the threshold sequence. Defaults to 2000. Only used if algo$vT is NULL.
q: The highest quantile for the threshold sequence. Defaults to 0.5. Only used if algo$vT is NULL. If q is zero, TAopt will run with algo$nT zero-thresholds (ie, like a Local Search).
x0: The initial solution. If this is a function, it will be called once without arguments to compute an initial solution, ie, x0 <- algo$x0(). This can be useful when the routine is called in a loop of restarts, and each restart is to have its own starting value.
vT: The thresholds. A numeric vector. If NULL (the default), TAopt will compute algo$nT thresholds. Passing threshold can be useful when similar problems are handled. Then the time to sample the objective function to compute the thresholds can be saved (ie, we save algo$nD function evaluations). If the thresholds are computed and algo$printDetail is TRUE, the time required to evaluate the objective function will be measured and an estimate for the remaining computing time is issued. This estimate is often very crude.
neighbour: The neighbourhood function, called as neighbour(x, ...). Its first argument must be a solution x; it must return a changed solution.
printDetail: If TRUE (the default), information is printed. If an integer i greater then one, information is printed at very ith iteration.
printBar: If TRUE (default is FALSE), a txtProgressBar (from package utils) is printed. The progress bar is not shown if printDetail is an integer greater than 1.
scale: The thresholds are multiplied by scale. Default is 1.
drop0: When thresholds are computed, should zero values be dropped from the sample of objective-function values? Default is FALSE.
stepUp: Defaults to 0. If an integer greater than zero, then the thresholds are recycled, ie, vT is replaced by rep(vT, algo$stepUp + 1) (and the number of thresholds will be increased by algo$nT times algo$stepUp). This option works for supplied as well as computed thresholds. Practically, this will have the same effect as restarting from a returned solution. (In Simulated Annealing, this strategy goes by the name of ‘reheating’.)
thresholds.only: Defaults to FALSE. If TRUE, compute only threshold sequence, but do not actually run TA.
storeF: if TRUE (the default), the objective function values for every solution in every generation are stored and returned as matrix Fmat.
storeSolutions: Default is FALSE. If TRUE, the solutions (ie, decision variables) in every generation are stored and returned in list xlist (see Value section below). To check, for instance, the current solution at the end of the ith generation, retrieve xlist[[c(2L, i)]].
classify: Logical; default is FALSE. If TRUE, the result will have a class attribute TAopt attached. This feature is experimental: the supported methods (plot, summary) may change without warning.
OF.target: Numeric; when specified, the algorithm will stop when an objective-function value as low as OF.target (or lower) is achieved. This is useful when an optimal objective-function value is known: the algorithm will then stop and not waste time searching for a better solution.

At the minimum, algo needs to contain an initial solution x0 and a neighbour function.

The total number of iterations equals algo$nT times (algo$stepUp + 1) times algo$nS (plus possibly algo$nD).

References

Dueck, G. and Scheuer, T. (1990) Threshold Accepting. A General Purpose Optimization Algorithm Superior to Simulated Annealing. Journal of Computational Physics. 90 (1), 161--175.

Dueck, G. and Winker, P. (1992) New Concepts and Algorithms for Portfolio Choice. Applied Stochastic Models and Data Analysis. 8 (3), 159--178.

Gilli, M., Këllezi, E. and Hysi, H. (2006) A Data-Driven Optimization Heuristic for Downside Risk Minimization. Journal of Risk. 8 (3), 1--18.

Gilli, M., Maringer, D. and Schumann, E. (2019) Numerical Methods and Optimization in Finance. 2nd edition. Elsevier. tools:::Rd_expr_doi("10.1016/C2017-0-01621-X")

Moscato, P. and Fontanari, J.F. (1990). Stochastic Versus Deterministic Update in Simulated Annealing. Physics Letters A. 146 (4), 204--208.

Schumann, E. (2012) Remarks on 'A comparison of some heuristic optimization methods'. https://enricoschumann.net/R/remarks.htm

Schumann, E. (2023) Financial Optimisation with R (NMOF Manual). https://enricoschumann.net/NMOF.htm#NMOFmanual

Winker, P. (2001). Optimization Heuristics in Econometrics: Applications of Threshold Accepting. Wiley.

Examples

Run this code

## Aim: given a matrix x with n rows and 2 columns,
##      divide the rows of x into two subsets such that
##      in one subset the columns are highly correlated,
##      and in the other lowly (negatively) correlated.
##      constraint: a single subset should have at least 40 rows

## create data with specified correlation
n <- 100L
rho <- 0.7
C <- matrix(rho, 2L, 2L); diag(C) <- 1
x <- matrix(rnorm(n * 2L), n, 2L) %*% chol(C)

## collect data
data <- list(x = x, n = n, nmin = 40L)

## a random initial solution
x0 <- runif(n) > 0.5

## a neighbourhood function
neighbour <- function(xc, data) {
    xn <- xc
    p <- sample.int(data$n, size = 1L)
    xn[p] <- abs(xn[p] - 1L)
    # reject infeasible solution
    c1 <- sum(xn) >= data$nmin
    c2 <- sum(xn) <= (data$n - data$nmin)
    if (c1 && c2) res <- xn else res <- xc
    as.logical(res)
}

## check (should be 1 FALSE and n-1 TRUE)
x0 == neighbour(x0, data)

## objective function
OF <- function(xc, data)
    -abs(cor(data$x[xc, ])[1L, 2L] - cor(data$x[!xc, ])[1L, 2L])

## check
OF(x0, data)
## check
OF(neighbour(x0, data), data)

## plot data
par(mfrow = c(1,3), bty = "n")
plot(data$x,
     xlim = c(-3,3), ylim = c(-3,3),
     main = "all data", col = "darkgreen")

## *Local Search*
algo <- list(nS = 3000L,
             neighbour = neighbour,
             x0 = x0,
             printBar = FALSE)
sol1 <- LSopt(OF, algo = algo, data=data)
sol1$OFvalue

## *Threshold Accepting*
algo$nT <- 10L
algo$nS <- ceiling(algo$nS/algo$nT)
sol <- TAopt(OF, algo = algo, data = data)
sol$OFvalue

c1 <- cor(data$x[ sol$xbest, ])[1L, 2L]
c2 <- cor(data$x[!sol$xbest, ])[1L, 2L]

lines(data$x[ sol$xbest, ], type = "p", col = "blue")

plot(data$x[ sol$xbest, ], col = "blue",
     xlim = c(-3,3), ylim = c(-3,3),
     main = paste("subset 1, corr.", format(c1, digits = 3)))

plot(data$x[!sol$xbest, ], col = "darkgreen",
     xlim = c(-3,3), ylim = c(-3,3),
     main = paste("subset 2, corr.", format(c2, digits = 3)))

## compare LS/TA
par(mfrow = c(1,1), bty = "n")
plot(sol1$Fmat[ ,2L],type="l", ylim=c(-1.5,0.5),
     ylab = "OF", xlab = "iterations")
lines(sol$Fmat[ ,2L],type = "l", col = "blue")
legend(x = "topright",legend = c("LS", "TA"),
       lty = 1, lwd = 2,col = c("black", "blue"))

Run the code above in your browser using DataLab