Chapter 7
Scientists and engineers often have to solve for the maximum or minimum of a multi-
dimensional function, sometimes with constraints on the values of some or all of the
variables. This is known as optimization, and is a rich, highly developed, and often
difficult problem. Generally the problem is phrased as a minimization, which shows
its kinship to the least-squares data fitting procedures discussed in a subsequent chap-
ter. If a maximum is desired, one simply solves for the minimum of the negative of
the function. The greatest difficulties typically arise if the multi-dimensional surface
has local minima in addition to the global minimum, because there is no way to show
that the minimum is local except by trial and error.
R has three functions in the base installation: optimize() for one-dimensional
problems, optim() for multi-dimensional problems, and constrOptim() for opti-
mization with linear constraints. (Note that even “unconstrained” optimization nor-
mally is constrained by the limits on the search range.) We shall consider each of
these with suitable examples, and introduce several add-on packages that expand the
power of the basic functions. Optimization is a sufficiently large and important topic
to deserve its own task view in R, at
In addition to the packages considered in this chapter, the interested reader
should become acquainted with the nloptr package, which is considered one of
the strongest and most comprehensive optimization packages in R. According to its
synopsis,“nloptr is an R interface to NLopt. NLopt is a free/open-source library for
nonlinear optimization, providing a common interface for a number of different free
optimization routines available online as well as original implementations of various
other algorithms.
7.1 One-dimensional optimization
The base R function for finding minima (the default) or maxima of functions of a
single variable is optimize(). As a concrete example, given a 16 ×20 sheet of
cardboard, find the size x of the squares to be cut from the corners that maximizes
the volume of the open box formed when the sides are folded up. x must be in the
range (0,8).
> optimize(function(x) x*(20-2*x)*(16-2*x), c(0,8), maximum=T)
0.0 1.0 2.0 3.0
-2 -1 0 1 2
Figure 7.1: Plot of function f (x) = x sin(4x) showing several maxima and minima.
[1] 2.944935
[1] 420.1104
Consider next the use of optimize with the function
> f = function(x) x*sin(4*x)
which plotted looks like this (Figure 7.1):
> curve(f,0,3)
It has two minima in the x = 0 3 range, with the global minimum near 2.8,
and two maxima, with the global maximum near 2.0. Applying optimize() in the
simplest way yields
> optimize(f,c(0,3))
[1] 1.228297
[1] -1.203617
which gives the local minimum because it is the first minimum encountered by the
search algorithm (Brent’s method, which combines root bracketing, bisection, and
inverse quadratic interpolation). Because we have a plot of the function, we can see
that we must exclude the local minimum from the lower and upper endpoints of the
search interval.
> optimize(f,c(1.5,3))
[1] 2.771403
[1] -2.760177
To find the global maximum we enter
> optimize(f,c(1,3),maximum=TRUE)
[1] 1.994684
[1] 1.979182
We could have obtained the same result by minimizing the negative of the function
> optimize(function(x) -f(x),c(1,3))
[1] 1.994684
[1] -1.979182
which finds the maximum in the right place but, of course, yields the negative of the
function value at the maximum.
If necessary, the desired accuracy can be adjusted with the tol option in the
function call.
The pracma package contains the function findmins(), which finds the posi-
tions of all the minima in the search interval by dividing it n times (default n = 100)
and applying optimize in each interval. To find the values at those minima, evaluate
the function.
> require(pracma)
> f.mins = findmins(f,0,3)
> f.mins # x values at the minima
[1] 1.228312 2.771382
> f(f.mins[1:2]) # function evaluated at the minima
[1] -1.203617 -2.760177
The Examples section of the help page for optimize() shows how to include
a parameter in the function call. It also shows how, for a function with a very flat
minimum, the wrong solution can be obtained if the search interval is not properly
chosen. Unfortunately, there is no clear way to choose the search interval in such
cases, so if the results are not as expected from inspecting the graph of the function,
other intervals should be explored.
The function f (x) = |x
8| yields some interesting behavior (Figure 7.2).
> f = function(x) abs(x^2-8)
> curve(f,-4,4)
Straightforward solving for the maximum over the entire range yields the result at
x = 0.
> optimize(f,c(-4,4),maximum=T)
[1] -1.110223e-16
[1] 8
Excluding the middle from the search interval finds the maxima at the extremes.

Get Using R for Numerical Analysis in Science and Engineering now with O’Reilly online learning.

O’Reilly members experience live online training, plus books, videos, and digital content from 200+ publishers.