矩陣連乘

矩陣連乘詳解

                               --crystal yi

  既然這篇文章叫作矩陣連乘詳解，那麼我就不能辜負詳解這兩個字，只有把一個原來不懂的的人弄懂了，才叫詳解。

 

  言歸正傳，首先讓咱們複習一下矩陣連乘的有關知識。對於矩陣知識很瞭解的人能夠跳過矩陣知識這塊內容，不過筆者建議最好複習一下。

 

矩陣知識:

矩陣的乘法：

 

 左面的矩陣的行數決定結果的行數，右面的矩陣的列數決定結果的行數，兩個矩陣在相乘時必定要有相同的行數或者列數。

矩陣相乘只有在第一個矩陣的列數和第二個矩陣的行數相同時纔有定義。假如A爲m×n矩陣，B爲n×p矩陣，則他們的乘AB(有時記作A·B)會是一個m×p矩陣。

例如一個2x3的矩陣與一個3x2的矩陣的乘會是一個2x2的矩陣 。

 

例子：

 

 

矩陣連乘：

 

設有矩陣M1,M2,M3,M4，
其維數分別是10×20, 20×50, 50×1 和1×100，現要求出這4個矩陣相乘的結果。咱們知道，若矩陣A的維數是p×q，矩陣B的維數是q×r，則A與B相乘後所得矩陣AB的維數是p×r。按照矩陣相乘的定義，求出矩陣AB中的一個元素須要作q次乘法（及q-1次加法）。這樣，要計算出AB就須要作p×q×r次乘法。爲簡單起見，且因爲加法比一樣數量的乘法所用時間要少得多，故這裏咱們暫不考慮加法的計算量。因爲矩陣連乘知足結合律，故計算矩陣連乘的方式能夠有多種。

 

例如，咱們能夠按M1(M2(M3M4))的方式去計算，

也能夠按(M1(M2M3))M4的方式去計算，所得結果是相同的。

可是值得注意的是，

按前一方式計算須要作125,000次乘法，

而按後一方式計算只須要作2,200次乘法。

因而可知，矩陣連乘的運算次序對於所須要的計算量

（所需乘法次數）有着極大的影響。

M3M4：50*1*100=5,000；M2(M3M4)：20*50*100=100,000

M1(M2(M3M4))：10*20*100=20,000

(M2M3)：20*50*1=1000；(M1(M2M3))：10*20*1=200 ；

(M1(M2M3))M4：10*1*100=1000

 

如何解決矩陣連乘問題：

 

分析：

 

設要求出矩陣連乘MiMi+1……Mj-1Mj（i<j）所需的最少乘法次數。
因共有j-i+1個矩陣，故稱這個矩陣連乘的規模是j-i+1

 

按照作最後一次乘法的位置進行劃分，該矩陣連乘一共可分爲j-i種狀況即有(j-i)種斷開方式：Mi(Mi+1┅Mj)，(MiMi+1)(Mi+2┅Mj)，┅，(MiMi+1┅Mj-1)Mj。其中任一對括號內的矩陣個數（即規模）不超過j-i。若咱們已知任一個規模不超過j-i的矩陣連乘所需的最少乘法次數，咱們就能夠很容易地計算出矩陣連乘MiMi+1┅Mj-1Mj（i<j）所需的最少乘法次數，其方法以下。將上述的j-i種狀況表示爲通式：(Mi┅Mk) (Mk+1┅Mj)（i≤k<j）。

 

記第t個矩陣Mt的列數爲rt，並令rt-1爲矩陣Mt的行數。
則Mi┅Mk連乘所得是ri-1×rk維矩陣，
Mk+1┅Mj連乘所得是rk×rj維矩陣，
故這兩個矩陣相乘須要作ri-1×rk×rj次乘法

 

因爲在此以前咱們已知
任一個規模不超過j-i的矩陣連乘所需的最少乘法次數，故(Mi┅Mk)和(Mk+1┅Mj)所需的最少乘法次數已知，將它們分別記之爲mi,k和mk+1,j。
形爲(Mi┅Mk)(Mk+1┅Mj)的矩陣連乘所需的最少乘法次數爲：

mi,k   +   mk+1,j   +   ri-1×rk×rj。

 

對知足i≤k<j 的共j-i種狀況逐一進行比較，咱們就能夠獲得
矩陣連乘MiMi+1┅Mj-1Mj（i<j）所需的最少乘法次數mi,j爲：
mi,j=min {   mi,k  +  mk+1,j   +   ri-1×rk×rj   }     （i≤k<j）

 

因而在初始時咱們定義mi,i=0（至關於單個矩陣的狀況）

求m1,2：即i=1, j=2，

就是2個矩陣，無需劃分 (k=1,由於i<=k<j)

m1,2=min{ m1,1  +  m2,2  +  ri-1×ri×ri+1

=ri-1×ri×ri+1=r0×r1×r2= 10×20×50=10000，
求m23：即i=2, j=3，故ri-1×ri×ri+1=r1×r2×r3=20×50×1=1000，

求m13：即i=1, j=3，min（i≤k<j）{mik+mk+1,j+ri-1×rk×rj}=
min{ (m11+m23+r0×r1×r3)(k=1)，(m12+m33+r0×r2×r3)(k=2)}=
min{(0+1000+200), (10000+0+500)}= min{1200,10500}=1200

 

 

矩陣在C語言中的表示方法：

 

A[i][j]     (i爲行，j爲列)

 

 

擁有了上面的準備知識看下面的內容就會輕鬆不少。

 

注：下文所用的mi,j或者m[i][j]表示MiMi+1┅Mj-1Mj（i<j）所需的最少乘法次數

 

首先，讓咱們本身思考一下怎麼解決一個n個矩陣連乘的問題，也許有了上面的知識你會想到下面的遞歸的方法：

 

假設已經算好了Mi┅Mk和Mk+1┅Mj這兩個子矩陣各自相乘的最小的次數，剩下的只須要算最後一次乘積的次數就能夠了，也就是這裏的兩個子矩陣相乘所涉及的乘法次數ri-1×rk×rj而後加上前面所得的兩個子矩陣各自所需的最小的相乘次數就能夠了。這樣就獲得了上面所說的這個式子：mi,j = mi,k  +  mk+1,j   +   ri-1×rk×rj  而後依次判斷在k等於某個數的時候所得的mi,j的值，剩下的只須要對全部的mi,j的值取一個最小值就ok了，因而咱們就獲得了上面所說的這個式子：  mi,j=min {   mi,k  +  mk+1,j   +   ri-1×rk×rj   }     （i≤k<j）

然而這樣作並不能求出mi,j的最小值，由於咱們不知道 mi,k 和 mk+1,j的值。可是，這並不會難倒咱們，由於咱們知道對於 mi,k 和 mk+1,j 的求解能夠基於遞歸思想。方法和求解mi,j的最小值同樣。這樣咱們就獲得了遞歸的式子：

 

下面代碼的p[i]表示的是第i個矩陣的列數，固然他的行數就用i-1表示，呵呵。

 1 下面代碼的p[i]表示的是第i個矩陣的列數，固然他的行數就用i-1表示，呵呵。
 2 
 3 int min(int a,intb){
 4 
 5     return a<b?a:b;
 6 
 7 }
 8 
 9 int d(int i,int j){
10 
11     int i,j,k,t,u;
12 
13     if(i==j)   return 0;  //i==j時就是mi,i，也就是單個矩陣，固然，單個矩陣不存在乘積因此爲0
14 
15     if(i==j-1)   return p[i-1]*p[i]*p[j];   //明顯，這是相鄰的兩個矩陣的乘積
16 
17      u=d(i,i)+d(i+1,j)+p[i-1]*p[i]*p[j];   //1
18 
19     for(k=i+1;k<j;k++){                   //2
20 
21       t=d(i,k)+d(k+1,j)+p[i-1]*p[k]*p[j];  //2
22 
23         u=min(t,u);                      //3
24 
25     }
26 
27 }

1           對於從矩陣i到矩陣j的連乘，當k==i時所得的mi,j

2           對於從矩陣i到矩陣j的連乘，當i<k<j時所得的 mi,j

注：1 ， 2 所說的結合起來就是 3 也就是上面藍字部分的分析

3   就是mi,j=min {   mi,k  +  mk+1,j   +   ri-1×rk×rj   }     （i≤k<j）

 

 

 

但是，有興趣的話你能夠算一下用遞歸作的效率，你會發現，複雜度是O(2^n)，這個複雜度太大了，稍大一點的數據就爆了不只timelimitexceeded並且也會memorylimitexceeded，因此咱們就必須換一種方法改用dp(動態規劃)，可是，咱們的思路仍是同樣的，換句話說，只不過把遞歸的形式改成dp的形式而已，沒什麼新的東西，若是說有的話也只是記錄了中間狀態。請看下面的代碼：

 1 int dp_matrix_multiply(){
 2 
 3     int i,j,r,k,tem;
 4 
 5     n-=1;
 6 
 7     for(i=1;i<=n;i++)    //1
 8 
 9         m[i][i]=0;      //1
10 
11     for(r=2;r<=n;r++){    //2
12 
13         for(i=1;i<=n-r+1;i++){  //3
14 
15             j=i+r-1;           //3
16 
17             m[i][j]=m[i+1][j]+p[i-1]*p[i]*p[j];   //4
18 
19             for(k=i+1;k<j;k++){                   //4
20 
21                 tem=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];    //4
22 
23                 m[i][j]=tem<m[i][j]?tem:m[i][j];            //4
24 
25             }
26 
27         }
28 
29     }
30 
31     return m[1][n];
32 
33 }

初看，是否是很頭痛？？我剛開始接觸的時候看到別人這麼寫也很頭大，可是思考了近3個小時後終於理解了。

 

首先咱們看一下注釋1，這部應該很好理解吧，就是初始化單個矩陣連乘的最小次數，固然單個矩陣不可能與其自身相乘因此賦值爲0

 

而後是註釋2，這個我想應該是初學者最難以理解的地方之一了吧(包括當時的我也同樣,呵呵)，其實這裏的r就是 the matrix train length，連乘的矩陣個數。這裏從2開始，有兩個緣由：首先單個矩陣的乘積(就是其自身，這一步已經在1中闡明瞭)已經初始化好了，也就是說咱們已經求得了r==1時的若干單個矩陣最小的乘積次數，接下去固然要求r==2時，若干個由兩個矩陣組成的的子矩陣乘積的最小的次數，而後依次類推。其次，對於這裏的r的解釋咱們能夠看一下這個狀態方程：  mi,j=min {   mi,k  +  mk+1,j   +   ri-1×rk×rj   }     （i≤k<j）

若是咱們要知道mi,j的值，咱們就必須知道  mi,k  和  mk+1,j的值，而mi,k和mk+1,j的求解又是經過更小的子模塊來實現，以此類推最後在mi,i，也就是單個矩陣的地方終止子問題的進一步求解，由於沒有比單個矩陣更小的子問題了，講到這裏咱們不妨抽點時間再看一下上面講的遞推的方法，你會發現遞歸的實質就是求解更小的子模塊，這不是和這裏的dp求解子模塊的方法同樣嗎。並且你會發現，註釋4的部分和遞歸求解的2，3部分幾乎如出一轍，那是否是一樣的意思呢，相信本身吧，你理解的是正確的，可是，你要注意一下dp的方法優越於遞歸的方法的主要緣由就是將遞歸中的t=d(i,k)+d(k+1,j)+p[i-1]*p[k]*p[j];改爲了tem=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];這樣用數組記錄狀態節約了大量的時間，由於遞歸方法中d(i,k)這個函數老是被屢次調用，不能用遞歸的方法解決大數據的緣由就在於此。

 

理解了上面的內容註釋3的意思也不難看懂了，i和j就是MiMi+1┅Mj-1Mj（i<j）中的i和j，至於i的值爲何從1到n-r+1以及j的值爲何是i+r-1不用我說你也應該理解，i從1到n-r+1就是r個矩陣連乘的第一個矩陣的下標(注：這裏矩陣的下標從1開始),既然這樣那麼j就是r個矩陣相乘最後那個矩陣的下標。

 

呵呵，經過上面的詳解應該理解了吧。

打鐵趁熱趕快作一下下面的題目：

http://poj.org/problem?id=1651               pku 1651

 1 #include <stdio.h>
 2 
 3 #include <string.h>
 4 
 5 #include <stdlib.h>
 6 
 7 #define len 105
 8 
 9 int m[len][len],p[len];
10 
11 int n;
12 
13 bool input(){
14 
15     int i;
16 
17     if(scanf("%d",&n)==1&&n>0){
18 
19         for(i=0;i<n;i++)
20 
21             scanf("%d",&p[i]);
22 
23         return true;
24 
25     }
26 
27     else return 0;
28 
29 }
30 
31 int dp_matrix_multiply(){
32 
33     int i,j,r,k,tem;
34 
35     n-=1;
36 
37     for(i=1;i<=n;i++)
38 
39         m[i][i]=0;
40 
41     for(r=2;r<=n;r++){
42 
43         for(i=1;i<=n-r+1;i++){
44 
45             j=i+r-1;
46 
47             //i==k
48 
49             m[i][j]=m[i+1][j]+p[i-1]*p[i]*p[j];
50 
51             for(k=i+1;k<j;k++){
52 
53                 tem=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
54 
55                 m[i][j]=tem<m[i][j]?tem:m[i][j];
56 
57             }
58 
59         }
60 
61     }
62 
63     return m[1][n];
64 
65 }
66 
67 int main(void){
68 
69     int ans;
70 
71     while(input()){
72 
73         ans=dp_matrix_multiply();
74 
75         printf("%d\n",ans);
76 
77     }
78 
79     return 0;
80 
81 }