\documentclass[nojss]{jss} \usepackage{amsmath} \usepackage{dcolumn} %%My declarations \newcommand{\hlink}{\htmladdnormallink} \newcommand{\MyPerp}{\perp \! \! \! \perp} \newcommand{\e}{\text{e-}} \newcolumntype{C}{>{$}c<{$}} \newcolumntype{L}{>{$}l<{$}} \newcolumntype{R}{>{$}r<{$}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% declarations for jss.cls %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% almost as usual \author{Walter R. Mebane, Jr. \\ University of Michigan \And Jasjeet S. Sekhon \\ UC Berkeley} %% for pretty printing and a nice hypersummary also set: \title{Genetic Optimization Using Derivatives: The \pkg{rgenoud} Package for \proglang{R}} \Plaintitle{Genetic Optimization Using Derivatives: The rgenoud Package for R} \Plainauthor{Walter R. Mebane, Jr., Jasjeet S. Sekhon} %% comma-separated \Shorttitle{\pkg{rgenoud}: Genetic Optimization Using Derivatives in \proglang{R}} %% an abstract and keywords \Abstract{ This introduction to the \proglang{R} package \pkg{rgenoud} is a modified version of \cite{MebaneJr.+Sekhon:2011}, published in the \textit{Journal of Statistical Software}. That version of the introduction contains higher resolution figures. \code{genoud} is an \proglang{R} function that combines evolutionary algorithm methods with a derivative-based (quasi-Newton) method to solve difficult optimization problems. \code{genoud} may also be used for optimization problems for which derivatives do not exist. \code{genoud} solves problems that are nonlinear or perhaps even discontinuous in the parameters of the function to be optimized. When the function to be optimized (for example, a log-likelihood) is nonlinear in the model's parameters, the function will generally not be globally concave and may have irregularities such as saddlepoints or discontinuities. Optimization methods that rely on derivatives of the objective function may be unable to find any optimum at all. Multiple local optima may exist, so that there is no guarantee that a derivative-based method will converge to the global optimum. On the other hand, algorithms that do not use derivative information (such as pure genetic algorithms) are for many problems needlessly poor at local hill climbing. Most statistical problems are regular in a neighborhood of the solution. Therefore, for some portion of the search space, derivative information is useful. The function supports parallel processing on multiple CPUs on a single machine or a cluster of computers.} \Keywords{genetic algorithm, evolutionary program, optimization, parallel computing, \proglang{R}} \Plainkeywords{genetic algorithm, evolutionary program, optimization, parallel computing, R} %% without formatting %% at least one keyword must be supplied %% publication information %% NOTE: This needs to filled out ONLY IF THE PAPER WAS ACCEPTED. %% If it was not (yet) accepted, leave them commented. \Volume{42} \Issue{11} %\Month{September} \Year{2011} \Submitdate{2007-02-11} \Acceptdate{2007-11-21} %% The address of (at least) one author should be given %% in the following format: \Address{ Walter R. Mebane, Jr.\\ % Professor\\ Department of Political Science\\ Department of Statistics\\ University of Michigan\\ Ann Arbor, MI 48109-1045\\ E-mail: \email{wmebane@umich.edu}\\ URL: \url{http://www.umich.edu/~wmebane} \\ [2ex] Jasjeet S. Sekhon\\ Department of Political Science\\ Department of Statistics\\ 210 Barrows Hall \#1950\\ UC Berkeley\\ Berkeley, CA 94720-1950\\ E-mail: \email{sekhon@berkeley.edu}\\ URL: \url{http://sekhon.berkeley.edu} } %% It is also possible to add a telephone and fax number %% before the e-mail in the following format: %% Telephone: +43/1/31336-5053 %% Fax: +43/1/31336-734 %% for those who use Sweave please include the following line (with % symbols): %% need no \usepackage{Sweave.sty} %% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \graphicspath{{Figures/}} \begin{document} \SweaveOpts{engine=R} %\VignetteIndexEntry{Using genoud} %\VignetteDepends{rgenoud} %\VignetteKeywords{rgenoud} %\VignettePackage{rgenoud} \section{Introduction} We developed the \proglang{R} package \pkg{rgenoud} to solve difficult optimization problems such as often arise when estimating nonlinear statistical models or solving complicated nonlinear, nonsmooth and even discontinuous functions.\footnote{The \pkg{rgenoud} software package is available from the Comprehensive \proglang{R} \citep{Rcore} Archive Network at \url{http://CRAN.R-project.org/package=rgenoud}.} Optimization difficulties often arise when the objective function (for instance, the log-likelihood) is a nonlinear function of the parameters. In such cases the function to be optimized is usually not globally concave. An objective function that is not globally concave may have multiple local optima, saddle points, boundary solutions or discontinuities. While the objective function for a statistical model is often concave in a neighborhood of the optimal solution, that neighborhood is often a small proportion of the parameter space of potential interest, and outside that neighborhood the function may be irregular. In such cases, methods of optimization that depend entirely on derivatives can be unreliable and often are virtually unusable. Newton-Raphson and quasi-Newton methods are among the commonly used optimization methods that rely completely on derivatives. Such methods work well when the function to be optimized is regular and smooth over the domain of parameter values that is of interest, but otherwise the methods often fail \citep{GillMurrayWright1981}. Even in models where such methods can be expected to work most of the time, resampling techniques such as the bootstrap \citep{EfronTibshirani1994} can generate resamples in which derivative-based optimization algorithms encounter severe difficulties. This is unfortunate because the methods most frequently used for optimization in problems of statistical estimation are entirely based on derivatives. The \pkg{rgenoud} package has been used by wide variety of users and developers. More than twelve \proglang{R} packages currently rely upon \pkg{rgenoud} including: \pkg{anchors} \citep[analyzing survey data with anchoring vignettes,][]{Wand+King:2011,WandKingLau2007,KingWand2007}; \pkg{BARD} \citep[automated redistricting,][]{Altman+McDonald:2011}; \pkg{boolean} \citep[boolean binary response models,][]{braumoeller2003,BraumoellerGoodrichKline2006}; \pkg{DiceOptim} \citep[kriging-based optimization for computer experiments,][]{DiceOptim}; \pkg{FAiR} \citep[factor analysis,][]{goodrich.fair}; \pkg{JOP} \citep[simultaneous optimization of multiple responses,][]{JOP}; \pkg{KrigInv} \citep[kriging-based inversion,][]{KrigInv}; \pkg{PKfit} \citep[data analysis in pharmacokinetics,][]{LeeLee2006}; \pkg{Matching} \citep[propensity score and multivariate matching,][]{SekhonMatching,Sekhon:2011}; \pkg{ivivc} \citep[in vitro-in vivo correlation modeling,][]{ivivc}; \pkg{multinomRob} \citep[robust multinomial models,][]{MebaneSekhon.multinomRob,MebaneSekhon2004}; and \pkg{Synth} \citep[synthetic control group method for comparative case studies,][]{AbadieGardeazabal2003,synth2008,Abadie+Diamond+Hainmueller:2011}. We present in Section \ref{sec:examples} an example using benchmark functions taken from \citet{yao.liu.lin1999}, followed by an example motivated by the \code{multinomRob} package. The benchmark suite includes functions that are high-dimensional, discontinuous or that have many local optima. The \code{multinomRob} package robustly estimates overdispersed multinomial regression models and uses \pkg{rgenoud} to solve a least quartile difference (LQD) generalized S-estimator \citep{MebaneSekhon2004}. The LQD is not a smooth function of the regression model parameters. The function is continuous, but the parameter space is finely partitioned by nondifferentiable boundaries. In another paper in this volume, the \pkg{Matching} package and its use of \pkg{rgenoud} are described in detail \citep{Sekhon:2011}. \pkg{Matching} provides functions for multivariate and propensity score matching and for finding optimal covariate balance based on a genetic search algorithm implemented in \pkg{rgenoud}. The search over optimal matches is discontinuous so no derivatives are used.\footnote{The \code{BFGS} option of \code{genoud} is set to \code{FALSE}, and the ninth operator which depends on derivatives is not used.} The search also involves lexical optimization which is a unique feature implemented in \pkg{rgenoud}.\footnote{Lexical optimization is useful when there are multiple fitness criteria; the parameters are chosen so as to maximize fitness values in lexical order---i.e., the second fit criterion is only considered if the parameters have the same fit for the first. See the \code{lexical} option and \citet{Sekhon:2011} for details. All of \code{genoud}'s options are described in the \proglang{R} help file for the function.} % If one wishes to use such a model, %significant optimization difficulties may need to be overcome. For %example, the estimating functions of all of the following models are %not in general globally concave: models with multivariate qualitative %variables, censoring or endogenous switching %\citep{maddala1983,amemiya1985}; linear covariance structures %\citep{bollen1989}; models for durations or transitions with %unobserved heterogeneity \citep{HeckmanSinger1985,lancaster1990}; %quasi-likelihood and extended quasi-likelihood models %\citep{McCullaghNelder1989} and systems based on such models %\citep{FahrmeirTutz1994}; and nonlinear simultaneous equation models %\citep{gallant1987}. For such models, optimization methods driven %entirely by derivatives are often unreliable. For some versions of %some of these models, algorithms such as EM that involve data %augmentation have been successfully applied to produce optimization %problems that derivative-based methods can solve %\citep[e.g.,][197]{lancaster1990}. The \pkg{rgenoud} package implements an updated and extended version of the \proglang{C} program \code{GENOUD} \citep{MebaneSekhon1997} described in \citep{SekhonMebane1998}. The many improvements include among other things the interface with \proglang{R}, which includes the ability to optimize functions written in \proglang{R}, options to optimize both floating point and integer-valued parameters, the ability to optimize loss functions which return multiple fitness values (lexical optimization), the ability to call \code{genoud} recursively, the ability to have the optimizer evaluate fits only for new parameter values, and the ability to use multiple computers, CPUs or cores to perform parallel computations. The \pkg{rgenoud} program combines an evolutionary algorithm with a quasi-Newton method. The quasi-Newton method is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method \citep[][119]{GillMurrayWright1981} implemented in \proglang{R}'s \code{optim} function. When the BFGS is being used, our program offers the option of using either \pkg{rgenoud}'s built-in numerical derivatives (which are based on code taken from \citealt[][337--344]{GillMurrayWright1981}) or user-supplied analytical derivatives.\footnote{User supplied derivatives may be provides via the \code{gr} option.} Our program can also work without the BFGS, in which case no derivatives are needed and the optimizer will work even when the function is discontinuous. The primary benefit from using derivatives is that the algorithm will then quickly find a local optimum when a current set of trial solution parameter values is in a smooth neighborhood of the local optimum point. Appropriate use of the BFGS can make the algorithm converge to the global optimum much more quickly. But premature or excessive use of the BFGS can prevent convergence to the global optimum.\footnote{The user can control whether \code{genoud} uses the BFGS at all (via the \code{BFGS} option), and if operators that use the BFGS are used (via the \code{P9} option), how often they are used.} As always, it is hazardous to rely on an optimizer's default settings. Our program does not eliminate the need for judgment, testing and patience. As Gill, Murray and Wright observe, ``there is no guaranteed strategy that will resolve every difficulty'' \citep[285]{GillMurrayWright1981}. In this article, we very briefly review the theory of random search algorithms that supports the assertion that \pkg{rgenoud} has a high probability of finding the global optimum when such exists. And we present three examples of how to use the \code{genoud} function: to optimize a simple but fiendish scalar Normal mixture model; to minimize a suite of benchmark functions that have previously been used to test evolutionary programming optimization algorithms; and to optimize a version of the only intermittently differentiable LQD estimator. Additional details on both the theory and performance of \code{genoud} can be found in our article that describes \code{GENOUD} \citep{SekhonMebane1998}. \section{Background on Genetic Optimization} An evolutionary algorithm (EA) uses a collection of heuristic rules to modify a population of trial solutions in such a way that each generation of trial values tends to be, on average, better than its predecessor. The measure for whether one trial solution is better than another is the trial solution's fitness value. In statistical applications, the fitness is a function of the summary statistic being optimized (e.g., the log-likelihood). \pkg{rgenoud} works for cases in which a solution is a vector of floating-point or integer numbers that serve as the parameters of a function to be optimized. The search for a solution proceeds via a set of heuristic rules, or \textit{operators}, each of which acts on one or more trial solutions from the current population to produce one or more trial solutions to be included in the new population. EAs do not require derivatives to exist or the function to be continuous in order find the global optimum. The EA in \pkg{rgenoud} is fundamentally a genetic algorithm (GA) in which the code-strings are vectors of numbers rather than bit strings, and the GA operators take special forms tuned for the floating-point or integer vector representation. A GA uses a set of randomized genetic operators to evolve a finite population of finite code-strings over a series of generations \citep{holland1975,goldberg1989,GrefenstetteBaker1989}. The operators used in GA implementations vary \citep{davis1991,FilhoTreleavenAlippi1994}, but in an analytical sense the basic set of operators can be defined as reproduction, mutation, crossover and inversion. The variation in these operators across different GA implementations reflects the variety of codes best suited for different applications. Reproduction entails selecting a code-string with a probability that increases with the code-string's fitness value. Crossover and inversion use pairs or larger sets of the selected code-strings to create new code-strings. Mutation randomly changes the values of elements of a single selected code-string. Used in suitable combinations, the genetic operators tend to improve average fitness of each successive generation, though there is no guarantee that average fitness will improve between every pair of successive generations. Average fitness may well decline. But theorems exist to prove that parts of the trial solutions that have above average fitness values in the current population are sampled at an exponential rate for inclusion in the subsequent population \citep[][139--140]{holland1975}. Each generation's population contains a biased sample of code-strings, so that a substring's performance in that population is a biased estimate of its average performance over all possible populations \citep{dejong1993,grefenstette1993}. The long-run properties of a GA may be understood by thinking of the GA as a Markov chain. A state of the chain is a code-string population of the size used in the GA. For code-strings of finite length and GA populations of finite size, the state space is finite. If such a GA uses random reproduction and random mutation, all states always have a positive probability of occurring. A finite GA with random reproduction and mutation is therefore approximately a finite and irreducible Markov chain.\footnote{\citet[][372--419]{feller1970} and \citet[][107--142]{billingsley1986} review the relevant properties of Markov chains. The randomness in an actual GA depends on the performance of pseudorandom number generators. This and the limitations of floating point arithmetic mean it is not literally true that an actual GA has a positive probability of reaching any state from any other state, and some states may in fact not be reachable from a given set of initial conditions.} An irreducible, finite Markov chain converges at an exponential rate to a unique stationary distribution \citep[][128]{billingsley1986}. This means that the probability that each population occurs rapidly converges to a constant, positive value. \citet{NixVose1992} and \citet{vose1993} use a Markov chain model to show that in a GA where the probability that each code-string is selected to reproduce is proportional to its observed fitness, the stationary distribution strongly emphasizes populations that contain code-strings that have high fitness values. They show that asymptotic in the population size---i.e., in the limit for a series of GAs with successively larger populations---populations that have suboptimal average fitness have probabilities approaching zero in the stationary distribution, while the probability for the population that has optimal average fitness approaches one. If $k>1$ populations have optimal average fitness, then in the limiting stationary distribution the probability for each approaches $1/k$. The theoretical results of Nix and Vose imply that a GA's success as an optimizer depends on having a sufficiently large population of code-strings. If the GA population is not sufficiently large, then the Markov chain that the GA approximately implements is converging to a stationary distribution in which the probabilities of optimal and suboptimal states are not sharply distinguished. Suboptimal populations can be as likely to occur as optimal ones. If the stationary distribution is not favorable, the run time in terms of generations needed to produce an optimal code-string will be excessive. For all but trivially small state spaces, an unfavorable stationary distribution can easily imply an expected running time in the millions of generations. But if the stationary distribution strongly emphasizes optimal populations, relatively few generations may be needed to find an optimal code-string. In general, the probability of producing an optimum in a fixed number of generations tends to increase with the GA population size. \begin{table} \begin{center} \line(1,0){450} \end{center} \begin{enumerate} \item[\code{P1}] Cloning. Copy $\mathbf{X}_{t}$ into the next generation, $\mathbf{X}_{t+1}$. \item[\code{P2}] Uniform Mutation. At random choose $i\in\mathbf{N}$. Select a value $\tilde{x}_i \sim U(\underline{x}_i, \overline{x}_i)$. Set $X_i=\tilde{x}_i$. \item[\code{P3}] Boundary Mutation. At random choose $i\in\mathbf{N}$. Set either $X_i=\underline{x}_i$ or $X_i=\overline{x}_i$, with probability $1/2$ of using each value. \item[\code{P4}] Non-uniform Mutation. At random choose $i\in\mathbf{N}$. Compute $p = (1-t/T)^B u$, where $t$ is the current generation number, $T$ is the maximum number of generations, $B>0$ is a tuning parameter and $u \sim U(0,1)$. Set either $X_i=(1-p)x_i + p\underline{x}_i$ or $X_i=(1-p)x_i + p\overline{x}_i$, with probability $1/2$ of using each value. \item[\code{P5}] Polytope Crossover. Using $m=\max(2,n)$ vectors $\mathbf{x}$ from the current population and $m$ random numbers $p_j\in(0,1)$ such that $\sum_{j=1}^m p_j = 1$, set $\mathbf{X} = \sum_{j=1}^m p_j \mathbf{x}_j$. \item[\code{P6}] Simple Crossover. Choose a single integer $i$ from $\mathbf{N}$. Using two parameter vectors, $\mathbf{x}$ and $\mathbf{y}$, set $X_i = p x_i + (1-p) y_i$ and $Y_i = p y_i + (1-p) x_i$, where $p\in(0,1)$ is a fixed number. \item[\code{P7}] Whole Non-uniform Mutation. Do non-uniform mutation for all the elements of $\mathbf{X}$. \item[\code{P8}] Heuristic Crossover. Choose $p \sim U(0,1)$. Using two parameter vectors, $\mathbf{x}$ and $\mathbf{y}$, compute $\mathbf{z} = p(\mathbf{x} - \mathbf{y}) + \mathbf{x}$. If $\mathbf{z}$ satisfies all constraints, use it. Otherwise choose another $p$ value and repeat. Set $\mathbf{z}$ equal to the better of $\mathbf{x}$ and $\mathbf{y}$ if a satisfactory mixed $\mathbf{z}$ is not found by a preset number of attempts. In this fashion produce two $\mathbf{z}$ vectors. \item[\code{P9}] Local-minimum Crossover. Choose $p \sim U(0,1)$. Starting with $\mathbf{x}$, run BFGS optimization up to a preset number of iterations to produce $\mathbf{\tilde{x}}$. Compute $\mathbf{z} = p\mathbf{\tilde{x}} + (1-p)\mathbf{x}$. If $\mathbf{z}$ satisfies boundary constraints, use it. Otherwise shrink $p$ by setting $p = p/2$ and recompute $\mathbf{z}$. If a satisfactory $\mathbf{z}$ is not found by a preset number of attempts, return $\mathbf{x}$. This operators is extremely computationally intensive, use sparingly. \end{enumerate} \begin{flushleft} Notation: \\ [2ex] $\mathbf{X} =\begin{bmatrix}X_1,\dots,X_n\end{bmatrix}$ is the vector of $n$ parameters $X_i$. $\underline{x}_i$ is the lower bound and $\overline{x}_i$ is the upper bound on values for $X_i$. $x_i$ is the current value of $X_i$, and $\mathbf{x}$ is the current value of $\mathbf{X}$. $\mathbf{N} = \{1,\dots,n\}$. $p \sim U(0,1)$ means that $p$ is drawn from the uniform distribution on the $[0,1]$ interval. \end{flushleft} \caption{`genoud' Operators. Adapted from \citet{SekhonMebane1998}.} \label{tab:operators} \begin{center} \line(1,0){450} \end{center} \end{table} The evolutionary algorithm in \pkg{rgenoud} uses nine operators that are listed in Table~\ref{tab:operators}. The operators extend and modify a set of operators used in \code{GENOCOP} \citep{MichalewiczSwaminathanLogan1993, Michalewicz1992}. The operators are numbered using syntax matching that used to refer to them by \pkg{rgenoud}. The \textit{cloning} operator simply makes copies of the best trial solution in the current generation (independent of this operator, \pkg{rgenoud} always retains one instance of the best trial solution). The \textit{uniform mutation}, \textit{boundary mutation} and \textit{non-uniform mutation} operators act on a single trial solution. Uniform mutation changes one parameter in the trial solution to a random value uniformly distributed on the domain specified for the parameter. Boundary mutation replaces one parameter with one of the bounds of its domain. Non-uniform mutation shrinks one parameter toward one of the bounds, with the amount of shrinkage decreasing as the generation count approaches the specified maximum number of generations. \textit{Whole non-uniform mutation} does non-uniform mutation for all the parameters in the trial solution. \textit{Heuristic crossover} uses two trial solutions to produce a new trial solution located along a vector that starts at one trial solution and points away from the other one. \textit{Polytope crossover} (inspired by simplex search, \citealt[][94--95]{GillMurrayWright1981}) computes a trial solution that is a convex combination of as many trial solutions as there are parameters. \textit{Simple crossover} computes two trial solutions from two input trial solutions by swapping parameter values between the solutions after splitting the solutions at a randomly selected point. This operator can be especially effective if the ordering of the parameters in each trial solution is consequential. \textit{Local-minimum crossover} computes a new trial solution in two steps: first it does a preset number of BFGS iterations starting from the input trial solution; then it computes a convex combination of the input solutions and the BFGS iterate. \section{Examples} \label{sec:examples} The only function in the \pkg{rgenoud} package is \code{genoud}. The interface of this function is similar to that of the \code{optim} function in \proglang{R}. But the function has many additional arguments that control the behavior of the evolutionary algorithm. \subsection{Asymmetric Double Claw:} Our first example, which we also studied in \citet{SekhonMebane1998}, is a normal mixture called the Asymmetric Double Claw (ADC). We plot the function in Figure~\ref{fig:mixture1}. \begin{figure} \caption{Normal Mixture: Asymmetric Double Claw} \vspace{-.5in} \label{fig:mixture1} \begin{center} \includegraphics{fig1.pdf} \end{center} \vspace{-.7in} \end{figure} Mathematically, this mixture is defined as \begin{eqnarray} \label{eq:adclaw} f_{\mathrm{ADC}} & = \sum_{m=0}^{1} \frac{46}{100} \mathrm{N}\left(2m-1,\frac{2}{3}\right) + \sum_{m=1}^{3} \frac{1}{300} \mathrm{N}\left(\frac{-m}{2},\frac{1}{100}\right) + \sum_{m=1}^{3} \frac{7}{300} \mathrm{N}\left(\frac{m}{2},\frac{7}{100}\right), \end{eqnarray} where N is the normal density. The asymmetric double claw is difficult to maximize because there are many local solutions. There are five local maxima in Figure \ref{fig:mixture1}. Standard derivative-based optimizers would simply climb up the hill closest to the starting value. To optimize this normal mixture we must first create a function for it <>= claw <- function(xx) { x <- xx[1] y <- (0.46*(dnorm(x,-1.0,2.0/3.0) + dnorm(x,1.0,2.0/3.0)) + (1.0/300.0)*(dnorm(x,-0.5,.01) + dnorm(x,-1.0,.01) + dnorm(x,-1.5,.01)) + (7.0/300.0)*(dnorm(x,0.5,.07) + dnorm(x,1.0,.07) + dnorm(x,1.5,.07))) return(y) } @ And we now make a call to \pkg{rgenoud} using this function: <>= library("rgenoud") claw1 <- genoud(claw, nvars=1, max=TRUE, pop.size=3000) @ The first argument of \code{genoud} is the function to be optimized. The first argument of that function must be the vector of parameters over which optimizing is to occur. Generally, the function should return a scalar result.\footnote{The function to be optimized may return a vector if one wishes to do lexical optimization. Please see the \code{lexical} option to \code{genoud}.} The second argument of \code{genoud} in this example (\code{nvars}) is the number of variables the function to be optimized takes. The third argument, \code{max=TRUE}, tells \code{genoud} to maximize the function instead of its default behavior which is to minimize. The fourth option \code{pop.size} controls the most important part of the evolutionary algorithm, the population size. This is the number of individuals \code{genoud} uses to solve the optimization problem. As noted in the theoretical discussion, the theorems related to evolutionary algorithms are asymptotic in the population size so larger is generally better but obviously takes longer. The maximum solution of the double claw density is reliably found by \code{genoud} even using the default value of \code{pop.size=1000}. Reliability does increase as the \code{pop.size} is made larger. Unfortunately, because of the stochastic nature of the algorithm, it is impossible to generally answer the question of what is the best population size to use. Other options determine the maximum number of generations the evolutionary algorithm computes. These options are \code{max.generations}, \code{wait.generations} and \code{hard.generation.limit}. The specified termination point also affects how some of the operators perform: the two non-uniform mutation operators introduce smaller ranges of variation in the trial solutions as the generation count approaches the specified \code{max.generations} value. There are many more options that can be used to fine-tune the behavior of the algorithm. The output printed by \code{genoud} is controlled by a \code{print.level} argument. The default value, \code{print.level=2}, produces relatively verbose output that gives extensive information about the set of operators being used and the progress of the optimization. Normally \proglang{R} conventions would suggest setting the default to be \code{print.level=0}, which would suppress output to the screen, but because \code{genoud} runs may take a long time, it can be important for the user to receive some feedback to see the program has not died and to be able to see where the program got stuck if it eventually fails to make adequate progress. The output printed by the preceding invocation of \code{genoud}, which uses the default value for a \code{print.level} argument, is as follows. \begin{CodeChunk} \begin{CodeOutput} Fri Feb 9 19:33:42 2007 Domains: -1.000000e+01 <= X1 <= 1.000000e+01 Data Type: Floating Point Operators (code number, name, population) (1) Cloning........................... 372 (2) Uniform Mutation.................. 375 (3) Boundary Mutation................. 375 (4) Non-Uniform Mutation.............. 375 (5) Polytope Crossover................ 375 (6) Simple Crossover.................. 376 (7) Whole Non-Uniform Mutation........ 375 (8) Heuristic Crossover............... 376 (9) Local-Minimum Crossover........... 0 HARD Maximum Number of Generations: 100 Maximum Nonchanging Generations: 10 Population size : 3000 Convergence Tolerance: 1.000000e-03 Using the BFGS Derivative Based Optimizer on the Best Individual Each Generation. Checking Gradients before Stopping. Using Out of Bounds Individuals. Maximization Problem. GENERATION: 0 (initializing the population) Fitness value... 4.112017e-01 mean............ 4.990165e-02 variance........ 9.708147e-03 #unique......... 3000, #Total UniqueCount: 3000 var 1: best............ 9.966758e-01 mean............ 3.453097e-02 variance........ 3.373681e+01 GENERATION: 1 Fitness value... 4.113123e-01 mean............ 2.237095e-01 variance........ 2.566140e-02 #unique......... 1858, #Total UniqueCount: 4858 var 1: best............ 9.995043e-01 mean............ 4.615946e-01 variance........ 7.447887e+00 [...] GENERATION: 10 Fitness value... 4.113123e-01 mean............ 2.953888e-01 variance........ 2.590842e-02 #unique......... 1831, #Total UniqueCount: 21708 var 1: best............ 9.995033e-01 mean............ 8.403935e-01 variance........ 5.363241e+00 GENERATION: 11 Fitness value... 4.113123e-01 mean............ 2.908561e-01 variance........ 2.733896e-02 #unique......... 1835, #Total UniqueCount: 23543 var 1: best............ 9.995033e-01 mean............ 8.084638e-01 variance........ 6.007372e+00 'wait.generations' limit reached. No significant improvement in 10 generations. Solution Fitness Value: 4.113123e-01 Parameters at the Solution (parameter, gradient): X[ 1] : 9.995033e-01 G[ 1] : -6.343841e-09 Solution Found Generation 1 Number of Generations Run 11 Fri Feb 9 19:33:45 2007 Total run time : 0 hours 0 minutes and 3 seconds \end{CodeOutput} \end{CodeChunk} %DESCRIBE OUTPUT After printing the date and time of the run, the program prints the domain of values it is allowing for each parameter of the function being optimized. In this case the default domain values are being used. Naturally it is important to specify domains wide enough to include the solution. In practice with highly nonlinear functions it is often better to specify domains that are relatively wide than to have domains that narrowly and perhaps even correctly bound the solution. This surprising behavior reflects the fact with a highly nonlinear function, a point that is close to the solution in the sense of simple numerical proximity may not be all that close in the sense of there being a short feasible path to get to the solution. Next the program prints the Data Type. This indicates whether the parameters of the function to be optimized are being treated as floating point numbers or integers. For more information about this, see the \code{data.type.int} argument. The program then displays the number of operators being used, followed by the values that describe the other characteristics set for this particular run: the maximum number of generations, the population size and the tolerance value to be used to determine when the parameter values will be deemed to have converged. The output then reports whether BFGS optimization will be applied to the best trial solution produced by the operators in each generation. For problems that are smooth and concave in a neighborhood of the global optimum, using the BFGS in this way can help \code{genoud} quickly converge once the best trial solution is in the correct neighborhood. This run of \code{genoud} will also compute the gradient at the best trial solution before stopping. In fact this gradient checking is used as a convergence check. The algorithm will not start counting its final set of generations (the \code{wait.generations}) until each element of the gradient is smaller in magnitude than the convergence tolerance. Gradients are never used and BFGS optimization is not used when the parameters of the function to be optimized are integers. The next message describes how strictly \code{genoud} is enforcing the boundary constraints specified by the domain values. By default (\code{boundary.enforcement=0}), the trial solutions are allowed to wander freely outside the boundaries. The boundaries are used only to define domains for those operators that use the boundary information. Other settings of the \code{boundary.enforcement} argument induce either more stringent or completely strict enforcement of the boundary constraints. Notice that the boundary constraints apply to the parameters one at a time. To enforce constraints that are defined by more complicated functional or data-dependent relationships, one can include an appropriate penalty function as part of the definition of the function to be optimized, letting that function define an extremely high (if minimizing) or low (if maximizing) value to be returned if the desired conditions are violated. After reporting whether it is solving a minimization or a maximization problem, \code{genoud} reports summary statistics that describe the distribution of the fitness value and the parameter values across the trial solutions at the end of each generation. In the default case where \code{genoud} is keeping track of all the distinct solutions considered over the whole course of the optimizing run, these generational reports also include a report of the number of unique trial solutions in the current population and the number of unique solutions ever considered. The benefit of keeping track of the solutions is to avoid repeatedly evaluating the function being optimized for the identical set of parameter values. This can be an important efficiency when function evaluations are expensive, as they can be in statistical applications where the data are extensive. This tracking behavior is controlled by the \code{MemoryMatrix} argument. Upon convergence, or when the hard maximum generation limit is reached, the program prints the fitness value at the best trial solution and that solution's parameter values. In this case the solution was found after one generation. While the Asymmetric Double Claw might present a difficult challenge for a gradient-based optimizer that uses only local hill climbing, it is an almost trivially simple problem for \code{genoud}. \subsection{Benchmark Functions:} The second example is a suite of 23 benchmark nonlinear functions used in \citet{yao.liu.lin1999} to study a pair of evolutionary programming optimization algorithms. Function definitions are in Table~\ref{tab:23funcs}. Because it includes a random component and so lacks a reproducible minimum value, we ignore the function numbered function 7 in their sequence.\footnote{The omitted function is $\sum_{i=1}^n ix_i^4 +U(0,1)$, where $U(0,1)$ is a uniformly distributed random variable on the unit interval that takes a new value whenever the function is evaluated. This stochastic aspect means that even given the set of parameters that minimize the nonstochastic component $\sum_{i=1}^n ix_i^4$, i.e., $x_i=0$, the value of the function virtually never attains the minimum possible value of zero. An optimizer that evaluated the function at $x_i=0$ would not in general obtain a function value smaller than the function value obtained for a wide range of different parameter values. Hence we do not consider this function to be a good test for function optimization algorithms.} Implementations of these functions are available in the supplemental \proglang{R} file provided with this article.\footnote{This supplemental file is available at \url{http://sekhon.berkeley.edu/rgenoud/func23.R}.} These \proglang{R} definitions include the values of the constants used in functions 14, 15 and 19 through 23. The function argument domains are restricted to the specific domains used by \citet{yao.liu.lin1999} via bounds that are stated in the list named \code{testbounds} in the supplemental file. \begin{table} \begin{center} \line(1,0){450} \end{center} \begin{tabular}{LLLL} \text{func.} & \text{definition} & n & \text{minimum}^a \\ 1 & \sum_{i=1}^n x_i^2 & 30 & 0 \\ 2 & \sum_{i=1}^n |x_i| + \prod_{i=1}^n |x_i| & 30 & 0 \\ 3 & \sum_{i=1}^n (\sum_{j=1}^i x_j)^2 & 30 & 0 \\ 4 & \max_i \{|x_i|, 1\leq i\leq n\} & 30 & 0 \\ 5 & \sum_{i=1}^{n-1} [100(x_{i+1}-x_i^2)^2 +(x_i-1)^2], & 30 & 0 \\ 6 & \sum_{i=1}^n (\lfloor x_i + 0.5\rfloor)^2 & 30 & 0 \\ 7 & \sum_{i=1}^n ix_i^4 +U(0,1) & 30 & 0 \\ 8 & \sum_{i=1}^n -x_i\sin(\sqrt{|x_i|}) & 30 & -12569.5 \\ 9 & \sum_{i=1}^n [x_i^2-10\cos(2\pi x_i)+10] & 30 & 0 \\ 10 & -20\exp\left(-0.2\sqrt{n^{-1}\sum_{i=1}^n x_i^2}\,\right) \\ & \qquad -\exp(n^{-1}\sum_{i=1}^n\cos 2\pi x_i)+20+e & 30 & 0 \\ 11 & (1/1000) \sum_{i=1}^n x_i^2 - \prod_{i=1}^n\cos\left(x_i/\sqrt{i}\right) + 1 & 30 & 0 \\ 12 & n^{-1}\pi \left\{10 \sin^2(\pi y_1)+\sum_{i=1}^{n-1}(y_i-1)^2 [1+10\sin^2(\pi y_{i+1})]+(y_n-1)^2\right\} \\ & \quad +\sum_{i=1}^n u(x_i,10,100,4), & 30 & 0 \\ & y_i=1+(x_i+1)/4, \quad u(x_i,a,k,m)=\begin{cases} k(x_i-a)^m, & x_i>a, \\ 0, & -a\leq x_i\leq a,\\ k(-x_i-a)^m, & x_i<-a, \end{cases} & \\ 13 & \left\{ \sin^2(3\pi x_1)+\sum_{i=1}^{n-1}(x_i-1)^2[1+\sin^2(3\pi x_{i+1})]\right\}/10 \\ & \quad +(x_n-1)[1+\sin^2(2\pi x_n)]/10 +\sum_{i=1}^n u(x_i,5,100,4) & 30 & 0 \\ 14 & \left\{1/500+ \sum_{j=1}^{25} 1/\left[j+\sum_{i=1}^2(x_i-a_{ij})^6\right]\right\}^{-1} & 2 & 1 \\ 15 & \sum_{i=1}^{11} [a_i-x_1(b_i^2+b_ix_2)/(b_i^2+b_ix_3+x_4)]^2 & 4 & 0.0003075 \\ 16 & 4x_1^2-2.1x_1^4+x_1^6/3+x_1x_2-4x_2^2+4x_2^4 & 2 & -1.0316285 \\ 17 & [x_2-5.1x_1^2/(4\pi^2) + 5x_1/\pi - 6]^2 + 10[1-1/(8\pi)]\cos(x_1) + 10 & 2 & 0.398 \\ 18 & [1+(x_1+x_2+1)^2(19-14 x_1+3 x_1^2-14 x_2+6 x_1 x_2+3 x_2^2)] \\ & \quad \times [30+(2 x_1-3 x_2)^2 (18-32 x_1+12 x_1^2+48 x_2-36 x_1 x_2+27 x_2^2)] & 2 & 3 \\ 19 & -\sum_{i=1}^4 c_i\exp\left[-\sum_{j=1}^4 a_{ij}(x_i-p_{ij})^2\right] & 4 & -3.86 \\ 20 & -\sum_{i=1}^4 c_i\exp\left[-\sum_{j=1}^6 a_{ij}(x_i-p_{ij})^2\right] & 6 & -3.32 \\ 21 & \sum_{i=1}^5 [(x-a_i)^{\prime}(x-a_i)+c_i]^{-1} & 4 & -10 \\ 22 & \sum_{i=1}^7 [(x-a_i)^{\prime}(x-a_i)+c_i]^{-1} & 4 & -10 \\ 23 & \sum_{i=1}^{10} [(x-a_i)^{\prime}(x-a_i)+c_i]^{-1} & 4 & -10 \end{tabular} \begin{flushleft} Notes: $^a$Minimum function value within specified bounds as given by \citet[85]{yao.liu.lin1999}. \end{flushleft} \vspace*{.1in} \caption{23 Benchmark Functions} \label{tab:23funcs} \begin{center} \line(1,0){450} \end{center} \end{table} As \citet{yao.liu.lin1999} describe, optimizing each of functions 1--13 presents a high-dimensional problem. These functions each have $n=30$ free parameters. Functions 1--5 are unimodal, with function 5 being a 30-dimensional version of the banana-shaped Rosenbrock function. Function 6 is a step function. Function 6 has one minimum value that occurs when all arguments are in the interval $x_i\in[0,.5)$, and the function is discontinuous. Functions 8--13 are multimodal, defined such that the number of local minima increases exponentially with the number of arguments. \citet[84]{yao.liu.lin1999} describe these functions as among ``the most difficult class of problems for many optimization algorithms (including [evolutionary programming]).'' Functions 14--23, which have between two and six free parameters each, each have only a few local minima. Nonetheless the evolutionary programming algorithms considered by \citet{yao.liu.lin1999} have trouble optimizing functions 21--23. Although \citet[85]{yao.liu.lin1999} state that each of these functions has a minimum value of $-10$, over 50 replications the two algorithms they consider achieve solutions averaging between $-5.52$ and $-9.10$ \citep[88, Table~IV]{yao.liu.lin1999}. We use these benchmark functions to illustrate not only how effective \code{genoud} can be with a range of difficult problems, but also to emphasize an important aspect of how one should think about trying to tune the optimizer's performance. Theory regarding genetic algorithms suggests that optimal solutions are more likely to appear as both the population size of candidate solutions and the number of generations increase. In \code{genoud} two arguments determine the number of generations. One is \code{max.generations}: if \code{hard.generation.limit=TRUE}, then the value specified for \code{max.generations} is a binding upper limit. The other is \code{wait.generations}, which determines when the algorithm terminates if \code{hard.generation.limit=FALSE}. But even if \code{hard.generation.limit=TRUE}, then \code{wait.generations} determines for how many generations the algorithm continues once the best parameter vector and the value of the function being optimized appear to have settled down. The fact that the current best solution is not changing should not be treated as decisive, because this solution may be merely a local optimum or a saddlepoint. If the population size is sufficient, the algorithm tends to build a population of trial solutions that contain parameters in neighborhoods of all the competitive local optima in the domain defined by the parameter boundaries. Even while the current best solution is stable, the algorithm is improving the solutions near other local optima. So having a higher \code{wait.generations} value never worsens the algorithm's efficacy. Increasing \code{max.generations} may or may not in itself improve optimization. The value of \code{max.generations} sets the value of $T$ used in the mutation operators---operators 4 and 7 in Table~\ref{tab:operators}. These mutation operators perform random search near a trial solution that has been selected for mutation only when the current generation count is an appreciable fraction of $T$. So increasing \code{max.generations} without changing \code{wait.generations} increases the period during which random search is occurring over a wider domain. For multimodal functions such a wider search may be helpful, but sometimes failing to search more densely near the current trial solutions is not good. We use \code{genoud} to minimize the 23 functions using two values for \code{pop.size} (5000 and 10000) and two values for \code{max.generations} (30 and 100). Following \citet{yao.liu.lin1999}, we replicate each optimization 50 times. The following code describes the computations. The list \code{testfuncs}, vector \code{testNparms} and list \code{testbounds} are defined in the supplemental \proglang{R} file, and it is assume that this file is loaded with the \code{source("supplemental.R")} command. The vector \code{gradcheck} is true for all elements except the ones corresponding to functions 6 and 7. <>= source("supplemental.R") gradcheck <- rep(TRUE,23) gradcheck[6:7] <- FALSE sizeset <- c(5000,10000) genset <- c(30,100) nreps <- 50 gsarray <- array(NA, dim=c(length(sizeset), length(genset), 23, nreps)) asc <- function(x) { as.character(x) } dimnames(gsarray) <- list(asc(sizeset), asc(genset), NULL, NULL); for (gsize in sizeset) { for (ngens in genset) { for (i in 1:23) { for (j in 1:nreps) { gsarray[as.character(gsize), as.character(ngens),i,j] <- genoud(testfuncs[[i]], nvars=testNparms[i], pop.size=gsize, max.gen=ngens, hard.gen=TRUE, Domains=testbounds[[i]], solution.tol=1e-6, boundary=1, gradient.check=gradcheck[i], print=0)$value } } } } @ Using \code{genoud} to minimize the benchmark functions produces excellent results, at least when the \code{pop.size} and \code{max.generations} arguments are sufficiently large. Table~\ref{tab:23mean} reports the mean function values for each configuration of the arguments. These values may be compared both to the true function minima given by \citet{yao.liu.lin1999} (see the rightmost column in Table~\ref{tab:23funcs}) and to the average minimum function values \citet{yao.liu.lin1999} report for their ``fast'' evolutionary programming (FEP) algorithm, which appear in the last column of Table~\ref{tab:23mean}. The \code{genoud} averages for the \code{max.gen=100} configurations equal or are close to the true minima for all the functions except function 13. One can reasonably argue that the average solutions for function 5 are not as close to zero as might be desired: these averages are close to $10^{-7}$, while the averages for other functions that have a true minimum of zero are $10^{-15}$ or smaller. And the averages for functions 6, 12 and 15 in the \code{pop.size=5000} case are off. The effect of increasing \code{pop.size} is apparent with respect to both those three functions and also functions 13 and 20--23: the average minima are smaller with \code{pop.size=10000} than with \code{pop.size=5000}. Except for functions 6 and 12 in the \code{pop.size=5000} case and function 13, all the \code{genoud} solution averages for \code{max.gen=100} are either slightly or substantially better than the corresponding FEP solution averages. The results in Table~\ref{tab:23mean} clearly illustrate the potential consequences of not allowing \code{genoud} to run for a sufficient number of generations. While some of the \code{genoud} solutions for \code{max.gen=30} have competitive means, several of the means are not good at all. \begin{table} \begin{center} \line(1,0){450} \end{center} \begin{tabular}{llllll} & \multicolumn{4}{c}{\code{genoud}} & \\ \cline{2-5} & \multicolumn{2}{c}{\code{pop.size=5000}$^a$} & \multicolumn{2}{c}{\code{pop.size=10000}$^a$} & \\ \text{func.} & \code{max=30} & \code{max=100} & \code{max=30} & \code{max=100} & \text{FEP}$^b$ \\ 1 & 1.453\e 32 & 1.658\e 32 & 2.416\e 32 & 6.134\e 32 & 5.7\e 4 \\ 2 & 6.55\e 16 & 7.212\e 16 & 9.652\e 16 & 1.043\e 15 & 8.1\e 3 \\ 3 & 4.633\e 18 & 3.918\e 18 & 3.787\e 18 & 4.032\e 18 & 1.6\e 2 \\ 4 & 6.203\e 17 & 6.542\e 17 & 9.453\e 17 & 7.85\e 17 & 0.3 \\ 5 & 0.07973 & 5.887\e 08 & 8.133\e 08 & 8.917\e 08 & 5.06 \\ 6 & 18.58 & 0.08 & 9.38 & 0 & 0 \\ % 7 & 0.1114 & 0.07473 & 0.1114 & 0.06466 & 7.6\e 3 \\ 8 & -12569.49 & -12569.49 & -12569.49 & -12569.49 & -12554.5 \\ 9 & 2.786 & 0 & 0.9353 & 0 & 4.6\e 2 \\ 10 & 2.849 & 3.997\e 15 & 2.199 & 3.997\e 15 & 1.8\e 2 \\ 11 & 7.994\e 17 & 7.105\e 17 & 9.548\e 17 & 6.439\e 17 & 1.6\e 2 \\ 12 & 5.371\e 19 & 0.004147 & 0.002073 & 1.178\e 19 & 9.2\e 6 \\ 13 & 0.02095 & 0.006543 & 0.006629 & 0.003011 & 1.6\e 4 \\ 14 & 0.998 & 0.998 & 0.998 & 0.998 & 1.22 \\ 15 & 0.0003441 & 0.0004746 & 0.0003807 & 0.0003807 & 5.0\e 4 \\ 16 & -1.0316285 & -1.0316285 & -1.0316285 & -1.0316285 & -1.03 \\ 17 & 0.3979 & 0.3979 & 0.3979 & 0.3979 & 0.398 \\ 18 & 3 & 3 & 3 & 3 & 3.02 \\ 19 & -3.863 & -3.863 & -3.863 & -3.863 & -3.86 \\ 20 & -3.274 & -3.277 & -3.279 & -3.286 & -3.27 \\ 21 & -9.85 & -9.444 & -9.95 & -10.05 & -5.52 \\ 22 & -9.771 & -10.09 & -10.19 & -10.3 & -5.52 \\ 23 & -10.1 & -9.997 & -10.32 & -10.21 & -6.57 \end{tabular} \begin{flushleft} Note: $^b$Average minimum function values (over 50 replications) obtained using \code{genoud}. $^a$Mean best function values (over 50 replications) reported for the ''fast'' evolutionary programming algorithm, from \citet[85 and 88, Tables II--IV]{yao.liu.lin1999}. \end{flushleft} \vspace*{.1in} \caption{Mean Values of 22 Optimized Functions} \label{tab:23mean} \begin{center} \line(1,0){450} \end{center} \end{table} The effect of increasing \code{pop.size} are even more clearly apparent in Table~\ref{tab:23sdev}, which reports the standard deviations of the respective minima across the 50 replications. With the exceptions of functions 6, 12, 13, 15 and 21 with \code{pop.size=5000}, the \code{genoud} solutions for \code{max.gen=100} vary much less than the corresponding FEP solutions. For those functions and also for functions 20, 22 and 23, the \code{max.gen=100} solutions with \code{pop.size=10000} vary noticeably less than the solutions with \code{pop.size=5000}. \begin{table} \begin{center} \line(1,0){450} \end{center} \begin{tabular}{llllll} & \multicolumn{4}{c}{\code{genoud}} & \\ \cline{2-5} & \multicolumn{2}{c}{\code{pop.size=5000}$^a$} & \multicolumn{2}{c}{\code{pop.size=10000}$^a$} & \\ \text{func.} & \code{max=30} & \code{max=100} & \code{max=30} & \code{max=100} & \text{FEP}$^b$ \\ 1 & 9.997\e 32 & 7.059\e 32 & 9.562\e 32 & 2.1\e 31 & 1.3\e 4 \\ 2 & 1.668\e 15 & 1.621\e 15 & 2.116\e 15 & 2.102\e 15 & 7.7\e 4 \\ 3 & 4.568\e 18 & 3.342\e 18 & 4.38\e 18 & 5.136\e 18 & 1.4\e 2 \\ 4 & 1.793\e 16 & 1.758\e 16 & 2.055\e 16 & 2.002\e 16 & 0.5 \\ 5 & 0.5638 & 4.921\e 08 & 5.573\e 08 & 4.955\e 08 & 5.87 \\ 6 & 5.65 & 0.274 & 3.528 & 0 & 0 \\ % 7 & 0.0837 & 0.05153 & 0.08042 & 0.04883 & 2.6\e 3 \\ 8 & 3.749\e 10 & 1.071\e 12 & 8.948\e 09 & 6.365\e 13 & 52.6 \\ 9 & 1.864 & 0 & 1.179 & 0 & 1.2\e 2 \\ 10 & 0.7146 & 0 & 0.702 & 0 & 2.1\e 3 \\ 11 & 1.209\e 16 & 1.582\e 16 & 1.289\e 16 & 8.713\e 17 & 2.2\e 2 \\ 12 & 2.336\e 18 & 0.02052 & 0.01466 & 7.423\e 19 & 3.6\e 6 \\ 13 & 0.03427 & 0.006867 & 0.006903 & 0.001508 & 7.3\e 5 \\ 14 & 5.638\e 12 & 8.894\e 11 & 1.029\e 12 & 4.35\e 12 & 0.56 \\ 15 & 0.0001813 & 0.0003546 & 0.000251 & 0.0002509 & 3.2\e 4 \\ 16 & 1.315\e 14 & 9.751\e 15 & 1.233\e 14 & 1.054\e 14 & 4.9\e 7 \\ 17 & 5.422\e 15 & 5.51\e 15 & 4.925\e 15 & 1.392\e 14 & 1.5\e 7 \\ 18 & 1.509\e 13 & 3.477\e 14 & 6.18\e 14 & 2.907\e 14 & 0.11 \\ 19 & 7.349\e 15 & 1.521\e 15 & 1.344\e 14 & 7.255\e 15 & 1.4\e 5 \\ 20 & 0.05884 & 0.0583 & 0.05765 & 0.05504 & 5.9\e 2 \\ 21 & 1.212 & 1.776 & 1.005 & 0.7145 & 1.59 \\ 22 & 1.937 & 1.269 & 1.052 & 0.7459 & 2.12 \\ 23 & 1.479 & 1.636 & 1.066 & 1.29 & 3.14 \end{tabular} \begin{flushleft} Note: $^a$Standard deviation of the minimum function values (over 50 replications) obtained using \code{genoud}. $^b$Standard deviation of the best function values (over 50 replications) reported for the ''fast'' evolutionary programming algorithm, from \citet[85 and 88, Tables II--IV]{yao.liu.lin1999}. \end{flushleft} \vspace*{.1in} \caption{Standard Deviations of Values of 22 Optimized Functions} \label{tab:23sdev} \begin{center} \line(1,0){450} \end{center} \end{table} \subsection{A Logistic Least Quartile Difference Estimator:} Our third example is a version of the LQD estimator used in \code{multinomRob}. Using the \proglang{R} function \code{IQR} to compute the interquartile range, the function to be minimized may be defined as follows.\footnote{The LQD problem solved in \code{multinomRob} is somewhat different. There the problem is to minimize the $\binom{h_{K}}{2}$ order statistic of the set of absolute differences among the standardized residuals, where $h_K$ is a function of the sample size and the number of unknown regression model coefficients \citep{MebaneSekhon2004}. The problem considered in the current example is simpler but exhibits similar estimation difficulties.} <>= LQDxmpl <- function(b) { logistic <- function(x) { 1/(1+exp(-x)) } sIQR <- function(y, yhat, n) { IQR((y-yhat)/sqrt(yhat*(n-yhat)), na.rm=TRUE) } sIQR(y, m*logistic(x %*% b), m) } @ For this example we define \code{LQDxmpl} after we compute the simulated data, so the data vector \code{y}, matrix \code{x} and scalar \code{m} are in scope to evaluate to have the values we simulate: <>= m <- 100 x <- cbind(1,rnorm(1000),rnorm(1000)) b1 <- c(.5, 1, -1) b2 <- c(0, -1, 1) logistic <- function(x) { 1/(1+exp(-x)) } y <- rbinom(1000, m, logistic(c(x[1:900,] %*% b1, x[901:1000,] %*% b2))) @ The data simulate substantial contamination. The first 900 observations are generated by one binomial regression model while the last 100 observations come from a very different model. Presuming we are interested in the model that generated the bulk of the data, ignoring the contamination in a generalized linear model with the binomial family produces very poor results: <>= summary(glm1 <- glm(cbind(y,m-y) ~ x[,2] + x[,3], family="binomial")) @ \begin{CodeChunk} \begin{CodeOutput} Call: glm(formula = cbind(y, m - y) ~ x[, 2] + x[, 3], family = "binomial") Deviance Residuals: Min 1Q Median 3Q Max -22.9168 -1.1693 0.3975 1.5895 24.6439 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.439492 0.007097 61.93 <2e-16 *** x[, 2] 0.679847 0.007985 85.14 <2e-16 *** x[, 3] -0.716963 0.007852 -91.31 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \end{CodeOutput} \end{CodeChunk} Of course, if we knew which observations to omit the results would be much better: <>= suby <- y[1:900] subx <- x[1:900,] summary(glm2 <- glm(cbind(suby,m-suby) ~ subx[,2] + subx[,3], family="binomial")) @ \begin{CodeChunk} \begin{CodeOutput} Call: glm(formula = cbind(suby, m - suby) ~ subx[, 2] + subx[, 3], family = "binomial") Deviance Residuals: Min 1Q Median 3Q Max -3.21478 -0.71699 0.03528 0.67867 2.88314 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.501880 0.008036 62.46 <2e-16 *** subx[, 2] 1.003592 0.009779 102.63 <2e-16 *** subx[, 3] -0.984295 0.009437 -104.30 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \end{CodeOutput} \end{CodeChunk} But in practical applications it is unknown apriori which observations should be treated as outliers. As the definition of \code{LQDxmpl} indicates, the LQD is based on minimizing the interquartile range (IQR) of the standardized residuals. Because the quartiles correspond to different data points for different values of the regression coefficients, the fitness function is not smooth, which is to say it is not everywhere differentiable. In general, at every point in the parameter space where the identity of the first or third quartile point changes, the function is not differentiable. Figure~\ref{fig:lqd1} illustrates this. A higher resolution version of this figure is available in \cite{MebaneJr.+Sekhon:2011}---also see \url{http://sekhon.berkeley.edu/papers/rgenoudJSS.pdf}. The IQR clearly has a minimum in the vicinity of the coefficient values used to generate most of the data. But contour plots for the numerically evaluated partial derivative with respect to the second coefficient parameter testify to the function's local irregularity. The function we use to evaluate this numerical derivative is defined as follows. <>= dLQDxmpl <- function(b) { eps <- 1e-10 logistic <- function(x) { 1/(1+exp(-x)) } sIQR <- function(y, yhat, n) { IQR((y-yhat)/sqrt(yhat*(n-yhat)), na.rm=TRUE) } dsIQR <- vector() for (i in 1:length(b)) { beps <- b beps[i] <- b[i]+eps dsIQR <- c(dsIQR, (sIQR(y, m*logistic(x %*% beps), m)- sIQR(y, m*logistic(x %*% b), m))/eps) } return(dsIQR) } @ Setting the intercept equal to 0.5, the code to generate the plotted values is % -rw-rw---- 1 wrm1 research 1878 Feb 10 2007 lqdxmpl2plotB.R <>= blen <- 3 lenbseq <- length(bseq <- seq(-2,2,length=200)) bseq3 <- seq(-1.2,-.9,length=200) bseq2 <- seq(.89,1.1,length=200) IQRarr <- IQRarrA <- array(NA, dim=c((1+blen), lenbseq, lenbseq)) dimnames(IQRarrA) <- list(NULL, as.character(bseq), as.character(bseq)) dimnames(IQRarr) <- list(NULL, as.character(bseq2), as.character(bseq3)) for (i in 1:lenbseq) { for (j in 1:lenbseq) { IQRarrA[1,i,j] <- LQDxmpl(c(.5, bseq[i], bseq[j])) IQRarrA[-1,i,j] <- dLQDxmpl(c(.5, bseq[i], bseq[j])) IQRarr[1,i,j] <- LQDxmpl(c(.5, bseq2[i], bseq3[j])) IQRarr[-1,i,j] <- dLQDxmpl(c(.5, bseq2[i], bseq3[j])) } } @ The following code produces the plots: <>= par(mfrow=c(2,2), lwd=.1) contour(bseq,bseq, IQRarrA[1,,], main="IQR", xlab="b[2]", ylab="b[3]") contour(bseq,bseq, IQRarrA[3,,], main="partial derivative w/r/t b[2]", xlab="b[2]", ylab="b[3]") loc2 <- (150:160)-5 loc3 <- (40:50)+5 contour(bseq[loc2],bseq[loc3], IQRarrA[3,loc2,loc3], main="partial derivative w/r/t b[2]", xlab="b[2]", ylab="b[3]") contour(bseq2,bseq3, IQRarr[3,,], main="partial derivative w/r/t b[2]", xlab="b[2]", ylab="b[3]") @ If the IQR function were smooth, we would see clearly separated, slightly curved contour lines, reflecting the nonlinearity of the \code{logistic} function, but there is nothing like that. Instead, looking over the domain $[-2,2]^2$ for the second and third regression coefficient parameters, with 200 points evaluated along each axis (the upper right plot), there is a splotchy cloud. This reflects the fact that the derivative changes sign very frequently over the domain: of the 40,000 points at which the derivative is evaluated, it is positive at 12,460 points and negative at 27,540 points. \begin{figure}[t] \caption{Contour Plots of the LQD Fitness Function and of its Partial Derivatives} \label{fig:lqd1} \begin{center} \includegraphics{lqdxmpl2B.pdf} % -rw-rw---- 1 wrm1 research 2745 Feb 10 23:15 lqdxmpl2plotB.Rout \end{center} \vspace{-.3in} \end{figure} The LQD fitness function is not appreciably smoother close to the true values. The bottom two plots show the partial derivatives with respect to the second coefficient parameter evaluated over the domain $[.89,1.1]\times[-1.2,-.9]$. The bottom left plot evaluates the derivatives at 11 points along each axis while the bottom right plot uses 200 points along each axis. In the left plot it is easier to see the intricacy of the partitioning of the parameter space as the identity of the first or third quartile point changes. The bottom right plot shows the intricacy in fact replicates at the finer grain of the more detailed evaluations (to see this it is probably necessary to magnify the plot while viewing online the \textit{Journal of Statistical Software} version of this documentation).\footnote{Please see \cite{MebaneJr.+Sekhon:2011} or \url{http://sekhon.berkeley.edu/papers/rgenoudJSS.pdf}.} Within this smaller domain the sign of the derivative changes even more frequently than it does over the domain $[-2,2]^2$: the derivative is positive at 18,098 points and negative at 21,902 points. While the LQD fitness function may be differentiable in a neighborhood of the global solution, that neighborhood, if it exists, is clearly not very big. As likely is that the global solution is located at a point where the function is not differentiable. Hence a numerically evaluated gradient is not meaningful for evaluating whether the solution has been found. At the true solution, numerical gradient values may differ substantially from zero. To use the LQD to estimate the binomial regression model parameters we use the following call to \code{genoud}. Because gradient information is of questionable relevance for this problem, we turn off the termination condition that the gradient at the solution be smaller than the value specified for the \code{solution.tolerance} argument. We retain the default setting \code{BFGS=TRUE} because, in principle, gradient-driven optimization may help in each of the many differentiable neighborhoods, even if it is useless across the nondifferentiable boundaries. Our experience optimizing the LQD \citep{MebaneSekhon2004} shows that using the BFGS in this way improves performance, even though the gradient is not useful for evaluating whether the solution has been found. <>= LQD1 <- genoud(LQDxmpl, nvars=3, max=FALSE, pop.size=2000, max.generations=300, wait.generations=100, gradient.check=FALSE, print=1) @ This invocation of the \code{LQDxmpl} function matches the behavior of \code{multinomRob} in that it produces an estimate for the intercept parameter along with the other coefficients. In a linear regression context, the interquartile range statistic contains no information about the intercept, so the LQD is not an estimator for that parameter. With a binomial regression model there is some information about the intercept due to the nonlinearity of the \code{logistic} function. The LQD estimate for the intercept should nonetheless be expected not to be very good. Results from the preceding estimation are as follows. % macht:/home/xchg/jss07/rgenoud: % -rw-rw---- 1 wrm1 research 5080 Aug 7 02:31 lqdxmpl1b.Rout \begin{CodeChunk} \begin{CodeOutput} Tue Aug 7 02:27:08 2007 Domains: -1.000000e+01 <= X1 <= 1.000000e+01 -1.000000e+01 <= X2 <= 1.000000e+01 -1.000000e+01 <= X3 <= 1.000000e+01 [...] HARD Maximum Number of Generations: 300 Maximum Nonchanging Generations: 100 Population size : 2000 Convergence Tolerance: 1.000000e-06 Using the BFGS Derivative Based Optimizer on the Best Individual Each Generation. Not Checking Gradients before Stopping. Using Out of Bounds Individuals. Minimization Problem. Generation# Solution Value 0 4.951849e-01 56 1.298922e-01 57 1.298891e-01 59 1.298848e-01 60 1.298820e-01 61 1.298793e-01 62 1.298768e-01 63 1.298744e-01 'wait.generations' limit reached. No significant improvement in 100 generations. Solution Fitness Value: 1.298743e-01 Parameters at the Solution (parameter, gradient): X[ 1] : 8.130357e-02 G[ 1] : 9.616492e-03 X[ 2] : 8.889485e-01 G[ 2] : -1.167897e-01 X[ 3] : -9.327966e-01 G[ 3] : -3.090130e-02 Solution Found Generation 63 Number of Generations Run 164 Tue Aug 7 02:31:09 2007 Total run time : 0 hours 4 minutes and 1 seconds \end{CodeOutput} \end{CodeChunk} Recall that the gradient is not reliably informative at the solution. To check whether this solution is believable, we might try reestimating the model using a larger population and larger specified number of generations: <>= LQD1 <- genoud(LQDxmpl, nvars=3, max=FALSE, pop.size=10000, max.generations=1000, wait.generations=300, gradient.check=FALSE, print=1) @ At the price of a greatly increased running time (from four minutes up to one hour 53 minutes), the results are better than the first run (even though the summary measure of fit is slightly worse): % macht:/home/xchg/jss07/rgenoud: % -rw-rw---- 1 wrm1 research 5236 Aug 7 05:16 lqdxmpl1d.Rout \begin{CodeChunk} \begin{CodeOutput} Minimization Problem. Generation# Solution Value 0 2.238865e-01 2 1.301149e-01 3 1.300544e-01 4 1.300482e-01 6 1.300375e-01 7 1.300343e-01 8 1.300323e-01 134 1.299662e-01 135 1.299099e-01 136 1.298867e-01 137 1.298843e-01 138 1.298822e-01 139 1.298791e-01 141 1.298774e-01 'wait.generations' limit reached. No significant improvement in 300 generations. Solution Fitness Value: 1.298770e-01 Parameters at the Solution (parameter, gradient): X[ 1] : 2.013748e-01 G[ 1] : -7.394125e-02 X[ 2] : 9.526390e-01 G[ 2] : 7.807607e-02 X[ 3] : -9.642458e-01 G[ 3] : 3.834052e-02 Solution Found Generation 141 Number of Generations Run 442 Tue Aug 7 05:16:37 2007 Total run time : 1 hours 53 minutes and 45 seconds \end{CodeOutput} \end{CodeChunk} %$ This example demonstrates a key difficulty that arises when optimizing irregular functions in the absence of gradients. It is difficult to assess when or whether an optimum has been found. The estimated coefficient values are close to the ones used to generate most of the data, except as expected the estimate for the intercept is not good. The estimates are better than if we had ignored the possibility of contamination. But whether these are the best possible estimates is not clear. If we were to use an even larger population and specify that an even greater number of generations be run, perhaps a better solution would be found. %http://sekhon.berkeley.edu/rgenoud/R/genoud_cluster_manual.R %http://sekhon.berkeley.edu/rgenoud/R/genoud_cluster_manual_tunnel.R %\code{genoud} can be tricked into doing constrained optimization by %passing in a function which contains a penalty function. Even for less irregular problems convergence is difficult to determine. Nonlinear optimizers often report false convergence, and users should not simply trust whatever convergence criteria an optimizer uses. \citet{McCulloughVinod2003} offer four criteria for verifying the solution of a nonlinear solver. These criteria are meaningful only for problems that meet regularity conditions at the solution, notably differentiability, and as such are not useful for the LQD example offered above. The four criteria are: (1) making sure the gradients are zero; (2) inspecting the solution path (i.e., the trace) to make sure it follows the expected rate of convergence; (3) evaluating the Hessian to make sure it is well-conditioned;\footnote{Following an exchange with \citet{DrukkerWiggins2004}, \citet{McCulloughVinod2004b} modify their third suggestion to note that determining if the Hessian is well-conditioned in part depends on how the data are scaled. That is, a Hessian that appears to be ill-conditioned may be made well-conditioned by rescaling. So if an Hessian appears to be ill-conditioned, \citet{McCulloughVinod2004b} recommend that the analyst attempt to determine if rescaling the data can result in a well-conditioned Hessian.} and (4) profiling the likelihood to ascertain if a quadratic approximation is adequate. One may need to take the results from \code{genoud} and pass them to \code{optim} to conduct some of these diagnostic tests such as to profile the likelihood. It is also good practice to use more than one optimizer to verify the solution \citep{stokes2004}. Note that \code{genoud} uses its own random number seeds and internal pseudorandom number generators to insure backward compatibility with the original C version of the software and to make cluster behavior more predictable. These seeds are set by the \code{unif.seed} and \code{int.seed} options. The \proglang{R} \code{set.seed} command is ignored by \code{genoud}. %user supplied gradients and lexical search? \section{Conclusion} \label{sec:conclusion} The \code{genoud} function provides many more options than can be reviewed in this brief paper. These options are described in the \proglang{R} help file. The most important option influencing how well the evolutionary algorithm works is the \code{pop.size} argument. This argument controls the population size---i.e., it is the number of individuals \code{genoud} uses to solve the optimization problem. As noted above, the theorems proving that genetic algorithms find good solutions are asymptotic in both population size and the number of generations. It is therefore important that the \code{pop.size} value not be small. On the other hand, computational time is finite so obvious trade-offs must be made. As the LQD example illustrates, a larger population size is not necessarily demonstrably better. The most important options to ensure that a good solution is found, aside from \code{pop.size}, are \code{wait.generations}, \code{max.generations} and \code{hard.generation.limit}. Many statistical models have objective functions that are nonlinear functions of the parameters, and optimizing such functions is tricky business \citep{AltmanGillMcDonald2003}. Optimization difficulties often arise even for problems that are generally considered to be simple. For a cautionary tale on how optimization can be a difficult task even for such problems see Stokes' (\citeyear{stokes2004}) effort to replicate a probit model estimated by \citet[pp. 335]{maddala1992}. A recent controversy over estimating a nonlinear model estimated by maximum likelihood offers another cautionary tale \citep{DrukkerWiggins2004,McCulloughVinod2003,McCulloughVinod2004a,McCulloughVinod2004b,SchacharNalebuff2004}. End users are generally surprised to learn that such optimization issues can arise, and that results can substantially vary across optimizers and software implementations. The \pkg{rgenoud} package provides a powerful and flexible global optimizer. When compared with traditional derivative-based optimization methods, \pkg{rgenoud} performs well \citep{SekhonMebane1998}. Optimization of irregular functions is, however, as much of an art as science. And an optimizer cannot be used without thought even for simple surfaces, let alone spaces that require a genetic algorithm. We hope that the availability of a scalable global optimizer will allow users to work with difficult functions more effectively. \bibliography{rgenoud} \end{document} % LocalWords: optima GillMurrayWright EfronTibshirani GENOUD MebaneSekhon LQD % LocalWords: SekhonMebane multinomRob overdispersed nondifferentiable Broyden % LocalWords: Goldfarb Shanno BFGS optim genoud EA EAs GA holland goldberg GAs % LocalWords: GrefenstetteBaker davis FilhoTreleavenAlippi generation's dejong % LocalWords: substring's grefenstette billingsley Vose NixVose vose GA's ADC % LocalWords: GENOCOP MichalewiczSwaminathanLogan Michalewicz Polytope dnorm % LocalWords: nvars Nonchanging UniqueCount var MemoryMatrix IQR LQDxmpl exp % LocalWords: sIQR yhat na rm cbind rnorm rbinom Min Std Pr Signif suby subx % LocalWords: apriori Raphson saddlepoints fn Genoud's gr hessian unif int ssh % LocalWords: Tausworthe stdout tempdir polytope musil deckard makeSOCKcluster % LocalWords: usernames username PVM MPI peakgeneration parm UC Outlier seq mu % LocalWords: lapply xlab ylab col biclaw mNd ALTMAN rgenoud WandKingLau PKfit % LocalWords: KingWand braumoeller BraumoellerGoodrichKline pharmacokinetics % LocalWords: LeeLee pseudorandom pdf glm rw wrm lqdxmpl cl makeCluster probit % LocalWords: localhost clusterExport stopCluster AltmanGillMcDonald GENetic % LocalWords: McCulloughVinod SchacharNalebuff GenMatch GUI URL boolean Sekhon % LocalWords: CPUs nonlinearity FAiR goodrich SekhonMatching ivivc vitro vivo % LocalWords: multinomial Synth AbadieGardeazabal synth nonsmooth Rcore maxima % LocalWords: resampling resamples optimizer's eval nonstochastic testbounds % LocalWords: LLLL func ij unimodal Rosenbrock multimodal minima saddlepoint % LocalWords: testfuncs testNparms gradcheck sizeset genset nreps gsarray tol % LocalWords: gsize ngens FEP LLLLLL sqrt dLQDxmpl eps dsIQR beps blen lenbseq % LocalWords: bseq IQRarr IQRarrA dimnames mfrow lwd loc online DrukkerWiggins % LocalWords: maddala funcs sapply runif ifelse rbind crossprod byrow vXXi % LocalWords: interquartile