动态规划之矩阵连乘问题

问题描述：

矩阵连乘问题是通过给矩阵连乘时加括号，使得总的计算量最小。

考虑3个矩阵相乘的例子， A1,A2,A3,假设这3个矩阵的维数分别为 10x100,100x50,5x50

若按照（（A1A2）A3）方式计算，需要的数乘次数为10x100x5+10x5x50 = 7500

若按照（A1(A2A3)）方式计算，需要的数乘次数为100x5x50+10x100x50 =75000

问题分析：

1.分析最优解的结构 

    设计求解具体问题的动态规划算法的第一步是刻画该问题的最优解的结构特征。我们将矩阵连乘积AiAi+1....Aj简记为A[ i : j ]。考察计算A[ 1: n]的最优计算次序。设这个计算次序在矩阵Ak和Ak+1之间将矩阵链断开，1<=k<n，则其相应的完全加括号形式为((A1...Ak)(Ak+1...An))。以此次序，总的计算量为A[ 1 : k ]的计算量加上A[ k+1 : n ]的计算量， 再加上A[ 1 : k ]和A[ k+1 : n ]相称的计算量。
    这个问题的关键特征是：计算A[ 1 :n ]的最优次序所包含的计算矩阵子链a[ 1 : k ]和A[ k+1 : n ]的次序也是最优的。因此，矩阵连乘积计算次序问题的最优解包含着其子问题的最优解。这种性质称为最优子结构性质。问题的最优子结构性质是该问题可以用动态规划算法求解的显著特征。

2.建立递归关系

   

  设计动态规划算法的第二步就是递归地定义最优值。对于矩阵连乘积的最有计算次序问题，设计算A[i:j], 1<=i<=j<=n，所需的最少数乘次数为m[i][j]，则原问题的最优值为m[1][n]。

   当i=j时，A[i;j ]=Ai，为单一矩阵，无需计算，因此m[i][i]=0。

   当i < j时，可以利用最优子结构的性质来计算m[i][j]。事实上，若计算A[i:j]的最优次序在Ak和Ak+1之间断开，i<=k<j，则m[i][j]=m[i][k]+m[k+1][j]+Pi-1*Pk*Pj。其中Pi表示第i个矩阵的列数，也是第i-1个矩阵的行数，P0表示第一个矩阵的行数。由于在计算时并不知道断开点k的位置，所以k还未定。不过k的位置只有j-i个可能。从而m[i][j]可以递归地定义为

                      当i=j   m[i][j] = 0

                      当i<j   m[i][j] = min{ m[i][k]+m[k+1][j]+Pi-1*Pk*Pj }

    m[i][j]给出了最优值，即计算A[i:j]所需的最少数乘次数。同时还确定了计算A[i:j]的最优次序中的断开位置k，也就是说，对于这个k有

                       m[i][j]=m[i[k]+m[k+1][j] + Pi-1*Pk*Pj

   若将对应于m[i][j]的断开位置k记为s[i][j]，在计算最优值m[i][j]后，可以递归地有s[i][j]构造出相应的最优解。

 3. 计算最优值

    根据计算m[ i ][ j ]的递归式，容易写一个递归算法计算m[ 1 ][ n ]。但是简单地递归将好费指数计算时间。在递归计算时，许多子问题被重复计算多次。这也是该问题可以用动态规划算法求解的又一显著特征。

   用动态规划算法解决此问题，可依据其递归是以自底向上的方式进行计算。在计算的过程中，保存以解决的子问题答案。每个子问题只计算一次，而在后面需要时只要简单查一下，从而避免大量的重复计算。

4. 代码实现

   下面的算法中，输入参数p0,p1,……p2存储在数组P中。算法除了输出最优值数组外还输出记录最优断开置的数组s.

   算法首先计算出m[i][i] = 0,i=1,2,3……,n,然后再根据递归式，按矩阵链增长的方式依次计算m[i][i+1],i=1,2,……n-1,(矩阵链长度为2）；m[i][i+2],i=1,2,……n-2,(矩阵链长度为3)；在计算m[i][j]时，只用到已经计算出的m[i][k]和m[k+1][j].

#include<iostream>#define LEN 6 // 连乘的矩阵个数using namespace std;//矩阵的维数分别是 30x35,35x15,15x5,5x10,10x20,20x25int p[LEN + 1] ={ 30,35,15,5,10,20,25};int m[LEN + 1][LEN + 1] = { 0,0,0};int s[LEN + 1][LEN + 1] = { 0,0,0};void MatrixChain(int *p,int m[][LEN + 1], int s[][LEN + 1]){    for(int i = 0; i <= LEN; ++i){        m[i][i] = 0;    }//初始化，即链长为1    for(int r = 2; r <= LEN; ++r){//链长从2到LEN        for(int i = 1; i <= LEN - r + 1; ++i){对i位置的矩阵计算链长为r的最优解            int j = i + r - 1;            m[i][j] = m[i + 1][j] + p[i-1]*p[i]*p[j];            s[i][j] = i;            for(int k = i + 1; k < j; ++k){//k 是i,j 之间的断点                int t = m[i][k] + m[k + 1][j] + p[i-1]*p[k]*p[j];//m[i][k] + m[k + 1][j]链长短于当前链，因此已经计算出                if(m[i][j] > t){                    s[i][j] = k;//断开点                    m[i][j] = t;//最小计算量                }            }        }    }}void TrackBack(int i, int j,int s[][LEN + 1]){    if(i == j || i == LEN){        cout<<"A"<<i;        return;    }     cout<<"(";    TrackBack(i,s[i][j],s);    cout<<")(";    TrackBack(s[i][j] + 1,j,s);    cout<<")";}int main(){    MatrixChain(p,m,s);     TrackBack(1,LEN,s);    cout<<endl;}