算法基础:动态规划（二）

关于动态规划，这里再举一个书上的经典例子。

求最长公共子序列（LCS）

问题描述：

令序列$X=x_1{x}_2{x}_3...x_m$$X = x_1x_2x_3...x_m$ ，序列 $Y=y_1{y}_2...y_k$$Y = y_1y_2...y_k$ 是$X$$X$ 的子序列，存在$X$$X$ 的一个严格递增下标序列<$i_1,i_2,...,i_k$$i_1,i_2,...,i_k$>，使得对于所有的$j=1,2,...,k$$j = 1,2,...,k$有$x_{i_j}=y_j$$x_{i_j} = y_j$。例如，$X=ABCBDAB,Y=BCDB$$X = ABCBDAB , Y = BCDB$，则$Y$$Y$是$X$$X$的一个子序列。

给定两个序列$X$$X$和$Y$$Y$，称序列$Z$$Z$是$X\text{和}Y$$X和Y$的公共子序列，是指$Z$$Z$同时是$X$$X$和$Y$$Y$的组序列。最长公共子序列问题定义为：给定序列$X=x_1{x}_2...x_m$$X = x_1x_2...x_m$和$Y=y_1{y}_2...y_n$$Y = y_1y_2...y_n$，求这两个序列的最长公共子序列 。以$X=ABCBDAB,Y=BCDB$$X = ABCBDAB , Y = BCDB$为例，两者的公共子序列为$Z=BCDB$$Z = BCDB$。

问题分析：

如果是匹配字符串问题的话，可以用经典KMP算法，可惜不是，求不得。如果使用暴力算法，穷举$X$$X$的所有子序列，并与$Y$$Y$序列比较，其时间复杂度可以达到$O(2^{m}n)$$O(2^mn)$，对于长序列来说，这个代价是不可接受的。那我们使用动态规划按照以下步骤来解决这个问题:

1）刻画最长公共子序列问题的最优子结构。

设$X=x_1{x}_2...x_m$$X=x_1x_2...x_m$和$Y=y_1{y}_2...y_n$$Y=y_1y_2...y_n$是两个序列，$Z=z_1{z}_2...z_k$$Z=z_1z_2...z_k$是$X$$X$和$Y$$Y$的一个最长公共子序列。

（1）如果$x_m=y_n$$x_m=y_n$，那么$z_k=x_m=y_n$$z_k=x_m=y_n$，且$Z_{k-1}$$Z_{k-1}$是$X_{m-1}$$X_{m-1}$和$Y_{n-1}$$Y_{n-1}$的一个最长公共子序列。

（2）如果$x_m\ne y_n$$x_m{\ne}y_n$，那么$z_k\ne x_m$$z_k{\ne}x_m$，蕴含着$Z$$Z$是$X_{m-1}$$X_{m-1}$和$Y$$Y$的一个最长公共子序列。

（3）如果$x_m\ne y_n$$x_m{\ne}y_n$，那么$z_k\ne y_n$$z_k{\ne}y_n$，蕴含着$Z$$Z$是$X$$X$和$Y_{n-1}$$Y_{n-1}$的一个最长公共子序列。

2） 递归定义最优解的值

3） 代码演示

C++

​x
94
1
#include <algorithm>
2
#include <iostream>
3
#include <vector>
4

5
using std::cout;
6
using std::endl;
7
using std::string;
8
using std::vector;
9

10

11

12
//算法主体
13
void Match(string x,string y,int lenx,int leny){
14

15
    //建立图表
16
    vector<vector<std::pair<int,int>>> diagram(lenx+1);
17
    for(auto & v:diagram){
18
        v.resize(leny+1,{0,0}); // 初始化
19
    }
20
    
21

22
    // pair的第二位用于表示当前字符是否属于X和Y的公共子序列
23
    // 1，2分别表示“当前X的字符不在公共子序列中，且当前X_i与Y_j的LCS和X_(i-1)与Y_j的LCS一样",
24
    // ”当前Y的字符不在公共子序列中，且当前X_i与Y_j的LCS和X_i与Y_(j-1)的LCS一样"
25
    // 3表示当前X与Y字符相等，且X_i与Y_j的LCS等于X_(i-1)与Y_(j-1)的LCS+1
26
    for(int i = 1; i < lenx+1 ; ++i){
27
        for(int j = 1; j < leny+1 ; ++j){
28
            if(x[i-1] == y[j-1]){
29
                diagram[i][j].first = diagram[i-1][j-1].first+1;
30
                diagram[i][j].second = 3;
31
            }
32
            else{
33
                if(diagram[i-1][j].first > diagram[i][j-1].first){
34
                    diagram[i][j].first = diagram[i-1][j].first;
35
                    diagram[i][j].second = 1;
36
                } 
37
                else{
38
                    diagram[i][j].first = diagram[i][j-1].first;
39
                    diagram[i][j].second = 2;
40
                }
41
            }
42
        }
43
    }
44
    
45
    // 回溯得到最长公共序列 
46
    string lcs;
47
    int a = lenx;
48
    int b = leny;
49
    while(a > 0 && b > 0){
50
        switch(diagram[a][b].second)
51
        {
52
        case 1:
53
            --a;
54
            break;
55
        case 2:
56
            --b;
57
            break;
58
        case 3:
59
            lcs+=x[a-1];
60
            --a;
61
            --b;
62
            break;
63
        }
64
    }
65

66
    // 打印结果
67
    for(int i = 1; i < lenx+1; ++i){
68
        for(int j = 1 ; j < leny+1 ; ++j){
69
            printf("{%d,%d} ",diagram[i][j].first,diagram[i][j].second);
70
        }
71
        printf("\n");
72
    } 
73
    cout << "lcs: ";
74
    std::reverse(lcs.begin(),lcs.end());
75
    for(const auto & it : lcs){
76
        cout << it;
77
    }
78
    cout << endl;
79
    
80
    cout << "size of lcs: " << diagram[lenx][leny].first << endl;
81

82
}
83

84

85

86
int main(){
87
    string X = "BFCDCVB";
88
    string Y = "EFBFUCB";
89

90

91
    Match(X,Y,X.length(),Y.length());
92
    return 0;
93
}
94

Python

xxxxxxxxxx
50
1
#!/usr/bin/env python
2
# coding=utf-8
3

4

5
def Match(x,y):
6
   
7
    # 初始化图
8
    diagram = []
9
    for i in range(len(x)+1):
10
        diagram.append([]); 
11
        for j in range(len(y)+1):
12
            diagram[i].append([0,0])
13
    
14
    # 填表
15
    for i in range(1,len(x)+1):
16
        for j in range(1,len(y)+1):
17
            if x[i-1] == y[j-1] :
18
                diagram[i][j][0] = diagram[i-1][j-1][0] + 1
19
                diagram[i][j][1] = 3
20
            else:
21
                if diagram[i-1][j] > diagram[i][j-1]:
22
                    diagram[i][j][0] = diagram[i-1][j][0]
23
                    diagram[i][j][1] = 1
24
                else:
25
                    diagram[i][j][0] = diagram[i][j-1][0]
26
                    diagram[i][j][1] = 2
27

28
    # 回溯
29
    LCS = []
30
    xi,yi = len(x),len(y) 
31
    while (xi > 0) and (yi > 0):
32
        if diagram[xi][yi][1] == 1:
33
            xi-=1
34
        elif diagram[xi][yi][1] == 2:
35
            yi-=1
36
        elif diagram[xi][yi][1] == 3:
37
            LCS.insert(0,x[xi-1])
38
            xi-=1
39
            yi-=1
40
    
41
    # 打印结果
42
    print("LCS: {}".format("".join(LCS)))
43
    print(f"length: {diagram[-1][-1][0]}")
44

45
if __name__ == "__main__":
46
    X = "BFCDCVB"
47
    Y = "EFBFUCB"
48

49
    Match(X,Y)
50

这么看来动态规划算法的有两个比较难的点：写出最优解结构的递归式和找到合适的回溯算法。