动态规划--矩阵连乘的最优乘法顺序

载自：蓝飞技术部落格

最优矩阵连乘积

Accepted: 10 Total Submit: 18

Time Limit: 1000ms Memony Limit: 32768KB

Description

在科学计算中经常要计算矩阵的乘积。矩阵A和B可乘的条件是矩阵A的列数等于矩阵B的行数。若A是一个p×q的矩阵，B是一个q×r的矩阵，则其乘积C=AB是一个p×r的矩阵。其标准计算公式为：





由该公式知计算C=AB总共需要pqr次的数乘。

为了说明在计算矩阵连乘积时加括号方式对整个计算量的影响，我们来看一个计算3个矩阵{A1，A2，A3}的连乘积的例子。设这3个矩阵的维数分别为10×100，100×5和5×50。若按第一种加括号方式((A1A2)A3)来计算，总共需要10×100×5+10×5×50=7500次的数乘。若按第二种加括号方式(A1(A2A3))来计算，则需要的数乘次数为100×5×50+10×100×50=75000。第二种加括号方式的计算量是第一种加括号方式的计算量的10倍。由此可见，在计算矩阵连乘积时，加括号方式，即计算次序对计算量有很大影响。

于是，人们自然会提出矩阵连乘积的最优计算次序问题，即对于给定的相继n个矩阵{A1,A2,…,An}(其中Ai的维数为pi-1×pi ，i=1,2,…,n)，如何确定计算矩阵连乘积A1A2…An的一个计算次序(完全加括号方式)，使得依此次序计算矩阵连乘积需要的数乘次数最少。

Input

有若干种案例，每种两行，第一行是一个非负整数n表示矩阵的个数，n=0表示结束。接着有n行，每行两个正整数，表示矩阵的维数。

Ouput

对应输出最小的乘法次数。

Sample Input

3

10 100

100 5

5 50

6

30 35

35 15

15 5

5 10

10 20

20 25

0

Sample Output

7500

15125

Source

思路是采用动态规划：

设a
[m]为第n个矩阵到第m个矩阵连乘的最小乘法数（n, m >= 1），b[i], b[i + 1]为第i个矩阵的行数和列数（i >= 0），那么：

a
[n + 1]易求，为相邻两个矩阵相乘的乘法数，即b[n - 1] * b
* b[n + 1]；
An - Am可以任意拆分为An - Ai及Ai+1 - Am两部分相乘（n <= i <= m），则a
[m]为所有拆分情况中乘法次数最少的一种，即a
[m] = min(a
[i] + a
+ b[n - 1] * b[i] * b[m])，注意b的下标是从0开始，比a的下标少1。

这样说可能比较抽象，我们举个例子：

A1 * A2 * A3 * A4可以拆分为以下几种形式，且最小乘法数必是以下情况中的最乘法数量：

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

而A1 * A2 * A3与A2 * A3 * A4又可以拆分为以下两种形式：
A1 * (A2 * A3)
(A1 * A2) * A3

以此类推，一直到相邻两个矩阵的相乘，而两个矩阵的最小乘法数即是两个矩阵相乘所需的乘法数，易求。

而计算顺序则与分析相反，由两两相邻的矩阵开始，一直推算到所有矩阵的最优结果。以第二个测试数据为例，a
[m]计算结果如下：

a	1	2	3	4	5	6
1	0	15750	7875	9375	11875	15125
2	0	0	2625	4375	7125	10500
3	0	0	0	750	2500	5375
4	0	0	0	0	1000	3500
5	0	0	0	0	0	5000
6	0	0	0	0	0	0

易知计算顺序为由左下往右上，最终得出的a[1][6]（表示从第1个矩阵到第6个矩阵的最小乘法数量）便是所求答案。

一种思路是三层循环嵌套，第一层为所求矩阵个数，第二层为所求矩阵开始位置，第三层为拆分的所有子情况，代码如下：

#include<iostream>

#include<vector>

using namespace std;

inline int min(int a, int b){

        return a > b ? b : a;

}

int main(){

        int t, i, j, r, a[101][101], b[101];

        while(cin >> t && t){

                for(i = 0; i < t; ++i){

                        cin >> b[i] >> b[i + 1];

                }

                for(i = 0; i <= 100; ++i){

                        for(j = 0; j <= 100; ++j){

                                a[i][j] = 0;    //全部归零

                        }

                }

                for(r = 2; r <= t; ++r){        //求相邻r个矩阵的最优值

                        for(i = 1; i + r - 1 <= t; ++i){        //从第i个矩阵开始的r个矩阵

                                int num = 0xffffff;     //定义一个较大的值

                                for(j = i; j < i + r - 1; ++j){ //从第i个矩阵开始的r个矩阵拆分为i到j和j+1到i+r-1两个部分

                                        num = min(num, a[i][j] + a[j + 1][i + r - 1] + b[i - 1] * b[j] * b[i + r - 1]);

                                }

                                a[i][i + r - 1] = num;

                                //cout << i << ":" << i + r - 1 << " " << a[i][i + r - 1] << " ";

                        }

                        //cout <<  endl;

                }

                cout << a[1][t] << endl;

        }

        return 0;

}

三层循环嵌套可能比较难理解，当时在实验室也是因为三层循环三个变量以及两个数组下标之间的关系没有理清楚所以搞的一片混乱，像这种关联数据非常多的情况应该先把关系分析好再来写代码，否则改来改去很难找清楚到底是哪里出的问题，后来直接把三个循环全删了重新写才把这道题做出来。

其实对于这种有点类似斐波那契数列的递推问题我们也可以通过递归来解决，从而免去两层循环：

#include<iostream>

#include<vector>

using namespace std;

int b[101], sum = 0;

inline int min(int a, int b){

        return a > b ? b : a;

}

int GetMuls(int n, int m){

        if(n == m)

                return 0;

        ++sum;

        int num = 0xffffff;     //足够大的数

        for(int i = n; i < m; ++i){     //将n - m的矩阵拆分为n到i及i + 1到m两个部分

                num = min(num, GetMuls(n, i) + GetMuls(i + 1, m) + b[n - 1] * b[i] * b[m]);     //将最小值保存下来

        }

        return num;

}

int main(){

        int t, i;

        while(cin >> t && t){

                for(i = 0; i < t; ++i){

                        cin >> b[i] >> b[i + 1];

                }

                cout << GetMuls(1, t) << endl;

                //cout << sum << endl;  //统计递归次数

        }

        return 0;

}

使用递归代码顿时简洁了许多，甚至连a
[m]都省了。虽然可以OJ上通过，但这种做法并不推荐，因为这是一种效率非常低的做法，类似斐波那契数列，该问题可以分解为两个子问题，而每个子问题又都可以分解为更小的两个子问题，递归的次数呈几何倍数增长，若测试数据较大的话这种算法必然超时。同样，我们也可以利用斐波那契数列问题的解决思想，将计算过的值先储存起来，再次用到的时候直接返回，用空间换时间。剪枝后的代码如下：

#include<iostream>

#include<vector>

using namespace std;

int a[101][101], b[101], sum = 0;

inline int min(int a, int b){

        return a > b ? b : a;

}

int GetMuls(int n, int m){

        if(n == m)

                return 0;

        if(a[n][m] > 0) //如果a
[m]已经求出

                return a[n][m]; //直接返回a
[m]的值

        ++sum;

        int num = 0xffffff;

        for(int i = n; i < m; ++i){

                num = min(num, GetMuls(n, i) + GetMuls(i + 1, m) + b[n - 1] * b[i] * b[m]);

        }

        a[n][m] = num;  //保存进数组中

        return num;

}

int main(){

        int t, i, j;

        while(cin >> t && t){

                for(i = 0; i < t; ++i){

                        cin >> b[i] >> b[i + 1];

                }

                for(i = 0; i <= 100; ++i){

                        for(j = 0; j <= 100; ++j){

                                a[i][j] = 0;

                        }

                }

                cout << GetMuls(1, t) << endl;

                //cout << sum << endl;

        }

        return 0;

}

大家可以把上面统计递归次数的sum输出前的注释去掉，对比一下就可以明显地看出时间复杂度的区别。

上面所求的都是最优值，最后再说一说最优解的求法：

我们把每一步求得最优值时的拆分方法记录下来，保存于数组c
[m]中，c
[m]表示将An - Am的最佳拆分方案为An - Ac
[m]及Ac
[m] + 1 - Am两个部分，仍以第二个测试数据为例：

c	1	2	3	4	5	6
1	0	1	1	3	3	3
2	0	0	2	3	3	3
3	0	0	0	3	3	3
4	0	0	0	0	4	5
5	0	0	0	0	0	5
6	0	0	0	0	0	0

得到最佳拆分方案后的问题是，如何利用c
[m]为矩阵连乘式子加上括号，用字符串保存矩阵式子再在字符串中插入括号是不切实际的，因为插入括号后代表矩阵的字母位置会改变，很难确定下一个括号插入的位置。
我所用的方法是使用二维数组来储存每一个矩阵间隙间（间隙数量比矩阵数量多一个）的两种括号数量，d[i][0]为第i个空隙中左括号的数量，d[i][1]为第i个空隙中右括号的数量：

#include<iostream>

#include<vector>

using namespace std;

int a[101][101], b[101], sum = 0;

int c[101][101], d[101][2];

inline int min(int a, int b){

        return a > b ? b : a;

}

int GetMuls(int n, int m){

        if(n == m)

                return 0;

        if(a[n][m] > 0)

                return a[n][m];

        ++sum;

        int num = 0xffffff;

        int tmp = 0xffffff;

        for(int i = n; i < m; ++i){

                num = min(num, GetMuls(n, i) + GetMuls(i + 1, m) + b[n - 1] * b[i] * b[m]);

                if(num < tmp){

                        tmp = num;

                        c[n][m] = i;    //记录最佳分解方案时i的值

                        //cout << i << endl;

                }

        }

        a[n][m] = num;

        return num;

}

void Solve(int n, int m){       //求出n到m最佳分解方案的括号数量及位置

        if(m - n <= 1)

                return;

        int i = c[n][m];

        if(i - n > 0){  //在不止一个矩阵的情况下才加括号，防止单个矩阵直接被括号包围

                ++d[n][0];      //第n个空隙处左括号数量加1

                ++d[i + 1][1];  //第i + 1个空隙处右括号数量加1

        }

        if(m - i - 1 > 0){      //在不止一个矩阵的情况下才加括号，防止单个矩阵直接被括号包围

                ++d[i + 1][0];  //第i + 1个空隙处左括号数量加1

                ++d[m + 1][1];  //第m + 1个空隙处右括号数量加1

        }

        Solve(n, i);    //计算n到i的括号数量及位置

        Solve(i + 1, m);        //计算i + 1到m的括号数量及位置

}

int main(){

        int t, i, j;

        while(cin >> t && t){

                for(i = 0; i < t; ++i){

                        cin >> b[i] >> b[i + 1];

                }

                for(i = 0; i <= 100; ++i){

                        for(j = 0; j <= 100; ++j){

                                a[i][j] = 0;

                                c[i][j] = 0;

                        }

                        d[i][0] = 0;

                        d[i][1] = 0;

                }

                cout << GetMuls(1, t) << endl;

                /*

                for(i = 1; i <= t; ++i){

                        for(j = 1; j <= t; ++j){

                                cout << c[i][j] << " ";

                        }

                        cout << endl;

                }

                */

                Solve(1, t);    //计算1到t的括号数量及位置

                for(i = 1; i <= t + 1; ++i){    //输出最优解

                        for(j = 0; j < d[i][1]; ++j){

                                cout << ")";    //由于不可能出现左括号右括号相邻的情况，因此右括号先输出

                        }

                        if(i > 1 && i < t + 1)

                                cout << "*";    //矩阵之间输出乘号

                        for(j = 0; j < d[i][0]; ++j){

                                cout << "(";    //输出左括号

                        }

                        if(i < t + 1)

                                cout << "A" << i;       //输出代表矩阵的字母

                }

                cout << endl;

                //cout << sum << endl;

        }

        return 0;

}

测试案例输出结果：

7500

(A1*A2)*A3

15125

(A1*(A2*A3))*((A4*A5)*A6)