经典动态规划引例--矩阵链相乘

转自http://blog.youkuaiyun.com/hearthougan/article/details/25834141

这个感觉有必要说一下，因为很多经典的问题都是以它为根基扩展的，譬如：石子合并类型的。    

  给定由n个要相乘的矩阵构成的序列：<A1, A2, A3···An>，要计算乘积：A1*A2*A3*····An  ---- <1>为了计算<1>式乘积，我们知道矩阵相乘是满足结合律的，故无论怎么添加括号，都会产生相同的结果。例如：矩阵链<A1, A2, A3,  A4>乘积A1*A2*A3*A4可用五种不同的方式添加括号：

(A1*(A2*(A3*A4)))

(A1*((A2*A3)*A4))

((A1*A2)*(A3*A4))

((A1*(A2*A3))*A4)

(((A1*A2)*A3)*A4)

    矩阵链添加括号的方式，对矩阵的运算有很大的影响，首先看一下两个矩阵相乘的代价：

Matrix Multi_Matrix(Matrix a, Matrix b, int n)
{
1    Matrix c;
2   int i, j, k;
3    for(i = 0; i < row(A); ++i)
4   {
5      for(j = 0; j < col(B); ++j)
6        {
7           c.iMatrix[i][j] = 0;
8           for(k = 0; k < col(A); ++k)
9            {
10                c.iMatrix[i][j] ^= (a.iMatrix[i][k] & b.iMatrix[k][j]); 
11           }
12         }
13    }
14    return c;
}

      我们都直到如果两个矩阵是相容的（即A的列数等于矩阵B的行数）时，才可以进行运算，设矩阵A,B分别为p*q矩阵,q*r矩阵，则结果矩阵C为p*r的，计算C的时间由第10行所决定的，运算次数为：p*q*r

到此，为说明矩阵链不同的加括号方式所带来的矩阵乘积代价是不同的，现以下为例：

     设有四个矩阵A1,A2,A3,A4，它们的维数分别是:

   

             A1:p1×p2=50×10，A2:p2×p3=10×40

           A3:p3×p4=40×30，A4:p4×p5=30×5

    因此最优的添加括号顺序，会为运算带来更多的回报。那么如何添加括号呢？即最优的断开位置(在断开位置添加括号，问题划归为两个子问题)在哪儿呢？

如果枚举每一个加括号的位置：

    设P(n)表示n个矩阵链相乘的方案数，则当n = 1时, 是一个平凡问题，只有一个矩阵，则P(n) = 1；当n>=2时，是一个非平凡问题，当在第k个位置添加括号时，则有P(k)*P(n-k)种方案数，由于k的位置是在[1, n-1]之间因此有如下递推公式：

这个解是一个Catalan数序列，解的个数是关于n的指数形式,因此暴力不是一个好的方法。

步骤一：

    寻找最优子结构，利用子问题的最优解来求解原问题的最优解，对于矩阵链相乘，我们记：A[i··j]表示A[i]*A[i+1]*····A[j]的乘积结果，其中Ai的维数表示为p[i]*p[i+1]；

    由于矩阵相乘满足结合律，那么计算A[i]*····A[j]，可以在A[k]和A[k+1]之间分开,即首先计算A[i···Ak] = A[i]*A[i+1]*····A[k]和A[k+1···j] = A[k+1]*A[k+2]*···A[j]然后计算两者相乘，这样就可以求出A[i]*A[i+1]*···A[j], 其中k的可能位置：i <= k < j;如此，A[i]*A[i+1]*···A[j]的计算代价，就是计算A[i]*A[i+1]*···A[k]与A[k+1]*A[k+2]*····A[j]代价之和，然后再加上两者相乘的代价。

    那么如果A[i]*A[i+1]*···A[j]的一个最优加括号方式是在A[k]和A[k+1]处分开，那么对于它的“前缀”A[i]*A[i+1]*···A[k]的最优加括号方式也一定是最优的，为什么呢？如果A[i]*A[i+1]*···A[k]存在一个更优的加括号顺序，那么把它带入到A[i]*A[i+1]*···A[j]中，则A[i]*A[i+1]*···A[j]的计算代价一定比最优加括号顺序的代价小，矛盾了！同理[k+1]*A[k+2]*······A[j]的加括号顺序也一定是最优的。这样矩阵链相乘就满足了最优子结构性质。

步骤二：

    找出递推公式:首先交代一下三个数组的含义：m[i][j]记录的是Ai乘到Aj的最小代价，s[i][j] = k表示A[i]*A[i+1]*···A[j]的最佳断开位置为k，p[i]表示矩阵Ai的行数，p[i+1]表示它的列数。

    假设A[i]*A[i+1]*···A[j]的最优断开位置为k,那么A[i]*A[i+1]*····A[j]记为：m[i][j]；则：A[i]*A[i+1]*····A[k]的最小计算代价为：m[i][k]；[k+1]*A[k+2]*····A[j]的计算代价记为：m[k+1][j].而A[i···k]*A[k+1····j]可表示为：p[i]*p[k+1]*p[j+1]。则m[i][j] = m[i][k] + m[k+1][j] + p[i]*p[k+1]*p[j+1]。

    由于实际上k的位置是在[i, j-1]之间，因此k的位置是不确定的,那么可以得出递推公式：

步骤三：

    如果利用递归求解，那么可能会处理很多相同的子问题，那么处理重叠子问题，也是动态规划一个特性，例如：

考虑四个矩阵的连乘积A1A2A3A4，其中

A1:p1×p2,  A2:p2×p3,  A3:p3×p4,  A4:p4×p5

    上图是利用自顶向下的方式求解，会有很多重叠的子问题，多次处理，从而降低了效率，如：m[2][3],m[3][4]处理了两次，等等，如果规模更大的话，重叠的子问题也会随之而增加。

    因此我们利用自底向上的方式来求解问题，首先处理长度为1的矩阵链，即只有一个矩阵，有n个，那么相乘的次数为0，即：m[i][i] = 0，1<=i<=n；

    长度为2的矩阵链，即只有两个矩阵相乘例如：A1*A2， A2*A3，····A[i][j]····，A[n-1]*A[n]，共有n-1 = n+2-1个，其中的j如何求呢？显然它等于j = i+r-1。

·····································

    长度为r的矩阵链，有n-r+1个。

····················

    显然长度为n的只有一个即：m[1][n]，也为最终的求解。

到此，理论完事。

代码：

[cpp]  view plain copy

1. #include <iostream>  
2. #include <cstdio>  
3. #include <cstring>  
4. #include <climits>  
5.   
6. using namespace std;  
7.   
8. const int MAXN = 110;  
9.   
10. int m[MAXN][MAXN], s[MAXN][MAXN], p[MAXN];  
11.   
12. void Multi_Matrix(int n)  
13. {  
14.     int i, j, r, k, t;  
15.     for(i = 0; i < MAXN; ++i)  
16.     {  
17.         for(j = 0; j < MAXN; ++j)  
18.         {  
19.             if(i == j)  
20.                 m[i][j] = 0;  
21.             else  
22.                 m[i][j] = INT_MAX;  
23.             s[i][j] = 0;  
24.         }  
25.     }  
26.     for(r = 2; r <= n; ++r)//查找长度为r的矩阵连乘  
27.     {  
28.         for(i = 1; i <= n-r+1; ++i)  
29.         {  
30.             j = i+r-1;  
31.             for(k = i; k < j; ++k)  
32.             {  
33.                 t = m[i][k] + m[k+1][j] + p[i]*p[k+1]*p[j+1];  
34.                 if(t < m[i][j])  
35.                 {  
36.                     m[i][j] = t;  
37.                     s[i][j] = k;  
38.                 }  
39.             }  
40.         }  
41.     }  
42. }  
43.   
44.   
45. void Add(int i, int j)  
46. {  
47.     if(i == j)  
48.     {  
49.         printf("A%d", i);  
50.         return ;  
51.     }  
52.     printf("(");  
53.     Add(i, s[i][j]);  
54.     printf("*");  
55.     Add(s[i][j]+1, j);  
56.     printf(")");  
57. }  
58.   
59. int main()  
60. {  
61.     int n, i;  
62.     while(scanf("%d", &n))  
63.     {  
64.         memset(p, 0, sizeof(p));  
65.         for(i = 1; i <= n+1; ++i)  
66.             scanf("%d", &p[i]);  
67.         Multi_Matrix(n);  
68.         printf("%d\n", m[1][n]);  
69.         Add(1, n);  
70.         printf("\n");  
71.     }  
72.     return 0;  
73. }  
74.   
75. /** 
76. 6 
77. 6 7 3 1 2 4 5 
78.  
79. 4 
80. 50 10 40 30 5 
81. */