算法导论(4)动态规划
动态规划通常用来解决最优化问题。与分治方法相似,都是通过组合子问题的解来求解原问题。不同之处在于,分治方法将问题划分为互不相交的子问题,递归地求解子问题,此时如果子问题出现重叠的话就会反复求解,而动态规划对每个子问题只求解一次,将其解保存再一个表格中,从而不用重复计算。
我们通常按如下4步来设计一个动态规划算法:
- 刻画一个最优解的结构特征
- 递归地定义一个最优解的值
- 计算最优解的值,通常采用自底向上的方法
- 利用计算出的信息构造一个最优解
接下来通过几个示例来说明怎么用动态规划来设计算法。
1.钢条切割
问题描述
钢条切割:给定一条长度为n的钢条和一个价格表,求钢条的最佳切割方案(使销售收益 r n r_n rn最大)。价格表可以如下:
长度 i i 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 |
n取1到10的最佳方案我们可以直接看出:
r1=1 | r2=5 | r3=8 | r4=10 | r5=13 | r6=17 | r7=18 | r8=22 | r9=25 | r10=30 |
---|---|---|---|---|---|---|---|---|---|
1=1 | 2=2 | 3=3 | 4=2+2 | 5=2+3 | 6=6 | 7=1+6//2+2+3 | 8=2+6 | 9=3+6 | 10=10 |
我们考虑一下容易发现一个结论:
r
n
=
max
(
p
n
,
r
1
+
r
n
−
1
,
.
.
.
,
r
n
−
1
+
r
1
)
r_n=\max(p_n,r_1+r_{n-1},...,r_{n-1}+r_1)
rn=max(pn,r1+rn−1,...,rn−1+r1)
简写为:
r
n
=
max
1
≤
i
≤
n
(
p
i
+
r
n
−
i
)
r_n=\max_{1\le i\le n}(p_i+r_{n-i})
rn=1≤i≤nmax(pi+rn−i)
由此我们可以采取一种自顶向下的递归方法来求解这个问题:
CUT-ROD(p,n)
if n==0
return 0
q=-infty
for i=1 to n
q=max(q,p[i]+CUT-ROT(p,n-i))
return q
这个算法其实有点问题的,有些细节上的问题,这是书上的版本,其中一个明显的错误是,数组p只有10个元素,如果n大于10的话就会越界,所以上面的算法理应改为
if n==0
return 0
q=-infty
for i=1 to n
if i<=10
q=max(q,p[i]+CUT-ROT(p,n-i))
else if i!=n
q=max(q,CUT-ROD(p,i)+CUT-ROD(p,n-i))
return q
(这是我自己的版本,可能不符合伪代码的语言规则)
这个算法思路简单,但由于它通过递归调用会反复求解相同的子问题,所以随着n的增加工作量会爆炸性的增加。
根据上面的算法可以容易的给出C++版本的算法实现
#include<iostream>
#include<vector>
#include<climits>
using namespace std;
int cut_rod(vector<int>p, int n) {
if (n == 0)
return 0;
int q = -INT_MAX;
vector<int>vq;
vq.push_back(q);
for (int i = 0; i < n; ++i) {
if(i<10)
vq.push_back(p[i] + cut_rod(p, n - i - 1));
else if(i!=n-1)
vq.push_back(cut_rod(p,i+1) + cut_rod(p, n - i - 1));
}
int max = -INT_MAX;
for (auto c : vq) {
max = c > max ? c : max;
}
return max;
}
int main() {
vector<int>p = { 1,5,8,9,10,17,17,20,24,30 };
int m = cut_rod(p, 12);
cout << m << endl;
}
分析可知CUT-ROD的运行时间为 Θ ( 2 n ) \Theta(2^n) Θ(2n),为n的指数函数。对于长度为n的钢条,该算法需要考察 2 n − 1 2^{n-1} 2n−1种切割方案。。
用动态规划方法求解最优钢条切割问题
动态规划会将每个子问题求解后将结果保存下来,它是消耗额外的内存空间来节省运行时间,它时间上的节省是非常大的,可能将一个指数时间的解转化为一个多项式时间的解。动态规划有两种等效的实现方法:带备忘的自定向下法和自底向上法。
下面是加入备忘机制的自顶向下CUT-ROD过程
MEMOIZED-CUT-ROD(p,n)
let r[0,n]be a new array
for i=0 to n
r[i]=-infty
return MEMOIZED-CUT-ROD-AUX(p,n,r)
MEMOIZED-CUT-ROD-AUX(p,n,r)
if r[n]>=0
return r[n]
if n==0
q==0
else q=-infty
for i=1 to n
q=max(q,p[i]+MEMOIZED-CUT-ROD-AUX(p,n-i,r))
r[n]=q
return q
C++版本实现上面的算法,只需要在之前的代码的基础上稍微修改一下就行了
#include<iostream>
#include<vector>
#include<climits>
using namespace std;
int cut_rod(vector<int>p, int n, vector<int>& r);
//必须先声明,否则会导致找不到标识符错误
int m_cut(vector<int>p,int n) {
vector<int>r;
for (int i = 0; i <= n; ++i)
r.push_back(-INT_MAX);
return cut_rod(p, n, r);
}
int cut_rod(vector<int>p, int n,vector<int>& r) {
if (r[n] >= 0)
return r[n];
int q;
if (n == 0)
q=0;
else {
q = -INT_MAX;
vector<int>vq;
vq.push_back(q);
for (int i = 0; i < n; ++i) {
if (i < 10)
vq.push_back(p[i] + cut_rod(p, n - i - 1,r));
else if (i != n - 1)
vq.push_back(cut_rod(p, i + 1,r) + cut_rod(p, n - i - 1,r));
}
for (auto c : vq) {
q = c > q ? c : q;
}
}
r[n] = q;
return q;
}
int main() {
vector<int>p = { 1,5,8,9,10,17,17,20,24,30 };
int m = m_cut(p, 17);
cout << m << endl;
}
而自顶向上的版本就比较简单:
BOTTOM-UO-CUT-ROD(p,n)
let r[0..n] be a new array
r[0]=0
for j =1 to n
q=-infty
for i=1 to j
q=max(q,p[i]+r[j-i])
r[j]=q
return r[n]
C++版本如下,本身的算法很简单,注意自顶向上版本不是递归的,它本身思考问题的角度与前面是不一样的,正好反过来思考。这个算法是简单有效的。
int bucr(vector<int>p, int n) {
vector<int>r = { 0 };
for (int j = 1; j <= n; ++j) {
int q = -INT_MAX;
for (int i = 1; i <= j; ++i) {
if (i <= 10) {
q=(p[i - 1] + r[j - i])>q? (p[i - 1] + r[j - i]) :q;
}
else if (i != j) {
q=(r[i] + r[j - i])>q? (r[i] + r[j - i]):q;
}
}
r.push_back(q);
}
return r[n];
}
int main() {
vector<int>p = { 1,5,8,9,10,17,17,20,24,30 };
//int m = m_cut(p, 25);
int n = bucr(p, 25);
cout << "---"<<n<<endl;
}
另外我们这里要求最佳方案,我们现在仅仅求出了最佳方案下的可出售的金额,我们需要稍微修改 一下输出最佳最佳切割方案。
所以下面是改进版,会返回最佳方案
#include<iostream>
#include<vector>
#include<climits>
using namespace std;
int bucr(vector<int>p, int n, vector<int>& s) {
vector<int>r = { 0 };
vector<vector<int>>ss = { {0} };
for (int j = 1; j <= n; ++j) {
int q = -INT_MAX;
for (int i = 1; i <= j; ++i) {
if (i <= 10) {
if (p[i - 1] + r[j - i] > q){
q = p[i - 1] + r[j - i];
s.clear();
s.push_back(i);
for (auto c : ss[j - i])
if(c!=0)
s.push_back(c);
}
}
else if (i != j) {
if ((r[i] + r[j - i]) > q) {
q = (r[i] + r[j - i]);
s.clear();
for (auto c : ss[i])
if(c!=0)
s.push_back(c);
for (auto c : ss[j - i])
if (c!=0)
s.push_back(c);
}
}
}
r.push_back(q);
ss.push_back(s);
}
return r[n];
}
int main() {
vector<int>p = { 1,5,8,9,10,17,17,20,24,30 };
vector<int>s;
int n = bucr(p, 25,s);
cout <<n<<endl;
for (auto c : s)
cout << c << endl;
}
2.矩阵链乘法
问题描述
计算n个矩阵的乘积时,由于矩阵乘法满足结合律,我们可以用括号来明确计算顺序。然而不同的计算顺序下的乘积运算的代价(标量的乘法次数)也是不同的。对于两个矩阵的乘法的计算过程如下
MATRIX-MULTIPLY(A,B)
if A.columns!=B.rows
error"incompatible dimensions"
else let C be a new A.rows * B.columns matrix
for i =1 to A.rows
for j= 1 to B.columns
cij=0
for k=1 to A.columns
cij=cij+a_ik*b_kj
return C
如果A是 p × q p\times q p×q的矩阵,B是 q × r q \times r q×r的矩阵,那么计算他们的乘积所需的运行时间由他们的标量乘法的次数pqr决定。我们用标量乘积的次数表示计算代价。
打个比方,对于矩阵链 ( A 1 , A 2 , A 3 ) (A_1,A_2,A_3) (A1,A2,A3)相乘,假设三个矩阵的规模分别为 10 × 100 , 100 × 5 , 5 × 50 10 \times 100, 100\times 5, 5\times 50 10×100,100×5,5×50,如果按顺序 ( ( A 1 A 2 ) A 3 ) ((A_1 A_2)A_3) ((A1A2)A3)计算一共需要 10 × 100 × 5 + 10 × 5 × 50 = 7500 10\times 100\times 5+10\times 5\times 50=7500 10×100×5+10×5×50=7500次标量乘法,而按顺序 ( A 1 ( A 2 A 3 ) ) (A_1(A_2A_3)) (A1(A2A3))计算的话需要 10 × 100 × 50 + 100 × 5 × 50 = 75000 10\times 100 \times 50+100\times 5 \times 50=75000 10×100×50+100×5×50=75000次标量乘法,所以第一种顺序应该比第二种顺序快大概10倍。
矩阵链乘法问题:给定n个矩阵的链 ( A 1 , . . . , A n ) (A_1,...,A_n) (A1,...,An),矩阵 A i A_i Ai的规模为 p i − 1 × p i ( 1 ≤ i ≤ n ) p_{i-1}\times p_i(1\le i\le n) pi−1×pi(1≤i≤n),求完全括号化方案使得乘积所需标量乘法的次数最少。
这个问题不难理解。我们先理论上分析这个问题。
首先通过分析发现括号化方案的数量如下:
P
(
n
)
=
{
1
n
=
1
∑
k
=
1
n
−
1
P
(
k
)
P
(
n
−
k
)
n
≥
2
P(n)=\begin{cases}1 & n=1\\ \sum_{k=1}^{n-1}P(k)P(n-k) &n\ge2 \end{cases}
P(n)={1∑k=1n−1P(k)P(n−k)n=1n≥2
这是一个递归公式(递推公式)。它的增长速度为
Ω
(
4
n
/
n
3
/
2
)
\Omega(4^n/n^{3/2})
Ω(4n/n3/2),所以括号化方案的数量是与n成指数关系的。
应用动态规划方案
我们还是四步走,按照前面给出的四个步骤来做。
- 最优括号化方案的结构特征
为了寻找这个问题的最优解,我们可以将问题划分为两个子问题,(矩阵链划分为两个子的矩阵链( A i A i + 1 . . . A k A_iA_{i+1}...A_k AiAi+1...Ak)和( A k + 1 . . . A j A_{k+1}...A_j Ak+1...Aj)),求出子问题的解然和将子问题组合起来。一定要考察所有的划分点。
- 一个递归求解方案
A
i
.
.
.
A
j
A_i...A_j
Ai...Aj的最小括号化方案的递归求解公式为
m
[
i
,
j
]
=
{
0
i
=
j
min
i
≤
k
≤
j
{
m
[
i
,
k
]
+
m
[
k
+
1
,
j
]
+
p
i
−
1
p
k
p
j
}
i
<
j
m[i,j]=\begin{cases} 0 & i=j\\ \min_{i\le k\le j}\{m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\} & i<j\end{cases}
m[i,j]={0mini≤k≤j{m[i,k]+m[k+1,j]+pi−1pkpj}i=ji<j
- 计算最优代价
我们前面提到了由两种等效的方法用动态规划构造算法。我们下面用自顶向上的方法来构造算法求最优解。
MATRIX-CHAIN-ORDER§
n=p.length-1
let m[1..n,1..n]and s[1..n-1,2..n]be new tables
for i=1 to n
m[i,i]=0
for l=2 to n
for i=1 to n-l+1
j=i+l-1
m[i,j]=infty
for k=i to j-1
q=m[i,k]+m[k+1,j]+p_i-1p_kp_j
if q<m[i,j]
m[i,j]=q
s[i,j]=k
return m and s
仔细去体会它这个计算过程,理解了以后就可以按照这个过程用程序来实现了。
C++来按照上面的伪代码实现这类问题的求解
#include<iostream>
#include<vector>
#include<climits>
using namespace std;
vector<vector<int>> mco(vector<int>p) {
int n = p.size() - 1;
vector<vector<int>>m, s;
for (int i = 0; i < n; ++i) {
vector<int>vi;
for (int j = 0; j < n; ++j) {
vi.push_back(i == j ? 0 : INT_MAX);
}
m.push_back(vi);
s.push_back(vi);
}
for (int l= 2; l <= n; ++l) {
for (int i = 1; i <= n - l + 1; ++i) {
int j = i + l-1;
for (int k = i; k <= j - 1; ++k) {
int q = m[i - 1][ k - 1] + m[k][ j - 1] + p[i - 1] * p[k] * p[j];
if (q < m[i-1][j-1]) {
m[i - 1][j - 1] = q;
s[i - 1][ j - 2] = k;
}
}
}
}
return m;
}
int main() {
vector<int>p = { 30,35,15,5,10,20,25 };
vector<vector<int>> m = mco(p);
for (auto c : m) {
for (auto cc : c)
cout << cc << "--";
cout << endl;
}
}
- 构造最优解
然后我们这一步就是根据上一步的求得的最小代价输出最优括号化方案。这里需要点技巧了。我们上一步的构造的二维数组s中保存了分割点k,我们要利用它来输出我们求得的最有括号化方案。过程如下
PRINT-OPTIMAL-PARENS(s,i,j)
if i==j
print"A"
else
print"("
PRINT-OPTIMAL-PARENS(s,i,s[i,j])
PRINT-OPTIMAL-PARENS(s,s[i,j]+1,j)
print")"
添加这个过程我们便得到了最终的完成版,直接输出我们要求的最佳方案。
#include<iostream>
#include<vector>
#include<climits>
using namespace std;
vector<vector<int>> mco(vector<int>p,vector<vector<int>>&s) {
int n = p.size() - 1;
vector<vector<int>>m;
for (int i = 0; i < n; ++i) {
vector<int>vi;
for (int j = 0; j < n; ++j) {
vi.push_back(i == j ? 0 : INT_MAX);
}
m.push_back(vi);
s.push_back(vi);
}
for (int l= 2; l <= n; ++l) {
for (int i = 1; i <= n - l + 1; ++i) {
int j = i + l-1;
for (int k = i; k <= j - 1; ++k) {
int q = m[i - 1][ k - 1] + m[k][ j - 1] + p[i - 1] * p[k] * p[j];
if (q < m[i-1][j-1]) {
m[i - 1][j - 1] = q;
s[i - 1][ j - 1] = k;
}
}
}
}
return m;
}
void pop(vector<vector<int>>s, int i, int j) {
if (i == j)
cout << "A_" << i;
else {
cout << "(";
pop(s, i, s[i-1][j-1]);
cout << "*";
pop(s, s[i-1][j-1] + 1, j);
cout << ")";
}
}
int main() {
vector<int>p = { 30,35,15,5,10,20,25 };
vector<vector<int>>s;
vector<vector<int>> m = mco(p,s);
for (auto c : m) {
for (auto cc : c)
cout << cc << "--";
cout << endl;
}
pop(s, 1,6);
}
3.动态规划原理
尽管我们通过前面两个实例对于动态规划已经有一些了解了,但我们还是不知道动态规划在解决哪些问题时有用。
通常来说,适合用动态规划求解的最优化问题需要具备两个要素:最优子结构和子问题重叠。
最优子结构
如果一个问题的最优解包含其子问题的最优解,则称其具有最优子结构。用动态规划时我们用子问题的最优解来构造原问题的最优解。
重叠子问题
问题的递归算法会反复求解相同的子问题,而不是生成新的子问题。如果递归的每一步都生成新的子问题,那么这是就适合用分治方法来求了。
备忘
另外如我们前面的钢条切割问题,我们可以通过加入备忘机制来使自顶向下的低效的递归算法变得同自顶向上方法一样可以记录子问题的解。