动态规划之 矩阵链乘法

承接着上一篇的装配线调度问题，下一个要解决的问题是矩阵链乘法问题：给定n个要相乘的矩阵构成序列<A1,A2,...,An>,计算乘积A1A2...An。为了得到计算的结果，可以将两个矩阵相乘的算法作为一个子程序，根据括号的顺序做全部的矩阵乘法，矩阵的乘法满足结合律，无论怎么加括号都会得到相同的结果，但注意并不满足交换律。例如，4个矩阵的乘积可用5种不同的括号顺序表示：

但是要注意，尽管矩阵括号的顺序并不影响最后的结果，但对运算的代价有很大的影响，先观察矩阵的标准算法：

如果A是p*q矩阵，B是q*r矩阵，那么结果的矩阵C就是p*r矩阵，计算C的时间由第七行的标量乘法次数决定，这个次数等于p*q*r。为了说明加上不同的括号计算的时间是不同的，假设是3个矩阵进行相乘<A1,A2,A3>,维度分别是10*100,100*5,5*50，假设（A1*A2）*A3进行计算的话，那么次数将是10*100*5+10*5*50=7500，但假如是按照A1*(A2*A3)顺序进行计算的话，那么总的次数将是100*5*50+10*100*50=50000,显然使用前一种括号的顺序可以使得速度快了近10倍。

所以，矩阵链乘法问题可以表述如下：给定n个矩阵构成的一个链<A1,A2,...,An>，其中i=1,2,...,n，矩阵Ai的维数为P(i-1)*Pi,对乘积A1A2...An以一种最小化标量乘法次数进行加括号。

计算全部括号的重数

同样的，穷举所有的加括号的方案并不能得到一个较好的算法，假设P(n)表示一串n个矩阵可能的加的括号的数目，当n=1时，只有一个矩阵，只有一种方案，当n>=2时，一个加全部括号的矩阵乘积，就是两个加全部括号的矩阵子乘积的乘积，而且这两个矩阵的分裂可能发生在第k个和第k+1个矩阵之间，得到递归式：

其增长形势仍旧是指数形式，穷举策略并不是一个好的方法。

步骤一：最优加全部括号的结构

和之前提到的一样，动态规划的第一步都是寻找最优的子结构，用子问题的最优解构造出原问题的最优解。对于矩阵链的乘法问题，假设A(i..j)表示乘积AiAi+1...Aj相乘的结果，其中i<=j，如果这个问题是非平凡的，即i<j，那么对乘积AiAi+1...Aj的任何加全部括号形式都将乘积在Ak和Ak+1之间分开，这样计划计算这个乘积就是先计算Ai...Ak,然后计算Ak+1...Aj,最后将这两个结果相乘。

所以，这个问题的最优子结构如下：假设AiAi+1...Aj的一个最优节哀全部括号把乘积在Ak和Ak+1之间分开，则对AiAi+1...Aj最优价全部括号的前缀子链AiAi+1...Ak的加全部括号必须是AiAi+1...Ak的一个最优加全部括号。所以现在就可以通过对子问题的最优解来构造原问题的最优解，一个矩阵链的乘积的最优解法需要分割乘积，且任何最优解都包含了子问题的最优解，所以就可以把问题分割为两个子问题，寻找问题实例的最优解，来构造原问题的最优解。当然在分割某个位置的时候必须保证已经考虑过了所有可能的位置，并检查了最优的一个。

步骤二：一个递归解

第二步就是根据子问题的最优解递归定义一个最优解的代价，计m[i,j]为计算Ai..j所需要的标量乘法的最小值，对于矩阵链，可以像上文中所述的以k位置进行分解，分解为左右两个部分，这样就得到了递归式：

为了跟踪如何构造一个最优解，令s[i,j]为这样的一个k值：在该处分割可以得到一个最优加全部括号。

步骤三：计算最优代价

当然，可以很容易的根据来写最小代价m[1,n]的递归算法，但这个算法也是具有指数时间，它与检查每一种加全部括号乘积的强力法差不多。但是可以注意到，原问题只有相当少的子问题，一个递归算法在递归树的不同分支可能会多次遇到同一个问题，重叠子问题也就是动态规划的第二个特征。

现在，我们通过动态规划的方式进行求解，使用自底向上的表格法来计算最优代价。

假设矩阵Ai的维数是Pi-1*Pi，输入是一个序列P=<P0,P1,...,Pn>，其中length[p]=n+1.程序使用一个辅助表m[1..n,1..n]来保存m[i,j]的代价，用一个辅助表s[1..n,1..n]来记录m[i,j]取得最优代价时k的值。最后利用s构造最优解。

为了正确实现自底向上的方法，必须要确定哪些表项要被用来计算m[i,j]，当计算j-i+1个矩阵相乘时，其代价m[i,j]仅依赖于计算一个有少于j-i+1个矩阵乘积的代价。下面是伪代码：

该算法首先对2~3行中的i=1,2...n置m[i,i]=0，这是长度为1的矩阵链的代价，在4~12行循环的第一次执行中，利用递归式来计算m[i,i+1]，循环的第二次计算m[i,i+2]，然后一直进行下去。

下图显示了算法执行中m和s的形态，注意图中只记录了一半。

步骤四：构造一个最优解

通过之前给出的函数中的s计算，在每一个表项s[i,j]中，记录了对乘积Ai*Ai+1...*Aj进行分裂的最优k值。因此，按照最优方式计算Ai..n时，最终矩阵的相乘顺序是A1..s[1,n]*As[1,n]+1..n，之前的乘法可以递归给出。

最后是程序的代码：

#include <stdio.h>#include <stdlib.h>int p[7]={30,35,15,5,10,20,25};int m[6+1][6+1];int s[6+1][6+1];void MATRIX_CHAIN_ORDER(){    int n=7-1;    int i,l,j,k,q;    for(i=1;i<=n;i++)        m[i][i]=0;    for(l=2;l<=n;l++)    {        for(i=1;i<=n-l+1;i++)        {            j=i+l-1;            m[i][j]=1000000;            for(k=i;k<=j-1;k++)            {                q=m[i][k]+m[k+1][j]+p[i-1]*p[k]*p[j];                if(q<m[i][j])                {                    m[i][j]=q;                    s[i][j]=k;                }            }        }    }}void PRINT_OPTIMAL_PARENS(int i,int j){    if(i==j)        printf("A%d",i);    else    {        printf("(");        PRINT_OPTIMAL_PARENS(i,s[i][j]);        PRINT_OPTIMAL_PARENS(s[i][j]+1,j);        printf(")");    }}int main(){    int i,j;    printf("输入的矩阵维数如下：\n");    for(i=0;i<6;i++)        printf("A%d=%d*%d ",i+1,p[i],p[i+1]);    printf("\n");    printf("得到的矩阵为：\n");    MATRIX_CHAIN_ORDER();    for(i=1;i<=6;i++)    {        for(j=1;j<=6;j++)            printf("%6d ",m[i][j]);        printf("\n");    }    printf("\n加括号的顺序是\n");    PRINT_OPTIMAL_PARENS(1,6);    return 0;}