Longest Common Subsequence, in depth: the recurrence, traceback, linear-space Hirschberg and how diff uses it

The problem, and why brute force fails

Given X of length n and Y of length m, find a longest sequence Z that is a subsequence of both. Note the article: there can be several. For X = ABCBDAB and Y = BDCABA, both BCBA and BDAB have length 4, and no common subsequence is longer.

Do not confuse subsequence with substring. A substring is contiguous; a subsequence is not. The longest common substring of the same pair is only length 2 (for example AB), and it needs a different recurrence. Mixing up the two is the single most common LCS bug. Pattern search over contiguous text is covered in string matching algorithms.

Brute force enumerates the 2ⁿ subsequences of X and checks each against Y in O(m). At n = 40 that is over a trillion candidates. The way out is the observation every dynamic program rests on: the answer for the whole problem can be assembled from answers to smaller problems of the same shape, and there are only polynomially many of those. If the general method is new to you, read dynamic programming in depth first.

Deriving the recurrence

Define L[i][j] as the LCS length of the prefixes X[0:i] and Y[0:j]. There are (n+1)(m+1) such subproblems. Now look only at the last characters, X[i-1] and Y[j-1].

• They match. Some LCS of the two prefixes ends with that character. Any common subsequence that does not use it can be changed to use it without getting shorter, by swapping its last element for this pair. So L[i][j] = L[i-1][j-1] + 1.
• They differ. No common subsequence can end with both of them, so at least one is unused. Drop X's last character or drop Y's, and keep the better result: L[i][j] = max(L[i-1][j], L[i][j-1]).
• Base case. An empty prefix has nothing in common with anything: L[0][j] = L[i][0] = 0.

Each cell depends only on the cell above, the cell to the left and the diagonal. So filling rows top to bottom and left to right always has its inputs ready. That gives O(nm) time and, for the full table, O(nm) space. The matching case is the one to check carefully in any proof or code review. Taking the diagonal when characters match is always safe, which is what makes the recurrence greedy-safe at that one step and saves a third comparison.

Worked example: filling the table

Here is the table for the running example, produced by the code in the next section. Row labels are X, column labels are Y, and the extra row and column of zeros are the empty-prefix base case. The shaded cells are the traceback path from the bottom-right corner; green cells are where the traceback took a character.

Read a few cells by hand to build intuition. At row A (i = 1), column A (j = 4), the characters match, so the cell is the diagonal (0) plus 1. At row C (i = 3), column C (j = 3), they match again: the diagonal holds 1 (from AB against BD, which share B), so the cell is 2. At row D (i = 5), column A (j = 4), they differ. The cell above holds 2 and the cell to the left holds 2, so the cell is 2. The final answer, 4, sits in the bottom-right corner.

Code: table and traceback

The table gives the length. To get an actual subsequence, walk back from the corner. On a match, emit the character and move diagonally. Otherwise, move to whichever neighbour holds the larger value. The tie-break, up or left, decides which of several optimal answers you get. Moving up on ties yields BCBA for the example. Moving left on ties can yield a different optimal subsequence of the same length.

def lcs_table(a, b):
    n, m = len(a), len(b)
    L = [[0] * (m + 1) for _ in range(n + 1)]
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            if a[i - 1] == b[j - 1]:
                L[i][j] = L[i - 1][j - 1] + 1
            else:
                L[i][j] = max(L[i - 1][j], L[i][j - 1])
    return L


def lcs_string(a, b):
    L = lcs_table(a, b)
    i, j, out = len(a), len(b), []
    while i > 0 and j > 0:
        if a[i - 1] == b[j - 1]:
            out.append(a[i - 1])
            i, j = i - 1, j - 1
        elif L[i - 1][j] >= L[i][j - 1]:   # tie: prefer moving up
            i -= 1
        else:
            j -= 1
    return "".join(reversed(out))


assert lcs_string("ABCBDAB", "BDCABA") == "BCBA"

The traceback is O(n + m), because every step decrements i, j or both. The same functions work on any sequences whose elements support equality, not just characters. Pass lists of lines and the result is a list of lines in common (drop the final join), which is exactly what a diff needs.

Less memory: two rows and Hirschberg

The full table costs nm cells. For two 100,000-line files that is ten billion integers, which is not going to fit. If you only need the length, notice that row i reads only row i-1, so two rows suffice. Iterate over the longer sequence and keep rows the size of the shorter one, giving O(min(n, m)) space.

def lcs_length(a, b):
    if len(b) > len(a):
        a, b = b, a                     # rows sized by the shorter input
    prev = [0] * (len(b) + 1)
    for x in a:
        cur = [0] * (len(b) + 1)
        for j, y in enumerate(b, 1):
            cur[j] = prev[j - 1] + 1 if x == y else max(prev[j], cur[j - 1])
        prev = cur
    return prev[-1]

Two rows lose the traceback, though. Hirschberg's algorithm (1975) recovers the subsequence in linear space with a divide-and-conquer trick. Split X in half. Compute the last row of LCS lengths for the top half of X against every prefix of Y, and the same for the reversed bottom half against every reversed suffix of Y. The split point k of Y that maximises left[k] + right[m-k] is where an optimal path crosses the middle row. Recurse on the two halves.

def _last_row(a, b):
    prev = [0] * (len(b) + 1)
    for x in a:
        cur = [0] * (len(b) + 1)
        for j, y in enumerate(b, 1):
            cur[j] = prev[j - 1] + 1 if x == y else max(prev[j], cur[j - 1])
        prev = cur
    return prev


def hirschberg(a, b):
    if not a or not b:
        return ""
    if len(a) == 1:
        return a if a[0] in b else ""
    mid = len(a) // 2
    left = _last_row(a[:mid], b)
    right = _last_row(a[mid:][::-1], b[::-1])
    k = max(range(len(b) + 1), key=lambda s: left[s] + right[len(b) - s])
    return hirschberg(a[:mid], b[:k]) + hirschberg(a[mid:], b[k:])

The work at each level of recursion sums to at most about nm, and the levels shrink geometrically, so total time stays O(nm), roughly double the plain table. Space is O(n + m). On the example it returns BDAB rather than BCBA. That is a different, equally long answer, which is a useful reminder that tests should check length and validity rather than one specific string.

From LCS to diff

Treat each file as a sequence of lines. Lines in the LCS are unchanged. Lines of the old file not in it are deletions, and lines of the new file not in it are insertions. With only insertions and deletions allowed, the minimum number of edits is n + m - 2 * LCS, so a shortest diff and a longest common subsequence are two views of the same answer.

def line_diff(old, new):
    L = lcs_table(old, new)
    i, j, ops = len(old), len(new), []
    while i > 0 or j > 0:
        if i > 0 and j > 0 and old[i - 1] == new[j - 1]:
            ops.append("  " + old[i - 1]); i, j = i - 1, j - 1
        elif j > 0 and (i == 0 or L[i][j - 1] >= L[i - 1][j]):
            ops.append("+ " + new[j - 1]); j -= 1
        else:
            ops.append("- " + old[i - 1]); i -= 1
    return ops[::-1]

Production diff tools do not fill an nm table. Myers' 1986 algorithm runs in O((n + m) D) time, where D is the size of the edit script, so near-identical files diff almost linearly. Implementations also strip common prefixes and suffixes first, and Git's patience diff anchors on lines that occur once in each file (histogram diff extends that to rare lines) to produce more readable hunks. Those refinements change which of several optimal or near-optimal alignments you see, not the underlying idea.

Variants worth knowing

• Longest common substring. Contiguous version: on a match S[i][j] = S[i-1][j-1] + 1, otherwise 0, and the answer is the maximum cell anywhere, not the corner.
• Longest palindromic subsequence. The LCS of a string and its reverse. See longest palindromic subsequence.
• Edit distance. Same table shape with a substitution option and costs instead of matches. Levenshtein distance is the classic version.
• Permutations. If both inputs are permutations of the same distinct symbols, map each symbol to its position in Y. The LCS becomes the longest increasing subsequence of the mapped X, solvable in O(n log n).
• Many sequences. For k sequences the table has k dimensions. LCS of an arbitrary number of sequences is NP-hard (Maier, 1978), so multiple alignment tools rely on heuristics.
• Weighted matches. Score matches by importance, for example rare tokens higher. The table shape is the same as in knapsack-style DPs: replace +1 with a weight.

Failure modes

• Off-by-one indexing. The table is (n+1) by (m+1) and compares a[i-1] with b[j-1]. Indexing the strings with i and j reads past the end or skips the first character.
• Substring recurrence by mistake. Resetting to 0 on a mismatch computes longest common substring. Check a pair like AXB and AB, whose LCS is 2 and longest common substring is 1.
• Asserting one specific answer. Ties make several outputs correct. Test that the result is a subsequence of both inputs and has the optimal length, preferably against a brute force on small random inputs.
• Quadratic memory on big inputs. nm cells of 4 to 8 bytes run out of memory long before nm operations run out of time. Use two rows for the length and Hirschberg, or a Myers-style diff, for alignments.
• Recursive memoisation in Python. A top-down version recurses up to n + m deep and hits the default recursion limit around a thousand. Prefer the bottom-up loop.
• Unnormalised input for diffs. Trailing whitespace and line-ending differences make every line unique and the LCS collapses. Normalise before comparing, or report the whitespace difference separately.

What to do next

1. Implement lcs_table and lcs_string from memory, then test them against a brute-force subsequence search on random strings over a small alphabet.
2. Fill the ABCBDAB / BDCABA table by hand once, and check it against the figure.
3. Change the tie-break and confirm you get a different subsequence with the same length.
4. Convert your solver to the two-row length version, then implement Hirschberg and compare its output length to the full-table version.
5. Write line_diff over two versions of a source file and compare its output with git diff; note where the hunks differ and why.
6. Solve longest common substring and longest palindromic subsequence as variations on the same table to make the recurrence choices stick.