几道常见的动态规划题
通常暴力穷举的方式是一种糟糕的策略,动态规划正是一种解决类似问题的思想,如果一个问题满足最优子结构,就可以通过把原问题可以分解为几个子问题来解决,即全局的最优解一定是某个局部的最优解,我们需要一张表来保存前一次计算的结果,以便递推出原问题的解,这样可以避免重复计算。
0-1背包
固定容量的背包,希望能够装入价值最大的物品。假设 m ( i , j ) m(i,j) m(i,j)是背包容量为 j j j,可选择物品为0~i时0-1 背包问题的最优值。 w i w_i wi为第i个物品的容量, v i v_i vi为第i个物品的价值。
- 当i=0时,表示没有可选择的物品,总价值为0。
- 当j=0时,表示背包容量为0,不好装东西,总价值为0。
- 当 j ≥ w i j\geq w_i j≥wi时,可以选择放入物品i或不放入物品i,若不放入物品i,则只剩下i-1个物品可以选择,背包容量仍为j,若放入物品i,也只剩下i-1个物品可供选择,背包剩余容量为 j − w i j-w_i j−wi。
- 当 0 ≤ j < w i 0 \leq j < w_i 0≤j<wi时,背包放不下物品i,只能在剩余i-1个物品中选择。
递推式:
m
(
i
,
j
)
=
{
0
i
f
i
=
0
o
r
j
=
0
max
{
m
(
i
−
1
,
j
)
,
m
(
i
−
1
,
j
−
w
i
)
+
v
i
}
i
f
j
≥
w
i
m
(
i
−
1
,
j
)
i
f
0
≤
j
<
w
i
m(i,j)=\left\{\begin{matrix} 0 & if \ i=0 \ or \ j=0 \\ \max \{m(i-1, j), m(i-1, j-w_i)+v_i\} & if \ j \geq w_i\\ m(i-1,j) & if \ 0 \leq j < w_i \end{matrix}\right.
m(i,j)=⎩⎨⎧0max{m(i−1,j),m(i−1,j−wi)+vi}m(i−1,j)if i=0 or j=0if j≥wiif 0≤j<wi
代码:
import numpy as np
def knapsack(c, w, v):
assert len(w)==len(v)
m = np.full((len(w) + 1, c + 1), -1)
m[0,:] = 0 # i=0
m[:, 0] = 0 # j=0
for i in range(1, len(w) + 1):
for j in range(1, c + 1):
if w[i - 1] > j:
m[i, j] = m[i-1, j]
else:
m[i, j] = max(m[i-1, j], m[i-1, j-w[i - 1]] + v[i - 1])
print(m)
# solution
j = c
x = np.full(len(w), -1)
for i in range(len(w)):
if m[i+1, j] == m[i, j]:
x[i] = 0
else:
x[i] = 1
j = j - w[i]
return x
c = 10 # 背包容量
w = [2,2,6,5,4] # 物品容量
v = [6,3,5,4,6] # 物品价值
s = knapsack(c, w, v)
print(s)
[[ 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 6 6 6 6 6 6 6 6 6]
[ 0 0 6 6 9 9 9 9 9 9 9]
[ 0 0 6 6 9 9 9 9 11 11 14]
[ 0 0 6 6 9 9 9 10 11 13 14]
[ 0 0 6 6 9 9 12 12 15 15 15]]
[1 1 0 0 1]
可见,最大价值为15,选择第0个、第1个以及第4个物品。
最小编辑距离
给出两个单词source和target,求出删除、替换或插入某个字符使得source变为target的最少次数。
例如:source=‘intention’, target=‘execution’,最少次数为5,总共有5步。
- intention -> inention (删除 ‘t’)
- inention -> enention (将 ‘i’ 替换为 ‘e’)
- enention -> exention (将 ‘n’ 替换为 ‘x’)
- exention -> exection (将 ‘n’ 替换为 ‘c’)
- exection -> execution (插入 ‘u’)
假设 e d i t [ i , j ] edit[i,j] edit[i,j]为source的第0到i个字符与target的第0到j个字符分别组成的字符串的最小编辑距离。
- 当 i = 0 a n d j = 0 \ i=0 \ and \ j=0 i=0 and j=0 时,两个都为空串,距离为0。
- 当 i = 0 a n d j > 0 \ i=0 \ and \ j>0 i=0 and j>0 时,插入j次就行,距离为j。
- 当 j = 0 a n d i > 0 \ j=0 \ and \ i>0 j=0 and i>0 时,删除i次就行,距离为i。
- 当 i ≥ 1 a n d j ≥ 1 \ i\geq 1 \ and \ j\geq1 i≥1 and j≥1 时,看是否添加source的第i个字符或删除target的第j个字符能否使得剩余部分相同,即继续计算 e d i t [ i − 1 , j ] edit[i-1,j] edit[i−1,j]或 e d i t [ i , j − 1 ] edit[i,j-1] edit[i,j−1],如果添加或删除后还不等,则需要继续计算 e d i t [ i − 1 , j − 1 ] edit[i-1,j-1] edit[i−1,j−1],我们取最小的情况。
递推式:
e
d
i
t
[
i
,
j
]
=
{
0
i
f
i
=
0
a
n
d
j
=
0
j
i
f
i
=
0
a
n
d
j
>
0
i
i
f
j
=
0
a
n
d
i
>
0
min
{
e
d
i
t
[
i
−
1
,
j
]
+
i
n
s
c
o
s
t
(
s
o
u
r
c
e
i
)
,
e
d
i
t
[
i
−
1
,
j
−
1
]
+
s
u
b
s
t
c
o
s
t
(
s
o
u
r
c
e
i
,
t
a
r
g
e
t
j
)
,
e
d
i
t
[
i
,
j
−
1
]
+
d
e
l
c
o
s
t
(
t
a
r
g
e
t
j
)
)
}
i
f
i
≥
1
a
n
d
j
≥
1
edit[i,j]=\left\{\begin{matrix} 0 & if \ i=0 \ and \ j=0\\ j & if \ i=0 \ and \ j>0\\ i & if \ j=0 \ and \ i>0\\ \min\{edit[i-1,j]+inscost(source_i),\\edit[i-1,j-1]+substcost(source_i,target_j),\\edit[i,j-1]+delcost(target_j))\} & if \ i\geq 1 \ and \ j\geq1 \end{matrix}\right.
edit[i,j]=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧0jimin{edit[i−1,j]+inscost(sourcei),edit[i−1,j−1]+substcost(sourcei,targetj),edit[i,j−1]+delcost(targetj))}if i=0 and j=0if i=0 and j>0if j=0 and i>0if i≥1 and j≥1
其中,
i
n
s
c
o
s
t
(
s
o
u
r
c
e
i
)
inscost(source_i)
inscost(sourcei) 为插入的损失,一般为1,
d
e
l
c
o
s
t
(
t
a
r
g
e
t
j
)
delcost(target_j)
delcost(targetj)为删除的损失,一般为1,
s
u
b
s
t
c
o
s
t
(
s
o
u
r
c
e
i
,
t
a
r
g
e
t
j
)
substcost(source_i,target_j)
substcost(sourcei,targetj) 为替换的损失,可以为1或2。
代码:
import numpy as np
def min_edit_distance(target, source):
m = len(target)
n = len(source)
distance = np.full((m+1, n+1), np.Inf)
for i in range(m + 1):
distance[i, 0] = i
for j in range(n + 1):
distance[0, j] = j
for i in range(m):
for j in range(n):
if target[i]==source[j]:
distance[i+1, j+1] = distance[i,j]
else:
# ins-cost=1, subst-cost=2, del-cost=1
distance[i+1, j+1] = min(distance[i,j+1] + 1, distance[i,j] + 1, distance[i+1,j] + 1)
print(distance)
return distance[m, n]
target = 'intention'
source = 'execution'
print(min_edit_distance(target, source))
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[1. 1. 2. 3. 4. 5. 6. 6. 7. 8.]
[2. 2. 2. 3. 4. 5. 6. 7. 7. 7.]
[3. 3. 3. 3. 4. 5. 5. 6. 7. 8.]
[4. 3. 4. 3. 4. 5. 6. 6. 7. 8.]
[5. 4. 4. 4. 4. 5. 6. 7. 7. 7.]
[6. 5. 5. 5. 5. 5. 5. 6. 7. 8.]
[7. 6. 6. 6. 6. 6. 6. 5. 6. 7.]
[8. 7. 7. 7. 7. 7. 7. 6. 5. 6.]
[9. 8. 8. 8. 8. 8. 8. 7. 6. 5.]]
5.0
最长公共子序列
所谓子序列即在一段序列中删除任意元素后剩余的序列,子序列的元素可以不相邻,但要维持原本先后次序。例如ABCBDAB和BDCABA的最长公共子序列为BCBA,长度为4。
假设
X
=
<
x
1
,
x
2
,
.
.
,
x
m
>
X=<x_1, x_2, .., x_m>
X=<x1,x2,..,xm>和
Y
=
<
y
1
,
y
2
,
.
.
.
,
y
n
>
Y=<y_1,y_2,...,y_n>
Y=<y1,y2,...,yn>代表两个序列,令
Z
=
<
z
1
,
z
2
,
.
.
,
z
k
>
Z=<z_1,z_2,..,z_k>
Z=<z1,z2,..,zk>为
X
X
X和
Y
Y
Y的最长公共子序列。
- 如果 x m = y n x_m=y_n xm=yn,那么 z k = x m = y n z_k=x_m=y_n zk=xm=yn并且 Z k − 1 Z_{k-1} Zk−1是 X m − 1 X_{m-1} Xm−1与 Y n − 1 Y_{n-1} Yn−1的最长公共子序列。
- 如果 x m ≠ y n x_m\neq y_n xm̸=yn,那么 z k ≠ x m z_k\neq x_m zk̸=xm意味着 Z Z Z是 X m − 1 X_{m-1} Xm−1与 Y Y Y的最长公共子序列。
- 如果 x m ≠ y n x_m\neq y_n xm̸=yn,那么 z k ≠ y n z_k\neq y_n zk̸=yn意味着 Z Z Z是 X X X与 Y n − 1 Y_{n-1} Yn−1的最长公共子序列。
上面三条可用反证法证明,说明最长公共子序列问题是满足最优子结构的。
递推式:
c
[
i
,
j
]
=
{
0
i
f
i
=
0
o
r
j
=
0
c
[
i
−
1
,
j
−
1
]
+
1
i
f
i
,
j
>
0
a
n
d
x
i
=
y
i
max
(
c
[
i
,
j
−
1
]
,
c
[
i
−
1
,
j
]
)
i
f
i
,
j
>
0
a
n
d
x
i
≠
y
i
c[i,j]=\left\{\begin{matrix} 0 & if \ i=0 \ or \ j=0\\ c[i-1,j-1]+1 & if \ i,j>0 \ and \ x_i=y_i\\ \max(c[i,j-1],c[i-1,j]) & if \ i,j>0 \ and \ x_i\neq y_i \end{matrix}\right.
c[i,j]=⎩⎨⎧0c[i−1,j−1]+1max(c[i,j−1],c[i−1,j])if i=0 or j=0if i,j>0 and xi=yiif i,j>0 and xi̸=yi
其中,
c
[
i
,
j
]
c[i,j]
c[i,j]为
X
i
X_i
Xi和
Y
j
Y_j
Yj的最长公共子序列的长度。
代码:
import numpy as np
import sys
def lcs_length(x, y):
m = len(x)
n = len(y)
b = np.empty((m, n), dtype='str')
c = np.zeros((m+1, n+1))
for i in range(1, m+1):
for j in range(1, n+1):
if x[i-1]==y[j-1]: # 序列长度是从0开始
c[i,j] = c[i-1, j-1] + 1
b[i-1, j-1] = '↖'
elif c[i-1, j]>=c[i, j-1]:
c[i,j] = c[i-1, j]
b[i-1, j-1] = '↑'
else:
c[i,j] = c[i, j-1]
b[i-1, j-1] = '←'
return c[i, j], c, b
其中,表
c
[
0..
m
,
0..
n
]
c[0..m,0..n]
c[0..m,0..n]保存了最长公共子序列的长度,它是按行的顺序填表依次计算的。表
b
[
1..
m
,
1..
n
]
b[1..m,1..n]
b[1..m,1..n]则用来构造最优解,
b
[
i
,
j
]
b[i,j]
b[i,j]指向了构造最优子问题的选择方向,如下图所示。
然后,只要按着箭头指向的方向就可寻找到最优解了。
def print_lcs(b, x, len_x, len_y):
i = len_x
j = len_y
if i==0 or j==0:
return
if b[i-1, j-1]=='↖':
print_lcs(b, x, i-1, j-1)
sys.stdout.write(x[i-1])
elif b[i-1, j-1]=='↑':
print_lcs(b, x, i-1, j)
else:
print_lcs(b, x, i, j-1)
x = 'ABCBDAB'
y = 'BDCABA'
# x = '10010101'
# y = '010110110'
l, c, b=lcs_length(x, y)
print(l)
print(c)
print(b)
print_lcs(b, x, len(x), len(y))
递归打印如下:
4.0
[[0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1. 1. 1.]
[0. 1. 1. 1. 1. 2. 2.]
[0. 1. 1. 2. 2. 2. 2.]
[0. 1. 1. 2. 2. 3. 3.]
[0. 1. 2. 2. 2. 3. 3.]
[0. 1. 2. 2. 3. 3. 4.]
[0. 1. 2. 2. 3. 4. 4.]]
[['↑' '↑' '↑' '↖' '←' '↖']
['↖' '←' '←' '↑' '↖' '←']
['↑' '↑' '↖' '←' '↑' '↑']
['↖' '↑' '↑' '↑' '↖' '←']
['↑' '↖' '↑' '↑' '↑' '↑']
['↑' '↑' '↑' '↖' '↑' '↖']
['↖' '↑' '↑' '↑' '↖' '↑']]
BCBA
其实,该算法可以省去表 b b b,如果只需计算最终长度在空间上还可以继续优化,具体参见《算法导论》该小节后的习题。
矩阵链乘
给定若干个矩阵的序列 < A 1 , . . , A n > <A_1,..,A_n> <A1,..,An>,假设它们能够按顺序相乘,因为矩阵乘法满足结合律,不同加括号的方式可能会对乘法的性能有不同的影响。例如, < A 1 , A 2 , A 3 > <A_1,A_2,A_3> <A1,A2,A3>相乘,三个矩阵的规模分别为 10 × 100 10\times100 10×100、 100 × 5 100\times5 100×5、 5 × 50 5\times50 5×50,如果按照 ( ( A 1 A 2 ) A 3 ) ((A_1A_2)A_3) ((A1A2)A3)计算,计算 A 1 A 2 A_1A_2 A1A2需要做 10 × 100 × 5 = 5000 10\times100\times5=5000 10×100×5=5000次标量乘法,再与 A 3 A_3 A3相乘又需要做 10 × 5 × 50 = 2500 10\times5\times50=2500 10×5×50=2500次标量乘法,共7500次。如果按照 ( A 1 ( A 2 A 3 ) ) (A_1(A_2A_3)) (A1(A2A3))的顺序计算,计算 A 2 A 3 A_2A_3 A2A3需 100 × 5 × 50 = 25000 100\times5\times50=25000 100×5×50=25000次, A 1 A_1 A1再与之相乘需 10 × 100 × 50 = 50000 10\times100\times50=50000 10×100×50=50000次,共75000次。因此第一种比第二种快了10倍。
令 m [ i , j ] m[i,j] m[i,j]表示计算矩阵 A i . . j A_{i..j} Ai..j即 A i A i + 1 . . . A j A_iA_{i+1}...A_j AiAi+1...Aj所需标量乘法次数的最小值,假设 A i . . j A_{i..j} Ai..j的最优分割点 k k k在 A k A_k Ak和 A k + 1 A_{k+1} Ak+1之间,其中 i ≤ k < j i\leq k < j i≤k<j,我们用 s [ i , j ] s[i,j] s[i,j]来保存最优分割点 k k k。假定 A i A_i Ai的规模为 p i − 1 × p i ( i = 1 , 2 , . . , n ) p_{i-1}\times p_i (i=1,2,..,n) pi−1×pi(i=1,2,..,n)。
- 如果 i = j i=j i=j,说明矩阵链只包含唯一矩阵,不需做运算,代价为0。
- 如果 i < j i<j i<j, m [ i , j ] m[i,j] m[i,j]等于计算 A i . . k A_{i..k} Ai..k与 A k + 1.. j A_{k+1..j} Ak+1..j的代价加上两者相乘代价 p i − 1 p k p j p_{i-1}p_kp_j pi−1pkpj对于所有 k ( i ≤ k < j ) k (i \leq k<j) k(i≤k<j)的最小值。
递推式:
m
[
i
,
j
]
=
{
0
i
f
i
=
j
max
i
≤
k
<
j
{
m
[
i
,
k
]
+
m
[
k
+
1
,
j
]
+
p
i
−
1
p
k
p
j
}
i
f
i
<
j
m[i,j]=\left\{\begin{matrix} 0 & if \ i=j\\ \max\limits_{i\leq k<j}\{m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\} & if \ i<j \end{matrix}\right.
m[i,j]={0i≤k<jmax{m[i,k]+m[k+1,j]+pi−1pkpj}if i=jif i<j
代码:
import sys
import numpy as np
def matrix_chain_order(p):
n = len(p) - 1
m = np.mat(np.zeros((n, n)), dtype=np.int64)
s = np.mat(np.zeros((n, n)), dtype=np.int64)
for l in range(2, n + 1):
for i in range(0, n - l + 1):
j = i + l - 1 # 链长l=j-i+1
m[i, j] = sys.maxsize # 初始化为正无穷 2**63 - 1
for k in range(i, j):
# 注意:下标是从0开始的
q = m[i, k] + m[k + 1, j] + p[i] * p[k + 1] * p[j + 1]
if q < m[i, j]:
m[i, j] = q
s[i, j] = k
return m[0, n - 1], m, s
同样可以递归打印最优括号化方案:
def print_optimal_parens(s, i, j):
if i == j:
sys.stdout.write('A' + str(i))
else:
sys.stdout.write('(')
print_optimal_parens(s, i, s[i, j])
print_optimal_parens(s, s[i, j] + 1, j)
sys.stdout.write(')')
p = [5, 20, 50, 1, 100] # p是维度矩阵
v, m, s = matrix_chain_order(p)
print(v)
print(m)
print(s)
print_optimal_parens(s, 0, 3)
打印如下:
1600
[[ 0 5000 1100 1600]
[ 0 0 1000 3000]
[ 0 0 0 5000]
[ 0 0 0 0]]
[[0 0 0 2]
[0 0 1 2]
[0 0 0 2]
[0 0 0 0]]
((A0(A1A2))A3)
除此之外,我们还可以将其改为递归形式:
def matrix_chain_order2(p):
n = len(p) - 1
m = np.mat(np.zeros((n, n)), dtype=np.int64)
s = np.mat(np.zeros((n, n)), dtype=np.int64)
for i in range(0, n):
for j in range(i, n):
m[i, j] = sys.maxsize
return lookup_chain(m, s, p, 0, n - 1), m, s
def lookup_chain(m, s, p, i, j):
if m[i, j] < sys.maxsize:
return m[i, j]
if i == j:
m[i, j] = 0
else:
for k in range(i, j):
# 注意:下标是从0开始的
q = lookup_chain(m, s, p, i, k) + lookup_chain(m, s, p, k + 1, j) \
+ p[i] * p[k + 1] * p[j + 1]
if q < m[i, j]:
m[i, j] = q
s[i, j] = k
return m[i, j]