Longest Common Subsequence - ricorrenza, tabella e ricostruzione

La LCS è il secondo grande esempio di programmazione dinamicaSchema in cinque passi: sottostruttura ottima, ricorrenza sui costi, informazione addizionale, calcolo bottom-up o memoizzato, ricostruzioneProgrammazione dinamica - Fibonacci, memoizzazione e schema del paradigma → ed è il modello di tutti i problemi su stringhe che seguono.

Stringhe, prefissi, sottostringhe, sottosequenze

Sia $\Sigma$ un alfabeto finito. Una stringa è una sequenza finita $X = \langle x_1, \dots, x_m \rangle$ di simboli di $\Sigma$; $\varepsilon$ è la stringa vuota e $\Sigma^*$ l'insieme di tutte le stringhe finite (stella di Kleene).

• Prefisso $X_i = \langle x_1, \dots, x_i \rangle$, per $0 \le i \le m$, con $X_0 = \varepsilon$.
• Suffisso $X^i = \langle x_i, \dots, x_m \rangle$, per $1 \le i \le m + 1$, con $X^{m+1} = \varepsilon$.
• Sottostringa $X_{i..j} = \langle x_i, \dots, x_j \rangle$ per $1 \le i \le j \le m$ (simboli contigui); $X_{i..j} = \varepsilon$ se $i > j$.
• Sottosequenza: una stringa $Z = \langle z_1, \dots, z_k \rangle$ per cui esiste una successione di indici strettamente crescente $1 \le j_1 < j_2 < \dots < j_k \le m$ con $z_s = x_{j_s}$. Si dice che gli indici realizzano $Z$ in $X$. Una stessa $Z$ può avere più successioni di indici.

Esempio. $X = \langle A, B, C, B, B, D \rangle$ e $Z = \langle A, C, B \rangle$: $Z$ è sottosequenza di $X$, realizzata da $\langle 1, 3, 4 \rangle$ e anche da $\langle 1, 3, 5 \rangle$. Non è una sottostringa. Ogni sottostringa è una sottosequenza (con indici $i, i+1, \dots, j$), non viceversa.

Quante sono. Le sottostringhe sono $\Theta(m^2)$: ce ne sono $m$ di lunghezza $1$, $m - 1$ di lunghezza $2$, $\dots$, una di lunghezza $m$, più $\varepsilon$:

$$1 + \sum_{j=1}^{m}(m - j + 1) = 1 + \sum_{j'=1}^{m} j' = 1 + \frac{m(m+1)}{2}.$$

Le sottosequenze sono molte di più. Se i caratteri di $X$ sono tutti distinti, ogni sottoinsieme di $k$ indici di $\{1, \dots, m\}$ dà una sottosequenza distinta (la si ottiene ordinando gli indici), e i sottoinsiemi di taglia $k$ sono $\binom{m}{k}$ (Fattoriale e coefficienti binomialiFattoriale, permutazioni, disposizioni, combinazioni e coefficiente binomiale n su k, con il triangolo di Tartaglia.Fattoriale e coefficienti binomiali →). In totale

$$\sum_{k=0}^{m} \binom{m}{k} = 2^m$$

per il binomio di Newton (Binomio di NewtonLa formula per sviluppare (a+b)^n con i coefficienti binomiali.Binomio di Newton →); contando anche $\varepsilon$ (per $k = 0$) sono $2^m$. Per questo un approccio che enumera le sottosequenze di $X$ e le verifica su $Y$ è esponenziale.

Il problema

Definizione (LCS). Date $X = \langle x_1, \dots, x_m \rangle$ e $Y = \langle y_1, \dots, y_n \rangle$, determinare una stringa $Z = \langle z_1, \dots, z_k \rangle$ di lunghezza massima che sia sottosequenza sia di $X$ sia di $Y$.

Esempio. $X = \langle A, B, C, B, B, D \rangle$, $Y = \langle A, D, C, C, B, D \rangle$: la LCS è $Z = \langle A, C, B, D \rangle$, lunga $4$, realizzata in $X$ da $\langle 1, 3, 4, 6 \rangle$ e in $Y$ da $\langle 1, 4, 5, 6 \rangle$.

La lunghezza della LCS misura la somiglianza tra due stringhe; si usa in genomica (confronto di sequenze di DNA) e nella correzione automatica dei testi.

La LCS può non essere unica: $X = \langle A, B, C \rangle$ e $Y = \langle B, A, C \rangle$ hanno due LCS di lunghezza $2$, $\langle A, C \rangle$ e $\langle B, C \rangle$. Si cerca una LCS.

Proprietà di sottostruttura ottima

Idea: guardare l'ultimo simbolo. Se $X$ e $Y$ finiscono con lo stesso simbolo, quel simbolo si può sempre accoppiare; altrimenti uno dei due simboli finali non è usato. Questo riconduce un'istanza a istanze formate da coppie di prefissi, che sono $(m + 1)(n + 1) = \Theta(mn)$: uno spazio quadratico, non esponenziale come le sottosequenze. La generica sottoistanza è $(X_i, Y_j)$ con $0 \le i \le m$ e $0 \le j \le n$, e si vuole $\text{LCS}(X_m, Y_n) = \text{LCS}(X, Y)$.

Teorema (sottostruttura della LCS). Siano $X_i$, $Y_j$ prefissi di $X$ e $Y$ e sia $Z = \langle z_1, \dots, z_k \rangle$ una LCS di $X_i$ e $Y_j$.

1. Se $i = 0$ o $j = 0$, allora $Z = \varepsilon$.
2. Se $i > 0$ e $j > 0$:
   • (a) se $x_i = y_j$, allora $z_k = x_i = y_j$ e $Z_{k-1} = \langle z_1, \dots, z_{k-1} \rangle$ è una LCS di $X_{i-1}$ e $Y_{j-1}$;
   • (b) se $x_i \ne y_j$, allora $Z$ è una LCS di $X_{i-1}, Y_j$ oppure di $X_i, Y_{j-1}$, e quindi è la più lunga tra le due.

Dimostrazione. Il punto 1 è banale: un prefisso vuoto non ha sottosequenze non vuote.

Caso (a), $x_i = y_j$. Siano $i_1 < \dots < i_k \le i$ e $j_1 < \dots < j_k \le j$ gli indici che realizzano $Z$ in $X_i$ e in $Y_j$.

• $z_k = x_i$. Per assurdo sia $z_k \ne x_i = y_j$ (oppure $k = 0$). Allora $i_k \ne i$ e $j_k \ne j$ (altrimenti $z_k$ sarebbe $x_i$ o $y_j$), cioè $i_k < i$ e $j_k < j$. Ma allora $\langle Z, x_i \rangle$ è una sottosequenza comune, realizzata dagli indici $i_1 < \dots < i_k < i$ in $X_i$ e $j_1 < \dots < j_k < j$ in $Y_j$, ed è più lunga di $Z$: contraddizione con la massimalità di $Z$.
• $Z_{k-1}$ è una LCS di $X_{i-1}$ e $Y_{j-1}$. Gli indici $i_1 < \dots < i_{k-1} < i_k \le i$ sono tutti $\le i - 1$, e lo stesso vale in $Y$: quindi $Z_{k-1}$ è una sottosequenza comune di $X_{i-1}$ e $Y_{j-1}$. Se esistesse una sottosequenza comune $W$ di $X_{i-1}$ e $Y_{j-1}$ più lunga di $Z_{k-1}$, allora $\langle W, x_i \rangle$ sarebbe una sottosequenza comune di $X_i$ e $Y_j$ più lunga di $Z$ (è la tecnica del cut-and-paste: si taglia $Z_{k-1}$ e si incolla $W$), contro la massimalità di $Z$.

Caso (b), $x_i \ne y_j$. Il simbolo $z_k$ non può essere uguale sia a $x_i$ sia a $y_j$. Se $z_k \ne x_i$ (o $k = 0$), l'indice $i$ non è usato per realizzare $Z$ in $X_i$, quindi $Z$ è una sottosequenza comune di $X_{i-1}$ e $Y_j$. Se $z_k \ne y_j$, per lo stesso motivo $Z$ è sottosequenza comune di $X_i$ e $Y_{j-1}$. Almeno una delle due vale. Nel primo caso $Z$ deve essere una LCS di $(X_{i-1}, Y_j)$ (una sottosequenza comune più lunga, che è anche sottosequenza comune di $X_i$ e $Y_j$, batterebbe $Z$). Nel secondo, di $(X_i, Y_{j-1})$. Infine le LCS di $(X_{i-1}, Y_j)$ e di $(X_i, Y_{j-1})$ sono entrambe sottosequenze comuni di $X_i$ e $Y_j$, quindi nessuna è più lunga di $Z$: $Z$ è la più lunga delle due. $\square$

Ricorrenza sui costi e informazione addizionale

Sia $\ell(i, j) = |\text{LCS}(X_i, Y_j)|$. Dal teorema, per $0 \le i \le m$ e $0 \le j \le n$:

$$\ell(i, j) = \begin{cases} 0 & i = 0 \ \lor \ j = 0 \\ 1 + \ell(i-1, j-1) & i, j > 0 \ \land \ x_i = y_j \\ \max\{\ell(i-1, j),\ \ell(i, j-1)\} & i, j > 0 \ \land \ x_i \ne y_j \end{cases}$$

Informazione addizionale. Per i casi base la soluzione è sempre $\varepsilon$: non si salva nulla. Per gli altri ci sono solo tre modi di costruire $Z = \text{LCS}(X_i, Y_j)$ da sottoistanze: $\langle \text{LCS}(X_{i-1}, Y_{j-1}), x_i \rangle$, oppure $\text{LCS}(X_i, Y_{j-1})$, oppure $\text{LCS}(X_{i-1}, Y_j)$. Basta ricordare quale: una freccia

$$b(i, j) = \begin{cases} \nwarrow & x_i = y_j \\ \leftarrow & x_i \ne y_j \ \land \ Z = \text{LCS}(X_i, Y_{j-1}) \\ \uparrow & x_i \ne y_j \ \land \ Z = \text{LCS}(X_{i-1}, Y_j) \end{cases}$$

Il codice bottom-up

Si usano due tabelle $\ell[0..m, 0..n]$ e $b[0..m, 0..n]$. La cella $(i, j)$ dipende solo da $(i-1, j-1)$, $(i-1, j)$ e $(i, j-1)$. Per un ordine coerente basta che al momento di calcolare $(i, j)$ queste tre siano già pronte (si inizializzano riga $0$ e colonna $0$, che contengono i casi base). Funzionano sia la scansione per righe (riga $1$ da sinistra a destra, poi riga $2$, ...) sia quella per colonne: in entrambe, le tre celle da cui si dipende precedono $(i, j)$. Per righe è preferibile in linguaggi che memorizzano le matrici per riga, per l'uso della cache.

LCS(X, Y)
  m <- |X|; n <- |Y|
  for i <- 0 to m do l[i, 0] <- 0
  for j <- 1 to n do l[0, j] <- 0
  for i <- 1 to m do
    for j <- 1 to n do
      if x_i = y_j then
        l[i, j] <- 1 + l[i-1, j-1]
        b[i, j] <- '↖'
      else if l[i-1, j] >= l[i, j-1] then
        l[i, j] <- l[i-1, j]
        b[i, j] <- '↑'
      else
        l[i, j] <- l[i, j-1]
        b[i, j] <- '←'
  return (l[m, n], b)

Costo. Contando solo i confronti tra simboli (modello di costo usuale per le stringhe), ce n'è uno per iterazione:

$$T(m, n) = \sum_{i=1}^{m} \sum_{j=1}^{n} 1 = m \cdot n = \Theta(mn),$$

quadratico (con $m = n$). L'algoritmo restituisce la lunghezza; per ricostruire la LCS bisogna restituire anche tutta la tabella $b$.

Ricostruzione. Si parte da $(m, n)$ e si segue la freccia: a ogni passo almeno una delle due coordinate scende, e $i + j$ diminuisce di almeno $1$. Quindi le chiamate sono $O(m + n)$.

PRINT_LCS(X, b, i, j)               { prima chiamata: PRINT_LCS(X, b, m, n) }
  if i = 0 or j = 0 then return
  if b[i, j] = '↖' then
    PRINT_LCS(X, b, i-1, j-1)       { prima la LCS dei prefissi... }
    print(x_i)                      { ...poi il simbolo comune }
  else if b[i, j] = '←' then PRINT_LCS(X, b, i, j-1)
  else PRINT_LCS(X, b, i-1, j)

Esempio eseguito

$X = \langle A, B, C, B, B, D \rangle$ e $Y = \langle A, D, C, C, B, D \rangle$. In grassetto il percorso della ricostruzione (partendo da $(6, 6)$). Ogni cella ha il valore $\ell(i, j)$ e la freccia $b(i, j)$.

Alcune celle a mano:

• $(1, 1)$: $x_1 = A = y_1$, quindi $\ell = 1 + \ell(0, 0) = 1$, freccia $\nwarrow$.
• $(2, 5)$: $x_2 = B = y_5$, quindi $\ell = 1 + \ell(1, 4) = 1 + 1 = 2$, freccia $\nwarrow$.
• $(3, 4)$: $x_3 = C = y_4$, quindi $\ell = 1 + \ell(2, 3) = 1 + 1 = 2$ (la cella $(3, 3)$ ha lo stesso valore, perché anche $y_3 = C$).
• $(4, 4)$: $x_4 = B \ne C = y_4$. Si confrontano $\ell(3, 4) = 2$ e $\ell(4, 3) = 2$: sono uguali, e per il $\ge$ del codice vince $\uparrow$.

Ricostruzione: $(6,6)$ $\nwarrow$ stampa $D$ (dopo), poi $(5,5)$ $\nwarrow$ $B$, poi $(4,4)$ $\uparrow$ sale a $(3,4)$ $\nwarrow$ $C$, poi $(2,3)$ $\uparrow$ e $(1,3)$ $\leftarrow$, $(1,2)$ $\leftarrow$, $(1,1)$ $\nwarrow$ $A$. Letti dall'alto in basso i simboli stampati sono $A, C, B, D$: $\ell(6, 6) = 4$ e la LCS è $\langle A, C, B, D \rangle$. Tutta la tabella è stata confrontata con una forza bruta (enumerazione di tutte le sottosequenze) su $400$ coppie casuali di stringhe corte.

Una LCS non è unica e il risultato dipende dal criterio nei casi di parità. Per $X = \langle A, B, C \rangle$ e $Y = \langle B, A, C \rangle$ nella cella $(2, 2)$ si ha $\ell(1, 2) = \ell(2, 1) = 1$: con il $\ge$ del codice si sceglie $\uparrow$ e la ricostruzione dà $\langle A, C \rangle$; scegliendo $\leftarrow$ nella parità si otterrebbe $\langle B, C \rangle$.

def lcs(X, Y):
    m, n = len(X), len(Y)
    l = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if X[i - 1] == Y[j - 1]:
                l[i][j] = 1 + l[i - 1][j - 1]
            else:
                l[i][j] = max(l[i - 1][j], l[i][j - 1])
    i, j, z = m, n, []
    while i > 0 and j > 0:                  # ricostruzione senza la tabella b
        if X[i - 1] == Y[j - 1]:
            z.append(X[i - 1]); i -= 1; j -= 1
        elif l[i - 1][j] >= l[i][j - 1]:
            i -= 1
        else:
            j -= 1
    return l[m][n], "".join(reversed(z))

assert lcs("ABCBBD", "ADCCBD") == (4, "ACBD")

(La tabella $b$ si può anche non memorizzare: la freccia si ricava a posteriori confrontando le celle vicine, come fa il ciclo qui sopra.)

Il divide and conquer sulla stessa ricorrenza

Un algoritmo ricorsivo sulla stessa sottostruttura, al caso peggiore (per esempio $X$ e $Y$ senza simboli in comune), ha un nodo interno con lavoro unitario per ogni coppia con $m, n > 0$ e due chiamate:

$$T(m, n) = \begin{cases} 0 & m = 0 \ \lor \ n = 0 \\ T(m-1, n) + T(m, n-1) + 1 & \text{altrimenti} \end{cases}$$

Minorazione (srotolamento). Con $m = n$: $T(n, n) = T(n-1, n) + T(n, n-1) + 1 \ge 2T(n-1, n-1) + 1$ (perché $T$ è non decrescente in ciascun argomento). Srotolando $i$ volte, $T(n, n) \ge 2^i T(n-i, n-i) + 2^i - 1$, e con $i = n$ si ottiene $T(n, n) \ge 2^n - 1$.

Valore esatto: $T(m, n) = \binom{m+n}{m} - 1$. Per induzione sulla somma $m + n$: se $m = 0$ o $n = 0$, $\binom{m+n}{m} - 1 = 0$. Altrimenti, con la regola di Pascal,

$$T(m, n) = \Big[\tbinom{m+n-1}{m-1} - 1\Big] + \Big[\tbinom{m+n-1}{m} - 1\Big] + 1 = \tbinom{m+n}{m} - 1.$$

Per $m = n = 20$ sono $\binom{40}{20} - 1 = 137\,846\,528\,819$ nodi, contro le $400$ celle della tabella. La causa sono le sottoistanze ripetute: la stessa coppia di prefissi $(i, j)$ viene raggiunta da molti cammini.

Risparmiare spazio

Se interessa solo la lunghezza della LCS non serve la tabella $b$, e per $\ell$ bastano due righe: la riga $i$ dipende solo dalla riga $i - 1$ e dalla riga $i$ stessa. Scambiando i ruoli in modo che $|Y| \le |X|$ si usa spazio $\Theta(\min\{m, n\})$ invece di $\Theta(mn)$.

LCS_LENGTH(X, Y)                    { |Y| <= |X| }
  m <- |X|; n <- |Y|
  for j <- 0 to n do prev[j] <- 0
  for i <- 1 to m do
    curr[0] <- 0
    for j <- 1 to n do
      if x_i = y_j then curr[j] <- 1 + prev[j-1]
      else curr[j] <- max(prev[j], curr[j-1])
    prev <- curr
  return prev[n]

(Con la ricostruzione bisogna conservare tutta la tabella, o ricalcolare a pezzi.) Il confronto con la versione a tabella intera è stato verificato su stringhe casuali.

L'algoritmo memoizzato

Si mantiene la ricorsione e si aggiunge la tabella globale. La sentinella è $-1$: è corretta perché $\ell \ge 0$. Le celle dei casi base valgono $0$.

INIT_LCS(X, Y)
  m <- |X|; n <- |Y|
  if m = 0 or n = 0 then return 0
  for i <- 0 to m do l[i, 0] <- 0
  for j <- 0 to n do l[0, j] <- 0
  for i <- 1 to m do
    for j <- 1 to n do l[i, j] <- -1
  return MEMO_LCS(X, Y, m, n)

MEMO_LCS(X, Y, i, j)
  if l[i, j] = -1 then
    if x_i = y_j then
      l[i, j] <- 1 + MEMO_LCS(X, Y, i-1, j-1)
    else
      a <- MEMO_LCS(X, Y, i-1, j)
      c <- MEMO_LCS(X, Y, i, j-1)
      l[i, j] <- max(a, c)
  return l[i, j]

• Correttezza: dalla sottostruttura e dalla scelta della sentinella.
• Costo al caso peggiore: si analizza l'albero delle chiamate. Per $1 \le i \le m$ e $1 \le j \le n$, la chiamata $(i, j)$ è un nodo interno al più una volta, a costo unitario (un confronto); le altre occorrenze restituiscono il valore in tabella. Quindi $T(m, n) \le m \cdot n$, uguale al bottom-up.
• Non sempre tutte le celle. Se le due stringhe sono simili (hanno lunghi tratti in comune), la discesa top-down tocca poche celle. Esempio: $X = \langle A, A, B, B, C \rangle$ e $Y = \langle B, B, C \rangle$. Le chiamate sono $(5,3) \to (4,2) \to (3,1) \to (2,0)$: ultimi simboli uguali ($C = C$, poi $B = B$, poi $B = B$), si scende in diagonale fino al caso base $(2, 0)$. Solo $3$ nodi interni (e un caso base) contro le $5 \times 3 = 15$ celle del bottom-up. Nelle applicazioni di genomica, dove le sequenze sono quasi uguali, il vantaggio è notevole. Controllato anche con Python: la discesa memoizzata tocca esattamente $(5,3), (4,2), (3,1)$.
• Per la coppia $ABCBBD$ e $ADCCBD$ dell'esempio eseguito il memoizzato calcola $16$ celle su $36$ (stringhe poco simili: il vantaggio è piccolo).

Domande d'esame

1. (Risposta numerica.) Siano $X = \langle a, f, r, a \rangle$ e $Y = \langle a, f, a \rangle$ e sia $l(i, j)$ la lunghezza della LCS di $X_i$ e $Y_j$. Calcolare $l(2, 2) + l(4, 3)$. È proprio uno dei quesiti della prima parte delle simulazioni: svolto in Esercizio - domande di prima parte su programmazione dinamica (simulazioni d'esame). Risposta: $2 + 3 = 5$.
2. (Risposta numerica.) Per $m = n = 10$, quanti confronti tra simboli esegue LCS? E quanti nodi interni al massimo ha l'albero del D&C puro? Traccia: $m \cdot n = 100$ confronti; $\binom{20}{10} - 1 = 184\,755$ nodi interni.
3. (Dimostrazione.) Se $x_i = y_j$ perché $\text{LCS}(X_i, Y_j)$ termina con questo simbolo? Traccia: altrimenti si potrebbe aggiungere il simbolo comune alla fine, ottenendo una sottosequenza comune più lunga (indici $i_k < i$ e $j_k < j$, perché $z_k \ne x_i, y_j$).