Phase 06

Heuristics & Optimization

When exact search is too slow: bounds, local search, annealing, populations, black-box budgets and game bots.

Phase 6: Heuristics & Optimization

Phase 5 ended with backtracking, and with a warning: it is exponential. This phase starts where that warning bites. Many of the most useful problems in the world are like that: plan the shortest delivery route, pack parcels into the fewest vans, build a timetable without clashes, choose the settings of a machine you cannot see inside. For none of them is there a known algorithm that is both fast and always exact, and most computer scientists believe there never will be.

So we change the question. Instead of "what is the best answer?" we ask "how good an answer can I find in the time I have?". The tools for that are heuristics: methods that usually find very good answers quickly, without a proof that they are the best. You will learn to build a solution greedily, improve it with small changes, escape the traps that small changes fall into, search with a whole population of candidates, optimise functions you can only evaluate, and write a bot that plans a few moves ahead against an opponent.

This phase is also scored differently. Most problems have no single right answer; the next section explains how they are graded.

When exact search is too slow

A closed route through n stops can be driven in (n - 1)! / 2 different orders (fix the start, and each loop can be driven in two directions). That number explodes:

import math
for n in (5, 10, 15, 20, 25):
    print(n, math.factorial(n - 1) // 2)
5 12
10 181440
15 43589145600
20 60822550204416000
25 310224200866619719680000

Trying every order is fine for 10 stops and hopeless for 20. Problems like this are called NP-hard: nobody knows an algorithm whose running time grows only polynomially with the input size, and finding one (or proving there is none) is the most famous open question in computer science, "P versus NP". In practice NP-hard means: exact methods work for small inputs, and for large inputs you use heuristics.

Exact methods can still be pushed a long way. Dynamic programming over subsets remembers, for every set of visited stops and every last stop, the shortest way to get there. That is 2^n · n states instead of n! orders: about 5 million for 20 stops instead of 10^17.

import math
from functools import lru_cache

def shortest_loop(d):
    n = len(d)
    @lru_cache(maxsize=None)
    def best(mask, last):
        # shortest path from 0 through exactly the stops in mask, ending at last
        if mask == (1 << last) | 1:
            return d[0][last]
        prev = mask & ~(1 << last)
        return min(best(prev, k) + d[k][last] for k in range(1, n) if prev >> k & 1)
    full = (1 << n) - 1
    return min(best(full, k) + d[k][0] for k in range(1, n))

pts = [(0, 0), (4, 0), (4, 3), (0, 3), (2, 5)]
d = [[math.dist(a, b) for b in pts] for a in pts]
print(round(shortest_loop(d), 2))
15.66

The other classic exact technique is branch and bound: backtracking plus a bound, an optimistic estimate of the best answer a partial solution could still reach. When even the optimistic estimate cannot beat the best answer found so far, the whole branch is skipped. For a knapsack, a good bound fills the remaining room greedily by value per weight and is allowed to take a fraction of the last item:

def knapsack_bb(items, cap):
    # items: (value, weight); sorted by value per weight for a good bound
    items = sorted(items, key=lambda it: it[0] / it[1], reverse=True)
    best = 0
    def bound(i, value, room):
        # fractional relaxation: fill the rest greedily, splitting the last item
        for v, w in items[i:]:
            if w <= room:
                value, room = value + v, room - w
            else:
                return value + v * room / w
        return value
    def search(i, value, room):
        nonlocal best
        best = max(best, value)
        if i == len(items) or bound(i, value, room) <= best:
            return                      # prune: even the optimistic bound cannot win
        v, w = items[i]
        if w <= room:
            search(i + 1, value + v, room - w)
        search(i + 1, value, room)
    search(0, 0, cap)
    return best

print(knapsack_bb([(60, 10), (100, 20), (120, 30)], 50))
220

Bounds are useful even when you do not search: a lower bound on a minimisation problem tells you how far from the best answer you can possibly be. The weight of a minimum spanning tree is a lower bound for any closed route (remove one road from a route and you have a spanning tree), and the total size of all parcels divided by the van capacity, rounded up, is a lower bound for the number of vans. If your heuristic lands on the bound, you know it is optimal.

Common mistakes: a bound that is not really optimistic (it underestimates what a branch can reach) prunes good branches and gives wrong answers. Bitmask DP with n above about 20 runs out of memory. And do not assume greedy is optimal because it looks reasonable: always look for a counterexample first.

Practice: Try Every Order, No Loop Can Be Shorter Than This, Fifteen Stops, One-Way Streets, Cargo Too Heavy to Tabulate

How these problems are scored

A heuristic problem has many acceptable answers of different quality, so it is graded in two parts.

  • Pass or fail. The answer must be feasible (a real route that visits every stop once, bins that are not overfull, a legal move) and at least as good as a simple baseline that the grader computes itself, such as the route that always drives to the nearest stop, or the best of as many random guesses as you were allowed evaluations. Each problem's description says what its baseline is.
  • Quality from 0 to 100. How much of the gap between the baseline and the best known answer you close. 0 means baseline-level, 100 means you matched (or beat) the best answer we know. For continuous functions the scale is logarithmic (every factor of ten closer to the optimum counts the same); for games it is the share of points you took.

The problem header shows the par quality: the quality of the reference solution. One star for a pass, two stars at three quarters of par, three stars for matching par without hints. In the arena you can rank a board by Quality.

Two rules make this fair. Every test starts from the same random state, so a randomised method gives the same result on every run and every device (the arena leader re-runs your code and must see the same quality). And you must limit your loops by a number of rounds, not by the clock: "run for one second" gives different answers on a fast laptop and a slow phone, and the tests have time limits.

Construct, then improve

Most heuristics have two phases. First construct a complete, reasonable solution quickly, usually greedily: always drive to the nearest unvisited stop, put each parcel into the first van it fits, give each node the smallest colour its neighbours are not using. Then improve it with small changes, a neighbourhood: swap two items, move one item elsewhere, reverse a stretch of a route, flip one node to the other side. Hill climbing applies improving moves until none is left.

For routes the classic move is 2-opt: take two roads (a, b) and (c, e), replace them with (a, c) and (b, e), which means reversing the stretch between them. Two roads that cross can always be uncrossed this way, and the route gets shorter.

import math, random

def length(pts, t):
    return sum(math.dist(pts[t[i]], pts[t[i - 1]]) for i in range(len(t)))

def two_opt(pts, tour):
    n = len(tour)
    improved = True
    while improved:
        improved = False
        for i in range(n - 1):
            for j in range(i + 2, n if i else n - 1):
                a, b, c, e = tour[i], tour[i + 1], tour[j], tour[(j + 1) % n]
                old = math.dist(pts[a], pts[b]) + math.dist(pts[c], pts[e])
                new = math.dist(pts[a], pts[c]) + math.dist(pts[b], pts[e])
                if new < old - 1e-9:
                    tour[i + 1:j + 1] = reversed(tour[i + 1:j + 1])
                    improved = True
    return tour

rng = random.Random(4)
pts = [(rng.randint(0, 100), rng.randint(0, 100)) for _ in range(12)]
tour = list(range(12))
print(round(length(pts, tour), 1))
print(round(length(pts, two_opt(pts, tour)), 1))
632.9
353.9

Three design choices matter more than anything else. Evaluate moves incrementally: 2-opt compares four road lengths instead of recomputing the whole route, which is the difference between O(1) and O(n) per move. Precompute what you look up constantly (a distance table). And pick a neighbourhood that is small enough to scan but rich enough to escape bad layouts.

When hill climbing stops, the solution is a local optimum: no single move improves it. That is not the same as the best solution. Different starts lead to different local optima, and the next section is about getting out of them.

Common mistakes: comparing floating-point lengths without a small tolerance, so the loop swaps back and forth forever. Recomputing the whole cost after every candidate move. Counting the same road twice in 2-opt (when i is 0, j must stop before the last position).

Practice: Closest Shelf First, Shorter Delivery Loop, Fewer Crates for the Parcels, Fewest Exam Slots

Escaping local optima

The simplest escape is random restarts: run hill climbing from many different starting points and keep the best result. It works surprisingly often, and it is the baseline every cleverer method must beat.

Simulated annealing accepts some worse moves on purpose. A move that makes the cost worse by delta is accepted with probability exp(-delta / T), where the temperature T starts high (almost everything is accepted, the search wanders) and cools slowly (only small setbacks are accepted, then none). The name comes from metallurgy: cooling metal slowly lets its atoms settle into a low-energy arrangement.

import math, random

def anneal(state, cost, neighbour, steps=20000, t0=50.0, t1=0.01):
    cur, cur_cost = state, cost(state)
    best, best_cost = cur, cur_cost
    for k in range(steps):
        t = t0 * (t1 / t0) ** (k / steps)          # geometric cooling
        cand = neighbour(cur)
        delta = cost(cand) - cur_cost
        if delta <= 0 or random.random() < math.exp(-delta / t):
            cur, cur_cost = cand, cur_cost + delta
            if cur_cost < best_cost:
                best, best_cost = cur, cur_cost
    return best, best_cost

# split numbers into two groups with sums as close as possible
nums = [71, 38, 95, 12, 64, 27, 50, 83, 9, 46]
def cost(side):
    return abs(sum(x if s else -x for x, s in zip(nums, side)))
def flip_one(side):
    i = random.randrange(len(side))
    return side[:i] + (not side[i],) + side[i + 1:]

random.seed(1)
start = tuple([True] * len(nums))
print(cost(start), anneal(start, cost, flip_one)[1])
495 1

The starting temperature should be about the size of a typical bad move: here the numbers are around 50, so t0 = 50. Start much colder and annealing is just hill climbing; much hotter and it wastes its time wandering. Always keep the best state seen, because the current state may drift away from it.

Tabu search takes a different route: it always makes the best move in the neighbourhood, even a worsening one, but remembers the last few moves in a tabu list and forbids undoing them for a while, so it cannot fall straight back into the optimum it just left. Min-conflicts is the version for puzzles with constraints (n queens, timetables): repeatedly pick something that is in conflict and move it to the place with the fewest conflicts.

Common mistakes: computing exp(-delta / t) with t reaching 0 (divide by zero); cool to a small positive value. Returning the final state instead of the best one. A neighbour function that sometimes returns an infeasible state.

Practice: Summit in the Fog, Queens on a Cracked Board, A Seating Plan Everyone Can Live With, House Rules for the Smart Home

Populations

The methods so far keep one candidate. Population methods keep many and let them share information.

A genetic algorithm (GA) imitates evolution. Each candidate is scored by its fitness. Each generation, selection picks parents that tend to be fit (a tournament: draw a few at random, keep the best), crossover combines two parents into a child, mutation makes a small random change, and elitism copies the best few unchanged into the next generation, so the best answer is never lost. For bit strings, one-point crossover takes the start of one parent and the end of the other. For permutations (routes, schedules) that would create duplicates, so an order crossover keeps a slice of one parent and fills the rest in the other parent's order:

import random

def order_crossover(p1, p2, rng):
    # keep a slice of parent 1, fill the rest in parent 2's order
    n = len(p1)
    i, j = sorted(rng.sample(range(n), 2))
    middle = p1[i:j]
    rest = [x for x in p2 if x not in middle]
    return rest[:i] + middle + rest[i:]

rng = random.Random(3)
print(order_crossover([0, 1, 2, 3, 4, 5, 6, 7], [7, 6, 5, 4, 3, 2, 1, 0], rng))
[7, 6, 5, 3, 4, 2, 1, 0]

Particle swarm optimisation (PSO) is for continuous problems: each particle has a position and a velocity, and each step its velocity is pulled towards the best point it has seen and the best point any particle has seen, plus some inertia. Ant colony optimisation (ACO) builds routes step by step: ants choose the next stop with a probability that grows with the pheromone on that road and shrinks with its length; after each round good routes deposit pheromone and all trails evaporate a little, so the colony gradually concentrates on short roads.

Population methods need more evaluations than hill climbing per improvement, so they shine when the landscape is rugged (many local optima) or when good solutions share building blocks that crossover can combine. On smooth problems a simple local search is often better, and it is worth checking.

Common mistakes: no elitism, so the best solution is lost to a bad crossover. A mutation rate so high that children are random, or so low that the population becomes identical after a few generations (loss of diversity). Evaluating the same candidate's fitness over and over instead of storing it.

Practice: One Generation of Seed Breeding, Hut-to-Hut Hiking Loop, Lowest Point on Rough Ground, Print Shop Queue with Ink Changes

Black-box optimization with a budget

Sometimes you cannot see the problem at all: you can only try an input and read how good it is. A simulator, a physical experiment, a machine learning model's training run. Each evaluation is expensive, so you get a budget: at most so many calls. The whole game is to spend them well.

Random search (try random points, keep the best) is the baseline. A far better simple method is the (1+1) evolution strategy: keep one point, try a random step from it, keep the step if it is better, and adapt the step size with the one-fifth success rule: if more than about one try in five succeeds, the steps are too timid, so grow them; if fewer, shrink them.

import random

def one_plus_one(f, x, budget, step=1.0):
    fx = f(x)
    wins = tries = 0
    for _ in range(budget - 1):
        y = [xi + random.gauss(0, step) for xi in x]
        fy = f(y)
        tries += 1
        if fy < fx:
            x, fx, wins = y, fy, wins + 1
        if tries == 10:                      # the one-fifth success rule
            step *= 1.5 if wins > 2 else 0.6
            wins = tries = 0
    return x, fx

random.seed(2)
bowl = lambda x: sum((xi - 3) ** 2 for xi in x)
x, fx = one_plus_one(bowl, [0.0, 0.0, 0.0], 300)
print(f"{fx:.2e}")
1.04e-08

Three hundred evaluations bring a 3-dimensional bowl within 10^-8 of its minimum; random search with the same budget would still be off by about 1. Curved valleys and landscapes with many dips are harder: there a swarm, the Nelder-Mead simplex method or covariance adaptation (which learns the shape of good steps, not just their size) pay off. The scores in this section are logarithmic because improvements come in factors of ten.

Common mistakes: calling the function again for a point you already evaluated (store the value). Losing track of how many evaluations you have used. Stepping outside the allowed box (clamp every coordinate).

Practice: Sliders That Do Not Interact, Tuning Knobs in the Dark, Readings You Cannot Trust

Adversarial search

In a game your "neighbourhood" is chosen alternately by you and an opponent who wants the opposite. Minimax assumes both play perfectly: the value of a position is the best of your moves if it is your turn, and the worst for you if it is the opponent's. For small games you can compute it exactly with memoised recursion. Here is a take-away game (remove 1, 2 or 3 stones; whoever takes the last stone wins), written in the compact negamax form: my value is the negative of my opponent's value after my move.

from functools import lru_cache

# a take-away game: remove 1, 2 or 3 stones; whoever takes the last stone wins
@lru_cache(maxsize=None)
def value(stones):
    # +1 if the player to move wins with best play, -1 if they lose
    if stones == 0:
        return -1                             # the previous player took the last stone
    return max(-value(stones - k) for k in (1, 2, 3) if k <= stones)

print([value(s) for s in range(1, 10)])
[1, 1, 1, -1, 1, 1, 1, -1, 1]

The pattern shows the game's secret: multiples of four lose for the player to move.

Real games are too big to search to the end, so you search a fixed number of moves ahead (the depth, measured in plies, single moves) and score the positions at the bottom with an evaluation function: a heuristic guess such as "open lines of three are good for me". Alpha-beta pruning makes the same decision as minimax while skipping branches that cannot matter: alpha is the best value you are already guaranteed, beta the best the opponent is guaranteed, and once alpha >= beta the rest of the branch is irrelevant.

import math

def alphabeta(node, depth, alpha, beta, maximising, children, score):
    kids = children(node)
    if depth == 0 or not kids:
        return score(node)
    if maximising:
        best = -math.inf
        for kid in kids:
            best = max(best, alphabeta(kid, depth - 1, alpha, beta, False, children, score))
            alpha = max(alpha, best)
            if alpha >= beta:
                break                         # the opponent will never allow this branch
        return best
    best = math.inf
    for kid in kids:
        best = min(best, alphabeta(kid, depth - 1, alpha, beta, True, children, score))
        beta = min(beta, best)
        if alpha >= beta:
            break
    return best

tree = {"root": ["a", "b"], "a": ["a1", "a2"], "b": ["b1", "b2"]}
leaf = {"a1": 3, "a2": 5, "b1": 2, "b2": 9}
print(alphabeta("root", 2, -math.inf, math.inf, True, lambda n: tree.get(n, []), lambda n: leaf.get(n, 0)))
3

In that tree the opponent answers a with a1 (value 3). After seeing b1 (value 2) under b, alpha-beta never looks at b2: the opponent can already hold b to 2, which is worse for you than the 3 you have. Pruning works best when good moves are tried first, so order moves by a quick guess (in four in a row: central columns first).

Common mistakes: scoring a won position the same no matter how soon the win comes, so the bot dawdles; add the remaining depth to a win. Evaluating from the wrong player's point of view in negamax. Returning no move when every move loses; always return a legal one.

Practice: Two Jars of Marbles, How Much Does the Scout Skip?, Four in a Row: Beat the House, Six-by-Six Reversi: Outflank the House

Checklist

Before moving on, make sure you can do each of these without looking anything up:

  • Explain why trying every order is hopeless for 20 stops, and what NP-hard means in practice.
  • Write a bitmask DP over (set of visited items, last item) and say when it is still feasible.
  • Add an optimistic bound to a backtracking search and prune with it; give a lower bound for a route or a packing.
  • Build a greedy solution, then improve it with a neighbourhood move evaluated incrementally.
  • Implement 2-opt and explain what a local optimum is.
  • Write simulated annealing with a cooling schedule, choose a starting temperature and keep the best state.
  • Describe tabu search and min-conflicts in one sentence each.
  • Implement a genetic algorithm with tournament selection, crossover, mutation and elitism; use order crossover for permutations.
  • Describe how a particle swarm and an ant colony share information.
  • Spend an evaluation budget with an adaptive step size, and never exceed it.
  • Write minimax (or negamax) with alpha-beta pruning, a depth limit and an evaluation function.
  • Keep every result reproducible: seeded randomness and loops bounded by rounds, not time.