Sequence Alignment DP: Longest Common Subsequence (LCS) & Wagner-Fischer Edit Distance

1. The Alignment Grid DAG: Wagner-Fischer Edit Distance & LCS

Given strings $A[0..n-1]$ and $B[0..m-1]$, prefix-pair Dynamic Programming defines subproblem $(i, j)$ on prefixes $A[0..i-1]$ and $B[0..j-1]$ across an $(n+1) \times (m+1)$ grid DAG with three incoming edges into $(i, j)$:

• Wagner–Fischer Levenshtein Edit Distance (1974): Minimum insertions, deletions, and substitutions to transform $A$ into $B$: $$dp[i][j] = \begin{cases} dp[i-1][j-1] & \text{if } A[i-1] = B[j-1] \\ 1 + \min\big(\underbrace{dp[i-1][j]}_{\text{Delete } A[i-1]},\; \underbrace{dp[i][j-1]}_{\text{Insert } B[j-1]},\; \underbrace{dp[i-1][j-1]}_{\text{Substitute}}\big) & \text{if } A[i-1] \neq B[j-1] \end{cases}$$ with base cases $dp[i][0] = i$ and $dp[0][j] = j$.
• Deletion-Insertion Edit Distance vs. LCS Duality: When substitutions are forbidden (or cost $2$), transforming $A$ into $B$ deletes $n - |\text{LCS}|$ characters and inserts $m - |\text{LCS}|$ characters: $$d_{\text{del/ins}}(A, B) = n + m - 2 \cdot |\text{LCS}(A, B)|$$

2. Needleman-Wunsch (Global) vs. Smith-Waterman (Local) Sequence Alignment

In computational genomics, mismatches carry penalty $-\mu$ and gaps carry penalty $-\delta$ while matches score $+\sigma$:

• Needleman–Wunsch (Global Alignment): Aligns full strings end-to-end; initializes $dp[i][0] = -i\delta, dp[0][j] = -j\delta$ and returns $dp[n][m]$.
• Smith–Waterman (Local Alignment): Finds the highest-scoring aligned substrings by adding a $0$ floor to every state transition ($dp[i][j] = \max(0, \dots)$)—allowing the alignment path to start fresh at any $(i, j)$ with score $0$—and returns $\max_{i,j} dp[i][j]$.

3. Complexity Matrix

4. Reference C++20 Implementation

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <string_view>
#include <unordered_map>
#include <vector>

// 1. Wagner-Fischer Levenshtein Edit Distance in O(N*M) Time and O(min(N, M)) Space
int edit_distance_space_optimized(std::string_view a, std::string_view b) {
    if (a.size() < b.size()) std::swap(a, b);
    const int n = static_cast<int>(a.size());
    const int m = static_cast<int>(b.size());

    std::vector<int> dp(m + 1);
    for (int j = 0; j <= m; ++j) dp[j] = j;

    for (int i = 1; i <= n; ++i) {
        int prev_diag = dp[0]; // Holds dp[i - 1][j - 1]
        dp[0] = i;             // Base case dp[i][0] = i deletions
        for (int j = 1; j <= m; ++j) {
            int saved_up = dp[j];
            if (a[i - 1] == b[j - 1]) {
                dp[j] = prev_diag;
            } else {
                dp[j] = 1 + std::min({ saved_up, dp[j - 1], prev_diag });
            }
            prev_diag = saved_up;
        }
    }
    return dp[m];
}

// 2. Ukkonen's O(N * K) Banded Edit Distance (Returns distance if <= k, else k + 1)
int banded_edit_distance(std::string_view a, std::string_view b, int k) {
    const int n = static_cast<int>(a.size());
    const int m = static_cast<int>(b.size());
    if (std::abs(n - m) > k) return k + 1;

    constexpr int INF = 1'000'000'000;
    std::vector<int> prev(m + 1, INF), curr(m + 1, INF);
    for (int j = 0; j <= std::min(m, k); ++j) prev[j] = j;

    for (int i = 1; i <= n; ++i) {
        int lo = std::max(1, i - k);
        int hi = std::min(m, i + k);
        curr[lo - 1] = (lo == 1) ? i : INF;

        int row_min = curr[lo - 1];
        for (int j = lo; j <= hi; ++j) {
            int cost = (a[i - 1] == b[j - 1]) ? 0 : 1;
            curr[j] = std::min({
                prev[j] + 1,         // Delete
                curr[j - 1] + 1,     // Insert
                prev[j - 1] + cost   // Match / Substitute
            });
            row_min = std::min(row_min, curr[j]);
        }
        if (hi < m) curr[hi + 1] = INF; // Guard boundary outside band
        if (row_min > k) return k + 1;   // Early exit: entire band exceeded threshold k!
        std::swap(prev, curr);
    }
    return std::min(k + 1, prev[m]);
}

// 3. Hunt-Szymanski O((N + M) log N) LCS when elements of a are distinct
int lcs_distinct_via_lis(const std::vector<int>& a, const std::vector<int>& b) {
    std::unordered_map<int, int> pos_in_a;
    pos_in_a.reserve(a.size() * 2);
    for (int i = 0; i < static_cast<int>(a.size()); ++i) {
        pos_in_a[a[i]] = i;
    }
    std::vector<int> tails;
    for (int x : b) {
        auto it = pos_in_a.find(x);
        if (it == pos_in_a.end()) continue;
        int idx = it->second;
        auto lb = std::lower_bound(tails.begin(), tails.end(), idx);
        if (lb == tails.end()) tails.push_back(idx);
        else *lb = idx;
    }
    return static_cast<int>(tails.size());
}