动态规划之矩阵连乘

以下内容参考（摘抄）《算法设计与分析》，王晓东编著，清华大学出版社2003年1月第1版。

给定n个矩阵{A₁,A₂,…,A_n}，其中A_i与A_i+1是可乘的，i=1,2,…,n-1。考察这n个矩阵的连乘积A₁A₂…A_n。由于矩阵乘法满足结合律，故计算矩阵的连乘积可以有许多不同的计算次序，这种计算次序可以用加括号的方式来确定。若一个矩阵连乘积的计算次序完全确定，则可以依此次序反复调用2个矩阵相乘的标准算法（有改进的方法，这里不考虑）计算出矩阵连乘积。若A是一个p×q矩阵，B是一个q×r矩阵，则计算其乘积C=AB的标准算法中，需要进行pqr次数乘。

矩阵连乘积的计算次序不同，计算量也不同，举例如下：

先考察3个矩阵{A₁,A₂,A₃}连乘，设这三个矩阵的维数分别为10×100，100×5，5×50。若按（（A₁A₂）A₃）方式需要的数乘次数为10×100×5＋10×5×50＝7500，若按（A₁（A₂A₃））方式需要的数乘次数为100×5×50＋10×100×50＝75000。

下面使用动态规划法找出矩阵连乘积的最优计算次序。

1，  设矩阵连乘积A_iA_i+1…A_j简记为A[i:j]，设最优计算次序在A_k和A_k+1之间断开，则加括号方式为：

（（A_iA_i+1…A_k）（A_k+1…A_j））

则依照这个次序，先计算A[i:k]和A[K+1:j]然后再将计算结果相乘，计算量是：

A[i:k]的计算量加上A[K+1:j]的计算量再加上它们相乘的计算量。

问题的一个关键是：计算A[i:j]的最优次序所包含的两个子过程（计算A[i:k]和A[K+1:j]）也是最优次序。

2，  设计算A[i:j]所需的最少数乘次数为m[i][j]。

i=j时为单一矩阵，则m[i][i]=0，

i<j时，设最优计算次序在A_k和A_k+1之间断开，则m[i][j]=m[i][k]+m[k+1][j]+p_ip_k+1p_j+1，其中p表示数组的维数，例如A₀到A₅共6个数组（为了C语言的描述方便，下标从0开始），他们表示如下：

//p[0]:第一个矩阵的行数

    //p[1]:第一个矩阵的列数，第二个矩阵的行数

    //p[2]:第二个矩阵的列数，第三个矩阵的行数

k此时并未确定，需要从i到j-1遍历以寻找一个最小的m[i][j]。我们把这个最小的k放在s[i][j]。

以下是完整实现代码，以一个具体的例子实现，稍加修改即可通用。

 

#include<iostream>

usingnamespace std;

//////////////////////////////////////////////////////////////////////////////

//MatrixChain计算m[i][j]所需的最少数乘次数

//并记录断开位置s[i][j]

//////////////////////////////////////////////////////////////////////////////

void MatrixChain(int *p,int n,int **m,int **s)

{

     for(int i=0;i<n;i++)

         m[i][i]=0;//单个矩阵相乘，所需数乘次数为0

 

     //以下两个循环是关键之一，以6个矩阵为例(为描述方便，m[i][j]用ij代替)

     //需按照如下次序计算

     //01 12 23 34 45

     //02 13 24 35

     //03 14 25

     //04 15

     //05

     //下面行的计算结果将会直接用到上面的结果。例如要计算14，就会用到12，24；或者13，34等等

     for(int r=1;r<n;r++)

     {

         for(int i=0;i<n-r;i++)

         {

              int j=i+r;

              //首先在i断开，即（Ai*（Ai+1...Aj））

              m[i][j]=m[i][i]+m[i+1][j]+p[i]*p[i+1]*p[j+1];

              s[i][j]=i;

              for(int k=i+1;k<j;k++)

              {

                   //然后在k（从i+1开始遍历到j-1）断开，即（（Ai...Ak)*（Ak+1...Aj））

                   int t=m[i][k]+m[k+1][j]+p[i]*p[k+1]*p[j+1];

                   if(t<m[i][j])//找到更好的断开方法

                   {

                       m[i][j]=t;//记录最少数乘次数

                       s[i][j]=k;//记录断开位置

                   }

              }

         }

     }

     //如果使用下面注释的循环，则是按照如下次序计算

     //01 02 03 04 05

     //12 13 14 15

     //23 24 25

     //34 35

     //45

     //当要计算时14，会用到12，24，而此时24并没有被计算出来。

/*

     for(int i=0;i<n;i++)

     {

         for( int j=i+1;j<n;j++)

         {

              m[i][j]=m[i][i]+m[i+1][j]+p[i]*p[i+1]*p[j+1];

              s[i][j]=i;

              for(int k=i+1;k<j;k++)

              {

                   int t=m[i][k]+m[k+1][j]+p[i]*p[k+1]*p[j+1];

                   if(t<m[i][j])

                   {

                       m[i][j]=t;

                       s[i][j]=k;

                   }

              }

         }

     }

     */

}

//////////////////////////////////////////////////////////////////////////////

//Traceback打印A[i:j]的加括号方式

//////////////////////////////////////////////////////////////////////////////

void Traceback(int i,int j,int **s)

{ 

     //s[i][j]记录了断开的位置，即计算A[i:j]的加括号方式为：

     //(A[i:s[i][j]])*(A[s[i][j]+1:j])

     if(i==j)return;

     Traceback(i,s[i][j],s);//递归打印A[i:s[i][j]]的加括号方式

     Traceback(s[i][j]+1,j,s);//递归打印A[s[i][j]+1:j]的加括号方式

 

     //能走到这里说明i等于s[i][j]，s[i][j]+1等于j

     //也就是说这里其实只剩下两个矩阵，不必再分了

     cout<<"A"<<i<<"和A"<<(s[i][j]+1)<<"相乘"<<endl;    

}

 

int _tmain(int argc, _TCHAR* argv[])

{

     int n=6;//矩阵的个数

     int *p=newint[n+1];

     //p[0]:第一个矩阵的行数

     //p[1]:第一个矩阵的列数，第二个矩阵的行数

     //p[2]:第二个矩阵的列数，第三个矩阵的行数

     p[0]=30;

     p[1]=35;

     p[2]=15;

     p[3]=5;

     p[4]=10;

     p[5]=20;

     p[6]=25;

 

     int **m,**s;

     m=newint*[n];

     for( int i=0;i<n;i++)

         m[i]=newint[n];

 

     s=newint*[n];

     for(int i=0;i<n;i++)

         s[i]=newint[n];   

 

     MatrixChain(p,n,m,s);

     Traceback(0,n-1,s);

for(int i=0;i<n;i++)   

         {   

              delete []m[i];

              m[i]=NULL; 

              delete []s[i];

              s[i]=NULL; 

         }   

         delete []m;   

         m=NULL;  

         delete []s;   

         s = NULL; 

         delete []p;   

         p = NULL;  

     return 0;

}

 

打印结果是：

A1和A2相乘

A0和A1相乘

A3和A4相乘

A3和A5相乘

A0和A3相乘

实际上要表达的是如下加括号方式：

（（A0（A1A2））（（A3A4）A5））

加了括号之后用第一个来代替，例如（A1A2）可看作A1，这个结果的数乘次数是15125。