算法设计技巧: 动态规划 (Dynamic Programming)

动态规划(Dynamic Programming)是一种求解优化问题的技巧. 它适用于这样一类优化问题: 原问题包含了一系列 重叠的 子问题. 换句话说, 父问题最优解的结构包含了其子问题的最优解. 例如斐波那契(Fibonacci)数列, 已知$F(n-1)$和$F(n-2)$的值, 我们可以在$O(1)$内计算$F(n)=F(n-1)+F(n-2)$. 因此$F(n)$的值可以通过计算所有子问题$F(i)$, $i=1,2,\dots ,n-1$得到.

基本步骤

给定一个优化问题, 我们可以按照如下思路设计动态规划算法.

1. 刻画最优解的结构. 首先定义子问题, 然后用反证法证明原问题的最优解包含子问题的最优解. 注意: 不是所有问题都可以用动态规划求解. 因此, 设计动态规划算法之前必须验证最优解结构的正确性.

2. 定义最优值的递归式. 根据最优解的结构定义 目标函数 的递归表达式. 注意: 函数初始值的定义.

3. 计算最优值并构造最优解. 计算递归表达式一般有两种方式: (a) 递归求解并记录所有子问题的计算结果(Top-down with memoization), 递归地利用子问题的解得到原问题的解; (b) 按问题规模从小到大求解(Bottom-up). 构造最优解的方式可以通过记录递归式的求解信息然后构造得到(详情可以参考下面的例子). (一般来说, 已知递归式构造最优解是比较直观的.)

最短路

以最短路问题为例(下图). 给定图$G=(V,E)$, $V$代表顶点的集合, $E$代表边的集合. $\forall (i,j)\in E$, 令$c_{i,j}\in {\mathbb{R}}^{+}$代表$(i,j)$的距离. 求任意$s,t\in V$两点之间的最短路.

给定任意两点$s,t$, 最短路问题可以用Dijkstra算法在$O(m+n\log n)$求解, 其中$m=|E|,n=|V|$. 如果对所有的$(s,t)\in V\times V$, 用Dijkstra算法求解每一个最短路问题, 总的时间复杂度问题$O(mn+n^2\log n)=O(n^5)$. 下面我们用动态规划的算法在$O(n^4)$内求解.

定义子问题

令 V k = { 1 , 2 , … , k } V_k = \{1, 2, \ldots, k\} Vk​={1,2,…,k}. 令 P k ( u , v ) = { u , i 1 , … , i s , v } P_k(u,v) = \{u, i_1, \ldots, i_s, v\} Pk​(u,v)={u,i1​,…,is​,v}代表$V_k$上的路(Path), 其中$i_1,\dots ,i_s\in V_k$, 称为内部点(Interior points). 设$P_k^{\ast }(u,v)$代表$V_k$上的最短路, 因此原问题是要计算$P_n^{\ast }(u,v),\forall (u,v)\in V\times V$.

最优解的结构

考虑一条最短路 P = { v 1 , . . . , v k } P = \{v_1, ..., v_k\} P={v1​,...,vk​}以及它的任意一条子路 P ′ = { v i , … , v j } P' = \{v_i, \ldots, v_j\} P′={vi​,…,vj​}, 那么$P^{\prime }$是从$v_i$到$v_j$的最短路.

证明. 假设存在一条比$P^{\prime }$更短的从$u$到$v$的路$P^{\prime \prime }$, 那么我们可以构造一条更短的从$v_1$到$v_k$的路: 沿着$P$从$v_1$出发, 然后沿着$P^{\prime \prime }$从$v_i$到$v_j$, 最后沿着$P$到达$v_k$(如下图所示). 与$P$是最短路的假设矛盾.

递归式

令$d_k(u,v)$代表$V_k$上从$u$到$v$的最短路距离. 为方便描述, 令$P$代表对应的最短路. 考虑两种情况:

1. $k$是$P$的内部点. 把$P$在$k$处截断成条子路$P_1,P_2$, 根据前文 最优解的结构, 我们知道$P_1$和$P_2$都是最短路. 我们有$d_k(u,v)=d_{k-1}(u,k)+d_{k-1}(k,v)$.
2. $k$不是$P$的内部点. 换句话说$k$不在最优解中, 因此$d_k(u,v)=d_{k-1}(u,v)$.

综上所述, 我们有
d k ( u , v ) = min ⁡ { d k − 1 ( u , v ) + d k − 1 ( k , v ) , d k − 1 ( u , v ) } , k ≥ 1. d_k(u,v) = \min \{d_{k-1}(u,v) + d_{k-1}(k,v), d_{k-1}(u,v)\}, \quad k\geq 1. dk​(u,v)=min{dk−1​(u,v)+dk−1​(k,v),dk−1​(u,v)},k≥1.

初始条件

$$d_0(u,v)=\{\begin{array}{ll}{c}_{u,v},& (u,v)\in E\\ 0,& u=v\\ \infty ,& u\ne v,(u,v)\notin E.\end{array}$$

时间复杂度

从两个角度来看: 1.宽度, 即根据子问题构造父问题的时间. 计算 min ⁡ { ⋅ } \min\{\cdot\} min{⋅}的时间为$O(1)$, 一共要计算$V\times V$个顶点对, 因此构造所有$d_k(u,v)$的时间是$O(n^2)$; 2. 深度, 即递归的深度. $k=1,2,\dots ,n$, 因此深度为$(n)$. 综上述所, 总的时间复杂度为$O(n^2\cdot n)=O(n^3)$.

求解

根据上述递归式, 我们可以采用Top-down或Bottom-up的方式求解. (个人比较喜欢Bottom-up的方式, 原因是代码比较干净, 容易分析).

# Floyd-Warshall Algorithm
d(u,v) = d0(u,v), for all u,v in V  # 初始化d0(u,v)
for k = 1 to n:
	for (u,v) in V * V:
		if d(u,v) > d(u,k) + d(k,v):
			d(u,v) = d(u,k) + d(k, v)

构造最优解

令$P_{u,v}\subseteq V$代表$u$到$v$的最短路. 回顾递归式 d k ( u , v ) = min ⁡ { d k − 1 ( u , v ) + d k − 1 ( k , v ) , d k − 1 ( u , v ) } d_k(u,v) = \min \{d_{k-1}(u,v) + d_{k-1}(k,v), d_{k-1}(u,v)\} dk​(u,v)=min{dk−1​(u,v)+dk−1​(k,v),dk−1​(u,v)}, 当$d_k(u,v)=d_{k-1}(u,v)+d_{k-1}(k,v)$时, 说明$k\in P_{u,v}$. 因此$P(u,v)=P(u,k)\cup P(k,u)$.

# Floyd-Warshall Algorithm
d(u,v) = d0(u,v), for all u,v in V  # 初始化d0(u,v)
P(u,v) = empty for all u,v in V  # 初始化u,v之间的最短路的内部点
for k = 1 to n:
	for (u,v) in V * V:
		if d(u,v) > d(u,k) + d(k,v):
			d(u,v) = d(u,k) + d(k, v)
			# 根据递归式构造原问题的解
			P(u,v) = P(u,k) + {k} +  P(k,v) 

总结

1. 构造递归表达式是动态规划的核心. 递归式定义好之后, 算法实现和最优解的构造都是比较直观的.
2. 定义递归式时需要考虑最优解的结构. 一般来说, 我们可以通过标准的反证法来证明结构的正确性. 特别需要注意的是: 不是所有的问题都可以用动态规划的方式求解. 以 最长路 为例, 上述最优解结构的证明还凑效吗? (留给读者思考)

Python实现

from copy import deepcopy


def shortest_paths(c):
    """ 计算Graph中任意两点的最短路
    :param c: 成本矩阵
        * c[i][j] = infinity if (i,j) not in E
        * c[i][j] = 0 if i = j
        * c[i][j] = cost from i to j, if (i,j) in E
    :return: i,j之间的最短路程以及最短路
    """
    n = len(c)
    d = deepcopy(c)  # 初始化i和j之间的最短路程
    paths = init_paths(n)  # 初始化i和j之间的最短路
    for k in range(n):
        for i in range(n):
            for j in range(n):
                if d[i][j] > d[i][k] + d[k][j]:
                    d[i][j] = d[i][k] + d[k][j]
                    # 记录最短路
                    paths[i][j] = paths[i][k] + [k] + paths[k][j]
    return d, paths


def init_paths(n):
    """ 初始化任意两点的最短路为空list.

    :param n: 顶点的个数
    :return: n*n的list
    """
    paths = [[]] * n
    for i in range(n):
        paths[i] = [[]] * n
    return paths

完整代码

最长公共子序列

在计算生物学中DNA可以表示成含4种氮碱基 { A , C , G , T } \{A, C, G, T\} {A,C,G,T}构成的序列. 考虑两个序列:
$$\begin{array}{rl}& X=ACCGGTCGAGTGCGCGGAAGCCGGCCGAA\\ & Y=GTCGTTCGGAATGCCGTTGCTCTGTAAA\end{array}$$
我们希望计算$X$和$Y$之间的相似性. 下面定义相似性: 给定两个字符串$s,s^{\prime }$, 如果$s$包含$s^{\prime }$的所有字符且保持相同的前后顺序, 我们说$s^{\prime }$是$s$的子序列. 例如acdf是abcdeef的子序列. 我们用$X$和$Y$中最长的公共子序列的长度来衡量$X$和$Y$的相似性.

问题描述

给定两个字符串$X,Y$. 计算$X,Y$的最长公共子序列$Z$.

定义子问题

• 令$X_i=⟨x_1,x_2,\dots ,x_i⟩$, $Y_j=⟨y_1,y_2,\dots ,y_j⟩$
• 假设$|X|=m,|Y|=n$, 目标是计算$X_m$和$Y_n$的最长公共子序列

递归式

令$c_{i,j}$代表$X_i$和$Y_j$最长公共子序列的长度. 考虑两种情况:

• $x_i=y_j$. 那么$x_i$(即$y_j$)一定在公共子序列中. 因此, 通过计算$X_{i-1}$和$Y_{j-1}$的最长公共子序列, 即可得到$X_i$和$Y_j$的最长公共子序列. 我们有$c_{i,j}=c_{i-1,j-1}+1$.
• $x_i\ne y_j$. 那么$x_i$或$y_j$不在最长的公共子序列中. 因此 c i , j = max ⁡ { c i − 1 , j , c i , j − 1 } c_{i,j} = \max\{c_{i-1,j}, c_{i,j-1}\} ci,j​=max{ci−1,j​,ci,j−1​}

综上所述, 我们有
c i , j = { c i − 1 , j − 1 + 1 i , j > 0 , x i = y j max ⁡ { c i , j − 1 , c i − 1 , j } i , j > 0 , x i ≠ x j 0 i = 0  or  j = 0 c_{i,j} = \begin{cases} c_{i-1,j-1} + 1 & i,j>0, x_i = y_j \\ \max \{c_{i,j-1}, c_{i-1,j}\} & i, j > 0, x_i\neq x_j\\ 0 & i=0 \text{ or } j = 0 \end{cases} ci,j​=⎩⎪⎨⎪⎧​ci−1,j−1​+1max{ci,j−1​,ci−1,j​}0​i,j>0,xi​=yj​i,j>0,xi​​=xj​i=0 or j=0​

Remark: 最优解的结构和时间复杂度分析留给读者.

构造最优解

在计算递归式的同时记录下标位置的变化, 用来追溯公共子序列. 回顾上面的递归式, 考虑三种情况:

1. $c_{i,j}=c_{i-1,j-1}+1$. 令$b_{i,j}=‘‘\nwarrow "$, 表示从$(i,j)$向左上移动到$(i-1,j-1)$;
2. $c_{i,j}=c_{i,j-1}$. 令$b_{i,j}=‘‘\uparrow "$, 表示从$(i,j)$向上移动到$(i,j-1)$;
3. $c_{i,j}=c_{i-1,j}$. 令$b_{i,j}=‘‘\leftarrow$, 表示从$(i,j)$向左移动到$(i-1,j)$.

从$(i,j)=(m,n)$开始, 按照$b_{i,j}$指示的方向判断$x_i,y_j$是否相同, 如果相同则把$x_i$加入到公共子序列中(如下图所示).

图片来自经典教材算法导论第三版(Introduction to Algorithms, Third Edition).

Python实现

import numpy as np


def longest_common_subsequence(x, y):
    """ 根据递归式计算最长公共子序列
    """
    m, n = len(x), len(y)
    # 把x,y中字符的下标向右移一位(方便计算)
    x1 = ' ' + x
    y1 = ' ' + y
    # c[m][n]为最大公共子序列的长度
    c = np.zeros((m+1, n+1))  # 初始化c
    # 记录下标改变的路径(用字符串表示)
    b = np.zeros((m+1, n+1)).tolist()
    for i in range(1, m+1):
        for j in range(1, n+1):
            if x1[i] == y1[j]:
                c[i][j] = c[i-1][j-1] + 1
                b[i][j] = 'lu'  # go left up
            elif c[i-1][j] >= c[i][j-1]:
                c[i][j] = c[i-1][j]
                b[i][j] = 'u'  # go up
            else:
                c[i][j] = c[i][j-1]
                b[i][j] = 'l'  # go left
    return get_common_subsequence(x1, y1, b)


def get_common_subsequence(x1, y1, b):
    """ 根据b还原最长公共子序列
    """
    i, j = len(x1) - 1, len(y1) - 1
    res = []
    while i and j:
        if x1[i] == y1[j]:
            res.insert(0, x1[i])
        if b[i][j] == 'lu':  # go left up
            i -= 1
            j -= 1
        elif b[i][j] == 'l':  # go left
            j -= 1
        elif b[i][j] == 'u':  # go up
            i -= 1
    return ''.join(res)

完整代码