动态规划概论

动态规划(dynamic programming)是运筹学的一个分支,是求解决策过程(decision process)最优化的数学方法。20世纪50年代初美国数学家R.E.Bellman等人在研究多阶段决策过程(multistep decision process)的优化问题时,提出了著名的最优化原理(principle of optimality),把多阶段过程转化为一系列单阶段问题,利用各阶段之间的关系,逐个求解,创立了解决这类过程优化问题的新方法——动态规划。

钢条切割问题

算法运行时间分析

如果采用最平凡的递推算法,由$r_n= \underset{1 \leq i \leq n}{\max} (p_i+r_{n-i})$
令$T(n)$表示算法第二个参数值为n的时候CUT-ROD的调用次数,可以推出:$T(n)=1+ \sum_{j=0}^{n-1}T(j)$

$T(n)=1+T(1)+T(2)+\cdots+T(n-1)$
$T(n-1)=1+T(1)+T(2)+\cdots+T(n-2)$

即$T(n)=2+2T(1)+\cdots+2T(n-2)=2T(n-1)$

$\frac{T(n)}{T(n-1)}=2$

由等比数列的递推公式,$T(n)=2^n$
所以,平凡的递归算法反复求解相同的子问题,耗费了大量的时间。

使用动态规划

01

求解子问题的两种方法:
使用memorized_cut_rod方法的时候,memorized_cut_rod_aux(p,n,r)的递归,最多递归到r[]的下标为n-1,因为在循环中执行的是for(int i=0;i<n;i++)

02

使用自底向上递归求解。

03

由上图所示,在递归过程中:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
r[0]=0;
for(int j=0;j<n;j++)
{
int result=-INFINITY;
for(int i=0;i<=j;i++)
{
//获取scale=j时候的max值,即最大收益result
//随着j++,我们获取最大收益的时候要返回r[j+1]的值,这样
//scale=scale+1时候的最大收益,即r[j+1],即为上一步求出的
//result值,即r[j+1]=result
result=max(result,p[i]+r[j-1]);
}
r[j+1]=result;

}
return r[n];
//返回的r[n]这里的n指问题的规模

算法实现过程:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
#include <iostream>
#define INFINITY 0x7fffffff
using namespace std;

int max(int a,int b)
{
return a>b?a:b;
}

int memorized_cut_rod_aux(int p[],int n,int r) //n表示问题的规模
{
//注意,这里p[i]的意义和r[i]的意义不一样
//p[i]表示对应的数组元素,i表示下标
//r[i]表示问题的规模为i所对应的解
//p[i]对应的问题是第i个,其规模为i+1,剩余的问题规模是n-(i+1)
//从另一个角度,0<=i<=n-1 0<=r<=n-1,所以递归的时候是(p,n-1-i,r)
int result;
if(r[n]>=0)
return r[n]; //表示问题已经被求解完毕了,返回该值

if(n==0)
result=0; //表示没有切割
else
{
result=-INFINITY;
for(int i=0;i<n;i++)
{
result=max(result,p[i]+memorized_cut_rod_aux(p,n-1-i,r));
}
}
//表示规模为n的问题求解完毕
r[n]=result;
return result;

}
int memorized_cut_rod(int p[],int n)
{
int *r=new int[n];
for(int i=0;i<n;i++)
{
r[i]=-INFINITY;
}

return memorized_cut_rod_aux(p,n,r); //这里(p,n,r)
//r[]的递归标最多为n-1
}

int BOTTOM_UP_CUT_ROD(int p[],int n)
{
int *r=new int [n];
r[0]=0;
for(int j=0;j<n;j++)
{
int result=-INFINITY;
for(int i=0;i<=j;i++)
{
//获取scale=j时候的max值,即最大收益result
//随着j++,我们获取最大收益的时候要返回r[j+1]的值,这样
//scale=scale+1时候的最大收益,即r[j+1],即为上一步求出的
//result值,即r[j+1]=result
result=max(result,p[i]+r[j-1]);
}
r[j+1]=result;

}
return r[n];
//返回的r[n]这里的n指问题的规模
}

主函数:

1
2
3
4
5
6
7
int main()
{
const int n=10;
int p[10]={1,5,8,9,10,17,17,20,24,30};
cout<<memorized_cut_rod(p,9)<<endl;
cout<<BOTTOM_UP_CUT_ROD(p,9)<<endl;
}

子问题图与重构解

这里特别说明一种表示方式

1
array *EXTENDED_BOTTOM_UP_CUT_ROD(int p[],int n)

这里p[]表示指向数组第一个元素的指针,注意区分以下两种表示方法:当表示第n个元素的时候

1
2
*(p+n)->element
p[n].element

这两种表示方式不同,一种是指针,一种是结构体表示。

具体的说明见下图:

04

动态规划解决钢条切割问题,问题的输出:注意让规模逐渐减小。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <iostream>
using namespace std;
#define INFINITY 0x7fffffff

struct rod
{
int result; //代表最大收益
int solution; //代表切割方案
};

rod *EXTENDED_BOTTOM_UP_CUT_ROD(int p[], int n)
{
rod *res=new rod[n];
res[0].result=0; //初始化,钢条规模为0,收益为0
int q;
for(int j=0;j<n;j++)
{
q=-INFINITY;
for(int i=0;i<=j;i++)
{
if(q<p[i]+res[j-i].result)
{
q=p[i]+res[j-i].result;
res[j+1].solution=i+1; //保存切割方式
} //这里的i表示切割下来的片段,规模有多大?
}
//此时已经完成了从1---j的循环
res[j+1].result=q;
}
return res+n; //注意:这里res表示指向第一个元素的指针
}
//res+n表示指针移动n位,新的指针指向n之后的那个位置
//rod *res=new rod[n] res[i].result (res+i)->result
//两种用法注意区分

//特别注意:这里res[0].result=0 最后res[j+1].result
//res[j+1].solution,意味这res[]的存储范围如下:
//res[1],res[2],......,res[n]
//res初始指针位置为res[0],res末指针的位置为res[n]

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
rod *res=EXTENDED_BOTTOM_UP_CUT_ROD(p,n);
//注意这里函数的返回值为res+n,即最后一个数组元素的地址
while(n>0)
{
cout<<(*res).solution<<" ";
n=n-(*res).solution;
//注意返回的res指的是res[n]的位置
//返回的是规模为n时候的切割方法,从最大规模往最小规模输出
//这里从大规模往小规模输出,指针也要调整
//从规模为n的位置逐渐减小到i,再逐渐减小到1,0
res=res-(*res).solution;

}
}

void main()
{
const int n=10;
int p[10]={1,5,8,9,10,17,17,20,24,30};
cout<<(*EXTENDED_BOTTOM_UP_CUT_ROD(p,4)).result<<endl;
PRINT_CUT_ROD_SOLUTION(p,4);
}

其他问题

15.1-2

反例如下图:

05

15.1-3 钢条切割的时候考虑固定成本为c

特别注意,这里的固定成本为c,在执行循环的时候,存在表达式的不同!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
//当j==i的时候有不同!
for(int j=0;j<n;j++)
{
//省略的code
for(int i=0;i<j;i++)
{
//i<j时候表示钢条一定会进行切割!
if(q<=p[i]+res[j-i].result-c)
{
q=p[i]+res[j-i].result-c;
res[j+1].solution=i+1;
cut_or_not=true; //表示进行切割了
}
}

//当i==j时候的情况不一样了!
//此时已经保存的最大收益是q
//很有可能
//p[i]+res[j-i].result-c<q<=p[i]+res[j-i].result
//这个时候j==i可能没有发生切割,没有发生切割和发生切割
//收益的表达式不同,一个减去成本,一个没有
if(j==i)
{
if(q<=p[i]+res[j-i].result && cut_or_not==false)
{
res[j+1].solution=i+1;
}
}

//j求解完毕
res[j+1].result=q;

}

return res+n;

具体实现的过程如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <iostream>
using namespace std;
#define INFINITY -0x7fffffff
#define cost 2

int max(int a, int b)
{
return a>b?a:b;
}

struct rod_cost
{
int result;
int solution;
bool cut_or_not;
};

rod_cost *EXTENDED_BOTTOM_UP_CUT_ROD(int p[], int n)
{
rod_cost *res=new rod_cost[n];
res[0].result=0;
int q;

for(int j=0;j<n;j++)
{
q=-INFINITY;
bool cut_or_not=false;
for(int i=0;i<j;i++)
{
if(q<=p[i]+res[j-i].result-cost)
{
q=p[i]+res[j-i].result-cost;
res[j+1].solution=i+1;
cut_or_not=true;
}
}

//每一次循环,j=0----j=n-1,表示问题的规模从1---n
//每一个规模,cut_or_not都初始化为false,最后看是否切割?
//如果发生了切割,在i的循环中,也就是j的子问题中
//让cut_or_not=true,最后判断是否切割?

if(j==i)
{
if(q<=p[i]+res[j-i].result && cut_or_not==false)
{
res[j+1].solution=i+1;
}
}

//求解完毕
res[j+1].result=q;
}

return res+n;
}

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
rod_cost *res=EXTENDED_BOTTOM_UP_CUT_ROD(p,n);
while(n)
{
cout<<(*res).solution<<" "
n=n-(*res).solution;
res=res-(*res).solution;
}
}

void main()
{
const int n=10;
int p[10]={1,5,8,9,10,17,17,20,24,30};
cout<<(*EXTENDED_BOTTOM_UP_CUT_ROD(p,10)).result<<endl;
PRINT_CUT_ROD_SOLUTION(p,10);
}

15.1-4修改memorized_cut_rod,使之返回最优收益值和切割方案

实现函数如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <iostream>
using namespace std;
#define INFINITY 0x7fffffff

struct rod
{
int result;
int solution;
};

int max(int a, int b)
{
return a>b?a:b;
}

rod *memorized_cut_rod_aux(int p[], int n, rod res[])
{
int q;
if(res[n].result>=0)
return res+n;

if(n==0)
q=0; //规模为0的时候,自然收益为0
else
{
q=-INFINITY;
for(int i=0;i<n;i++)
{
//会犯错误的地方:
//这里不能够直接写
//q=max(q,p[i]+(memorized_cut_rod(p,n-1-i,res))->solution)
//因为还要进行储存切割方案,在子问题规模为i+1的时候,有可能进行切割,也有可能不切割
int tmp=p[i]+memorized_cut_rod_aux(p,n-1-i,res)->solution;
//在p[i]处的最优方案是tmp,tmp和q进行比较判断后,推出在i+1处的最优方案

//这里把n分解成更小的子问题,怎么解决?通过递归解决
//memorized_cut_rod_aux(p,n-1-i,res)把n-1的子问题转化为n-1-i
if(q<tmp)
{
q=tmp;
res[n].solution=i+1;
//我们这里返回的就是问题规模为n时候的切割方案,这里的切割方案是:
//切割一次or无切割!
//i<n-1时候,q=tmp表示切割
//i==n-1的时候,tmp=p[n-1]+memorized_cut_rod(p,0,res)->solution
//如果此时还满足q<tmp,res[n].solution=n-1+1=n,表示无切割
//这一步的目的在于保存每一个循环过程中i的最优方案
}
}
}

res[n].result=q;
return res+n;
}

rod *memorized_cut_rod(int p[], int n)
{
rod *res=new rod[n];
for(int i=0;i<=n;i++)
{
res[i].result=-INFINITY;
}

return memorized_cut_rod_aux(p,n,res);
}

void PRINT_CUT_ROD_SOLUTION(int p[], int n)
{
rod *res=memorized_cut_rod(p,n);
cout<<"max result"<<res->result<<endl;
cout<<"solution: "<<endl;
while(n)
{
cout<<(*res).solution<<" ";
n=n-(*res).solution;
res=res-(*res).solution;
}
}

//重点:在于区分数组下标和子问题的规模,数组下标为i-1,子问题的规模为i

void main()
{
const int n=10;
int p[10]={1,5,8,9,10,17,17,20,24,30};
PRINT_CUT_ROD_SOLUTION(p,10);
cout<<endl;
}

15.1-5斐波那契数列

子问题图:

06

实现方法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <iostream>
using namespace std;

int fibonacci(int array[],int n)
{
if(n>=0)
array[0]=0;
if(n>=1)
array[1]=1;

for(int i=2;i<n;i++)
{
array[i]=array[i-1]+array[i-2];
}
return array[n];
}

特别说明:

1
2
3
4
5
6
7
8
9
10
11
for(int j=0;j<n;j++)
{
//code
for(int i=0;i<=j;i++)
{
q=max(q,p[i]+r[j-i]);
}
r[j]=q; //这样写也是可\ \;\;以的
}

return r[n-1]; //只是返回值会不一样

矩阵链乘法

矩阵规模序列

$A1 \qquad\ \;\;|A2 \qquad\ \;\;|A3 \qquad\ \;\;|$
$\underline{10}\times100\;\;|100\times5\;\;\;|\underline{5\times50}\quad\;\;\,|$

$A1A2 \qquad\times A3=10\times5\times50$

由此可以推出:
$A_{i \ldots k} \qquad \qquad \times \qquad \qquad A_{k+1 \ldots j}$
$=\underline{p_{i-1}} \times p_i \times p_{i+1} \times \qquad\quad \underline{p_k} \times p_{k+1}\times\underline{p_j}$
$=p_{i-1} \times p_k \times p_j$

由矩阵规模序列推出矩阵链乘法,如上图所示,相乘的代价是:把有下划线的项相乘的值。
可以得到递推式如下:

矩阵链乘法的运算过程

如下图所示:

08

可以看出,三角形处的min值取决于两个圆圈处的值,具体的实现过程如下:
line01:

line02:括号中表示划分方式

line03:

line04:

line05:

计算最优代价

09

具体的实现过程如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#include <iostream>
#include <vector>
#include <stdlib.h>
#include <ctime>
using namespace std;
#define n 7
#define INFINITY 0x7fffffff

struct Matrix
{
int rows; //表示行数
int columns; //表示列数
vector<vector<int> > matrix;
Matrix(int r,int c)
{
rows=r;
columns=c;
matrix.resize(rows);
for(int i=0;i<rows;i++)
matrix[i].resize(columns);
}
};

struct Matrix_Chain
{
vector<vector<int> > m; //表示运算次数
vector<vector<int> > s; //划分方式
Matrix_Chain()
{
m.resize(n-1); //有n-1个括号
for(int k=0;k<n-1;k++)
{
m[k].resize(n-1); //最后输出的是一个对角阵
//辅助表m[1...n-1, 1...n-1]来保存代价
//s[1...n-1,2...n]记录最优值m[i,j]对应的分割点k
}
s.resize(n-1);
for(int t=0;t<n-1;t++)
{
s[t].resize(n-1);
}
}
};

Matrix init(Matrix A)
{
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<A.columns;j++)
{
A.matrix[i][j]=rand()%10;
}
}
return A;
}

Matrix Matrix_Mutiply(Matrix A,Matrix B)
{
if(A.columns!=B.rows)
{
Matrix D(1,1);
D.matrix[0][0]=0;
cerr<<"incompatible dimensions"<<endl;
return D;
}
else
{
Matrix C(A.rows, B.columns);
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<B.columns;j++)
{
C.matrix[i][j]=0;
for(int k=0;k<A.columns;k++)
C.matrix[i][j]=C.matrix[i][j]+A.matrix[i][k]*B.matrix[k][j];
}
}
return C;
}
}

Matrix_Chain Matrix_Chain_Order(int p[])
{
Matrix_Chain T;
int len=n-1; //辅助表m[k,k]来保存代价,其中k=n-1,n为矩阵链长度
for(int i=0;i<len;i++) //s[k,k]用来保存最优值m[i,j]对应的分割点
{
T.m[i][i]=0;
}
for(int l=2;l<=len;l++) //l表示矩阵链长度
{
for(int i=1;i<=len-l+1;i++)
{
int j=i+l-1;
T.m[i-1][j-1]=INFINITY; //用来计算[i,i+1..k] [k+1..,j-1,j]这样划分的运算代价
for(int k=i;k<j;k++) //循环中的i表示第几个数,T.m[i]中的i表示数组下标
//千万注意,这里k从i开始!
{
int q=T.m[i-1][k-1]+T.m[k][j-1]+p[i-1]*p[k]*p[j];
if(q<T.m[i-1][j-1])
{
T.m[i-1][j-1]=q;
T.s[i-1][j-1]=k-1;
}
}
}
}
return T;
}

void print(Matrix A)
{
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<A.columns;j++)
cout<<A.matrix[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}

int main()
{
int p[n]={5,10,3,12,5,50,6};
Matrix_Chain T=Matrix_Chain_Order(p);
for(int i=0;i<n-1;i++)
{
for(int j=0;j<n-1;j++)
{
if(T.m[i][j]<0)
cout<<"-1"<<"\t";
else cout<<T.m[i][j]<<"\t";
}
cout<<endl;
}
}

运算结果:

10

矩阵链总乘法

矩阵链总乘法的递归表示:

11

递归运算的核心:递归进行到最底层的时候,算法的运行情况?
这里的运算是矩阵乘法:

1
2
3
4
if(j==i+1)
{
return Matrix_Mutiply(A[i],A[j]);
}

相当于二路归并排序,递归进行到最末端时候的情况。A[i]*A[i+1]
在每一次的递归过程中,划分的下标信息,是存储在S[1,…,n]中的,我们在函数中要传递该参数。

具体实现过程:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
#include <iostream>
#include <vector>
#include <stdlib.h>
#include <ctime>
#include <memory>
using namespace std;
#define n 7
#define INFINITY 0x7fffffff

struct Matrix
{
int rows; //表示行数
int columns; //表示列数
vector<vector<int> > matrix;
};

struct Matrix_Chain
{
vector<vector<int> > m; //表示运算次数
vector<vector<int> > s; //划分方式
Matrix_Chain()
{
m.resize(n-1); //有n-1个括号
for(int k=0;k<n-1;k++)
{
m[k].resize(n-1); //最后输出的是一个对角阵
//辅助表m[1...n-1, 1...n-1]来保存代价
//s[1...n-1,2...n]记录最优值m[i,j]对应的分割点k
}
s.resize(n-1);
for(int t=0;t<n-1;t++)
{
s[t].resize(n-1);
}
}
};

Matrix init(Matrix &A,int r,int c) //对A的矩阵具体值的初始化
{
A.rows=r;
A.columns=c;
A.matrix.resize(A.rows);
for(int i=0;i<A.rows;i++)
A.matrix[i].resize(A.columns);

srand((unsigned)time(0));
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<A.columns;j++)
{
A.matrix[i][j]=rand()%n+1;
}
}
return A;
}

void print(Matrix A)
{
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<A.columns;j++)
cout<<A.matrix[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}

Matrix Matrix_Mutiply(Matrix A,Matrix B)
{
if(A.columns!=B.rows)
{
Matrix D;
init(D,1,1);
D.matrix[0][0]=0;
cerr<<"incompatible dimensions"<<endl;
return D;
}
else
{
Matrix C;
init(C,A.rows,B.columns);
for(int i=0;i<A.rows;i++)
{
for(int j=0;j<B.columns;j++)
{
C.matrix[i][j]=0;
for(int k=0;k<A.columns;k++)
C.matrix[i][j]=C.matrix[i][j]+A.matrix[i][k]*B.matrix[k][j];
}
}
return C;
}
}

Matrix_Chain Matrix_Chain_Order(int p[]) //这里的i j统一表示是:数组中的第几个数?
{ //所以,数组下标是T.m[i-1][j-1]
Matrix_Chain T;
int len=n-1; //辅助表m[k,k]来保存代价,其中k=n-1,n为矩阵链长度
for(int i=0;i<len;i++) //s[k,k]用来保存最优值m[i,j]对应的分割点
{
T.m[i][i]=0;
}
for(int l=2;l<=len;l++) //l表示矩阵链长度
{
for(int i=1;i<=len-l+1;i++)
{
int j=i+l-1;
T.m[i-1][j-1]=INFINITY; //用来计算[i,i+1..k] [k+1..,j-1,j]这样划分的运算代价
for(int k=i;k<j;k++) //循环中的i表示第几个数,T.m[i]中的i表示数组下标
//千万注意,这里k从i开始!
{
int q=T.m[i-1][k-1]+T.m[k][j-1]+p[i-1]*p[k]*p[j];
if(q<T.m[i-1][j-1])
{
T.m[i-1][j-1]=q;
T.s[i-1][j-1]=k-1;
}
}
}
}
return T;
}


//矩阵链总乘法,相当于二路归并
//在二路归并排序中,递归处理,让区间的长度从l逐渐减小到0
//递归到尾端的时候,是两个数相乘A[i]*A[j] (j==i+1)

//函数参数:矩阵链乘法中,1、传递A[]矩阵的原始数据,2、括号的划分情况,存储在数组T.s[]中
Matrix Matrix_Chain_Multiply(Matrix A[],Matrix_Chain T,int i,int j)
{
if(j==i)
return A[i];
if(j==i+1)
return Matrix_Mutiply(A[i],A[j]);
Matrix t1=Matrix_Chain_Multiply(A,T,i,T.s[i][j]); //矩阵T.s[i][j]表示A[i...j]的划分方式
Matrix t2=Matrix_Chain_Multiply(A,T,T.s[i][j]+1,j);
return Matrix_Mutiply(t1,t2);
}

//同理,打印输出矩阵链乘法也是如此:
void Print_Optimal_Parents(Matrix_Chain T,int i,int j)
{
if(i==j)
cout<<"A"<<i<<" ";
else
{
cout<<"("<<" ";
Print_Optimal_Parents(T,i,T.s[i][j]);
Print_Optimal_Parents(T,T.s[i][j]+1,j);
cout<<")"<<" ";
}
}

int main()
{
Matrix test;
init(test,4,5);
print(test);
cout<<"Program begins: "<<endl;

int p[n]={5,10,3,12,5,50,6};
Matrix_Chain T=Matrix_Chain_Order(p);

for(int i=0;i<n-1;i++) //矩阵的dimension,等于矩阵链规模n-1
{
for(int j=0;j<n-1;j++) //原矩阵对应的值A[i....j]=A[0...n-2]
{
if(T.m[i][j]<0)
cout<<"-1"<<"\t";
else
cout<<T.m[i][j]<<"\t";
}
cout<<endl;
}

Print_Optimal_Parents(T,0,n-2);
cout<<endl;

cout<<"Concret Implement:"<<" "<<endl;

Matrix A[n]; //这里的A存放的是指向Matrix数组第一个元素的指针
//A[]的dimension的值,具体一点为p(i-1)*p(i)
for(int j=1;j<n;j++)
{
Matrix t;
init(t,p[j-1],p[j]); //A[rows*columns]--->p[j-1]*p[j]
A[j-1]=t; //p[]仅仅表示A[]的维度,具体的值用init初始化
cout<<endl;
print(A[j-1]);
cout<<endl;
}

Matrix result=Matrix_Chain_Multiply(A,T,0,n-2);
print(result);

return 0;
}

算法运行情况:

12

13

问题解答

15.2-3 用代入法证明递归公式(15.6)的结果是:$\Omega(2^n)$
$n=2\qquad\quad P(2)=p^2(1)=1$
$n=3\qquad\quad P(3)=2P(1)P(2)$
$n=4\qquad\quad P(4)=2^2P^2(1)P(2)+P^2(2)$
$n=5\qquad\quad P(5)=2^3P^3(1)P(2)+(2+2^2)P(1)P^2(2)$
$\cdots\cdots$

$P(n)=2^{[n/2]+1}P^{[n/2]+1}(1)P(2)+(2+2^2+\cdots+2^{[n/2]})P(1)P^2(2)$
由等比数列的求和公式,可以推出:
$P(n)=\Omega(2^n)$

15.2-4 对输入链长度为n的矩阵链乘法问题,包含多少个顶点?包含多少条边?这些边分别连接哪些顶点?
值得注意的是,矩阵链不等于二叉树,子问题就对应矩阵中的每个点。

如下图所示:

14

子问题图包含$n^2/2=O(n^2)$个顶点。
边数?第i层的边数为$i(i-1)$

总边数为$\sum_{i=1}^{i=n}i(i-1)=O(n^3)$

15.2-5 计算其他表项的时候,访问m[i,j],主要的时间代价是由
$q=m[i,k]+m[k+1,j]+p_{i-1}p_kp_j$ 产生的,每一次访问都对m[]调用2次。

总的遍历循环次数,为k可以取到的值的总数,即为矩阵链长度length
令$l=length$
$\sum_{i=1}^n\sum_{j=1}^nR(i,j)=\sum_{l=2}^{l=n}\sum_{i=1}^{n-l+1}2(l-1)=2\sum_{l=2}^n(n-(l-1))(l-1)=2\sum_{t=1}^{n-1}(n-t)t$
$=O(2\frac{n(n-1)(2n-1)}{6})=\frac{n^3-n}{3}$

其实本质上还是这张图:

08

图中红色线段长度表示$l$,即$\sum_{l=1}{l=n}$,每一个红色线段对应两个m[]点,在每一层,有$l-1$个三角形的点
所以,总代价是:
$\sum_{l=2}^{l=n}2 \times l(l-1)=2 \times \frac{n(n-1)(2n-1)}{6}=O(\frac{n^3-n}{3})$

15.2-6 对n个元素的表达式进行完全括号化,恰好需要n-1对括号

数学归纳法,k个元素要k-1对括号,当元素增加到k+1时,保持k-1对括号不变,再添加新的一个元素,括号数为k