Update vignette

1 files changed,134insertions(+),103deletions(-) M vignettes/neighbours.Rnw

M vignettes/neighbours.Rnw => vignettes/neighbours.Rnw +134 -103

@@ 40,7 40,7 @@ pv <- gsub("(.*)[.](.*)", "\\1-\\2", pv)\noindent \url{es@enricoschumann.net}\\ \bigskip \noindent Local-search algorithms proceed through the \noindent Local-search algorithms move through the solution space by randomly selecting new solutions that are close to the original solutions. Such a new solution is called a neighbour solution. The function@@ 57,8 57,9 @@ set.seed(34734)\noindent In the examples that follow, we will use a simple optimisation algorithm, a (stochastic) Local Search. If package \texttt{NMOF} is available, function \texttt{LSopt} is used; otherwise, we use a simple replacement (taken from \citealp[Chapter~13]{Gilli2019}). function \texttt{LSopt} from that package is used; otherwise, we use a simple replacement (taken from \citealp[Chapter~13]{Gilli2019}). <<LSopt>>= LSopt <- if (requireNamespace("NMOF")) {@@ 79,17 80,17 @@ LSopt <- if (requireNamespace("NMOF")) {} @ \section*{Selecting elements of a list} \section*{Example: Selecting elements of a list} We are given a numeric vector~$y$, and also a matrix~$X$, which has as many rows as $y$ has elements. matrix~$X$, which has as many rows as~$y$ has elements. The aim now is to find a subset of columns of~$X$ whose average is as weakly correlated with $y$ as possible. average is as lowly correlated with $y$ as possible. Let us create random data. <<data>>= ny <- 50 nx <- 1000 nx <- 2000 y <- runif(ny) X <- array(runif(ny * nx), dim = c(ny, nx)) @@@ 110,21 111,28 @@ x0[1:15] <- TRUEhead(x0, 20) @ We write an objective function to compute the actual \noindent We write an objective function to compute the actual quality of the solution \texttt{x0}. <<objective-fun>>= column_cor <- function(x, X, y) cor(rowMeans(X[, x]), y) @ We evaluate the solution. \noindent We evaluate the quality of~\texttt{x0}. <<>>= column_cor(x0, X, y) @ To run a Local Search, we need a neighbourhood function. Calling \texttt{neighbourfun} will create such a function, and it will bind the parameters to this function. \noindent To run a Local Search, we need a neighbourhood function. Calling \texttt{neighbourfun} will create such a function, and it will bind the parameters to this function.% \marginpar{\footnotesize\RaggedRight \citet{Gentleman2000} show that \textsf{R}'s scoping rules are particularly convenient for optimisation.} % <<nb>>= nb <- neighbourfun(type = "logical", kmin = 10, kmax = 20) @@@ 155,7 163,7 @@ par(mfrow = c(1, 2), las = 1, bty = "n",plot(y, rowMeans(X[, x0]), main = "Initial solution", pch = 19, cex = 0.5, ylim = c(0.2, 0.8), ylim = c(0.3, 0.7), ylab = "Linear combination of columns") par(yaxt = "n") plot(y, rowMeans(X[, sol.ls$xbest]),@@ 170,130 178,166 @@ axis(4)\section*{More restrictions} The neighbourhood function we used in the previous section included constraints: it would not include less than 10 or more than 20~\texttt{TRUE} values. Note that the neighbourhood function required a valid \texttt{x} as input.% section included constraints: it would include no fewer than 10 or no more than 20~\texttt{TRUE} values. Note that the neighbourhood function required a valid \texttt{x} as input.% \marginpar{\footnotesize\RaggedRight For invalid \texttt{x}, the result is undefined. Neighbourhood functions should not check the validity of their inputs, because of speed: } % inputs, because of speed: the functions are called thousands of times during an optimisation run, and every fraction of second matters.} % We could also set \texttt{kmin} and \texttt{kmax} to the same integer, so the number of \texttt{TRUE} values is fixed. (In this case, a slightly-different We may also set \texttt{kmin} and \texttt{kmax} to the same integer, so the number of \texttt{TRUE} values is fixed. (In this case, a slightly-different neighbourhood algorithm will be used.) <<fixed>>= x <- logical(8) x[1:3] <- TRUE nb <- neighbourfun(type = "logical", kmin = 3, kmax = 3) @ % \noindent Let us take a few random steps. % <<echo=false>>= % for (i in 1:10) { % if (i == 1) % cat(" 0 ", ifelse(x, "o", "."), % " | initial solution: o == TRUE, . == FALSE ", % sep = "", fill = TRUE) % x <- nb(x) % cat(format(i, width = 2), " ", % ifelse(x, "o", "."), sep = "", fill = TRUE) % } % @ \noindent We can also add a constraint about elements not to touch. Suppose the initial solution is the following: We can also add a constraint about elements not to touch. Suppose the initial solution is the following: <<>>= x <- logical(9) x[4:6] <- TRUE cat(ifelse(x, "o", "."), sep = "") compare_vectors(x) @ \noindent We restrict the changes that can be made to the solution: the first three elements must not be touched; only the remaining elements may change. (They are \texttt{active}.) <<restrict>>= active <- !logical(length(x)) active[1:3] <- FALSE active nb <- neighbourfun(type = "logical", kmin = 3, kmax = 3, active = active) @ are \texttt{active}.) <<restrict>>= active <- !logical(length(x)) active[1:3] <- FALSE active nb <- neighbourfun(type = "logical", kmin = 3, kmax = 3, active = active) @ \noindent Again, let us take a few random steps. The element in the middle remains \texttt{TRUE}, just as the first and last element remain \texttt{FALSE}. <<echo=false>>= xs <- list() xs[[1]] <- x for (i in 1:15) { if (i == 1) cat(" 0 ", ifelse(x, "o", "."), " | initial solution: o == TRUE, . == FALSE ", sep = "", fill = TRUE) x <- nb(x) cat(format(i, width = 2), " ", ifelse(x, "o", "."), sep = "", fill = TRUE) xs[[length(xs) + 1]] <- x <- nb(x) } do.call(compare_vectors, xs) @ \section*{Numeric solutions} \section*{Another example: minimising portfolio risk} When the argument \texttt{type} is numeric, the resulting neighbour function will randomly select elements and then add and subtract random numbers from those elements. When the argument \texttt{type} is set to \texttt{numeric}, the resulting neighbour function will randomly select elements, and then add and subtract typically-small numbers to and from those elements. The size of those numbers is controlled by argument \texttt{stepsize}, and with argument \texttt{random} set to TRUE (the default), the step sizes will vary randomly. <<>>= x <- rep(1, 5) nb <- neighbourfun(min = -Inf, max = Inf, type = "numeric", stepsize = 0.05, sum = FALSE) for (i in 1:5) { print(c(Sum = sum(x),x= x)) x <- nb(x) } @ \noindent The default for \texttt{sum} is actually TRUE. Suppose we are given a matrix \texttt{R} and aim to find a vector $x$ such that the variance of the elements in~$Rx$ is minimised. This is common problem in finance, in which $R$ could be a matrix of asset-price returns, with each column holding the returns of one asset, and $x$ a vector of portfolio weights. <<>>= x <- rep(1, 5) nb <- neighbourfun(min = -Inf, max = Inf, type = "numeric", stepsize = 0.05) For this particular goal -- the variance of returns --, we can first compute the variance--covariance matrix~$\Sigma$ of $R$, and then minimise $x' \Sigma x$. for (i in 1:5) { print(c(Sum = sum(x),x= x)) x <- nb(x) } <<>>= R <- NMOF::randomReturns(20, 120, 0.03, rho = 0.6) S <- cov(R) ## Sigma @ This setting is the default, because it is useful in many applications. Suppose we start with an equal-weight portfolio. <<>>= x0 <- rep(1/20, 20) @ \section*{An application: portfolio optimisation} The objective function. <<>>= variance <- function(x, S, ...) x %*% S %*% x variance(x0, S) @ We can solve this problem with Local Search. <<mv-ls1>>= nb <- neighbourfun(type = "numeric", min = 0, max = 0.10, stepsize = 0.005) sol.ls <- LSopt(variance, list(x0 = x0, neighbour = nb, nI = 1000, printBar = FALSE), S = S) Suppose we are given NMOF::minvar(S, wmin = 0, wmax = 0.1) @ Minimise variance An alternative, and perhaps more-straightforward way to write the objective function is the following. <<>>= variance2 <- function(x, R, ...) var(R %*% x) @ When we feed it to \texttt{LSopt}, we arrive at the same solution. <<mv-ls2>>= nb <- neighbourfun(type = "numeric", min = 0, max = 0.10, stepsize = 0.005) sol.ls <- LSopt(variance2, list(x0 = x0, neighbour = nb, nI = 1000, printBar = FALSE), R = R) NMOF::minvar(S, wmin = 0, wmax = 0.1) @ This second objective function is less efficient than the first. But it is much more flexible: only for minimising variance, could we take the shortcut via the variance--covariance matrix. But for other measures of risk, we cannot do that. One example is the so-called semi-variance, defined as \begin{equation} \label{eq:sv} \frac{1}{N}\sum_{i=1}^N \min(R_ix - m, 0)^2 \end{equation} All we have to do now is to exchange objective functions. <<>>= semivariance <- function(x, R, ...) { Rx <- R %*% x Rx.lower <- pmin(Rx - mean(Rx), 0) sum(Rx.lower) } @ \section*{Updating} In many applications, notably when doing data analysis, the solution \texttt{x} is used to compute a matrix/vector product $Ax$, in which $A$ is a $m$ times $n$~matrix, and $x$ is of length$n$. In many applications, such as the example in the previous section, but also in many other cases when doing data analysis, the solution \texttt{x} is used to compute a matrix/vector product $Ax$, in which $A$ is a $m$ times $n$~matrix, and $x$ is of length$n$. If we only change few elements in the solution, we can@@ 343,20 387,7 @@ data.frame(target[1:10], sol$xbest[1:10])@ <<>>= x <- runif(5)+0:4 nb <- neighbourfun(type = "numeric", stepsize = 0.1, min = 0, max = 50) plot(rep(1, 5), x, xlim = c(0, 100)) for (i in 2:100) { x <- nb(x) lines(rep(i, 5), x, type = "p") } @ \citet{Gentleman2000} \citet{Gilli2019} \bibliographystyle{plainnat}