Grasp

GRASP — Greedy Randomized Adaptive Search Procedures

You are an expert in metaheuristic optimization, specializing in GRASP and multi-start methods. This skill covers greedy randomized construction with restricted candidate lists (RCL), the construction/local-search iteration, alpha calibration, reactive GRASP, bias functions, and hybridization with path relinking. Use the framework below to decide whether GRASP fits a problem, implement it correctly in vectorized numpy, and report results that withstand reviewer scrutiny.

Initial Assessment

Establish the following before proposing or writing any GRASP code:

• Greedy structure. Identify the element-by-element construction: what is an "element" (a column, an assignment, an edge), and what incremental greedy score does adding it have? GRASP requires a cheap, meaningful greedy function. If none exists, reconsider the method.
• Problem class and objective. Minimization or maximization, single objective, and which constraints the construction must respect at every partial step.
• Hard vs soft constraints. Hard constraints shape the candidate set during construction; soft constraints belong in the objective or in penalties.
• Instance size. Number of elements per solution, candidates per construction step, and whether greedy scores can be updated incrementally instead of recomputed.
• Local search. Which neighborhood improves constructed solutions, and whether O(1)/O(n) delta evaluation exists. GRASP without a competent improvement phase is just semi-greedy sampling and usually underperforms.
• Evaluation cost. Time per construction and per local-search descent determines how many restarts fit the budget.
• Exact-vs-heuristic need. If instances are small enough for a MIP solver to close the gap in the available time, solve exactly and use GRASP only for warm starts.
• Time budget and stopping rule. Wall-clock limit, iteration limit, or stagnation-based stop; restarts are independent, so the budget maps linearly to iterations.
• Quality requirement. Target gap to best-known solutions or bounds; whether a single good solution or a distribution over seeds is required.
• Data format. How instances arrive (matrix files, benchmark formats, generated) and what an independent feasibility check looks like.
• Memory features in scope. Plain GRASP is memoryless; decide up front whether reactive alpha, hashing of duplicates, or path relinking is part of the deliverable.
• Parallel resources. Independent restarts parallelize with near-linear speedup; count available cores.
• Reproducibility. One np.random.default_rng(seed) per run; record seeds in all result tables.

Algorithm Anatomy

GRASP (Feo & Resende 1989, "A probabilistic heuristic for a computationally difficult set covering problem"; Feo & Resende 1995, "Greedy randomized adaptive search procedures") is a multi-start metaheuristic. Each iteration has two phases:

1. Greedy randomized construction. Build a complete solution element by element. At each step, score every feasible candidate with the greedy function, form a restricted candidate list (RCL), and pick one RCL member uniformly at random.
2. Local search. Descend from the constructed solution to a local optimum.

The incumbent over all iterations is returned. Iterations are independent: basic GRASP carries no state between restarts except the incumbent.

The RCL. Let $C$ be the current candidate set and $c(e)$ the incremental greedy cost of element $e$ (lower is better), with $c_{\min} = \min_{e \in C} c(e)$ and $c_{\max} = \max_{e \in C} c(e)$. The value-based RCL with parameter $\alpha \in [0, 1]$ is

$$ \mathrm{RCL}(\alpha) ;=; {, e \in C ;:; c(e) ,\le, c_{\min} + \alpha,(c_{\max} - c_{\min}) ,}. $$

The cardinality-based alternative keeps the $\lceil \alpha,|C| \rceil$ best-scored candidates. Semantics of the extremes:

• $\alpha = 0$: pure greedy construction (random only in tie-breaking).
• $\alpha = 1$: uniformly random construction.

Adaptive is the operative word: greedy scores are recomputed (or incrementally updated) after every selection, because the value of a candidate depends on the partial solution built so far. A non-adaptive variant that scores candidates once is cheaper but markedly weaker on problems with strong interaction between elements.

Why GRASP works. Pure greedy is a single, biased sample: it commits early to myopic choices and lands in one basin of attraction. Pure random restarts cover many basins but start local search from poor solutions, wasting descent time. Semi-greedy sampling (Hart & Shogan 1987, "Semi-greedy heuristics: an empirical study") produces starting points that are simultaneously good and diverse, so local search converges fast and to varied local optima. The minimum over iterations then improves with the number of restarts. Empirically, GRASP time-to-target values follow an approximately exponential distribution (Aiex, Resende & Ribeiro 2002, "Probability distribution of solution time in GRASP"), which has a practical consequence: running $p$ independent GRASPs in parallel divides expected time-to-target by roughly $p$.

Decision guidance.

• Use GRASP when a natural constructive heuristic exists (covering, location, assignment, scheduling, routing) and its score is cheap to evaluate or update.
• Prefer GRASP over iterated local search when high-quality local optima are scattered across the search space; prefer ILS when they cluster in a "big valley," because perturbing a good solution then beats rebuilding from scratch.
• Prefer GRASP over population methods when memory and diversity machinery would be overkill, or when an embarrassingly parallel method is wanted.
• Plain GRASP learns nothing across iterations. If the iteration budget is large, add memory: reactive alpha selection and path relinking are the two standard, well-tested upgrades (Resende & Ribeiro 2016, "Optimization by GRASP").

Complexity. One iteration costs the construction plus the descent. A construction of $k$ elements with $|C|$ candidates per step and incremental score updates costs $O(k,|C|)$ score evaluations; recomputing scores from scratch each step multiplies that by the update cost. Local search usually dominates; budget it accordingly and see local-search-and-neighborhoods for delta-evaluation patterns that keep it cheap.

Generic Framework

The skeleton every GRASP shares:

GRASP(alpha, max_iters):
    best <- None
    for k = 1 .. max_iters:                       # independent restarts
        x <- GreedyRandomizedConstruction(alpha)
        x <- LocalSearch(x)                       # descend to a local optimum
        if cost(x) < cost(best): best <- x
    return best

GreedyRandomizedConstruction(alpha):
    x <- empty partial solution
    while x is not complete:
        C    <- feasible candidate elements given x
        c(e) <- greedy score of e for all e in C  # adaptive: depends on x
        cmin, cmax <- min / max of c over C
        RCL  <- { e in C : c(e) <= cmin + alpha * (cmax - cmin) }
        e    <- uniform random element of RCL
        x    <- x + e                             # update scores incrementally
    return x

The RCL selection is the only piece that repeats across problems, so isolate it:

import numpy as np


def rcl_choice(scores: np.ndarray, alpha: float, rng: np.random.Generator,
               kind: str = "value") -> int:
    """Pick one candidate index from the restricted candidate list.

    scores: greedy score per candidate, lower = better, np.inf = infeasible.
    kind:   "value" keeps scores within alpha of the min-max range;
            "cardinality" keeps the best ceil(alpha * n_feasible) candidates.
    """
    feasible = np.flatnonzero(np.isfinite(scores))
    s = scores[feasible]
    if kind == "value":
        cutoff = s.min() + alpha * (s.max() - s.min())
        rcl = feasible[s <= cutoff + 1e-12]
    else:
        k = max(1, int(np.ceil(alpha * len(feasible))))
        rcl = feasible[np.argpartition(s, k - 1)[:k]]
    return int(rng.choice(rcl))

The driver is problem-independent. It takes the construction and the local search as callables, so the same loop serves every worked example below:

import numpy as np
from collections.abc import Callable


def grasp(construct: Callable[[float, np.random.Generator], np.ndarray],
          local_search: Callable[[np.ndarray], np.ndarray],
          cost: Callable[[np.ndarray], float],
          n_iters: int = 200, alpha: float = 0.2,
          seed: int = 0) -> tuple[np.ndarray, float, np.ndarray]:
    """Plain multi-start GRASP; returns (best, best_cost, best-so-far curve)."""
    rng = np.random.default_rng(seed)
    best, best_cost = None, np.inf
    curve = np.empty(n_iters)
    for k in range(n_iters):
        x = local_search(construct(alpha, rng))
        c = cost(x)
        if c < best_cost:
            best, best_cost = x.copy(), c
        curve[k] = best_cost
    return best, best_cost, curve


# Tiny wiring check: select exactly 2 of 5 items at minimum total cost.
_item_costs = np.array([5.0, 1.0, 4.0, 2.0, 9.0])

def _toy_construct(alpha: float, rng: np.random.Generator) -> np.ndarray:
    """Greedy randomized selection of 2 items (value-based RCL on item cost)."""
    chosen = np.zeros(5, dtype=bool)
    for _ in range(2):
        scores = np.where(chosen, np.inf, _item_costs)
        feasible = np.flatnonzero(np.isfinite(scores))
        s = scores[feasible]
        cutoff = s.min() + alpha * (s.max() - s.min())
        chosen[int(rng.choice(feasible[s <= cutoff + 1e-12]))] = True
    return chosen

best, best_cost, _ = grasp(_toy_construct, lambda x: x,
                           lambda x: float(_item_costs[x].sum()),
                           n_iters=50, alpha=0.3, seed=1)
print(best_cost)
# Expected: 3.0 -- the items costing 1 and 2, found within the 50 restarts.

Parameter guidance

Operator-level details (neighborhood catalogs, repair operators, delta tables) are covered by their own skills — see local-search-and-neighborhoods for moves and fitness-evaluation-and-caching for incremental evaluation; do not reinvent them inside the GRASP loop.

Worked Example 1: Set Covering GRASP

The set covering problem (SCP): given rows $1..m$, columns $j = 1..n$ with cost $c_j > 0$ and a boolean cover matrix $A \in {0,1}^{m \times n}$, choose a column set covering every row at minimum cost:

$$ \min \sum_{j=1}^{n} c_j x_j \quad \text{s.t.} \quad \sum_{j: a_{ij}=1} x_j \ge 1 ;; \forall i, \qquad x_j \in {0,1}. $$

The classic greedy score is cost per newly covered row, $c_j / |S_j \cap U|$ where $U$ is the uncovered set (Chvátal 1979, "A greedy heuristic for the set-covering problem"). This was GRASP's original application (Feo & Resende 1989). For models, bounds, and Lagrangian alternatives on this problem family, see set-covering-packing-partitioning.

Instance generator and an independent feasibility check:

import numpy as np


def generate_scp(m: int, n: int, density: float,
                 seed: int) -> tuple[np.ndarray, np.ndarray]:
    """Random SCP instance: boolean cover matrix A (m x n) and costs c (n,)."""
    rng = np.random.default_rng(seed)
    A = rng.random((m, n)) < density
    for i in np.flatnonzero(A.sum(axis=1) < 2):       # guarantee coverability
        A[i, rng.choice(n, size=2, replace=False)] = True
    c = 1.0 + rng.integers(0, 99, size=n).astype(float)
    return A, c


def scp_is_cover(x: np.ndarray, A: np.ndarray) -> bool:
    """Independent feasibility check: every row covered by a selected column."""
    return bool(A[:, x].any(axis=1).all())

Construction. Solutions are boolean column masks. The score vector is rebuilt after every selection from the shrinking uncovered set — the adaptive part:

import numpy as np


def construct_scp(A: np.ndarray, c: np.ndarray, alpha: float,
                  rng: np.random.Generator) -> np.ndarray:
    """Greedy randomized SCP construction; returns a boolean column mask.

    Greedy score: cost per newly covered row, recomputed after each pick.
    """
    m, n = A.shape
    x = np.zeros(n, dtype=bool)
    uncovered = np.ones(m, dtype=bool)
    while uncovered.any():
        gain = A[uncovered].sum(axis=0).astype(float)   # new rows per column
        score = np.where((gain > 0) & ~x, c / np.maximum(gain, 1.0), np.inf)
        feasible = np.flatnonzero(np.isfinite(score))
        s = score[feasible]
        cutoff = s.min() + alpha * (s.max() - s.min())
        j = int(rng.choice(feasible[s <= cutoff + 1e-12]))
        x[j] = True
        uncovered &= ~A[:, j]
    return x

Local search: first drop redundant columns (every covered row stays covered), then a first-improvement 1-out exchange — remove one column, repair greedily, keep the trial if it is strictly cheaper. The repair is a deterministic version of the construction:

import numpy as np


def drop_redundant(x: np.ndarray, A: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Remove columns whose rows are all covered twice, most expensive first."""
    x = x.copy()
    counts = A[:, x].sum(axis=1)                    # cover count per row
    for j in sorted(np.flatnonzero(x), key=lambda j: -c[j]):
        rows = A[:, j]
        if (counts[rows] >= 2).all():
            x[j] = False
            counts = counts - rows
    return x


def greedy_repair(x: np.ndarray, A: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Deterministic ratio-greedy repair of a partial (uncovering) solution."""
    x = x.copy()
    uncovered = ~A[:, x].any(axis=1)
    while uncovered.any():
        gain = A[uncovered].sum(axis=0).astype(float)
        score = np.where((gain > 0) & ~x, c / np.maximum(gain, 1.0), np.inf)
        j = int(np.argmin(score))
        x[j] = True
        uncovered &= ~A[:, j]
    return x


def scp_local_search(x: np.ndarray, A: np.ndarray, c: np.ndarray) -> np.ndarray:
    """Redundancy elimination + first-improvement 1-out exchange."""
    x = drop_redundant(x, A, c)
    improved = True
    while improved:
        improved = False
        for j in np.flatnonzero(x):
            trial = x.copy()
            trial[j] = False
            trial = drop_redundant(greedy_repair(trial, A, c), A, c)
            if c[trial].sum() < c[x].sum() - 1e-9:
                x, improved = trial, True
                break
    return x

The full algorithm is the generic grasp driver plus closures over the instance. Always benchmark against the alpha = 0 extreme — deterministic greedy plus the same local search — because that baseline is what randomization must beat:

import numpy as np

A, c = generate_scp(m=100, n=300, density=0.03, seed=7)

greedy_only = scp_local_search(
    construct_scp(A, c, alpha=0.0, rng=np.random.default_rng(0)), A, c)

best, best_cost, curve = grasp(
    construct=lambda alpha, rng: construct_scp(A, c, alpha, rng),
    local_search=lambda x: scp_local_search(x, A, c),
    cost=lambda x: float(c[x].sum()),
    n_iters=200, alpha=0.15, seed=42)

assert scp_is_cover(best, A)
print(f"greedy+LS: {c[greedy_only].sum():.0f}   GRASP: {best_cost:.0f}")
# Expected: both feasible covers; greedy+LS gives 648, GRASP 633 (~2% lower) --
# the 200 randomized restarts escape the deterministic greedy's bias.

Worked Example 2: QAP GRASP with Path Relinking

The quadratic assignment problem (QAP): given symmetric flow matrix $F$ and distance matrix $D$ (zero diagonals), find a permutation $\pi$ (facility $\to$ location) minimizing

$$ \min_{\pi} ; \sum_{a=1}^{n} \sum_{b=1}^{n} F_{ab} , D_{\pi(a)\pi(b)}. $$

GRASP was among the first competitive heuristics for QAP (Li, Pardalos & Resende 1994, "A greedy randomized adaptive search procedure for the quadratic assignment problem"). The construction places high-interaction facilities first; the local search is 2-exchange descent with the standard $O(n)$ swap delta — maintain a full delta table for larger instances, see fitness-evaluation-and-caching:

import numpy as np


def generate_qap(n: int, seed: int) -> tuple[np.ndarray, np.ndarray]:
    """Symmetric random QAP: integer flows F, Euclidean distances D, zero diagonals."""
    rng = np.random.default_rng(seed)
    F = np.triu(rng.integers(0, 10, (n, n)), k=1).astype(float)
    F = F + F.T
    pts = rng.random((n, 2)) * 100.0
    D = np.hypot(pts[:, None, 0] - pts[None, :, 0],
                 pts[:, None, 1] - pts[None, :, 1])
    return F, D


def qap_cost(perm: np.ndarray, F: np.ndarray, D: np.ndarray) -> float:
    """sum_{a,b} F[a,b] * D[perm[a], perm[b]]  (perm maps facility -> location)."""
    return float((F * D[np.ix_(perm, perm)]).sum())


def construct_qap(F: np.ndarray, D: np.ndarray, alpha: float,
                  rng: np.random.Generator) -> np.ndarray:
    """Greedy randomized placement, high-interaction facilities first.

    Greedy score of (facility f, free location l): the interaction cost
    between f at l and all already-placed facilities.
    """
    n = F.shape[0]
    perm = np.full(n, -1, dtype=int)
    free = np.ones(n, dtype=bool)                  # free locations
    placed: list[int] = []
    for f in np.argsort(-F.sum(axis=1)):           # facility insertion order
        locs = np.flatnonzero(free)
        if placed:
            pl = np.array(placed)
            inc = 2.0 * (D[np.ix_(locs, perm[pl])] @ F[f, pl])  # symmetric F, D
            cutoff = inc.min() + alpha * (inc.max() - inc.min())
            loc = int(rng.choice(locs[inc <= cutoff + 1e-12]))
        else:
            loc = int(rng.choice(locs))
        perm[f] = loc
        free[loc] = False
        placed.append(f)
    return perm


def swap_delta(perm: np.ndarray, F: np.ndarray, D: np.ndarray,
               a: int, b: int) -> float:
    """O(n) cost change for swapping the locations of facilities a and b."""
    mask = np.ones(len(perm), dtype=bool)
    mask[[a, b]] = False
    pk = perm[mask]
    u, v = perm[a], perm[b]
    return float(2.0 * ((F[a, mask] - F[b, mask]) * (D[v, pk] - D[u, pk])).sum())


def qap_local_search(perm: np.ndarray, F: np.ndarray,
                     D: np.ndarray) -> np.ndarray:
    """Best-improvement 2-exchange descent using O(n) swap deltas."""
    perm = perm.copy()
    n = len(perm)
    while True:
        best_delta, best_pair = -1e-9, None
        for a in range(n - 1):
            for b in range(a + 1, n):
                d = swap_delta(perm, F, D, a, b)
                if d < best_delta:
                    best_delta, best_pair = d, (a, b)
        if best_pair is None:
            return perm
        a, b = best_pair
        perm[a], perm[b] = perm[b], perm[a]

Path relinking turns memoryless GRASP into a method with intensification (Laguna & Martí 1999, "GRASP and path relinking for 2-layer straight line crossing minimization"). Keep a small elite set; after each GRASP iteration, walk step by step from the new local optimum toward a randomly chosen elite, always applying the best repairing swap, and keep the best intermediate solution on the path. The sketch below is the standard greedy forward relinking for permutations; a deeper treatment of reference-set management and relinking variants lives in scatter-search-path-relinking:

import numpy as np


def path_relink(start: np.ndarray, guide: np.ndarray, F: np.ndarray,
                D: np.ndarray) -> tuple[np.ndarray, float]:
    """Greedy walk from start toward guide; returns the best intermediate."""
    s = start.copy()
    best, best_cost = s.copy(), qap_cost(s, F, D)
    cur_cost = best_cost
    while True:
        diff = np.flatnonzero(s != guide)
        if len(diff) < 2:
            return best, best_cost
        moves = []
        for a in diff:
            b = int(np.flatnonzero(s == guide[a])[0])   # swap fixes position a
            moves.append((swap_delta(s, F, D, a, b), a, b))
        d, a, b = min(moves)                            # best repairing swap
        s[a], s[b] = s[b], s[a]
        cur_cost += d
        if cur_cost < best_cost - 1e-9:
            best, best_cost = s.copy(), cur_cost


def grasp_pr_qap(F: np.ndarray, D: np.ndarray, n_iters: int = 60,
                 alpha: float = 0.3, elite_size: int = 5,
                 seed: int = 0) -> tuple[np.ndarray, float]:
    """GRASP + path relinking to a random elite each iteration."""
    rng = np.random.default_rng(seed)
    elites: list[tuple[float, np.ndarray]] = []
    best, best_cost = None, np.inf

    def submit(perm: np.ndarray, cost: float) -> None:
        nonlocal best, best_cost
        if cost < best_cost:
            best, best_cost = perm.copy(), cost
        if all(np.any(perm != e) for _, e in elites):   # keep elites distinct
            elites.append((cost, perm.copy()))
            elites.sort(key=lambda t: t[0])
            del elites[elite_size:]

    for _ in range(n_iters):
        x = qap_local_search(construct_qap(F, D, alpha, rng), F, D)
        cx = qap_cost(x, F, D)
        if elites:
            _, guide = elites[int(rng.integers(len(elites)))]
            y, _ = path_relink(x, guide, F, D)
            y = qap_local_search(y, F, D)               # polish best intermediate
            cy = qap_cost(y, F, D)
            if cy < cx:
                submit(y, cy)
        submit(x, cx)
    return best, best_cost


F, D = generate_qap(n=12, seed=3)
best, best_cost = grasp_pr_qap(F, D, n_iters=40, alpha=0.3, elite_size=5, seed=1)
assert sorted(best.tolist()) == list(range(12))         # valid permutation
print(f"{best_cost:.1f}")
# Expected: 24650.1. Plain multi-start GRASP with the same seed and budget
# stalls at 24713.5; relinking against the elite set recovers the better
# value. Across seeds, GRASP+PR results here stay within 0.1% of each other.

Notes on the relinking step. Each swap fixes at least one disagreeing position, so a path between solutions at Hamming distance $h$ has at most $h - 1$ interior points; evaluating each step costs $O(h \cdot n)$ with the delta function. Apply full local search only to the best intermediate, not to every point on the path — relinking is a cheap trajectory scan, not a second descent phase.

Advanced Techniques

Reactive GRASP

Fixing alpha is the most common tuning mistake: the best value differs across instances and even across phases of a run. Reactive GRASP (Prais & Ribeiro 2000, "Reactive GRASP: an application to a matrix decomposition problem in TDMA traffic assignment") draws alpha from a discrete set with probabilities updated from observed solution quality — alphas that produced better averages get sampled more often:

import numpy as np


class ReactiveAlpha:
    """Reactive GRASP alpha selection (Prais & Ribeiro 2000).

    Samples alpha from a discrete set; every block_size iterations the
    sampling probabilities are reset proportional to
    (best_cost / avg_cost_of_alpha) ** beta.
    """

    def __init__(self, alphas: np.ndarray, block_size: int = 50,
                 beta: float = 10.0) -> None:
        self.alphas = alphas
        self.block_size = block_size
        self.beta = beta
        self.probs = np.full(len(alphas), 1.0 / len(alphas))
        self.sums = np.zeros(len(alphas))
        self.counts = np.zeros(len(alphas))
        self.best = np.inf
        self.iters = 0
        self.idx = 0

    def sample(self, rng: np.random.Generator) -> float:
        """Draw the alpha to use for the next GRASP iteration."""
        self.idx = int(rng.choice(len(self.alphas), p=self.probs))
        return float(self.alphas[self.idx])

    def report(self, cost: float) -> None:
        """Feed back the cost of the iteration run with the sampled alpha."""
        self.sums[self.idx] += cost
        self.counts[self.idx] += 1
        self.best = min(self.best, cost)
        self.iters += 1
        if self.iters % self.block_size == 0:
            avg = self.sums / np.maximum(self.counts, 1.0)
            avg[self.counts == 0] = self.best        # untried alphas stay attractive
            q = (self.best / avg) ** self.beta
            self.probs = q / q.sum()


# Reactive GRASP on the set covering functions from Worked Example 1.
A, c = generate_scp(m=100, n=300, density=0.03, seed=7)
ra = ReactiveAlpha(alphas=np.linspace(0.0, 0.9, 10), block_size=40, beta=50.0)
rng = np.random.default_rng(5)
best_cost = np.inf
for _ in range(200):
    alpha = ra.sample(rng)
    x = scp_local_search(construct_scp(A, c, alpha, rng), A, c)
    ra.report(float(c[x].sum()))
    best_cost = min(best_cost, float(c[x].sum()))
print(best_cost, np.round(ra.probs, 2))
# Expected: best_cost 633 (matches the fixed-alpha run above); probabilities
# concentrate on the low alphas, e.g. ~0.27 and ~0.24 on alpha 0.0 and 0.1.

Prais & Ribeiro used beta = 10. Note, however, that a strong local search compresses the cost differences between alphas — every start descends to a similar level — so a sharper exponent (30–100) is often needed before probability mass moves decisively. Report the final probability vector in experiments: it is free diagnostic output that tells you what alpha the instance class actually wants.

Bias Functions in RCL Sampling

Uniform sampling inside the RCL is the classic choice, but the selection can be biased toward better-ranked candidates without shrinking the list (Bresina 1996, "Heuristic- biased stochastic sampling"). This decouples diversity (RCL size, via alpha) from greediness (sampling sharpness, via the bias):

import numpy as np


def biased_rcl_choice(scores: np.ndarray, alpha: float, bias: str,
                      rng: np.random.Generator) -> int:
    """RCL sampling with a rank bias: 'uniform', 'linear' (1/rank),
    or 'exponential' (e^-rank). Lower score = better candidate."""
    feasible = np.flatnonzero(np.isfinite(scores))
    s = scores[feasible]
    rcl = feasible[s <= s.min() + alpha * (s.max() - s.min()) + 1e-12]
    rank = np.argsort(np.argsort(scores[rcl])) + 1.0
    w = {"uniform": np.ones_like(rank),
         "linear": 1.0 / rank,
         "exponential": np.exp(-rank)}[bias]
    return int(rng.choice(rcl, p=w / w.sum()))

A linear bias with a generous alpha (0.4–0.6) often matches a tuned uniform RCL while being less sensitive to the exact alpha value — a useful robustness trade.

Path Relinking Variants

Forward relinking (new solution toward elite, as implemented above) is the default. Backward relinking starts from the elite — usually better, since more of the path lies near the higher-quality endpoint. Mixed relinking alternates ends and meets in the middle. Truncated relinking explores only the first 20–40% of the path, where most improvements concentrate. Evolutionary path relinking periodically relinks all elite pairs as a post-processing intensification (Resende & Werneck 2004, "A hybrid heuristic for the p-median problem"). Implement forward + truncated first; measure before adopting anything fancier. Reference-set admission rules (quality plus minimum distance) are covered in scatter-search-path-relinking.

Avoiding Repeated Work Across Restarts

With small alpha, restarts collide: identical or near-identical constructions get re-descended at full cost. Hash each constructed solution (e.g. hash(x.tobytes()) for a canonical array form) and skip or re-randomize duplicates before local search. On problems where local search is the bottleneck this saves 10–30% of runtime at negligible cost. Memoization of full evaluations belongs to fitness-evaluation-and-caching.

Cost Perturbations and Sampled Greedy

Two construction variants matter in practice. Cost perturbation multiplies greedy scores by random factors (e.g. uniform on [1, 1.25]) and then runs a pure greedy construction — diversity comes from the data instead of the RCL; this often preserves more greedy quality than value RCLs (Resende & Werneck 2004). Sampled greedy evaluates only a random subset of candidates per step and picks the best of the sample — the method of choice when $|C|$ is huge and scoring every candidate per step is the bottleneck. Both plug into the same grasp driver unchanged.

Practical Challenges

Every iteration returns nearly the same solution. Alpha is too small, or greedy scores contain heavy ties broken deterministically. Verify the construction actually consumes randomness (log RCL sizes per step — size 1 means pure greedy), randomize tie-breaking, and raise alpha or switch to a cardinality RCL when score ranges collapse.

The incumbent stops improving after a small fraction of the budget. Expected behavior for memoryless sampling: the running minimum of i.i.d. draws improves ever more slowly. Add memory — reactive alpha, path relinking — or split the remaining budget across parallel seeds instead of extending a stagnant run.

Construction dominates runtime. Update greedy scores incrementally instead of recomputing them: in set covering, only columns intersecting newly covered rows change their gain. Use np.argpartition for cardinality RCLs, and switch to sampled greedy when the candidate set is very large.

Local search barely improves constructed solutions. Either alpha is so low that constructions are already locally optimal for the chosen neighborhood, or the neighborhood is too weak. Track mean relative improvement of the descent phase; if it is near zero, strengthen the neighborhood (see local-search-and-neighborhoods) before touching alpha.

Construction dead-ends on tightly constrained problems. When no candidate is feasible for a partial solution, do not silently return garbage: restart the construction (cheap), allow temporarily infeasible additions with penalty scores, or finish with a repair operator. Count dead-ends — a high rate means the greedy order itself needs rethinking.

Path relinking elite set collapses to near-duplicates. Admission by quality alone fills the set with clones and relinking paths of length zero. Admit a solution only if it differs from every elite by a minimum distance (e.g. Hamming distance ≥ n/10), and replace the closest worse elite rather than the worst one.

Alpha tuned on one instance fails on the rest. Tune on a training subset and report on a disjoint test subset, or sidestep the issue with reactive GRASP. Always report the alpha sensitivity curve — reviewers ask for it, and flat curves are a selling point.

Comparisons against ILS or simulated annealing look unfair. Give every method the same wall-clock or evaluation budget and the same local search where applicable, run 10+ seeds, and compare distributions, not single best values. GRASP's variance across seeds is usually low; exploit that in the write-up.

Tools & Libraries

Output Format

A complete GRASP deliverable contains:

1. Configuration summary — one table, so the run is reproducible:

2. Results table — per instance: best cost, mean ± std over seeds, gap to best-known or bound, time-to-best, iterations-to-best, and the alpha = 0 (greedy + LS) baseline cost in its own column.
3. Validation statement — every reported solution passed an independent feasibility check and an objective recomputation outside the algorithm code.
4. Convergence evidence — best-so-far curve (median over seeds with a min-max band) and, for runtime claims, a TTT plot at a stated target value.
5. Diagnostics — final reactive-alpha probabilities (if used), mean local-search improvement per iteration, duplicate-construction rate.
6. Code artifacts — instance generator with seeds, the algorithm, and the validator as separate functions, so reviewers can rerun everything.

Report iteration counts and wall-clock time: GRASP comparisons are meaningless under iteration counts alone when construction costs differ across methods.

Questions to Ask

• What greedy construction already exists for this problem, and what is its score?
• Are instances small enough that an exact solver closes the gap in the time budget?
• What is the wall-clock budget per instance, and how many cores are available?
• Is there a known neighborhood with cheap delta evaluation for the local search?
• Fixed alpha or reactive — is cross-instance robustness a requirement?
• Is path relinking in scope, or is plain multi-start enough for the deliverable?
• What baselines must the comparison include, and at what budgets and seed counts?
• How will solutions be validated independently of the algorithm code?
• Are there benchmark instances with best-known values to report gaps against?

Related Skills

• local-search-and-neighborhoods — when the improvement phase needs better neighborhoods, first-vs-best improvement decisions, or delta evaluation design.
• warm-starts-and-initial-solutions — when GRASP constructions feed a MIP start or another solver, or you need construction heuristics by problem class.
• set-covering-packing-partitioning — when the underlying problem is covering or partitioning and you need exact models, bounds, or Lagrangian alternatives.
• scatter-search-path-relinking — when path relinking becomes the main intensification tool and you need reference-set management and relinking variants in depth.