【算法导论】矩阵链乘法

1.前言

按照算法导论的编排，这一部分应该属于动态规划的范围。今天在这里我就介绍一下如何用动态规划来求解矩阵链乘问题。

矩阵链乘法对于接触过线性代数的人来讲都知道吧，是矩阵运算最基本的东西。一般而言，我们对于几个矩阵的链乘，手算都行，但是对于我这种经常算错数的人来讲，非得借助工具matlab不可。对于A1*A2*...*An这样的问题，我们一般手算通常也就是从左往右来，但是通常而言，效率是很低的。这里说的效率低不止是说手算速度慢，更多的是指我们选择了运算量比较大的方式来计算。

2.问题描述

给定n个矩阵构成的一个链<A1,A2,...,An>，矩阵Ai的维数为mi和ni，其中i=1,2,...n。对乘积<A1,A2,...,An>以一种最小化标量乘法次数的方式进行加全部括号。

注意：在矩阵链乘问题中，实际上并没有把矩阵相乘，目的是确定一个具有最小代价的矩阵相乘顺序。简单讲，即使找出如何打括号的方案，来使得相乘的运算代价最低。

举个例子：矩阵链<A1,A2,A3>，三个矩阵规模分别是(10*100),(100*5),(5*50)。对于这个矩阵链，计算方法有两种，一种是(A1A2)A3，另一种是A1(A2A3)。

计算量分别是：

(1)10*100*5+10*5*50=5000+2500=7500

(2)100*5*50+10*100*50=25000+50000=75000。

显然，第一种方法比第二种效率高了10倍。

3、动态规划分析过程

1）最优加全部括号的结构

动态规划第一步是寻找一个最优的子结构。假设现在要计算AiAi+1....Aj的值，计算AiAi+1....Aj过程当中肯定会存在某个k值（i<=k<j）将AiAi+1....Aj分成两部分，使得AiAi+1....Aj的计算量最小。找到最佳的k之后，分成两个子问题Ai...k和Ak+1...j，只需要继续递归寻找这两个子问题的最优解。

2）递归方案

我们令m[i,j]为AiAi+1....Aj 所需的标量乘法次数的最小值，那么，对于矩阵链<A1,A2,...,An>而言，其实就是求解m[1,n]。假如我们已经得到了最优的k，那么m[i,j]=m[i,k]+m[k+1,j]+m(i)*n(k)*n(j),那么我们的k是如何确定的呢？k的取值只可能有j-i种，所以我们就检查所有的情况，找出最优者就可以了。怎样才算最优？

当i=j时，就一个矩阵，k=i=j，不用算了，m[i,j]=0；

当i<j时，那么m[i,j]=m[i,k]+m[k+1,j]+m(i)*n(k)*n(j)。

3）计算最优代价

如果基于上述的方法进行递归计算，并不比检查所有括号化方案的暴力搜索方法更好。我们采用自底向上的方法来代替，并且对于重叠问题的结果进行记录，对于括号位置k也进行记录，生成两个表格。

4）构造最优解

我们根据第三步记录的表格k来进行最优括号化方案。因为表格里的每一个数实际上就记录的就是表格所在的位置。

5）源代码

#include <iostream>
using namespace std;
//全局变量len，记录了矩阵链的长度
int len;
//m和n表示矩阵为m行n列
struct matrix {
int m;
int n;
};
//带备忘的自底向上表格法
void matrix_chain_order(struct matrix A[],int **num,int **s)
{
int i,j,l,k;
int temp;

for (i=0;i<len;i++)
{
*(*(num+i)+i)=0;
}
for (l=2;l<=len;l++)
{
for(i=0;i<len-l+1;i++)
{
j=i+l-1;
*(*(num+i)+j)=-1;//-1代替正无穷

for(k=i;k<j;k++)
{
temp=*(*(num+i)+k)+*(*(num+k+1)+j)+A[i].m*A[j].n*A[k].n;
if(temp<*(*(num+i)+j)||*(*(num+i)+j)==-1)
{
*(*(num+i)+j)=temp;
*(*(s+i)+j)=k;
}
}
//输出了m[i,j]的最优方案和对应k
cout<<"m["<<i<<","<<j<<"] :      "<<*(*(num+i)+j)<<endl;
cout<<"最优解k为   "<<*(*(s+i)+j)<<endl;

}

}
}
//输出最优括号化方案
void print_optimal_parens(int **s,int i,int j)
{
if(i==j)
cout<<"A"<<i+1;
else
{
cout<<"(";
print_optimal_parens(s,i,*(*(s+i)+j));
print_optimal_parens(s,*(*(s+i)+j)+1,j);
cout<<")";
}
}

int main()
{
cout<<"输入要矩阵个数："<<endl;
int **num=0;
int **s=0;
int i;

cin>>len;

struct matrix *A;
A=new matrix [len];

//num表格表示m[i,j],s表格表示s表格
num=new int*[len];
for (i=0;i<len;i++)
num[i]= new int[len];

s=new int*[len];
for (i=0;i<len;i++)
s[i]= new int [len];

for (i=0;i<len;i++)
{
cout<<"矩阵"<<i+1<<"的行列：";
cin>>A[i].m>>A[i].n;
}

matrix_chain_order(A,num,s);
cout<<"最优括号化方案为 ："<<endl;
print_optimal_parens(s,0,len-1);

delete[] A;
delete[] num;
delete[] s;
while(1);
return 0;
}