《算法导论》学习（十八）----动态规划之矩阵链乘（C语言）

---

前言

本文主要讲解了动态规划中的矩阵链乘问题：给定一个矩阵链，得到它的最小代价计算次序。给出了动态规划方案的分析，并且给出了C语言实现。

---

一、矩阵链乘

1.问题描述

给定一个n个矩阵的序列（矩阵链）$<A_1,A_2,A_3,A_4,...,A_n>$，现在我们希望计算它的乘积
$$A_1{A}_2{A}_3{A}_4...A_n$$
对于矩阵链乘来说，我们可以通过加括号的手段来确定先让哪两个矩阵进行相乘。无论乘的次序如何，最终都不影响结果。
但是不同次序之间的矩阵相乘，计算机所要付出的代价完全不一样。
比如所$<A_1,A_2,A_3>,其中A_1为2\ast 3；A_2为3\ast 4；A_3为4\ast 1$
我们可以计算两种不同策略的切割代价：

• $((A_1{A}_2)A_3)=2\ast 3\ast 4+2\ast 4\ast 1=32$
• $(A_1(A_2{A}_3))=3\ast 4\ast 1+2\ast 3\ast 1=18$

显然，两种代价是完全不一样的。现在我们的问题就是给定一个矩阵链，得到它的最小代价计算次序。

二、问题解决

1.最优化的子问题结构

我们先直接给出相应的递归求解公式：
w [ i , j ] = { 0  if  i = j m i n { w [ i , k ] + w [ k + 1 , j ] + p i − 1 p k p j }  if  i < j w[i,j]=\begin{cases} 0& \text{ if }i=j \\ min\{ w[i,k]+w[k+1,j]+p_{i-1}p_kp_j\}& \text{ if } i<j \\ \end{cases} w[i,j]={0min{w[i,k]+w[k+1,j]+pi−1​pk​pj​}​ if i=j if i<j​

• w[i,j]表示第i个矩阵到第j个矩阵，路径上的矩阵链的最小代价
• $p_{i-1}{p}_k{p}_j$代表从最优子解向它的最优父解构造时需要的代价

2.动态规划

我们采用一个二维数组来记录子问题的最优解。比如w[2,4]的最优解就记录在数组w[1,3]种（因为数组从0开始）。
同时我们从长度为2的子矩阵链开始求最优解，之后长度不断加1，直到最后得到n矩阵链的最优解。
这个过程种采用三层循环

• 第一层循环从小到大遍历矩阵链的长度
• 第二层循环遍历第一层循环确定长度下的子矩阵链中第一个矩阵的位置
• 第三层循环遍历前两层循环确定好的矩阵链的全部”切割“方案

3.最优解构造

经过动态规划的三层循环，我们现在可以得到最小的代价的具体值，但是具体的矩阵相乘顺序是如何，我们需要构造最优解。
我们采用一个二维数s组来辅助我们构造出最优方案。

• 我们用s[0,6]来记录矩阵链w[1,7]第一次加括号的位置，即它的k值。
• 我们记录全部情况的矩阵链的最优k值。
• 然后通过递归调用算法，可以得到括号输出结果。具体方案见下方代码。各位读者也可以通过画递归树来迅速确定方案。

三、C代码

1.代码

#include<stdio.h>
#include<stdlib.h>
#include<time.h>



//由于有n个矩阵链乘
//需要n+1个存储单元来存储行列信息
//此处SIZE为存储单元的个数
//因此矩阵链乘的个数是SIZE-1 
#define SIZE 10
//此处可以规定矩阵的行列数为1~LIM的随机数 
#define LIM 100
//注意：
//由于计算机存储数大小的限制
//SIZE和LIM不宜过大
//过大肯定会使得程序崩溃 



int Matrix_chain_order(int *p,int s[][SIZE-1],int w[][SIZE-1],int size)
{
	int i=0;
	int j=0;
	int k=0;
	int q=0;
	int l=0;
	int r=0;
	for(i=0;i<size-1;i++)
	{
		w[i][i]=0;
		s[i][i]=0;
	}
	for(i=2;i<=(size-1);i++)
	{
		for(j=1;j<=(size-i);j++)
		{
			l=j;
			r=j+i-1;
			w[l-1][r-1]=10000000;
			for(k=l;k<r;k++)
			{
				q=w[l-1][k-1]+w[k][r-1]+p[l-1]*p[k]*p[r];
				if(q<w[l-1][r-1])
				{
					w[l-1][r-1]=q;
					s[l-1][r-1]=k;
				}
			}
		}
	}
	return w[0][size-2];
}



void print_priority(int s[][SIZE-1],int l,int r)
{
	if(l==r)
	{
		printf("A%d",l+1);
		return;
	}
	printf("(");
	print_priority(s,l,s[l][r]-1);
	print_priority(s,s[l][r],r);
	printf(")");		
}



int main()
{
	int p[SIZE];
	int s[SIZE-1][SIZE-1];
	int w[SIZE-1][SIZE-1];
	int min_pay;
	int i=0;
	int j=0;
	srand((unsigned)time(NULL));
	for(i=0;i<SIZE;i++)
	{
		p[i]=rand()%LIM+1;
	}
	for(i=0;i<SIZE;i++)
	{
		printf("%5d",p[i]);
	}
	printf("\n");
	min_pay=Matrix_chain_order(p,s,w,SIZE);
	printf("%5d\n",min_pay);
	print_priority(s,0,SIZE-2);
	return 0;
} 

2.结果

---

总结

文章的不妥之处请各位读者包涵指正。
动态规划的其它问题还可以参考：
《算法导论》学习（十七）----动态规划之钢条切割（C语言）