Local Search And Neighborhoods

Local Search and Neighborhoods

You are an expert in neighborhood-based search for combinatorial optimization. This skill covers neighborhood design (swap, insertion, 2-opt, Or-opt, exchange), constant- and linear-time delta evaluation, first- vs best-improvement pivoting, neighborhood scanning order, move data structures, and the limits of plain hill climbing. Use the framework below to build descent engines that are correct (deltas audited against full recomputation) and fast (no full objective evaluation inside the move loop), and that drop in unchanged as the inner loop of simulated annealing, tabu search, and iterated local search.

Initial Assessment

Establish these facts before proposing a neighborhood or writing any move code:

• Representation. Permutation (tour, schedule), binary vector (selection), integer

assignment array, or partition into sets (routes, machine loads)? The move catalog and the delta formulas depend entirely on this. See solution-encodings-style criteria: every move must map feasible representations to feasible representations, or you must plan a penalty scheme.

• Objective structure. Additive over solution elements (sum of edge lengths, sum of

assignment costs, linear penalties)? Additive objectives almost always admit O(1) or O(n) deltas. Max-type / critical-path objectives (flow-shop or job-shop makespan) do not; plan for an acceleration scheme (Taillard heads/tails, critical-path filtering) instead of a true delta.

• Hard vs soft constraints. Decide per constraint: keep moves feasibility-preserving,

or relax the constraint into the objective and evaluate the penalty change inside the delta. Mixing the two without a plan is the most common source of wrong deltas.

• Instance size. n ≤ 500: any O(n²) neighborhood with full sweeps is fine in plain

Python. n in 10³–10⁴: you need candidate lists, don't-look bits, or vectorized move scoring. n above 10⁴: you also need O(1) position queries and careful apply costs (segment reversal is O(n) per move).

• Evaluation cost. Time one full objective evaluation and one delta evaluation early.

If the delta is not at least ~n times cheaper than the full evaluation, the neighborhood implementation is wrong or the objective does not decompose.

• Role of the local search. Standalone descent, multistart, or inner loop of a

metaheuristic? The wrapper changes the acceptance logic; the move/delta code should not change at all. Keep them in separate functions from day one.

• Quality target. 2-opt local optima on random Euclidean TSP sit roughly 5% above

optimal from good starts (Johnson & McGeoch 1997, "The traveling salesman problem: a case study in local optimization"). If the target is tighter, plan a metaheuristic wrapper or a richer neighborhood, not more sweeps.

• Data format. Full distance matrix (O(n²) memory) vs coordinates plus k-nearest

neighbors via a KD-tree? At n = 20,000, a float64 matrix is 3.2 GB — candidate lists are not optional there.

• Determinism. Seed every random scan order and tie-break (np.random.default_rng)

so descents are reproducible run to run.

• Validation plan. An independent objective recompute and a delta audit (compare each

claimed delta against a from-scratch evaluation on random moves) must exist before any performance tuning.

Neighborhood Anatomy

A neighborhood is a function $N: S \to 2^S$ mapping each solution to the set of solutions reachable by one move. Local search descends until it reaches a local optimum:

$$ s^\ \text{ is a local optimum of } f \text{ w.r.t. } N \iff f(s^\) \le f(s) \quad \forall s \in N(s^\*). $$

A move $m$ applied to $s$ is written $s \oplus m$, and the entire craft of fast local search is computing the delta

$$ \Delta(s, m) = f(s \oplus m) - f(s) $$

without ever computing $f(s \oplus m)$ from scratch. Local optimality is always relative to N: a 2-opt optimum is usually not an Or-opt optimum, which is the basis of compound neighborhoods and variable neighborhood descent.

Move catalog

| Move | Acts on | Neighborhood size | Delta cost | Notes | |---|---|---|---|---| | Swap (exchange two elements) | permutation, assignment array | n(n−1)/2 | O(1) for additive costs; O(n) for QAP | Weak on tours (touches 4 edges but rarely improving); strong on assignments | | Insertion (shift one element) | permutation | (n−1)² | O(1) additive; O(nm) per job in flow shop via Taillard | The dominant scheduling neighborhood | | 2-opt (reverse a segment) | tours | n(n−1)/2 | O(1) symmetric distances; no O(1) delta if asymmetric | Apply costs O(n) (reversal); Croes (1958) | | Or-opt (relocate segment of 1–3) | tours | ≈ 3n² | O(1) | Orientation-preserving, so valid for asymmetric instances; Or (1976) | | k-flip | binary vectors | C(n,k) | O(degree) per flip | Knapsack, max-cut, SAT-style problems | | Inter-set exchange / relocate | partitions, routes | O(n²) | O(1)–O(route length) | CVRP inter-route moves; feasibility (load, time) checked in the delta |

Delta evaluation

For the symmetric TSP, reversing positions $i{+}1 \dots j$ of tour $\pi$ removes two edges and adds two edges; every interior edge keeps its (symmetric) length:

$$ \Delta{\text{2opt}}(i, j) = d{\pii \pij} + d{\pi{i+1} \pi{j+1}} - d{\pii \pi{i+1}} - d{\pij \pi_{j+1}}. $$

This O(1) formula is why 2-opt scans millions of moves per second. On asymmetric instances the reversed segment changes direction, so the formula is invalid — use orientation-preserving moves (Or-opt, 3-opt segment reinsertion) instead. The same discipline applies everywhere: write the delta as (cost of edges/terms created) − (cost of edges/terms destroyed), and prove to yourself that nothing else changes.

First vs best improvement

• Best improvement scans the whole neighborhood, applies the steepest move. Fewer,

larger steps; deterministic given the scan; pairs naturally with vectorized scoring of all moves at once.

• First improvement applies the first improving move found. Much cheaper sweeps early

in the descent when improving moves are dense; final quality is usually on par with best improvement, sometimes better, especially from random starts (Hansen & Mladenović 2006, "First vs. best improvement: an empirical study").

• Practical default: first improvement with a randomized scan order per sweep (removes

lexicographic bias), switching to best improvement only when moves are scored in bulk with numpy.

• The termination certificate is one full improving-move-free sweep. That last sweep costs

|N| delta evaluations no matter what — pruning rules and candidate lists are the only way to shrink it.

Scanning order and move data structures

Scan order matters because first improvement commits to whatever it sees first. Lexicographic order biases the search toward low indices; a seeded random permutation per sweep is the cheap fix. Queue-driven orders (re-examine elements whose surroundings just changed) are the systematic version and lead directly to don't-look bits (Bentley 1992, "Fast algorithms for geometric traveling salesman problems").

Keep these structures current under every apply:

• pos array with pos[element] = position, the inverse of the permutation — O(1) move

legality and delta lookups, updated in O(segment) on reversal.

• Cached per-element quantities used by deltas (route load, completion times, gains), each

with an O(1)/O(k) incremental update rule.

• For very large tours, array reversal is the bottleneck; doubly linked or two-level

list representations bound the apply cost (Fredman, Johnson, McGeoch & Ostheimer 1995, "Data structures for traveling salesmen").

Hill climbing and its limits

Plain hill climbing terminates at the first local optimum of $N$ — typically in O(n) to O(n log n) accepted moves empirically, and there it stops. Larger neighborhoods give better local optima at higher sweep cost; no fixed polynomial neighborhood removes local optima for NP-hard problems. Escaping is the wrapper's job: probabilistic acceptance (simulated annealing), memory (tabu search), perturbation and restarts (iterated local search). Build the descent so those wrappers reuse the identical move and delta code.

Generic Descent Engine

The reusable skeleton, independent of problem and move type:

LOCAL-SEARCH(s, N, rule):
    f_s  tuple[float, dict[str, int]]:
    """Generic descent. enumerate_moves() yields (move, delta) pairs for the
    CURRENT solution; apply_move mutates the solution in place."""
    stats = {"sweeps": 0, "moves": 0, "delta_evals": 0}
    while stats["moves"]  None:
    """Hill climbing for weighted max-cut: flip one vertex side, O(n) delta."""
    rng = np.random.default_rng(seed)
    w = rng.uniform(0.0, 1.0, (n, n))
    w = np.triu(w, 1) + np.triu(w, 1).T            # symmetric, zero diagonal
    side = rng.integers(0, 2, n)

    def cut_value(s: np.ndarray) -> float:
        return float(w[s[:, None] != s[None, :]].sum() / 2.0)

    def enumerate_moves() -> Iterable[tuple[int, float]]:
        same = side[:, None] == side[None, :]
        gain = (w * same).sum(axis=1) - (w * ~same).sum(axis=1)
        for i in np.argsort(-gain):                # promising flips first
            yield int(i), -float(gain[i])          # we minimize -cut

    def apply_move(i: int) -> None:
        side[i] ^= 1

    obj, stats = local_search(-cut_value(side), enumerate_moves, apply_move,
                              improvement="first")
    assert abs(obj + cut_value(side)) 20 rarely pays off |
| Or-opt segment lengths | 1–3 | Longer segments: bigger neighborhood, diminishing returns |
| Max sweeps / move budget | bound by time, not count, in experiments | A budget cap protects wrappers (ILS, SA) from runaway descents |
| Restarts (standalone use) | 5–50 seeded starts, keep best | Cheap variance reduction; superseded by ILS-style perturbation |
| Delta audit frequency | every run in tests; every ~10⁴ moves in production | Catches drift and wrong delta formulas at negligible cost |

## Worked Example: TSP 2-opt with O(1) Delta Evaluation

Tour as a numpy position array `tour` (city at each position). The move `(i, j)` reverses
`tour[i+1 : j+1]`; the delta is the four-edge formula above. First, the scalar
first-improvement descent with a randomized scan order:

```python
import numpy as np

def tour_length(tour: np.ndarray, dist: np.ndarray) -> float:
    """Length of the closed tour (array of city indices by position)."""
    return float(dist[tour, np.roll(tour, -1)].sum())

def two_opt_delta(tour: np.ndarray, dist: np.ndarray, i: int, j: int) -> float:
    """O(1) length change of reversing tour[i+1 : j+1], for positions i  tuple[np.ndarray, float]:
    """First-improvement 2-opt descent with a randomized scan order per sweep."""
    rng = np.random.default_rng(seed)
    tour = tour.copy()
    n = len(tour)
    length = tour_length(tour, dist)
    improved = True
    while improved:
        improved = False
        for i in rng.permutation(n - 2):
            i = int(i)
            for j in range(i + 2, n - (i == 0)):   # skip the wrap pair (0, n-1)
                delta = two_opt_delta(tour, dist, i, j)
                if delta  tuple[np.ndarray, np.ndarray]:
    """n points on the unit circle; the optimal tour visits them in angular order."""
    rng = np.random.default_rng(seed)
    theta = np.sort(rng.uniform(0.0, 2.0 * np.pi, n))
    pts = np.column_stack([np.cos(theta), np.sin(theta)])
    dist = np.linalg.norm(pts[:, None] - pts[None, :], axis=2)
    return pts, dist

_, dist = circle_instance(12, seed=42)
start = np.random.default_rng(7).permutation(12)
tour, length = two_opt_first(start, dist)
print(f"start={tour_length(start, dist):.4f} final={length:.4f} "
      f"optimal={tour_length(np.arange(12), dist):.4f}")
# Expected: start=14.6762 final=5.9567 optimal=5.9567. For points in convex
# position every 2-opt local optimum is crossing-free, and the only
# crossing-free tour through points in convex position is the hull tour.

Best improvement becomes attractive once the whole delta matrix is computed in one shot. With succ = np.roll(tour, -1) and edge[i] the cost of the edge leaving position i, the n×n delta matrix is four broadcasts — one (n, n) temporary per sweep (n = 5,000 → ~200 MB, so this variant is for n up to a few thousand):

import numpy as np

def two_opt_best_vectorized(tour: np.ndarray, dist: np.ndarray,
                            tol: float = 1e-9) -> tuple[np.ndarray, float]:
    """Best-improvement 2-opt; each sweep scores all O(n^2) moves at once."""
    tour = tour.copy()
    n = len(tour)
    length = float(dist[tour, np.roll(tour, -1)].sum())
    iu = np.triu_indices(n, k=2)                   # positions with j >= i+2
    while True:
        succ = np.roll(tour, -1)
        edge = dist[tour, succ]                    # cost of edge at each position
        delta = (dist[tour[:, None], tour[None, :]]
                 + dist[succ[:, None], succ[None, :]]
                 - edge[:, None] - edge[None, :])
        k = int(np.argmin(delta[iu]))
        i, j = int(iu[0][k]), int(iu[1][k])
        if delta[i, j] >= -tol:
            return tour, length
        tour[i + 1 : j + 1] = tour[i + 1 : j + 1][::-1]
        length += float(delta[i, j])

rng = np.random.default_rng(3)
pts = rng.random((100, 2))
dist = np.linalg.norm(pts[:, None] - pts[None, :], axis=2)
start = rng.permutation(100)
tour, length = two_opt_best_vectorized(start, dist)
print(f"start={float(dist[start, np.roll(start, -1)].sum()):.3f} final={length:.3f}")
# Expected: start=46.152 final=8.733 -- about 81% shorter; a best-improvement
# 2-opt local optimum for 100 uniform points (optimal is typically ~7.7-8.0,
# so this sits 5-12% above it, normal for 2-opt from a random start).

The wrap pair (i=0, j=n−1) stays in the index set because its delta is exactly zero (reversing everything after position 0 re-traverses the same closed tour), so it can never be selected as improving.

Every delta formula gets an audit before it gets trusted. This is non-negotiable: a wrong delta produces a descent that silently reports objectives that are not the objective.

import numpy as np

def audit_two_opt_deltas(n: int = 60, trials: int = 500, seed: int = 1) -> float:
    """Largest |claimed delta - true length change| over random 2-opt moves."""
    rng = np.random.default_rng(seed)
    pts = rng.random((n, 2))
    dist = np.linalg.norm(pts[:, None] - pts[None, :], axis=2)
    tour = rng.permutation(n)
    base = float(dist[tour, np.roll(tour, -1)].sum())
    worst = 0.0
    for _ in range(trials):
        i = int(rng.integers(0, n - 3))
        j = int(rng.integers(i + 2, n - (i == 0)))
        a, b, c, d = tour[i], tour[i + 1], tour[j], tour[(j + 1) % n]
        delta = float(dist[a, c] + dist[b, d] - dist[a, b] - dist[c, d])
        trial = tour.copy()
        trial[i + 1 : j + 1] = trial[i + 1 : j + 1][::-1]
        true_change = float(dist[trial, np.roll(trial, -1)].sum()) - base
        worst = max(worst, abs(true_change - delta))
    return worst

print(f"max |delta - true change| = {audit_two_opt_deltas():.2e}")
# Expected: below 1e-12 -- pure floating-point rounding; the formula is exact.

2-opt local optima still contain misplaced chains of cities that no segment reversal can fix. Or-opt relocates segments of length 1–3 (optionally reversed) elsewhere in the tour — an O(1) delta with three edges destroyed and three created, and since the segment keeps its orientation, the formula is also valid for asymmetric instances:

import numpy as np

def closed_tour_length(tour: list[int], dist: np.ndarray) -> float:
    """Length of a closed tour stored as a Python list of city indices."""
    n = len(tour)
    return float(sum(dist[tour[k], tour[(k + 1) % n]] for k in range(n)))

def first_or_opt_move(tour: list[int], 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.