Package adagio

May 29, 2016

Type Package
Title Discrete and Global Optimization Routines
Version 0.6.5
Date 2016-05-29
Author Hans Werner Borchers
Maintainer Hans W. Borchers <[email protected]>
Depends R (>= 3.0.0)
Imports graphics, stats
Description The R package 'adagio' will provide methods and algorithms for
discrete optimization and (evolutionary) global optimization.
License GPL (>= 3)
LazyLoad yes
LazyData yes
Repository CRAN
Repository/R-Forge/Project optimist
Repository/R-Forge/Revision 408
Repository/R-Forge/DateTimeStamp 2016-05-29 12:20:55
Date/Publication 2016-05-29 21:14:02
NeedsCompilation yes

R topics documented:
adagio-package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
CMAES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
fminviz,flineviz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
knapsack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
maxempty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
maxquad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1
2 adagio-package

maxsub . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
mknapsack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
nelmin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
SIAM test functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
simpleDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
simpleEA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
subsetsum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Testfunctions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
transfinite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Index 29

adagio-package Discrete and Global Optimization

Description
The R package adagio will provide methods and algorithms for discrete as well as (derivative-free)
global optimization.

Details

Package: adagio
Type: Package
Version: 0.6.4
Date: 2015-12-07
License: GPL (>= 3)
LazyLoad: yes

Author(s)
Hans W Borchers
Maintainer: HwB <[email protected]>

References
Arnold Neumaier (2010). Global Optimization. http://www.mat.univie.ac.at/~neum/glopt.
html

See Also
The R-Forge project Optimization and solving packages. https://r-forge.r-project.org/
projects/optimizer/
assignment 3

Examples
# See the functions in this package.

assignment Linear Sum Assignment Problem

Description
Linear (sum) assignment problem, or LSAP.

Usage
assignment(cmat)

Arguments
cmat quadratic integer matrix, the cost matrix.

Details
Solves the linear (sum) assignment problem for quadratic matrices with integer entries.

Value
List with components perm, the permutation that defines the minimum solution, min, the minimum
value, and err, which is -1 if an integer overflow occured.

Note
Faster than the Hungarian algorithm, but only applicable to quadratic cost matrices with integer
components.

Author(s)
Copyright(c) 1993 A. H. Morris, Jr., Naval Surface Warfare Center, using Fortran routines written
by S. Martello and P. Toth, University of Bologna. Released for free and general use, now under
GPL license, wrapped for R by Hans W Borchers <[email protected]>.

References
Burkard, R., M. DellAmico, and S. Martello (2009). Assignment Problems. Society for Industrial
and Applied Mathematics (SIAM).
Martello, S., and P. Toth (1990). Knapsack Problems: Algorithms and Computer Implementations.
John Wiley & Sons, Ltd.

See Also
clue::solve_LSAP
4 CMAES

Examples
## Example similar to clue::solve_LSAP
set.seed(8237)
x <- matrix(sample(1:100), nrow = 10)
y <- assignment(x)
# show permutation and check minimum sum
y$perm; y$t # 4 5 7 2 6 1 3 8 10 9
z <- cbind(1:10, y$perm) # 156
x[z] # 5 4 11 8 20 7 38 15 22 26
y$min == sum(x[z]) # TRUE

## Not run:
## Example: minimize sum of distances of complex points
n <- 100
x <- rt(n, df=3) + 1i * rt(n, df=3)
y <- runif(n) + 1i * runif(n)
cmat <- round(outer(x, y, FUN = function(x,y) Mod(x - y)), 2)
dmat <- round(100*cmat, 2)
system.time(T1 <- assignment(dmat)) # elapsed: 0.003
T1$min / 100 # 139.32

library("clue")
system.time(T2 <- solve_LSAP(cmat)) # elapsed: 0.034
sum(cmat[cbind(1:n, T2)]) # 139.32

## End(Not run)

CMAES Covariance Matrix Adaptation Evolution Strategy

Description
The CMA-ES (Covariance Matrix Adaptation Evolution Strategy) is an evolutionary algorithm for
difficult non-linear non-convex optimization problems in continuous domain. The CMA-ES is typ-
ically applied to unconstrained or bounded constraint optimization problems, and search space di-
mensions between three and fifty.

Usage
pureCMAES(par, fun, lower = NULL, upper = NULL, sigma = 0.5,
stopfitness = -Inf, stopeval = 1000*length(par)^2, ...)

Arguments
par objective variables initial point.
fun objective/target/fitness function.
lower,upper lower and upper bounds for the parameters.
sigma coordinate wise standard deviation (step size).
CMAES 5

stopfitness stop if fitness < stopfitness (minimization).

stopeval stop after stopeval number of function evaluations
... additional parameters to be passed to the function.

Details
The CMA-ES implements a stochastic variable-metric method. In the very particular case of a
convex-quadratic objective function the covariance matrix adapts to the inverse of the Hessian ma-
trix, up to a scalar factor and small random fluctuations. The update equations for mean and covari-
ance matrix maximize a likelihood while resembling an expectation-maximization algorithm.

Value
Returns a list with components xmin and fmin.
Be patient; for difficult problems or high dimensions the function may run for several minutes;
avoid problem dimensions of 30 and more!

Note
There are other implementations of Hansens CMAES in package cmaes (simplified form) and in
package parma as cmaes() (extended form).

Author(s)
Copyright (c) 2003-2010 Nikolas Hansen for Matlab code PURECMAES; converted to R by Hans
W Borchers.

References
Hansen, N., and A. Ostermeier (2001). Completely Derandomized Self-Adaptation in Evolution
Strategies. Evolutionary Computation 9(2), pp. 159-195.
URL: https://www.lri.fr/~hansen/cmaartic.pdf
Hansen, N., S.D. Mueller, P. Koumoutsakos (2003). Reducing the Time Complexity of the Deran-
domized Evolution Strategy with Covariance Matrix Adaptation (CMA-ES). Evolutionary Compu-
tation 11(1), pp. 1-18.
URL: https://www.lri.fr/~hansen/evco_11_1_1_0.pdf
Hansen, N. (2011). The CMA Evolution Strategy: A Tutorial.
URL: https://www.lri.fr/~hansen/cmatutorial.pdf
Hansen, N., D.V. Arnold, and A. Auger (2013). Evolution Strategies. To appear in Janusz Kacprzyk
and Witold Pedrycz (Eds.): Handbook of Computational Intelligence, Springer-Verlag (accepted for
publication).
URL: https://www.lri.fr/~hansen/es-overview-2014.pdf

See Also
cmaes::cmaes, parma::cmaes
6 fminviz,flineviz

Examples
## Not run:
## Polynomial minimax approximation of data points
## (see the Remez algorithm)
n <- 10; m <- 101 # polynomial of degree 10; no. of data points
xi <- seq(-1, 1, length = m)
yi <- 1 / (1 + (5*xi)^2) # Runge's function

pval <- function(p, x) # Horner scheme

outer(x, (length(p) - 1):0, "^") %*% p

pfit <- function(x, y, n) # polynomial fitting of degree n

qr.solve(outer(x, seq(n, 0), "^"), y)

fn1 <- function(p) # objective function

max(abs(pval(p, xi) - yi))

pf <- pfit(xi, yi, 10) # start with a least-squares fitting

sol1 <- pureCMAES(pf, fn1, rep(-200, 11), rep(200, 11))
zapsmall(sol1$xmin)
# [1] -50.24826 0.00000 135.85352 0.00000 -134.20107 0.00000
# [7] 59.19315 0.00000 -11.55888 0.00000 0.93453

print(sol1$fmin, digits = 10)

# [1] 0.06546780411

## Polynomial fitting in the L1 norm

## (or use LP or IRLS approaches)
fn2 <- function(p)
sum(abs(pval(p, xi) - yi))

sol2 <- pureCMAES(pf, fn2, rep(-100, 11), rep(100, 11))

zapsmall(sol2$xmin)
# [1] -21.93238 0.00000 62.91083 0.00000 -67.84847 0.00000
# [7] 34.14398 0.00000 -8.11899 0.00000 0.84533

print(sol2$fmin, digits = 10)

# [1] 3.061810639

## End(Not run)

fminviz,flineviz Visualize Function Minimum

Description

Visualizes multivariate functions around a point or along a line between two points in R^n.
hamiltonian 7

Usage
fminviz(fn, x0, nlines = 2*length(x0),
npoints = 51, scaled = 1.0)

flineviz(fn, x1, x2, npoints = 51, scaled = 0.1)

Arguments
fn multivariate function to be visualized.
x0,x1,x2 points in n-dimensional space.
nlines number of lines to plot.
npoints number of points used to plot a line.
scaled scale factor to extend the line(s).

Details
fminviz vizualizes the behavior of a multivariate function fn around a point x0. It randomly selects
nlines lines through x0 in R^n and draws the curves of the function along these lines in one graph.
Curves that have at least one point below fn(x0) are drawn in red, all others in blue. The scale on
the x-axis is the Euclidean distance in R^n. The scale factor can change it.
flineviz vizualizes the behavior of a multivariate function fn along the straight line between the
points x1 and x2. Points x1 and x2 are also plotted.

Value
Plots a line graph and returns NULL (invisibly).

Examples
## Not run:
f1 <- function(x) x[1]^2 - x[2]^2
fminviz(f1, c(0, 0), nlines = 10)

f2 <- function(x) (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2

flineviz(f2, c(0, 0), c(1, 1))

## End(Not run)

hamiltonian Finds a Hamiltonian path or cycle

Description
A Hamiltionian path or cycle (a.k.a. Hamiltonian circuit) is a path through a graph that visits each
vertex exactly once, resp. a closed path through the graph.
8 hamiltonian

Usage
hamiltonian(edges, start = 1, cycle = TRUE)

Arguments
edges an edge list describing an undirected graph.
start vertex number to start the path or cycle.
cycle Boolean, should a path or a full cycle be found.

Details
hamiltonian() applies a backtracking algorithm that is relatively efficient for graphs of up to 30
40 vertices. The edge list is first transformed to a list where the i-th component contains the list of
all vertices connected to vertex i.
The edge list must be of the form c(v1, v2, v3, v2, ...) meaning that there are edges
v1 --> v2, v3 --> v4, etc., connecting these vertices. Therefore, an edge list has an even
number of entries.
If the function returns NULL, there is no Hamiltonian path or cycle. The function does not check if
the graph is connected or not. And if cycle = TRUE is used, then there also exists an edge from the
last to the first entry in the resulting path.
Ifa Hamiltonian cycle exists in the graph it will be found whatever the starting vertex was. For a
Hamiltonian path this is different and a successful search may very well depend on the start.

Value
Returns a vector containing vertex number of a valid path or cycle, or NULL if no path or cycle has
been found (i.e., does not exist); If a cycle was requested, there exists an edge from the last to the
first vertex in this list of edges.

Note
See the igraph package for more information about handling graphs and defining them through
edge lists or other constructs.

Author(s)
Hans W. Borchers

References
Papadimitriou, Ch. H., and K. Steiglitz (1998). Optimization Problems: Algorithms and Complex-
ity. Prentice-Hall/Dover Publications.

See Also
Package igraph
hamiltonian 9

Examples
## Dodekaeder graph
D20_edges <- c(
1, 2, 1, 5, 1, 6, 2, 3, 2, 8, 3, 4, 3, 10, 4, 5, 4, 12,
5, 14, 6, 7, 6, 15, 7, 8, 7, 16, 8, 9, 9, 10, 9, 17, 10, 11,
11, 12, 11, 18, 12, 13, 13, 14, 13, 19, 14, 15, 15, 20, 16, 17, 16, 20,
17, 18, 18, 19, 19, 20)
hamiltonian(D20_edges, cycle = TRUE)
# [1] 1 2 3 4 5 14 13 12 11 10 9 8 7 16 17 18 19 20 15 6
hamiltonian(D20_edges, cycle = FALSE)
# [1] 1 2 3 4 5 14 13 12 11 10 9 8 7 6 15 20 16 17 18 19

## Herschel graph
# The Herschel graph the smallest non-Hamiltonian polyhedral graph.
H11_edges <- c(
1, 2, 1, 8, 1, 9, 1, 10, 2, 3, 2, 11, 3, 4, 3, 9, 4, 5,
4, 11, 5, 6, 5, 9, 5, 10, 6, 7, 6, 11, 7, 8, 7, 10, 8, 11)
hamiltonian(H11_edges, cycle = FALSE)
# NULL

## Not run:
## Example: Graph constructed from squares
N <- 45 # 23, 32, 45
Q <- (2:trunc(sqrt(2*N-1)))^2
sq_edges <- c()
for (i in 1:(N-1)) {
for (j in (i+1):N) {
if ((i+j)
sq_edges <- c(sq_edges, i, j)
}
}

require(igraph)
sq_graph <- make_graph(sq_edges, directed=FALSE)
plot(sq_graph)

if (N == 23) {
# does not find a path with start=1 ...
hamiltonian(sq_edges, start=18, cycle=FALSE)
# hamiltonian(sq_edges) # NULL
} else if (N == 32) {
# the first of these graphs that is Hamiltonian ...
# hamiltonian(sq_edges, cycle=FALSE)
hamiltonian(sq_edges)
} else if (N == 45) {
# takes much too long ...
# hamiltonian(sq_edges, cycle=FALSE)
hamiltonian(sq_edges)
}
## End(Not run)
10 knapsack

knapsack 0-1 Knapsack Problem

Description
Solves the 0-1 (binary) single knapsack problem.

Usage
knapsack(w, p, cap)

Arguments
w vector of weights.
p vector of profits.
cap maximal capacity of the knapsack.

Details
knapsack solves the 0-1, or: binary, single knapsack problem by using the synamic programming
approach. The problem can be formulated as:
Maximize sum(x*p) such that sum(x*w) <= cap, where x is a vector with x[i] == 0 or 1.

Value
A list with components capacity, profit, and indices.

Note
Will be replaced by a compiled version.

Author(s)
HwB email: <[email protected]>

References
Papadimitriou, C. H., and K. Steiglitz (1998). Combinatorial Optimization: Algorithms and Com-
plexity. Dover Publications 1982, 1998.
Horowitz, E., and S. Sahni (1978). Fundamentals of Computer Algorithms. Computer Science
Press, Rockville, ML.

See Also
knapsack::knapsack
maxempty 11

Examples
# Example 1
p <- c(15, 100, 90, 60, 40, 15, 10, 1)
w <- c( 2, 20, 20, 30, 40, 30, 60, 10)
cap <- 102
(is <- knapsack(w, p, cap))
# [1] 1 2 3 4 6 , capacity 102 and total profit 280

## Example 2
p <- c(70, 20, 39, 37, 7, 5, 10)
w <- c(31, 10, 20, 19, 4, 3, 6)
cap <- 50
(is <- knapsack(w, p, cap))
# [1] 1 4 , capacity 50 and total profit 107

maxempty Maximally Empty Rectangle Problem

Description
Find the largest/maximal empty rectangle, i.e. with largest area, not containing given points.

Usage
maxempty(x, y, ax = c(0, 1), ay = c(0, 1))

Arguments
x, y coordinates of points to be avoided.
ax, ay left and right resp. lower and upper constraints.

Details
Find the largest or maximal empty two-dimensional rectangle in a rectangular area. The edges
of this rectangle have to be parallel to the edges of the enclosing rectangle (and parallel to the
coordinate axes). Empty means that none of the points given are contained in the interior of the
found rectangle.

Value
List with area and rect the rectangle as a vector usable for the rect graphics function.

Note
The algorithm has a run-time of O(n^2) while there are run-times of O(n*log(n)) reported in the
literature, utilizing a more complex data structure. I dont know of any comparable algorithms for
the largest empty circle problem.
12 maxquad

Author(s)
HwB email: <[email protected]>

References
B. Chazelle, R. L. Drysdale, and D. T. Lee (1986). Computing the Largest Empty Rectangle. SIAM
Journal of Computing, Vol. 15(1), pp. 300315.
A. Naamad, D. T. Lee, and W.-L. Hsu (1984). On the Maximum Empty Rectangle Problem. Dis-
crete Applied Mathematics, Vol. 8, pp. 267277.

See Also
Hmisc::largest.empty with a Fortran implementation of this code.

Examples
N <- 100; set.seed(8237)
x <- runif(N); y <- runif(N)
R <- maxempty(x, y, c(0,1), c(0,1))
R
# $area
# [1] 0.08238793
# $rect
# [1] 0.7023670 0.1797339 0.8175771 0.8948442

## Not run:
plot(x, y, pch="+", xlim=c(0,1), ylim=c(0,1), col="darkgray",
main = "Maximally empty rectangle")
rect(0, 0, 1, 1, border = "red", lwd = 1, lty = "dashed")
do.call(rect, as.list(R$rect))
grid()
## End(Not run)

maxquad The MAXQUAD Test Function

Description
Lemarechals MAXQUAD optimization test function.

Usage
maxquad(n, m)

Arguments
n number of variables of the generated test function.
m number of functions to compete for the maximum.
maxsub 13

Details
MAXQUAD actually is a family of minimax functions, parametrized by the number n of variables
and the number m of functions whose maximum it is.

Value
Returns a list with components fn the generated test function of n variables, and gr the correspond-
ing (analytical) gradient function.

References
Kuntsevich, A., and F. Kappel (1997). SolvOpt The Solver for Local Nonlinear Optimization
Problems. Manual Version 1.1, Institute of Mathematics, University of Graz.
Lemarechal, C., and R. Mifflin, Eds. (1978). Nonsmooth Optimization. Pergamon Press, Oxford.
Shor, N. Z. (1985). Minimization Methods for Non-differentiable Functions. Series in Computa-
tional Mathematics, Springer-Verlag, Berlin.

Examples
# Test function of 5 variables, defined as maximum of 5 smooth functions
maxq <- maxquad(5, 5)
fnMaxquad <- maxq$fn
grMaxquad <- maxq$gr
# shor

maxsub Maximal Sum Subarray

Description
Find a subarray with maximal positive sum.

Usage
maxsub(x, inds = TRUE, compiled = TRUE)

maxsub2d(A)

Arguments
x numeric vector.
A numeric matrix
inds logical; shall the indices be returned?
compiled logical; shall the compiled version be used?
14 maxsub

Details
maxsub finds a contiguous subarray whose sum is maximally positive. This is sometimes called
Kadanes algorithm.
maxsub will use a compiled and very fast version with a running time of O(n) where n is the length
of the input vector x.
maxsub2d finds a (contiguous) submatrix whose sum of elements is maximally positive. The ap-
proach taken here is to apply the one-dimensional routine to summed arrays between all rows of A.
This has a run-time of O(n^3), though a run-time of O(n^2 log n) seems possible see the reference
below.
maxsub2d uses a Fortran workhorse and can solve a 1000-by-1000 matrix in a few secondsbut
beware of biggere ones

Value
Either just a maximal sum, or a list this sum as component sum plus the start and end indices as a
vector inds.

Note
In special cases, the matrix A may be sparse or (as in the example section) only have one nonzero
element in each row and column. Expectation is that there may exists a more efficient (say O(n^2))
algorithm in this extreme case.

Author(s)
HwB <[email protected]>

References
Bentley, Jon (1986). Programming Pearls, Column 7. Addison-Wesley Publ. Co., Reading, MA.
T. Takaoka (2002). Efficient Algorithms for the Maximum Subarray Problem by Distance Matrix
Multiplication. The Australasian Theory Symposion, CATS 2002.

Examples
## Find a maximal sum subvector
set.seed(8237)
x <- rnorm(1e6)
system.time(res <- maxsub(x, inds = TRUE, compiled = FALSE))
res

## Standard example: Find a maximal sum submatrix

A <- matrix(c(0,-2,-7,0, 9,2,-6,2, -4,1,-4,1, -1,8,0,2),
nrow = 4, ncol = 4, byrow =TRUE)
maxsub2d(A)
# $sum: 15
# $inds: 2 4 1 2 , i.e., rows = 2..4, columns = 1..2

## Not run:
## Application to points in the unit square:
mknapsack 15

set.seed(723)
N <- 50; w <- rnorm(N)
x <- runif(N); y <- runif(N)
clr <- ifelse (w >= 0, "blue", "red")
plot(x, y, pch = 20, col = clr, xlim = c(0, 1), ylim = c(0, 1))

xs <- unique(sort(x)); ns <- length(xs)

X <- c(0, ((xs[1:(ns-1)] + xs[2:ns])/2), 1)
ys <- unique(sort(y)); ms <- length(ys)
Y <- c(0, ((ys[1:(ns-1)] + ys[2:ns])/2), 1)
abline(v = X, col = "gray")
abline(h = Y, col = "gray")

A <- matrix(0, N, N)
xi <- findInterval(x, X); yi <- findInterval(y, Y)
for (i in 1:N) A[yi[i], xi[i]] <- w[i]

msr <- maxsub2d(A)

rect(X[msr$inds[3]], Y[msr$inds[1]], X[msr$inds[4]+1], Y[msr$inds[2]+1])

## End(Not run)

mknapsack Multiple 0-1 Knapsack Problem

Description
Solves the 0-1 (binary) multiple knapsack problem.

Usage
mknapsack(p, w, k, bck = -1)

Arguments
p integer vector of profits.
w integer vector of weights.
k integer vector of capacities of the knapsacks.
bck maximal number of backtrackings allowed; default: -1.

Details
Solves the 0-1 multiple knapsack problem for integer profits and weights
A multiple 0-1 knapsack problem can be formulated as:
maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + p(n)*(x(1,n) + ... + x
subject to w(1)*x(i,1) + ... + w(n)*x(i,n) <= k(i) for i=1,...,m
x(1,j) + ... + x(m,j) <= 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , j=1,...,n ,
The input problem description must satisfy the following conditions:
16 mknapsack

vs=-1 if n < 2 or m < 1

vs=-2 if some p(j) , w(j) or k(i) are not positive
vs=-3 if a knapsack cannot contain any item
vs=-4 if an item cannot fit into any knapsack
vs=-5 if knapsack m contains all the items
vs=-7 if array k is not correctly sorted
vs=-8 [should not happen]

Value
A list with compomnents, ksack the knapsack numbers the items are assigned to, value the total
value/profit of the solution found, and bs the number of backtracks used.

Note
With some care, this function can be used for the bounded and unbounded single knapsack problem
as well.

Author(s)
The Fortran source code is adapted from the free NSCW Library of Mathematical Subroutines.
The wrapping code has been written by yours package maintainer,
HwB email: <[email protected]>

References
Kellerer, H., U. Pferschy, and D. Pisinger (2004). Knapsack Problems. Springer-Verlag, Berlin
Heidelberg.
Martello, S., and P. Toth (1990). Knapsack Problems: Algorithms and Computer Implementations.
John Wiley & Sons, Ltd.

See Also
Other packages implementing knapsack routines.

Examples
## Example 1: single knapsack
p <- c(15, 100, 90, 60, 40, 15, 10, 1)
w <- c( 2, 20, 20, 30, 40, 30, 60, 10)
cap <- 102
(is <- mknapsack(p, w, cap))
which(is$ksack == 1)
# [1] 1 2 3 4 6 , capacity 102 and total profit 280

## Example 2: multiple knapsack

p <- c(110, 150, 70, 80, 30, 5)
w <- c( 40, 60, 30, 40, 20, 5)
nelmin 17

k <- c(65, 85)

is <- mknapsack(p, w, k)
# kps 1: 2,6; kps 2: 1,4; value: 345; backtracks: 14

## Example 3: multiple knapsack

p <- c(78, 35, 89, 36, 94, 75, 74, 79, 80, 16)
w <- c(18, 9, 23, 20, 59, 61, 70, 75, 76, 30)
k <- c(103, 156)
is <- mknapsack(p, w, k)
# kps 1: 1,3,6; kps 2: 4,5,9; value: 452; backtracks: 4

## Example 4: subset sum

p <- seq(2, 44, by = 2)^2
w <- p
is <- mknapsack(p, w, 2012)
sum((2 * which(is$ksack == 1))^2)

## Example 5: maximize number of items

w <- seq(2, 44, by = 2)^2
p <- numeric(22) + 1
is <- mknapsack(p, w, 2012)

nelmin Nelder-Mead Minimization Method

Description
An implementation of the Nelder-Mead algorithm for derivative-free optimization / function mini-
mization.

Usage
nelmin(fn, x0, tol = 1e-10, ...,
maxfeval = 1e4, step = rep(1.0, length(x0)))

nelminb(fn, x0, lower, upper, tol = 1e-10, ...,

maxfeval = 10000, step = rep(1, length(x0)))

Arguments
fn nonlinear function to be minimized.
x0 starting point for the iteration.
tol terminating limit for the variance of function values; can be made *very* small,
like tol=1e-50.
maxfeval maximum number of function evaluations.
step size and shape of initial simplex; relative magnitudes of its elements should
reflect the units of the variables.
... additional arguments to be passed to the function.
lower, upper lower and upper bounds.
18 nelmin

Details
Also called a simplex method for finding the local minimum of a function of several variables.
The method is a pattern search that compares function values at the vertices of the simplex. The
process generates a sequence of simplices with ever reducing sizes.
The simplex function minimisation procedure due to Nelder and Mead (1965), as implemented by
ONeill (1971), with subsequent comments by Chambers and Ertel 1974, Benyon 1976, and Hill
1978. For another elaborate implementation of Nelder-Mead in R based on Matlab code by Kelley
see package dfoptim.
nelminb uses transfinite to define the function on all of R^n and to retransform the solution to
the bounded domain. The starting value is not allowed to lie on the boundary.

Value
List with following components:

xmin minimum solution found.

fmin value of f at minimum.
fcount number of iterations performed.
restarts number of restarts.
errmess error message

Note
Original FORTRAN77 version by R ONeill; MATLAB version by John Burkardt under LGPL
license. Re-implemented in R by Hans W. Borchers.

References
Nelder, J., and R. Mead (1965). A simplex method for function minimization. Computer Journal,
Volume 7, pp. 308-313.
ONeill, R. (1971). Algorithm AS 47: Function Minimization Using a Simplex Procedure. Applied
Statistics, Volume 20(3), pp. 338-345.

See Also
dfoptim::nmk

Examples
## Classical tests as in the article by Nelder and Mead
# Rosenbrock's parabolic valley
rpv <- function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2
x0 <- c(-2, 1)
nelmin(rpv, x0) # 1 1

# Fletcher and Powell's helic valley

fphv <- function(x)
100*(x[3] - 10*atan2(x[2], x[1])/(2*pi))^2 +
SIAM test functions 19

(sqrt(x[1]^2 + x[2]^2) - 1)^2 +x[3]^2

x0 <- c(-1, 0, 0)
nelmin(fphv, x0) # 1 0 0

# Powell's Singular Function (PSF)

psf <- function(x) (x[1] + 10*x[2])^2 + 5*(x[3] - x[4])^2 +
(x[2] - 2*x[3])^4 + 10*(x[1] - x[4])^4
x0 <- c(3, -1, 0, 1)
nelmin(psf, x0) # 0 0 0 0, needs maximum number of function calls

# Bounded version of Nelder-Mead

lower <- c(-Inf, 0, 0)
upper <- c( Inf, 0.5, 1)
x0 <- c(0, 0.1, 0.1)
nelminb(fnRosenbrock, c(0, 0.1, 0.1), lower, upper)
# $xmin = c(0.7085595, 0.5000000, 0.2500000)
# $fmin = 0.3353605

## Not run:
# Can run Rosenbrock's function in 30 dimensions in one and a half minutes:
nelmin(fnRosenbrock, rep(0, 30), tol=1e-20, maxfeval=10^7)
# $xmin
# [1] 0.9999998 1.0000004 1.0000000 1.0000001 1.0000000 1.0000001
# [7] 1.0000002 1.0000001 0.9999997 0.9999999 0.9999997 1.0000000
# [13] 0.9999999 0.9999994 0.9999998 0.9999999 0.9999999 0.9999999
# [19] 0.9999999 1.0000001 0.9999998 1.0000000 1.0000003 0.9999999
# [25] 1.0000000 0.9999996 0.9999995 0.9999990 0.9999973 0.9999947
# $fmin
# [1] 5.617352e-10
# $fcount
# [1] 1426085
# elapsed time is 96.008000 seconds
## End(Not run)

SIAM test functions Trefethen and Wagon Test Functions

Description
Test functions for global optimization posed for the SIAM 100-digit challenge in 2002 by Nick
Trefethen, Oxford University, UK.

Usage
fnTrefethen(p2)
fnWagon(p3)

Arguments
p2 Numerical vector of length 2.
p3 Numerical vector of length 3.
20 simpleDE

Details
These are highly nonlinear and oscillating functions in two and three dimensions with thousands of
local mimima inside the unit square resp. cube (i.e., [-1, 1] x [-1, 1] or [-1, 1] x [-1, 1] x [-1, 1]).

Value
Function value is a single real number.

Author(s)
HwB <[email protected]>

References
F. Bornemann, D. Laurie, S. Wagon, and J. Waldvogel (2004). The SIAM 100-Digit Challenge: A
Study in High-Accuracy Numerical Computing. Society for Industrial and Applied Mathematics.

Examples
x <- 2*runif(5) - 1
fnTrefethen(x)
fnWagon(x)

## Not run:
T <- matrix(NA, nrow=1001, ncol=1001)
for (i in 1:1001) {
for (j in 1:1001) {
T[i, j] <- fnTrefethen(c(x[i], y[j]))
}
}
image(x, y, T)
contour(x, y, T, add=TRUE)

## End(Not run)

simpleDE Simple Differential Evolution Algorithm

Description
Simple Differential Evolution for Minimization.

Usage
simpleDE(fun, lower, upper, N = 64, nmax = 256, r = 0.4,
confined = TRUE, log = FALSE)
simpleDE 21

Arguments
fun the objective function to be minimized.
lower vector of lower bounds for all coordinates.
upper vector of upper bounds for all coordinates.
N population size.
nmax bound on the number of generations.
r amplification factor.
confined logical; stay confined within bounds.
log logical; shall a trace be printed.

Details
Evolutionary search to minimize a function: For points in the current generation, children are
formed by taking a linear combination of parents, i.e., each member of the next generation has
the form
p1 + r(p2 p3 )
where the pi are members of the current generation and r is an amplification factor.

Value
List with the following components:

fmin function value at the minimum found.

xmin numeric vector representing the minimum.

Note
Original Mathematica version by Dirk Laurie in the SIAM textbook. Translated to R by Hans W
Borchers.

Author(s)
HwB <[email protected]>

References
Dirk Laurie. A Complex Optimization". Chapter 5 In: F. Bornemann, D. Laurie, S. Wagon, and J.
Waldvogel (Eds.). The SIAM 100-Digit Challenge. Society of Industrial and Applied Mathematics,
2004.

See Also
simpleEA, DEoptim in the DEoptim package.
22 simpleEA

Examples
simpleDE(fnTrefethen, lower = c(-1,-1), upper = c(1,1))
# $fmin
# [1] -3.306869
# $xmin
# [1] -0.02440308 0.21061243 # this is the true global optimum!

simpleEA Simple Evolutionary Algorithm

Description
Simple Evolutionary Algorithm for Minimization.

Usage
simpleEA(fn, lower, upper, N = 100, ..., con = 0.1, new = 0.05,
tol = 1e-10, eps = 1e-07, scl = 1/2, confined = FALSE, log = FALSE)

Arguments
fn the objective function to be minimized.
lower vector of lower bounds for all coordinates.
upper vector of upper bounds for all coordinates.
N number of children per parent.
... additional parameters to be passed to the function.
con percentage of individuals concentrating to the best parents.
new percentage of new individuals not focussed on existing parents.
tol tolerance; if in the last three loops no better individuals were found up to this
tolerance, stop.
eps grid size bound to be reached.
scl scaling factor for shrinking the grid.
confined logical; shall the set of individuals be strictly respect the boundary? Default:
FALSE.
log logical, should best solution found be printed per step.

Details
Evolutionary search to minimize a function: For each point in the current generation, n random
points are introduced and the n best results of each generation (and its parents) are used to form the
next generation.
The scale shrinks the generation of new points as the algorithm proceeds. It is possible for some
children to lie outside the given rectangle, and therefore the final result may lie outside the unit
rectangle well. (TO DO: Make this an option.)
subsetsum 23

Value
List with the following components:

par numeric vector representing the minimum found.

val function value at the minimum found.
fun.calls number of function calls made.
rel.scl last scaling factor indicating grid size in last step.
rel.tol relative tolerance within the last three minima found.

Note
Original Mathematica Version by Stan Wagon in the SIAM textbook. Translated to R by Hans W
Borchers.

Author(s)
HwB <[email protected]>

References
Stan Wagon. Think Globally, Act Locally". Chapter 4 In: F. Bornemann, D. Laurie, S. Wagon,
and J. Waldvogel (Eds.). The SIAM 100-Digit Challenge. Society of Industrial and Applied Math-
ematics, 2004.

See Also
DEoptim in the DEoptim package.

Examples
simpleEA(fnTrefethen, lower=c(-1,-1), upper=c(1,1), log=FALSE)
# $par
# [1] -0.02440310 0.21061243 # this is the true global optimum!
# $val
# [1] -3.306869

subsetsum Subset Sum Problem

Description
Subset sum routine for positive integers.

Usage
subsetsum(S, t, method = "greedy")
24 subsetsum

Arguments
S vector of positive integers.
t target value.
method can be greedy or dynamic, where dynamic stands for the dynamic pro-
gramming approach.

Details
Searching for a set of elements in S that sum up to t by continuously adding more elements of S.
The first components will be preferred, i.e., if S is decreasing, the sum with larger elements will be
found, if increasing, the sum with smaller elements.
The dynamic method may be faster for large sets, but will also require much more memory if the
target value is large.

Value
List with the target value, if reached, and vector of indices of elements in S that sum up to t.
If no solution is found, the dynamic method will return indices for the largest value below the target,
the greedy method witll return NULL.

Note
Will be replaced by a compiled version.

Author(s)
HwB email: <[email protected]>

References
Horowitz, E., and S. Sahni (1978). Fundamentals of Computer Algorithms. Computer Science
Press, Rockville, ML.

See Also
maxsub

Examples
## Not run:
amount <- 4748652
products <-
c(30500,30500,30500,30500,42000,42000,42000,42000,
42000,42000,42000,42000,42000,42000,71040,90900,
76950,35100,71190,53730,456000,70740,70740,533600,
83800,59500,27465,28000,28000,28000,28000,28000,
26140,49600,77000,123289,27000,27000,27000,27000,
27000,27000,80000,33000,33000,55000,77382,48048,
51186,40000,35000,21716,63051,15025,15025,15025,
Testfunctions 25

15025,800000,1110000,59700,25908,829350,1198000,1031655)

# prepare set
prods <- products[products <= amount] # no elements > amount
prods <- sort(prods, decreasing=TRUE) # decreasing order

# now find one solution

system.time(is <- subsetsum(prods, amount))
# user system elapsed
# 0.320 0.032 0.359

prods[is]
# [1] 70740 70740 71190 76950 77382 80000 83800
# [8] 90900 456000 533600 829350 1110000 1198000

sum(prods[is]) == amount
# [1] TRUE
## End(Not run)

Testfunctions Optimization Test Functions

Description

Simple and often used test function defined in higher dimensions and with analytical gradients,
especially suited for performance tests. Analytical gradients, where existing, are provided with the
gr prefix. The dimension is determined by the length of the input vector.

Usage

fnRosenbrock(x)
grRosenbrock(x)
fnRastrigin(x)
grRastrigin(x)
fnNesterov(x)
grNesterov(x)
fnNesterov1(x)
fnHald(x)
grHald(x)
fnShor(x)
grShor(x)

Arguments

x numeric vector of a certain length.

26 Testfunctions

Details
Rosenbrock Rosenbrocks famous valley function from 1960. It can also be regarded as a least-
squares problem:
n1
X
(1 xi )2 + 100(xi+1 x2i )2
i=1

No. of Vars.: n >= 2

Bounds: -5.12 <= xi <= 5.12
Local minima: at f(-1, 1, . . . , 1) for n >= 4
Minimum: 0.0
Solution: xi = 1, i = 1:n

Nesterov Nesterovs smooth adaptation of Rosenbrock, based on the idea of Chebyshev polyno-
mials. This function is even more difficult to optimize than Rosenbrocks:
n1
X
(x1 1)2 /4 + (1 + xi+1 2x2i )
i=1

No. of Vars.: n >= 2

Bounds: -5.12 <= xi <= 5.12
Local minima: ?
Minimum: 0.0
Solution: xi = 1, i = 1:n

Rastrigin Rastrigins function is a famous, non-convex example from 1989 for global optimiza-
tion. It is a typical example of a multimodal function with many local minima:
n
X
10n + (x2i 10 cos(2xi ))
1

No. of Vars.: n >= 2

Bounds: -5.12 <= xi <= 5.12
Local minima: many
Minimum: 0.0
Solution: xi = 0, i = 1:n

Hald Halds function is a typical example of a non-smooth test function, from Hald and Madsen
in 1981.
x1 + x2 t i
max
1in 1 + x3 ti + x4 t2 3 exp(ti )
i + x5 t i

where ti = 1 + (i 1)/10 for 1 i 21.

transfinite 27

No. of Vars.: n =5
Bounds: -1 <= xi <= 1
Local minima: ?
Minimum: 0.0001223713
Solution: (0.99987763, 0.25358844, -0.74660757, 0.24520150, -0.03749029)

Shor Shors function is another typical example of a non-smooth test function, a benchmark for
Shors R-algorithm.

Value

Returns the values of the test function resp. its gradient at that point. If an analytical gradient is not
available, a function computing the gradient numerically will be provided.

References

Search the Internet.

Examples
x <- runif(5)
fnHald(x); grHald(x)

# Compare analytical and numerical gradient

shor_gr <- function(x) adagio:::ns.grad(fnShor, x) # internal gradient
grShor(x); shor_gr(x)

transfinite Boxed Region Transformation

Description

Transformation of a box/bound constrained region to an unconstrained one.

Usage

transfinite(lower, upper, n = length(lower))

Arguments

lower, upper lower and upper box/bound constraints.

n length of upper, lower if both are scalars, to which they get repeated.
28 transfinite

Details
Transforms a constraint region in n-dimensional space bijectively to the unconstrained Rn space,
applying a atanh resp. exp transformation to each single variable that is bound constraint.
It provides two functions, h: B = []x...x[] --> R^n and its inverse hinv. These functions can,
for example, be used to add box/bound constraints to a constrained optimization problem that is to
be solved with a (nonlinear) solver not allowing constraints.

Value
Returns to functions as components h and hinv of a list.

Note
Based on an idea of Ravi Varadhan, intrinsically used in his implementation of Nelder-Mead in the
dfoptim package.
For positivity constraints, x>=0, this approach is considered to be numerically more stable than
x-->exp(x) or x-->x^2.

Examples
lower <- c(-Inf, 0, 0)
upper <- c( Inf, 0.5, 1)
Tf <- transfinite(lower, upper)
h <- Tf$h; hinv <- Tf$hinv

## Not run:
## Solve Rosenbrock with one variable restricted
rosen <- function(x) {
n <- length(x)
x1 <- x[2:n]; x2 <- x[1:(n-1)]
sum(100*(x1-x2^2)^2 + (1-x2)^2)
}
f <- function(x) rosen(hinv(x)) # f must be defined on all of R^n
x0 <- c(0.1, 0.1, 0.1) # starting point not on the boundary!
nm <- nelder_mead(h(x0), f) # unconstraint Nelder-Mead
hinv(nm$xmin); nm$fmin # box/bound constraint solution
# [1] 0.7085596 0.5000000 0.2500004
# [1] 0.3353605

## End(Not run)
Index

Topic graphs grNesterov (Testfunctions), 25

hamiltonian, 7 grRastrigin (Testfunctions), 25
Topic manip grRosenbrock (Testfunctions), 25
transfinite, 27 grShor (Testfunctions), 25
Topic optimize
CMAES, 4 hamiltonian, 7
knapsack, 10
maxsub, 13 knapsack, 10
nelmin, 17
maxempty, 11
simpleDE, 20
maxquad, 12
simpleEA, 22
maxsub, 13, 24
subsetsum, 23
maxsub2d (maxsub), 13
Topic package
mknapsack, 15
adagio-package, 2
Topic testfunctions nelmin, 17
maxquad, 12 nelminb (nelmin), 17
SIAM test functions, 19
Testfunctions, 25 pureCMAES (CMAES), 4
Topic visualize
fminviz,flineviz, 6 SIAM test functions, 19
simpleDE, 20
adagio (adagio-package), 2 simpleEA, 21, 22
adagio-package, 2 subsetsum, 23
assignment, 3
Testfunctions, 25
CMAES, 4 transfinite, 27

flineviz (fminviz,flineviz), 6
fminviz (fminviz,flineviz), 6
fminviz,flineviz, 6
fnHald (Testfunctions), 25
fnNesterov (Testfunctions), 25
fnNesterov1 (Testfunctions), 25
fnRastrigin (Testfunctions), 25
fnRosenbrock (Testfunctions), 25
fnShor (Testfunctions), 25
fnTrefethen (SIAM test functions), 19
fnWagon (SIAM test functions), 19

grHald (Testfunctions), 25

29