9.矩陣鏈乘法
其實這個過程不是很理解，書裡的考究是使其能滿足最優解，就是我們之前一直說的動態規劃情況下滿足的最優解。
給定一個n個矩陣的序列（矩陣鏈）<A1,A2,...,An>，我們希望計算它們的乘積 A1A2...An，為了計算表示式，我們可以先用括號明確計算次序，然後利用標準的矩陣相乘演算法進行計算。 例如如果有矩陣鏈為<A1,A2,A3,A4>,則共有5種完全括號化的矩陣乘積鏈。
(A1(A2(A3A4)))、(A1((A2A3)A4))、((A1A2)(A3A4))、((A1(A2A3))A4)、((A1A2)A3)A4) 對矩陣鏈加括號的方式會對乘積運算的代價產生巨大影響。我們先來分析兩個矩陣相乘的代價。下面的虛擬碼的給出了兩個矩陣相乘的標準演算法，屬性rows和columns是矩陣的行數和列數。
演算法虛擬碼：
MATRIX-MULTIPLY(A,B)
1. if columns[A] !=row[B]
2. then error"矩陣乘法不相容"
3. else for i<-1 to rows[A]
4. do for j<- 1 to columns[B]
5. do C[i,j] <-0
6. for k<- 1 to columns[A]
7. do C[i,j]<-C[i,j]+A[i,k]*B[k,j]
8. return C
兩個矩陣A和B只有相容(compatible)，即A的列數等於B的行數時，才能相乘。如果A是p×q的矩陣，B是q×r的矩陣，那麼乘積C是p×r的矩陣。計算C所需要時間由第8行的標量乘法的次數決定的，即pqr。以矩陣鏈<A1,A2,A3>為例，來說明不同的加括號方式會導致不同的計算代價。假設三個矩陣的規模分別為10×100、100×5和5×50。如果按照((A1A2)A3)的順序計算，為計算A1A2(規模10×5)，需要做10*100*5=5000次標量乘法，再與A3相乘又需要做10*5*50=2500次標量乘法，共需7500次標量乘法。 如果按照(A1(A2A3))的順序計算，為計算A2A3(規模100×50)，需100*5*50=25000次標量乘法，再與A1相乘又需10*100*50=50000次標量乘法，共需75000
矩陣鏈乘法問題(matrix-chain multiplication problem)可描述如下：給定n個矩陣的鏈<A1,A2,...,An>，矩陣Ai的規模為p(i-1)×p(i) (1<=i<=n)，求完全括號化方案，使得計算乘積A1A2...An所需標量乘法次數最少。 步驟1：最優括號化方案的結構特徵 動態規劃的第一步是尋找最優子結構，然後就可以利用這種子結構從子問題的最優解構造出原問題的最優解。在矩陣鏈乘法問題中，我們假設A(i)A(i+1)...A(j)的最優括號方案的分割點在A(k)和A(k+1)之間。那麼，繼續對“字首”子鏈A(i)A(i+1)..A(k)進行括號化時，我們應該直接採用獨立求解它時所得的最優方案。 我們已經看到，一個非平凡(i≠j)的矩陣鏈乘法問題例項的任何解都需要劃分鏈，而任何最優解都是由子問題例項的最優解構成的。為了構造一個矩陣鏈乘法問題例項的最優解，我們可以將問題劃分為兩個子問題(A(i)A(i+1)...A(k)和A(k+1)A(k+2)..A(j)的最優括號化問題)，求出子問題例項的最優解，然後將子問題的最優解組合起來。我們必須保證在確定分割點時，已經考察了所有可能的劃分點，這樣就可以保證不會遺漏最優解。
步驟2：一個遞迴求解方案 下面用子問題的最優解來遞迴地定義原問題最優解的代價。對於矩陣鏈乘法問題，我們可以將對所有1<=i<=j<=n確定A(i)A(i+1)...A(j)的最小代價括號化方案作為子問題。令m[i,j]表示計算矩陣A(i..j)所需標量乘法次數的最小值，那麼，原問題的最優解—計算A(1..n)所需的最低代價就是m[1,n]。 我們可以遞迴定義m[i,j]如下。對於i=j時的平凡問題，矩陣鏈只包含唯一的矩陣A(i..j)=A(i)，因此不需要做任何標量乘法運算。所以，對所有i=1,2,...,n，m[i,i]=0。若i<j，我們利用步驟1中得到的最優子結構來計算m[i,j]。我們假設A(i)A(i+1)...A(j)的最優括號化方案的分割點在矩陣A(k)和A(k+1)之間，其中i<=k<j。那麼，m[i,j]就等於計算A(i..k)和A(k+1..j)的代價加上兩者相乘的代價的最小值。由於矩陣Ai的大小為p(i-1)*pi，易知A(i..k)和A(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) 此遞迴公式假定最優分割點k是已知的，但實際上我們是不知道。不過，k只有j-i種可能的取值，即k=i,i+1,...,j-1。由於最優分割點必在其中，我們只需檢查所有可能情況，找到最優者即可。 因此，A(i)A(i+1)...A(j)的最小代價括號化方案的遞迴求解公式變為： ①如果i=j，m[i,j]=0
②如果i<j，m[i,j]=min{m[i,k]+m[k+1,j]+p(i-1)p(k)p(j)} i<=k<j
m[i,j]的值給出了子問題最優解的代價，但它並未提供足夠的資訊來構造最優解。為此，我們用s[i,j]儲存最優括號化方案的分割點位置k，即使得m[i,j]=m[i,k]+[k+1,j]+p(i-1)p(k)p(j)成立的k值。
步驟3：計算最優代價 現在，我們可以很容易地基於遞迴公式寫出一個遞迴演算法，但遞迴演算法是指數時間的，並不必檢查若有括號化方案的暴力搜尋方法更好。注意到，我們需要求解的不同子問題的數目是相對較少的：每對滿足1<=i<=j<=n 的i和j對應一個唯一的子問題，共有n^2(最少)。遞迴演算法會在遞迴呼叫樹的不同分支中多次遇到同一個子問題。這種子問題重疊的性質是應用動態規劃的另一標識(第一個標識是最優子結構)。 我們採用自底向上表格法代替遞迴演算法來計算最優代價。此過程假定矩陣Ai的規模為p(i-1)×pi(i=1,2,...,n)。它的輸入是一個序列p=<p0,p1,...,pn>，其長度為p.length=n+1。過程用一個輔助表m[1..n,1..n]來儲存代價m[i,j]，用另一個輔助表s[1..n-1,2..n](s[1,2]..s[n-1,n]這裡i<j)記錄最優值m[i,j]對應的分割點k。我們就可以利用表s構造最優解。
對於矩陣A(i)A(i+1)...A(j)最優括號化的子問題，我們認為其規模為鏈的長度j-i+1。因為j-i+1個矩陣鏈相乘的最優計算代價m[i,j]只依賴於那麼少於j-i+1個矩陣鏈相乘的最優計算代價。因此，演算法應該按長度遞增的順序求解矩陣鏈括號化問題，並按對應的順序填寫表m。
演算法虛擬碼：
MATRIX-CHAIN-ORDER(p)
1. n<-length[p]-1
2. for i<-1 to n
3. do m[i,i]<-0
4. for l<-2 to n // l 為矩陣鏈的長度
5. do for i<-1 to n-l+1
6. do j<-i+l-1
7. m[i,j]<-∞
8. for k<-i to j-1
9. do q<-m[i,k]+m[k+1,j]+p(i-1)p(k)p(j)
10. if q<m[i,j]
11. then m[i,j]<-q //標量乘法運算數
12. s[i,j]<-k //分割點位置
13. return m and s
下圖展示了對一個長度為6的矩陣鏈執行此演算法的過程：
此圖展示了將主對角線置於水平狀態時的效果，計運算元鏈AiAi+1...Aj的乘積的最小代價m[i,j]可在Ai行與Aj行的交匯處找到。即圖中的頂點15125.
步驟4：構造最優解
雖然MATRIX_CHAIN_ORDER求出了計算矩陣鏈乘積所需的最少標量乘法運算次數，但它並未直接指出如何進行這種最優代價的矩陣鏈乘法計算。表s[i,j]記錄了一個k值，指出A(i)A(i+1)...A(j)的最優括號化方案的分割點應在A(k)和A(k+1)之間。
因此，我們A(1..n)的最優計算方案中最後一次矩陣乘法運算應該是以s[1,n]為分界的A(1..s[1,n])*A(s[1,n]+1..n)。我們可以用相同的方法遞迴地求出更早的矩陣乘法的具體計算過程，因為s[1,s[1,n]]指出了計算A(1..s[1,n])時應進行的最後一次矩陣乘法執行；s[s[1,n]+1,n]指出了計算A(s[1,n]+1..n)時應進行的最後一次矩陣乘法運算。下面給出的遞迴過程可以輸出<A(i),A(i+1),...,A(j)>的最優括號化方案。
演算法虛擬碼：
PRINT-OPTIMAL-PARENS(s,i,j)
1. if i=j
2. then print "A"i
3. else print "("
4. PRINT-OPTIMAL-PARENS(s,i,s[i,j])
5. PRINT-OPTIMAL-PARENS(s,s[i,j]+1,j)
6. print ")"
C++:
matrixchain.h
#define _matrixchain_h #include <stdlib.h> #include<limits.h> #include <stdio.h> #include"pair.h" pair matrixchainorder(int *p,int n){ int i,l,j,k,q; int *m=(int*)malloc((n+1)*(n+1)*sizeof(int)); int *s=(int*)malloc((n+1)*(n+1)*sizeof(int)); for(i=0;i<=n;i++) m[i*(n+1)+i]=0; for(l=2;l<=n;l++) for(i=1;i<=n-l+1;i++){ j=i+l-1; m[i*(n+1)+j]=INT_MAX; for(k=i;k<=j-1;k++){ q=m[i*(n+1)+k]+m[(k+1)*(n+1)+j]+p[i-1]*p[k]*p[j]; if(q<m[i*(n+1)+j]){ m[i*(n+1)+j]=q; s[i*(n+1)+j]=k; } } } return make_pair(m,s); } void printoptimalparens(int *s,int i,int j,int n){ if(i==j) printf("A%d",i); else{ printf("("); printoptimalparens(s,i,s[i*(n+1)+j],n); printoptimalparens(s,s[i*(n+1)+j]+1,j,n); printf(")"); } }
pair.h
#define _pair_h //定義了一個結構體用來存放生成的表m以及s typedef struct { void *first; void *second; }pair; pair make_pair(void *f,void *d){ pair p={f,d}; return p; }
main.cpp
#include<stdio.h> #include"matrixchain.h" int main(){ int p[]={30,35,15,5,10,20,25};//p表示矩陣的行列數，如30行35列，35行15列，15行5列。。。 int *s,*m; pair r=matrixchainorder(p,6); m=(int*)r.first; s=(int*)r.second; printoptimalparens(s,1,6,6); printf("\n%d\n",m[13]);//之所以m的下標是13的原因是，前面定義的i,j是加在一起的，[i*(n+1)+j]其中i為1，n為元素個數，j為6個矩陣鏈 free(s);//free()釋放calloc, malloc, realloc動態分配的空間 free(m); }
輸出結果為：
C++：
Matrixchaincpp.h
#define _Matrixchaincpp_h #include<iostream> using namespace std; pair<int*,int*> matrixchainorder(int *p,int n){ long i,l,j,k,q; int *m=new int[(n+1)*(n+1)]; int *s=new int[(n+1)*(n+1)]; for(i=0;i<=n;i++) m[i*(n+1)+i]=0; for(l=2;l<=n;l++) for(i=1;i<=n-l+1;i++){ j=i+l-1; m[i*(n+1)+j]=numeric_limits<int>::max(); for(k=i;k<=j-1;k++){ q=m[i*(n+1)+k]+m[(k+1)*(n+1)+j]+p[i-1]*p[k]*p[j]; if(q<m[i*(n+1)+j]){ m[i*(n+1)+j]=q; s[i*(n+1)+j]=k; } } } return make_pair(m,s); } void printoptimalparens(int *s,int i,int j,int n){ if (i==j) cout<<'A'<<i; else{ cout<<'('; printoptimalparens(s,i,s[i*(n+1)+j],n); printoptimalparens(s,s[i*(n+1)+j]+1,j,n); cout<<')'; } }
main.cpp
#include"Matrixchaincpp.h" int main(){ int p[]={30,35,15,5,10,20,25};//p表示矩陣的行列數，如30行35列，35行15列，15行5列。。。 int *s,*m; pair<int*,int*> r=matrixchainorder(p,6);//C++中預設預定義了資料結構類pair，可以直接使用 m=r.first; s=r.second; printoptimalparens(s,1,6,6); cout<<endl<<m[13]<<endl;//之所以m的下標是13的原因是，前面定義的i,j是加在一起的，[i*(n+1)+j]其中i為1，n為元素個數，j為6個矩陣鏈 delete []s;//free()釋放calloc, malloc, realloc動態分配的空間 delete []m; }
輸出結果如上圖一樣！
JAVA:
Matrixchain.java
package Jamin;
public class Matrixchain { public static Pair matrixchainorder(int[] p) {//p為輸入的矩陣行列數鏈 int n=p.length-1,i,l,j,k,q; int [][] m=new int[n+1][n+1],s=new int[n+1][n+1]; for(i=0;i<=n;i++) m[i][i]=0; for(l=2;l<=n;l++) for(i=1;i<=n-l+1;i++) { j=i+l-1; m[i][j]=Integer.MAX_VALUE; for(k=i;k<=j-1;k++) { q=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j]; if(q<m[i][j]) { m[i][j]=q; s[i][j]=k; } } } return Pair.make(m, s); } public static void printoptimalparens(int[][] s,int i,int j) { if(i==j) System.out.printf("A%d",i); else { System.out.print("("); printoptimalparens(s,i,s[i][j]); printoptimalparens(s,s[i][j]+1,j); System.out.print(")"); } } }
Pair.java
package Jamin;
public class Pair {//用於對m表和s表進行儲存的結構變數 public Object first; public Object second; public Pair() { first=new Object(); second=new Object(); } public static Pair make(Object a,Object b) { Pair p=new Pair(); p.first=a; p.second=b; return p; } }
Test.java
package Jamin;
public class Test {
public static void main(String[] args) { // TODO Auto-generated method stub int p[]= {30,35,15,5,10,20,25}; int[][] m,s; Pair r=Matrixchain.matrixchainorder(p); m=(int[][]) r.first; s=(int[][]) r.second; Matrixchain.printoptimalparens(s, 1, 6); System.out.println(); System.out.println(m[1][6]); }
} 輸出結果：
((A1(A2A3))((A4A5)A6)) 15125