Longest common subsequence

The longest common subsequence (LCS) of two strings is the longest string you can get from both of them by crossing out letters. Dynamic programming finds it with a grid over the beginnings (prefixes) of the two strings in $O(nm)$ time. If the common part must be one unbroken piece of both strings, you need the longest common substring, a different problem.

Explanation

Two strings, ABCD and ACBD, are written one above the other. Cross out as few letters as possible from both strings so that what is left is the same in both.

Let's join each kept letter of the top string to its partner in the bottom string by a line. The lines join equal letters, a letter has at most one line, and two lines can't cross, because the kept letters must come in the same order in both strings. The A's and the D's can be joined freely, but the B's and the C's cross, so only one of those two lines can stay, and the answer is 3 letters: ABD or ACD.

Fig. 1. The top row is the first string, the bottom row the second. A line joins two kept letters. The red line crosses a blue one.

What is left after crossing out letters, with the order kept, is a subsequence of a string. A string that is a subsequence of both strings is a common subsequence.

Repeated letters make it harder, as in the strings BCABD and ABCBDA. B appears twice in each, and A twice in the second, so a letter can have several partners, and it is not clear which one to choose.

The natural thing to try is to go through the first string from the left and join each letter to the first equal letter of the second string to the right of the last line. B goes to the B in place 2, C to the C in place 3, and A to the A at the very end of the second string. Now nothing is left to the right of that line, so the B and the D of the first string get no partner. That is 3 letters, the red lines in Fig. 2. The blue lines keep 4 letters, BCBD, by crossing out the A of the first string, so joining the A cost two letters later. Starting from the second string doesn't help either: it keeps ABD, again 3 letters.

Fig. 2. Left: joining each letter to the first partner that fits keeps 3 letters (red). Right: the blue lines keep 4.

Whether a line is a good choice depends on the letters after it. The last letters have nothing after them, so there the choice is simple. Here they are D in the first string and A in the second, and they are different. A line from the D must go to a letter before the A, and a line from the A must go to a letter before the D, so two such lines would cross. So at least one of them is crossed out. We don't know which, so we try both. Either we cross out the D and solve the question for BCAB and ABCBDA, or we cross out the A and solve it for BCABD and ABCBD. Each is the same question on shorter strings, and the answer is the better of the two.

When the last letters are equal, we join them to each other. That line is the last one in both strings, so it crosses nothing. What is left is the same question for the two strings without their last letters, and the answer is 1 more than its answer.

Every question that comes up this way is about the first few letters of one string against the first few letters of the other, a prefix of each. The first string has 6 prefixes (from 0 to 5 letters) and the second has 7, so there are only 6 × 7 = 42 such questions, and their answers fit in a grid.

This is dynamic programming over pairs of prefixes. The first string runs down the left of the grid and the second across the top. The row and the column marked ∅ are the empty prefix, which has nothing in common with anything, so they hold 0. Crossing out the last letter of the first string moves one row up, crossing out the last letter of the second string moves one column left, and crossing out both moves to the up-left cell. So the corner cell, for BCABD and ABCBDA, is the larger of the cell above it (BCAB and ABCBDA, 3) and the cell to its left (BCABD and ABCBD, 4), which is 4.

The player fills the whole grid, one cell per step.

The empty prefix has nothing in common with anything, so row ∅ and column ∅ hold 0.

B and A are different, so the cell takes the larger of the cell above and the cell to the left: 0.

B and B are equal, so the cell gets the up-left number plus 1: 0 + 1 = 1.

B and C are different: the larger of 0 and 1 is 1.

B and B are equal: 0 + 1 = 1.

B and D are different: the larger of 0 and 1 is 1.

B and A are different: the larger of 0 and 1 is 1.

C and A are different: the larger of 0 and 0 is 0.

C and B are different: the larger of 1 and 0 is 1.

C and C are equal: 1 + 1 = 2.

C and B are different: the larger of 1 and 2 is 2.

C and D are different: the larger of 1 and 2 is 2.

C and A are different: the larger of 1 and 2 is 2.

A and A are equal: 0 + 1 = 1.

A and B are different: the larger of 1 and 1 is 1.

A and C are different: the larger of 2 and 1 is 2.

A and B are different: the larger of 2 and 2 is 2.

A and D are different: the larger of 2 and 2 is 2.

A and A are equal: 2 + 1 = 3.

B and A are different: the larger of 1 and 0 is 1.

B and B are equal: 1 + 1 = 2.

B and C are different: the larger of 2 and 2 is 2.

B and B are equal: 2 + 1 = 3.

B and D are different: the larger of 2 and 3 is 3.

B and A are different: the larger of 3 and 3 is 3.

D and A are different: the larger of 1 and 0 is 1.

D and B are different: the larger of 2 and 1 is 2.

D and C are different: the larger of 2 and 2 is 2.

D and B are different: the larger of 3 and 2 is 3.

D and D are equal: 3 + 1 = 4.

D and A are different: the larger of 3 and 4 is 4. The corner holds the answer, 4.

1/31

Fig. 3. The grid filled row by row. Each cell holds the length of the LCS of two prefixes.

Step through the grid with the buttons or the arrow keys. The blue arrows point from the cells the new cell reads, and the framed one is the larger.

The grid holds the answer for every pair of prefixes, not only for the full strings: for example, BCA and ABCBDA have 3 letters in common.

Walking back to the letters

The corner says 4 but not which letters. To get the LCS itself, walk back from the corner. Where the letters of the row and the column are equal, the letter is in the LCS, and the walk goes up and to the left. Where they are different, the walk goes to the neighbour, above or to the left, that holds the same number as the current cell, since that is the cell the number came from. The walk stops at row ∅ or column ∅. Fig. 4 shows it on the example. The letters come out from the end, so the LCS read forwards is BCBD.

Fig. 4. The walk back from the corner. A thick arrow takes a letter. The letters are collected from the end, so they are read backwards.

Sometimes both neighbours hold the same number as the current cell. Then both moves are right: the strings have several longest common subsequences, and the choice decides which one the walk finds.

Shortest common supersequence

The shortest string that contains both strings as subsequences is called their shortest common supersequence. It has n + m - L letters, where $n$ and $m$ are the lengths of the strings and $L$ is the length of their LCS. Each letter of an LCS is written once and serves both strings, and every other letter is written once. No shorter string works: in any such string, the letters used by both strings form a common subsequence, so at most $L$ letters are shared. For BCABD and ABCBDA that is 5 + 6 − 4 = 7 letters, for example ABCABDA.

Fewest deletions to a palindrome

The fewest letters to cross out of a string to leave a palindrome is $n$ minus the length of the LCS of the string and its reverse. A palindrome reads the same both ways, so what is left is a common subsequence of the string and its reverse. For BCABD and its reverse DBACB the LCS has 3 letters, BCB, so 2 letters are crossed out.

The other direction also holds, though an LCS of a string and its reverse need not be a palindrome itself. Each line of such an LCS joins two places of the string with the same letter: a place $p$, counted from the left, and a place $q$, counted from the left in the string but found in the reverse. Going along the lines, $p$ grows and $q$ shrinks. Take the $t$ lines where $p$ is still before $q$. Their letters at the places $p$, followed by their letters at the places $q$, form a palindrome of $2t$ letters. The other L - t lines give a palindrome of 2(L - t) letters the same way, or one letter fewer if p = q on the first of them. One of the two has at least $L$ letters, so the count is exact.

Implementation

The first string is a with n letters, the second is b with m letters. The grid is a vector of vectors f with n + 1 rows and m + 1 columns, and the five steps of dynamic programming are:

1. State: f[i][j] is the length of the LCS of the first i letters of a and the first j letters of b. For example, f[3][3] = 2 is the cell for BCA and ABC.

2. Base case: f[0][j] = 0 and f[i][0] = 0, because an empty prefix has nothing in common with anything.

3. Order: i from 1 to n, and inside it j from 1 to m, so each cell comes after the cells above, to the left and up-left.

4. Transition: if a[i - 1] == b[j - 1], then f[i][j] = f[i - 1][j - 1] + 1; otherwise f[i][j] = max(f[i - 1][j], f[i][j - 1]). Letter i of a is a[i - 1], because a string counts places from 0.

5. Answer: f[n][m].

After the answer, the program walks back from (n, m) as in Fig. 4. When the two neighbours hold the same number, it goes up.

The input is two lines: the first string and the second string, without spaces. The program prints the length of the LCS and then one LCS. For the main example the input is:

BCABD
ABCBDA

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
using namespace std;

int main() {
    string a, b;   // for arrays of numbers use vector<int> here and read the lengths first
    cin >> a >> b;
    int n = a.size(), m = b.size();

    // state: f[i][j] = LCS length of the first i letters of a and the first j letters of b
    // base case: row 0 and column 0 stay 0
    vector<vector<int>> f(n + 1, vector<int>(m + 1, 0));

    // order: row by row, each row left to right
    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= m; j++) {
            // transition: letter i of a is a[i - 1], letter j of b is b[j - 1]
            if (a[i - 1] == b[j - 1]) {
                f[i][j] = f[i - 1][j - 1] + 1;
            } else {
                f[i][j] = max(f[i - 1][j], f[i][j - 1]);
            }
        }
    }

    // answer
    cout << f[n][m] << "\n";

    // walk back from (n, m); remove this block if you only need the length
    string s;
    int i = n, j = m;
    while (i > 0 && j > 0) {
        if (a[i - 1] == b[j - 1]) {
            s += a[i - 1];   // this letter is in the LCS
            i--;
            j--;
        } else if (f[i - 1][j] >= f[i][j - 1]) {
            i--;             // up, also when the two are equal
        } else {
            j--;             // left
        }
    }
    reverse(s.begin(), s.end());   // the letters were collected from the end
    cout << s << "\n";
}

On the main example the program prints 4 and BCBD. On the strings ABCD and ACBD it prints ABD: at the cell for ABC and ACB both neighbours hold 2, and the walk goes up. A walk that goes left there gives ACD.

The grid has (n + 1)(m + 1) cells, and each is filled in constant time, so the program takes $O(nm)$ time and memory (see time complexity). The walk back takes O(n + m) steps more. For two strings of 5000 letters that is 2.5 \cdot 10^7 cells, well within a second by the rule of thumb of about 10^8 simple steps per second.

Only the length, with two rows

While it fills row i, the program reads only row i - 1 (f[i - 1][j - 1] and f[i - 1][j]) and row i itself (f[i][j - 1]). So when only the length is needed, prev can hold row i - 1 and cur row i. After each row, cur becomes prev. The memory drops to $O(m)$, and the time stays $O(nm)$. Put the shorter string in b to keep the rows short.

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
using namespace std;

int main() {
    string a, b;
    cin >> a >> b;
    int n = a.size(), m = b.size();

    vector<int> prev(m + 1, 0), cur(m + 1, 0);   // rows i - 1 and i
    for (int i = 1; i <= n; i++) {
        cur[0] = 0;
        for (int j = 1; j <= m; j++) {
            if (a[i - 1] == b[j - 1]) {
                cur[j] = prev[j - 1] + 1;
            } else {
                cur[j] = max(prev[j], cur[j - 1]);
            }
        }
        swap(prev, cur);   // row i becomes the row above
    }
    cout << prev[m] << "\n";
}

After the last swap, the last row is in prev, so the answer is prev[m]. On the main example the program prints 4.

Correctness proof

Take a cell (i, j) and assume every cell above it and to its left holds the right length. Here a_i is letter $i$ of the first string (a[i - 1] in the code), b_j is letter $j$ of the second, and f_{i,j} is f[i][j], so a_i and b_j are the last letters of the two prefixes. Every value the transition writes is the length of a real common subsequence, so it is never too big. We show it is never too small: take a longest common subsequence of the two prefixes, drawn as lines.

If a_i \ne b_j, the two last letters can't both have lines. Say a_i is joined to some letter of $b$ before b_j, and b_j to some letter of $a$ before a_i. Then the first line goes from the end of $a$ to earlier in $b$, the second from earlier in $a$ to the end of $b$, and they cross (Fig. 5). So one of the last letters has no line, and the whole subsequence fits in the prefixes without it: its length is at most f_{i-1,j} or f_{i,j-1}.

Fig. 5. The last letters D and A are different. If each had a line to another letter, the lines would cross.

If a_i = b_j, we can change the subsequence so that these two letters are joined to each other, without making it shorter. If neither has a line, add the line between them: it is the last in both strings, so it crosses nothing, and the subsequence gets longer, which is impossible for a longest one. If only a_i has a line, to an earlier equal letter of $b$, move its end to b_j: the letter is the same, and nothing is to the right of it. The same works if only b_j has a line. Both can't have lines to other letters, because those lines would cross, as above. Once the last letters are joined, the rest is a common subsequence of the two shorter prefixes, so the length is at most f_{i-1,j-1} + 1.

Common mistakes

• Joining each letter to the first equal letter that fits (Fig. 2). It keeps 3 letters where 4 are possible. Fill the grid instead.

• Adding 1 to the larger neighbour when the letters are equal, f[i][j] = max(f[i - 1][j], f[i][j - 1]) + 1, which gives one letter two lines. For B and BB it prints 2 instead of 1. Equal letters read the up-left cell.

• Keeping one number per prefix of the first string: f[i] is the LCS of the first i letters of a with all of b, and f[i] = f[i - 1] + 1 when letter i occurs in b. This state doesn't remember where in b the last line ended, so lines can cross. For AB and BA it gives 2, but the answer is 1. The second index records how much of b is still free.

• Using LCS when the problem asks for a common substring. For AXBYC and ABC the LCS has 3 letters, but the longest common substring has 1. For a substring, a cell resets to 0 when the letters differ, and the answer is the largest cell, not the corner. For long strings, use a suffix array or hashing.

Similar algorithms

Longest increasing subsequence keeps an increasing subsequence of one sequence. When both sequences are permutations of the numbers 1 to $n$, rename each number of the first sequence by its place in it and apply the same renaming to the second. Their LCS is then the longest increasing subsequence of the renamed second sequence, found in O(n \log n).

Edit distance uses the same grid over prefixes, but counts the fewest inserted, deleted and replaced letters that turn one string into the other. With only inserting and deleting, it equals n + m - 2L: every letter outside the LCS is deleted from the first string or inserted from the second, once.

If letters may only be added at the end of a string to make a palindrome, it is not an LCS problem: the question is the longest palindrome at the end of the string, which the Z-function finds.

Practice