A global optimization paradigm based on change of measures.

A global optimization framework, COMBEO (Change Of Measure Based Evolutionary Optimization), is proposed. An important aspect in the development is a set of derivative-free additive directional terms, obtainable through a change of measures en route to the imposition of any stipulated conditions aimed at driving the realized design variables (particles) to the global optimum. The generalized setting offered by the new approach also enables several basic ideas, used with other global search methods such as the particle swarm or the differential evolution, to be rationally incorporated in the proposed set-up via a change of measures. The global search may be further aided by imparting to the directional update terms additional layers of random perturbations such as 'scrambling' and 'selection'. Depending on the precise choice of the optimality conditions and the extent of random perturbation, the search can be readily rendered either greedy or more exploratory. As numerically demonstrated, the new proposal appears to provide for a more rational, more accurate and, in some cases, a faster alternative to many available evolutionary optimization schemes.

Introduction

Potential applications of global optimization are ubiquitous across a large spectrum of problems in science and engineering, from very-large-scale integration circuit design, to optimization of transportation routes to determination of protein structures. An essential aspect of the solution to any such problem is the extremization of an objective functional, possibly subject to a prescribed set of constraints. The solution is deemed to have been achieved if the global extremum of the objective functional is reached. For many practical applications, the objective functional could be non-convex, non-separable and even non-smooth. This precludes gradient-based local optimization approaches to treat these problems, an aspect that has driven extensive research in devising global optimization strategies, many of which employ stochastic (i.e. random evolutionary) search of heuristic or meta-heuristic origin [1]. Indeed, finding the global optimum within a search space of the parameters or design variables is much more challenging than arriving at a local extremum. In the absence of generic directional information (equivalent to the Gateaux or directional derivative in the search for a local extremum of a sufficiently smooth cost functional), a direct search seems to be the only option left. Some of the notable schemes in this genre include variants of the genetic algorithm (GA) [2], differential evolution (DE) [3], particle swarm optimization (PSO) [4], ant colony optimization [5] and covariance matrix adaptation evolution strategy (CMA-ES) [6]. In a search for the global extremum, stochastic schemes typically score over their gradient-based deterministic counterparts [7,8], even for cases involving sufficiently smooth cost functionals with well-defined directional derivatives. Nevertheless, thanks to the proper directional information guiding the update, gradient-based methods possess the benefit of faster convergence to the nearest extremum, as long as the objective functional is sufficiently smooth. To the authors' knowledge, the existing evolutionary global optimization schemes do not offer the benefit of well-founded directional information. Most evolutionary schemes depend on a random (diffusion-type) scatter applied to the available candidate solutions or particles and some criteria for selecting the new particles. Unfortunately, a ‘greedy’ approach for the global optimum (e.g. selecting particles of higher ‘fitness’ with higher probability), though tempting from the perspective of faster convergence, faces the hazard of getting trapped in local extrema. As a fix to a premature and possibly wrong convergence, many evolutionarily global search schemes adopt safeguards. Despite the wide acceptability of a few evolutionary methods of the heuristic/meta-heuristic origin [9], the underlying justification is often based on sociological or biological metaphors [1-6,10] that are hardly founded on a sound probabilistic basis even though a random search forms a key ingredient of the algorithm. In this context, some interesting comments on the facile usage of natural metaphors in the development of many global optimization algorithms, without a serious engagement with scientific ratiocination, may be found in [10]. Be that as it may, the popular adoption of these schemes is not only owing to the algorithmic simplicity, but mainly because of their effectiveness in treating many non-deterministic polynomial-time hard optimization problems, which may be contrasted with a relatively poorer performance of some of the well-grounded stochastic schemes, e.g. simulated annealing [11], stochastic tunnelling [12], etc. Whether a slow convergence to the global optimum, not infrequently encountered with some of the popular evolutionary schemes, could be fixed by reorienting them with an alternative stochastic framework, however, remains a moot point, which is addressed in this work to an extent. Another related question is regarding the precise number of parameters in the algorithm that the end-user must tune for an appropriate exploration–exploitation trade-off [13,14], a phrase used to indicate the relative balance of random scatter of the particles vis-a-vis their selection based on the ‘fitness’. Admittedly, the notion of random evolution, built upon Monte Carlo (MC) sampling of the particles, efficiently explores the search space, though at the cost of possibly slower convergence [15] in contrast to a gradient-based method. For better exploration, many such schemes, e.g. the GA, adopt steps like ‘crossover’, ‘mutation’, ‘selection’, etc. While ‘crossover’ and ‘mutation’ impart variations in the particles, by the ‘selection’ step each particle is assigned a ‘weight’ or fitness value (a measure of importance) determined as the ratio of the realized value of the objective functional to the available maximum of the same. In the selection step, the fitness values, used to update the particles via selection of the best-fit individuals for subsequent crossover and mutation, may be considered functionally analogous to the derivatives used in gradient-based approaches. Analogy may also be drawn between the fitness values and the likelihood ratios commonly employed in nonlinear stochastic filtering, e.g. the family of sequential importance sampling filters [16]. Even though a weight-based approach reduces a measure of misfit between the available best and the rest within a finite ensemble of particles, they bring with them the curse of degeneracy, wherein all particles but one tend to receive zero weights as the iterations progress [17]. This problem, also referred to as ‘particle collapse’ in the stochastic filtering parlance, can be arrested, in principle, by exponentially increasing the ensemble size (number of particles), a requirement that can hardly be met in practice [18]. Evolutionary schemes that replace the weight-based multiplicative approach by adopting additive updates for the particles, often succeed in eliminating particle degeneracy. Evolutionary schemes like DE, PSO, etc., that are known to perform better than GA, use such additive particle update strategies without adversely affecting the ‘craziness’ or randomness in the evolution. On reaching the optimum, the additive correction shrinks to zero. Unfortunately, none of these methods obtain the additive term, which may be looked upon as a stochastic equivalent to the derivative-based directional information, in a rigorous or optimal manner—a fact that is perhaps responsible for a painstakingly large number of functional evaluations in some cases.

In framing a rigorous stochastic basis leading to directional updates of the additive type, several possible schemes are suggested in this work. First, it is shown that the problem of optimization may be generically posed as a martingale problem [19,20], which must however be randomly perturbed to facilitate the global search. For instance, one may perform a greedy search for an extremum of an objective functional by solving a martingale problem, which is in the form of an integro-differential equation herein referred to as the extremal equation. In line with the Freidlin–Wentzell theory [21], the martingale problem may further be perturbed randomly so as to enable the greedy scheme to converge to the global extremum as the strength of the perturbation vanishes asymptotically. The solution through this approach, which is viewed as a stochastic process parametered via the iterations, depends crucially on an error or ‘innovation’ function and is deemed to have been realized when the innovation function, interpreted as a stochastic process, is driven to a zero-mean noise or martingale [19]. The martingale structure of the innovation essentially implies that the mean of the computed cost functional in future iterations remains constant and invariant to zero-mean random perturbations of the argument vector (the design variables). The argument vector should then correspond to an extremum. In the evolutionary random set-up, both the cost functional and its argument vector are treated as stochastic diffusion processes even though the original extremization problem is deterministically posed. Realizing a zero-mean martingale structure for the innovation requires a change of measures that in turn leads to additive gain-type updates on the particles. The gain coefficient matrix, which is a replacement for and generalization over the Frechet derivative of a smooth cost functional, provides an efficacious directional search without requiring the functional to be differentiable. Additionally, in order to facilitate the global search, a general random perturbation strategy (using steps such as ‘scrambling’, ‘relaxation’ and ‘selection’) is adopted to insure against a possible trapping of particles in local extrema.

An important feature of the present approach is the flexibility with which the innovation process may be designed in order to meet a set of conflicting demands en route the detection of the global extremum. The efficiency of the global search basically relies upon the ability of the algorithm to explore the search space while preserving some directionality that helps to quickly resolve the nearest extremum. The development of the proposed set-up, referred to as COMBEO (Change Of Measure Based Evolutionary Optimization), recognizes the near impossibility of a specific optimization scheme performing uniformly well across a large class of problems. In this context, one may refer to variants of the famous ‘no free lunch theorem’, where it is proved that the average performance of any two algorithms across all possible problems is identical [22]. Accommodation of this fact had earlier led to a proposal of an evolutionary scheme that simultaneously ran different optimization methods for a given problem with some communications built among the updates by the different methods [23]. A major related contribution of this work is to bring the basic ideas for global search used in a few well-known optimization schemes under a single unified framework propped up by a sound probabilistic basis. In a way explained later, the ideas (or their possible generalizations) behind some of the existing optimization methods may sometimes be readily incorporated in the present setting by just tweaking the innovation process.

The rest of the paper is organized as follows. In §2, the problem of finding an extremum of an objective functional is posed as a martingale problem represented through an integro-differential equation, whose solution only realizes a local optimum. The integro-differential equation is discretized and weakly solved within an MC setting so as to circumvent its inherent circularity. Section 3 discusses several random perturbation schemes to arrive at the global optimum efficiently without getting stuck in local traps. By combining these tools in different ways, a few pseudo-codes are presented in §4. In §5, the performance of COMBEO is compared with DE, PSO and CMA-ES in the context of several benchmark problems. In the same section, the new optimization approach is also applied for the quantitative recovery of boundary layer density variation in a high-speed flow given the light travel time data from an experiment. Finally, the conclusions of this study are drawn in §6.

Search for an extremum through a martingale problem

Here, it is demonstrated that functional extremization, subject to a generic set of equality constraints, can be posed as a martingale problem upon a proper characterization of the design variables within a probabilistic setting. However, before elaborating on this, a few fundamental features of the new evolutionary approach are listed below:

(i) at a given iterate, the solution is a random variable (taking values in the search space);

(ii) hence the solution process along the iteration axis may be considered a stochastic process and the iterations must be so designed that the optimal solution is asymptotically approached by the mean (i.e. the first moment) of the solution process; and

(iii) since upon convergence the mean should be iteration-invariant, random fluctuations about the converged mean must assume the structure of a zero-mean stochastic process (along the iteration axis), thereby allowing us to identify the fluctuations as a noise process, e.g. a Brownian motion, or, more generally, a zero-mean martingale.

In posing the global optimization problem within a stochastic framework, a complete probability space [19] is considered, within which the solution to the optimization problem must exist as an -measurable random variable. Here Ω, known as the population set (sample space), necessarily includes all possible candidate solutions from which a finite set of randomly chosen realizations or particles is evolved along the iteration axis τ. The introduction of τ as a positive monotonically increasing function on is required to qualify the evolution of the solution as a continuously parametered stochastic process. A necessary aspect of the evolution is a random initial scatter imparted to the particles so as to search the sample space for solutions that extremize the cost functional while satisfying the posed constraints, if any. Since the extremal values, global or otherwise, of the objective functional may not be known a priori, the particles are updated iteratively (i.e. along τ) by conditioning on the evolution history of a so-called extremal cost process based on the available particles. The extremal cost process is defined as a stochastic process such that, for any ξ∈[0,τ], its mean is given by the available best objective functional based on its evaluations across the particles at the iteration step denoted by ξ. Now denoting by the filtration (i.e. statistical information pertaining to the evolution history) containing the increasing family of sub σ-algebras generated by the extremal cost process, the aim is to determine the update by conditioning the guessed or predicted solution on . The guessed solution could be a Brownian motion or a random walk along τ. In case there are equality constraints, also contains history for the best realized constraint until τ. For a multivariate multimodal nonlinear objective functional in , a point x* needs to be found such that . Since the design variable x is evolved in τ as a stochastic process, it is parametrized as x:=x(τ). It may be noted that there may not be any inherent physical dynamics in its evolution, and in such cases x may be evolved as a zero-mean Brownian motion or random walk in Ω (even though, as will be shown later, this step is not absolutely necessary). Since only a finite number of iterations can be performed in practice, τ is discretized as τ0<τ1<⋯<τ so that x is evolved over every τ increment. Thus, for the ith iteration with τ∈(τ,τ], x may be thought of as being governed by the following stochastic differential equation (SDE): Here is a zero-mean vector Brownian process with the covariance matrix , where is the noise intensity matrix with its (j,k)th element denoted as g. The discrete τ-marching map for equation (2.1) has the form: where (⋅) stands for (⋅). Let the extremal cost process, generating the filtration , be denoted as . Within a strictly deterministic setting, a smooth functional f(x) tends to become stationary, in the sense of a vanishing first variation, as x approaches an extremal value x*. In the stochastic setting, a counterpart of this scenario could be that any future conditioning of the process x on , which is the filtration generated by , s<ξ, until a past iteration (i.e. ξ<τ), identically yields the random variable x itself, i.e. , where E represents the expectation operator. Via the Markov structure of x one may equivalently postulate that a necessary condition for the extremization of f:=f(x) is to require that E(x|)=x. This characterization renders the process x a martingale with respect to the extremal cost filtration , i.e. once extremized, any future conditional mean of x on the cost filtration remains iteration invariant. See [24] for an introductory treatise on the theory of martingales. An equivalent way of stating this condition is based on an error or ‘innovation’ process defined as −f. An extremization would require updating x to statistically match the extremal process so that the innovation is driven to a zero mean martingale, e.g. a zero-mean Brownian motion or, more generally, a stochastic integral in Ito's sense (or, in a τ-discrete setting, a zero-mean random walk). This enables one to pose the determination of an extremum as a martingale problem, originally conceptualized by Stroock & Varadhan [20] to provide an improved setting for solutions to SDEs beyond Ito's theory. If the update scheme is effective, one anticipates asymptotically, i.e. as . In other words, and, as the norm of the noise intensity associated with x approaches zero, . However, since strictly zero noise intensity is infeasible within the stochastic set-up, the solution in a deterministic setting may be thought of as a degenerate version (with a Dirac measure) of that obtained through the stochastic scheme.

Note that a ready extension of this approach is possible with multi-objective functions, which would merely require building up the innovation as a vector-valued process. Similarly, the approach may also be adapted for treating equality constraints, wherein zero-valued constraint functions, perturbed by zero-mean noise processes, may be used to construct the innovation processes. The easy adaptability of the current set-up for multiple objective functions or constraints may be viewed as an advantage over most existing evolutionary optimization schemes, where a single cost functional needs to be constructed based on all the constraints, a feature that may possibly engender instability for a class of large dimensional problems.

Driving x to an -martingale, or forcing the innovation process to a zero-mean martingale, may be readily achieved through a change of measures. In evolutionary algorithms, e.g. the GA, the accompanying change of measures requires the assignment of weights or fitness values to the current set of particles and a subsequent selection of those with higher fitness, e.g. by rejection sampling. For a better exploration of the sample space, the GA uses steps like ‘crossover’ and ‘mutation’. These steps are biologically inspired safeguards against an intrinsic limitation of the weight-based update, which tends to diminish all the weights but one to zero, thereby leaving the realized particles nearly identical after a few iterations. This problem is referred to as that of weight collapse or particle impoverishment and precipitates premature convergence to a wrong solution. As a way out of this degeneracy in the realized population set (i.e. the ensemble), the present scheme employs a purely additive update, derived following a Girsanov change of measures P→Q (see [19] for a detailed exposition on the Girsanov transformation), which ensures that the innovation process, originally described under measure P, becomes a zero-mean martingale under the new measure Q. A basic tool in effecting this change of measure is the Radon–Nikodym derivative Λ:=dP/dQ (assuming absolute continuity of Q w.r.t. P and vice versa), a scalar valued random variable also called the likelihood ratio (or the fitness or the weight) that multiplicatively updates the design process x as xΛ. In order to bypass the degeneracy problem noted above, one may obtain an additive update by expanding xΛ using Ito's formula [19]. Within an MC set-up, an immediate consequence of the additive update is that particles with lower weights are never eliminated, but rather corrected to have more fitness by being driven closer to an extremum. Drawing analogy with Taylor's expansion of a smooth function(al) that obtains its first-order term based on Newton's directional derivative, the present version of the additive correction may be interpreted as a non-Newton directional term that drives the innovation process to a zero-mean martingale and a precise form of this term is derived later in this section. In general, the innovation process may be a vector. One such case occurs when, by way of addressing possible numerical instability, a single cost functional is split into several so that each of the correspondent innovations is driven to a zero-mean martingale. This in turn ensures that the innovation corresponding to the original cost functional is also driven to a zero-mean martingale. For a simpler exposition of the basic ideas, the method below is presented using a single cost functional only. The formulation can be trivially extended for vector innovation processes of any finite dimension.

Within a τ-discrete framework with τ∈(τ,τ], x is evolved following equation (2.1). An innovation constraint that may be satisfied for an accelerated (greedy) search in a neighbourhood of an extremum is given by where is the extremal cost process, f the objective functional and a P-Brownian increment representing the diffusive fluctuations. Deriving the subsequent integro-differential equation for the search is best accomplished by recasting equation (2.3) as an SDE. Towards this, an -measureable process := may be constructed to arrive at the following incremental form: Here Δτ:=τ−τ is taken as a ‘small’ increment. Since the τ-axis is fictitious, Δτ is replaced by Δτ=τ−τ, so that equation (2.4) is modified to which is essentially correspondent to the SDE:

Note that the replacement of Δτ by Δτ merely modifies the intensity of the noise process Δη in equation (2.3). The form of SDE (2.6) has the desirable feature that the fictitious diffusion coefficient Δτ is an order ‘smaller’ relative to the drift coefficient f(x). However, since η is not standard Brownian, it is more convenient to rewrite equation (2.6) as where W is standard P-Brownian and ρ a more general form of (scalar-valued) noise intensity that may be made a function of τ. For multi-objective optimization problems (where is a vector stochastic process), ρ would be the intensity matrix and the expressions below are so written that they are valid for both scalar and vector cases. Equation (2.7) may be rewritten as where . While it is possible, and even desirable, to replace the Brownian noise term above by one whose quadratic variation is zero (e.g. a Poisson martingale), such a modification is not central to the basic idea and will be considered in future extensions of the work. SDE (2.7) is assumed to satisfy the standard conditions [19] so as to ensure the existence of a weak solution. The -measurable locally optimal (extremal) solution may now be identified with the conditional mean . Considering a new measure Q under which x from equation (2.1) satisfies equation (2.7), the conditional mean may be expressed via the generalized Bayes' formula as where the expectation E(⋅) is taken with respect to the new measure Q and Λ is the scalar fitness given by

As shown in appendix A (a theorem and its corollary) using Ito's expansions of xΛ and (Λ)−1, the incrementally additive updates to arrive at an extremum must follow from the differential equation: and

Note that the directionality of a search process provided by the second (integral) term on the right-hand side (r.h.s.) of the equation above may also be gauged from the fact the integrand in this term can be interpreted as a Malliavin-type derivative using Clarke–Ocone theorem [25]. But the appearance of the unknown term π(f) on the r.h.s. prevents equations (2.11) or (2.12) to be qualified as an SDE in π(x). Indeed, a necessarily nonlinear dependence of f on x and the consequent non-Gaussianity of π(f) would prevent writing the latter in terms of π(x). This results in the so-called closure problem in solving for π(x).

While a direct solution of equations (2.11) or (2.12) yields an extremum (i.e. a local extremum only) in principle, exact/analytical solutions are infeasible owing to the circularity inherent in the closure problem. This has a parallel in nonlinear stochastic filtering, wherein the Kushner–Stratonovich equation [26] (an equivalent of equation (2.11)) also suffers from a similar circularity problem. Motivated by the MC filters used to solve nonlinear filtering problems [27,28], an MC scheme may be developed for a numerical treatment of equations (2.11) or (2.12). A two-stage strategy, viz. prediction and update, may apparently be considered, even though, as we will soon see, the prediction step could be entirely eliminated in the final scheme. The prediction step, as in most evolutionary optimization schemes, aims at an initial random exploration based on equation (2.1). Consider the ith iteration, i.e. τ∈(τ,τ] and let N denote the ensemble size so that one realizes N predicted particles or MC candidates, that must be updated via equation (2.12). For an MC-based numerical solution to equation (2.12), a sample-averaged form of the equation is first written as is the ensemble-averaged approximation to the conditional mean π(.). A particle-wise representation of equation (2.13) is given by where , , , and . is an N-dimensional row vector with entries 1. Note that the second term on the r.h.s. of equation (2.14) is the additive update/correction term. For solving equation (2.14), the MC approximation to equation (2.12), a τ-discrete numerical scheme is required. Such a scheme would typically involve the following two steps.

Prediction

The predicted particle set at τ is generated using an Euler–Maruyama (EM) discretized map where , , etc.

Additive update

The update to the predicted particles is through an EM approximation to the integral in equation (2.14). Higher order integration schemes could also be considered, especially for evaluating the correction integral. The discrete update equation is presently given by where Δ=Δr, and the predicted solution, , is used to compute . Moreover, . Equation (2.16) may be recast as Recall from equation (2.4) that Δ:=Δτ so that equation (2.17) may be rearranged as where . When the particles are somewhat away from a local extremum (e.g. during the initial stages of evolution), the (norm of the) correction term is large. Hence the particles traverse more in the search space. In such cases, the innovation process would not behave as a zero-mean martingale as it would have a significant drift component. Since evolutions in this regime may have sharper τ-gradients, it is appropriate to modify the coefficient matrix in equation (2.18) so as to incorporate information on these gradients through previous estimates. Thus, and are replaced, respectively, by the following approximations: and

Note that we have used Ito's formula while approximating the r.h.s. in equation (2.19). Using equations (2.19) and (2.20), equation (2.18) may be modified as

Once the converged estimate is obtained, the innovation noise covariance should satisfy the identity:

Prior to convergence to an extremum, the (norm of the) r.h.s. of the equation above would typically be relatively large. Thus, one could impart higher diffusion to the search in the initial stages by replacing in equation (2.21) by Here α∈(0,1]. It is typically taken as 0.8 in the numerical illustrations based on pseudo-codes 1, 2 and 3. Equation (2.21) thus takes the final form:

A more concise form of the update equation is where the gain-like update coefficient matrix is given by

In the update strategy of equation (2.25), one may still improve on the search space exploration by multiplying the gain-weighted innovation term by a scalar factor β<1. This is equivalent to increasing the noise covariance and hence allowing the particles to be more diffusive or ‘explorative’. Thus, the update equation becomes

Coalescence, scrambling and relaxation: schemes for global search

The possible trapping of particles in local extrema is a major challenge to any global optimization scheme. For β≥1, the scheme described in equation (2.26) yields a greedy search that may often end up in a local extremum. Although it seems possible to avoid the local traps through the innovation noise, whose intensity could be tuned by the scalar factor β, such an approach may not be quite effective. As a local extremum is approached, the ‘strength’ (or norm) of the update term would be small and consequently its sensitivity to variations in β would also be poor. This makes the choice of β difficult (e.g. necessitating β to be too small for the search scheme to be efficient) and less effective for the global search. Moreover, smaller β also implies larger diffusion and hence poorer convergence. A more effective way out is a random perturbation applied to the particles so as to force out the ones trapped in the local wells. A general perturbation scheme, combining three basic approaches referred to as ‘coalescence’, ‘scrambling’ and ‘relaxation’, is now described. While the coalescence component is implemented through yet another martingale problem, the scrambling part requires a perturbation kernel. The relaxation component, on the other hand, requires accepting improvements with some probability. The importance of the random perturbation steps, in the present context, may also be gauged from the following fact. For a class of problems considered later, the faster convergence to a local extremum, engendered by the greedy scheme, equation (2.26), may force wrong convergence despite the application of random perturbations in conjunction with the greedy local search. Hence it would be worth exploring alternative search schemes, which completely eliminate scheme (2.26) and rely on the random perturbation steps alone.

Since the aim is to pose ‘coalescence’ as a martingale problem, the associated update should share the same generic structure as equation (2.26), which is recast below for the jth candidate where . The basic idea here is to provide to the evolving solution layers of random perturbation, whose inverse intensity may be formally indexed by a positive integer l such that the perturbation vanishes as . Within the τ-discrete setting, we start with the prediction and denote by the randomly perturbed, l-indexed increment so that as the random perturbations vanish asymptotically. During the ith iteration, the perturbed increment is arrived at using two transitional increments, u and v determined by two perturbation operators, say T1 and T2, respectively. While the operator T1 corresponds to the ‘local search’ and/or ‘coalescence’ operations, T2 encapsulates the ‘scrambling’ and ‘relaxation’ operations. The transitions may be indicated as , where T3 is a selection operator, commonly used with most evolutionary optimization schemes in some form or the other. Here is the final increment at τ, i.e. , which is input to the next iteration. denotes the predicted increment. Ideally, one may start the iterations with small l (i.e. high perturbation intensity), which is then gradually increased with progressing iterations. However, in the current implementations of COMBEO, the perturbation intensity is kept uniform all through the iterations. This is possible as, upon convergence to the global optimum, the applied perturbations merely yield zero-mean random fluctuations about the optimum, which are averaged out when the sample expectation operation is performed in the MC simulation. In the view of this and for notational ease, the left superscript l is removed from the variables in the discussion to follow and the perturbed nature of the variables should be clear from the context. The operators are now defined below.

Operator T1: local search and coalescence

The operation for local search has been described in §2 and through equation (3.1). The operation of ‘coalescence’ is now outlined. This perturbation is motivated by the observation that the probability density function (PDF, if it exists) associated with the converged measure π(.) should be unimodal, with its only peak located at the global extremum. Thus, when convergence to the global extremum occurs, all the particles should coalesce at the optimum point, except for a zero-mean noisy scatter around the latter. Ideally, for the sake of optimization accuracy, the noisy scatter should also have a low intensity so as to keep sample fluctuations under control. Once the global optimization scheme converges, the noisy scatter should then behave as a zero-mean martingale as a function of τ and with a unimodal transitional PDF. A zero-mean Brownian motion, which has a Gaussian PDF, is one such martingale. Clearly, such a property does not hold away from the global optimum, where the PDF would be multi-modal with a peak located at every local extremum detected by the algorithm.

Now, the aim is to obtain a scheme (or rather a family of schemes) to force the coalescence of particles. Consider the update of the jth particle such that coalescence of particles can be enforced. The random scatter around could be quantified by , where 1(j) denotes a random permutation on the indexing set {1,N}∖{j} based on a uniform measure. The goal is then to drive to a zero-mean vector Brownian increment with intensity matrix (typically assumed to be diagonal with the entries chosen uniformly for all j). In this case, ‘coalescence’ then implies that, in the limit of the noise intensity in δ(j) approaching zero, all the particles would tend to coalesce into a single location at the global extremum. Thus, similar to the innovation −f(x) on the left-hand side of equation (2.3), one treats as yet another innovation process so that, upon convergence, the identity holds almost surely. Here, is responsible for an additional layer of randomness that gives every particle the structure of a stochastic process. Accordingly, the extremal cost filtration must be suitably expanded/modified to include the sub-filtration generated by for s≤τ. If the coalescence innovation, in the form as indicated above, is included within our search process, equation (3.1) must be modified as where Recall that over-tildes indicate either the predicted particles or functions evaluated using the predicted particles, as appropriate. With a convenient abuse, the same notation for the gain-like update coefficient matrix is retained in equation (3.3) (used earlier in the local update equation (2.25)). may now be computed as and Here, the jth column of is given as . Let a be a positive real number and define the nearest integer smaller than a by ⌊a⌋. Then the integer valued perturbation parameter for the local extremization cum coalescence step may be identified as .

A real strength of the coalescence step is in the non-unique choice of the innovation vector—a feature that enables one to design powerful global search schemes. Indeed, one may also borrow some of the basic concepts from different existing global optimization schemes and adapt them within the present coalescence framework. For example, the personal and global information used in the PSO may be incorporated here by constructing the innovation as given below: Here is the personal best location corresponding to the jth particle along its evolution history until τ. is the available best location among all the particles in the population (ensemble) until τ, i.e. during the same evolution history. In contrast to the PSO, the gain-type coefficient matrix here is founded on a sound probabilistic basis and hence its matrix structure enables iteration-dependent differential weighting of the scalar components of the innovation vector, thereby yielding faster convergence without losing the exploratory efficacy of the global search. Indeed, it is well recognized that the use of three parameters (e.g. the cognitive and social learning factors and the inertia weight) in the original PSO may lead to solutions that sensitively depend upon the choice of these parameters. Even though it is possible to incorporate within the current set-up, the basic ideas behind some of the augmented PSO schemes, which attempt at removing some of the shortfalls of the original PSO [29], such a detailed exercise is kept outside the ambit of the current work. While we do not provide a proof here, the convergence and uniqueness of the iterative increment through this step, for a non-decreasing sequence of l converging to a limit point l* which may be large, could be shown based on the seminal work of Stroock & Varadhan [20] on martingale problems.

Operator T2: scrambling and relaxation

Similar to T1, the operator T2 also corresponds to random perturbations of the particles in the population. Identifying the gain-like coefficient matrix in the update equation (3.1) as a derivative-free stochastic counterpart to the Frechet derivative, the term may be considered the equivalent of a directional derivative term responsible for updating the jth particle. Consequently, around any extremum, the L2(P) norm is likely to be small. This has the effect of rendering further updates of the jth particle small, leading to a possible stalling of the optimization scheme. In order to move out of these local traps, a possible way is to swap the gain-weighted directional information for the jth particle with that of another randomly selected one. This form of random perturbation at τ may be implemented by replacing the update equation, e.g. equation (3.2), by any one of the following two perturbed equations: or where is update vector originally computed for the jth particle via equation (3.2) and 2 a random permutation on the integer set {1,N}. Out of the two perturbed equations above to implement scrambling, equation (3.8) is presently adopted. Formally, this perturbation may be described by a probability kernel p on [1,N]×[1,N] such that: . Clearly, as (or as , i.e. as convergence to the global optimum takes place), the matrix p(i,j) should ideally approach the identity matrix, i.e. p(i,j)→δ, where δ is the Kronecker delta.

By borrowing a basic idea from the DE and probably at the cost of a somewhat slower convergence for a class of problems, a more effective modification of the above scrambling strategy may be contemplated with a view to ensuring that the particles explore even more in the search space. Noting that every particle is an n-dimensional vector, one may execute the scrambling operation separately (n times) for individual scalar components of the particles instead of swapping the update terms for the particles as a whole. This may further be followed up by a so-called relaxation strategy, wherein the resulting modification is accepted with some probability γ<1. This improvement, enabling element crossovers across different particles, allows the particles to assume larger variations and prompt them to explore more. The resulting update equation is then given by where, for instance, denotes the updated mth scalar component of the jth particle at τ. Since, upon convergence (and possibly owing to the coalescence step, if applied), all particles crowd around the peak of a unimodal PDF with progressing iterations, directional scrambling with relaxation across such converged particles should not, in any way, affect the numerical accuracy of the estimated global extremum. Hence, in practical implementations of our scheme, p(i,j) need not strictly approach the identity matrix for large l.

Operator T3: selection

Use of diffusion-based random perturbations during exploration might sometimes result in ‘bad’ candidates. This may necessitate a selection step wherein candidates for the next iteration (say, the ith iteration) are chosen based on some selection criteria quantified by a selection or fitness function g(|x), ∈Ω. A general construction of the function, which is a Markov transition kernel on Ω and is conditioned on the ensemble of particles at the last iteration, should satisfy the following properties:

(i)  if j≠k; ∀ j,k∈[1,N],

(ii)  where ς∈(0,1], and

(iii) .

The updated jth particle appearing in the above clauses is computed using equation (3.7) (or equation (3.8)), which may incorporate different combinations of the three operations, e.g. local extremization, coalescence and scrambling-cum-relaxation. For the operator T3, the integer-valued inverse-perturbation parameter l may be identified with l=⌊1/(1−ς)⌋. In the current numerical implementations of COMBEO as described in §4, ς=1 is consistently adopted. This corresponds to l being infinity across all iterations and implies that the selection procedure is deterministic.

This work has primarily been aimed at the proposal for a class of new evolutionary optimization schemes and a verification of its performance purely through the numerical route. Accordingly, a detailed convergence analysis of the asymptotic dynamics, based on a combination of the martingale theory of Stroock–Varadhan and the random perturbation theory of Freidlin–Wentzell, will be considered in a separate study.

Algorithm development

The local search and random perturbation approaches, as described in §§2 and 3, merely provide a set of general tools whose different combinations could lead to different evolutionary schemes for optimization. For instance, incorporation of the innovation term in equation (2.3) yields a greedy algorithm that, despite possibly faster convergence, is likely to be a poor performer in the global search. On the other hand, to facilitate a more exhaustive global search at the cost of a substantively slower convergence, only a PSO-type innovation as in equation (3.6) could be adopted. For better clarity and a more objective assessment of the strengths and demerits of different search tools, three pseudo-codes are presented in this section. A more ambitious alternative would have been to combine all the presented ideas in a single pseudo-code that could automatically offer a right mix of convergence speed with exploratory efficiency depending on the nature of the problem at hand. This exercise however needs a non-trivial extension of the current work and hence constitutes an interesting future problem at this stage. pseudo-code 1 uses the innovation term of equation (2.3) for a greedy local search while employing scrambling of particles as a whole (not element wise) for the global search. Importantly, it may be noted that while the prediction step of equation (2.1) appears to be helpful in exploration, this step is not practically useful with our schemes, especially as it does not exploit any directional information and as more efficient tools for explorations have already been laid out. Hence, given a random initial scatter to the particles provided at the beginning of the iterations, the prediction step has been completely eliminated from all subsequent iterations. This would typically mean that the number of evaluations of the objective functional would be reduced by half. Note that, as there is no prediction step, the over-tildes in the notations of the variables are henceforth removed.

COMBEO pseudo-code 1

1. Discretize the τ-axis, say , using a partition such that τ0<⋯<τ and τ−τ=Δτ (=Δτ:=1/M if a uniform step-size is chosen for i=0,…,M−1). Assign τ0=1 and adopt Δτ small (approx. 10−7). Choose an ensemble size N.

2. From the domain of definition randomly generate, following a uniform distribution, the ensemble of initial particles (the initial population) for the solution vector. For each discrete τ, i=1,…,M−1, execute the following steps.

3. Additive update

Choose α∈(0,1); a typical prescription could be α≈0.8, even though the method generally performs well for other values in the interval indicated. Update each particle as where is the ith correction vector. 2(j) is the jth random permutation based on a uniformly distributed probability measure over the integer set{1,…,N}. For convenience, the expression for is reproduced below: ( is a scalar constant and ζ a uniform random number between 0 and 1). 1(j) is defined as the jth candidate from another (independent) random permutation of the integer set {1,…,N}. Finally,

4. If then retain as the updated particle; else set .

5. If i

In the next pseudo-code, the DE-like scalar element-wise scrambling scheme for particles is implemented and the resulting modification accepted with some relaxation criteria. The greedy local search based on the innovation term is not included in this pseudo-code as it might negatively bear on the global exploration. It may be possible to modify the pseudo-code to augment the innovation vector with additional innovation terms like after properly weighing the different terms. However, such a variation is presently not attempted.

COMBEO pseudo-code 2

1. Replicate steps 1 and 2 from pseudo-code 1.

2. Additive update

Choose α∈(0,1) as in pseudo-code 1 and update each particle as described below. Draw a random number r uniformly from the set {1,…,n}. As in pseudo-code 1, define σ2(j) as the jth candidate from a uniformly random permutation over the integer set {1,…,N}.

Initialize: , j=1,…,N.

If ζ1]), then with as in pseudo-code 1 and , where σ1(j) is a uniformly random permutation on the index set {1,N}∖{j};

else Set r→r+1. Run the loop on r till it reaches n.

3. Replicate steps 4 and 5 from pseudo-code 1.

By using the available personal and global best information as in the PSO, one may propose yet another global optimization pseudo-code, as suggested earlier while discussing the coalescence step. Thus, the search here is based purely on a martingale problem without requiring a random perturbation, especially the scrambling step. This is an interesting proposal in the sense that the search for the global optimum is attempted by just constructing an innovation vector and driving the same to a zero-mean martingale. In the previous two algorithms (pseudo-codes 1 and 2), the update terms were not constructed using past information beyond one step. In pseudo-code 3 presented below, the update term is however based on the evolution history until the last iteration in order to explore the search space more exhaustively. Without a loss of generality, the specific algorithm given below is for a global maximization problem.

COMBEO pseudo-code 3

1. Same as in pseudo-code 1.

2. Generate the initial population .

Initialize: and .

Initialize the local and global best particles and . Initialize the correction term For each discrete τ, i=1,…,M−1, execute the following steps.

3. Additive update

Choose α∈(0,1).

Compute .

Update each particle as: is the jth update vector. The expression for is also reproduced below: ( is a scalar constant and ζ a uniform random number in [0,1]).

4. Update the personal and global information as:

if then retain as the updated particle;

if then ;

else x=x.

Before concluding this section, a word of caution regarding the claims made on a relatively superior performance of the present approach vis-a-vis a few others, e.g. the DE and the PSO, should be in order. Optimization schemes like the DE or the PSO have been extensively improved over the last 15 years or so and it would be wrong to qualify the current method (as captured through the pseudo-codes above) as being the most competent. Nevertheless, the inherent flexibility of our proposal should permit the incorporation of the basic concepts behind some of the recent improvements in the DE or the PSO [30-33], thereby rendering the latter even more effective. Such details however need to be worked out and do not form part of this study.

Numerical illustrations

Benchmark problems

While the efficiency of a global optimization strategy is largely dependent on its ability to search all the promising search regions, the complexity level of a given optimization problem could exponentially increase with increasing number of design variables. As an example, even though Rosenbrock's function is unimodal in two dimensions (n=2), it is highly multimodal in still higher dimensional search spaces. Hence, a global optimization scheme, which works successfully for lower dimensional problems, may very well fail as the dimension of the problem increases. Moreover, the level of difficulty in finding the global optimum may also be dependent on the specific characteristics of the objective function, which goes to explain the varied degree of difficulty in solving different optimization problems of the same dimension. If the problem is separable or partially separable, i.e. if the objective function can be additively split in terms of component functions, each of which is expressible in terms of just one element or a small subset of elements of the design variable vector x, the original problem may actually be decoupled into a set of sub-problems. Each sub-problem, involving one scalar variable or a small set of scalar variables, may be solved separately, e.g. by a simple line search for one scalar variable. Depending on the degree of complexity, optimization problems could be categorized as separable, m-non-separable and non-separable [34]. In between the two extreme cases, i.e. separable and non-separable, are the m-non-separable functionals (m being the maximum number of scalar design variables appearing in the descriptions of the component functions) that correspond to partially separable problems. The nomenclature separability bears a similar meaning as epistasis in biology. In this work, in testing the effectiveness of the proposed schemes to some extent, the following problems have been considered.

I. Separable functions:

(a) F1: shifted elliptic function

(b) F2: shifted Rastrigin's function

(c) F3: shifted Ackley's function

II. Single-group m-non-separable functions:

(d) F4: single-group shifted and m-rotated elliptic function

(e) F5: single-group shifted and m-rotated Rastrigin's function

(f) F6: single-group shifted m-dimensional Schwefel's Problem 1.2

(g) F7: single-group shifted m-dimensional Rosenbrock's function

III. n/2m-group m-non-separable functions:

(h) F8: n/2m-group shifted and m-rotated Rastrigin's function

(i) F9: n/2m-group shifted and m-rotated Schwefel's Problem 1.2

IV. n/m-group m-non-separable functions:

(j) F10: n/m-group shifted and m-rotated Schwefel's Problem 1.2

V. Non-separable functions:

(k) F11: shifted Schwefel's Problem 1.2

The above functions are constructed using the following basic functions.

B1 The sphere function:

B2 The elliptic function:

B3 The rotated elliptic function: Frot_elliptic(x)=Felliptic(z), z=x*M (M being an orthogonal matrix)

B4 Schwefel's Problem 1.2:

B5 Rosenbrock's function:

B6 Rastrigin's function:

B7 Rotated Rastrigin's function: Frot_Rastrigin(x)=FRastrigin(z), z=x*M (M being a orthogonal matrix)

B8 Ackley's function:

B9 Rotated Ackley's function: Frot_Ackley(x)=FAckley(z), z=x*M

Explicit expressions for the functions F1–F11 are given below. Let z=x−o. denote a shifted global optimum and σ the random permutation of the vector {1,…,n}.

(a) F1(x)=Felliptic(z)

(b) F2(x)=FRastrigin(z)

(c) F3(x)=FAckley(z)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k) F11(x)=FSchwefel(z).

In solving optimization problems involving functions F1–F11, the performance of pseudo-code 2, as compared with the parent DE, is given in table 1. The population set in both the cases consists of 2000 particles. The scalar parameter c, which is similar to the acceptance criteria given by the crossover constant (CR) in DE [3] is taken as 0.1. In implementing DE, CR is also taken as 0.1. For both the schemes, the programme is terminated if the number of iterations exceeds the maximum allowed threshold, which is presently given by it_max=1×105. In all the tables, ε denotes an error norm ≤10−5. By the way of investigating the consistent reproducibility of the results reported herein, the proposed optimization scheme (pseudo-code 2) and the DE are run five times for F3. It is observed that in all the cases pseudo-code 2 provides the optimum in roughly the same number of iterations (same range) mentioned in table 2, whereas the DE fails to converge in all the runs. In the literature, the above benchmark problems are typically considered to be ‘difficult’ to solve by most of the existing evolutionary optimization schemes [34].

Table 1.

Comparative performance of pseudo-code 2 and DE against objective functions F1–F11 (n=40). (OF, objective functional; NI, number of iterations; MEN, mean error norm.)

	pseudo-code 2		DE
OF	NI	MEN	NI	MEN
F1	1.69×103	ε	2.42×104	ε
F2	4.68×103	ε	2.16×104	ε
F3	3.73×103	ε	it_max	20
F4	1.34×104	ε	it_max	1.2×104
F5	3×103	ε	2.37×104	ε
F6	2.51×103	ε	2.12×104	ε
F7	1.11×104	ε	2.81×104	ε
F8	4.15×104	ε	5.84×104	ε
F9	2.29×103	ε	3.62×104	ε
F10	8.84×103	ε	it_max	575
F11	it_max	7.7×10−4	it_max	5.78×103

Table 2.

Comparative performance of pseudo-code 3, pseudo-code 1 and PSO against objective functions B1–B9 (n=40). (OF, objective functional; NI, number of iterations; MEN, mean error norm.)

	pseudo-code 3 (pseudo-code 1)		PSO
OF	NI	MEN	NI	MEN
B1	146 (37)	ε	330	ε
B2	172 (55)	ε	353	ε
B3	212 (64)	ε	630	ε
B4	190 (49)	ε	456	ε
B5	154 (35)	ε	280	ε
B6	155 (42)	ε	353	ε
B7	172 (46)	ε	330	ε
B8	195 (61)	ε	485	ε
B9	201 (59)	ε	440	ε

Indeed, for these problems, the PSO, pseudo-codes 1 and 3 could not converge to the global optimum within the presently allowed number of iterations. However, in order to demonstrate that pseudo-codes 1 and 3 could indeed be very accurate with faster convergence for many global optimization problems, these codes are tested with benchmark global optimization problems B1–B9. While these problems are also challenging in their own right and often used in the literature to assess the performance of optimization algorithms, they do pose a comparatively lesser degree of complexity in reaching the global optimum vis-a-vis F1–F11. The results, reported in table 2, are also compared against the PSO under a similar environment. In all the cases, the population set is taken to consist of 50 particles. As anticipated, pseudo-code 3 and especially pseudo-code 1 show conspicuously faster convergence than the PSO while maintaining the same or higher level of accuracy. In any case, there is clearly scope for improving these two pseudo-codes (1 and 3) so that they can be rendered as competitive as pseudo-code 2. In exploring modification over the present algorithms so as to render them more robust and applicable to higher dimensional optimization problems, one possible way, for instance, would be intelligent augmentation of multiple adaptive strategies [35]. Moreover, since there are no specific formulae for selecting the parameter values in the present algorithms, parameter control mechanisms within the martingale problem set-up may also have to be suitably devised [35].

In order to gain further understanding on the performance of COMBEO, pseudo-code 2 is tested against yet another set of benchmark problems with the cost functionals, herein denoted as F(I)1–24, as considered in (http://coco.gforge.inria.fr/doku.php). In evolutionary algorithms, tuning of some of the parameters could be a tricky issue [36,37]. Since the performance of a scheme may crucially depend upon the choice of such parameters [36,37], which are often left to be selected at the user-end, a robust scheme should have the number of such tunable parameters kept minimum. Note that only one tuning parameter (CR) is introduced in the proposed scheme (pseudo-code 2). A detailed comparison of the performance of CMA-ES vis-a-vis COMBEO with varying CR values (e.g. 0.1, 0.9 and a number randomly drawn from a uniform distribution over 0 and 1) in the context of the new set of benchmark problems as considered in (http://coco.gforge.inria.fr/doku.php) has been undertaken and the results are reported in tables 3 and 4, only for the case of n=20. In order to limit the length of the article, a more detailed comparison for different values of n is provided in the electronic supplementary material. It may be noted that, in order to maintain a fair comparison between COMBEO and CMA-ES, the realized state with the minimum cost functional value in the available ensemble is considered as the solution. In these tables, EFE, error 1, error 2, error 3 and error 4, respectively, denote the average number of function evaluations among independent restarts, mean error, best error, worst error and median error. Tables 3 and 4 and electronic supplementary material, tables SI–SX, also include simulation results via CMA-ES [6]. The source code for CMA-ES is taken from https://www.lri.fr/~hansen/cmaes_inmatlab.html.

Table 3.

Comparative performance of pseudo-code 2 (with varying CR) and CMAES against objective functions F(I)1–F(I)14 (n=20).

F(I) index	method	EF(I)E	error 1	error 2	error 3	error 4
F(I)1	pseudo-code 2(CR=0.1)	32 909	8.70×10−9	8.70×10−9	8.70×10−9	8.70×10−9
F(I)1	pseudo-code 2(CR=rand)	50 909	7.98×10−9	7.98×10−9	7.98×10−9	7.98×10−9
F(I)1	pseudo-code 2(CR=0.9)	1 000 000	4.95×10−5	3.24×10−5	5.98×10−5	4.44×10−5
F(I)1	CMA-ES	4537	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)2	pseudo-code 2(CR=0.1)	46 818	8.31×10−9	8.31×10−9	8.31×10−9	8.31×10−9
F(I)2	pseudo-code 2(CR=rand)	70 455	6.67×10−9	6.67×10−9	6.67×10−9	6.67×10−9
F(I)2	pseudo-code 2(CR=0.9)	1 000 000	4.15×10−1	3.30×10−1	6.39×10−1	4.09×10−1
F(I)2	CMA-ES	21 147	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)3	pseudo-code 2(CR=0.1)	72 909	9.50×10−9	9.50×10−9	9.50×10−9	9.50×10−9
F(I)3	pseudo-code 2(CR=rand)	1 000 000	5.68×100	5.32×100	9.32×100	5.32×100
F(I)3	pseudo-code 2(CR=0.9)	1 000 000	9.49×101	8.92×101	1.07×102	9.71×101
F(I)3	CMA-ES	7822	2.33×101	1.99×101	2.59×101	2.39×101
F(I)4	pseudo-code 2(CR=0.1)	79 727	6.87×10−9	6.87×10−9	6.87×10−9	6.87×10−9
F(I)4	pseudo-code 2(CR=rand)	1 000 000	9.32×100	8.89×100	1.20×101	8.89×100
F(I)4	pseudo-code 2(CR=0.9)	1 000 000	1.04×102	1.01×102	1.06×102	1.06×102
F(I)4	CMA-ES	7421	3.62×101	2.69×101	5.27×101	2.89×101
F(I)5	pseudo-code 2(CR=0.1)	4727	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)5	pseudo-code 2(CR=rand)	3455	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)5	pseudo-code 2(CR=0.9)	2818	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)5	CMA-ES	340	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)6	pseudo-code 2(CR=0.1)	1 000 000	1.49×10−2	1.42×10−2	2.08×10−2	1.42×10−2
F(I)6	pseudo-code 2(CR=rand)	1 000 000	2.06×10−1	1.58×10−1	2.65×10−1	1.58×10−1
F(I)6	pseudo-code 2(CR=0.9)	1 000 000	5.54×102	4.62×102	6.44×102	5.65×102
F(I)6	CMA-ES	16 332	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)7	pseudo-code 2(CR=0.1)	1 000 000	1.15×100	1.05×100	1.50×100	1.11×100
F(I)7	pseudo-code 2(CR=rand)	1 000 000	1.46×100	1.05×100	1.82×100	1.74×100
F(I)7	pseudo-code 2(CR=0.9)	1 000 000	3.97×101	3.81×101	4.02×101	4.02×101
F(I)7	CMA-ES	3126	4.99×100	4.07×100	6.84×100	4.68×100
F(I)8	pseudo-code 2(CR=0.1)	1 000 000	1.20×101	1.17×101	1.26×101	1.17×101
F(I)8	pseudo-code 2(CR=rand)	1 000 000	1.52×101	1.51×101	1.56×101	1.51×101
F(I)8	pseudo-code 2(CR=0.9)	1 000 000	1.92×101	1.88×101	2.00×101	1.92×101
F(I)8	CMA-ES	21 876	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)9	pseudo-code 2(CR=0.1)	1 000 000	1.73×101	1.72×101	1.73×101	1.73×101
F(I)9	pseudo-code 2(CR=rand)	1 000 000	1.66×101	1.65×101	1.67×101	1.67×101
F(I)9	pseudo-code 2(CR=0.9)	1 000 000	1.83×101	1.81×101	1.85×101	1.84×101
F(I)9	CMA-ES	22 236	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)10	pseudo-code 2(CR=0.1)	1 000 000	2.83×104	2.65×104	3.14×104	2.65×104
F(I)10	pseudo-code 2(CR=rand)	1 000 000	1.01×105	8.36×104	1.55×105	8.36×104
F(I)10	pseudo-code 2(CR=0.9)	1 000 000	3.79×105	3.53×105	6.37×105	3.53×105
F(I)10	CMA-ES	21 391	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)11	pseudo-code 2(CR=0.1)	1 000 000	6.29×101	6.29×101	6.29×101	6.29×101
F(I)11	pseudo-code 2(CR=rand)	1 000 000	1.69×102	1.44×102	2.19×102	1.60×102
F(I)11	pseudo-code 2(CR=0.9)	1 000 000	1.31×102	1.15×102	2.23×102	1.15×102
F(I)11	CMA-ES	17 781	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)12	pseudo-code 2(CR=0.1)	1 000 000	1.13×10−1	2.05×10−2	1.33×10−1	1.33×10−1
F(I)12	pseudo-code 2(CR=rand)	1 000 000	9.25×10−2	1.52×10−2	3.05×10−1	1.52×10−2
F(I)12	pseudo-code 2(CR=0.9)	1 000 000	1.72×103	1.42×103	3.27×103	1.42×103
F(I)12	CMA-ES	36 406	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8
F(I)13	pseudo-code 2(CR=0.1)	1 000 000	5.08×10−1	1.60×10−1	1.31×100	4.58×10−1
F(I)13	pseudo-code 2(CR=rand)	1 000 000	5.11×10−2	3.30×10−2	8.28×10−2	3.30×10−2
F(I)13	pseudo-code 2(CR=0.9)	1 000 000	4.38×101	4.31×101	4.38×101	4.38×101
F(I)13	CMA-ES	8332	2.96×10−7	2.96×10−7	2.96×10−7	2.96×10−7
F(I)14	pseudo-code 2(CR=0.1)	1 000 000	2.91×10−3	2.46×10−3	4.67×10−3	2.46×10−3
F(I)14	pseudo-code 2(CR=rand)	1 000 000	5.81×10−3	4.56×10−3	8.44×10−3	4.56×10−3
F(I)14	pseudo-code 2(CR=0.9)	1 000 000	8.22×10−2	8.22×10−2	8.22×10−2	8.22×10−2
F(I)14	CMA-ES	36 078	−1.00×10−8	−1.00×10−8	−1.00×10−8	−1.00×10−8

Table 4.

Comparative performance of pseudo-code 2 (with varying CR) and CMAES against objective functions F(I)1–F(I)14 (n=20).

F(I) index	method	EF(I)E	error 1	error 2	error 3	error 4
F(I)15	pseudo-code 2(CR=0.1)	1 000 000	8.66×101	8.11×101	9.36×101	8.54×101
F(I)15	pseudo-code 2(CR=rand)	1 000 000	8.94×101	8.48×101	1.05×102	8.48×101
F(I)15	pseudo-code 2(CR=0.9)	1 000 000	1.22×102	1.21×102	1.33×102	1.21×102
F(I)15	CMA-ES	7657	1.69×101	1.69×101	1.69×101	1.69×101
F(I)16	pseudo-code 2(CR=0.1)	1 000 000	1.40×101	1.09×101	1.64×101	1.41×101
F(I)16	pseudo-code 2(CR=rand)	1 000 000	1.55×101	1.52×101	1.71×101	1.52×101
F(I)16	pseudo-code 2(CR=0.9)	1 000 000	1.94×101	1.83×101	2.20×101	1.83×101
F(I)16	CMA-ES	11 041	7.93×10−1	6.31×10−1	2.41×100	6.31×10−1
F(I)17	pseudo-code 2(CR=0.1)	1 000 000	1.18×10−2	8.39×10−3	1.70×10−2	1.20×10−2
F(I)17	pseudo-code 2(CR=rand)	1 000 000	9.43×10−3	8.35×10−3	1.34×10−2	8.35×10−3
F(I)17	pseudo-code 2(CR=0.9)	1 000 000	5.46×10−1	5.03×10−1	6.17×10−1	5.44×10−1
F(I)17	CMA-ES	12 505	7.83×10−3	3.43×10−3	1.40×10−2	4.95×10−3
F(I)18	pseudo-code 2(CR=0.1)	1 000 000	3.45×10−1	2.93×10−1	4.10×10−1	3.39×10−1
F(I)18	pseudo-code 2(CR=rand)	1 000 000	1.65×10−1	1.50×10−1	1.90×10−1	1.50×10−1
F(I)18	pseudo-code 2(CR=0.9)	1 000 000	3.22×100	3.18×100	3.59×100	3.18×100
F(I)18	CMA-ES	12 964	4.95×10−2	4.12×10−2	7.00×10−2	4.12×10−2
F(I)19	pseudo-code 2(CR=0.1)	1 000 000	3.07×100	2.75×100	3.30×100	3.12×100
F(I)19	pseudo-code 2(CR=rand)	1 000 000	2.65×100	2.53×100	3.89×100	2.53×100
F(I)19	pseudo-code 2(CR=0.9)	1 000 000	4.29×100	4.02×100	4.95×100	4.02×100
F(I)19	CMA-ES	147 687	1.62×100	7.59×10−1	4.16×100	7.59×10−1
F(I)20	pseudo-code 2(CR=0.1)	1 000 000	5.07×10−2	4.96×10−2	6.16×10−2	4.96×10−2
F(I)20	pseudo-code 2(CR=rand)	1 000 000	6.86×10−1	6.43×10−1	6.98×10−1	6.98×10−1
F(I)20	pseudo-code 2(CR=0.9)	1 000 000	2.41×100	2.38×100	2.44×100	2.44×100
F(I)20	CMA-ES	5963	1.57×100	1.39×100	1.81×100	1.55×100
F(I)21	pseudo-code 2(CR=0.1)	1 000 000	1.31×10−2	9.34×10−4	1.44×10−2	1.44×10−2
F(I)21	pseudo-code 2(CR=rand)	87 364	9.12×10−9	9.12×10−9	9.12×10−9	9.12×10−9
F(I)21	pseudo-code 2(CR=0.9)	1 000 000	2.07×100	1.52×100	2.54×100	2.51×100
F(I)21	CMA-ES	5113	7.29×100	6.92×10−1	6.44×101	1.80×100
F(I)22	pseudo-code 2(CR=0.1)	1 000 000	1.95×100	1.95×100	1.95×100	1.95×100
F(I)22	pseudo-code 2(CR=rand)	1 000 000	1.95×100	1.95×100	1.95×100	1.95×100
F(I)22	pseudo-code 2(CR=0.9)	1 000 000	1.96×100	1.96×100	1.96×100	1.96×100
F(I)22	CMA-ES	4753	2.52×100	6.92×10−1	1.83×101	6.92×10−1
F(I)23	pseudo-code 2(CR=0.1)	1 000 000	1.41×100	1.12×100	1.58×100	1.58×100
F(I)23	pseudo-code 2(CR=rand)	1 000 000	1.28×100	1.16×100	1.83×100	1.16×100
F(I)23	pseudo-code 2(CR=0.9)	1 000 000	1.88×100	1.73×100	2.34×100	1.73×100
F(I)23	CMA-ES	3823	2.90×100	2.64×100	3.15×100	3.04×100
F(I)24	pseudo-code 2(CR=0.1)	1 000 000	1.05×102	1.03×102	1.12×102	1.04×102
F(I)24	pseudo-code 2(CR=rand)	1 000 000	9.77×101	9.27×101	1.12×102	9.27×101
F(I)24	pseudo-code 2(CR=0.9)	1 000 000	1.30×102	1.25×102	1.36×102	1.30×102
F(I)24	CMA-ES	6175	4.67×101	3.89×101	4.96×101	4.96×101

As we see from tables 3 and 4 and the ones included in the electronic supplementary material, while COMBEO generally scores in achieving higher accuracy while approaching the global optimum, CMA-ES does have a distinct computational advantage in that it can perform with substantively lower ensemble sizes, N. It would thus be worthwhile to incorporate some of the adaptive features of CMA-ES within COMBEO in a future work.

An application to quantitative estimation of density variation in high-speed flows through inversion of light travel-time data

Having assessed the proposed global optimization framework through a variety of benchmark problems, we now consider a considerably larger dimensional optimization problem that involves the quantitative estimation of density variation in a high-speed flow obstructed by a blunt-nose aerodynamic vehicle. This example is taken up to demonstrate that, while pseudo-code 2 often provides an effective search tool for the global extremum of complex objective functionals, a greedier search scheme (as in pseudo-code 1 or its possible variants) could nevertheless be more appropriate for many practically useful, yet large dimensional, problems where a scheme like pseudo-code 2 (or most available global schemes such as the DE or the PSO) might prove computationally prohibitive.

The experimental set-up is shown in figure 1, see [38] for more details. Here, the region of interest (ROI), the flow around the object (in this case, a blunt-nosed missile model), is illuminated by a plane wave. The distorted plane wave trans-illuminates a random dot pattern (RDP). The RDP in effect serves the purpose of the expensive lenslet split Shack–Hartmann sensor. The geometric shadow cast at an adjacent plane is imaged by a high-speed camera. The shadow, in comparison to the original RDP, is space-shifted. Local mean shifts =(d,d) and a related quantity, the slopes of the wavefront (d/d,d/d), are estimated by cross-correlating local sub-region of the distorted RDP with the corresponding original and finding the shift of the cross-correlation peak. The estimated local slopes are integrated to find the smooth wavefront Φ [39]. The recovered unwrapped Φ, which carries the distortion owing to its passage through the ROI, is shown in figure 2. The distorted wavefront Φ (or equivalently the phase delay) evaluated at the ith detector, is related to the refractive index distribution ϱ(r) in the ROI, r=(x,y,z) being the position vector, through the equation: Here L is the length of the ith ray and ds is an element from the Fermat's path of the light-ray through ϱ(r), which obeys the Eikonel equation:

Figure 1.

A schematic diagram of the experimental set-up.

Figure 2.

A three-dimensional view of the reconstruction via the proposed framework.

Therefore, the equation connecting Φ to ϱ(r) is nonlinear. The collection of Φ at all the detectors is denoted using the vector :={Φ}. The problem of recovery of ϱ(r) from the experimentally measured data, denoted by exp, obtained from Φ corresponding to the locations of detectors, is posed as an optimization problem which is solved within the proposed framework. The corresponding objective functional may be written as However, given the large system dimension (n=1012) and with the algorithm requiring repetitive numerical inversions for Φ, most available global optimization schemes (e.g. the DE or the PSO) are rendered impracticable using the commonly available computing systems. An acceleration of the computing speed is however possible in the present framework by so modifying the innovation processes as to render the search greedier while retaining some features of the global search, especially the scrambling step. One way could be to start with pseudo-code 1. In order to make the algorithm faster, we might replace the innovation used in pseudo-code 1 by the difference of the objective functional fexp computed over two succeeding iterations. Unfortunately, this leads to instability in the algorithm owing to the fact that the innovation becomes one-dimensional in contrast to the number of design variables which is of order O(104). This may be remedied by splitting the original objective functional into 1012 individual ones (i.e. one corresponding to each design variable) and thus constructing 1012 separate scalar innovation processes for the new objective function vector. This enables solving the problem in only a few iterations (as few as just five iterations) with only 30 MC particles. Accordingly, we now redefine , where . The present modification requires the realized innovation I at the ith iteration in pseudo-code 1 (specifically in step 3) to be given as It is recalled that this modified algorithm uses only scrambling for effecting the global search. This is an interesting example of the practical advantages in being able to suitably design the innovation process within the present framework. It is worth noting that our efforts at rendering the DE and the PSO greedier by tuning their appropriate algorithmic parameters have failed to make them effective in solving the current problem.

The recovered ϱ(r) is converted into density distribution using the Gladstone–Dale equation. A cross section of the recovered density profile is shown in figures 2 and 3. The above reconstruction is verified (figure 3) by computing the flow density distribution close to the boundary layer by solving Navier–Stoke's equation using the commercially available CFD software Fluent. The comparison of the density variations near the boundary layer, shown in figure 3, is seen to be quite close.

Figure 3.

A projection of the reconstruction via the proposed framework.

Concluding remarks

Despite the widely reported success and appealing simplicity of the methodological aspects of many existing evolutionary schemes for global optimization, often inspired by biological or social observations, the specific forms of their update scheme in approaching the global extremum are generally not rigorously grounded. While a major aspect of this work is the attempt at finding an underlying logical thread that perhaps relates the basic ideas behind many such update strategies, the notion of a non-uniquely specifiable martingale approach in deriving the update term and its functional similarity with the directional derivative, used in most deterministic optimization schemes, has by far been the most fundamental aspect in this article. Such a directional search, free from the notion of classical derivatives, has indeed been responsible for superior convergence features of the COMBEO family of schemes in contrast to some of the well-known evolutionary methods, e.g. the PSO, the DE and the CMA-ES. An effective global exploration of the search space however requires several random perturbation strategies, of which the coalescence step is treated as part of the martingale problem. It is shown that a variant of the coalescence step may mimic a PSO-like search. A DE-like search, on the other hand, may be organized through a variant of scrambling, which is yet another random perturbation strategy involving directional information swaps. The framework of COMBEO, founded on the theory of stochastic processes, promises to yield highly competitive update schemes in order to meet the conflicting demands of faster convergence and efficacious exploration. As demonstrated via numerical work on a host of higher dimensional benchmark objective functions, the COMBEO schemes show consistently superior performance in comparison to several existing evolutionary methods, prominently the DE and the PSO. The flexibility of the proposed optimization paradigm is also demonstrated by evolving a greedier version of the search scheme and applying the same to the reconstruction of boundary layer density variation in a high-speed flow experiment.

The generality of the framework afforded by COMBEO should motivate further research in improving the specific numerical schemes considered in this work. For instance, the COMBEO algorithms could possibly be rendered computationally more competitive by incorporating the adaptive features of the CMA-ES or some of the many available modifications or augmentations of the DE.