2. Stochastic Global Optimization

2.1. Why stochastic global optimization

DSGE and VAR likelihoods are routinely multimodal, with shallow ridges, flat plateaus, and parameter-region dependence on regime probabilities. Local optimizers (fmincon, csminwel, newrat) converge to the basin they happen to start in. Stochastic global optimizers trade per- iteration cost for coverage: they sample the bounded box [lb, ub] with a population of candidate solutions, exchange information across candidates, and contract toward the best region only as the iteration budget is consumed. They are the right tool when:

  • the model has more than a handful of parameters;

  • the user does not have a strong prior on a starting region;

  • a posterior maximization is intended to seed an MCMC sampler that must explore beyond a single mode.

RISE ships a unified family of nine global optimizers (plus the legacy @bee and the deterministic mads baseline), all sharing the same calling convention and outer-loop machinery.

2.2. The unified shell

Every algorithm in the +globalopt package is a thin entry point that wires an algorithm-specific step function into globalopt.generic_loop. The shell handles:

  • initial-population construction;

  • iteration counting and budget checks (MaxIter, MaxFunEvals, MaxTime);

  • progress display (Display = 'iter' | 'on' | 'off');

  • manual stopping via a ManualStoppingFile.txt sentinel in the current directory (see utils.optim.manual_stopping);

  • optional finite-difference Hessian at the optimum (returned as the 4th output).

This keeps every algorithm’s implementation small: the step function implements only the algorithm-specific update; everything else is shared.

2.3. Common calling convention

All nine algorithms expose the same signature, both via the package and via the rise_<name> shortcut:

[xfinal, ffinal, exitflag]    = rise_lshade(fh, x0, lb, ub)
[xfinal, ffinal, exitflag]    = rise_lshade(fh, x0, lb, ub, options)
[xfinal, ffinal, exitflag]    = rise_lshade(fh, x0, lb, ub, options, hyper)
[xfinal, ffinal, exitflag, H] = rise_lshade(...)

Inputs:

  • fh – function handle, fval = fh(x).

  • x0 – initial guess column vector. Used by some algorithms to seed the initial population (e.g. CMA-ES’s initial mean); ignored by others.

  • lb, ub – column vectors with lower and upper bounds.

  • options (optional struct) – the MATLAB-style fields below.

  • hyper (optional struct) – algorithm-specific hyperparameters (see each algorithm’s entry).

Outputs:

  • xfinal – best parameter vector found.

  • ffinal – best objective value.

  • exitflag – 1 on success.

  • H (optional 4th output) – finite-difference Hessian at xfinal (computed lazily on request).

2.3.1. Common options (every algorithm)

field

meaning

default

MaxIter

outer-loop iteration cap

200

MaxFunEvals

total function-eval cap

inf

MaxTime

wall-clock cap (seconds)

inf

Display

'iter' / 'on' / 'off'

'off'

2.4. The algorithm catalogue

The nine algorithms group naturally into four families. Each entry gives the philosophy, the seminal reference, and the algorithm-specific hyper fields.

2.4.1. Differential Evolution family

These three algorithms share the DE skeleton – a population of candidates, mutation by weighted vector differences, binomial crossover, greedy selection – and add successive refinements.

DE (rise_de, +globalopt.de)

The original Storn-Price algorithm (1997). Strategy selectable via hyper.strategy: 'rand/1', 'best/1', 'current-to-best/1', 'rand/2', 'best/2'. Hyperparameters: NP (population, default 20 + 5*dim), F (mutation factor, 0.8), CR (crossover probability, 0.9), strategy ('rand/1').

L-SHADE (rise_lshade, +globalopt.lshade)

Tanabe & Fukunaga (2014), CEC 2014 winner. Three additions: JADE-style current-to-pbest/1/bin mutation with an external archive; per-individual F_i (Cauchy) and CR_i (Normal) sampled from a success-history memory (SHADE); linear population-size reduction from NP_init to NP_min as the budget is consumed. Hyperparameters: NP_init (min(18*dim, 50)), NP_min (4), H (history memory size, 6), p_best (greediness, 0.11), arc_rate (archive size, 2.6).

jSO (rise_jso, +globalopt.jso)

Brest, Maucec & Boskovic (2017), CEC 2017 winner. Refinement of L-SHADE: weighted mutation F_w * (pbest - x) + F * (r1 - r2) with F_w depending on the fraction of the budget consumed (suppresses overly large early steps); F capped at 0.7 before 60% of the budget; generation-dependent greediness; smaller archive (arc_rate = 1.0). Hyperparameters: NP_init (min(round(25*log(dim)*sqrt(dim)), 50)), NP_min (4), H (5), p_max (0.25).

2.4.2. Population / swarm metaheuristics

ABC (rise_abc, +globalopt.abc)

Artificial Bee Colony (Karaboga 2005, with the Akay-Karaboga 2012 per-dimension modification rate). Three roles per cycle: employed bees do one mutation per food source; onlookers run NP roulette-wheel mutations weighted by fitness; scouts uniformly replace food sources whose trial counter exceeds limit. Hyperparameters: NP (food sources, max(20, 5*dim)), limit (abandonment threshold, NP*dim), MR (per-dimension modification rate, 0.5), phi_low, phi_high (perturbation bounds, [-1, 1]).

ACO (rise_aco, +globalopt.aco)

Continuous Ant Colony Optimization, ACO_R (Socha & Dorigo 2008). Keeps a ranked archive of k solutions with Gaussian-kernel weights. Each iteration samples m new ants by picking a guide from the archive (roulette), sampling each coordinate from N(guide_d, sigma_d) with sigma_d set by the mean absolute distance in dimension d to the other archived solutions, then merging into the archive. Hyperparameters: k (archive, max(20, 5*dim)), m (ants per iteration, k), q (elitism, 0.1), xi (convergence, 0.85).

BBO (rise_bbo, +globalopt.bbo)

Biogeography-Based Optimization (Simon 2008). Habitats are sorted best-first; the top n_elite are preserved; each remaining habitat receives features (decision-variable values) by migration from randomly chosen “donor” habitats with rates set by I (immigration) and E (emigration). Features are then perturbed by Gaussian creep with std mut_sigma (as a fraction of ub - lb). Hyperparameters: NP (max(20, 5*dim)), I (1), E (1), mut_prob (0.05), mut_sigma (0.02), n_elite (2).

2.4.3. Knowledge-based

AGSK (rise_agsk, +globalopt.agsk)

Adaptive Gaining-Sharing Knowledge (Mohamed, Hadi & Mohamed 2020, CEC 2020 winner). At every generation the population is sorted by fitness; for each individual a per-dimension partition splits the n coordinates into a junior subset (local exploitation using rank-neighbours) and a senior subset (global exploration using top/middle/bottom slices). The junior fraction shrinks as the budget is consumed. Update factors K_F, K_R are sampled per individual from SHADE-style success-history memories. Hyperparameters are mostly self-tuning; the implementation here is the single-memory variant.

2.4.4. Evolution strategy

CMA-ES (rise_cma_es, +globalopt.cma_es)

Covariance Matrix Adaptation Evolution Strategy (Hansen 2016). The distinct one in the family: model-based rather than mutation-based. Maintains a multivariate Gaussian sampling distribution; updates its mean by weighted recombination of the best offspring, and its covariance via rank-1 and rank-mu updates that adapt to the local geometry. The default adaptation constants are derived from lambda and dim. Hyperparameters: sigma0 (initial step size, 0.3), lambda (population, 4 + floor(3*log(dim))), mu (parents, floor(lambda/2)).

2.4.5. Direct-search baseline

MADS (rise_mads, +globalopt.mads)

Mesh Adaptive Direct Search (Audet & Dennis 2006), in the coordinate-aligned variant closer to Generalized Pattern Search (GPS) than to OrthoMADS. Deterministic given x0 and the bounds: no internal randomness. Coordinate-aligned poll directions (+/- e_i), opportunistic polling (accept the first improving direction), per-dimension step sizes that contract on failure (contract = 0.5) and expand on success (expand = 2.0). Useful as a non-stochastic reference when comparing the population-based optimizers. Hyperparameters: init_step (0.1*(ub - lb)), tol_step (1e-8), expand (2.0), contract (0.5), opportunistic (true).

2.5. Calling from estimate

The optimizer named argument of estimate(...) accepts a function handle, a string, or a cell {name, options...}. All rise_<name> shortcuts are valid drop-ins:

m = estimate(m, data = mydata, optimizer = @rise_lshade);

m = estimate(m, data = mydata, ...
             optimizer = {'rise_lshade', 'MaxFunEvals', 5e4});

To override the algorithm hyperparameters at estimation time, write a small wrapper that closes over the desired hyper:

my_lshade = @(fh, x0, lb, ub, opts, varargin) ...
    rise_lshade(fh, x0, lb, ub, opts, struct('NP_init', 100));

m = estimate(m, optimizer = my_lshade);

2.6. Calling directly

Outside estimate, the algorithms are standalone optimizers and can solve any objective:

rosenbrock = @(x) sum(100*(x(2:end) - x(1:end-1).^2).^2 ...
                      + (1 - x(1:end-1)).^2);

n  = 10;
lb = -5*ones(n,1); ub = 5*ones(n,1);
x0 = lb + (ub - lb).*rand(n,1);

options = struct('MaxFunEvals', 1e4, 'Display', 'iter');

[xstar, fstar] = rise_lshade(rosenbrock, x0, lb, ub, options);

2.7. Manual stopping

Long runs can be stopped cleanly without losing the current best: create a file named ManualStoppingFile.txt (any content) in the working directory. utils.optim.manual_stopping is polled inside generic_loop; the next iteration boundary terminates the run and returns the current best.

2.8. The legacy @bee class

RISE also retains the older @bee class (m/optimizers/@bee/bee.m with m/optimizers/bee_gate.m as the function-style entry point). It exposes a different option set and remained the default for RISE’s MCMC posterior-mode search until the +globalopt family was introduced. New work should prefer one of the nine rise_<name> shortcuts above; @bee is kept for backward compatibility with existing scripts.

2.9. Algorithms not currently included

Two algorithms sometimes asked about – firefly (Yang 2008) and glowworm swarm optimization (Krishnanand & Ghose 2009) – are not in the +globalopt family. The family is opinionated toward algorithms that have been independently winners or runners-up of the CEC real-parameter optimization competition (DE, L-SHADE, jSO, AGSK, CMA-ES), plus the most-cited swarm metaheuristics (ABC, ACO, BBO), plus a deterministic direct-search baseline (MADS). Firefly and glowworm slot in cleanly if needed: drop a +globalopt/+firefly/ package with firefly.m (thin entry point wiring hyper.algorithm_step) and firefly_step.m (init / iterate modes), matching the existing 9 algorithms’ shape.