I have been inspired by Varadi’s latest series of posts in http://cssanalytics.wordpress.com regarding Heuristic optimization techniques (particularly Genetic Algorithms, GAs, and Particle Swarm Optimization, PSO) to perform a quick analysis comparing said techniques. I highly recommend you read the http://cssanalytics.wordpress.com blog posts from Sept. 3-6, 2013, for a detailed explanation of GAs and PSO before viewing my analysis.
A Quick Analysis of Heuristic Optimization by Stochastic Genetic Algorithms and Particle Swarm
We investigate the optimization performance of both Genetic Algorithms (GAs) and Particle Swarm Optimization (PSO). Particularly we are interested in comparing the efficiency and accuracy of these techniques in optimizing a variety of well known test functions. The test functions used are outlined below.
We consider three classes of GAs:
- Matlab’s Built in Traditional GA (MatlabGA).
- A traditional GA as outlined by Varadi in http://cssanalytics.wordpress.com/ (GA)
- A simplistic Micro-GA (mGA) as outlined by Varadi in http://cssanalytics.wordpress.com/.
A Brief Review of Traditional GAs
Recall, like most evolutionary algorithms, a GA is a stochastic sampling method that repeatedly refines a population or set of candidate solutions. The method is based on the Darwinian notion of “Survival of the Fittest” in which only the elite solutions, (those that return the most optimal objective function values) are kept, while sub-optimal solutions are discarded.
A typical GA begins by randomly generating an initial population candidate solutions, called the parents. The initial population of each generation is then enlarged through production of a combination of cross-over children and/or mutation children. Cross-over children are derived by mixing and matching individual components of the parents, in what Varadi describes as Mating. Formally, at each component 
Next we evaluate the corresponding objective function value at each of the potential solution in the current population. Only a fraction of the individuals in the current population that provide the most optimal objective function value will continue onward to form the parent population of the next generation.
The mGA, on the other hand, is a simplified version of the traditional GA and involves a Tournament phase for determining which parent solutions are worthy to mate. Those parents that are worthy to mate produce cross-over children, as described above and in http://cssanalytics.wordpress.com/.
For consistency amongst the varying GAs we have kept each population size and maximum number of generations fixed at 100 and 100d, respectively, where d is the number of input dimensions to the test function. Convergence is established either by achieving the maximum number of generations or if the average of the relative difference between a generation’s solution and best solution found over the most recent 50 generations is less than a tolerance, 
A Brief Review of PSO
PSO is a gradient free, stochastic optimization algorithm. Like GAs and mGAs, PSO begins by first generating an initial set of candidate solutions or particles, collectively called the swarm. At each iteration, individual particles of the swarm navigate the search-space based on each particles personal best position and the overall best position found by the swarm. Therefore, previous locations that minimize the objective function, 
where 
For optimization of each test function I have chosen to use a swarm size of 20. For a fair comparison with the GAs, the stopping criteria for PSO is that the average relative difference of the best solution found over the most recent 50 iterations is less than a tolerance,
Lastly, note that since both the GA and PSO are stochastic optimization routines, the results presented below have each been averaged over 100 simulations.
Numerical Results
We analyze three quantities:
- Accuracy: Average Percent Error =
, where
is the true, known, global optimum value.
- Efficiency: Average number of Function Evaluations (FEs) required for convergence.
- General Performance Assessment:} GPA= FEs
in which accuracy holds twice the weight as efficiency and smaller GPA values indicate optimal performance.
Simulation results for each test function can be found in Table 1. From Table 1 we can draw the following conclusions:
- In general MatlabGA outperforms both my own traditional GA and mGA, with the exception of optimizing the 2-D Hump Function and 6-D Hartmann function. This is not surprising given the level of sophistication of Matlab’s built in optimization algorithms as compared to the 1-day coding blitz of my GAs.
- The mGA is indeed the most efficient of the GAs, as measured by the number of FEs required to establish convergence. Of course improved efficiency comes at the cost of decreased accuracy – ultimately the user’s optimization goals will play a huge factor in algorithm selection.
- As expected, PSO outperforms all GAs for objective functions that are strictly convex, with minimal noise and a single global optimum.
- As Varadi suggests, PSO (or a hybrid PSO method) may be the preferred technique for several portfolio optimization problems.
- Objective surfaces such as the Schwefel Function and the Rastrigin Function are noisy and multi-modal, containing several local optima. For global optimization of both these functions, MatlabGA proves to outperform PSO.
For full test function details see:
http://www-optima.amp.i.kyoto-u.ac.jp/member/student/hedar/Hedar_files/TestGO.htm
. Therefore, we will use DIRECT to find a sub-region within the feasible region that contains the global optimum, while simultaneously reducing the chance of searching in regions that contain local optima. 
. DIRECT then divides the hypercube into smaller hyper-rectangles by evaluating the objective function at one-third of the distance, 
, from the center in all coordinate directions ,
be a positive constant and let 
be the current best function value. A hyper-rectangle j is said to be potentially optimal if there exists 
such that}
is the center of hyper-rectangle j and 
is the distance from the center to the corresponding vertices. In general, a value of 
is used to ensure that 
exceeds 
for all hyper-rectangles of equal size, 
.
for all other hyper-rectangles k, and 
for all other hyper-rectangles k, and 
for all longest dimensions i. Of course, if the hyper-rectangle is a hypercube, then the objective function is evaluated in all coordinate directions 
, which is identical to the initial sampling step. We then divide the hyper-rectangle, with order determined by the objective function values returned during the sampling phase.
where 
is a weighting parameter specified by the used, 
are random numbers on the interval 
denotes the best previous position of that particle, and 
denotes the best previous position of the swarm. Qualitatively, the velocity component of PSO is specified by three parameters:
. Unless otherwise specified, the first iteration searches with a mesh size of 1 in all n directions, where n is the dimension of the objective function. The objective function is then computed at each of the mesh points, until it finds a value that is less than the value of the current objective function. Once it has done so, the algorithm sets to that point and begins polling again. After a successful poll, the mesh size is multiplied by a factor of 2. If a poll is unsuccessful, the algorithm remains at it’s current position, and polls once again, this time with a mesh size that is 0.5 the former mesh size. GPS can handle both bound and general constraints optimization problems. Humphries’ hybridization technique involves implementing PSO until we encounter an iteration that does not improve our solution. From there we perform one step of GPS, to poll for a better solution. If no better solution is found, we reduce the stencil size and poll again.
. The linear approximation, as measured by the tangent to the curve at the point 
such that 
. Therefore, assuming the objective function is twice differentiable, Newton’s Method uses 
in order calculate the roots of 
, and the iterative process becomes
and 
with the gradient, 
and the Hessian, 
, respectively. The main advantages of Newton’s method is its robustness and rapid local quadratic rate of convergence. Unfortunately, this method requires computation and storage of the Hessian matrix 
with a finite differencing scheme
The BFGS quasi-Newton approach can therefore be thought of as a generalization of the secant method. 
by approximating the Hessian matrix 
at $x_k$. We can then obtain a search direction 
Like Newton’s Method, we follow along this search direction to obtain a new guess
where 
is the step-size, often determined through a scalar optimization or a line search procedure.
, and is that of “most” descent. This simplification, however, decreases the rate of convergence such that the method of steepest descent converges linearly.
and 
for search direction 
.
using the BFGS update.
under the medium-scale implementation. Therefore, the first iteration of fminunc follows the method of steepest descent. We define 
, where 
is the magnitude of the gradient of the likelihood function, 
. The gradient function, 
, is derived using a forward finite differencing scheme. If 
. 
is a scalar and
, and therefore simplifies to
, along the line of steepest descent in order to obtain a second point 
. The two points along with the initial gradient approximation are used to obtain a quadratic interpolation. A third point, 
, is obtained through extrapolation of the quadratic interpolant, and brackets the minimum such that 
. Using each of 
, for 
and the first gradient approximation, allows us to obtain a cubic interpolant, whose minimum, 
, estimates the true minimum. We then proceed to bisect the bracket interval using two of 
along with the minimum estimate 
, will sample the objective function in all coordinate directions via 
. Here 
, where 
and 
represent the components of the lower and upper bounds of the feasible region, respectively, and 
. Therefore, if a sampling phase is unsuccessful in finding a more optimal position than the current incumbent point, the scale, and thus the stencil, is reduced by a factor of 0.5, and the sampling phase is repeated. By default, the optimization will terminate after the scale has been reduced by a factor of 
. It is, however, the quasi-Newton iteration that provides IF with a clear advantage over a general pattern search algorithm, and has shown to improve overall efficiency and ability of global convergence. 
, determined at the end of each sampling step. Therefore, IF must be given two main stopping criterion; one for each loop.
. Here, 
, where 
is a scaling factor that attempts to eliminate discrepancies in the step size of the quasi-Newton iteration, which arise when the values of the objective function are either very small or very large. The default scale value is 
, where 
is the tolerance level for terminating the quasi-Newton iteration. Highly sensitive stopping tolerances are required for optimization of 
, and the corresponding simulator outputs, 
, where 
, n is the number of user supplied design points, and d is the dimension of the problem. We denote 
as the set of n design points, and 
as the vector that holds the corresponding simulator output. The simulator output is thus modelled as
where 
is the overall mean, and 
is a collection of random variables, which have a joint normal distribution, and thus is a Gaussian Process. For example, in such a case where the overall mean, 
, then the random variables, ![E[z(x_i)]=0](https://s0.wp.com/latex.php?latex=E%5Bz%28x_i%29%5D%3D0&bg=efefef&fg=555555&s=0&c=20201002)
, variance ![\text{Var}[z(x_i)]=\sigma^2](https://s0.wp.com/latex.php?latex=%5Ctext%7BVar%7D%5Bz%28x_i%29%5D%3D%5Csigma%5E2&bg=efefef&fg=555555&s=0&c=20201002)
, and covariance function ![\text{Cov}[(z(x_i),z(x_j)]=\sigma^2 R](https://s0.wp.com/latex.php?latex=%5Ctext%7BCov%7D%5B%28z%28x_i%29%2Cz%28x_j%29%5D%3D%5Csigma%5E2+R&bg=efefef&fg=555555&s=0&c=20201002)
. 
is an 
vector of scaling parameters. Determining the optimal 
scaling parameter is crucial for ensuring the integrity of the GP model fit. Originally, Dr. Ranjan proposed the use of the scaling parameter 
in place of 
, which confines the parameter space to 
. However, they noticed that optimal 
, making the optimization process challenging. We therefore proposed a re-parameterization of the 
, such that 
\cite{P2}. Such a re-parameterization is preferred as it smoothes the parameter space, such that it is centered about 
, and ensures that potential 
to the spatial distribution 
for all 
and 
, given by
is an 
vector of all ones. Calculation of the mean and variance estimators, along with the repeated update of the inverse and determinant of the correlation matrix, 
, and 
respectively, allow us to derive the profile log-likelihood![\mathcal{L}_\beta = -2\log{(L_\beta)} \ \alpha \ \log{(|R|)} + n\log{[(Y-1_n\hat{\mu}(\beta))'R^{-1}(Y-1_n \hat{\mu}(\beta))]}.](https://s0.wp.com/latex.php?latex=%5Cmathcal%7BL%7D_%5Cbeta+%3D+-2%5Clog%7B%28L_%5Cbeta%29%7D+%5C++%5Calpha+%5C+%5Clog%7B%28%7CR%7C%29%7D+%2B+n%5Clog%7B%5B%28Y-1_n%5Chat%7B%5Cmu%7D%28%5Cbeta%29%29%27R%5E%7B-1%7D%28Y-1_n+%5Chat%7B%5Cmu%7D%28%5Cbeta%29%29%5D%7D.+&bg=efefef&fg=555555&s=0&c=20201002)
, for a desired predictor point ![\hat{y}(x^*) = \hat{\mu}+r'R^{-1}(Y-1_n\hat{\mu}) = \left[\frac{(1_n-r'R^{-1}1_n)}{{1_n}'R^{-1}1_n} {1_n}' + r'\right]R^{-1}Y = C'Y,](https://s0.wp.com/latex.php?latex=%5Chat%7By%7D%28x%5E%2A%29+%3D+%5Chat%7B%5Cmu%7D%2Br%27R%5E%7B-1%7D%28Y-1_n%5Chat%7B%5Cmu%7D%29+%3D+%5Cleft%5B%5Cfrac%7B%281_n-r%27R%5E%7B-1%7D1_n%29%7D%7B%7B1_n%7D%27R%5E%7B-1%7D1_n%7D+%7B1_n%7D%27+%2B+r%27%5Cright%5DR%5E%7B-1%7DY+%3D+C%27Y%2C&bg=efefef&fg=555555&s=0&c=20201002)
![r=[r_1(x^*),...,r_n(x^*)]](https://s0.wp.com/latex.php?latex=r%3D%5Br_1%28x%5E%2A%29%2C...%2Cr_n%28x%5E%2A%29%5D&bg=efefef&fg=555555&s=0&c=20201002)
, and ![r_i(x^*)=\text{corr}[z(x^*),z(x_i)],](https://s0.wp.com/latex.php?latex=r_i%28x%5E%2A%29%3D%5Ctext%7Bcorr%7D%5Bz%28x%5E%2A%29%2Cz%28x_i%29%5D%2C&bg=efefef&fg=555555&s=0&c=20201002)
the correlation between 
and 
, as measured by the mean squared error (MSE),![s^2(x^*)=E \left[ ( \hat{y}(x^*)- {y}(x^*))^2\right] = \hat{\sigma}^2 (1_n-2C'r+C'RC).](https://s0.wp.com/latex.php?latex=s%5E2%28x%5E%2A%29%3DE+%5Cleft%5B+%28+%5Chat%7By%7D%28x%5E%2A%29-+%7By%7D%28x%5E%2A%29%29%5E2%5Cright%5D+%3D+%5Chat%7B%5Csigma%7D%5E2+%281_n-2C%27r%2BC%27RC%29.+&bg=efefef&fg=555555&s=0&c=20201002)
. The addition of the nugget parameter smoothes the model predictor and therefore the GP does not produce an exact interpolant of our design points. To avoid unnecessary over-smoothing of the model predictions a lower bound on the nugget parameter is introduced, and is given by
, where 
are the eigenvalues of R, and 
is the desired threshold for the $
. We define an nXn matrix R to be near-singular if its condition number 
exceeds the threshold value 
denotes the matrix L2-norm. A value of a=25 is recommended for a space-filling Latin hypercube designs (LHDs) that is used in determining the design points. In terms of optimization, it is important to note that 
is a function of both the design points and the optimization parameter 
is dependent upon several components, namely the mean and variance estimators 
and $\hat{\sigma}^2(\beta)$, as well as the inverse and determinant of 
, all of which are dependent upon the hyper-parameter 
is a 
vector of parameters, with individual components 
, for all $i=\{1,…d\}$. Additionally, the nugget parameter 
are of no benefit. We can, however, provide numerical approximations to 
Butler's Math 