前言：今天接着學習動態規划算法，學習如何用動態規划來分析解決矩陣鏈乘問題。首先回顧一下矩陣乘法運算法，並給出C++語言實現過程。然后采用動態規划算法分析矩陣鏈乘問題並給出C語言實現過程。
1 #include <iostream> 2 using namespace std; 3 #define A_ROWS 3 4 #define A_COLUMNS 2 5 #define B_ROWS 2 6 #define B_COLUMNS 3 7 void matrix_multiply(int A[A_ROWS][A_COLUMNS],int B[B_ROWS][B_COLUMNS],int C[A_ROWS][B_COLUMNS]); 8 int main() 9 { 10 int A[A_ROWS][A_COLUMNS] = {1,0, 11 1,2, 12 1,1}; 13 int B[B_ROWS][B_COLUMNS] = {1,1,2, 14 2,1,2}; 15 int C[A_ROWS][B_COLUMNS] = {0}; 16 matrix_multiply(A,B,C); 17 for(int i=0;i<A_ROWS;i++) 18 { 19 for(int j=0;j<B_COLUMNS;j++) 20 cout<<C[i][j]<<" "; 21 cout<<endl; 22 } 23 return 0; 24 } 25 void matrix_multiply(int A[A_ROWS][A_COLUMNS],int B[B_ROWS][B_COLUMNS],int C[A_ROWS][B_COLUMNS]) 26 { 27 if(A_COLUMNS != B_ROWS) 28 cout<<"error: incompatible dimensions."<<endl; 29 else 30 { 31 int i,j,k; 32 for(i=0;i<A_ROWS;i++) 33 for(j=0;j<B_COLUMNS;j++) 34 { 35 C[i][j] = 0; 36 for(k=0;k<A_COLUMNS;k++) 37 C[i][j] += A[i][k] * B[k][j]; //將A的每一行的每一列與B的每一列的每一行的乘積求和 38 } 39 } 40 }
程序測試結果如下所示：
2、矩陣鏈乘問題描述
給定n個矩陣構成的一個鏈<A1,A2,A3,.......An>，其中i=1,2,...n，矩陣A的維數為pi-1pi，對乘積 A1A2...An 以一種最小化標量乘法次數的方式進行加全部括號。
注意：在矩陣鏈乘問題中，實際上並沒有把矩陣相乘，目的是確定一個具有最小代價的矩陣相乘順序。找出這樣一個結合順序使得相乘的代價最低。
3、動態規划分析過程
1）最優加全部括號的結構
動態規划第一步是尋找一個最優的子結構。假設現在要計算AiAi+1....Aj的值，計算Ai...j過程當中肯定會存在某個k值（i<=k<j）將Ai...j分成兩部分，使得Ai...j的計算量最小。分成兩個子問題Ai...k和Ak+1...j，需要繼續遞歸尋找這兩個子問題的最優解。
有分析可以到最優子結構為：假設AiAi+1....Aj的一個最優加全括號把乘積在Ak和Ak+1之間分開，則Ai..k和Ak+1..j也都是最優加全括號的。
2）一個遞歸解
設m[i,j]為計算機矩陣Ai...j所需的標量乘法運算次數的最小值，對此計算A1..n的最小代價就是m[1,n]。現在需要來遞歸定義m[i,j]，分兩種情況進行討論如下：
當i==j時：m[i,j] = 0，（此時只包含一個矩陣）
當i<j 時：從步驟1中需要尋找一個k（i≤k＜j）值，使得m[i,j] ＝min{m[i,k]+m[k+1,j]+pi-1pkpj} （i≤k＜j）。
3）計算最優代價
雖然給出了遞歸解的過程，但是在實現的時候不采用遞歸實現，而是借助輔助空間，使用自底向上的表格進行實現。設矩陣Ai的維數為pi-1pi，i=1,2.....n。輸入序列為：p=<p0，p1,...pn>，length[p] = n+1。使用m[n][n]保存m[i,j]的代價，s[n][n]保存計算m[i,j]時取得最優代價處k的值，最后可以用s中的記錄構造一個最優解。書中給出了計算過程的偽代碼，摘錄如下：
1 MAXTRIX_CHAIN_ORDER(p) 2 n = length[p]-1; 3 for i=1 to n 4 do m[i][i] = 0; 5 for t = 2 to n //t is the chain length 6 do for i=1 to n-t+1 7 j=i+t-1; 8 m[i][j] = MAXLIMIT; 9 for k=i to j-1 10 q = m[i][k] + m[k+1][i] + qi-1qkqj; 11 if q < m[i][j] 12 then m[i][j] = q; 13 s[i][j] = k; 14 return m and s;
MATRIX_CHAIN_ORDER具有循環嵌套，深度為3層，運行時間為O(n3)。如果采用遞歸進行實現，則需要指數級時間Ω(2n)，因為中間有些重復計算。遞歸是完全按照第二步得到的遞歸公式進行計算，遞歸實現如下所示：
1 int recursive_matrix_chain(int *p,int i,int j,int m[N+1][N+1],int s[N+1][N+1]) 2 { 3 if(i==j) 4 m[i][j] = 0; 5 else 6 { 7 int k; 8 m[i][j] = MAXVALUE; 9 for(k=i;k<j;k++) 10 { 11 int temp = recursive_matrix_chain(p,i,k,m,s) +recursive_matrix_chain(p,k+1,j,m,s) + p[i-1]*p[k]*p[j]; 12 if(temp < m[i][j]) 13 { 14 m[i][j] = temp; 15 s[i][j] = k; 16 } 17 } 18 } 19 return m[i][j]; 20 }
對遞歸算計的改進，可以引入備忘錄，采用自頂向下的策略，維護一個記錄了子問題的表，控制結構像遞歸算法。完整程序如下所示：
1 int memoized_matrix_chain(int *p,int m[N+1][N+1],int s[N+1][N+1]) 2 { 3 int i,j; 4 for(i=1;i<=N;++i) 5 for(j=1;j<=N;++j) 6 { 7 m[i][j] = MAXVALUE; 8 } 9 return lookup_chain(p,1,N,m,s); 10 } 11 12 int lookup_chain(int *p,int i,int j,int m[N+1][N+1],int s[N+1][N+1]) 13 { 14 if(m[i][j] < MAXVALUE) 15 return m[i][j]; //直接返回，相當於查表 16 if(i == j) 17 m[i][j] = 0; 18 else 19 { 20 int k; 21 for(k=i;k<j;++k) 22 { 23 int temp = lookup_chain(p,i,k,m,s)+lookup_chain(p,k+1,j,m,s) + p[i-1]*p[k]*p[j]; //通過遞歸的形式計算，只計算一次，第二次查表得到 24 if(temp < m[i][j]) 25 { 26 m[i][j] = temp; 27 s[i][j] = k; 28 } 29 } 30 } 31 return m[i][j]; 32 }
4)構造一個最優解
第三步中已經計算出來最小代價，並保存了相關的記錄信息。因此只需對s表格進行遞歸調用展開既可以得到一個最優解。書中給出了偽代碼，摘錄如下：
1 PRINT_OPTIMAL_PARENS(s,i,j) 2 if i== j 3 then print "Ai" 4 else 5 print "("; 6 PRINT_OPTIMAL_PARENS(s,i,s[i][j]); 7 PRINT_OPTIMAL_PARENS(s,s[i][j]+1,j); 8 print")";
4、編程實現
采用C++語言實現這個過程，現有矩陣A1(30×35)、A2(35×15)、A3(15×5)、A4(5×10)、A5(10×20)、A6(20×25)，得到p=<30,35,15,5,10,20,25>。實現過程定義兩個二維數組m和s，為了方便計算其第一行和第一列都忽略，行標和列標都是1開始。完整的程序如下所示：
1 #include <iostream> 2 using namespace std; 3 4 #define N 6 5 #define MAXVALUE 1000000 6 7 void matrix_chain_order(int *p,int len,int m[N+1][N+1],int s[N+1][N+1]); 8 void print_optimal_parents(int s[N+1][N+1],int i,int j); 9 10 int main() 11 { 12 int p[N+1] = {30,35,15,5,10,20,25}; 13 int m[N+1][N+1]={0}; 14 int s[N+1][N+1]={0}; 15 int i,j; 16 matrix_chain_order(p,N+1,m,s); 17 cout<<"m value is: "<<endl; 18 for(i=1;i<=N;++i) 19 { 20 for(j=1;j<=N;++j) 21 cout<<m[i][j]<<" "; 22 cout<<endl; 23 } 24 cout<<"s value is: "<<endl; 25 for(i=1;i<=N;++i) 26 { 27 for(j=1;j<=N;++j) 28 cout<<s[i][j]<<" "; 29 cout<<endl; 30 } 31 cout<<"The result is:"<<endl; 32 print_optimal_parents(s,1,N); 33 return 0; 34 } 35 36 void matrix_chain_order(int *p,int len,int m[N+1][N+1],int s[N+1][N+1]) 37 { 38 int i,j,k,t; 39 for(i=0;i<=N;++i) 40 m[i][i] = 0; 41 for(t=2;t<=N;t++) //當前鏈乘矩陣的長度 42 { 43 for(i=1;i<=N-t+1;i++) //從第一矩陣開始算起，計算長度為t的最少代價 44 { 45 j=i+t-1;//長度為t時候的最后一個元素 46 m[i][j] = MAXVALUE; //初始化為最大代價 47 for(k=i;k<=j-1;k++) //尋找最優的k值，使得分成兩部分k在i與j-1之間 48 { 49 int temp = m[i][k]+m[k+1][j] + p[i-1]*p[k]*p[j]; 50 if(temp < m[i][j]) 51 { 52 m[i][j] = temp; //記錄下當前的最小代價 53 s[i][j] = k; //記錄當前的括號位置，即矩陣的編號 54 } 55 } 56 } 57 } 58 } 59 60 //s中存放着括號當前的位置 61 void print_optimal_parents(int s[N+1][N+1],int i,int j) 62 { 63 if( i == j) 64 cout<<"A"<<i; 65 else 66 { 67 cout<<"("; 68 print_optimal_parents(s,i,s[i][j]); 69 print_optimal_parents(s,s[i][j]+1,j); 70 cout<<")"; 71 } 72 73 }
程序測試結果如下所示：
5、總結
動態規划解決問題關鍵是分析過程，難度在於如何發現其子問題的結構及子問題的遞歸解。這個需要多多思考，不是短時間內能明白。在實現過程中遇到問題就是數組，數組的下標問題是個比較麻煩的事情，如何能夠過合理的去處理，需要一定的技巧。