Genetic Optimization Using Derivatives: The rgenoud Package for R

This introduction to the R package rgenoud is a modiﬁed version of Mebane and Sekhon (2011), published in the Journal of Statistical Software . That version of the introduction contains higher resolution ﬁgures. genoud is an R function that combines evolutionary algorithm methods with a derivative-based (quasi-Newton) method to solve diﬃcult optimization problems. genoud may also be used for optimization problems for which derivatives do not exist. 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 ﬁnd 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 need-lessly 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.


Introduction
We developed the R package rgenoud to solve difficult optimization problems such as often arise when estimating nonlinear statistical models or solving complicated nonlinear, nonsmooth and even discontinuous functions. 1Optimization 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 (Gill, Murray, and Wright 1981).Even in models where such methods can be expected to work most of the time, resampling techniques such as the bootstrap (Efron and Tibshirani 1994) 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.
We present in Section 3 an example using benchmark functions taken from Yao, Liu, and Lin (1999), followed by an example motivated by the multinomRob package.The benchmark suite includes functions that are high-dimensional, discontinuous or that have many local optima.The multinomRob package robustly estimates overdispersed multinomial regression models and uses rgenoud to solve a least quartile difference (LQD) generalized S-estimator (Mebane and Sekhon 2004).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 Matching package and its use of rgenoud are described in detail (Sekhon 2011a).Matching provides functions for multivariate and propensity score matching and for finding optimal covariate balance based on a genetic search algorithm implemented in rgenoud.The search over optimal matches is discontinuous so no derivatives are used. 2The search also involves lexical optimization which is a unique feature implemented in rgenoud. 3he rgenoud package implements an updated and extended version of the C program GENOUD (Mebane and Sekhon 1997) described in (Sekhon and Mebane 1998).The many improvements include among other things the interface with R, which includes the ability to optimize functions written in 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 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 rgenoud program combines an evolutionary algorithm with a quasi-Newton method.The quasi-Newton method is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (Gill et al. 1981, 119) implemented in R's optim function.When the BFGS is being used, our program offers the option of using either rgenoud's built-in numerical derivatives (which are based on code taken from Gill et al. 1981, 337-344) or user-supplied analytical derivatives. 4Our 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.5As 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" (Gill et al. 1981, 285).In this article, we very briefly review the theory of random search algorithms that supports the assertion that rgenoud has a high probability of finding the global optimum when such exists.And we present three examples of how to use the 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 genoud can be found in our article that describes GENOUD (Sekhon and Mebane 1998).

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).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 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 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 (Holland 1975;Goldberg 1989;Grefenstette and Baker 1989).The operators used in GA implementations vary (Davis 1991;Filho and Alippi 1994), 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 (Holland 1975, 139-140).Each generation's population contains a biased sample of codestrings, so that a substring's performance in that population is a biased estimate of its average performance over all possible populations (De Jong 1993;Grefenstette 1993).
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.6An irreducible, finite Markov chain converges at an exponential rate to a unique stationary distribution (Billingsley 1986, 128).This means that the probability that each population occurs rapidly converges to a constant, positive value.Nix and Vose (1992) and Vose (1993) 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.
The evolutionary algorithm in rgenoud uses nine operators that are listed in Table 1.The operators extend and modify a set of operators used in GENOCOP (Michalewicz, Swaminathan, and Logan 1993;Michalewicz 1992).The operators are numbered using syntax matching that used to refer to them by rgenoud.The cloning operator simply makes copies of the best trial solution in the current generation (independent of this operator, rgenoud always retains one instance of the best trial solution).The uniform mutation, boundary mutation and 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.Whole non-uniform mutation does non-uniform mutation for all the parameters in the trial solution.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.Polytope crossover (inspired by simplex search, Gill et al. 1981, 94-95) computes a trial solution that is a convex combination of as many trial solutions as there are parameters.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.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.

Examples
The only function in the rgenoud package is genoud.The interface of this function is similar to that of the optim function in R.But the function has many additional arguments that control the behavior of the evolutionary algorithm.

Asymmetric Double Claw:
Our first example, which we also studied in Sekhon and Mebane (1998), is a normal mixture called the Asymmetric Double Claw (ADC).We plot the function in Figure 1.Mathematically, this mixture is defined as 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 1.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 P1 Cloning.Copy X t into the next generation, X t+1 .
probability 1/2 of using each value.
P4 Non-uniform Mutation.At random choose i ∈ 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 ∼ U (0, 1).Set either with probability 1/2 of using each value.
P5 Polytope Crossover.Using m = max(2, n) vectors x from the current population and m random numbers p j ∈ (0, 1) such that m j=1 p j = 1, set X = m j=1 p j x j .
P6 Simple Crossover.Choose a single integer i from N. Using two parameter vectors, x and y, set P7 Whole Non-uniform Mutation.Do non-uniform mutation for all the elements of X.
P8 Heuristic Crossover.Choose p ∼ U (0, 1).Using two parameter vectors, x and y, compute z = p(x − y) + x.If z satisfies all constraints, use it.Otherwise choose another p value and repeat.Set z equal to the better of x and y if a satisfactory mixed z is not found by a preset number of attempts.In this fashion produce two z vectors.
P9 Local-minimum Crossover.Choose p ∼ U (0, 1).Starting with x, run BFGS optimization up to a preset number of iterations to produce x.
If z satisfies boundary constraints, use it.Otherwise shrink p by setting p = p/2 and recompute z.If a satisfactory z is not found by a preset number of attempts, return x.This operators is extremely computationally intensive, use sparingly. Notation: x i is the lower bound and x i is the upper bound on values for X i .x i is the current value of X i , and x is the current value of X. N = {1, . . ., n}. p ∼ U (0, 1) means that p is drawn from the uniform distribution on the [0, 1] interval.
The fourth option pop.size controls the most important part of the evolutionary algorithm, the population size.This is the number of individuals 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 genoud even using the default value of pop.size=1000.Reliability does increase as the 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 max.generations,wait.generations and hard.generation.limit.
The specified termination point also affects how some of the operators perform: the two nonuniform mutation operators introduce smaller ranges of variation in the trial solutions as the generation count approaches the specified max.generations value.There are many more options that can be used to fine-tune the behavior of the algorithm.
The output printed by genoud is controlled by a print.levelargument.The default value, 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 R conventions would suggest setting the default to be print.level=0,which would suppress output to the screen, but because 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 genoud, which uses the default value for a print.levelargument, is as follows.
Fri 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 data.type.intargument.
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 genoud quickly converge once the best trial solution is in the correct neighborhood.This run of 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 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 genoud is enforcing the boundary constraints specified by the domain values.By default (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 boundary.enforcementargument 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, 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 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 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 genoud.

Benchmark Functions:
The second example is a suite of 23 benchmark nonlinear functions used in Yao et al. (1999) to study a pair of evolutionary programming optimization algorithms.Function definitions are in Table 2.Because it includes a random component and so lacks a reproducible minimum value, we ignore the function numbered function 7 in their sequence. 8Implementations of these functions are available in the supplemental R file provided with this article.9These 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 Yao et al. (1999) via bounds that are stated in the list named testbounds in the supplemental file.
As Yao et al. (1999) 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 ∈ [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.Yao et al. (1999, 84) 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 Yao et al. (1999) have trouble optimizing functions 21-23.Although Yao et al. (1999, 85) 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 (Yao The omitted function is n i=1 ix 4 i + 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 n i=1 ix 4 i , i.e., xi = 0, the value of the function virtually never attains the minimum possible value of zero.An optimizer that evaluated the function at xi = 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.
We use these benchmark functions to illustrate not only how effective 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 genoud two arguments determine the number of generations.One is max.generations: if hard.generation.limit=TRUE,then the value specified for max.generations is a binding upper limit.The other is wait.generations,which determines when the algorithm terminates if hard.generation.limit=FALSE.But even if hard.generation.limit=TRUE,then wait.generationsdetermines 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 wait.generationsvalue never worsens the algorithm's efficacy.
Increasing max.generations may or may not in itself improve optimization.The value of max.generations sets the value of T used in the mutation operators-operators 4 and 7 in Table 1.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 max.generations without changing wait.generationsincreases 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 genoud to minimize the 23 functions using two values for pop.size (5000 and 10000) and two values for max.generations (30 and 100).Following Yao et al. (1999), we replicate each optimization 50 times.The following code describes the computations.The list testfuncs, vector testNparms and list testbounds are defined in the supplemental R file, and it is assume that this file is loaded with the source("supplemental.R") command.The vector gradcheck is true for all elements except the ones corresponding to functions 6 and 7.
> source("supplemental.R") > gradcheck <-rep( Using genoud to minimize the benchmark functions produces excellent results, at least when the pop.size and max.generations arguments are sufficiently large.Table 3 reports the mean function values for each configuration of the arguments.These values may be compared both to the true function minima given by Yao et al. (1999) (see the rightmost column in Table 2) and to the average minimum function values Yao et al. (1999) report for their "fast" evolutionary programming (FEP) algorithm, which appear in the last column of Table 3.The genoud averages for the 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 pop.size=5000 case are off.The effect of increasing pop.size is apparent with respect to both those three functions and also functions 13 and 20-23: the average minima are smaller with pop.size=10000 than with pop.size=5000.Except for functions 6 and 12 in the pop.size=5000 case and function 13, all the genoud solution averages for max.gen=100 are either slightly or substantially better than the corresponding FEP solution averages.
The results in Table 3 clearly illustrate the potential consequences of not allowing genoud to run for a sufficient number of generations.While some of the genoud solutions for max.gen=30 have competitive means, several of the means are not good at all.The effect of increasing pop.size are even more clearly apparent in Table 4, 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 pop.size=5000, the genoud solutions for max.gen=100 vary much less than the corresponding FEP solutions.For those functions and also for functions 20, 22 and 23, the max.gen=100 solutions with pop.size=10000 vary noticeably less than the solutions with pop.size=5000.
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: Call: glm(formula = cbind(y, m -y) ~x[, 2] + x[, 3], family = "binomial") But in practical applications it is unknown apriori which observations should be treated as outliers.
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] × [−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 Journal of Statistical Software version of this documentation).11Within this smaller domain the sign of the derivative changes even 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 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 solution.toleranceargument.We retain the default setting 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 (Mebane and Sekhon 2004) 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 LQDxmpl function matches the behavior of 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 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.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.
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.McCullough and Vinod (2003) 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; 12 and (4) profiling the likelihood to ascertain if a quadratic approximation is adequate.One may need to take the results from genoud and pass them to 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 (Stokes 2004).
Note that 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 unif.seedand int.seed options.The R set.seed command is ignored by genoud.

Conclusion
12 Following an exchange with Drukker andWiggins (2004), McCullough andVinod (2004a) 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, McCullough and Vinod (2004a) recommend that the analyst attempt to determine if rescaling the data can result in a well-conditioned Hessian.
The genoud function provides many more options than can be reviewed in this brief paper.These options are described in the R help file.The most important option influencing how well the evolutionary algorithm works is the pop.size argument.This argument controls the population size-i.e., it is the number of individuals 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 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 pop.size, are wait.generations,max.generations and hard.generation.limit.Many statistical models have objective functions that are nonlinear functions of the parameters, and optimizing such functions is tricky business (Altman, Gill, and McDonald 2003).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' (2004) effort to replicate a probit model estimated by Maddala (1992, pp. 335).A recent controversy over estimating a nonlinear model estimated by maximum likelihood offers another cautionary tale (Drukker and Wiggins 2004;McCullough andVinod 2003, 2004b,a;Shachar and Nalebuff 2004).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 rgenoud package provides a powerful and flexible global optimizer.When compared with traditional derivative-based optimization methods, rgenoud performs well (Sekhon and Mebane 1998).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.

Figure 2 :
Figure 2: Contour Plots of the LQD Fitness Function and of its Partial Derivatives 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.

Table 3 :
Mean Values of 22 Optimized Functions rgenoud: Genetic Optimization Using Derivatives in R Standard deviation of the minimum function values (over 50 replications) obtained using genoud.bStandarddeviation of the best function values (over 50 replications) reported for the "fast" evolutionary programming algorithm, fromYao et al. (1999, 85 and  88, TablesII-IV).
Of course, if we knew which observations to omit the results would be much better: 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: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):