The longest common subsequence

A subsequence of a string is what is left after deleting any characters you like, without changing the order of the rest: ACE is a subsequence of ABCDE, but EA is not. Unlike a substring, a subsequence doesn't have to be contiguous. The longest common subsequence (LCS) of two strings is the longest string that is a subsequence of both. For ABCB and BDCAB it is BCB, of length 3.

S1 = ABCB above S2 = BDCAB; lines join the B, C and B of S1 to the B, C and B of S2, in the same order with gaps allowed, giving BCB
A common subsequence keeps letters in the same order in both strings but may skip letters; BCB is the longest here.

Checking every subsequence of one string against the other is hopeless: a string of length m has 2m subsequences. Dynamic programming solves the problem in time proportional to the product of the two lengths. Type two strings (up to 15 characters each) into the S1 and S2 fields and press LCS Recursive, LCS Memoized or LCS Table.

The recurrence

The page's code works on prefixes, identified by the index of their last character: LCS(S1, S2, x, y) is the LCS length of S1[0..x] and S2[0..y], and an index of -1 means the empty prefix. Compare the last characters of the two prefixes:

  • If one prefix is empty, the LCS is empty: length 0.
  • If S1[x] == S2[y], that shared last character can end a longest common subsequence, so the answer is 1 plus the LCS of both prefixes with it removed.
  • Otherwise at least one of the two last characters is not in the LCS. Drop the last character of S1, or the last of S2, solve both and take the better one.
def LCS(S1, S2, x, y):
    if (x == -1) or (y == -1):
        return 0
    else if (S1[x] == S2[y]):
        return 1 + LCS(S1, S2, x-1, y-1)
    else:
        return max(LCS(S1, S2, x-1, y), LCS(S1, S2, x, y-1))

The answer for the whole strings is LCS(S1, S2, len(S1)-1, len(S2)-1).

Version 1: plain recursion (LCS Recursive)

Running the recurrence as it stands works, but in the mismatch case it makes two calls that each shrink the problem by only one character, and those two calls share most of their subproblems: LCS(x-1, y) and LCS(x, y-1) can both go on to call LCS(x-1, y-1). These are overlapping subproblems. In the worst case, when the strings have no character in common, two strings of length n cause 2·C(2n, n) − 1 calls: 25,739 calls for two 8-letter strings and 369,511 for two 10-letter strings, while there are at most (n+1)2 different pairs (x, y) (81 and 121). That grows only a little more slowly than 4n. On the canvas every call is a line LCS(S1, S2, x, y) [LCS(prefix1, prefix2)], indented by its depth, and is replaced by its value when it returns.

Version 2: memoization (LCS Memoized)

Memoization keeps the recursion but stores every result in a table T indexed by (x+1, y+1), so that the -1 indices fit. Each call first looks in the table and only recurses if the entry is still empty:

def LCSMem(S1, S2, x, y):             # T[i][j] == -1 means "not computed yet"
    if T[x+1][y+1] != -1:
        return T[x+1][y+1]            # already known: just read it
    if (x == -1) or (y == -1):
        result = 0
    else if (S1[x] == S2[y]):
        result = 1 + LCSMem(S1, S2, x-1, y-1)
    else:
        result = max(LCSMem(S1, S2, x-1, y), LCSMem(S1, S2, x, y-1))
    T[x+1][y+1] = result
    return result

Each pair (x, y) is now computed at most once. On the canvas you see each result fly into its table cell, and on a later request fly back out of it. Only the cells the recursion actually needs get filled, so some cells may stay empty.

Version 3: bottom-up table (LCS Table)

A cell only depends on the cell diagonally before it, the cell to its left and the cell above, so the whole table can be filled in order without recursion. On this page S1 runs across the top and S2 down the left side, with the blue indices -1, 0, 1, ... next to them. The first row and column (index -1, the empty prefix) are all 0; then the table is filled row by row:

for i in 0 .. len(S1):  T[i][0] = 0          # empty prefix of S2
for j in 0 .. len(S2):  T[0][j] = 0          # empty prefix of S1
for j in 1 .. len(S2):                       # one row per character of S2
    for i in 1 .. len(S1):                   # one column per character of S1
        if S1[i-1] == S2[j-1]:
            T[i][j] = T[i-1][j-1] + 1        # diagonal + 1
        else:
            T[i][j] = max(T[i-1][j], T[i][j-1])   # left, above

While a cell is being filled, its two characters are highlighted; on a match the diagonal value plus one moves into the cell, on a mismatch the left and upper neighbours are highlighted and the larger value moves in. The length of the LCS ends up in the bottom-right cell.

Recovering the subsequence. The table stores lengths, not the string itself. To read off an actual LCS, start in the bottom-right cell and walk back. If the two characters of the current cell match, that character belongs to the LCS: record it and move diagonally up-left. Otherwise move to the neighbour that holds the larger value (left or up; on a tie this page moves up). Stop at the border. The characters are found from last to first, and the canvas collects them under Sequence. Both the Table and the Memoized buttons end with this walk.

Worked example: S1 = ABCB, S2 = BDCAB

                   S1:     A   B   C   B
          index      -1    0   1   2   3
  S2   -1             0    0   0   0   0
  B     0             0    0   1   1   1
  D     1             0    0   1   1   1
  C     2             0    0   1   2   2
  A     3             0    1   1   2   2
  B     4             0    1   2   2   3

Some of the cells in detail:

  • Row B (S2[0]), column B (S1[1]): the characters match, so the cell is the diagonal 0 plus 1 = 1. Every cell to its right and below is at least 1.
  • Row D, column C: D and C differ, so the cell is max(left 1, above 1) = 1.
  • Row C, column C: a match; the diagonal (row D, column B) holds 1, so the cell is 2 — the common subsequence BC.
  • Row A, column A: a match with diagonal 0, so 1 (the subsequence A).
  • Row B (S2[4]), column B (S1[3]): a match; the diagonal (row A, column C) holds 2, so the answer is 3.

Traceback, written as recursion indices (S1 index, S2 index):

(3,4)  S1[3]=B, S2[4]=B  match        -> take B, go diagonal to (2,3)
(2,3)  S1[2]=C, S2[3]=A  mismatch     left = 1, up = 2  -> go up to (2,2)
(2,2)  S1[2]=C, S2[2]=C  match        -> take C, go diagonal to (1,1)
(1,1)  S1[1]=B, S2[1]=D  mismatch     left = 0, up = 1  -> go up to (1,0)
(1,0)  S1[1]=B, S2[0]=B  match        -> take B, go diagonal to the border: stop
LCS table with S1 = ABCB across and S2 = BDCAB down, bottom-right cell 3. The walk back goes from that 3 diagonally (match B) to row A column C, up to row C column C, diagonally (match C) to row D column B, up to row B column B, diagonally (match B) to the border. A side panel shows the cell rule: on a match the diagonal plus 1, otherwise the max of left and up
Each cell is diagonal + 1 on a match or max(left, up) otherwise; walking back from the corner and taking the matched letters gives BCB.

Read backwards, the characters taken are B, C, B: the LCS is BCB, and here it is the only common subsequence of length 3. For this small input the plain recursion makes 30 calls; memoization makes 19 (15 different subproblems plus 4 lookups). The classic textbook pair ABCBDAB / BDCABA gives length 4; with this page's tie rule the traceback finds BDAB.

Why it is correct

Let Z be a longest common subsequence of the prefixes ending at x and y. If S1[x] == S2[y], we may assume Z ends with that character (if it didn't, appending that character to Z would give a longer common subsequence), and the rest of Z is a common subsequence of the shorter prefixes, which must itself be longest. If the characters differ, Z can't end with both, so it is a common subsequence of one of the two shortened pairs, and the larger of the two answers is exactly right. This optimal substructure plus induction on x + y proves the recurrence; the three versions just evaluate it in different orders.

Running time and space

  • Plain recursion: exponential, up to 2·C(2n, n) − 1 calls for two strings of length n. The recursion depth is at most m + n.
  • Memoized and table: (m+1)(n+1) cells, each computed in constant time: Θ(mn) time and space. The traceback adds only O(m + n) steps.
  • If you only need the length, keep just two rows: O(min(m, n)) space. Hirschberg's algorithm recovers the subsequence itself in linear space too.

Common mistakes and variants

  • Subsequence vs. substring: the longest common substring (contiguous) uses a similar table, but a mismatch resets the cell to 0 instead of copying the larger neighbour.
  • Off-by-one: the table has one extra row and column for the empty prefix, so cell T[i][j] belongs to the characters S1[i-1] and S2[j-1]. Mixing up the two index systems is the most frequent bug.
  • The LCS is not unique: different tie rules in the traceback give different, equally long answers.
  • Recovering the sequence needs the full table (or stored "arrows"); the two-row trick only gives the length.
  • The comparison is case-sensitive: a and A are different characters.
  • Edit distance (Levenshtein) and sequence alignment (Needleman-Wunsch) are close relatives: same kind of table, but with costs for insertions, deletions and substitutions.

Where it is used

File comparison tools such as Unix diff and the diff views of version control systems find a longest common subsequence of the lines of two files: lines in the LCS are unchanged, the rest are shown as deleted or added. In bioinformatics, LCS and its weighted cousins align DNA and protein sequences to measure how related they are. The same idea is used in plagiarism detection, merging edits, and spell checkers.