动态规划——矩阵链相乘问题（C代码实现）C3-A-Zexal的矩阵链乘

引入

矩阵连乘的开销

给出一个有n个矩阵(任何相邻的两个矩阵均相容)的矩阵序列$（<A_1,A_2,\dots ,A_n>)$求它们的乘积$(A_1{A}_2\dots A_n)$。

一般的做法是按照从左到右的顺序一个一个求矩阵的积，但是矩阵乘法是满足结合率的，我们可以通过给序列加上括号规定矩阵连乘的求解顺序。
矩阵连乘全括号：仅有一个矩阵，或者两个“矩阵连乘全括号”的乘积且外层包括一个括号，例如：$(A_1((A_2{A}_3)A_4))$。

一对规模分别为$p\times q$，$q\times r$的矩阵相乘的乘法开销为$pqr$。对于矩阵连乘问题，采用不同的加括号方式，可能导致不同的甚至差异巨大的乘法开销。

给定三个矩阵$A_1{A}_2{A}_3$，他们的规模分别是$10\times 100$，$100\times 5$，$5\times 50$。
按照$((A_1{A}_2)A_3)$的顺序运算的乘法开销是$10\times 100\times 5+10\times 50\times 5=7500$。
但按照$(A_1(A_2{A}_3))$的顺序运算的乘法开销是$100\times 5\times 50+10\times 100\times 50=75000$。
后者是前者的$10$倍。

因此为了使矩阵乘法运算速度更快，可以预先计算出最好的给矩阵链加括号的方式以减小总乘法开销。

加括号的方式总数

令$P(n)$表示$n$个矩阵连乘时所有可能的全括号的个数。
如果$n=1$，只有一种方法加括号
如果$n>=2$，则在这些矩阵里面可以找到一个分割位置，$A_k$与$A_{k+1}$中间，将所有的分割方式累加。
$$P(n)=\{\begin{array}{ll}1,& \text{if n =1 }\\ {\sum }_{k=1}^{n-1}P(k)P(n-k),& \text{if n ≥ odd}\end{array}$$
这个递归结果将是 $\Omega (2^n)$ （这个计算也可以从另一个角度考虑：$n$个位置可以选择插入或者不插入括号，那么总结果就是$2^n$的）的。如果采用暴力枚举来计算最优括号分割方式的话将消耗大量的时间

用动态规划法解决矩阵链相乘问题

优化子结构

在计算$A_i{A}_{i+1}\dots A_k{A}_{k+1}\dots A_j$的最优括号划分方案的时候，我们先找到一个最佳分割方案把这个矩阵列以$A_k{A}_{k+1}$为边界划分成两块，看成计算$A_i{A}_{i+1}\dots A_k$和$A_{k+1}\dots A_j$。之后再分别在这两块矩阵里寻找最佳边界继续划分，不断形成更小的子问题。
下面证明每一个子问题的最优括号划分方案合并起来就是整个问题的最优括号划分方案。

证明： 对于矩阵链$A_i{A}_{i+1}\dots A_k{A}_{k+1}\dots A_j$，寻找到一个最佳分割将其划分成$(A_i{A}_{i+1}\dots A_k)(A_{k+1}\dots A_j)$，记为$M\cdot N$（M，N，P均代表的是乘法计算次数）。
如果对于$A_i{A}_{i+1}\dots A_k$，能找到一个更佳的括号划分方案的话，得到的结果应该是$P$，且$P<M$。这样原矩阵链的最佳分割方案就应该是$P\cdot N$而不是$M\cdot N$了。
所以每一个子矩阵链的最佳划分方案合并在一起就是整个矩阵链的最佳划分方案。

基于最佳子结构，我们可以从子问题的最优解构造原问题的最优解。

递归求解方案

现在要计算矩阵链$A_i{A}_{i+1}\dots A_k{A}_{k+1}\dots A_j$的最优括号划分方案，最优括号划分即使乘法次数最少的方案，那么我们可以用一个二维数组$m[i,j]$来存储从$i$到$j$的这$j-i+1$个矩阵消耗的最少乘法次数，比如$m[1,n]$存放的就是$A_1$连乘到$A_n$所消耗的最少乘法次数，下面来看看如何计算$m[i,j]$。
在优化子结构里面我们提到了，一个问题的最优括号划分方案等于两个子问题括号划分问题的合并。同样的，消耗的乘法次数也是两个子问题的合并。不过这里要注意的是，分别计算完两个子矩阵链之后还需要把前后两个乘在一起。下面我们规定$p$是读入的每个矩阵的行列数，下标从$0$开始。$n$个矩阵只会读入$n+1$个数据，每相邻两个数代表一个矩阵的行列。

比如要求3个矩阵的连乘，读入4个数据，p₀=10，p₁=20，p₂=30，p₄=40。那么第一个矩阵是10×20的，第二个矩阵是20×30的，第三个矩阵是30×40的。我们可以看到矩阵A_i的行数=p_i-1,列数=p_i，这个结论后面就会用到。

这里假设我们已经找到了最佳分割点$k$位于$A_k$和$A_{k+1}$之间，则$m[i,j]$的计算公式是：
$$m[i,j]=m[i,k]+m[k+1,j]+p_{i-1}{p}_k{p}_j$$
这里$m[i,j]$是要求的最后的总最少的乘法次数，$A_k$和$A_{k+1}$中间是最佳的分割点，$m[i,k]$是分割点前面的矩阵连乘的最少次数，$m[k+1,j]$是分割点后面的矩阵连乘的最少次数，$p_{i-1}{p}_k{p}_j$ 是将前后两个矩阵链的结果合成最终矩阵时消耗的乘法次数。

为什么是$p_{i-1}{p}_k{p}_j$ ？
两个矩阵相乘的时候得到的结果矩阵取得是第一个矩阵的行和第二个矩阵的列。
那么前i到k个矩阵相乘的结果取得就是第$i$个矩阵的行和第$k$个矩阵的列，应用上面的结论的话就是$p_{i-1}$和$p_k$。同理第$k+1$到$j$个矩阵相乘的结果取得就是第$k+1$个矩阵的行$p_k$和第$j$个矩阵的列$p_j$，
一对规模分别为$p\times q$，$q\times r$的矩阵相乘的乘法开销为$pqr$，因此这两个矩阵（$p_{i-1}$×$p_k$，$p_k$×$p_j$）相乘的乘法开销就是$p_{i-1}{p}_k{p}_j$ 。

现在我们已经得到了每一个$m[i,j]$的计算公式，但这是基于最佳分割点$k$已经找到了的情况。实际上$k$是我们需要寻找的东西，而这个分割点一定在$i$和$j$中间，因此$m[i,j]$是对$k$取$i$到$j$中的任意值得到的结果中取最小，它的最终计算公式应该为：
m [ i , j ] = { 0 , if  i  = j   min ⁡ i ≤ k ＜ j   { m [ i , k ] + m [ k + 1 , j ] + p i − 1 p k p j } , if  i < j m[i,j]= \begin{cases} 0, & \text {if $i$ =$j$ } \\\min_{ i≤k＜j } \ \lbrace m[i,k]+m[k+1,j]+p_{i-1}p_kp_j\rbrace , & \text{if $i$<$j$} \end{cases} m[i,j]={0,mini≤k＜j​ {m[i,k]+m[k+1,j]+pi−1​pk​pj​},​if i =j if i<j​
对此我们可以写出用递归求解的c代码：

#define INFINITY 2147483647
int RECURSIVE_MATRIX_CHAIN(int p[],int i,int j)
{
	int k;
	if(i==j)
		return 0;
	m[i][j]=INFINITY;
	for(k=i;k<j;j++)
	{
		q=RECURSIVE_MATRIX_CHAIN(p,i,k)+RECURSIVE_MATRIX_CHAIN(p,k+1,j)+p[i-1]*p[k]*p[j];
		if (q<m[i][j]) m[i][j]=q;
	}
	return m[i][j];
}

时间复杂度分析：
$T(1)\ge 1,$
$T(n)\ge 1+{\sum }_{k=1}^{n-1}(T(k)+T(n-k)+1)=2{\sum }_{i=1}^{n-1}T(i)+n\ge 2^{n-1}$
时间复杂度并没有好于暴力枚举法，这是因为递归的过程中重叠的子问题被反复计算浪费了大量时间。

重叠子问题&构造最优解

基于上面的思路我们发现可以将已经计算出来的$m[i,j]$结果直接填入表中，减少冗余的递归计算。同时我们额外定义一个二维数组$s$用以记录每次分割的位置。具体的实现原因在代码中备注。

比如$5$个矩阵连乘，$s[1,5]$存的就是这一个序列的第一个分割点，假设是$3$，那么接下来的分割点就要在$s[1,3]$和$s[4,5]$里找，一直找到$s$的两个下标相等为止。

据此我们可以写出C代码：

#define INFINITY 2147483647
void MATRIX_CHAIN_ORDER(int p[],int length)//用于计算m和s的函数
{
	int n;
	n=length-1; 
	int i,j,k,q,l;
	for(i=1;i<=n;i++) 
	{
		m[i][i]=0;//只有一个矩阵的时候乘法代价为0
	} 
	for(l=2;l<=n;l++)//求解的顺序是自底向上的，先从规模最小的2个矩阵链求起一直求到矩阵链长度为n为止
	{
		for(i=1;i<=n-l+1;i++)//i用于决定一个矩阵链的起点，这个起点的范围就在1~n-l+1（总长度-当前长度+当前长度中多减的代表最靠后的起点的1）中间
		{
			j=i+l-1;//j决定一个矩阵链的终点，终点就在i+j-1（起点+长度-长度中起点的重复计数）
			m[i][j]=INFINITY;
			for(k=i;k<j;k++)//在i和j中寻找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_OPTIML_PARENS(int i,int j)//用于打印加括号方式的函数
{
	if(i==j)	printf("A%d",i);//当序列里只有一项的时候直接打出
	else
	{
		printf("(");//打出左括号
		PRINT_OPTIML_PARENS(i,s[i][j]);//打出分割点前的元素
		PRINT_OPTIML_PARENS(s[i][j]+1,j);//打出分割点后的元素
		printf(")");//打出右括号
	}
}

一道题目

题面

题目链接
用加括号的方式给出最优的矩阵相乘方案
输入
多组数据输入
第一行一个整数 n，表示矩阵链的长度（1<=n<=300）
接下来一行n+1个数表示这些矩阵的行数和列数
别问我为什么只有n+1个数，每相邻的两个数表示一个矩阵的大小
输出
对于每组数据，输出两行，第一行为计算次数，第二行为计算方案，用加括号的方式给出最优的矩阵相乘方案
如果不幸最优方案不唯一，选择优先计算左边的矩阵
输入样例

3
10 30 5 60
3
10 20 5 4

输出样例

4500
((A1A2)A3)
1200
((A1A2)A3)

Hint
在第二组样例中，选择((A1A2)A3)时，结果为10×20×5+10×5×4=1200
选择A1(A2A3)时，结果为20×5×4 + 10×20×4 = 1200
这时选择第一种，优先计算左边的！

AC代码

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
int N;
int p[305];
int m[305][305];
int s[305][305];

void MATRIX_CHAIN_ORDER(int p[])
{
	int n;
	n=N; 
	int i,j,k,q,l;
	for(i=1;i<=n;i++) 
	{
		m[i][i]=0;
	} 
	for(l=2;l<=n;l++)
	{
		for(i=n-l+1;i>=1;i--)//这里唯一不同的是，因为要求先算左边的，所以改变了i和j的选取顺序
		{
			j=i+l-1;
			m[i][j]=2147483647;
			for(k=j-1;k>=i;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_OPTIML_PARENS(int i,int j)
{
	if(i==j)	printf("A%d",i);
	else
	{
		printf("(");
		PRINT_OPTIML_PARENS(i,s[i][j]);
		PRINT_OPTIML_PARENS(s[i][j]+1,j);
		printf(")");
	}
}
int main()
{
	int i;
	while(scanf("%d",&N)!=EOF)
	{
    	for(i=0;i<N+1;i++) 
		{
			scanf("%d",&p[i]);
		} 
		MATRIX_CHAIN_ORDER(p);
		printf("%d\n",m[1][N]);
		PRINT_OPTIML_PARENS(1,N);
		printf("\n");
		memset(m,0,sizeof(m)); 
		memset(s,0,sizeof(s));
	}
	return 0;
}