Dynamic Programming

We solve the problem of finding the longest common subsequence of \(x=x_{1...m}\) and \(y=y_{1...n}\) by taking the best of the three possible cases:

1. The longest common subsequence of the strings \(x_{1...m-1}\) and \(y_{1...n}\)
   
2. The longest common subsequence of the strings \(x_{1...m}\) and \(y_{1...n-1}\)
   
3. If \(x_{m}\) is the same as \(y_{n}\), the longest common subsequence of the strings \(x_{1...m-1}\) and \(y_{1...n-1}\), followed by the common last character.
   

用 \(c[i][j]\) 记录序列 \(x_{1...i}\) 和 \(y_{1...j}\) 的最长公共子序列的长度，则有递归关系：\[c[i][j] = \begin{cases} 0 & \text{if} \; i=0 \; \text{or} \; j=0 \\ c[i-1][j-1] +1 & \text{if} \; i,j>0; x_i = y_j \\ \max \{c[i][j-1], c[i-1][j] \} & \text{if} \; i,j>0; x_i \neq y_j\end{cases}\]
利用这个递归式能算出 \(c[m][n]\) ，这就是最长公共子序列的长度。

可直接用程序求解上面的递归公式：

#include<stdio.h>
#include<string.h>

#define MAX(a,b) (((a)>(b))?(a):(b))

/* call it like this lcs_len(s, t, strlen(s)-1, strlen(t)-1) */
int lcs_len(char *s, char *t, int s_len, int t_len) {
    if(s_len < 0 || t_len < 0) {
        return 0;
    }

    if(s[s_len] == t[t_len]) {
        return lcs_len(s, t, s_len-1, t_len-1)+1;
    } else {
        return MAX(lcs_len(s, t, s_len-1, t_len), lcs_len(s, t, s_len, t_len-1));
    }
}

int main(void) {
    char a[]="ABCBDAB";
    char b[]="BDCABA";

    printf("length: %d\n", lcs_len(a, b, strlen(a)-1, strlen(b)-1));    // 输出 length: 4
    return 0;
}

上面这种递归解法的效率不高，会有过多的函数调用，而且有重复的计算。

前面的递归解法效率不高。下面我们设计一个表格化的、自底向上基于前面递归式的动态规化算法。其伪代码如下（摘自《算法导论（第 2 版）》）：

LCS-LENGTH(X, Y)
m ← length[X]
n ← length[Y]
for i ← 1 to m
  do c[i, 0] ← 0
for j ← 0 to n
  do c[0, j] ← 0
for i ← 1 to m
  do for j ← 1 to n
       do if xi = yj then
               c[i, j] ← c[i - 1, j - 1] + 1
               b[i, j] ← "↖"
          else if c[i - 1, j] ≥ c[i, j - 1] then
               c[i, j] ← c[i - 1, j]
               b[i, j] ← "↑"
          else c[i, j] ← c[i, j - 1]
               b[i, j] ← "←"
return c and b

上面伪代码中， \(c[i][j]\) 记录序列 \(x_{1...i}\) 和 \(y_{1...j}\) 的最长公共子序列的长度。 \(b[i][j]\) 记录着用来输出最优解的辅助信息，当 \(b[i][j]\) 等于“↖”时，意味着 \(x_i = y_j\) 是 LCS 中的一个元素，当 \(b[i][j]\) 等于“↑”时，意味着 \(x_i\) 不是 LCS 中的元素，当 \(b[i][j]\) 等于“←”时，意味着， \(y_j\) 不是 LCS 中的元素。

如果只要求最长公共子序列的长度，则可以不记录 \(b[i][j]\) 。
下面利用 \(b[i][j]\) 输出最长公共子序列。
打印最长公共子序列的伪代码如下（摘自《算法导论（第 2 版）》）：

PRINT-LCS(b, X, i, j)
if i = 0 or j = 0
    then return
if b[i, j] = "↖"
    then PRINT-LCS(b, X, i - 1, j - 1)
         print xi
else if b[i, j] = "↑"
    then PRINT-LCS(b, X, i - 1, j)
else PRINT-LCS(b, X, i, j - 1)

遇到“↖”就打印出 \(x_i\) ，同时把 X、Y 都减去一个字符，再递归调用打印程序；遇到“↑”就把 X 减去一个字符，再递归调用打印程序；遇到“←”就把 Y 减去一个字符，再递归调用打印程序。

以 X=BDCABA，Y=ABCBDAB 为例，完成计算后， \(c[i][j]\) （下图中每个方格中的数字）和 \(b[i][j]\) （下图中每个方格中的小箭头）如图 1 所示：

图 1 中，从左上角到右下角的大箭头表示 \(c[i][j]\) 和 \(b[i][j]\) 中元素的计算顺序。

用动态规划解决最长公共子序列问题，其空间复杂度比较高，为 \(O(mn)\) 。
如果当两个序列达到上万个节点时，内存消耗就成为了大问题。如确定一种致病基因的基本方法就是将该疾病患者的 DNA 链与健康者的 DNA 链相比较，找出其中的相同部分与不同部分，再进行分析。这种基因链的比较就可以用 LCS 算法来解决，但 DNA 序列一般都很长。所以动态规划方法不太实用。

针对动态规划空间复杂度高的缺陷，求解 LCS 问题有几种线性空间复杂度的其它算法，如 Hirschberg’s algorithm 和 Nakatsu algorithm 。

这是一个很有名的一个问题，在 Mark Allen Weiss 的《数据结构与算法分析——C 语言描述(原书第 2 版)》2.4.3 节中介绍了四种方法，《编程之美》2.14 节也有介绍，《编程珠玑(第二版)》第 8 章也介绍了该问题。

动态归划的思路：设 \(b[j]\) 是这些以 \(A[j]\) 结尾的众多子序列中和最大的一个，则原始问题的解就是 \(b[j], 1≤j≤n\) 中最大的一个。

//用动态规划算法解决最大子序和问题。
#include<iostream>

using std::cout;
using std::endl;

//返回值为子数组的最大和。
//A[]为数组，N为数组元素个数，start和end分别返回子数组开始、结束位置。
int max_sub_array(const int A[], int N, int& start, int& end) {
    int b = A[0];
    int sum = A[0];
    int start_tmp = 0;
    int end_tmp = 0;    //start_tmp和end_tmp用来辅助定位start、end。
    for (int j = 1; j < N; j++) {
        if (b > 0) {
            b = b + A[j];
            end_tmp = j;
        } else {
            b = A[j];
            start_tmp = end_tmp = j;
        }
        if (sum < b) {
            sum = b;
            start = start_tmp;
            end = end_tmp;
        }
    }
    return sum;
}

int main() {
    int A[] = { -2, 11, -4, 13, -5, -2 };
    int start, end;
    cout << max_sub_array(A, 6, start, end) << endl;
    cout << "start: A" << start << endl;  // 这个例子会输出 start: A1
    cout << "end: A" << end << endl;      // 这个例子会输出 end: A3
}

问题描述：
给定一个 \(M \times N\) 的矩阵，请找出此矩阵的一个子矩阵，并且此矩阵的各个元素的和最大，输出这个子矩阵和最大值。
例如，在如下这个矩阵中：

 0 -2 -7  0
 9  2 -6  2
-4  1 -4  1
-1  8  0 -2

拥有最大和的子矩阵为：

 9 2
-4 1
-1 8

其和为 15。

POJ（北京大学程序在线评测系统）上有这个题：http://poj.org/problem?id=1050

说明：最大子矩阵问题可以看做是最大子序和问题的推广，它也是动态规划的经典实例。具体分析略。

用 ed(i,j) 表示字符串 A[0,i] 和 B[0,j] 的编辑距离。假设 A[0,i] 和 B[0,j] 的形式分别为 str1a 和 str2b 。

如果 a==b，则 ed(i,j) = ed(i-1,j-1)。
如果 a != b，则 ed(i,j) 是下面三个可能值中最小的那个：
(1) 如果 a 替换为 b 时，可使 A 变为 B，则 ed(i,j) = ed(i-1,j-1) + 1；
(2) 如果在 a 后面插入字符 d，可使 A 变为 B，则 ed(i,j) = ed(i,j-1) + 1；
(3) 如果将 a 删除，可使 A 变为 B，则 ed(i,j) = ed(i-1,j) + 1。

显然，当 i=0 时，ed(i,j) = j；当 j=0 时，ed(i,j)=i。

总结，可得编辑距离的求解公式：
\[ed(i,j) = \begin{cases} i & \text{if} \; j = 0 \\ j & \text{if} \; i = 0 \\ ed(i-1,j-1) & \text{if} \; i,j > 0; a_i = b_j \\ \min\{ed(i-1,j-1), ed(i,j-1), ed(i-1,j)\} + 1 & \text{if} \; i,j > 0; a_i \ne b_j \\ \end{cases}\]

用如下递归程序可直观地求解上面公式。

// len_s and len_t are the number of characters in string s and t respectively
int EditDistance(char *s, int len_s, char *t, int len_t)
{
  int cost;

  /* base case: empty strings */
  if (len_s == 0) return len_t;
  if (len_t == 0) return len_s;

  /* test if last characters of the strings match */
  if (s[len_s-1] == t[len_t-1]) {
    return EditDistance(s, len_s - 1, t, len_t - 1);
  } else {
    int minimum = min(EditDistance(s, len_s - 1, t, len_t), EditDistance(s, len_s, t, len_t - 1));
    return min(minimum, EditDistance(s, len_s - 1, t, len_t - 1)) + 1;
  }
}

上面程序效率不高，会重复计算之前计算过的值。

递归解法效率太低，动态规则解法如下：

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

int EditDistance(const string & word1, const string & word2) {
  int n = word1.length();
  int m = word2.length();

  int dp[n+1][m+1];

  for(size_t i=0; i<=n; i++){
    dp[i][0] = i;
  }
  for(size_t j=0; j<=m; j++){
    dp[0][j] = j;
  }

  for(size_t i=1; i<=n; i++){
    for(size_t j=1; j<=m; j++){
      if(word1[i-1] == word2[j-1]){
        dp[i][j] = dp[i-1][j-1];
      }else{
        dp[i][j] = min(dp[i-1][j-1], min(dp[i-1][j],dp[i][j-1])) + 1;
      }
    }
  }

  return dp[n][m];
}

int main() {
  string A = "GGATCGA";
  string B = "GAATTCAGTTA";
  cout << EditDistance(A,B) << endl;  // 输出 5
}

说明，上面动态规划算法的时间复杂度为 \(O(nm)\) ，空间复杂度也为 \(O(nm)\) 。如果仅是求解编辑距离，不需要求解从字符串 A 到字符串 B 的具体转换过程，则可以对空间作进一步优化，优化过程可参考：https://en.wikipedia.org/wiki/Levenshtein_distance#Iterative_with_two_matrix_rows

这个问题，用穷举法其时间复杂度太高，我们尝试用动态规划方法求解。

(1)、寻找最优子结构。
假设 \(A_iA_{i+1} \cdots A_j\) 的最优加全部括号方案的分割点在 \(A_k\) 与 \(A_{k+1}\) 之间。那么，继续对“前缀”子链 \(A_iA_{i+1} \cdots A_k\) 进行括号化时，我们应该直接采用独立求解它时所得的最优方案。类似地，对于子链 \(A_{k+1}A_{k+2} \cdots A_j\) 我们也有相似的结论：在原问题 \(A_iA_{i+1} \cdots A_j\) 的最优括号化方案中，对子链 \(A_{k+1}A_{k+2} \cdots A_j\) 进行括号化的方法，就是它自身的最优括号化方案。

(2)、建立递归关系。
设 \(m[i][j]\) 为 \(A_iA_{i+1} \cdots A_j, (1≤i≤j≤n)\) 的加全部括号所需的标量乘法次数的最小值，则原问题的最优值为 \(m[1][n]\) 。
当 \(i = j\) 时，为单一矩阵，无需计算，因此 \(m[i][i]=0\) 。
当 \(i < j\) 时，我们假设 \(A_iA_{i+1} \cdots A_j\) 的最优加全部括号方案的分割点在 \(A_k\) 与 \(A_{k+1}\) 之间，那么 \(m[i][j]\) 就等于 \(A_iA_{i+1} \cdots A_k\) 的代价加上 \(A_{k+1}A_{k+2} \cdots A_j\) 的代价，再加上两者相乘的代价。其中两者相乘的代价为 \(p_{i-1}p_kp_j\) （设矩阵 \(A_i\) 的大小为 \(p_{i-1} \times p_{i}\) ）。因此，我们得到：
\[m[i][j] = m[i][k] + m[k+1][j] + p_{i-1}p_kp_j\]
上面递归公式假设最优分割点 \(k\) 是已知的，但实际上我们是不知道的。不过， \(k\) 的取值只可能是 \(k=i, i+1, \cdots, j-1\) ，由于最优分割点必在其中，我们只需要检查所有可能情况，找到最优者即可。
综上所述，有：
\[m[i][j] = \begin{cases} 0 & \text{if} \; i=j \\ \min_{i \le k < j} \{m[i][k] + m[k+1][j] + p_{i-1}p_kp_j \} & \text{if} \; i< j \end{cases}\]

不过，上面式子中，并未提供足够的信息来构造最优解。为了有助于跟踪如何构造一个最优解，定义 \(s[i][j]\) 为使 \(m[i][j] = m[i][k] + m[k+1][j] + p_{i-1} p_kp_j\) 成立时的 \(k\) 值，也就是分裂处的值。

//矩阵链乘法的动态规划算法实现。
//m[i][j]和s[i][j]的下标都从1开始，这是为了和前面的描述尽量一致。
//代码参考了《计算机算法设计与分析(第3版, 王晓东编)》3.1节。
#include<iostream>
#include<iomanip>    //使用setw()时需要
using std::cout;
using std::endl;
using std::setw;

#define N 6  //矩阵个数

//p为矩阵维数数组，n为矩阵个数。
void MatrixChain(int *p, int n, int m[][N + 1], int s[][N + 1]) {
    for (int i = 1; i <= n; i++) {
        m[i][i] = 0; //i=j时，m[i][j]为0。
    }
    for (int l = 2; l <= n; l++) {  //子序列长度由2到n递增
        for (int i = 1; i <= n - l + 1; i++) {  //注意i的上限
            int j = i + l - 1;
            m[i][j] = m[i + 1][j] + p[i - 1] * p[i] * p[j];
            s[i][j] = i;
            for (int k = i + 1; k < j; k++) {
                if (m[i][j] > m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j]) { //根据递归表达式计算m[i][j]
                    m[i][j] = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
                    s[i][j] = k;  //记录s[i][j]
                }
            }
        }
    }
}

//递归打印加上括号的矩阵链。
//s为保存着辅助信息的数组。i,j为矩阵链范围。
void Print_parens(int s[][N + 1], int i, int j) {  //递归地构造最优解
    int k = s[i][j];
    if (i != j) {
        if (!(i == 1 && j == N)) {
            cout << "(";
        }
        Print_parens(s, i, k);
        Print_parens(s, k + 1, j);
        if (!(i == 1 && j == N)) {
            cout << ")";
        }
    } else {
        cout << "A" << i;
    }
}

int main() {
    int p[N + 1] = { 30, 35, 15, 5, 10, 20, 25 };
    int m[N + 1][N + 1] = { 0 };  //不初赋值不会影响结果，但m的“下三角数组”会为随机数。
    int s[N + 1][N + 1] = { 0 };  //不初赋值不会影响结果，但s的“下三角数组”会为随机数。

    MatrixChain(p, N, m, s);
    Print_parens(s, 1, N);

    //打印m数组和s数组。
    cout << endl << endl;
    cout << "m=" << endl;
    for (int i = 1; i <= N; i++) { //由于前面约定数组从1开始。
        for (int j = 1; j <= N; j++) {
            cout << setw(8) << m[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
    cout << "s=" << endl;
    for (int i = 1; i <= N; i++) {
        for (int j = 1; j <= N; j++) {
            cout << setw(8) << s[i][j] << " ";
        }
        cout << endl;
    }
    return 0;
}

上面代码的输出：

(A1(A2A3))((A4A5)A6)

m=
       0    15750     7875     9375    11875    15125
       0        0     2625     4375     7125    10500
       0        0        0      750     2500     5375
       0        0        0        0     1000     3500
       0        0        0        0        0     5000
       0        0        0        0        0        0

s=
       0        1        1        3        3        3
       0        0        2        3        3        3
       0        0        0        3        3        3
       0        0        0        0        4        5
       0        0        0        0        0        5
       0        0        0        0        0        0