【算法导论】矩阵链乘法

1.前言       

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

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

2.问题描述

    给定n个矩阵构成的一个链<A1,A2,...,An>，矩阵A_i的维数为m_i和n_i，其中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）最优加全部括号的结构

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

2）递归方案

    我们令m[i,j]为A_iA_i+1....A_{j 所需的标量乘法次数的最小值，那么，对于矩阵链<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;
}