【動態規劃】矩陣鏈乘法
阿新 • • 發佈：2019-02-02
矩陣鏈乘法
求解矩陣鏈相乘問題時動態規劃演算法的另一個例子。給定一個n個矩陣的序列（矩陣鏈）<A1,A2,...,An>，我們希望計算它們的乘積 A1A2...An
為了計算表示式，我們可以先用括號明確計算次序，然後利用標準的矩陣相乘演算法進行計算。完全括號化(fully parenthesized)：它是單一矩陣，或者是兩個完全括號化的矩陣乘積鏈的積。
例如如果有矩陣鏈為<A1,A2,A3,A4>,則共有5種完全括號化的矩陣乘積鏈。
(A1(A2(A3A4)))、(A1((A2A3)A4))、((A1A2)(A3A4))、((A1(A2A3))A4)、((A1(A2A3))A4)
對矩陣鏈加括號的方式會對乘積運算的代價產生巨大影響。我們先來分析兩個矩陣相乘的代價。下面的虛擬碼的給出了兩個矩陣相乘的標準演算法，屬性rows和columns是矩陣的行數和列數。
以矩陣鏈<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所需標量乘法次數最少。
因為括號方案的數量與n呈指數關係，所以通過暴力搜尋窮盡所有可能的括號化方案來尋找最優方案是一個糟糕策略。
應用動態規劃方法
下面用動態規劃方法來求解矩陣鏈的最優括號方案，我們還是按照之前提出的4個步驟進行：
1.刻畫一個最優解的結構特徵
2.遞迴地定義最優解的值
3.計算最優解的值，通常採用自底向上的方法
4.利用計算出的資訊構造一個最優解
接下來按順序進行這幾個步驟，清楚地展示針對本問題每個步驟應如何做。
步驟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
步驟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構造最優解。
下圖展示了對一個長度為6的矩陣鏈執行此演算法的過程。由於定義m[i,j]僅在i<=j時有意義，因此表m只使用主對角線之上的部分，圖中的表是經過旋轉的，主對角線已經旋轉到水平方向，MATRIX_CHAIN_ORDER按自下而上、自左至右的順序計算所有行。當計算表m[i,j]時，會用到乘積p(i-1)p(k)p(j)(k=i,i+1,...j-1)，語句m[i,j]西南方向和東南方向上所有表項。
步驟4：構造最優解
求解矩陣鏈相乘問題時動態規劃演算法的另一個例子。給定一個n個矩陣的序列（矩陣鏈）<A1,A2,...,An>，我們希望計算它們的乘積 A1A2...An
為了計算表示式，我們可以先用括號明確計算次序，然後利用標準的矩陣相乘演算法進行計算。完全括號化(fully parenthesized)：它是單一矩陣，或者是兩個完全括號化的矩陣乘積鏈的積。
例如如果有矩陣鏈為<A1,A2,A3,A4>,則共有5種完全括號化的矩陣乘積鏈。
(A1(A2(A3A4)))、(A1((A2A3)A4))、((A1A2)(A3A4))、((A1(A2A3))A4)、((A1(A2A3))A4)
對矩陣鏈加括號的方式會對乘積運算的代價產生巨大影響。我們先來分析兩個矩陣相乘的代價。下面的虛擬碼的給出了兩個矩陣相乘的標準演算法，屬性rows和columns是矩陣的行數和列數。
兩個矩陣A和B只有相容(compatible)，即A的列數等於B的行數時，才能相乘。如果A是p×q的矩陣，B是q×r的矩陣，那麼乘積C是p×r的矩陣。計算C所需要時間由第8行的標量乘法的次數決定的，即pqr。MATRIX-MULTIPKLY(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 c(ij)=0 for k = 1 to A.columns c(ij)=c(ij)+a(ik)*b(kj) return C
以矩陣鏈<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所需標量乘法次數最少。
因為括號方案的數量與n呈指數關係，所以通過暴力搜尋窮盡所有可能的括號化方案來尋找最優方案是一個糟糕策略。
應用動態規劃方法
下面用動態規劃方法來求解矩陣鏈的最優括號方案，我們還是按照之前提出的4個步驟進行：
1.刻畫一個最優解的結構特徵
2.遞迴地定義最優解的值
3.計算最優解的值，通常採用自底向上的方法
4.利用計算出的資訊構造一個最優解
接下來按順序進行這幾個步驟，清楚地展示針對本問題每個步驟應如何做。
步驟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。(C++實現)
void MATRIX_CHAIN_ORDER(int *p,int Length,int m[][M],int s[][M])
{
int q,n=Length-1;
for(int i=1;i<=n;i++) m[i][i]=0;
for(int l=2;l<=n;l++) /* 矩陣鏈的長度 */
{
for(int i=1;i<=n-l+1;i++)
{
int j=i+l-1; /* 等價於 l=j-i+1 */
m[i][j]=INT_MAX;
for(int 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;
}
}
}
}
}下圖展示了對一個長度為6的矩陣鏈執行此演算法的過程。由於定義m[i,j]僅在i<=j時有意義，因此表m只使用主對角線之上的部分，圖中的表是經過旋轉的，主對角線已經旋轉到水平方向，MATRIX_CHAIN_ORDER按自下而上、自左至右的順序計算所有行。當計算表m[i,j]時，會用到乘積p(i-1)p(k)p(j)(k=i,i+1,...j-1)，語句m[i,j]西南方向和東南方向上所有表項。
n=6和矩陣規模如下表時，MATRIX_CHAIN_ORDER計算出m表和s表
6個矩陣相乘所需的最少標量乘法運算次數為m[1,6]=15125。表中有些表項被標記了深色陰影，相同的陰影表示在第13行中計算m[2,5]時同時訪問了這些表項：
步驟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)>的最優括號化方案。
void PRINT_OPTIMAL_PARENS(int s[][M],int i,int j)
{
if(i == j) cout<<"A"<<i;
else
{
cout<<"(";
PRINT_OPTIMAL_PARENS(s,i,s[i][j]);
PRINT_OPTIMAL_PARENS(s,s[i][j]+1,j);
cout<<")";
}
}對上圖中，n=6時，呼叫PRINT_OPTIMAL_PARENS(s,1,6)輸出括號化方案 ((A1(A2A3))((A4A5)A6))
下面給出完整的程式碼，當n=6時的最優解和括號方案。
#include<iostream>
using namespace std;
const int INT_MAX=2147483647;
int const M=7;
void MATRIX_CHAIN_ORDER(int *p,int Length,int m[][M],int s[][M])
{
int q,n=Length-1;
for(int i=1;i<=n;i++) m[i][i]=0;
for(int l=2;l<=n;l++) /* 矩陣鏈的長度 */
{
for(int i=1;i<=n-l+1;i++)
{
int j=i+l-1; /* 等價於 l=j-i+1 */
m[i][j]=INT_MAX;
for(int 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;
}
}
}
}
}
void PRINT_OPTIMAL_PARENS(int s[][M],int i,int j)
{
if(i == j) cout<<"A"<<i;
else
{
cout<<"(";
PRINT_OPTIMAL_PARENS(s,i,s[i][j]);
PRINT_OPTIMAL_PARENS(s,s[i][j]+1,j);
cout<<")";
}
}
int main()
{
int p[M]={30,35,15,5,10,20,25};
int m[M][M],s[M][M];
MATRIX_CHAIN_ORDER(p,M,m,s);
cout<<"當n=6時最優解為: \n"<<m[1][6];
cout<<"\n括號化方案為：\n";
PRINT_OPTIMAL_PARENS(s,1,6);
return 0;
}