Memetic Algorithms

Memetic Algorithms

You are an expert in memetic algorithms (MAs): population-based metaheuristics that hybridize a genetic algorithm with local search so that every individual in the population is a local optimum (or near one). This skill covers the canonical MA loop, Lamarckian vs Baldwinian learning, budgeting local-search frequency and depth, restart management, and keeping a population diverse when strong local search keeps collapsing it. Use the framework below to design, implement, and tune an MA, with complete worked implementations for the quadratic assignment problem (QAP) and the traveling salesman problem (TSP).

Initial Assessment

Establish the following before designing, implementing, or debugging a memetic algorithm:

• Problem class and encoding. Permutation (TSP, QAP, flow shop), binary selection (knapsack, set covering), integer assignment, or decoder-based? The encoding fixes which crossover, mutation, and neighborhood moves are legal. For representation choice see solution-encodings; for operator catalogs see crossover-operators and mutation-and-perturbation-operators.
• Local-search ingredients. Does a neighborhood with delta (incremental) evaluation already exist? What is the cost of one full descent in evaluations and milliseconds? An MA without fast delta evaluation is usually a mistake — fix that first (see local-search-and-neighborhoods and fitness-evaluation-and-caching).
• Evaluation budget. Wall-clock limit, evaluation-count limit, or both? Expect local search to consume 90%+ of all evaluations; the budget split between evolution and learning is the central design decision.
• Lamarckian feasibility. Can an improved phenotype be written back into the genotype? Direct encodings: yes. Decoder-based encodings (random keys, priority rules): often no — the improved schedule may have no preimage, which forces Baldwinian or repair-style designs.
• Constraint handling. Are all neighborhood moves feasibility-preserving, or do you need repair after crossover/mutation? Decide where infeasibility is allowed to exist (never, only pre-repair, or penalized throughout).
• Quality requirement. Gap to best-known values on benchmarks, or "good solution in 5 minutes"? MAs are the state of the art for QAP and TSP benchmarks but are heavier machinery than iterated local search (ILS).
• Baselines. Has anyone run plain GA, multistart local search, or ILS on this problem with the same budget? An MA must beat its own components, or the hybrid is not earning its complexity.
• Instance scale and protocol. Sizes (n), number of instances, tuning/test split, number of seeds per configuration.
• Reproducibility. Single np.random.default_rng(seed) per run; seed recorded with every result row.
• Compute model. Pure numpy on one core, multiprocessing for parallel descents, or batch (vectorized) fitness evaluation?

Algorithm Anatomy

A memetic algorithm (Moscato 1989, "On Evolution, Search, Optimization, Genetic Algorithms and Martial Arts: Towards Memetic Algorithms") layers individual learning — local search — on top of populational evolution — selection, crossover, mutation. The term "cultural algorithm" is often used loosely for the same hybrid; strictly, cultural algorithms (Reynolds 1994) add a shared belief space that biases variation, but the practical design questions (who learns, how much, and what gets inherited) are identical, and the framework below covers both readings.

The local-optimum subspace. For minimization of $f$ over a finite set $X$ with neighborhood $N$, define the set of local optima

$$L_N = \{\, x \in X : f(x) \le f(y)\ \ \forall y \in N(x) \,\}.$$

After Lamarckian local search, every population member lies in $LN$, so the MA effectively evolves over $LN$ — typically orders of magnitude smaller than $X$ and, on many problems, structured: good local optima share components (a "big valley", Boese 1995 for TSP; high fitness-distance correlation for QAP in Merz & Freisleben 2000). Crossover exploits exactly that sharing: it recombines components of two local optima, mutation supplies escape energy, and local search projects the child back onto $L_N$:

$$x{\text{child}} = \mathrm{LS}\big(\mathrm{mut}(\mathrm{cx}(x{p1}, x{p2}))\big) \in LN.$$

Learning models. What happens after local search improves a child from $x$ to $x^\* = \mathrm{LS}(x)$:

| Model | Genotype stored | Fitness used by selection | Use when | |---|---|---|---| | Lamarckian | $x^\$ (writeback) | $f(x^\)$ | Direct encodings. Default for combinatorial optimization. | | Baldwinian | $x$ (unchanged) | $f(x^\)$ | Writeback impossible (decoder-based encodings) or genotypic diversity must be protected. | | Partial Lamarckian | $x^\$ with probability $q$ | $f(x^\*)$ | Compromise; $q \approx 0.5$ can beat both extremes on rugged landscapes (Houck, Joines, Kay & Wilson 1997). |

Whitley, Gordon & Mathias (1994, "Lamarckian evolution, the Baldwin effect and function optimization") showed Baldwinian search can win on deceptive landscapes, but for the combinatorial problems in this repository Lamarckian writeback is almost always faster and is the recommended default.

Budget accounting. With $\lambda$ offspring per generation, $G$ generations, local-search probability $p{ls}$, and $\bar e{ls}$ average (delta-)evaluations per descent:

$$E{\text{total}} \;=\; \underbrace{\lambda G}{\text{children}} \;+\; \underbrace{p{ls}\,\lambda G\,\bar e{ls}}_{\text{local search}} .$$

In practice the second term dominates — well over 90% of evaluations happen inside descents. MA design is therefore mostly the art of spending the local-search budget well: which children learn ($p_{ls}$ or a selection rule), and how deeply (ls_max_moves).

When to use an MA versus alternatives:

| Situation | Recommendation | |---|---| | Fast local search with delta evaluation exists, and good solutions share components | MA. Crossover recombines building blocks across basins (shared edges in TSP, shared assignments in QAP). | | Solutions do not share exploitable structure, or no respectful crossover exists | ILS or tabu search — a population adds bookkeeping but little search power. | | Objective is expensive and has no delta evaluation | Plain GA with surrogate evaluation, or an MA with highly selective, shallow local search. | | Feasibility is hard to maintain under moves | Decoder- or repair-based GA first; add local search restricted to feasible moves later. | | Best-known-quality results needed on classic benchmarks | MA is the state of the art for QAP (Merz & Freisleben 2000) and TSP (Nagata & Kobayashi 2013, the EAX genetic algorithm). |

Complexity per generation. $O(\lambda \cdot c{\text{var}} + p{ls}\lambda \cdot \bar e{ls} \cdot c{\text{move}})$ where $c{\text{var}}$ is crossover+mutation cost ($O(n)$ for the permutation operators used below) and $c{\text{move}}$ is the cost of one neighborhood move evaluation — $O(n)$ for a QAP swap delta, $O(1)$ for a TSP 2-opt delta. Population sizes are deliberately small (10–50): each member is an expensive, high-quality local optimum, not a cheap random sample.

Generic Memetic Framework

The skeleton every MA below instantiates:

MEMETIC_ALGORITHM(f, config):
    P ← INIT_POPULATION(pop_size)             # random or construction-heuristic seeds
    for x in P:
        x ← LOCAL_SEARCH(x, ls_max_moves)     # start the evolution from local optima
    best ← argmin_{x in P} f(x)
    while budget remains:
        O ← ∅
        while |O|  float:
    """Mean pairwise Hamming distance between rows, normalized to [0, 1]."""
    m, n = pop.shape
    diff = pop[:, None, :] != pop[None, :, :]
    return float(diff.sum() / (m * (m - 1) * n))

def memetic_algorithm(
    init_population: Callable[[np.random.Generator], np.ndarray],
    evaluate: Callable[[np.ndarray], float],
    crossover: Callable[[np.ndarray, np.ndarray, np.random.Generator], np.ndarray],
    mutate: Callable[[np.ndarray, np.random.Generator], np.ndarray],
    local_search: Callable[[np.ndarray, float, int], tuple[np.ndarray, float]],
    config: MAConfig,
) -> tuple[np.ndarray, float, list[float]]:
    """Generic memetic algorithm for minimization over array-encoded solutions.

    `init_population(rng)` must return an int array of shape (pop_size, n).
    `local_search(x, f_x, max_moves)` must return an improved pair (x', f').
    Lamarckian mode stores x' in the population; Baldwinian mode keeps x but
    lets f' drive selection and replacement.
    """
    rng = np.random.default_rng(config.seed)

    def seeded_population(fresh: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, float]:
        """Local-search every individual; also return the best improved pair."""
        fits = np.empty(len(fresh), dtype=float)
        top_x, top_f = fresh[0].copy(), float("inf")
        for i in range(len(fresh)):
            x1, f1 = local_search(fresh[i], evaluate(fresh[i]), config.ls_max_moves)
            if config.lamarckian:
                fresh[i] = x1
            fits[i] = f1
            if f1  np.ndarray:
    """Random permutation population of shape (POP, N)."""
    return np.array([rng.permutation(N) for _ in range(POP)])

def displacement(x: np.ndarray) -> float:
    """Toy objective: sum_i |x[i] - i| (0 exactly at the identity)."""
    return float(np.abs(x - np.arange(len(x))).sum())

def pmx(p1: np.ndarray, p2: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Partially mapped crossover: copy a slice of p1, repair conflicts by mapping."""
    n = len(p1)
    a, b = np.sort(rng.choice(n, size=2, replace=False))
    child = p2.copy()
    child[a:b + 1] = p1[a:b + 1]
    mapping = dict(zip(p1[a:b + 1].tolist(), p2[a:b + 1].tolist()))
    for i in list(range(a)) + list(range(b + 1, n)):
        v = int(child[i])
        while v in mapping:
            v = mapping[v]
        child[i] = v
    return child

def swap_mutation(x: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """Exchange two random positions."""
    y = x.copy()
    i, j = rng.choice(len(y), size=2, replace=False)
    y[[i, j]] = y[[j, i]]
    return y

def swap_descent(x: np.ndarray, f_x: float, max_moves: int) -> tuple[np.ndarray, float]:
    """First-improvement swap descent with O(1) move deltas."""
    x, f, n, moves = x.copy(), f_x, len(x), 0
    improved = True
    while improved and moves  tuple[np.ndarray, np.ndarray]:
    """Random symmetric QAP: integer flows F and rounded Euclidean distances D."""
    rng = np.random.default_rng(seed)
    F = np.triu(rng.integers(1, 10, size=(n, n)), 1)
    F = F + F.T
    pts = rng.uniform(0.0, 100.0, size=(n, 2))
    D = np.rint(np.hypot(pts[:, None, 0] - pts[None, :, 0],
                         pts[:, None, 1] - pts[None, :, 1])).astype(int)
    return F, D

def qap_cost(p: np.ndarray, F: np.ndarray, D: np.ndarray) -> float:
    """Objective sum_ij F[i, j] * D[p[i], p[j]] (facility i sits at location p[i])."""
    return float((F * D[np.ix_(p, p)]).sum())

def swap_delta(p: np.ndarray, r: int, s: int, F: np.ndarray, D: np.ndarray) -> float:
    """O(n) cost change of swapping the locations of facilities r and s.

    Valid for symmetric F and D with zero diagonals (Taillard 1991,
    'Robust taboo search for the quadratic assignment problem').
    """
    k = np.delete(np.arange(len(p)), [r, s])
    return float(2.0 * ((F[r, k] - F[s, k]) * (D[p[s], p[k]] - D[p[r], p[k]])).sum())

def two_exchange_ls(p: np.ndarray, f_p: float, max_moves: int,
                    F: np.ndarray, D: np.ndarray) -> tuple[np.ndarray, float]:
    """First-improvement 2-exchange descent using O(n) swap deltas."""
    p, f, n, moves = p.copy(), f_p, len(p), 0
    improved = True
    while improved and moves  np.ndarray:
    """CX: decompose positions into cycles, alternate the source parent per cycle."""
    n = len(p1)
    child = np.full(n, -1)
    pos_in_p1 = np.empty(n, dtype=int)
    pos_in_p1[p1] = np.arange(n)
    take_p1 = bool(rng.integers(0, 2))
    assigned = np.zeros(n, dtype=bool)
    for start in range(n):
        if assigned[start]:
            continue
        i, cycle = start, []
        while not assigned[i]:
            assigned[i] = True
            cycle.append(i)
            i = int(pos_in_p1[p2[i]])
        src = p1 if take_p1 else p2
        child[cycle] = src[cycle]
        take_p1 = not take_p1
    return child

def qap_memetic(F: np.ndarray, D: np.ndarray, pop_size: int = 20,
                n_generations: int = 40, ls_max_moves: int = 5000,
                mutation_rate: float = 0.3, elite_keep: int = 2,
                seed: int = 0) -> tuple[np.ndarray, float]:
    """Lamarckian memetic algorithm for the symmetric QAP."""
    rng = np.random.default_rng(seed)
    n = F.shape[0]
    pop = np.array([rng.permutation(n) for _ in range(pop_size)])
    fit = np.empty(pop_size)
    for i in range(pop_size):
        pop[i], fit[i] = two_exchange_ls(pop[i], qap_cost(pop[i], F, D),
                                         ls_max_moves, F, D)
    for _ in range(n_generations):
        order = np.argsort(fit)
        new_pop = [pop[i].copy() for i in order[:elite_keep]]
        new_fit = [float(fit[i]) for i in order[:elite_keep]]
        while len(new_pop)  np.ndarray:
    """Random Euclidean TSP: n points in the unit square -> distance matrix."""
    rng = np.random.default_rng(seed)
    pts = rng.uniform(0.0, 1.0, size=(n, 2))
    return np.hypot(pts[:, None, 0] - pts[None, :, 0],
                    pts[:, None, 1] - pts[None, :, 1])

def tour_length(tour: np.ndarray, D: np.ndarray) -> float:
    """Total length of the closed tour."""
    return float(D[tour, np.roll(tour, -1)].sum())

def two_opt_ls(tour: np.ndarray, f_tour: float, max_moves: int,
               D: np.ndarray) -> tuple[np.ndarray, float]:
    """2-opt descent; for each first edge, the second edge is scanned vectorized.

    Removing edges (a,b) and (c,d) and reconnecting as (a,c),(b,d) reverses the
    segment between them; the delta D[a,c] + D[b,d] - D[a,b] - D[c,d] is O(1)
    per candidate and computed for all candidates j at once.
    """
    tour, f, n, moves = tour.copy(), f_tour, len(tour), 0
    improved = True
    while improved and moves  0 else n - 2   # skip the closing-edge pairing at i == 0
            j = np.arange(i + 2, j_hi + 1)
            if len(j) == 0:
                continue
            c, d = tour[j], tour[(j + 1) % n]
            deltas = D[a, c] + D[b, d] - D[a, b] - D[c, d]
            k = int(np.argmin(deltas))
            if deltas[k]  np.ndarray:
    """OX: copy a slice of p1, fill remaining cities in p2's cyclic order."""
    n = len(p1)
    a, b = np.sort(rng.choice(n, size=2, replace=False))
    child = np.empty(n, dtype=int)
    child[a:b + 1] = p1[a:b + 1]
    used = np.zeros(n, dtype=bool)
    used[p1[a:b + 1]] = True
    order = np.concatenate([p2[b + 1:], p2[:b + 1]])
    child[np.concatenate([np.arange(b + 1, n), np.arange(a)])] = order[~used[order]]
    return child

def double_bridge(tour: np.ndarray, rng: np.random.Generator) -> np.ndarray:
    """4-opt double bridge: cut the tour at three points and reorder the segments."""
    n = len(tour)
    i, j, k = np.sort(rng.choice(np.arange(1, n), size=3, replace=False))
    return np.concatenate([tour[:i], tour[j:k], tour[i:j], tour[k:]])

def tsp_memetic(D: np.ndarray, pop_size: int = 20, n_generations: int = 30,
                ls_max_moves: int = 20_000, mutation_rate: float = 0.4,
                elite_keep: int = 2, seed: int = 0) -> tuple[np.ndarray, float]:
    """Lamarckian memetic algorithm for the symmetric TSP (2-opt local search)."""
    rng = np.random.default_rng(seed)
    n = D.shape[0]
    pop = np.array([rng.permutation(n) for _ in range(pop_size)])
    fit = np.empty(pop_size)
    for i in range(pop_size):
        pop[i], fit[i] = two_opt_ls(pop[i], tour_length(pop[i], D),

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [hajibabaie](https://github.com/hajibabaie)
- **Source:** [hajibabaie/combinatorial-optimization-skills](https://github.com/hajibabaie/combinatorial-optimization-skills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.