矩阵链乘法问题

给定一个n个矩阵的序列$\langle A_1,A_2,A_3...A_n\rangle$,我们要计算他们的乘积：$A_1A_2A_3...A_n$，由于矩阵乘法满足结合律，加括号不会影响结果，但是不同的加括号方法，算法复杂度有很大的差别：
考虑矩阵链$：\langle A_1,A_2,A_3\rangle$，三个矩阵规模分别为$10\times100、100\times5、5\times50$
如果按$((A_1A_2)A_3)$方式，需要做$10*100*5=5000$次，再与$A_3$相乘，又需要$10*5*50=2500$，共需要7500次运算：
如果按$(A_1（A_2A_3）)$方式计算，共需要$100*5*50+10*100*50 =75000$次标量乘法，具有10倍的差别。可见一个好的加括号方式，对计算效率有很大影响。
为了得到所需乘法次数最少的方案，需要计算所有种方案的代价。
对一个n个矩阵的链，令P(n) 表示可供选择的括号化方案的数量。
完全括号化方案的矩阵乘积可以描述为两个完全括号化的部分相乘的形式，
$P(n) = 1$ , $n= 1$
$P(n) =\sum_{k=1}^{n-1}P(k)P(n-k)$, $n\geq 2$
k为分割点，即第k个矩阵和第k+1个矩阵之间
可以看出括号化方案的数量与n呈指数关系$Ω(2^n)$，若采用暴力搜索比较所有括号化方案的代价，效率很差！

动态规划

1.最优括号化方案的结构特征

由于要求得矩阵链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$的最优括号化方案,我们可以将问题划分为两个子问题$\langle A_i,A_{i+1}...A_k\rangle$和$\langle A_{k+1},A_{k+2}...A_j\rangle$的最优括号化方案的组合，这也是可以采用动态规划的一个重要标示。即一个大的问题的解是其子问题的组合。我们需要遍历所有的k值$i \leq k \leq j-1$即考查所用的划分点。

2.递归求解方案

令$m[i,j] ,i \leq j$标示矩阵链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$的最优括号化方案所需乘法次数的最小值。
当$i = j$时，$m[i,j] = 0$，只有一个矩阵不涉及乘法运算
当$i < j$时，假设在矩阵链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$分割点位置为$A_k和A_{k+1}$之间，如上分析，$m[i,j]为 m[i,k]与m[k+1,j]$的代价和，还要加上两者最后相乘涉及的运算次数。假如$A_{i}$的大小为$p_{i-1} \times p_i$,则子矩阵链$m[i,k]$乘积后的矩阵大小为$p_{i-1} \times p_k$, $m[k+1,j]$乘积后的大小为$p_{k} \times p_j$,所以最后一次乘积做的乘法运算次数为$p_{i-1}p_kp_j$。
即：
$m[i,j] = 0, (i=j)$
$m[i,j] = min\{m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\},i<j，i \leq k \leq j-1$

3.计算最优代价

递归算法会多次遇到同一个子问题，与钢铁切割很类似，每一次高层的运算，都会调用底层结果，越是底层，被调用的次数越多。所以可以采用自底向上的方法，先对底层逐个求解，当上层需要调用底层时，底层已经被求解完毕。
用m[i][j]二维矩阵保存对应链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$长度为j-i+1的最优计算代价q。
用s[i][j]二维矩阵保存对应链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$长度为j-i+1的最优划分位置k。

void Matrix_Chain_Order(int p[],int n)
{
   int i,j,L,k,q;
   for (i=1;i<=n;i++)      //先对单个矩阵的链，求解，即所有m[i][i] =0;
   {
     m[i][i]=0;          
   }
   for(L=2;L<=n;L++)     //从两个矩阵链的长度开始，逐次增加矩阵链的长度
       for(i=1;i<=n-L+1;i++)  //在给定p[]中的矩阵链中，对所有种长度为L的情况计算
       {
           j = i+L-1;
           m[i][j] = -1;
           for(k=i;k<=j-1;k++)   //遍历所有可能的划分点k，计算出最优的划分方案
           {
             q = m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];//计算划分的代价
             if ( q < m[i][j] || m[i][j] == -1)  
             {
                m[i][j] = q;     //最优的代价q保存在m[i][j]中
                s[i][j] = k;     //最优的划分位置k保存在s[i][j]中
             }
           }
       }

}

构造最优解

矩阵链$\langle A_i,A_{i+1},A_{i+2}...A_j\rangle$，由于二维矩阵s[i][j]记录了对应划分位置k，指出了应该在$A_k和A_{k+1}$之间，同样在矩阵链\left[中最优划分位置一定保存在数组$s \left[i\right]\left[s[i,j]\right]$内，矩阵；链$\langle A_{k+1},A_{k+2}...A_j\rangle$的最优划分位置一定保存在$s \left[s[i][j]+1]\right]\left[j\right]$数组内，可以不断递归出最优解。

例程

/************************************************************************
CSDN 勿在浮沙筑高台 
http://blog.csdn.net/luoshixian099
算法导论--动态规划（矩阵链乘法）
2015年6月3日                    
************************************************************************/
#include <STDIO.H>
#include <STDLIB.H>
int m[7][7]={0};
int s[7][7]={0};
void Print_Optimal_Parens(int s[][7],int i,int j)  //构造最优解
{
   if ( i ==j)
   {
       printf("A%d",i);
   }
   else
   {
       printf("(");
       Print_Optimal_Parens(s,i,s[i][j]);
       Print_Optimal_Parens(s,s[i][j]+1,j);
       printf(")");
   }
}
void Matrix_Chain_Order(int p[],int n)
{
   int i,j,L,k,q;
   for (i=1;i<=n;i++)      //先对单个矩阵的链，求解，即所有m[i][i] =0;
   {
     m[i][i]=0;          
   }
   for(L=2;L<=n;L++)     //从两个矩阵链开始，逐次增加矩阵链的长度
       for(i=1;i<=n-L+1;i++)  //在给定p[]中的矩阵链中，对所有种长度为L的情况计算
       {
           j = i+L-1;
           m[i][j] = -1;
           for(k=i;k<=j-1;k++)   //遍历所有可能的划分点k，计算出最优的划分方案，
           {
             q = m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];
             if ( q < m[i][j] || m[i][j] == -1)  
             {
                m[i][j] = q;     //最优的代价q保存在m[i][j]中
                s[i][j] = k;     //最优的划分位置k保存在s[i][j]中
             }
           }
       }

}
void main()
{
    int p[]={30,35,15,5,10,20,25};    //矩阵的输入
    int length = sizeof(p)/sizeof(p[0])-1;   //矩阵长度
    int i,j;
    Matrix_Chain_Order(p,length);

    for(i =1;i<=6;i++)
    {
        for (j=1;j<=6;j++)
        {
            printf("%8d",m[i][j]);
        }
        printf("\n");
    }

    Print_Optimal_Parens(s,1,6);
    printf("\n");
}