算法导论--动态规划（矩阵链乘法）

Posted 2020-10-01 ljbguanli

矩阵链乘法问题

给定一个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");
}