通常我们需要按照如下4个步骤设计一个动态规划的算法
- 刻画一个最优解的结构特征
- 递归地定义最优解的值
- 计算最优解的值，通常采用置底向上的方法
- 利用计算出的信息构造一个最有解
下面开始理解动态规划，从动态规划的钢条切割例子开始
1.钢条切割
| 长度i | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| 价格 p i p_{i} pi | 1 | 5 | 8 | 9 | 10 | 17 | 17 | 20 | 24 | 30 |
给定如上的的价格表，和一个长度为p的钢条，求出这个长度钢条以怎样的切割方案才可以获得最大利润？
下面给出切割钢条问题的递归解：
r
i
(
i
=
1
,
2
,
.
.
.
10
)
r_{i}(i=1,2,...10)
ri(i=1,2,...10)表示长度为i的最大利润。
更一般的，我们构造如下的递归方程，对于
r
n
(
n
≥
1
)
r_{n}(n\geq1)
rn(n≥1),我们可以用更短的钢条切割收益来描述它：
r
n
=
m
a
x
(
p
n
,
r
1
+
r
n
−
1
,
r
2
+
r
n
−
2
,
.
.
.
,
r
n
−
1
+
r
1
)
r_{n} = max(p_{n},r_{1}+r_{n-1},r_{2}+r_{n-2},...,r_{n-1}+r_{1})
rn=max(pn,r1+rn−1,r2+rn−2,...,rn−1+r1)
p
n
p_{n}
pn:表示对于长度为n的钢条不切割可以获得的利润
其他你n-1个参数对应的另外n-1种方案，例如：对于
r
i
+
r
i
−
1
r_{i}+r_{i-1}
ri+ri−1表示将钢条切割成i和i+1两段，接着求解这两段的最佳切割收益
r
i
r_{i}
ri和
r
n
−
i
r_{n-i}
rn−i.
我们发现求解
r
n
r_{n}
rn可以变成求解
r
i
r_{i}
ri和求解
r
i
+
1
r_{i+1}
ri+1这两个子问题的组合所以说钢条切割问题满足最优子结构*性质：
问题的最优解由相关子问题的最优解组合而成，而这些子问题可以独立求解.
下面我们介绍另一种比较简单的递归求解方法：
将钢条分成左边的钢条和右边的钢条，左边钢条不再切割，而对右边的钢条进行递归切割。这样不做任何切割的方案可以描述为：左边长度为n，收益为
p
n
p_{n}
pn,剩余部分长度为0，对应的收益
r
0
=
0
r_{0} = 0
r0=0。我们可以用如下的公式来表示：
r
n
=
m
a
x
l
≤
i
≤
n
(
p
i
+
r
n
−
i
)
r_{n} = \underset{l\leq i\leq n}{max}(p_{i}+r_{n-i})
rn=l≤i≤nmax(pi+rn−i)
# 自顶向下的递归实现
p = [1,5,8,9,10,17,17,20,24,30]
n = len(p)
def cutRod(p,n):
if n == 0:
return 0
q = -float('inf')
for i in range(0,n):
q = max(q,p[i]+cutRod(p,n-(i+1)))
return q
cutRod(p,10)
30
CUT—ROD函数会反复的计算相同的子问题，例如对于cutRod(4)来说，我们可以生成如下的递归树：
4
/|\ \
/ | | \
/ \ \ \
3 2 1 0
/|\ /\ \
/ | \ 1 0 0
2 1 0 |
/\ 0
1 0
|
0
当这个过程递归展开时，它所做的工作量会爆炸性的增长。
利用动态规划求解钢条切割问题。
如何利用动态规划优化钢条切割问题，即所有的子问题我们都只求解一遍，如果后来需要用到这个子问题的解，只需要查找保存的结果不需要再次计算。
动态规划有以下另种实现方法：
- 带备忘的自顶向下的方法：这个方法仍然按照自然的递归形式编写过程，这个过程会将求解得到的子问题的解保存在数组或者在哈希表中。如果需要到一个子问题的解是，只需要查找在数组或者哈希表中是否已经保存了这个子问题的解，如果有直接调用，如果没有，再递归的计算这个子问题。
- 自底向上的方法：这个方法一般需要定义子问题的规模。因为任何的子问题都需要计算更小的子问题，我们按照子问题的大小从小到大计算，当计算到某个子问题的时候，它所依赖的子问题都已经计算完成了。这样就可以保证，当我们第一次遇到这个子问题的时候，就是需要求解它的时候。
# 自顶向下的CUT-ROD过程的Python代码
def memoizedCutRod(p,n):
r = [-float('inf')]*n
return memoziedCutRodAux(p,n,r)
def memoziedCutRodAux(p,n,r):
if r[n-1] >= 0:
return r[n-1]
if n == 0:
q = 0
else:
q = -float('inf')
for i in range(n):
q = max(q,p[i]+memoziedCutRodAux(p,n-(i+1),r))
r[n-1] = q
return q
memoizedCutRod(p,10)
30
# 自底向上的方法
def bottomUpCutRod(p,n):
r = [0]*n
for j in range(n):
q = -float('inf')
for i in range(j+1):
q = max(q,p[i]+r[j-(i+1)])
r[j] = q
return r[-1]
bottomUpCutRod(p,10)
30
求斐波那契数列的第n项：
F
0
=
0
F
1
=
1
F
n
=
F
n
−
1
+
F
n
−
2
F_{0} = 0 \\ F_{1} = 1 \\ F_{n} = F_{n-1}+F_{n-2}
F0=0F1=1Fn=Fn−1+Fn−2
def focb(n):
res = [0]*(n+1)
res[1] = 1
for i in range(2,n+1):
res[i] = res[i-1]+res[i-2]
return res[-1]
focb(4)
3
2.矩阵链乘法
首先给出一个矩阵相乘的方法。
属性rows和columns是矩阵的行数和列数。
import numpy as np
def matrixMut(A,B):
rowsA = len(A)
colsA = len(A[0])
rowsB = len(B)
colsB = len(B[0])
if colsA != rowsB:
raise Exception("incompatible dimensions")
else:
C = [[0 for _ in range(colsB)] for _ in range(rowsA)]
for i in range(rowsA):
for j in range(colsB):
for k in range(colsA):
C[i][j] += A[i][k] * B[k][j]
return C
A = np.random.randint(0,10,(4,3))
B = np.random.randint(0,10,(3,4))
C = matrixMut(A,B)
C1 = np.dot(np.mat(A),np.mat(B))
C
[[24, 84, 84, 115], [46, 56, 96, 80], [43, 112, 127, 152], [37, 84, 101, 119]]
C1
matrix([[ 24, 84, 84, 115],
[ 46, 56, 96, 80],
[ 43, 112, 127, 152],
[ 37, 84, 101, 119]])
上面是矩阵相乘的程序，并且用numpy的矩阵相乘函数np.dot()检验。
我们发现对于一个矩阵相乘例如：
<
A
1
,
A
2
>
<A_{1}, A_{2}>
<A1,A2>,他们的规模分别为10×100、 100*5.那么我们要做的标量乘法就是
10
⋅
100
⋅
5
10\cdot 100 \cdot 5
10⋅100⋅5次。
对于一个矩阵链乘，会有不同的括号化方案，对应的不同的标量乘法次数。我们要找到一种括号化方案，似的矩阵的链乘积的标量乘法次数最少。
利用暴力法分析问题规模。
如果利用暴力法，那么我们就需要计算所有的括号化方案，来找到一个标量乘法次数最少的一个括号化方案，那么用P(n)来表示，那矩阵链的长度为n时，的
括号化的方案数。那么我们就有如下的递归方程。
P
(
n
)
=
{
1
i
f
  
n
=
=
1
∑
k
=
1
n
−
1
P
(
k
)
P
(
n
−
k
)
i
f
  
n
≥
2
P(n) = \left\{\begin{matrix} 1 & if\; n==1 \\ \sum ^{n-1}_{k=1}P(k)P(n-k) & if\; n\geq 2 \end{matrix}\right.
P(n)={1∑k=1n−1P(k)P(n−k)ifn==1ifn≥2
这个递归式的结果为
ω
(
2
n
)
\omega (2^n)
ω(2n),这是一种糟糕的策略。
应用动态规划的方法
利用我们一开头就提到的应用动态规划的步骤。
- 刻画一个最优的子结构
- 构造递归式
- 计算最有解通常使用自底向上的方法
- (如必要)根据计算的信息构造一个最优解。
第一步那就是找到一个最优的子结构了。
对于
A
i
A
i
+
1
.
.
.
A
j
A_{i}A_{i+1}...A{j}
AiAi+1...Aj,对其进行括号化，加入我们在矩阵
A
k
A_{k}
Ak的后边添加一个")"那么我们下一步就要计算前边半条矩阵链相乘得到矩阵
A
i
,
k
A_{i,k}
Ai,k的代价，加上后边矩阵链相乘得到矩阵
A
k
+
1
,
j
A_{k+1,j}
Ak+1,j的代价，然后加上这两个矩阵相乘的代价，这就是全部的代价。
那么最优子结构就是对于
A
i
A
i
+
1
.
.
.
A
j
A_{i}A_{i+1}...A{j}
AiAi+1...Aj，找到他的最优去括号化方案，将其分成两段，其中前段是
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak,那么单独求这半个的链的最优括号化方案，如果将这个方案放在整个链中也是最优的括号化方案。这个原因可以用剪切–粘贴技术证明：
假如我们有一个矩阵链
A
i
A
i
+
1
.
.
.
A
j
A_{i}A_{i+1}...A{j}
AiAi+1...Aj：它的最优括号化方案的分割点在k
方案一： 将其分成两段，
A
i
A
i
+
1
.
.
.
A
k
,
A
k
+
1
A
k
+
2
.
.
.
A
j
A_{i}A_{i+1}...A{k},A_{k+1}A_{k+2}...A{j}
AiAi+1...Ak,Ak+1Ak+2...Aj,我们分别单独求其的最优括号化方案，得到了:
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak,
A
k
+
1
A
k
+
2
.
.
.
A
j
A_{k+1}A_{k+2}...A{j}
Ak+1Ak+2...Aj (这里我们用不同的颜色来表示最优的括号化方案。)
将这两段拼接起来就是整体链的最佳括号化方案。
方案二：如果我们直接求这个整体链的最佳括号化方法，得到如下的最佳方案。
A i A i + 1 . . . A k A_{i}A_{i+1}...A{k} AiAi+1...Ak A k + 1 A k + 2 . . . A j A_{k+1}A_{k+2}...A{j} Ak+1Ak+2...Aj
接下来我们将这个方案二的紫色的前一段
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak提取出来，现在我们假设这一段不是
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak的最优解，也就是说
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak这个方案，比方案
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak,代价更高。那么我们将
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak这一段添加到方案二中得到：
A
i
A
i
+
1
.
.
.
A
k
A_{i}A_{i+1}...A{k}
AiAi+1...Ak
A
k
+
1
A
k
+
2
.
.
.
A
j
A_{k+1}A_{k+2}...A{j}
Ak+1Ak+2...Aj
显然这个是更优的解，所以就证明了方案二是不成立的。
构造递归式
用m[i,j]来表示计算矩阵
A
i
.
.
j
A_{i..j}
Ai..j所需要的最小标量乘法。
对于i==j的情况只有一个矩阵必然有m[i,j]=0,而对于i<j的情况，假设最佳分割点在矩阵
A
k
A_{k}
Ak和
A
k
+
1
A_{k+1}
Ak+1之间，那么我们就可以写出如下的递归式：
m
[
i
,
j
]
=
m
[
i
,
k
]
+
m
[
k
+
1
,
j
]
+
p
i
−
1
p
k
p
j
m[i,j] = m[i,k] + m[k+1,j] + p_{i-1}p{k}p{j}
m[i,j]=m[i,k]+m[k+1,j]+pi−1pkpj
而实际中我们不知道k是多少，这里k的取值一共有j-i种可能，即
k
=
i
,
i
+
1
,
.
.
.
,
j
k=i,i+1,...,j
k=i,i+1,...,j,我们遍历这些情况找到最佳的k值。所以递归式为：
m
[
i
,
j
]
=
{
0
i
f
  
i
=
=
j
m
i
n
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\\ \underset {i\leq k < j} {min} \{ m[i,k], m[k+1,j] + p_{i-1}p_{k}p_{j}\} & if\; i < j \end{matrix}\right.
m[i,j]={0i≤k<jmin{m[i,k],m[k+1,j]+pi−1pkpj}ifi==jifi<j
计算最优代价
假定矩阵
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)，他的输入是一个维度序列，
<
p
0
,
p
1
,
.
.
.
p
n
>
<p_{0},p_{1},...p_{n}>
<p0,p1,...pn>,这个过程用一个辅助表m[i,j]保存计算代价，而用s[1,…n-1,2…n]记录切割位置，以此来构造最优解。
为了采用自底向上的方法，我们必须确定计算m[i,j],会用到那些子问题，j-i+1个链相乘只依赖于那些比它短的链相乘的代价。
例如计算
A
1
,
3
A_{1,3}
A1,3的代价，他的长度为3，我们会用到链中所有长度小于3的链乘积代价，那么计算
A
1
,
4
A_{1,4}
A1,4,我们就会用到所有长度小于4的链的乘积代价，我们可以发现，我们可能会计算相同的子问题，这种子问题的重叠性也是动态规划的另一个标识。
那么我们采用自底向上的方法，我先解决规模最小的子问题，当我们需要解规模较大的子问题时，它所以依赖的较小子问题我们都已经解出来了，直接调用即可。所以我们就应该按照链的长度的递增顺序解决这个问题。
# 矩阵链乘积
def matrixChainOrder(p):
n = len(p)-1
m = [[0]*(n) for _ in range(n)]
s = [[0]*(n-1) for _ in range(n-1)]
for l in range(1,n):
for i in range(n-l):
j = i + l
m[i][j] = float('inf')
for k in range(i,j):
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-1] = k+1
return m,s
p = [30,35,15,5,10,20,25]
m,s = matrixChainOrder(p)
m,s
([[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]],
[[1, 1, 3, 3, 3],
[0, 2, 3, 3, 3],
[0, 0, 3, 3, 3],
[0, 0, 0, 4, 5],
[0, 0, 0, 0, 5]])
最后构造最优解
目前我们可以求出一个矩阵链乘的最小代价，但是我们还不能显示出最优的括号化方案，求得的s,s[i,j]保存的就是从矩阵 A i A_{i} Ai到矩阵 A j A_{j} Aj加入“)”的位置。我们可以递归的加入"()"
def printOptParens(s,i,j):
if i==j:
print("A{}".format(str(i+1)),end="")
else:
print("(",end="")
printOptParens(s,i,s[i][j-1]-1)
printOptParens(s,s[i][j-1],j)
print(")",end="")
printOptParens(s,0,5)
((A1(A2A3))((A4A5)A6))
1,对矩阵规模序列<5,10,3,12,5,50,6>,求矩阵链最优括号化方案
p = [5,10,3,12,5,50,6]
m,c = matrixChainOrder(p)
m,c
([[0, 150, 330, 405, 1655, 2010],
[0, 0, 360, 330, 2430, 1950],
[0, 0, 0, 180, 930, 1770],
[0, 0, 0, 0, 3000, 1860],
[0, 0, 0, 0, 0, 1500],
[0, 0, 0, 0, 0, 0]],
[[1, 2, 2, 4, 2],
[0, 2, 2, 2, 2],
[0, 0, 3, 4, 4],
[0, 0, 0, 4, 4],
[0, 0, 0, 0, 5]])
printOptParens(c,0,5)
((A1A2)((A3A4)(A5A6)))
- 设计递归算法MATRIX-CHAIN-MULTIPLY(A,s,i,j),实现矩阵链乘的算法。
def matrixChainMut(A,s,i,j):
if i==j:
return A[i]
if i==j+1:
return matrixMut(A[i],A[j])
else:
B1 = matrixChainMut(A,s,i,s[i][j-1]-1)
B2 = matrixChainMut(A,s,s[i][j-1],j)
return matrixMut(B1,B2)
# 根据维度链生成一系列随机矩阵
def randomMatrixs(p):
matrixs = list()
for i in range(len(p)-1):
matrixs.append(np.random.randint(0,10,(p[i],p[i+1])))
return matrixs
matrixs = randomMatrixs(p)
matrixChainMut(matrixs,s,0,5)
[[685009380, 681729124, 741674364, 597968918, 622221146, 737326806],
[574114480, 571367532, 621146920, 500886252, 521416180, 617575464],
[741957590, 738384282, 803052608, 647465550, 673855082, 798375320],
[920582443, 916179757, 996816790, 803671532, 836228448, 990964623],
[959841592, 955230604, 1038517360, 837402704, 871703472, 1032528964]]
694
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
