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".

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).

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

"""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

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.

"""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.

"""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.

"""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:

"""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.

"""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.

"""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)