Edit Distance, Alignment & Approximate Matching

Edit Metrics and Dynamic Programming

A string distance assigns a cost to transforming one sequence into another. The chosen operations determine both the meaning of the distance and the recurrence needed to compute it. Hamming distance compares aligned positions and therefore requires equal-length strings. Levenshtein distance permits insertion, deletion, and substitution, while weighted variants assign application-specific costs to those operations. Damerau–Levenshtein distance additionally permits transposing adjacent symbols.

For strings $A$ and $B$, define $D[i][j]$ as the minimum cost of transforming $A[0 \dots i)$ into $B[0 \dots j)$. Under unit-cost Levenshtein distance, $D[i][0]=i$, $D[0][j]=j$, and

$$ D[i][j]=\min\left(D[i-1][j]+1,\ D[i][j-1]+1,\ D[i-1][j-1]+[A[i-1]\ne B[j-1]]\right). $$

Weighted distance replaces unit costs with operation-dependent functions. Metric guarantees require non-negative symmetric costs, identity of indiscernibles, and triangle-compatible costs. Restricted optimal-string-alignment distance adds $D[i-2][j-2]$ for a reversed adjacent pair but cannot edit the same substring twice; full Damerau–Levenshtein uses last-occurrence information to support interacting transpositions.

Alignment Reconstruction and Acceleration

The full $O(nm)$ table supports traceback. From $(n,m)$, a diagonal move aligns two symbols, an upward move deletes a source symbol, and a leftward move inserts a target symbol. Reversing these choices produces an optimal alignment. When only the distance is needed, rolling rows reduce auxiliary space to $O(\min(n,m))$.

If the distance is at most $k$, a valid path cannot leave the band $|i-j|\le k$. Banded DP evaluates $O(k)$ cells per row, while thresholded Ukkonen search stops when no state can satisfy the bound. Myers' bit-vector algorithm packs many DP states into machine words and is especially effective for short patterns.

from math import inf


def edit_distance_with_traceback(source: str, target: str):
    """Return distance, one aligned pair, and the edit operations."""
    n, m = len(source), len(target)
    dp = [[0] * (m + 1) for _ in range(n + 1)]
    parent = [[None] * (m + 1) for _ in range(n + 1)]

    for i in range(1, n + 1):
        dp[i][0], parent[i][0] = i, "delete"
    for j in range(1, m + 1):
        dp[0][j], parent[0][j] = j, "insert"

    for i in range(1, n + 1):
        for j in range(1, m + 1):
            candidates = [
                (dp[i-1][j-1] + (source[i-1] != target[j-1]),
                 "match" if source[i-1] == target[j-1] else "substitute"),
                (dp[i-1][j] + 1, "delete"),
                (dp[i][j-1] + 1, "insert"),
            ]
            dp[i][j], parent[i][j] = min(candidates, key=lambda item: item[0])

    top, bottom, operations = [], [], []
    i, j = n, m
    while i > 0 or j > 0:
        operation = parent[i][j]
        if operation in {"match", "substitute"}:
            top.append(source[i-1]); bottom.append(target[j-1])
            i -= 1; j -= 1
        elif operation == "delete":
            top.append(source[i-1]); bottom.append("-"); i -= 1
        elif operation == "insert":
            top.append("-"); bottom.append(target[j-1]); j -= 1
        else:
            raise RuntimeError("invalid traceback state")
        operations.append(operation)

    return (dp[n][m], "".join(reversed(top)), "".join(reversed(bottom)),
            list(reversed(operations)))


def banded_edit_distance(source: str, target: str, k: int) -> int | None:
    """Return the distance if it is at most k; otherwise return None."""
    if k < 0:
        raise ValueError("k must be non-negative")
    n, m = len(source), len(target)
    if abs(n - m) > k:
        return None

    previous = {j: j for j in range(min(m, k) + 1)}
    for i in range(1, n + 1):
        lower, upper = max(0, i-k), min(m, i+k)
        current = {0: i} if lower == 0 else {}
        for j in range(max(1, lower), upper + 1):
            value = min(
                previous.get(j, inf) + 1,
                current.get(j-1, inf) + 1,
                previous.get(j-1, inf) + (source[i-1] != target[j-1]),
            )
            if value <= k:
                current[j] = int(value)
        if not current:
            return None
        previous = current

    result = previous.get(m)
    return result if result is not None and result <= k else None


if __name__ == "__main__":
    distance, top, bottom, operations = edit_distance_with_traceback("kitten", "sitting")
    print(distance, top, bottom, operations, sep="\n")
    assert distance == 3
    assert banded_edit_distance("kitten", "sitting", 3) == 3
    assert banded_edit_distance("kitten", "sitting", 2) is None

The LCS Family and Diff

Let $L[i][j]$ be the LCS length of $A[0 \dots i)$ and $B[0 \dots j)$. Matching final symbols extend $L[i-1][j-1]$; otherwise take $\max(L[i-1][j],L[i][j-1])$. If only insertions and deletions cost one, an LCS of length $\ell$ yields distance $|A|+|B|-2\ell$ and shortest-common-supersequence length $|A|+|B|-\ell$.

Standard DP costs $O(nm)$ time and space. Hirschberg reconstructs an LCS in linear auxiliary space by combining forward and reverse rows. Hunt–Szymanski processes matching-position lists and excels when the number of matching pairs is small. Bitset algorithms update many states per machine operation. Diff tools turn an LCS backbone into insertions and deletions; tie-breaking among optimal witnesses affects readability.

flowchart TD
    INPUT["Strings A and B"]
    LCS["LCS length = l"]
    KEEP["Keep l shared symbols"]
    EDITS["Insert/delete distance = |A| + |B| - 2l"]
    SCS["Shortest common supersequence = |A| + |B| - l"]
    INPUT --> LCS
    LCS --> KEEP
    KEEP -->|"delete |A| - l; insert |B| - l"| EDITS
    KEEP -->|"merge unmatched symbols"| SCS

Global and Local Sequence Alignment

Needleman–Wunsch computes a global alignment, forcing both sequences to participate from beginning to end. Smith–Waterman computes a local alignment by allowing states to restart at zero; traceback begins at the maximum cell and stops at zero. Semi-global alignment makes selected terminal gaps free, which suits fitting a short read into a longer reference.

With affine gaps, a length-$k$ gap costs $g_o+(k-1)g_e$. Gotoh's algorithm retains $O(nm)$ time using three states: $M$ ends in a symbol pair, $X$ ends in a gap in the second sequence, and $Y$ ends in a gap in the first.

Approximate Search Filters and Indexes

Approximate Bitap extends Shift-Or with one bit-vector level per error and works well for short patterns. Myers' bit-vector method similarly accelerates Levenshtein matching. A $q$-gram filter rejects candidates with insufficient fragment overlap; pigeonhole filtering splits a pattern into $k+1$ pieces so that an alignment with at most $k$ errors must preserve an appropriate seed. Seed-and-extend searches seeds first and verifies only nearby candidates.

BK-trees organize strings under a metric: after measuring distance $d$ to a node, a radius-$k$ query visits only child edges in $[d-k,d+k]$. Levenshtein automata recognize all strings within distance $k$ of a fixed pattern and can be intersected with a trie or finite-state dictionary. Filters must be followed by exact verification unless their construction itself proves the final guarantee.