動態規划---例題1.矩陣連乘問題
一.問題描述
矩陣A和B可乘的條件是矩陣A的列數等於矩陣B的行數.若A是一個p×q的矩陣，B是一個q×r的矩陣，則其乘積C=AB是一個p×r的矩陣.其標准計算公式為：
計算C=AB總共需要pqr次的數乘.
給定n個矩陣{A1,A2,…,An}.其中Ai與Ai+1是可乘的，i=1,2,…,n-1.要求計算出這n個矩陣的連乘積A1A2…An.
二.解題思路
矩陣乘法滿足結合律，故連乘積的計算可以有許多不同的計算次序.這種計算次序可以用加括號的方式來確定.若一個矩陣連乘積的計算次序已完全確定，也就是說該連乘積已完全加括號，則我們可以通過反復調用兩個矩陣相乘的標准算法計算出矩陣連乘積.完全加括號的矩陣連乘積可遞歸地定義為：
1.單個矩陣是完全加括號的；
2.若矩陣連乘積A是完全加括號的，則A可表示為兩個完全加括號的矩陣連乘積B和C的乘積並加括號，即A=(BC).
例如，矩陣連乘積A1A2A3 A4可以有以下5種不同的完全加括號方式：
- (A1(A2(A3A4)))
- (A1((A2A3)A4))
- ((A1A2)(A3A4))
- ((A1(A2A3))A4)
- (((A1A2)A3)A4)
每一種完全加括號方式對應於一種矩陣連乘積的計算次序，而這種計算次序與計算矩陣連乘積的計算量有着密切的關系.
為了說明在計算矩陣連乘積時加括號方式對整個計算量的影響，我們來看一個計算3個矩陣{A1，A2，A3}的連乘積的例子.
設這3個矩陣的維數分別為10×100，100×5和5×50.若按第一種加括號方式((A1A2)A3)來計算，總共需要
10×100×5+10×5×50=7500次的數乘.
若按第二種加括號方式(A1(A2A3))來計算，則需要的數乘次數為
100×5×50+10×100×50=75000
后者的計算量是前者的10倍.可見，在計算矩陣連乘積時，計算次序對計算量影響很大.§矩陣連乘積的最優計算次序問題，即對於給定的相繼n個矩陣{A1,A2,…,An}(其中Ai的維數為pi-1×p****i ，i=1,2,…,n)，如何確定計算矩陣連乘積A1A2…An的一個計算次序(完全加括號方式)，使得依此次序計算矩陣連乘積需要的數乘次數最少.
解這個問題的最容易想到的方法是窮舉搜索法.也就是列出所有可能的計算次序，並計算出每一種計算次序相應需要的計算量，然后找出最小者.
然而，這樣做計算量太大.
事實上，對於n個矩陣的連乘積，設有P(n)個不同的計算次序.由於可以首先在第k個和第k+1個矩陣之間將原矩陣序列分為兩個矩陣子序列，k=1,2,…,n-1；然后分別對這兩個矩陣子序列完全加括號；最后對所得的結果加括號，得到原矩陣序列的一種完全加括號方式.所以關於P(n)，我們有遞推式如下：
解此遞歸方程可得,P(n)實際上是Catalan數，
即P(n)=C(n-1)，其中，
也就是說，P(n)隨着n的增長是指數增長的.因此，窮舉搜索法不是一個有效算法.
下面我們來考慮用動態規划法解矩陣連乘積的最優計算次序問題.此問題是動態規划的典型應用之一.
1.分析最優解的結構
將矩陣連乘積AiAi+1…Aj簡記為Aij.先看計算A1 n的一個最優次序.設這個計算次序在矩陣Ak和Ak+1之間將矩陣鏈斷開，1<=k<n，則完全加括號方式為((A1…Ak)(Ak+1…An)).照此，我們要先計算A1…k和Ak+1…n，然后，將所得的結果相乘才得到A1 n.顯然其總計算量為計算A1…k的計算量加上計算Ak+1…n的計算量，再加上A1…k與Ak+1…n相乘的計算量.
問題的關鍵特征：計算A1 n的一個最優次序所包含的計算A1…k的次序也是最優的.事實上，若有一個計算A1…k的次序需要的計算量更少，則用此次序替換原來計算A1…k的次序，得到的計算A1…n的次序需要的計算量將比最優次序所需計算量更少，這是一個矛盾.同理，計算A1…n的一個最優次序所包含的計算矩陣子鏈Ak+1…n的次序也是最優的.
最優解包含其子問題的最優解，這種性質稱為最優子結構性質.另外，該問題顯然滿足無后向性，因為前面的矩陣鏈的計算方法和后面的矩陣鏈的計算方法無關.
2.建立遞歸關系
對於矩陣連乘積的最優計算次序問題，設計算Ai…j ,1≤i≤j≤n，所需的最少數乘次數為m[i,j]，原問題的最優值為m[1,n].
- 當i=j時，Ai…j=Ai為單一矩陣，無需計算，因此m[i,i]=0，i=1,2,…,n ；
- 當i<j時，可利用最優子結構性質來計算m[i,j].事實上，若計算Ai…j的最優次序在Ak和Ak+1之間斷開，i≤k<j，則：m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj
斷開點A的位置還未確定.不過k的位置只有j-i個可能，即k∈{i,i+1,…,j-1}.因此k是這j-i個位置中計算量達到最小的那個位置.從而m[i,j]可遞歸地定義為：
m[i,j]給出了最優值，即計算Ai…j所需的最少數乘次數.同時還確定了計算Ai…j的最優次序中的斷開位置k，也就是說，對於這個k有m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj .若將對應於m[i,j]的斷開位置k記錄在s[i,j]中，則相應的最優解便可遞歸地構造出來.
三.計算遞歸值
根據m[i,j]的遞歸定義(2.1)，容易寫一個遞歸程序來計算m[1,n].然而簡單地遞歸計算將耗費指數計算時間.注意：在遞歸計算過程中，不同的子問題個數只有θ(n2)個.事實上，對於1≤i≤j≤n不同的有序對(i,j)對應於不同的子問題.因此，不同子問題的個數最多只有
個.可見，許多子問題被重復計算多次.
用動態規划算法解此問題，可依據遞歸式(2.1)以自底向上的方式進行計算，在計算過程中，保存已解決的子問題答案，每個子問題只計算一次，而在后面需要時只要簡單查一下，從而避免大量的重復計算，最終得到多項式時間的算法.
下面所給出的計算m[i,j]動態規划算法中，輸入是序列P={p0,p1,…,pn}，輸出除了最優值m[i,j]外，還有使
m[i,j] = m[i,k] + m[k+1,j] + pi-1pkpj
達到最優的斷開位置k=s[i,j]，1≤i≤j≤n;
該算法按照
- m[1,1]
- m[2,2] m[1,2]
- m[3,3] m[2,3] m[1,3]
- ... ... ...
- m[n,n] m[n-1,n] ... .... ... m[1,n]
的順序根據公式(2.1)計算m[i, j].
該算法的計算時間上界為O(n3).算法所占用的空間顯然為O(n2).由此可見，動態規划算法比窮舉搜索法要有效得多.
4.構建最優解
算法MatrixChain只是計算出了最優值，並未給出最優解.也就是說，通過MatrixChain的計算，我們只知道計算給定的矩陣連乘積所需的最少數乘次數，還不知道具體應按什么次序來做矩陣乘法才能達到數乘次數最少.
然而，它己記錄了構造一個最優解所需要的全部信息.
因為s[i, j]中的數k表明計算矩陣鏈Ai…j的最佳方式應在矩陣Ak和Ak+1之間斷開，即最優的加括號方式應為(A1...k)(Ak+1…n).因此，從s[i,j]記錄的信息可知計算A1…n的最優加括號方式為 (A1…s[1,n])(As[1,n]+1…n).
而計算A1…s[1,n]的最優加括號方式為(A1…s[1,s[1,n]])(As[1,s[1,n]]+1…s[1,n]).同理可以確定計算As[1,n]+1…n的最優加括號方式在s[s[1,n]+1,n]處斷開.…照此遞推下去，最終可以確定As[1,n]+1…n的最優完全加括號方式，即構造出問題的一個最優解.
代碼如下:
// 矩陣乘法鏈動態規划算法
// 矩陣連乘,動態規划迭代實現(自底向上)
// A1 30*35 A2 35*15 A3 15*5 A4 5*10 A5 10*20 A6 20*25
// p[0-6] = {30,35,15,5,10,20,25}
#include<bits/stdc++.h>
using namespace std;
const int L = 7;
int MatrixChain(int *p, int n, int ** m, int **s);
void Traceback(int i, int j, int **s); //構造最優解
int main()
{
int p[L] = {30,35,15,5,10,20,25};
int **s = new int *[L];
int **m = new int *[L];
for(int i=0; i<L; i++)
{
s[i] = new int[L];
m[i] = new int[L];
}
cout<<"矩陣的最少計算次數為"<<MatrixChain(p, 6, m, s)<<endl;
cout<<"矩陣最優計算次序為:"<<endl;
Traceback(1, 6, s);
system("pause");
return 0;
}
int MatrixChain(int *p, int n, int ** m, int ** s) //T(n) = O(n^3)
{
for(int i=1; i<=n; i++) m[i][i] = 0;
for(int r=2; r<=n; r++) //r個矩陣連乘.r為當前計算的鏈長(子問題規模)
for(int i=1; i<=n-r+1; i++) //n-r+1為最后一個r鏈的前邊界
{
int j=i+r-1; //計算前邊界為i,鏈長為r的后邊界
m[i][j] = m[i+1][j] + p[i-1]*p[i]*p[j]; //將鏈ij划分為A(i)*(A[i+1:j]). ]m[i][i]=0省略.
s[i][j] = i;
for(int k=i+1; k<j; k++) //分割位置從i+1不斷向后移動,將鏈ij划分為(A[i:k]*A[k+1:j])
{
int t = m[i][k] +m[k+1][j] + p[i-1]*p[k]*p[j];
if(t < m[i][j])
{
m[i][j] = t;
s[i][j] = k;
}
}
}
return m[1][L-1];
}
void Traceback(int i, int j, int **s)
{
if(i==j) return ;
Traceback(i, s[i][j], s);
Traceback(s[i][j]+1, j, s);
cout<<"Mutiply A"<<i<<","<<s[i][j];
cout<<" and A"<<(s[i][j]+1)<<","<<j<<endl;
}
運行結果:
參考我的老師畢方明《算法設計與分析》課件.
如果覺得本篇文章對你有所幫助,歡迎大家來到我的個人博客網站---喬治的編程小屋逛一逛吧.