\section{Statistical search for optimization}
\label{sec:optim}

<<echo=false,results=hide>>=
library(tgp)
seed <- 0; set.seed(seed)
@ 

There has been considerable recent interest in the use of
statistically generated search patterns (i.e., locations of relatively
likely optima) for optimization.  A popular approach is to estimate a
statistical (surrogate) model, and use it to design a set of
well-chosen candidates for further evaluation by a direct optimization
routine.  Such statistically designed search patterns can be used
either to direct the optimization completely (e.g.,
\cite{JoneSchoWelc1998} or \cite{RommShoe2007}) or to work in hybrid
with local pattern search optimization (as in
\cite{TaddLeeGrayGrif2009}).  An bonus feature of the statistical
surrogate approach is that it may be used to tackle problems of
optimization under uncertainty, wherein the function being optimized
is observed with noise.  In this case the search is for input
configurations which optimize the response with high probability.
Direct-search methods would not apply in this scenario without
modification. However, a sensible hybrid could involve inverting the
relationship between the two approaches so that direct-search is used
on deterministic predictive surfaces from the statistical surrogate
model.  This search can be used to find promising candidates to
compliment space-filling ones at which some statistical
improvement criterion is evaluated.

Towards situating {\tt tgp} as a promising statistical surrogate model
for optimization (in both contexts) the approach developed by Taddy,
et al.~\cite{TaddLeeGrayGrif2009}, has been implemented to produce a
list of input locations that is ordered by a measure of the potential
for new optima.  The procedure uses samples from the posterior
predictive distribution of treed GP regression models to estimate
improvement statistics and build an ordered list of search locations
which maximize expected improvement.  The single location improvement
is defined $I(\mb{x}) = \mathrm{max}\{f_{min}-f(\mb{x}),0\}$, where
$f_{min}$ is the minimum evaluated response in the search (refer to
\cite{SchoWelcJone1998} for extensive discussion on general
improvement statistics and initial vignette~\cite{gramacy:2007} for
details of a base implementation in {\tt tgp}).  Thus, a high
improvement corresponds to an input location that is expected to be
much lower than the current minimum.  The criterion is easily changed
to a search for maximum values through negation of the response.  The
improvement is always non-negative, as points which do not turn out to
be new minimum points still provide valuable information about the
output surface.  Thus, in the expectation, candidate locations will be
rewarded for high response uncertainty (indicating a poorly explored
region of the input space), as well as for low mean predicted
response.  Our {\tt tgp} generated search pattern will consist of $m$
locations that recursively maximize (over a discrete candidate set) a
sequential version of the expected multi-location improvement
developed by Schonlau, et al.~\cite{SchoWelcJone1998}, defined as
$\mbb{E}\left[I^g(\mb{x}_1, \ldots, \mb{x}_m)\right]$ where
\begin{equation} 
\label{eqn:imult} I^g(\mb{x}_1,
\ldots, \mb{x}_m) = \left(\mathrm{max}\{(f_{min}-f(\mb{x}_1)), \ldots,
(f_{min}-f(\mb{x}_m)), 0 \}\right)^g.  
\end{equation}
Increasing $g \in \{0,1,2,3,\ldots\}$ increases the global scope of
the criteria by rewarding in the expectation extra variability at
$\mb{x}$.  For example, $g=0$ leads to $\mbb{E}[I^0(\mb{x})] =
\Pr(I(\mb{x})>0)$ (assuming the convention $0^0=0$), $g=1$ yields the
standard statistic, and $g=2$ explicitly rewards the improvement
variance since $\mbb{E}[I^2(\mb{x})] = \mr{var}[I(\mb{x})] +
\mbb{E}[I(\mb{x})]^2$.  For further discussion on the role of $g$, see
\cite{SchoWelcJone1998} .

Finding the maximum expectation of (\ref{eqn:imult}) is practically
impossible for the full posterior distribution of $I^g(\mb{x}_1,
\ldots, \mb{x}_m)$, and would require conditioning on a single fit for
the model parameters (for example, static imputation of predictive GP
means can be used to recursively build the improvement set
\cite{GinsLe-RCarr2009}).  However, {\tt tgp} just seeks to maximize
over a discrete list of predictive locations.  In fact, the default is
to return an ordering for the entire {\tt XX} matrix, thus defining a
ranking of predictive locations by order of decreasing expected
improvement.  There is no restriction on the form for {\tt
  XX}.\footnote{A full optimization routine would require that the
  search pattern is placed within an algorithm iterating towards
  convergence, as in \cite{TaddLeeGrayGrif2009}.  However, we
  concentrate here on the statistical problem of choosing the next
  samples optimally.  We shall touch on issues of convergence in
  Section \ref{sec:optimskel} where we describe a skeleton scheme for
  optimization extending {\sf R}'s internal {\tt optim}
  functionality.}  The structure of this scheme will dictate the form
for {\tt XX}.  If it is the case that we seek simply to explore the
input space and map a list of potential locations for improvement,
using LHS to choose {\tt XX} will suffice.

The discretization of decision space allows for a fast iterative
solution to the optimization of $\mbb{E}\left[I^g(\mb{x}_1, \ldots,
  \mb{x}_m)\right]$.  This begins with evaluation of the simple
improvement $I^g(\tilde{\mb{x}}_i)$ over $\tilde{\mb{x}}_i \in {\bf
  \tilde X}$ at each of $T=$ {\tt BTE[2]-BTE[1]} MCMC iterations (each
corresponding to a single posterior realization of {\tt tgp}
parameters and predicted response after burn-in) to obtain the
posterior sample
\begin{equation*}
 \mathcal{I} = \left\{ 
\begin{array}{rcl} I^g(
\tilde{\mb{x}}_1)_1& \ldots& I^g(\tilde{\mb{x}}_m)_1\\ 
&\vdots& \\ 
I^g( \tilde{\mb{x}}_1)_T& \ldots& I^g(\tilde{\mb{x}}_m)_T
 \end{array}\right\}.  
\end{equation*} 
Recall that in {\tt tgp} parlance, and as input to the {\tt b*}
functions: $\tilde{\mb{X}}\equiv $ {\tt XX}.

We then proceed iteratively to build an {\it ordered} collection of
$m$ locations according to an iteratively refined improvement:
Designate $\mb{x}_1 = \mathrm{argmax}_{\tilde{\mb{x}} \in {\bf \tilde
    X}} \mbb{E}\left[I^g( \tilde{\mb{x}})\right]$, and for
$j=2,\ldots,m$, given that $\mb{x}_1, \ldots, \mb{x}_{j-1}$ are
already included in the collection, the next member is
\begin{eqnarray*}
 \mb{x}_j &=& \mathrm{argmax}_{\tilde{\mb{x}} \in {\bf \tilde X}} \mbb{E}\left[
\mathrm{max}\{I^g( \mb{x}_1, \ldots, \mb{x}_{j-1}), I^g(\tilde{\mb{x}}) \} \right]\\ 
&=& \mathrm{argmax}_{\tilde{\mb{x}} \in {\bf \tilde X}}
\mbb{E}[\left(\mathrm{max}\{(f_{min}-f(\mb{x}_1)), \ldots,
(f_{min}-f(\mb{x}_{j-1})), (f_{min}-f(\tilde{\mb{x}})), 0\}\right)^g ] \\
&=& \mathrm{argmax}_{\tilde{\mb{x}} \in {\bf \tilde X}} \mbb{E}\left[I^g(\mb{x}_1,
\ldots, \mb{x}_{j-1},\tilde{\mb{x}})\right].
\end{eqnarray*} 
Thus, after each $j^{\rm th}$ additional point is added to the set, we have
the maximum expected $j$--location improvement conditional on the first
$j-1$ locations.  This is not necessarily the unconditionally maximal
expected $j$--location improvement; instead, point $\mb{x}_j$ is the
location which will cause the greatest increase in expected
improvement over the given $(j-1)$--location expected improvement.

The posterior sample $\mathcal{I}$ acts as a discrete approximation to
the true posterior distribution for improvement at locations within
the candidate set {\tt XX}.  Based upon this approximation, iterative
selection of the point set is possible without any re-fitting of the
{\tt tgp} model.  Conditional on the inclusion of
$\tilde{\mb{x}}_{i_1},\ldots,\tilde{\mb{x}}_{i_{l-1}}$ in the
collection, a posterior sample of the $l$--location improvement
statistics is calculated as
\begin{equation*}
\mathcal{I}_l = \left\{
\begin{array}{rcl}
I^g( \tilde{\mb{x}}_{i_1},\ldots,\tilde{\mb{x}}_{i_{l-1}}, \tilde{\mb{x}}_1)_1 
& \ldots& 
I^g(  \tilde{\mb{x}}_{i_1},\ldots,\tilde{\mb{x}}_{i_{l-1}}, \tilde{\mb{x}}_m)_1\\
&\vdots& \\
I^g(\tilde{\mb{x}}_{i_1},\ldots,\tilde{\mb{x}}_{i_{l-1}}, {\tilde x}_1)_T& 
\ldots& I^g(\tilde{\mb{x}}_{i_1},\ldots,\tilde{\mb{x}}_{i_{l-1}},\tilde{\mb{x}}_m)_T
\end{array}\right\},
\end{equation*}
where the element in the $t^{\rm th}$ row and $j^{\rm th}$ column of
this matrix is calculated as max$\{I^g(\tilde{\mb{x}}_{i_1}$, $\ldots,$
$\tilde{\mb{x}}_{i_{l-1}})_t$, $I^g(\tilde{\mb{x}}_j)_t\}$ and the
$l^{\rm th}$ location included in the collection corresponds to the
column of this matrix with maximum average.  Since the multi-location
improvement is always at least as high as the improvement at any
subset of those locations, the same points will not be chosen twice
for inclusion.  In practice, very few iterations (about 10\% of the
total candidate size under the default inference and regression
model(s)) through this ordering process can be performed before the
iteratively updated improvement statistics become essentially zero.
Increasing the number of MCMC iterations ({\tt BTE[2]-BTE[1]}) can
mitigate this to a large extent.\footnote{Once a zero (maximal)
  iterative improvement is attained the rest of the ranking is
  essentially arbitrary, at which point {\tt tgp} cuts off the process
  prematurely.}  We refer the reader to \cite{TaddLeeGrayGrif2009} for
further details on this approach to multi-location improvement search.



\subsection{A simple example}

We shall use the Rosenbrock function to illustrate the production of
an ordered collection of (possible) adaptive samples to maximize the
expected improvement within {\tt tgp}.  Specifically, the two
dimensional Rosenbrock function is defined as
<<>>=
rosenbrock <- function(x){ 
  x <- matrix(x, ncol=2)
  100*(x[,1]^2 - x[,2])^2 + (x[,1] - 1)^2 
}
@ 
and we shall bound the search space for adaptive samples to the
rectangle: $-1\le x_i \le 5$ for $i=1,2$.  The single global minimum
of the Rosenbrock function is at $(1,1)$.
<<>>=
rosenbrock(c(1,1))
@ 
This function involves a long steep valley with a gradually sloping
floor, and is considered to be a difficult problem for local
optimization routines.

We begin by drawing an LHS of 40 input locations within the 
bounding rectangle, and evaluating the function at these locations.
<<>>=
rect <- cbind(c(-1,-1),c(5,5))
X <- lhs(40, rect)
Z <- rosenbrock(X)
@ 

We will fit a {\tt bgp} model to this data to predict the Rosenbrock
response at unobserved (candidate) input locations in {\tt XX}.  The
{\tt improv} argument may be used to obtain an ordered list of places
where we should be looking for new minima.  In particular, specifying
{\tt improv=c(1,10)} will return the 10 locations which maximize the
iterative multi-location expected improvement function, with $g=1$
(i.e., Eq.~(\ref{eqn:imult})).  Note that {\tt improv=TRUE} is also
possible, in which case {\tt g} defaults to one and the entire list of
locations is ranked.  Our candidate set is just a space filling LHS
design.  In other situations, it may be useful to build an informative
LHS design (i.e., to specify {\tt shape} and {\tt mode} arguments for
the {\tt lhs} function) to reflect what is already known about the
location of optima.
<<>>=
XX <- lhs(200, rect)
rfit <- bgp(X,Z,XX,improv=c(1,10), verb=0)
@ 
Upon return, the \verb!"tgp"!-class object {\tt rfit} includes the matrix
{\tt improv}, which is a list of the expected single location
improvement for the 200 {\tt XX} locations, and the top 10 ranks.
Note that the {\tt rank}s for those points which are not included in
the top 10 are set to {\tt nrow(XX)=}\Sexpr{nrow(XX)}.  Here are the
top 10:
<<>>=
cbind(rfit$improv,XX)[rfit$improv$rank <= 10,]
@ 
This iterative algorithm may produce ranks that
differ significantly from a straightforward ordering of expected
improvement.  This leads to a list that better explores the input
space, since the expected improvement is naturally balanced against a
desire to search the domain.

We plot the results with the usual function, by setting {\tt
  as="improv"}, in Figure \ref{optim-fit1}.
\begin{figure}[htb!]
<<label=optim-fit1,fig=TRUE,echo=TRUE,width=10,height=6,include=FALSE>>= 
plot(rfit, as="improv")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=6.5in,trim=0 25 0 25]{tgp2-optim-fit1}
\caption{The {\em left} panel shows the mean predicted Rosenbrock function
  response, and on the {\em right} we have expected single location
  improvement with the top 10 points (labelled by rank) plotted on
  top.}
\label{optim-fit1}
\end{figure}
The banana--shaped region of higher expected improvement corresponds to
the true valley floor for the Rosenbrock function, indicating the that
{\tt bgp} model is doing a good job of prediction.  Also, we note that
the ordered input points are well dispersed throughout the valley---a
very desirable property for adaptive sampling candidates.

It is straightforward, with {\tt predict.tgp}, to obtain a new ordering
for the more global {\tt g=5} (or any new {\tt g}).
Figure \ref{optim-fit2} shows a more diffuse expected improvement
surface and a substantially different point ordering.  In practice, we
have found that {\tt g=2} provides a good compromise between local and
global search. 

\begin{figure}[htb!]
<<label=optim-fit2,fig=TRUE,echo=TRUE,width=5,height=6,include=FALSE>>= 
rfit2 <- predict(rfit, XX=XX, BTE=c(1,1000,1), improv=c(5,20), verb=0) 
plot(rfit2, layout="as", as="improv")
@
<<echo=false,results=hide>>=
graphics.off()
@
\centering
\includegraphics[width=3.25in,trim=0 25 0 25]{tgp2-optim-fit2}
\caption{The expected improvement surface and top 20 ordered locations,
for {\tt g=5}.}
\label{optim-fit2}
\end{figure}


\subsection{A skeleton optimization scheme}
\label{sec:optimskel}

%% The nature of global optimization demands that a fine balance be
%% struck between global and local search.  Therefore, designing a
%% one--size--fits--all approach would be a daunting task.  For one
%% thing, assessing convergence in any formal sense would be quite
%% difficult, although in practice it would be straightforward to
%% ``force'' convergence by (eventually) focusing the method on finding a
%% local solution.  In the case where the function evaluations are
%% deterministic, final convergence to a the local solution is always
%% possible through the use of {\tt R}'s {\tt optim} function, for
%% example.  A method using {\tt tgp} based on a similar, but more
%% formalized approach, using a direct/pattern search (in place of {\tt
%%   optim}) has been recently demonstrated in the context of
%% sequentially designing computer experiments to solve an optimization
%% problem \cite{TaddLeeGrayGrif2009}. Generally speaking, the result is
%% a sensible compromise between local and global search.  When the
%% function evaluations are noisy one can always create a deterministic
%% approximation, i.e., via the MAP predictive distribution (i.e., a
%% kriging surrogate), for use with {\tt optim} in order to obtain
%% convergence to a local optima.
%% 
%% It may be possible to base assessments of convergence on the
%% improvement statistic, which would naturally tend to zero as more
%% points are added into the design. But any such assessments would hinge
%% upon being able to drive the (Monte Carlo) method used to infer the
%% model parameters---on which the improvement statistic is based---to a
%% fixed point.  In the context of MCMC this is only guaranteed as the
%% number of samples gathered tends to infinity.  Even if obtaining
%% asymptotic convergence in this way is clearly a pipe dream, the
%% practical application of this idea, and those based on local
%% optimization mentioned above, can still bear fruit.  Insight into
%% convergence in practice is still a very tangible concept.  Moreover,
%% for many applications the considerations of convergence may even take
%% a back seat to other budgetary constraints where the efficient
%% allocation of an available resource (say computer cycles) is more
%% important than a bottom--line based upon convergence which may only be
%% achieved at all costs in the best of scenarios.

The capabilities outlined above are useful in their own right, as a
search list or candidate set ranked by expected improvement gain
provides concrete information about potential optima.  However, a full
optimization framework requires that the production of these sets of
search locations are nested within an iterative search scheme.  The
approach taken by Taddy, et al.~\cite{TaddLeeGrayGrif2009}, achieves
this by taking the {\tt tgp} generated sets of locations and using
them to augment a local optimization search algorithm.  In this way,
the authors are able to achieve robust solutions which balance the
convergence properties of the local methods with the global scope
provided by {\tt tgp}.  Indeed, any optimization routine capable of
evaluating points provided by an outside source could benefit from a
{\tt tgp} generated list of search locations.

In the absence of this sort of formal hybrid search algorithm, it is
still possible to devise robust optimization algorithms based around
{\tt tgp}.  A basic algorithm is as follows: first, use a LHS to
explore the input space (see the {\tt lhs} function included in {\tt
  tgp}).  Repeatedly fit one of the {\tt b*} models with {\tt
  improv!=FALSE} to the evaluated iterates to produce a search set,
then evaluate the objective function over this search set, as
described earlier.  Then evaluate the objective function over the
highest ranked locations in the search set.  Continue until you are
confident that the search has narrowed to a neighborhood around the
true optimum (a good indicator of this is when all of the top-ranked
points are in the same area).  At this point, the optimization may be
completed by {\tt optim}, {\sf R}'s general purpose local optimization
algorithm in order to guarentee convergence.  The {\tt optim} routine
may be initialized to the best input location (i.e. corresponding the
most optimal function evaluation) found thus far by {\tt tgp}.

Note that this approach is actually an extreme version of a template
proposed by Taddy, et al.~\cite{TaddLeeGrayGrif2009}, where the
influence of global (i.e. {\tt tgp}) search is downweighted over time
rather than cut off.  In either case, a drawback to such approaches
is that they do not apply when the function being optimized is
deterministic.  An alternative scheme is to employ both {\tt tgp}
search and a local optimization at each iteration.  The idea is that a
mix of local and global information is provided throughout the entire
optimization, but with an added twist.  Rather than apply {\tt optim}
on the stochastic function directly, which would not converge due to
the noise, it can be applied on a deterministic (MAP) kriging surface
provided by {\tt tgp}.  The local optima obtained can be used to
augment the candidate set of locations where the improvement statistic
is gathered---which would otherwise be simple LHS.  That way the search
pattern produced on output is likely to have a candidate with high
improvement.  

To fix ideas, and for the sake of demonstration, the {\tt tgp} package
includes a skeleton function for performing a single iteration in the
derivative--free optimization of noisy black--box functions.  The
function is called {\tt optim.step.tgp}, and the name is intended to
emphasize that it performs a single step in an optimization by trading
off local {\tt optim}--based search of {\tt tgp} predictive (kriging
surrogate) surfaces, with the expected posterior improvement.  In
other words, it is loosely based on some the techniques alluded to
above, but is designed to be augmented/adjusted as needed.  Given $N$
pairs of inputs and responses $(\mb{X}, \mb{Z})$, {\tt optim.step.tgp}
suggests new points at which the function being optimized should be
evaluated.  It also returns information that can be used to assess
convergence.  An outline follows.

The {\tt optim.step.tgp} function begins by constructing a set of
candidate locations, either as a space filling LHS over the input
space (the default) or from a treed $D$--optimal design, based on a
previously obtained \verb!"tgp"!-class model.  {\sf R}'s {\tt optim}
command is used on the MAP predictive surface contained within the
object to obtain an estimate of the current best guess $\mb{x}$-location of
the optimal solution.  A standalone subroutine called {\tt
  optim.ptgpf} is provided for this specific task, to be used within
{\tt optim.step.tgp} or otherwise.  Within {\tt optim.step.tgp}, {\tt
  optim.ptgpf} is initialized with the data location currently
predicted to be the best guess of the minimum.  The optimal
$x$-location found is then added into the set of candidates as it is
likely that the expected improvement would be high there.

Then, a new \verb!"tgp"!-class object is obtained by applying a {\tt
  b*} function to $(\mb{X}, \mb{Z})$ whilst sampling from the
posterior distribution of the improvement statistic.  The best one,
two, or several locations with highest improvement ranks are suggested
for addition into the design.  The values of the maximum improvement
statistic are also returned in order to track progress in future
iterations.  The \verb!"tgp"!-class object returned is used to
construct candidates and initialize the {\tt optim.ptgpf} function in
future rounds.

To illustrate, consider the 2-d exponential data from the
initial vignette \cite{gramacy:2007} as our noisy
function $f$.
<<>>=
f <- function(x) { exp2d.Z(x)$Z }
@ 
Recall that this data is characterized by a mean value of
\[
f(\mb{x}) =  x_1 \exp(-x_1^2 - x_2^2)
\]
which is observed with a small amount of Gaussian noise (with sd
$=0.001$).  Elementary calculus gives that the minimum of $f$ is
obtained at $\mb{x} = (-\sqrt{1/2},0)$.

The {\tt optim.step.tgp} function requires that the search domain be
defined by a bounding rectangle, and we require an initial design to
start things off.  Here we shall use $[-2,6]^2$ with an LHS design
therein.
<<>>=
rect <- rbind(c(-2,6), c(-2,6))
X <- lhs(20, rect)
Z <- f(X)
@ 
The following code proceeds with several rounds of sequential 
design towards finding the minimum of {\tt f}.
<<keep.source=TRUE>>=
out <- progress <- NULL
for(i in 1:20) {
  
  ## get recommendations for the next point to sample
  out <- optim.step.tgp(f, X=X, Z=Z, rect=rect, prev=out, verb=0)

  ## add in the inputs, and newly sampled outputs
  X <- rbind(X, out$X)
  Z <- c(Z, f(out$X))
  
  ## keep track of progress and best optimum
  progress <- rbind(progress, out$progress)
}
@ 
The {\tt progress} can be tracked through the rows of a {\tt
  data.frame}, as constructed above, containing a listing of the input
location of the current best guess of the minimum for each round,
together with the value of the objective at that point, as well as the
maximum of the improvement statistic.
\begin{figure}[ht!]
\centering
<<label=optim-progress,fig=TRUE,echo=TRUE,width=14,height=7.5,include=FALSE>>=
par(mfrow=c(1,2))
matplot(progress[,1:2], main="x progress",
        xlab="rounds", ylab="x[,1:2]", type="l", lwd=2)
legend("topright", c("x1", "x2"), lwd=2, col=1:2, lty=1:2)
plot(log(progress$improv), type="l", main="max log improv",
     xlab="rounds", ylab="max log(improv)")
@
<<echo=false,results=hide>>=
graphics.off()
@
\includegraphics[trim=40 20 0 0]{tgp2-optim-progress}
%\vspace{-0.5cm}
\caption{Progress in iterations of {\tt optim.step.tgp} shown by tracking
  the $\mb{x}$--locations of the best guess of the minimum ({\em left})
  and the logarithm of the maximum of the improvement statistics at
  the candidate locations ({\em right})}
\label{f:optim:progress}
\end{figure}
In addition to printing this data to the screen, plots such as the
ones in Figure \ref{f:optim:progress} can be valuable for assessing
convergence.  As can be seen in the figure, the final iteration
gives an $\mb{x}$-value that is very close to the correct result,
and is (in some loose sense) close to convergence.
<<>>=
out$progress[1:2]
@ 

As mentioned above, if it is known that the function evaluations are
deterministic then, at any time, {\sf R}'s {\tt optim} routine can be
invoked---perhaps initialized by the $\bm{x}$-location in
\verb!out$progress!---and convergence to a local optimum thus
guaranteed.  Otherwise, the quantities in \verb!out$progress! will
converge, in some sense, as long as the number of MCMC rounds used in
each round, above, ($T=$ {\tt BTE[2]-BTE[1]}) tends to infinity.  Such
arguments to the {\tt b*} functions can be set via the ellipses ({\tt
  ...})  arguments to {\tt optim.step.tgp}.\footnote{This runs
  contrary to how the ellipses are used by {\tt optim} in order to
  specify static arguments to {\tt f}.  If setting static arguments to
  {\tt f} is required within {\tt optim.step.tgp}, then they must be
  set in advance by adjusting the default arguments via {\tt
    formals}.}  A heuristic stopping criterion can be based on the
maximum improvement statistic obtained in each round as long as the
candidate locations become dense in the region as $T\rightarrow
\infty$.  This can be adjusted by increasing the {\tt NN} argument to
{\tt optim.step.tgp}.

The internal use of {\tt optim} within {\tt optim.step.tgp} on the
posterior predictive (kriging surrogate) surface via {\tt optim.ptgpf}
may proceed with any of the usual method arguments.  I.e.,
<<>>=
formals(optim)$method
@ 
however the default ordering is switched in {\tt optim.ptgpf}
and includes one extra method.
<<>>=
formals(optim.ptgpf)$method
@ 
Placing \verb!"L-BFGS-B"! in the default position is sensible since
this method enforces a rectangle of constraints as specified by {\tt
  rect}.  This guarentees that the additional candidate found by {\tt
  optim.ptfpf} will be valid.  However, the other {\tt optim} methods
generally work well despite that they do not enforce this constraint.
The final method, \verb!"optimize"!, applies only when the inputs to
{\tt f} are 1-d.  In this case, the documentation for {\tt optim}
suggests using the {\tt optimize} function instead.

%% Bobby says to Taddy: I also have a good example using the 1-d 
%% sinusoidal data, which is slightly more challenging since it
%% has two equivalent minima, but perhaps it is not ideal.  Do you
%% have any other good examples that we can test this on?  Do we 
%% need any more?
