# Parallel And Hybrid Metaheuristics > When the user wants to parallelize a metaheuristic or combine several searches — island models with migration, master-slave fitness evaluation, parallel multistart, cooperative search, or algorithm portfolios — with the Python practicalities (multiprocessing, joblib, vectorization-first). Also use when the user mentions "parallel metaheuristic," "island model," "migration," "master-slave," "cooperative search," "portfolio," or when single runs are too slow and CPU cores sit idle. For removing Python loops inside one process, see numpy-vectorization-for-optimization; for choosing the base algorithm, see metaheuristic-design-principles. - Skill: `hajibabaie/parallel-and-hybrid-metaheuristics` (Agent Skill) - Install (CLI): `npx skillmds@latest add hajibabaie/parallel-and-hybrid-metaheuristics` - Raw SKILL.md: https://api.skillmd.com/api/skills/hajibabaie/parallel-and-hybrid-metaheuristics/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: AI & ML - Author: hajibabaie (https://skillmd.com/u/hajibabaie) - Updated: 2026-09-17 - Page: https://skillmd.com/skills/hajibabaie/parallel-and-hybrid-metaheuristics --- # Parallel and Hybrid Metaheuristics You are an expert in parallel and hybrid metaheuristics for combinatorial optimization. This skill covers the architectures that distribute or combine searches — island models with migration, master-slave fitness evaluation, independent multistart, cooperative search with shared elite memory, and algorithm portfolios — together with the Python practicalities (multiprocessing, joblib, picklable workers, seed management) that decide whether a parallel design actually pays off. Use the framework below to pick an architecture from the bottleneck profile, implement it with process pools and numpy, and report speedups honestly. ## Initial Assessment Establish the following before recommending an architecture or writing any code: - **Profile first.** What fraction of wall-clock time is fitness evaluation? What fraction is operators, copying, bookkeeping? Parallelizing a part that takes 10% of runtime caps the speedup at 1.11x regardless of core count. - **Cost of one evaluation.** Microseconds (pure arithmetic), milliseconds (delta-evaluated neighborhood scans), or seconds (simulation, embedded LP/ML model)? This single number selects the architecture: see the decision guidance below. - **Vectorization status.** Is the base algorithm already numpy-vectorized at the population level? One vectorized process routinely beats eight processes running slow Python loops. Climb the vectorization rung before the parallelism rung. - **Base method class.** Population-based (GA, DE, ACO) maps naturally to islands and master-slave; single-solution methods (ILS, SA, tabu) map to multistart, cooperative search, and portfolios. - **Hardware.** Physical cores (not hyperthreads), memory per worker, single machine or cluster. Most research workloads end at one machine with 8-64 cores; design for that first. - **Operating system and start method.** Windows and macOS spawn worker processes (fresh interpreter, everything pickled); Linux can fork. Spawn imposes module-level worker functions and `if __name__ == "__main__"` guards. - **Picklability.** Can the objective and instance data cross a process boundary? Open solver handles, database connections, GUI objects, and lambdas cannot. - **Goal type.** Shorter wall-clock to reach the same quality, or better quality within the same wall-clock? Island models and cooperation target the second; master-slave targets the first. - **Reproducibility requirement.** Bit-for-bit reproducible runs rule out asynchronous designs whose result depends on OS scheduling; synchronous patterns with `SeedSequence.spawn` keep determinism. - **Experiment-level parallelism.** A study running 10 instances x 10 seeds already has 100 independent jobs. Parallelizing across runs is trivial, perfectly efficient, and often all that is needed — check this before touching the algorithm's internals. - **Time budget per run and per study.** Communication-heavy designs amortize poorly over short runs; a 5-second run cannot absorb 0.5 s of pool startup. ## Parallel Architectures Five architectures cover practice. The classification follows Crainic & Toulouse (2010), "Parallel Metaheuristics" (Handbook of Metaheuristics), and Alba (2005), "Parallel Metaheuristics: A New Class of Algorithms". | Architecture | Grain | What runs in parallel | Communication | Changes search behavior? | |---|---|---|---|---| | Master-slave (global) | fine | fitness evaluations of one population | every generation, master <-> workers | no — same trajectory, faster | | Island model (coarse) | coarse | whole subpopulations | rare migration events | yes — restricted mixing | | Cellular (fine-grained) | very fine | per-individual updates on a grid | implicit, overlapping neighborhoods | yes — slow diffusion | | Independent multistart | run-level | full runs with different seeds | none | no — same per-run distribution | | Cooperative search | coarse | heterogeneous searches | shared elite pool / adaptive memory | yes — by design | Algorithm **portfolios** are multistart generalized across algorithms: run several different methods (or configurations) on the same instance and keep the best result (Huberman, Lukose & Hogg 1997, "An Economics Approach to Hard Computational Problems"; Gomes & Selman 2001, "Algorithm Portfolios"). ### Speedup arithmetic Speedup and efficiency with $p$ workers, where $T_p$ is wall-clock time: $$ S(p) = \frac{T_1}{T_p}, \qquad E(p) = \frac{S(p)}{p}. $$ Amdahl (1967): if only a fraction $\phi$ of the work parallelizes, $$ S(p) \le \frac{1}{(1-\phi) + \phi/p}. $$ With $\phi = 0.9$ (evaluation is 90% of runtime), 16 workers give at most 6.4x. This is why master-slave designs must move *all* heavy work, not just part of it, into workers. For independent multistart there is a second, quality-based argument. If one run reaches a target quality with probability $q$ within the budget, $p$ independent runs succeed with probability $1-(1-q)^p$. Empirically, time-to-target of well-designed GRASP/ILS-type methods is approximately exponentially distributed (Aiex, Resende & Ribeiro 2002, "Probability distribution of solution time in GRASP"), and the minimum of $p$ i.i.d. exponentials is exponential with rate multiplied by $p$ — independent multistart then achieves near-linear speedup in expected time-to-target with zero communication. ### Decision guidance - **Use master-slave when** one evaluation costs >= 1-10 ms and dominates runtime (simulation-based objectives, embedded solver calls). It preserves the algorithm exactly; only wall-clock changes, so no re-tuning is needed. - **Use an island model when** the method is population-based, evaluation is cheap, and the goal is better quality at equal time. Restricted migration preserves diversity that a single panmictic population loses (Whitley, Rana & Heckendorn 1999, "The island model genetic algorithm: on separability, population size and convergence"). Islands also give near-perfect parallel efficiency because communication is rare. - **Use independent multistart when** the base method is single-solution with high run-to-run variance (ILS, SA). It is the easiest correct parallelization and a mandatory baseline for any fancier design: a cooperative scheme that cannot beat independent multistart at equal core-count is not cooperating usefully. - **Use cooperative search when** different searches can exploit each other's elites — e.g., construction-heavy and improvement-heavy methods, or differently-tuned copies of one method. The shared elite pool is the "adaptive memory" of Rochat & Taillard (1995), "Probabilistic diversification and intensification in local search for vehicle routing". - **Use a portfolio when** no single algorithm dominates across instances or runtime is heavy-tailed; racing configurations hedges the choice. Portfolios compose with everything above (each entry may itself be parallel). - **Skip cellular models in Python** unless targeting GPUs; their very fine grain suits SIMD hardware, not process pools. ### The vectorization-first ladder Climb in order; each rung is cheaper to implement and to debug than the next. 1. **Vectorize within one process** — numpy population operations, batch fitness, delta evaluation. Often 10-100x; no new failure modes. 2. **Parallelize evaluation** — master-slave over a process pool; only when one call is expensive enough to amortize pickling (rule of thumb: task >= 1 ms, better >= 10 ms). 3. **Parallelize runs** — multistart, islands, portfolios; communication rare or absent. 4. **Distribute** — ray / mpi4py / dask across machines, only when one machine is provably insufficient. The Python global interpreter lock (GIL) serializes pure-Python bytecode across threads, so CPU-bound metaheuristics need *processes*, not threads. numpy releases the GIL inside large array operations, which is exactly why rung 1 works so well. ## Reference Implementation: Synchronous Island Model The island model is an algorithmic structure first and a parallelization second: the restricted-migration dynamics improve search even when islands execute sequentially in one process. The framework below is problem-independent — `init` and `step` are the only problem-specific parts (swap in any population method; see genetic-algorithms for a complete GA `step`). ```text ISLAND_MODEL(init, step, k, interval, m, topology, G): for i in 1..k: (P_i, F_i) <- init(rng_i) # rng_i spawned from one SeedSequence for gen in 1..G: for i in 1..k: # parallelizable: islands are independent (P_i, F_i) <- step(P_i, F_i, rng_i) if gen mod interval == 0: # rare, structured communication for i in 1..k: E_i <- copies of the m best of island i # select emigrants FIRST for i in 1..k: for each source j of island i under topology: replace the m worst of island i with E_j return best individual over all islands ``` ```python """Problem-independent synchronous island model over swappable (init, step) callables.""" from dataclasses import dataclass from typing import Callable import numpy as np InitFn = Callable[[np.random.Generator], tuple[np.ndarray, np.ndarray]] StepFn = Callable[[np.ndarray, np.ndarray, np.random.Generator], tuple[np.ndarray, np.ndarray]] @dataclass class IslandConfig: """Topology and migration policy of a synchronous island model.""" n_islands: int = 8 migration_interval: int = 25 # generations between migration events n_migrants: int = 2 # emigrants per island per event topology: str = "ring" # "ring" or "complete" def migrate(pops: list[np.ndarray], fits: list[np.ndarray], cfg: IslandConfig) -> None: """Copy each island's best into its neighbors, replacing the targets' worst. Emigrants are selected before any replacement, so a solution overwritten on its home island still travels intact. Minimization; modifies in place. """ k = len(pops) sel = [np.argpartition(f, cfg.n_migrants)[:cfg.n_migrants] for f in fits] movers = [(p[s].copy(), f[s].copy()) for p, f, s in zip(pops, fits, sel)] for i in range(k): sources = [(i - 1) % k] if cfg.topology == "ring" else \ [j for j in range(k) if j != i] inc = np.vstack([movers[j][0] for j in sources]) inc_fit = np.concatenate([movers[j][1] for j in sources]) worst = np.argpartition(fits[i], -inc.shape[0])[-inc.shape[0]:] pops[i][worst] = inc fits[i][worst] = inc_fit def run_island_model(init: InitFn, step: StepFn, cfg: IslandConfig, n_generations: int, seed: int) -> tuple[np.ndarray, float]: """Evolve cfg.n_islands subpopulations with periodic migration; return the best.""" rngs = [np.random.default_rng(s) for s in np.random.SeedSequence(seed).spawn(cfg.n_islands)] pops, fits = map(list, zip(*(init(r) for r in rngs))) for gen in range(1, n_generations + 1): for i in range(cfg.n_islands): pops[i], fits[i] = step(pops[i], fits[i], rngs[i]) if gen % cfg.migration_interval == 0: migrate(pops, fits, cfg) i = int(np.argmin([f.min() for f in fits])) j = int(np.argmin(fits[i])) return pops[i][j].copy(), float(fits[i][j]) if __name__ == "__main__": def init(rng: np.random.Generator) -> tuple[np.ndarray, np.ndarray]: """40 random binary strings; objective: minimize the number of zeros.""" pop = rng.integers(0, 2, size=(40, 64)) return pop, (64 - pop.sum(axis=1)).astype(float) def step(pop: np.ndarray, fit: np.ndarray, rng: np.random.Generator) -> tuple[np.ndarray, np.ndarray]: """Elitist bit-flip step applied to every population slot at once.""" child = np.where(rng.random(pop.shape) < 1 / 64, 1 - pop, pop) child_fit = (64 - child.sum(axis=1)).astype(float) keep = child_fit < fit return np.where(keep[:, None], child, pop), np.where(keep, child_fit, fit) sol, val = run_island_model(init, step, IslandConfig(n_islands=4), 600, seed=3) print(val) # Expected: 0.0 (the all-ones string is found) on nearly all seeds within # 600 generations; 4 islands of 40 with ring migration every 25 generations. ``` ### Parameter guidance | Parameter | Typical range | What it trades off | |---|---|---| | `n_islands` | 4-16 on one machine (one process per island at the parallel rung) | more islands: more diversity and parallelism, but each subpopulation gets smaller and weaker | | subpopulation size | 30-200 per island | within-island selection strength vs total evaluation budget per generation | | `migration_interval` | 10-50 generations | frequent: islands merge into one effective panmictic population; rare: islands converge separately and synergy arrives too late | | `n_migrants` | 1-5 (<= 5% of the subpopulation) | takeover speed vs diversity; migration is itself a selection operator | | `topology` | ring (default), complete, random k-regular | takeover time: ring slowest (most diversity), complete fastest (Cantu-Paz 2000, "Efficient and Accurate Parallel Genetic Algorithms") | | migrant policy | best-replace-worst (default); random-replace-random for gentler pressure | the selection pressure added by migration on top of within-island selection | | sync vs async | synchronous (default in Python) | reproducibility and simplicity vs no barrier wait at segment ends | ## Worked Example: Island GA with Migration for TSP A complete, parallel island GA for the symmetric TSP. The pattern is **segment-synchronous**: each process-pool task evolves one island for `migration_interval` generations, the master gathers the islands, performs ring migration, and dispatches the next segment. This keeps communication at one gather/scatter per segment, so parallel efficiency stays high. Operator details (order crossover, inversion mutation, tournament selection) are standard; see crossover-operators and mutation-and-perturbation-operators for the full catalogs and genetic-algorithms for panmictic GA design. ```python """Island GA for the symmetric TSP: ring migration, islands in worker processes.""" from concurrent.futures import ProcessPoolExecutor import numpy as np def tour_lengths(pop: np.ndarray, dist: np.ndarray) -> np.ndarray: """Vectorized lengths of a population of closed tours; pop has shape (P, n).""" return dist[pop, np.roll(pop, -1, axis=1)].sum(axis=1) def nearest_neighbor(dist: np.ndarray, start: int) -> np.ndarray: """Nearest-neighbor construction; standard warm start for TSP populations.""" n = dist.shape[0] unvisited = np.ones(n, dtype=bool) unvisited[start] = False tour = [start] for _ in range(n - 1): cand = np.flatnonzero(unvisited) nxt = int(cand[np.argmin(dist[tour[-1], cand])]) tour.append(nxt) unvisited[nxt] = False return np.array(tour) def order_crossover(p1: np.ndarray, p2: np.ndarray, rng: np.random.Generator) -> np.ndarray: """OX: keep a random slice of p1, fill the rest in p2's order.""" n = p1.shape[0] a, b = np.sort(rng.choice(n, size=2, replace=False)) child = np.full(n, -1, dtype=p1.dtype) child[a:b + 1] = p1[a:b + 1] child[np.concatenate([np.arange(b + 1, n), np.arange(a)])] = ( p2[~np.isin(p2, p1[a:b + 1])]) return child def generation(pop: np.ndarray, fit: np.ndarray, dist: np.ndarray, rng: np.random.Generator, elite: int = 2, p_mut: float = 0.9) -> tuple[np.ndarray, np.ndarray]: """One elitist generation: tournament-3 parents, OX, inversion mutation.""" P, n = pop.shape new = np.empty_like(pop) new[:elite] = pop[np.argsort(fit)[:elite]] cand = rng.integers(0, P, size=(2, P - elite, 3)) winners = np.take_along_axis( cand, np.argmin(fit[cand], axis=2)[..., None], axis=2)[..., 0] for i in range(P - elite): child = order_crossover(pop[winners[0, i]], pop[winners[1, i]], rng) if rng.random() < p_mut: a, b = np.sort(rng.choice(n, size=2, replace=False)) child[a:b + 1] = child[a:b + 1][::-1] # inversion mutation new[elite + i] = child return new, tour_lengths(new, dist) def evolve_segment(args: tuple[np.ndarray, np.ndarray, np.ndarray, np.random.SeedSequence, int] ) -> tuple[np.ndarray, np.ndarray]: """Advance one island by `gens` generations (module-level so it pickles).""" pop, fit, dist, seed, gens = args rng = np.random.default_rng(seed) for _ in range(gens): pop, fit = generation(pop, fit, dist, rng) return pop, fit def ring_migration(pops: list[np.ndarray], fits: list[np.ndarray], n_migrants: int) -> None: """Island i-1 sends copies of its best to island i, replacing i's worst.""" k = len(pops) sel = [np.argpartition(f, n_migrants)[:n_migrants] for f in fits] movers = [(pops[i][sel[i]].copy(), fits[i][sel[i]].copy()) for i in range(k)] for i in range(k): sol, val = movers[(i - 1) % k] worst = np.argpartition(fits[i], -n_migrants)[-n_migrants:] pops[i][worst] = sol fits[i][worst] = val def island_ga_tsp(dist: np.ndarray, n_islands: int = 4, pop_size: int = 80, n_segments: int = 12, seg_gens: int = 25, n_migrants: int = 2, seed: int = 0, workers: int | None = None ) -> tuple[np.ndarray, float]: """Segment-synchronous island GA: evolve islands in parallel, migrate, repeat.""" n = dist.shape[0] ss = np.random.SeedSequence(seed) rng = np.random.default_rng(ss) pops = [np.array([nearest_neighbor(dist, int(rng.integers(n))) if j < 4 else rng.permutation(n) # 4 warm starts/island for j in range(pop_size)]) for _ in range(n_islands)] fits = [tour_lengths(p, dist) for p in pops] with ProcessPoolExecutor(max_workers=workers) as ex: for _ in range(n_segments): seeds = ss.spawn(n_islands) # fresh child streams per segment jobs = [(pops[i], fits[i], dist, seeds[i], seg_gens) for i in range(n_islands)] pops, fits = map(list, zip(*ex.map(evolve_segment, jobs))) ring_migration(pops, fits, n_migrants) i = int(np.argmin([f.min() for f in fits])) j = int(np.argmin(fits[i])) return pops[i][j].copy(), float(fits[i][j]) if __name__ == "__main__": rng = np.random.default_rng(42) pts = rng.random((40, 2)) dist = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2) tour, length = island_ga_tsp(dist, seed=7) assert sorted(tour.tolist()) == list(range(40)) # valid permutation print(round(length, 3)) # Expected: ~5.2-5.6 (seed 7 gives 5.437): the GA improves the best # nearest-neighbor start (5.71) but, without local search, will not close # the gap to the ~4.5-5.0 2-opt level — hybridize (memetic-algorithms) for # that. Wall-clock drops near-linearly up to n_islands workers. ``` Two implementation points matter. First, the per-segment `ss.spawn(n_islands)` gives every island a fresh, statistically independent stream each segment without shipping generator state back from the workers — simple and reproducible. Second, `dist` is re-pickled into every task; for instances above a few MB, move it into the workers once via the pool `initializer` pattern shown in the master-slave section. ## Worked Example: Parallel Multistart ILS with multiprocessing Independent multistart is the correct first parallelization for single-solution methods. Each worker runs a full ILS (2-opt descent plus double-bridge kicks; see local-search design notes in metaheuristic-design-principles and the ILS loop in its own skill) from its own seed; the master takes the minimum. No communication, perfect load-balance granularity at the run level, and the exponential time-to-target argument above predicts near-linear speedup to a fixed quality target. ```python """Parallel multistart ILS for the symmetric TSP with a process pool.""" from concurrent.futures import ProcessPoolExecutor import numpy as np def tour_length(tour: np.ndarray, dist: np.ndarray) -> float: """Length of one closed tour.""" return float(dist[tour, np.roll(tour, -1)].sum()) def two_opt(tour: np.ndarray, dist: np.ndarray) -> np.ndarray: """2-opt descent; the delta scan over all j for fixed i is vectorized.""" n = tour.shape[0] improved = True while improved: improved = False for i in range(n - 2): a, b = tour[i], tour[i + 1] j = np.arange(i + 2, n) c, d = tour[j], tour[(j + 1) % n] delta = dist[a, c] + dist[b, d] - dist[a, b] - dist[c, d] k = int(np.argmin(delta)) if delta[k] < -1e-10: jj = i + 2 + k tour[i + 1: jj + 1] = tour[i + 1: jj + 1][::-1] improved = True return tour def double_bridge(tour: np.ndarray, rng: np.random.Generator) -> np.ndarray: """Classic 4-opt kick: cut into four segments A|B|C|D, reorder to A C B D.""" n = tour.shape[0] a, b, c = np.sort(rng.choice(np.arange(1, n), size=3, replace=False)) return np.concatenate([tour[:a], tour[b:c], tour[a:b], tour[c:]]) def ils_run(args: tuple[np.ndarray, np.random.SeedSequence, int] ) -> tuple[float, np.ndarray]: """One independent ILS run (module-level so it pickles under spawn).""" dist, seed, n_kicks = args rng = np.random.default_rng(seed) cur = two_opt(rng.permutation(dist.shape[0]), dist) cur_len = tour_length(cur, dist) best, best_len = cur.copy(), cur_len for _ in range(n_kicks): cand = two_opt(double_bridge(cur, rng), dist) cand_len = tour_length(cand, dist) if cand_len < cur_len - 1e-10: # better-only acceptance cur, cur_len = cand, cand_len if cur_len < best_len: best, best_len = cur.copy(), cur_len return best_len, best def parallel_multistart_ils(dist: np.ndarray, n_starts: int, n_kicks: int, seed: int, workers: int | None = None ) -> tuple[np.ndarray, float]: """Run n_starts independent ILS in separate processes; return the best.""" seeds = np.random.SeedSequence(seed).spawn(n_starts) jobs = [(dist, s, n_kicks) for s in seeds] with ProcessPoolExecutor(max_workers=workers) as ex: results = list(ex.map(ils_run, jobs)) best_len, best = min(results, key=lambda r: r[0]) return best, best_len if __name__ == "__main__": rng = np.random.default_rng(11) pts = rng.random((60, 2)) dist = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2) best, best_len = parallel_multistart_ils(dist, n_starts=8, n_kicks=150, seed=0) print(round(best_len, 3)) # Expected: 6.305 — on this instance all 8 starts converge to the same # (almost surely optimal) tour, while plain 2-opt local optima average ~6.9. # With 8 cores, wall-clock is close to the cost of ONE run, not eight. ``` Practical notes for this pattern: - **Seeds:** never pass `seed + i` by hand; `SeedSequence.spawn` guarantees independent streams even across thousands of workers and avoids accidental stream overlap. - **Windows/macOS spawn:** the worker (`ils_run`) and everything it calls must be importable at module top level, and the driver must sit under `if __name__ == "__main__":` — otherwise each worker re-executes the script and the pool recurses. - **Reporting:** return per-start results (not just the min) so the run table records the distribution across starts; one row per (instance, seed, start) is the tidy format. ## Master-Slave Evaluation, Portfolios, and Cooperative Search ### Master-slave fitness evaluation The master runs the (unchanged) metaheuristic loop; workers evaluate candidates. Use it only when one evaluation is expensive — for cheap objectives, vectorized batch evaluation in one process wins (see fitness-evaluation-and-caching for delta evaluation and memoization, which often remove the bottleneck entirely). Two details carry the pattern: push large read-only data into workers **once** via `initializer`, and keep the executor alive across generations instead of rebuilding it. ```python """Master-slave parallel evaluation: workers hold read-only data, master evolves.""" from concurrent.futures import ProcessPoolExecutor import numpy as np _DATA: dict[str, np.ndarray] = {} def init_worker(dist: np.ndarray) -> None: """Run once per worker process: stash large read-only arrays in a global.""" _DATA["dist"] = dist def eval_one(tour: np.ndarray) -> float: """Evaluation stand-in: tour length computed the slow way on purpose. Replace the body with the real expensive objective (simulation, ML model, embedded LP). A process pool only pays off when one call costs >= ~1 ms. """ dist = _DATA["dist"] total = 0.0 for a, b in zip(tour, np.roll(tour, -1)): total += float(dist[a, b]) return total def evaluate_parallel(ex: ProcessPoolExecutor, pop: np.ndarray, chunksize: int = 8) -> np.ndarray: """Order-preserving parallel batch evaluation of a population.""" return np.fromiter(ex.map(eval_one, pop, chunksize=chunksize), dtype=float, count=pop.shape[0]) if __name__ == "__main__": rng = np.random.default_rng(2) pts = rng.random((50, 2)) dist = np.linalg.norm(pts[:, None, :] - pts[None, :, :], axis=2) pop = np.array([rng.permutation(50) for _ in range(64)]) with ProcessPoolExecutor(max_workers=4, initializer=init_worker, initargs=(dist,)) as ex: fit = evaluate_parallel(ex, pop) # reuse `ex` across generations serial = np.array([float(dist[t, np.roll(t, -1)].sum()) for t in pop]) print(bool(np.allclose(fit, serial))) # Expected: True — parallel evaluation must be a pure speed change, never a # semantic one. Verify this equivalence in a test before trusting speedups. ``` `chunksize` batches tasks to amortize inter-process round trips: with microsecond tasks even `chunksize=256` cannot save the pattern (vectorize instead); with multi-second tasks use `chunksize=1` so slow tasks do not strand fast workers. joblib wraps the same idea with worker reuse and automatic batching; its `loky` backend survives interactive sessions and reuses the pool between calls: ```python """Master-slave evaluation with joblib: reusable pool, automatic batching.""" import numpy as np from joblib import Parallel, delayed def eval_one(x: np.ndarray) -> float: """Costly black-box stand-in: a sequential recursion over the vector.""" state = 0.0 for v in x: state = 0.8 * state + float(np.tanh(v)) return state if __name__ == "__main__": rng = np.random.default_rng(0) pop = rng.normal(size=(200, 64)) with Parallel(n_jobs=-1, batch_size=16) as pool: # keep workers warm fit = np.asarray(pool(delayed(eval_one)(x) for x in pop)) print(fit.shape) # Expected: (200,) — same values and order as the serial list comprehension # [eval_one(x) for x in pop]; speedup appears only for genuinely costly calls. ``` ### Algorithm portfolios A portfolio races entries (algorithm, configuration, seed) on the same instance and keeps the best. It exploits run-time and run-quality variance instead of fighting it, and it removes the "which configuration do I trust" decision when tuning budget is short. Entries finish at different times, so collect with `as_completed`; note that Python cannot kill a *running* future — enforce budgets inside the workers. ```python """Parameter portfolio: race SA configurations on min x^T Q x, x binary.""" from concurrent.futures import ProcessPoolExecutor, as_completed import numpy as np def annealer(Q: np.ndarray, seed: np.random.SeedSequence, n_moves: int, t0_scale: float, cooling: float) -> tuple[float, np.ndarray]: """Single-flip simulated annealing with O(n) move deltas; one portfolio entry.""" rng = np.random.default_rng(seed) n = Q.shape[0] diag = np.diag(Q).copy() x = rng.integers(0, 2, size=n) val = float(x @ Q @ x) best_val, best_x = val, x.copy() temp = t0_scale * float(np.abs(Q).sum()) / n for _ in range(n_moves): i = int(rng.integers(n)) delta = (1 - 2 * x[i]) * (float((Q[i] + Q[:, i]) @ x) - 2 * diag[i] * x[i] + diag[i]) if delta < 0 or rng.random() < np.exp(-delta / max(temp, 1e-12)): x[i] ^= 1 val += delta if val < best_val: best_val, best_x = val, x.copy() temp *= cooling return best_val, best_x def run_entry(args: tuple[str, np.ndarray, np.random.SeedSequence, int, float, float]) -> tuple[str, float, np.ndarray]: """Portfolio worker: one labeled entry (module-level so it pickles).""" label, Q, seed, n_moves, t0_scale, cooling = args val, x = annealer(Q, seed, n_moves, t0_scale, cooling) return label, val, x def race(Q: np.ndarray, n_moves: int, seed: int, workers: int | None = None) -> list[tuple[str, float]]: """Race four SA configurations in parallel; return (label, value) sorted.""" configs = [("hot-slow", 2.0, 0.9999), ("hot-fast", 2.0, 0.999), ("cold-slow", 0.3, 0.9999), ("cold-fast", 0.3, 0.999)] seeds = np.random.SeedSequence(seed).spawn(len(configs)) results: list[tuple[str, float]] = [] with ProcessPoolExecutor(max_workers=workers) as ex: futs = [ex.submit(run_entry, (lab, Q, sd, n_moves, t0, c)) for (lab, t0, c), sd in zip(configs, seeds)] for fut in as_completed(futs): # report entries as they finish label, val, _ = fut.result() results.append((label, val)) return sorted(results, key=lambda r: r[1]) if __name__ == "__main__": rng = np.random.default_rng(5) Q = rng.normal(size=(60, 60)) table = race(Q, n_moves=120_000, seed=1) print(table[0][0], round(table[0][1], 2)) # Expected: "hot-slow -155.11" for these seeds (deterministic: the sort # removes completion-order effects). WHICH configuration wins varies with # seed and instance — exactly the variance a portfolio is built to exploit. ``` ### Cooperative search Cooperation adds communication to multistart: workers periodically deposit their best solutions into a shared elite pool and restart from (perturbed) pool members. This is the pool-and-restart form of adaptive memory (Rochat & Taillard 1995); it intensifies around the globally best material while the perturbation keeps workers from collapsing onto one basin. The segment-synchronous skeleton avoids shared-memory machinery entirely. ```python """Cooperative multistart: independent SA segments plus a shared elite pool.""" from concurrent.futures import ProcessPoolExecutor import numpy as np def sa_segment(args: tuple[np.ndarray, np.ndarray, np.random.SeedSequence, int, float]) -> tuple[float, np.ndarray]: """One SA segment on min x^T Q x from a given start (module-level, pickles).""" Q, x0, seed, n_moves, temp = args rng = np.random.default_rng(seed) diag = np.diag(Q).copy() x = x0.copy() val = float(x @ Q @ x) best_val, best_x = val, x.copy() for _ in range(n_moves): i = int(rng.integers(Q.shape[0])) delta = (1 - 2 * x[i]) * (float((Q[i] + Q[:, i]) @ x) - 2 * diag[i] * x[i] + diag[i]) if delta < 0 or rng.random() < np.exp(-delta / max(temp, 1e-12)): x[i] ^= 1 val += delta if val < best_val: best_val, best_x = val, x.copy() temp *= 0.999 return best_val, best_x def cooperative_search(Q: np.ndarray, n_workers: int = 4, n_rounds: int = 5, n_moves: int = 20_000, seed: int = 0 ) -> tuple[np.ndarray, float]: """Pool-and-restart cooperation: share the incumbent between SA rounds.""" ss = np.random.SeedSequence(seed) rng = np.random.default_rng(ss) n = Q.shape[0] starts = [rng.integers(0, 2, size=n) for _ in range(n_workers)] elite_x, elite_val = None, np.inf temp = float(np.abs(Q).sum()) / n with ProcessPoolExecutor(max_workers=n_workers) as ex: for _ in range(n_rounds): seeds = ss.spawn(n_workers) jobs = [(Q, starts[w], seeds[w], n_moves, temp) for w in range(n_workers)] for val, x in ex.map(sa_segment, jobs): if val < elite_val: elite_val, elite_x = val, x.copy() flips = rng.random((n_workers, n)) < 0.05 # diversified restarts starts = [np.where(flips[w], 1 - elite_x, elite_x) for w in range(n_workers)] temp *= 0.5 # cool across rounds return elite_x, float(elite_val) if __name__ == "__main__": rng = np.random.default_rng(9) Q = rng.normal(size=(80, 80)) x, val = cooperative_search(Q) print(round(val, 2)) # Expected: -294.42 for these seeds (the segment pattern is deterministic). # Benchmark requirement: at equal total moves, this should beat or match # n_workers fully independent SA runs on most seeds — verify, don't assume. ``` The 5% restart-flip rate is the cooperative analog of ILS kick strength: too small and all workers intensify the same basin; too large and cooperation degenerates to independent multistart. Tune it exactly like a perturbation strength. ## Advanced Techniques ### Asynchronous islands and topology effects Synchronous segments waste time at the barrier when islands run at different speeds. Asynchronous designs let each island push/pull migrants through queues (`multiprocessing.Queue`, ray actors, or MPI nonblocking sends) without a global generation clock. The cost is reproducibility: results depend on OS scheduling, so log every migrant event if you go asynchronous. Topology is a second lever on the same diversity dial: takeover time — the generations until one genotype occupies everything — grows from complete graphs (fastest takeover, least diversity) through random k-regular to rings (slowest takeover, most diversity); see Cantu-Paz (2000). When migration is rare and the topology sparse, an island model approximates independent multistart; when migration is frequent and dense, it approximates one big panmictic population. Place your design deliberately between those poles. ### Heterogeneous islands Nothing forces islands to share parameters or even operators: give each island a different mutation rate, crossover, or acceptance temperature, and migration lets the configurations rescue each other on instances where one setting stalls. This is a portfolio folded into an island model and is the cheapest insurance against mis-tuning — the spread of parameters substitutes for a tuning study when the budget is small. ### Cooperative-search design rules Cooperation changes the search trajectory; treat the communication policy as part of the algorithm, not as plumbing (Toulouse, Crainic & Gendreau 1996 showed that careless sharing can make cooperating searches *worse* than independent ones). Three rules of thumb: share **rarely** (every segment, not every iteration), share **selectively** (a bounded elite pool with duplicate elimination, not every improvement), and **diversify on receipt** (perturb what you adopt, as in the restart-flip above, so receivers do not pile onto the sender's basin). Bound the pool (10-50 entries) and evict by worst-objective-first with a distance filter so the memory stays diverse. ### Measuring parallel performance honestly Report speedup against the **best serial implementation** (the vectorized one), not against a deliberately slow baseline — this single rule eliminates most inflated claims. For quality-targeting designs, report time-to-target plots (Aiex, Resende & Ribeiro 2007, "TTT plots") and quality-vs-wall-clock curves over 10+ seeds rather than a single speedup number. Super-linear speedup is possible for island models because 8 islands of 50 is a *different algorithm* than one population of 400 (Alba 2002, "Parallel evolutionary algorithms can achieve super-linear performance") — when you observe it, say that the search dynamics changed; do not present it as a hardware miracle. ### Scaling beyond one machine Cross the machine boundary only when a profiled single-machine setup is insufficient. ray maps cleanly onto this skill's patterns: actors for asynchronous islands holding their own state, `ray.put` for sharing large instance data once per node, tasks for master-slave evaluation. mpi4py is the standard on academic clusters (synchronous segment exchange via `allgather`, asynchronous migration via `isend/irecv`). Keep the algorithm code identical and swap only the executor layer; if your code is structured like the examples above (module-level pure workers + a thin driver), the port is mechanical. ## Practical Challenges **"Can't pickle local object" or workers crash immediately on Windows.** Spawn-based start (Windows, macOS default) re-imports the main module in each worker: every worker function must live at module top level (no lambdas, no closures, no methods of unpicklable objects), and the driver must be guarded by `if __name__ == "__main__":`. Without the guard each worker re-runs the script and recursively spawns pools. On Linux the same code "works" under fork — test on the deployment OS, not just the development one. **The parallel version is slower than the serial one.** Process startup (hundreds of ms), per-task pickling, and result transfer are fixed costs; they dominate when tasks are tiny. Fix the granularity: batch work with `chunksize` / `batch_size`, move whole runs or segments (not single evaluations) into tasks, ship large read-only data once via `initializer` instead of per task, and reuse one executor for the whole run. If a task still finishes in under a millisecond, the correct fix is rung 1: vectorize, don't parallelize. **All workers produce identical results.** Classic seeding bug: the same seed (or the same global RNG state under fork) reaches every worker. Spawn child seeds with `np.random.SeedSequence(seed).spawn(k)` and build a fresh `default_rng` inside each worker from its own child. Never reseed from `time()` in workers — simultaneous starts collide. **The machine is oversubscribed and everything crawls.** numpy's BLAS may use all cores for matrix work; p worker processes times c BLAS threads equals p*c runnable threads thrashing the caches. In workers, pin BLAS to one thread (set `OMP_NUM_THREADS=1` / `MKL_NUM_THREADS=1` before numpy import, or use `threadpoolctl.threadpool_limits(1)` in the initializer) and give the pool one process per physical core. **Migration destroyed diversity and all islands converged to the same tour.** Migration is a selection operator; best-replace-worst on a complete topology with many migrants causes rapid takeover. Lengthen the interval, cut migrants to 1-2, switch to a ring, or migrate random individuals instead of the best. Diagnose by tracking per-island best and inter-island distance: if islands agree early, migration pressure is too high. **Results are not reproducible run-to-run even with fixed seeds.** Asynchronous collection (`as_completed`), unordered `imap_unordered`, …(truncated)