Dynamic Programming/Sequence Alignment & LCS

Lesson 8.24,400 words

Sequence Alignment & LCS

Two strings can be compared by how much of one appears inside the other. The longest common subsequence (LCS) and edit distance are the two classic measures, and they are the same dynamic program with different costs.

╌╌╌╌

How similar are two strings? algorithm and altruistic share the letters a, l, t, i, c in order, and that shared string, the longest common subsequence, is one of the most useful measures of similarity in computing. It underlies the diff utility, version-control merges, and (with a change of cost function) the alignment of DNA and protein sequences in computational biology.1 This lesson develops the LCS dynamic program in full, then shows that edit distance is the same dynamic program with the costs rearranged.

The problem

A subsequence of a string is what remains after deleting zero or more characters, keeping the rest in their original order. It need not be contiguous: ace is a subsequence of abcde, but aec is not. Given two strings and , a common subsequence is a string that is a subsequence of both, and we want a longest one.

A common subsequence is a set of order-preserving matches threading the two strings: the matched letters appear in both, left to right, though the gaps between them differ.

A common subsequence threads order-preserving matches through both strings. Here is the longest common subsequence of and ; matched letters connect, the rest are skipped.

A brute-force search is hopeless: has subsequences, and checking each against gives . We need the recipe, and it pays to follow it literally. The DP recipe runs through fixed steps: (0) simplify the goal, (1) define the subproblems and notation, (2) write the DP equations, (3) prove them correct by induction, (4) turn them into iterative pseudocode, (5) return to the original goal, then analyze. We walk LCS through those steps verbatim.

Step 0–1: Simplify the goal, define the subproblem

Simplify first. Computing the string drags around bookkeeping; computing its length is cleaner. So we first solve for the LCS length, then recover an actual subsequence in a cheap second pass (Step 5). With that simplification, the decisive move, which recurs across all sequence DPs, is to index subproblems by prefixes of the two strings. Write and for the two inputs.

The answer we want is . There are subproblems, one per pair of prefix lengths and .

Step 2: The DP equations

Now apply optimal substructure by looking at the last characters, and .2 Write the recurrence as a single over three cases, plus a base case:

Read the three branches as moves on the prefixes:

  • Case 1 drops : an LCS of and is a common subsequence of and , so .
  • Case 2 drops symmetrically, so .
  • Case 3 fires only when : the shared character can end an optimal LCS, so we append it to the best LCS of the strictly shorter prefixes, giving .

Why exactly these three? Because every common subsequence of and falls into one of three shapes: it ignores (case 1), ignores (case 2), or uses both as a final match (case 3, possible only if ). The over the applicable cases is the longest of them. When , only cases 1 and 2 apply, and the recurrence collapses to .

Each entry depends only on its left, upper, and upper-left neighbors, so filling the table row by row, left to right (or column by column) respects every dependency. The same three-neighbour stencil drives every sequence DP below: only the labels on the arrows change.

The shared stencil of every sequence DP: cell reads its left, upper, and upper-left neighbours; the diagonal fires on a match (), the other two on a drop or an edit.

Step 3: Correctness by induction on

Step 4: Iterative pseudocode

The equations are non-circular, since every entry reads strictly smaller prefixes, so they convert directly into a bottom-up table fill. Allocate , zero the border, then sweep:

Algorithm 1:LCS-Length(A[1..m],B[1..n])\textsc{LCS-Length}(A[1..m], B[1..n]) — fill the DP table
  1. 1
    for i0i \gets 0 to mm do
  2. 2
    OPT[i][0]0\OPT[i][0] \gets 0
    empty BB prefix
  3. 3
    for j0j \gets 0 to nn do
  4. 4
    OPT[0][j]0\OPT[0][j] \gets 0
    empty AA prefix
  5. 5
    for i1i \gets 1 to mm do
  6. 6
    for j1j \gets 1 to nn do
  7. 7
    OPT[i][j]max(OPT[i1][j], OPT[i][j1])\OPT[i][j] \gets \max\parens{\OPT[i-1][j],\ \OPT[i][j-1]}
    cases 1, 2
  8. 8
    if A[i]=B[j]A[i] = B[j] then
  9. 9
    OPT[i][j]max(OPT[i][j], OPT[i1][j1]+1)\OPT[i][j] \gets \max\parens{\OPT[i][j],\ \OPT[i-1][j-1] + 1}
    case 3
  10. 10
    return OPT[m][n]\OPT[m][n]

Every cell costs , so the fill is .

The DP table, filled

Take and . We build the table of values. Row and column are all zero (empty prefix); every other cell is filled by the recurrence. The shaded diagonal steps mark the matches that build the answer, and the red arrows trace the reconstruction walk (Step 5) backwards from the corner.

Filled LCS table for and with the traceback path arrowed.

The bottom-right entry reads : the longest common subsequence of BDCAB and ABCB has length , namely BCB (the only one here, though in general there may be ties). The arrows enter each shaded match cell diagonally (emitting B, then C, then B, read from the corner upward) and step straight up or left through non-match cells, exactly as the reconstruction below prescribes.

Step 5: Reconstructing the subsequence

The table gives the length; this is where the Step 0 simplification is paid back. In a second pass we walk backwards from , undoing the recurrence. At cell : if , that character belongs to the LCS — emit it and step diagonally to (case 3); otherwise move to whichever neighbor, up or left, holds the larger value (the one the of cases 1 and 2 chose).

Algorithm 2:LCS-Reconstruct(A,B,OPT,i,j)\textsc{LCS-Reconstruct}(A, B, \OPT, i, j) — recover the subsequence
  1. 1
    if i=0i = 0 or j=0j = 0 then
  2. 2
    return the empty string
    empty prefix
  3. 3
    if A[i]=B[j]A[i] = B[j] then
  4. 4
    return LCS-Reconstruct(A,B,OPT,i1,j1)\textsc{LCS-Reconstruct}(A, B, \OPT, i-1, j-1) followed by A[i]A[i]
    case 3 match
  5. 5
    else if OPT[i1][j]OPT[i][j1]\OPT[i-1][j] \ge \OPT[i][j-1] then
  6. 6
    return LCS-Reconstruct(A,B,OPT,i1,j)\textsc{LCS-Reconstruct}(A, B, \OPT, i-1, j)
    from above (case 1)
  7. 7
    else
  8. 8
    return LCS-Reconstruct(A,B,OPT,i,j1)\textsc{LCS-Reconstruct}(A, B, \OPT, i, j-1)
    from left (case 2)

The walk takes one step toward the origin each call, so it runs in time, cheap compared with building the table.

Trace it on the worked table above, starting at with value . The red arrows in the figure are this walk:

  • : , a match. Emit B, step to .
  • : . Compare the up neighbour against the left neighbour ; the tie breaks upward, step to .
  • : , a match. Emit C, step to .
  • : . Up neighbour ties left neighbour ; step up to .
  • : , a match. Emit B, step to .
  • : , stop.

The emitted characters, corner-first, are B, C, B; reversed into forward order they read BCB — the length- LCS the corner promised. The tie-break rule (up before left) is arbitrary; the other choice would recover an equally long subsequence, and a longest common subsequence need not be unique.

Running time and space

The table has entries. Filling one costs a character comparison and a of at most three previously-computed neighbours, all . Summing over the fill,

and the border initialization adds only , which absorbs. The reconstruction pass is , also absorbed. LCS therefore runs in time — a decisive improvement over the brute force, and for two length- strings the difference is a million cell updates against roughly subsequence checks.

Space is for the full table, but the recurrence reads only the current row and the one directly above it. Keeping two length- rows and swapping them after each (a rolling array) drops the footprint to , and choosing the shorter string as the inner axis makes it .

Rolling-array space optimization. Cell reads only the previous row (upper and upper-left) and the current row's left neighbour, so two rows suffice; after finishing row the previous row is discarded and the buffers swap.

The catch is universal to this trick: collapsing the table erases the information the traceback needs, so the two-row version yields only the length, never the subsequence itself. Recovering the actual alignment in linear space is possible — Hirschberg's divide-and-conquer refinement does it in time and space3 — but that machinery is a topic for later.

longest_common_subsequence.pypython
from typing import Sequence, TypeVar

Element = TypeVar("Element")

def lcs_length(left: Sequence[Element], right: Sequence[Element]) -> int:
  """
    Length of a longest common subsequence of `left` and `right`.\n
    Fills the full (m+1) x (n+1) prefix table and reports the corner.\n
  """
  rows: int = len(left)
  cols: int = len(right)

  # table[i][j] = LCS length of left[:i] and right[:j]; border rows/cols are 0
  # because an empty prefix shares nothing.
  table: list[list[int]] = [[0 for _ in range(cols + 1)] for _ in range(rows + 1)]

  for i in range(1, rows + 1):
    for j in range(1, cols + 1):
      if left[i - 1] == right[j - 1]:
        # case 3: extend the diagonal by the shared final element.
        table[i][j] = table[i - 1][j - 1] + 1
      else:
        # cases 1, 2: drop one side, keep the longer survivor.
        table[i][j] = max(table[i - 1][j], table[i][j - 1])

  return table[rows][cols]

def lcs(left: Sequence[Element], right: Sequence[Element]) -> list[Element]:
  """
    A longest common subsequence itself, as a list of elements.\n
    Builds the same table, then walks backwards from the bottom-right corner:\n
    a match emits its element and steps diagonally; otherwise it follows the\n
    neighbour the max chose. Ties are broken toward the upper neighbour, so the\n
    result is one valid LCS (there may be several of equal length).\n
  """
  rows: int = len(left)
  cols: int = len(right)

  # empty-prefix borders are 0; interior cells are filled below.
  table: list[list[int]] = [[0 for _ in range(cols + 1)] for _ in range(rows + 1)]

  # fill the prefix table: diagonal +1 on a match, else the longer drop.
  for i in range(1, rows + 1):
    for j in range(1, cols + 1):
      if left[i - 1] == right[j - 1]:
        table[i][j] = table[i - 1][j - 1] + 1
      else:
        table[i][j] = max(table[i - 1][j], table[i][j - 1])

  # walk back from the corner, emitting matches and following the max choice.
  result: list[Element] = []
  row, col = rows, cols
  while row > 0 and col > 0:
    if left[row - 1] == right[col - 1]:
      result.append(left[row - 1])
      row -= 1
      col -= 1
    elif table[row - 1][col] >= table[row][col - 1]:
      row -= 1
    else:
      col -= 1

  # the walk emits corner-to-origin, so flip back to forward order.
  result.reverse()
  return result
longest_palindromic_subsequence.pypython
from typing import Sequence, TypeVar

Element = TypeVar("Element")

def lps_length(sequence: Sequence[Element]) -> int:
  """
    Length of a longest palindromic subsequence of `sequence`.\n
    Fills the substring table over spans of increasing length.\n
  """
  count: int = len(sequence)
  if count == 0:
    return 0

  # table[i][j] = LPS length of sequence[i..j]; single elements are length 1.
  table: list[list[int]] = [[0 for _ in range(count)] for _ in range(count)]
  for index in range(count):
    table[index][index] = 1

  for span in range(2, count + 1):
    for low in range(count - span + 1):
      high: int = low + span - 1
      if sequence[low] == sequence[high]:
        # endpoints pair off; span 2 has no interior, so guard the inner read.
        inner: int = table[low + 1][high - 1] if span > 2 else 0
        table[low][high] = inner + 2
      else:
        table[low][high] = max(table[low + 1][high], table[low][high - 1])

  return table[0][count - 1]

def lps(sequence: Sequence[Element]) -> list[Element]:
  """
    A longest palindromic subsequence itself, as a list of elements.\n
    Builds the same table, then walks inward from the (0, n-1) corner: matching\n
    endpoints are emitted on both sides; otherwise it follows the longer\n
    neighbour. The two halves are stitched around the (possibly empty) centre.\n
  """
  count: int = len(sequence)
  if count == 0:
    return []

  table: list[list[int]] = [[0 for _ in range(count)] for _ in range(count)]
  for index in range(count):
    table[index][index] = 1

  for span in range(2, count + 1):
    for low in range(count - span + 1):
      high: int = low + span - 1
      if sequence[low] == sequence[high]:
        inner: int = table[low + 1][high - 1] if span > 2 else 0
        table[low][high] = inner + 2
      else:
        table[low][high] = max(table[low + 1][high], table[low][high - 1])

  # Collect the left half forward and the right half reversed, then join.
  left_half: list[Element] = []
  right_half: list[Element] = []
  low: int = 0
  high: int = count - 1
  while low < high:
    if sequence[low] == sequence[high]:
      left_half.append(sequence[low])
      right_half.append(sequence[high])
      low += 1
      high -= 1
    elif table[low + 1][high] >= table[low][high - 1]:
      low += 1
    else:
      high -= 1

  centre: list[Element] = [sequence[low]] if low == high else []
  right_half.reverse()
  return left_half + centre + right_half
longest_increasing_subsequence.pypython
from bisect import bisect_left
from typing import Sequence, TypeVar

from comparable import Comparable

Element = TypeVar("Element", bound=Comparable)

def lis_length_quadratic(sequence: Sequence[Element]) -> int:
  """
    Length of a longest strictly increasing subsequence, in O(n^2).\n
    `best_ending[i]` is the LIS length finishing at index `i`; each entry looks\n
    back at every smaller predecessor and extends the best of them.\n
  """
  count: int = len(sequence)
  if count == 0:
    return 0

  # best_ending[i] = LIS length finishing at i; extend the best smaller run.
  best_ending: list[int] = [1 for _ in range(count)]
  for index in range(1, count):
    for earlier in range(index):
      if sequence[earlier] < sequence[index]:
        best_ending[index] = max(best_ending[index], best_ending[earlier] + 1)

  return max(best_ending)

def lis_length(sequence: Sequence[Element]) -> int:
  """
    Length of a longest strictly increasing subsequence, in O(n log n).\n
    `tails[k]` holds the smallest tail seen for an increasing run of length\n
    `k + 1`; binary search places each element, and the final length of\n
    `tails` is the answer.\n
  """
  tails: list[Element] = []
  for value in sequence:
    position: int = bisect_left(tails, value)
    if position == len(tails):
      tails.append(value)
    else:
      # value can end a length-(position+1) run with a smaller tail.
      tails[position] = value
  return len(tails)

def lis(sequence: Sequence[Element]) -> list[Element]:
  """
    A longest strictly increasing subsequence itself, in O(n log n).\n
    Runs the patience method while recording, for each placed element, the\n
    index of its predecessor in the chain, then walks the links back from the\n
    final element of the longest run.\n
  """
  count: int = len(sequence)
  if count == 0:
    return []

  # tail_index[k] = index in `sequence` of the current smallest tail of a
  # length-(k+1) run. predecessor[i] = index chained before i, or -1.
  tail_values: list[Element] = []
  tail_index: list[int] = []
  predecessor: list[int] = [-1 for _ in range(count)]

  for index in range(count):
    value: Element = sequence[index]
    position: int = bisect_left(tail_values, value)

    # link each element to the tail of the run it extends.
    if position > 0:
      predecessor[index] = tail_index[position - 1]

    # start a new longest run, or lower the tail of an existing length.
    if position == len(tail_values):
      tail_values.append(value)
      tail_index.append(index)
    else:
      tail_values[position] = value
      tail_index[position] = index

  # walk the predecessor chain back from the tail of the longest run.
  result: list[Element] = []
  cursor: int = tail_index[-1]
  while cursor != -1:
    result.append(sequence[cursor])
    cursor = predecessor[cursor]

  result.reverse()
  return result
distinct_subsequences.pypython
from typing import Sequence, TypeVar

Element = TypeVar("Element")

def distinct_subsequences(
  source: Sequence[Element],
  target: Sequence[Element],
) -> int:
  """
    Number of distinct subsequences of `source` that equal `target`.\n
    Counts placements by sweeping the prefix table; positions of repeated\n
    target elements are summed, so identical subsequences are not double\n
    counted.\n
  """
  rows: int = len(source)
  cols: int = len(target)

  # count[i][j] = subsequences of source[:i] equal to target[:j]. The empty
  # target is matched once by every source prefix (the empty pick).
  count: list[list[int]] = [[0 for _ in range(cols + 1)] for _ in range(rows + 1)]
  for i in range(rows + 1):
    count[i][0] = 1

  for i in range(1, rows + 1):
    for j in range(1, cols + 1):
      # ways that skip source[i-1] entirely.
      count[i][j] = count[i - 1][j]
      if source[i - 1] == target[j - 1]:
        # plus ways that use source[i-1] to cover target[j-1].
        count[i][j] += count[i - 1][j - 1]

  return count[rows][cols]
comparable.pypython
from typing import Any, Protocol, TypeVar


class Comparable(Protocol):
  """
    Anything orderable with `<` (int, float, str, tuple, date, …).\n
  """

  # `other` is position-only so built-ins (int, str, …), whose dunder
  # operands are position-only, structurally satisfy the protocol.
  def __lt__(self, other: Any, /) -> bool: ...
  def __gt__(self, other: Any, /) -> bool: ...
  def __le__(self, other: Any, /) -> bool: ...
  def __ge__(self, other: Any, /) -> bool: ...

The same machine: edit distance

Edit distance (the Levenshtein distance) asks the closely related question: what is the minimum number of single-character insertions, deletions, and substitutions that transform into ?4 It is the cost model behind spell-checkers and diff, and it is structurally identical to LCS.

Again we look at the last characters. If , they need no edit and we align them for free. Otherwise we make one of three moves (delete , insert , or substitute ), each at cost , and recurse on the correspondingly shorter prefixes:

The base cases say it: turning a length- prefix into the empty string costs deletions, and building a length- prefix from nothing costs insertions.

It is the same three-neighbour stencil as LCS — only the arrow labels change. The diagonal is free on a match and costs to substitute; a step down deletes , a step right inserts , each at cost . We take the cheapest incoming move instead of the longest:

Edit distance's three moves into cell : the diagonal is free on a match and costs to substitute, a step down deletes , and a step right inserts — the LCS stencil, now minimizing edits instead of maximizing matches.

Filled on a small pair, the table looks just like the LCS one but now minimizes. Take and : the shaded diagonal marks the free matches, and the corner reports the answer.

Edit-distance table for ; the corner counts one substitution () and one insertion (). The red arrows trace the alignment back: a diagonal step is a match or substitution, a horizontal step an insertion.

Reading the corner, . The traceback recovers the alignment itself, walking corner to origin and reading each step as the edit that produced it:

  • value : the left neighbour is one cheaper, so this step is an insertion of . Move left to .
  • value : , and the upper-left is the cheapest source, so this is a substitutionT R. Move diagonally to .
  • value : , a free match. Move diagonally to .
  • value : , a free match. Move to and stop.

Read forward, the alignment is: keep C, keep A, substitute T R, insert S — turning CAT into CARS in the promised two edits. The match cells on the diagonal ( and ) copy the upper-left value unchanged, the same LCS diagonal step, but counting saved edits instead of matched characters.

The fill is the LCS loop with in place of and the borders seeded to the prefix lengths rather than zeros:

Algorithm 3:Edit-Distance(A[1..m],B[1..n])\textsc{Edit-Distance}(A[1..m], B[1..n]) — fill the DP table
  1. 1
    for i0i \gets 0 to mm do
  2. 2
    D[i][0]iD[i][0] \gets i
    delete all of A[1..i]A[1..i]
  3. 3
    for j0j \gets 0 to nn do
  4. 4
    D[0][j]jD[0][j] \gets j
    insert all of B[1..j]B[1..j]
  5. 5
    for i1i \gets 1 to mm do
  6. 6
    for j1j \gets 1 to nn do
  7. 7
    if A[i]=B[j]A[i] = B[j] then
  8. 8
    D[i][j]D[i1][j1]D[i][j] \gets D[i-1][j-1]
    free match
  9. 9
    else
  10. 10
    D[i][j]1+min(D[i1][j], D[i][j1], D[i1][j1])D[i][j] \gets 1 + \min\parens{D[i-1][j],\ D[i][j-1],\ D[i-1][j-1]}
    delete, insert, substitute
  11. 11
    return D[m][n]D[m][n]

Every cell is still , so the fill is , and the alignment is recovered by the same corner-to-origin traceback as LCS.

edit_distance.pypython
from typing import NamedTuple, Sequence, TypeVar

Element = TypeVar("Element")

class Edit(NamedTuple):
  """
    One step of an alignment script: the operation and the elements it acts on.\n
    `operation` is one of "match", "substitute", "delete", "insert"; for a\n
    delete `target` is None, for an insert `source` is None.\n
  """
  operation: str
  source: object
  target: object

def edit_distance(source: Sequence[Element], target: Sequence[Element]) -> int:
  """
    Minimum number of insert/delete/substitute edits taking `source` to\n
    `target`. Fills the full prefix table and reports the corner.\n
  """
  rows: int = len(source)
  cols: int = len(target)

  # distance[i][j] = edits to turn source[:i] into target[:j]. Borders count
  # a run of deletions (empty target) or insertions (empty source).
  distance: list[list[int]] = [[0 for _ in range(cols + 1)] for _ in range(rows + 1)]
  for i in range(rows + 1):
    distance[i][0] = i
  for j in range(cols + 1):
    distance[0][j] = j

  # fill: free on a match, else 1 + cheapest of delete / insert / substitute.
  for i in range(1, rows + 1):
    for j in range(1, cols + 1):
      if source[i - 1] == target[j - 1]:
        distance[i][j] = distance[i - 1][j - 1]
      else:
        distance[i][j] = 1 + min(
          distance[i - 1][j],      # delete source[i-1]
          distance[i][j - 1],      # insert target[j-1]
          distance[i - 1][j - 1],  # substitute
        )

  return distance[rows][cols]

def edit_script(
  source: Sequence[Element],
  target: Sequence[Element],
) -> list[Edit]:
  """
    A minimum-length sequence of edits transforming `source` into `target`.\n
    Builds the same table, then walks backwards from the corner, preferring a\n
    free match, then substitute, then delete, then insert. The returned list is\n
    in forward (left-to-right) order; its non-match steps number exactly the\n
    edit distance.\n
  """
  rows: int = len(source)
  cols: int = len(target)

  # border rows/cols cost a run of deletions or insertions.
  distance: list[list[int]] = [[0 for _ in range(cols + 1)] for _ in range(rows + 1)]
  for i in range(rows + 1):
    distance[i][0] = i
  for j in range(cols + 1):
    distance[0][j] = j

  # fill: free on a match, else 1 + cheapest of delete / insert / substitute.
  for i in range(1, rows + 1):
    for j in range(1, cols + 1):
      if source[i - 1] == target[j - 1]:
        distance[i][j] = distance[i - 1][j - 1]
      else:
        distance[i][j] = 1 + min(
          distance[i - 1][j],
          distance[i][j - 1],
          distance[i - 1][j - 1],
        )

  # walk back from the corner, preferring match, substitute, delete, insert.
  script: list[Edit] = []
  row, col = rows, cols
  while row > 0 or col > 0:
    if (
      row > 0
      and col > 0
      and source[row - 1] == target[col - 1]
      and distance[row][col] == distance[row - 1][col - 1]
    ):
      script.append(Edit("match", source[row - 1], target[col - 1]))
      row -= 1
      col -= 1
    elif (
      row > 0
      and col > 0
      and distance[row][col] == distance[row - 1][col - 1] + 1
    ):
      script.append(Edit("substitute", source[row - 1], target[col - 1]))
      row -= 1
      col -= 1
    elif row > 0 and distance[row][col] == distance[row - 1][col] + 1:
      script.append(Edit("delete", source[row - 1], None))
      row -= 1
    else:
      script.append(Edit("insert", None, target[col - 1]))
      col -= 1

  script.reverse()
  return script

Compare this against LCS line by line. Both index subproblems by prefix pairs; both branch on whether the last characters match; both fill an table where each entry reads its left, upper, and upper-left neighbors; both run in . The only differences are the costs and the optimization direction: LCS maximizes matched characters, edit distance minimizes edits. Sequence DPs are a single template, parameterized by what a match earns and a mismatch costs. Recognize the template and a whole family of problems (LCS, edit distance, sequence alignment, longest common substring, and string matching) falls to the same code.

Alignment, bioinformatics, and the quadratic wall

The LCS/edit-distance template is the single most consequential dynamic program in applied computing, because it is sequence alignment. Needleman and Wunsch (1970) introduced the global-alignment DP for comparing protein and nucleotide sequences; Smith and Waterman (1981) adapted it to local alignment (the best-matching substring pair, by clamping the score at and tracking the global maximum, the same move the maximum-subarray DP makes). These two recurrences are the foundation of computational biology, and the diff utility, git's merge machinery, spell checkers, and DNA read-mapping all descend from the same table. Gotoh (1982) refined the model with affine gap penalties — charging a large cost to open a gap and a small cost to extend it, which needs three coupled tables but stays — because a single long insertion is biologically more plausible than many scattered ones.

The catch is scale. A table is fine for two short strings but ruinous for two human chromosomes, and the Hirschberg linear-space trick (noted above) fixes the memory but not the time. Whether the time can be beaten is now settled conditionally: Backurs and Indyk (2015) and Bringmann and Künnemann (2015) proved that edit distance and LCS admit no strongly subquadratic algorithm unless the Strong Exponential Time Hypothesis is false. So the quadratic table is a genuine wall, and the practical response has been to give up exactness: heuristic aligners like BLAST (Altschul et al., 1990) and FASTA seed on short exact matches and extend them, trading a guarantee of optimality for the speed that made genome-scale search possible. The abstract define the subproblem, fill the table discipline of this lesson is, in this one instance, a multi-billion-dollar tool.5

Takeaways

  • Index sequence subproblems by prefixes: is the key to LCS, after the Step 0 move of solving for length first.
  • The recurrence is a over three cases (drop in case 1, drop in case 2, or extend the diagonal by when in case 3) over a base case of for an empty prefix.
  • Prove it by induction on in two directions: (build a witness subsequence) and (every common subsequence fits one of the three cases).
  • Fill the table in ; reconstruct in a second pass, walking backwards from and emitting a character on every diagonal match.
  • Only the length is needed? Two rows give space, at the cost of losing the reconstruction.
  • Edit distance is the same dynamic program: same prefix subproblems, same table shape, same , minimizing edits instead of maximizing matches. Sequence DP is one reusable template.

Footnotes

  1. Skiena, §10 — Dynamic Programming: the longest common subsequence as a similarity measure underlying diff and sequence alignment.
  2. CLRS, Ch. 15 — Dynamic Programming: the LCS recurrence obtained by examining the last characters of each prefix.
  3. Erickson, Ch. 3 — Dynamic Programming: Hirschberg's divide-and-conquer computes an optimal alignment in linear space by recursing on the midpoint column, keeping the time bound.
  4. Erickson, Ch. 3 — Dynamic Programming: edit (Levenshtein) distance as the minimum-cost insert/delete/substitute alignment filling an table.
  5. Needleman & Wunsch (1970) global and Smith & Waterman (1981) local sequence alignment; Gotoh (1982) affine gaps. Backurs & Indyk (2015) and Bringmann & Künnemann (2015): no strongly subquadratic edit distance / LCS under SETH — the practical reason heuristic aligners like BLAST (Altschul et al., 1990) exist.
Practice

╌╌ END ╌╌