tabu-search

When the user wants to design or implement tabu search — tabu lists and tenure, move attributes, aspiration criteria, frequency-based diversification, candidate list strategies, and elite-restart inte

Files1

View on GitHub

SKILL.md

GitHub

Tabu Search

You are an expert in metaheuristic optimization, specifically in designing and implementing tabu search (Glover 1986, "Future paths for integer programming and links to artificial intelligence"; Glover & Laguna 1997, Tabu Search). This skill covers the memory structures that make tabu search work — tabu lists and tenure, move attributes, aspiration criteria, frequency-based long-term memory, candidate list strategies, and elite-based intensification — plus two reference implementations: Taillard's robust tabu search for the QAP and a critical-path tabu search for job-shop scheduling. Use the framework below to pick each memory component deliberately, implement against the verified code, and validate results against exact baselines on small instances.

Initial Assessment

Establish the following before writing any tabu search code. Each answer changes a design decision downstream.

  • Problem class and representation. Permutation, binary vector, assignment, sequencing per machine? The representation fixes which neighborhoods and which move attributes are available. See solution-representation guidance in local-search-and-neighborhoods before adding memory on top.
  • Neighborhood and its size. Tabu search scans many moves per iteration (often the whole neighborhood). An O(n²) swap neighborhood with O(1) delta lookup is the sweet spot; an O(n²) neighborhood with O(n) evaluation each already costs O(n³) per iteration and needs candidate lists.
  • Delta evaluation cost. Tabu search is only competitive when move evaluation is incremental. If each move needs a full objective recomputation, fix that first (fitness-evaluation-and-caching) — memory structures cannot compensate for a slow scan.
  • What should be forbidden? Decide the move attribute: the exact inverse move, an element-position assignment, a destroyed arc/edge. This is the single most important design choice (see the granularity table below).
  • Tenure regime. Fixed, randomized per move, or reactive? Randomized tenure (Taillard 1991) is the robust default; fixed tenure needs per-instance tuning.
  • Evaluation budget and time limit. Iterations × (scan cost) must fit the budget. Tabu search has no natural stopping point; pick max iterations, max iterations without improvement, or wall-clock time explicitly.
  • Deterministic or stochastic runs? Tabu search with fixed tenure and deterministic tie-breaking is deterministic given the start. With randomized tenure/tie-breaking, plan a multi-seed protocol for reporting.
  • Feasibility handling. Will the search stay inside the feasible region (feasibility-preserving moves), or oscillate across the boundary with penalties? Strategic oscillation needs a penalty schedule.
  • Known optima or best-known values. For QAPLIB, Taillard, OR-Library instances, best-known values exist — report gaps against them. For new problems, build a small-instance exact baseline (MIP/CP) first.
  • Cycling risk indicators. Plateaus and symmetric solutions raise cycling risk; plan solution hashing for cycle detection and a diversification mechanism from the start, not as an afterthought.
  • Single long run vs restarts. Decide whether the budget goes into one long run with diversification phases or several restarts from elite solutions; both need the elite pool machinery below.
  • Comparison baselines. At minimum: best-improvement local search with random restarts, and one alternative metaheuristic (e.g., simulated annealing). A tabu list must demonstrably beat memoryless descent on your instances.

Algorithm Anatomy

Tabu search is best-improvement local search plus memory. At iteration $t$, with current solution $x_t$, neighborhood $\mathcal{N}(x_t)$, active tabu attribute set $T_t$, and incumbent $x^*$:

$$ x_{t+1} = \operatorname*{argmin}_{y \in \mathcal{N}(x_t)} \left{ f(y) ;:; \text{attr}(x_t \to y) \notin T_t ;\text{ or }; f(y) < f(x^*) \right} $$

The move is executed even if $f(x_{t+1}) > f(x_t)$ — tabu search climbs out of local optima deterministically, unlike simulated annealing's probabilistic acceptance. The tabu list then forbids reversing recent moves for tenure $\tau$ iterations, so the search cannot fall straight back into the optimum it just left. The condition $f(y) < f(x^*)$ is the standard improved-best aspiration: a tabu status is overridden whenever the move reaches a new global best, because the tabu list exists to prevent cycling, and a strictly better solution cannot be a repeat.

The three memories

Memory Time scale Stores Drives
Recency (short-term) last $\tau$ iterations move/solution attributes recently changed cycle avoidance — the tabu list itself
Frequency (long-term) whole run transition counts (how often each move fired), residence counts (how long each attribute was present) diversification: penalize over-used moves, force rarely-used assignments
Quality (medium-term) whole run elite solutions and their common features intensification: restart from elites, fix shared attributes

A minimal tabu search uses only recency memory. A production tabu search uses all three: recency every iteration, frequency when the search stagnates, quality when the budget allows restarts. Glover & Laguna (1997) treat this layering as the defining feature of the method.

Move attributes: what exactly becomes tabu

The tabu list stores attributes of executed moves, not full solutions. The attribute granularity sets how aggressively the list prunes:

Attribute recorded A later move is tabu if it ... Restrictiveness Typical use
Full solution hash recreates an exact previous solution weakest — blocks only exact revisits cycle detection, not prevention (hashing via fitness-evaluation-and-caching)
Exact inverse move, e.g. "swap positions (i, j)" applies the inverse move weak — the same solution can be rebuilt by other move sequences very large neighborhoods where stronger attributes over-restrict
Element-position pair, e.g. "unit $r$ returns to location $\ell$" restores a recently destroyed assignment medium — the Ro-TS choice for QAP assignment and permutation problems
Element only, e.g. "element $r$ moves at all" touches a recently moved element strong — freezes elements scheduling moves with many positions per element
Solution feature, e.g. "arc $u \to v$ re-enters the tour/sequence" re-creates a recently removed feature strong, feature-aligned TSP edges, JSP disjunctive arcs

Rule of thumb: pick the attribute that matches what the objective depends on (edges for tours, assignments for QAP, arcs for schedules). Too-fine attributes let the search cycle through near-duplicates; too-coarse attributes forbid large unexplored regions and force constant aspiration overrides.

Tenure

Tenure $\tau$ is how long an attribute stays tabu, in iterations. The classic guidance:

  • Fixed: $\tau = 7$ was Glover's early constant; $\tau \approx \sqrt{|\mathcal{N}|}$ or $\tau \in [0.5\sqrt{n}, 2\sqrt{n}]$ are common starting points. Fixed tenure is fragile — one value rarely fits all instances.
  • Randomized per move: draw $\tau \sim U[\tau_{\min}, \tau_{\max}]$ each time an attribute is recorded. Taillard (1991, "Robust taboo search for the quadratic assignment problem") uses $U[0.9n,,1.1n]$ for QAP and shows the randomization itself prevents cycles that any fixed value permits. This is the recommended default.
  • Reactive: adapt $\tau$ from observed solution repetitions (Battiti & Tecchiolli 1994, "The reactive tabu search") — see Advanced Techniques.

Symptoms: tenure too short → the search returns to the same local optimum every few iterations (detect by hashing solutions); tenure too long → objective trace shows long flat ramps because most of the neighborhood is forbidden, and the all-tabu fallback fires often.

Aspiration criteria

  1. Improved-best (always implement): override tabu if the move yields $f < f(x^*)$.
  2. Default / least-tabu (always implement): if every move is tabu and none aspirates, execute the move whose tabu status expires earliest (or the best tabu move). Never skip an iteration — a frozen search wastes the remaining budget.
  3. Forced long-term aspiration (Ro-TS): execute a move regardless of its value if it creates an assignment not seen for $A$ iterations (e.g., $A = 5n^2$). This is diversification built into the aspiration test.
  4. Regional aspiration: override tabu if the move improves on the best solution within the current search region (requires defining regions; rarely worth it before the first three).

When to choose tabu search

Situation Verdict
Cheap O(1)–O(n) delta evaluation, moderate neighborhood (≤ ~10⁵ moves) Tabu search is a top contender — full-scan best-improvement exploits the delta structure SA wastes
QAP, job-shop, graph coloring (tabucol, Hertz & de Werra 1987), vehicle routing intra-route Tabu search variants are or were state of the art; start here
Very expensive evaluation, no delta structure Prefer methods that evaluate few candidates per iteration (ILS, SA) or fix the evaluation first
Huge neighborhoods (ruin-and-recreate scale) Use large neighborhood search; tabu memory can still guard the repair
Highly constrained, hard-to-stay-feasible problems Consider strategic oscillation or decoder-based methods before plain tabu search

Per-iteration complexity = (number of candidate moves scanned) × (delta lookup cost) + (memory update cost). The QAP implementation below achieves O(n²) per iteration: O(n²) scan with O(1) lookups plus O(n²) delta-matrix maintenance.

Generic Tabu Search Framework

The skeleton every variant shares:

TABU SEARCH (minimization)
Input: initial solution x0, neighborhood N(.), tenure rule tau(.),
       aspiration test A(.), iteration budget T
x <- x0;  x_best <- x0
tabu <- empty map {attribute -> iteration until which it is active}
for t = 1 .. T:
    C <- candidate moves from N(x)              # full scan or candidate list
    admissible <- { m in C : tabu[attr_check(m)] < t  or  A(m) holds }
    if admissible is empty:
        m* <- move in C whose tabu status expires earliest   # default aspiration
    else:
        m* <- argmin_{m in admissible} f(x (+) m)   # best admissible, even uphill
    x <- x (+) m*
    tabu[attr_record(m*)] <- t + tau(t)         # forbid reversal for tau iterations
    if f(x) < f(x_best): x_best <- x
    # optional: update frequency memory; trigger diversification/intensification
return x_best

The generic engine below separates the four problem-specific callables (objective, move enumeration, move execution, tenure rule) from the memory logic. Use it to prototype; the two worked applications then show why production implementations specialize the scan and the delta bookkeeping for speed.

import numpy as np
from dataclasses import dataclass
from collections.abc import Callable, Hashable, Iterable
from typing import Any

@dataclass(frozen=True)
class Move:
    """One candidate move: payload to execute, objective change, tabu attributes."""
    payload: Any
    delta: float
    check_key: Hashable    # move is tabu while this attribute is active
    record_key: Hashable   # attribute activated when the move is executed

def tabu_search(
    initial: Any,
    objective: Callable[[Any], float],
    enumerate_moves: Callable[[Any], Iterable[Move]],
    execute: Callable[[Any, Any], Any],
    tenure: Callable[[np.random.Generator], int],
    n_iters: int,
    seed: int = 0,
) -> tuple[Any, float, list[float]]:
    """Attribute-based best-improvement tabu search (minimization).

    Returns (best solution, best value, best-so-far trace per iteration)."""
    rng = np.random.default_rng(seed)
    current, cur_val = initial, objective(initial)
    best, best_val = current, cur_val
    tabu_until: dict[Hashable, int] = {}
    trace = [best_val]
    for t in range(1, n_iters + 1):
        chosen, chosen_val = None, np.inf
        fallback, fallback_exp = None, np.inf   # least-tabu default aspiration
        for mv in enumerate_moves(current):
            new_val = cur_val + mv.delta
            if tabu_until.get(mv.check_key, 0) >= t and new_val >= best_val:
                if tabu_until[mv.check_key] < fallback_exp:
                    fallback, fallback_exp = mv, tabu_until[mv.check_key]
                continue
            if new_val < chosen_val:
                chosen, chosen_val = mv, new_val
        if chosen is None:                      # every move tabu, none aspired
            if fallback is None:
                break                           # neighborhood empty: stop
            chosen, chosen_val = fallback, cur_val + fallback.delta
        current = execute(current, chosen.payload)
        cur_val = chosen_val
        tabu_until[chosen.record_key] = t + tenure(rng)
        if cur_val < best_val:
            best, best_val = current, cur_val
        trace.append(best_val)
    return best, best_val, trace

# Demo: QUBO  min x^T Q x  with bit-flip moves and O(n) vectorized deltas.
def qubo_flip_deltas(Q: np.ndarray, x: np.ndarray) -> np.ndarray:
    """Objective change of flipping each bit of x (Q symmetric)."""
    s, dg = Q @ x, np.diag(Q)
    return (1 - 2 * x) * (dg + 2 * (s - dg * x))

rng = np.random.default_rng(7)
n = 14
Q = rng.integers(-10, 11, (n, n))
Q = (Q + Q.T) // 2

def moves(x: np.ndarray) -> list[Move]:
    """All bit flips; attribute = flipped index (element attribute)."""
    d = qubo_flip_deltas(Q, x)
    return [Move(payload=i, delta=float(d[i]), check_key=i, record_key=i)
            for i in range(n)]

def do_flip(x: np.ndarray, i: int) -> np.ndarray:
    """Return a copy of x with bit i flipped."""
    y = x.copy()
    y[i] = 1 - y[i]
    return y

x0 = rng.integers(0, 2, n)
best, best_val, trace = tabu_search(
    x0, lambda x: float(x @ Q @ x), moves, do_flip,
    tenure=lambda g: int(g.integers(5, 9)), n_iters=200, seed=1)
bits = ((np.arange(2 ** n)[:, None] >> np.arange(n)) & 1).astype(np.int64)
exact = int(np.einsum("bi,ij,bj->b", bits, Q, bits).min())
print(best_val, exact)
# Expected: tabu search reaches -143.0, equal to the brute-force optimum -143

Parameter guidance

Parameter Typical range Trades off
Tenure $\tau$ $U[0.9n, 1.1n]$ (QAP swaps); 8–14 (JSP arcs); $\sqrt{ \mathcal{N}
Tenure regime randomized > reactive > fixed (in robustness order) randomized needs a range, not a value; reactive adds two meta-parameters
Iteration budget $10^3$–$10^6$; or $50n$–$1000n$ without improvement solution quality vs wall-clock time
Aspiration improved-best + least-tabu fallback always; forced after $A = 5n^2$ idle larger $A$ = less forced diversification
Candidate list size 5–20 % of $\mathcal{N}$ when full scan is too slow scan speed vs move quality per iteration
Diversification trigger $0.5n$–$5n$ iterations without incumbent improvement too eager destroys convergence; too lazy wastes budget
Elite pool size 5–20 distinct solutions memory + restart variety vs dilution with mediocre elites
Frequency penalty weight scale to the median |delta| of accepted moves too small = no effect; too large = frequency dictates the search

Operator-level details (swap vs insertion neighborhoods, delta derivations per move type, perturbation design) belong to local-search-and-neighborhoods and fitness-evaluation-and-caching — this skill assumes a working neighborhood and adds memory on top.

Worked Application 1: Robust Tabu Search for the QAP

The quadratic assignment problem: given flow matrix $F \in \mathbb{Z}^{n \times n}$ and distance matrix $D \in \mathbb{Z}^{n \times n}$, find a permutation $\pi$ (unit $u$ placed at location $\pi(u)$) minimizing

$$ C(\pi) = \sum_{u=1}^{n} \sum_{v=1}^{n} F_{uv} , D_{\pi(u)\pi(v)}. $$

Exact methods stall around $n \approx 30$; Taillard's robust tabu search (Ro-TS) remains the reference heuristic and the standard baseline on QAPLIB (see quadratic-assignment-problem for the problem-side view). The swap move exchanges the locations of units $r$ and $s$. With $p = \pi(r)$, $q = \pi(s)$, the cost change is computable in O(n):

$$ \Delta(\pi, r, s) = F_{rr}(D_{qq}{-}D_{pp}) + F_{rs}(D_{qp}{-}D_{pq}) + F_{sr}(D_{pq}{-}D_{qp}) + F_{ss}(D_{pp}{-}D_{qq}) + \sum_{k \neq r,s} \big[ F_{kr}(D_{\pi(k)q}{-}D_{\pi(k)p}) + F_{ks}(D_{\pi(k)p}{-}D_{\pi(k)q}) + F_{rk}(D_{q\pi(k)}{-}D_{p\pi(k)}) + F_{sk}(D_{p\pi(k)}{-}D_{q\pi(k)}) \big] $$

Ro-TS maintains the full $\Delta$ matrix. After executing swap $(r, s)$, every pair ${u, v}$ disjoint from ${r, s}$ updates in O(1) (Taillard 1991), with $\pi$ taken before the swap:

$$ \Delta'(u,v) = \Delta(u,v) + (F_{ru}{-}F_{rv}{+}F_{sv}{-}F_{su})(D_{\pi(r)\pi(u)}{-}D_{\pi(r)\pi(v)}{+}D_{\pi(s)\pi(v)}{-}D_{\pi(s)\pi(u)}) + (F_{ur}{-}F_{vr}{+}F_{vs}{-}F_{us})(D_{\pi(u)\pi(r)}{-}D_{\pi(v)\pi(r)}{+}D_{\pi(v)\pi(s)}{-}D_{\pi(u)\pi(s)}) $$

Only the 2n−3 pairs touching $r$ or $s$ need the O(n) formula again, so one iteration costs O(n²) total — scan, choice, and maintenance included — versus O(n³) for naive re-evaluation. Both formulas below handle general (asymmetric) matrices and are verified against brute-force recomputation.

import numpy as np

def qap_cost(F: np.ndarray, D: np.ndarray, perm: np.ndarray) -> int:
    """Full objective: sum_{u,v} F[u,v] * D[perm[u], perm[v]]."""
    return int(np.sum(F * D[np.ix_(perm, perm)]))

def swap_delta(F: np.ndarray, D: np.ndarray, perm: np.ndarray,
               r: int, s: int) -> int:
    """O(n) cost change of swapping perm[r] and perm[s] (general matrices)."""
    pr, ps = perm[r], perm[s]
    d = (F[r, r] * (D[ps, ps] - D[pr, pr])
         + F[r, s] * (D[ps, pr] - D[pr, ps])
         + F[s, r] * (D[pr, ps] - D[ps, pr])
         + F[s, s] * (D[pr, pr] - D[ps, ps]))
    idx = np.arange(len(perm))
    k = np.flatnonzero((idx != r) & (idx != s))
    pk = perm[k]
    d += int(np.sum(F[k, r] * (D[pk, ps] - D[pk, pr])
                    + F[k, s] * (D[pk, pr] - D[pk, ps])
                    + F[r, k] * (D[ps, pk] - D[pr, pk])
                    + F[s, k] * (D[pr, pk] - D[ps, pk])))
    return int(d)

def full_delta_matrix(F: np.ndarray, D: np.ndarray,
                      perm: np.ndarray) -> np.ndarray:
    """All swap deltas; O(n^3) but computed once at the start of the search."""
    n = len(perm)
    delta = np.zeros((n, n), dtype=np.int64)
    for r in range(n - 1):
        for s in range(r + 1, n):
            delta[r, s] = swap_delta(F, D, perm, r, s)
            delta[s, r] = delta[r, s]
    return delta

The search itself. Tabu attribute = (unit, location) pair: after moving unit $r$ away from location $\pi(r)$, the assignment "$r$ back at $\pi(r)$" is tabu for a tenure drawn uniformly from $[0.9n, 1.1n]$ — redrawn for every recording, which is what makes the method robust to instance-specific tuning. A swap is tabu only if both units would re-enter recently-left locations. Two aspiration tests: improved-best, and the forced long-term test (an assignment absent for $5n^2$ iterations forces the move) — the latter is Ro-TS's built-in frequency-style diversification.

import numpy as np
# Continues the same module: uses qap_cost, swap_delta, full_delta_matrix
# from the previous block.

def robust_tabu_search(F: np.ndarray, D: np.ndarray, n_iters: int = 10_000,
                       seed: int = 0,
                       tenure: tuple[int, int] | None = None,
                       aspiration: int | None = None) -> tuple[np.ndarray, int]:
    """Taillard (1991) robust tabu search for the QAP; O(n^2) per iteration.

    Returns (best permutation, best cost)."""
    n = F.shape[0]
    rng = np.random.default_rng(seed)
    t_lo, t_hi = tenure if tenure else (int(0.9 * n), int(np.ceil(1.1 * n)))
    asp = aspiration if aspiration else 5 * n * n
    perm = rng.permutation(n)
    cost = qap_cost(F, D, perm)
    best_perm, best_cost = perm.copy(), cost
    delta = full_delta_matrix(F, D, perm)
    tabu_until = np.zeros((n, n), dtype=np.int64)      # [unit, location]
    iu, ju = np.triu_indices(n, 1)                     # all swaps, vectorized
    for t in range(1, n_iters + 1):
        d = delta[iu, ju]
        ti = tabu_until[iu, perm[ju]]   # unit iu entering location perm[ju]
        tj = tabu_until[ju, perm[iu]]   # unit ju entering location perm[iu]
        tabu = (ti >= t) & (tj >= t)    # tabu only if BOTH re-entries are recent
        aspired = (cost + d < best_cost) | (ti < t - asp) | (tj < t - asp)
        admissible = ~tabu | aspired
        if admissible.any():
            cand = np.flatnonzero(admissible)
            k = int(cand[np.argmin(d[cand])])
        else:
            k = int(np.argmin(d))       # all moves tabu: take the best anyway
        r, s = int(iu[k]), int(ju[k])
        cost += int(delta[r, s])
        a = F[r, :] - F[s, :]           # O(1)-per-pair update, pre-swap perm
        b = D[perm[r], perm] - D[perm[s], perm]
        c = F[:, r] - F[:, s]
        e = D[perm, perm[r]] - D[perm, perm[s]]
        delta += (np.subtract.outer(a, a) * np.subtract.outer(b, b)
                  + np.subtract.outer(c, c) * np.subtract.outer(e, e))
        tabu_until[r, perm[r]] = t + rng.integers(t_lo, t_hi + 1)
        tabu_until[s, perm[s]] = t + rng.integers(t_lo, t_hi + 1)
        perm[r], perm[s] = perm[s], perm[r]
        for u in range(n):              # exact O(n) refresh of pairs touching r, s
            for v in (r, s):
                if u != v:
                    dv = swap_delta(F, D, perm, u, v)
                    delta[u, v] = dv
                    delta[v, u] = dv
        if cost < best_cost:
            best_cost, best_perm = cost, perm.copy()
    return best_perm, best_cost

# Tiny synthetic instance, checked against brute force.
from itertools import permutations
rng = np.random.default_rng(42)
n = 8
F = rng.integers(0, 10, (n, n)); np.fill_diagonal(F, 0)
D = rng.integers(1, 10, (n, n)); np.fill_diagonal(D, 0)
opt = min(qap_cost(F, D, np.array(p)) for p in permutations(range(n)))
best_perm, best_cost = robust_tabu_search(F, D, n_iters=2000, seed=1)
print(opt, best_cost, qap_cost(F, D, best_perm))
# Expected: 1155 1155 1155 — Ro-TS matches the brute-force optimum, and the
# independently recomputed cost confirms the incremental bookkeeping.

Implementation notes that matter in practice:

  • The incremental cost must be advanced with delta[r, s] before the outer-product update touches the matrix, and the update term must use the pre-swap permutation. Both orderings are easy to get wrong; the final qap_cost recomputation in the demo is the regression check that catches them.
  • The "both re-entries recent" tabu test (logical AND) is deliberately weaker than OR. With OR, after a few iterations almost every swap is tabu and the search runs on aspiration overrides alone.
  • For symmetric instances with zero diagonals the delta formula collapses to $2\sum_{k\neq r,s}(F_{sk}-F_{rk})(D_{\pi(r)\pi(k)}-D_{\pi(s)\pi(k)})$; keep the general form unless profiling says otherwise.

Worked Application 2: Job-Shop Scheduling with a Critical-Path Neighborhood

The job-shop problem (JSP): $n$ jobs, each a fixed sequence of operations over $m$ machines, one operation per machine per job; minimize makespan $C_{\max}$. A solution is an operation order for each machine. In the disjunctive-graph view, fixing all machine orders turns the schedule into a longest-path computation: each operation's start time is the longest path from the source, and $C_{\max}$ is the longest path overall. A critical path is a chain of operations realizing $C_{\max}$, linked by job-precedence or machine-precedence arcs.

Two classic results make the neighborhood (van Laarhoven, Aarts & Lenstra 1992, "Job shop scheduling by simulated annealing"):

  1. Swapping two adjacent critical operations on the same machine never creates a cycle — every such move is feasible without checking.
  2. Swapping adjacent operations that are not both critical can never reduce the makespan — so restricting moves to critical machine-arcs loses nothing. This is a candidate list strategy that is exact, not approximate.

The tabu attribute is the destroyed disjunctive arc: after swapping $u$ before $v$ into $v$ before $u$ on machine $M$, re-creating arc $u \to v$ on $M$ is tabu. This is the attribute family behind the best JSP tabu searches (Taillard 1994, "Parallel taboo search techniques for the job shop scheduling problem"; Nowicki & Smutnicki 1996, "A fast taboo search algorithm for the job shop problem"). Graph machinery first:

import numpy as np

def makespan(machines: np.ndarray, ptimes: np.ndarray,
             seq: list[list[int]]) -> tuple[int, np.ndarray, np.ndarray,
                                            np.ndarray, np.ndarray] | None:
    """Longest-path schedule for fixed machine sequences.

    machines[j, k] / ptimes[j, k]: machine and duration of job j's k-th op.
    seq[m]: job order on machine m. Operation id = j * n_mach + k.
    Returns (cmax, start, completion, job_pred, mach_pred); None if the
    sequences induce a cycle (only possible for externally supplied seq)."""
    n_jobs, n_mach = machines.shape
    n_ops = n_jobs * n_mach
    op_of = np.zeros((n_jobs, n_mach), dtype=np.int64)   # (job, machine) -> op
    for j in range(n_jobs):
        for k in range(n_mach):
            op_of[j, machines[j, k]] = j * n_mach + k
    p = ptimes.ravel()
    job_pred = np.full(n_ops, -1, dtype=np.int64)
    for j in range(n_jobs):
        job_pred[j * n_mach + 1: (j + 1) * n_mach] = \
            np.arange(j * n_mach, (j + 1) * n_mach - 1)
    mach_pred = np.full(n_ops, -1, dtype=np.int64)
    for m in range(n_mach):
        ops = op_of[seq[m], m]
        mach_pred[ops[1:]] = ops[:-1]
    succ_job = np.full(n_ops, -1, dtype=np.int64)
    succ_mach = np.full(n_ops, -1, dtype=np.int64)
    for o in range(n_ops):
        if job_pred[o] >= 0:
            succ_job[job_pred[o]] = o
        if mach_pred[o] >= 0:
            succ_mach[mach_pred[o]] = o
    indeg = (job_pred >= 0).astype(np.int64) + (mach_pred >= 0).astype(np.int64)
    start = np.zeros(n_ops, dtype=np.int64)
    comp = np.zeros(n_ops, dtype=np.int64)
    stack, seen = list(np.flatnonzero(indeg == 0)), 0
    while stack:                       # topological longest-path pass
        o = stack.pop()
        seen += 1
        st = 0
        if job_pred[o] >= 0:
            st = max(st, comp[job_pred[o]])
        if mach_pred[o] >= 0:
            st = max(st, comp[mach_pred[o]])
        start[o], comp[o] = st, st + p[o]
        for nx in (succ_job[o], succ_mach[o]):
            if nx >= 0:
                indeg[nx] -= 1
                if indeg[nx] == 0:
                    stack.append(nx)
    if seen < n_ops:
        return None
    return int(comp.max()), start, comp, job_pred, mach_pred

def critical_path(machines: np.ndarray, ptimes: np.ndarray,
                  seq: list[list[int]]) -> tuple[int, list[int], np.ndarray]:
    """One critical path (source to sink) and the machine-predecessor array."""
    res = makespan(machines, ptimes, seq)
    if res is None:
        raise ValueError("machine sequences induce a cycle")
    cmax, start, comp, job_pred, mach_pred = res
    o = int(np.argmax(comp))
    path = [o]
    while True:
        o = path[-1]
        if mach_pred[o] >= 0 and comp[mach_pred[o]] == start[o]:
            path.append(int(mach_pred[o]))
        elif job_pred[o] >= 0 and comp[job_pred[o]] == start[o]:
            path.append(int(job_pred[o]))
        else:
            break
    path.reverse()
    return cmax, path, mach_pred

def critical_machine_pairs(machines: np.ndarray, path: list[int],
                           mach_pred: np.ndarray) -> list[tuple[int, int, int]]:
    """Adjacent same-machine critical pairs as (machine, first_job, second_job)."""
    n_mach = machines.shape[1]
    pairs = []
    for a, b in zip(path[:-1], path[1:]):
        if mach_pred[b] == a:
            ja, jb = a // n_mach, b // n_mach
            pairs.append((int(machines[ja, a % n_mach]), int(ja), int(jb)))
    return pairs

The tabu search. Initial solution: order jobs on every machine by the operation's position inside its job (ties by job id) — this "level order" is provably acyclic and needs no scheduler. Tenure is randomized on a small range; the JSP arc attribute is much stronger than the QAP assignment attribute, so tenures of 8–14 work where QAP needs ~$n$. Each candidate move is evaluated by a full O(|ops|) recomputation, which is honest and simple; Taillard (1994) shows how head/tail values give an O(1) makespan lower bound per move when the scan itself becomes the bottleneck — see job-shop-scheduling for that refinement and for the Nowicki–Smutnicki block neighborhood.

import numpy as np
# Continues the same module: uses makespan, critical_path,
# critical_machine_pairs from the previous block.

def initial_sequences(machines: np.ndarray) -> list[list[int]]:
    """Acyclic start: order jobs on each machine by op position in the job."""
    n_jobs, n_mach = machines.shape
    pos = np.zeros((n_jobs, n_mach), dtype=np.int64)
    for j in range(n_jobs):
        for k in range(n_mach):
            pos[j, machines[j, k]] = k
    return [list(np.lexsort((np.arange(n_jobs), pos[:, m])))
            for m in range(n_mach)]

def apply_swap(seq: list[list[int]], m: int, ja: int, jb: int) -> list[list[int]]:
    """Copy of seq with jobs ja, jb (adjacent on machine m, ja first) swapped."""
    new = [list(s) for s in seq]
    i = new[m].index(ja)
    new[m][i], new[m][i + 1] = jb, ja
    return new

def jsp_tabu_search(machines: np.ndarray, ptimes: np.ndarray,
                    n_iters: int = 500, tenure: tuple[int, int] = (8, 14),
                    seed: int = 0) -> tuple[list[list[int]], int]:
    """Critical-path tabu search for the JSP (N1 neighborhood).

    Tabu attribute: the destroyed machine arc (machine, first_job, second_job).
    Returns (best machine sequences, best makespan)."""
    rng = np.random.default_rng(seed)
    seq = initial_sequences(machines)
    cmax, path, mach_pred = critical_path(machines, ptimes, seq)
    best_seq, best_cmax = [list(s) for s in seq], cmax
    tabu: dict[tuple[int, int, int], int] = {}
    for t in range(1, n_iters + 1):
        pairs = critical_machine_pairs(machines, path, mach_pred)
        if not pairs:
            break                                  # cmax = max job length: optimal
        chosen, chosen_val = None, None            # (move, seq, schedule result)
        fallback, fallback_exp = None, None
        for m, ja, jb in pairs:
            cand_seq = apply_swap(seq, m, ja, jb)  # feasible by construction
            res = critical_path(machines, ptimes, cand_seq)
            val = res[0]
            if tabu.get((m, jb, ja), 0) >= t and val >= best_cmax:
                exp = tabu[(m, jb, ja)]            # tabu, no aspiration
                if fallback_exp is None or exp < fallback_exp:
                    fallback, fallback_exp = ((m, ja, jb), cand_seq, res), exp
                continue
            if chosen_val is None or val < chosen_val:
                chosen, chosen_val = ((m, ja, jb), cand_seq, res), val
        if chosen is None:                         # all tabu: least-tabu fallback
            chosen = fallback
        (m, ja, jb), seq, (cmax, path, mach_pred) = chosen
        tabu[(m, ja, jb)] = t + int(rng.integers(tenure[0], tenure[1] + 1))
        if cmax < best_cmax:
            best_cmax, best_seq = cmax, [list(s) for s in seq]
    return best_seq, best_cmax

# ft06 (Fisher & Thompson 1963), 6 jobs x 6 machines, known optimum 55.
machines = np.array([[2, 0, 1, 3, 5, 4], [1, 2, 4, 5, 0, 3],
                     [2, 3, 5, 0, 1, 4], [1, 0, 2, 3, 4, 5],
                     [2, 1, 4, 5, 0, 3], [1, 3, 5, 0, 4, 2]])
ptimes = np.array([[1, 3, 6, 7, 3, 6], [8, 5, 10, 10, 10, 4],
                   [5, 4, 8, 9, 1, 7], [5, 5, 5, 3, 8, 9],
                   [9, 3, 5, 4, 3, 1], [3, 3, 9, 10, 4, 1]])
start_cmax = critical_path(machines, ptimes, initial_sequences(machines))[0]
best_seq, best_cmax = jsp_tabu_search(machines, ptimes, n_iters=500, seed=0)
print(start_cmax, best_cmax)
# Expected: 60 55 — from the level-order start (makespan 60) the critical-path
# tabu search reaches 55, the proven optimum of ft06, within 500 iterations.

Two observations worth internalizing. First, the neighborhood here is the candidate list: instead of all O(n·m) adjacent swaps, only critical-arc swaps are scanned (typically 3–10 on ft06-sized instances), and result 2 above guarantees this filter discards no improving move. Second, the uphill steps are essential: from many schedules every critical swap increases the makespan, and a memoryless descent stops exactly there; the tabu list is what lets the search traverse those ridges without immediately undoing them.

Advanced Techniques

Reactive tenure

Battiti & Tecchiolli (1994) replace tenure tuning with feedback: hash every visited solution; if a solution repeats within a window, the search is cycling — increase tenure multiplicatively. If no repetition occurs for a long stretch, decay tenure so the search is no more restricted than necessary. Canonical hashing of solutions is covered in fitness-evaluation-and-caching.

import numpy as np

class ReactiveTenure:
    """Battiti & Tecchiolli (1994): adapt tenure from solution repetitions."""

    def __init__(self, start: int, cap: int, window: int = 50,
                 up: float = 1.1, down: float = 0.9) -> None:
        self.tenure, self.cap, self.window = start, cap, window
        self.up, self.down = up, down
        self.seen: dict[int, int] = {}     # solution hash -> last iteration
        self.last_change = 0

    def update(self, solution_hash: int, t: int) -> int:
        """Call once per iteration with the current solution's hash."""
        last = self.seen.get(solution_hash)
        self.seen[solution_hash] = t
        if last is not None and t - last < self.window:        # cycling
            self.tenure = min(self.cap, max(self.tenure + 1,
                                            int(self.tenure * self.up)))
            self.last_change = t
        elif t - self.last_change > self.window:               # quiet: relax
            self.tenure = max(1, int(self.tenure * self.down))
            self.last_change = t
        return self.tenure

rt = ReactiveTenure(start=8, cap=40)
perm = np.array([3, 1, 0, 2])
h = hash(tuple(perm))
print(rt.update(h, t=1), rt.update(h, t=5))
# Expected: 8 9 — the repeat at t=5 (within the window) raises the tenure.

Frequency-based long-term memory and diversification

Recency memory forgets after $\tau$ iterations; frequency memory never forgets. Track transition counts (how often each move attribute fired) and residence counts (how long each attribute stayed in the solution). When the incumbent stalls, enter a diversification phase: add a penalty $w \cdot \text{freq}(m)/t$ to each move's delta so the scan prefers rarely used moves, or force the least-used assignment outright (the Ro-TS forced aspiration is exactly this, embedded in the aspiration test). Scale $w$ to the median |delta| of recently accepted moves so the penalty competes with, but does not drown, the objective.

import numpy as np

class FrequencyMemory:
    """Long-term transition/residence counts for a permutation problem."""

    def __init__(self, n: int) -> None:
        self.transition = np.zeros((n, n), dtype=np.int64)   # swap (i, j) count
        self.residence = np.zeros((n, n), dtype=np.int64)    # iters of i at pos
        self.iters = 0

    def record(self, perm: np.ndarray, r: int, s: int) -> None:
        """Call once per iteration, after executing swap (r, s)."""
        self.transition[r, s] += 1
        self.transition[s, r] += 1
        self.residence[np.arange(len(perm)), perm] += 1
        self.iters += 1

    def penalties(self, weight: float) -> np.ndarray:
        """Penalty matrix to ADD to the delta matrix while diversifying."""
        if self.iters == 0:
            return np.zeros_like(self.transition, dtype=np.float64)
        return weight * self.transition / self.iters

mem = FrequencyMemory(5)
perm = np.array([2, 0, 1, 4, 3])
mem.record(perm, 0, 1)
mem.record(perm, 0, 1)
mem.record(perm, 2, 3)
print(round(mem.penalties(10.0)[0, 1], 2), round(mem.penalties(10.0)[2, 3], 2))
# Expected: 6.67 3.33 — the over-used swap (0, 1) is penalized twice as hard.

Intensification by restarting from elite solutions

Quality memory: keep a bounded pool of the best distinct solutions seen. When diversification has not produced a new incumbent for a long phase, restart the search from a random elite with a cleared (or shortened) tabu list and a fresh tenure draw — the second visit explores different uphill exits because the randomized tenure changes which reversals are forbidden. Nowicki & Smutnicki's TSAB owes much of its strength to exactly this back-jump-to-elite mechanism. A stronger variant fixes the attributes shared by all elites (common assignments, common arcs) and searches only the residual problem; path relinking between elites is the systematic version of that idea. Restart policy comparisons and diversity measurement belong to diversity-and-population-management.

import numpy as np

class ElitePool:
    """Bounded pool of best distinct solutions for intensification restarts."""

    def __init__(self, capacity: int) -> None:
        self.capacity = capacity
        self.pool: list[tuple[float, tuple[int, ...]]] = []

    def offer(self, value: float, solution: np.ndarray) -> None:
        """Insert unless duplicate; keep the `capacity` best by value."""
        key = tuple(int(v) for v in solution)
        if any(k == key for _, k in self.pool):
            return
        self.pool.append((value, key))
        self.pool.sort(key=lambda e: e[0])
        del self.pool[self.capacity:]

    def restart_from(self, rng: np.random.Generator) -> np.ndarray:
        """Pick an elite uniformly at random as the new current solution."""
        _, sol = self.pool[rng.integers(len(self.pool))]
        return np.array(sol)

pool = ElitePool(capacity=3)
for v, p in [(10, [0, 1, 2]), (8, [1, 0, 2]), (12, [2, 1, 0]),
             (8, [1, 0, 2]), (7, [2, 0, 1])]:
    pool.offer(float(v), np.array(p))
print([v for v, _ in pool.pool])
# Expected: [7.0, 8.0, 10.0] — duplicate rejected, worst (12) evicted.

Candidate list strategies

When the full neighborhood scan is too slow, scan a subset — but choose the subset with structure, not blindly (Glover & Laguna 1997, ch. 3):

  • Exact filters: problem theory proves some moves cannot improve — the JSP critical-arc restriction above is the model case. Always look for one first.
  • Elite candidate list: scan the full neighborhood occasionally, keep the best $k$ moves; in between, scan only that list, refreshing it when its best move degrades below a threshold.
  • Successive filtering: cheap test first (e.g., lower bound on delta), exact evaluation only for survivors.
  • Random sampling: evaluate a random fraction per iteration; simplest, weakest, but unbiased — combine with first-improvement-style early exit.
  • Neighborhood decomposition: partition moves (by region, machine, vehicle) and rotate through partitions; pairs naturally with frequency memory to decide which partition deserves attention.

Strategic oscillation

For tightly constrained problems, let the search cross the feasibility boundary on purpose: allow infeasible moves with a violation penalty, increase the penalty while infeasible and decrease it while feasible (multiplicative factors ~1.2/0.85 every few iterations work as a start), and record boundary crossings in the tabu list so the oscillation does not thrash on one constraint. This turns the boundary region — where optima of capacity-style problems live — into the most-searched region instead of a wall.

Practical Challenges

The search still cycles despite the tabu list. The attribute is too fine (inverse-move attributes let other move orders rebuild the same solution) or the tenure too short. Diagnose by hashing visited solutions and logging repeats (fitness-evaluation-and-caching); fix by widening the attribute (element-position or feature level), randomizing tenure over a range, or switching to reactive tenure. Cycling with period $p$ needs tenure comfortably above $p$.

Every move is tabu and the search freezes. Tenure is too long relative to the number of distinct attributes (a tenure of $n$ with only $2n$ attributes freezes half the neighborhood), or the attribute is too coarse. Always implement the least-tabu fallback so an iteration executes something; if the fallback fires more than a few percent of iterations, shorten tenure or refine the attribute.

Incremental deltas drift away from the true objective. Order-of-update bugs (cost updated after the delta matrix, update term computed with the post-swap permutation) corrupt the search silently — it still runs, just on wrong numbers. In tests, recompute the objective from scratch every $k$ iterations and compare exactly; the worked examples above end with exactly this check. Keep the independent recomputation function separate from the incremental code path.

The neighborhood scan dominates runtime. Profile before optimizing: if delta lookup is O(1) but the scan is a Python loop, vectorize the scan over index arrays as in the Ro-TS implementation (mask + argmin over np.triu_indices). If evaluation itself is the cost, install candidate lists, preferring exact filters. Only then consider numba on the inner loop.

Aspiration never fires, or fires constantly. Never fires: the improved-best test is too strict for long runs on large instances — fine, that is its job; check the least-tabu fallback instead. Fires constantly: tenure or attribute is over-restrictive and aspiration has become the de-facto selection rule, which silently degrades the method to plain best-improvement; fix the memory, do not loosen the aspiration further.

Strong basins recapture the search after every escape. The tabu list handles small basins; deep, wide basins need long-term memory. Trigger a diversification phase after $0.5n$–$5n$ improvement-free iterations: frequency-penalized selection, a forced sequence of least-used moves, or a restart from an elite plus a perturbation. If diversification phases dominate the run, the neighborhood itself is too myopic — revisit local-search-and-neighborhoods before adding more memory.

Results are not reproducible or vary wildly across seeds. Pass explicit np.random.default_rng(seed) generators everywhere (tenure draws, tie-breaking, restarts); never use global random state. Report best/mean/worst over ≥10 seeds with fixed budgets. If the across-seed spread is large, the search is start-sensitive: strengthen diversification or lengthen runs rather than cherry-picking the best seed.

Tenure tuning consumes the project. Stop tuning constants. Adopt randomized tenure on a range tied to instance size ($[0.9n, 1.1n]$ for QAP-like attributes, small constants for strong arc attributes), or reactive tenure. Tune the range, not a point value, and only on a training split of instances.

Tools & Libraries

Library / resource When to use Note
numpy always — delta matrices, vectorized scans, RNG default_rng(seed); vectorize the scan before reaching for compilation
numba inner move loops that resist vectorization (JSP graph passes) @njit the longest-path evaluation; keep a pure-Python reference for tests
networkx independent validation of longest-path / DAG claims in scheduling use dag_longest_path_length as the slow oracle in unit tests
OR-Tools CP-SAT exact baselines for JSP/timetabling on small-medium instances confirms optima your tabu search should hit on test instances
gurobipy exact MIP baselines; linearized QAP for tiny instances exact methods stall ~n=30 on QAP — small instances only
QAPLIB / OR-Library / Taillard sets standard benchmark instances with best-known values report gaps to best-known, not just raw objectives
Optuna tuning tenure ranges, diversification triggers, penalty weights tune on a training instance split; hold out test instances
pandas run logs and result tables (one row per run) aggregate per instance × seed; feed statistical comparison

Output Format

A complete tabu search deliverable contains:

  1. Design summary table — every memory decision made explicit:
Component Choice
Neighborhood all pairwise swaps, |N| = n(n−1)/2, O(1) delta lookup
Tabu attribute (unit, location) assignments destroyed by the move
Tabu test both re-entries recent (AND)
Tenure uniform on [0.9n, 1.1n], redrawn per recording
Aspiration improved-best; forced after 5n² absent iterations; least-tabu fallback
Long-term memory transition frequencies; diversification after 2n idle iterations
Intensification elite pool of 10; restart on diversification failure
Budget / stopping 10⁵ iterations or 60 s, whichever first

  1. Run report — one row per (instance, seed): instance name, $n$, seed, best objective, best-known value, gap %, iteration of last improvement, total iterations, wall-clock seconds. Aggregate per instance: best/mean/worst over seeds.

  2. Convergence evidence — best-so-far trace per run (the trace list in the generic engine); plot mean ± range over seeds. Flag the fraction of budget after the last improvement: above ~50 % suggests the budget is wasted or diversification is broken.

  3. Validation statement — best solutions re-checked by an independent objective recomputation (never the incremental value), feasibility verified by a separate checker, and agreement with exact baselines on all small test instances (as in both worked examples above).

  4. Artifacts — solutions as plain arrays/JSON with instance and seed metadata; traces as CSV (iteration, current, best); the exact config (all parameters + git hash) embedded in the result file.

Questions to Ask

  • What problem and representation — permutation, assignment, sequences per machine, binary vector? Is a neighborhood with cheap delta evaluation already available?
  • What instance sizes now, and at the target scale? What per-run time budget?
  • Is there a known-optimum or best-known reference set for validation and gap reporting?
  • Must the search stay feasible at all times, or is controlled infeasibility (penalties, oscillation) acceptable?
  • Is there an existing local search to upgrade — or are we choosing among metaheuristics, in which case why tabu search over SA/ILS for this problem?
  • What should be reported: single best solution, or a statistical comparison across seeds and instances?
  • Any symmetry in solutions (relabeling, rotations) that makes attribute design and cycle detection harder?
  • Will the implementation later need parallel runs or hybridization (e.g., tabu search inside a memetic loop) that constrain the code structure?

Related Skills

  • local-search-and-neighborhoods — when the neighborhood, move data structures, or delta evaluation needs design work before any memory is added; tabu search inherits all of it.
  • fitness-evaluation-and-caching — when move evaluation is the bottleneck or cycling must be detected: O(n) QAP swap updates, canonical solution hashing, memoized evaluators.
  • quadratic-assignment-problem — when the user works on QAP itself: formulations, linearizations, QAPLIB benchmarks, and why robust tabu search is the method of choice.
  • job-shop-scheduling — when the user needs the full JSP toolbox: disjunctive MIP, CP-SAT interval models, Nowicki–Smutnicki block neighborhoods, shifting bottleneck.
  • diversity-and-population-management — when stagnation persists and restart policies, diversity measurement, or population-level memory need more machinery than frequency penalties.

Related in development

update-docs-vercel-next-js

This skill should be used when the user asks to "update documentation for my changes", "check docs for this PR", "what d...

by vercel
137,316

docstring-pytorch-pytorch

Write docstrings for PyTorch functions and methods following PyTorch conventions. Use when writing or updating docstring...

by pytorch
96,875

at-dispatch-v2

Convert PyTorch AT_DISPATCH macros to AT_DISPATCH_V2 format in ATen C++ code. Use when porting AT_DISPATCH_ALL_TYPES_AND...

dispatch
by pytorch
95,362