《算法设计与分析》学习笔记_第六章_矩阵链乘法问题

1、矩阵链乘法问题的描述

　　给定n个矩阵构成的一个链<A₁,A₂,A₃,.......A_n>，其中i=1,2,...n，矩阵A的维数为p_i-1p_i，对乘积 A₁A₂...A_n 以一种最小化标量乘法次数的方式进行加全部括号。

　　注意：在矩阵链乘问题中，实际上并没有把矩阵相乘，目的是确定一个具有最小代价的矩阵相乘顺序。找出这样一个结合顺序使得相乘的代价最低。

2、动态规划分析过程

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

    假设AiAi+1....Aj的一个最优加全括号把乘积在Ak和Ak+1之间分开，则Ai..k和Ak+1..j也都是最优加全括号的。

设m[i,j]为计算机矩阵A_i...j所需的标量乘法运算次数的最小值，对此计算A_1..n的最小代价就是m[1,n]。现在需要来递归定义m[i,j]，分两种情况进行讨论如下：

当i==j时：m[i,j] = 0，（此时只包含一个矩阵）

当i<j时：从步骤1中需要寻找一个k（i≤k＜j）值，使得m[i,j]＝min{m[i,k]+m[k+1,j]+p_i-1p_kp_j}（i≤k＜j）。

设矩阵Ai的维数为pi-1pi，i=1,2.....n。输入序列为：p=<p0，p1,...pn>，length[p] = n+1。使用m[n][n]保存m[i,j]的代价，s[n][n]保存计算m[i,j]时取得最优代价处k的值，最后可以用s中的记录构造一个最优解。

动态规划算法的伪代码如下：

 DPMaxtix_Chain (p)
 2   n = length[p]-1;
 3   for i=1 to n do
 4        m[i][i] = 0;
 5   for t = 2 to n do //t 是矩阵链的长度
 6        for i=1 to n-t+1 do
 7        		j=i+t-1;
 8              m[i][j] = MAXLIMIT;
 9       		for k=i to j-1 do
10             		q = m[i][k] + m[k+1][i] + qi-1qkqj;
11                  if q < m[i][j] then
12                  	m[i][j] = q;
13                      s[i][j] = k;
14   return m and s;

DPMaxtix_Chain具有循环嵌套，深度为3层，运行时间为O(n^3)

3、递归

int Recursive_Matrix_Chain(int *p, int i, int j, int m[N+1][N+1], int s[N+1][N+1])
{
	if(i == j)
		m[i][j] = 0;
	else
	{
		int k;
		m[i][j] = MaxValue;
		for(k=i;k<j;k++)
		{
			int temp = Recursive_Matrix_Chain(p,i,k,m,s)+Recursive_Matrix_Chain(p,k+1,j,m,s)+ p[i-1]*p[k]*p[j];
			if(temp < m[i][j])
			{
				m[i][j] = temp;
				s[i][j] = k;
			}
		}
	}
	return m[i][j];
}

如果采用递归进行实现，则需要指数级时间Ω(2^n)，因为中间有些重复计算。

4、备忘录 对递归算计的改进，可以引入备忘录，采用自顶向下的策略，维护一个记录了子问题的表，控制结构像递归算法。

int Memoized_Matrix_Chain(int *p, int m[N+1][N+1], int s[N+1][N+1])
{
	int i,j;

	for(i = 1; i <= N; i++)
		for(j = 1; j <= N; j++)
		{
			m[i][j] = MaxValue;
		}

		return Lookup_Chain(p, 1, N, m, s);
}

int Lookup_Chain(int *p, int i, int j, int m[N+1][N+1], int s[N+1][N+1])
{
	if(m[i][j] < MaxValue)
	{
		return m[i][j]; //直接返回，相当于查表
	}

	if(i == j)
	{	
		m[i][j] = 0;
	}

	else
	{
		int k;

		for(k = i; k < j; k++)
		{
			int temp = Lookup_Chain(p, i, k, m, s) + Lookup_Chain(p, k+1, j,m, s) + p[i-1] * p[k] * p[j];  //通过递归的形式计算，只计算一次，第二次查表得到
			
			if(temp < m[i][j])
			{
				m[i][j] = temp;
				s[i][j] = k;
			}
		}
	}

	return m[i][j];
}

我们可以根据s[i][k]和s[k+1][j]的值递归加括号，算法如下

Print_Optimal_Parens(s, i, j)
	if i== j then 
		print "Ai"
	else
		print "(";
		Print_Optimal_Parens(s, i, s[i][j]);
		Print_Optimal_Parens(s, s[i][j]+1, j);
		print")";

采用C++语言实现这个过程，现有矩阵A ₁ (30×35)、A ₂ (35×15) _、 A3(15×5)、A4(5×10)、A5(10×20)、A6(20×25)，得到p=<30,35,15,5,10,20,25>。实现过程定义两个二维数组m和s，为了方便计算其第一行和第一列都忽略，行标和列标都是1开始。完整的程序如下所示：

#include <iostream>
using namespace std;

#define N 6
#define MaxValue 1000000

void Matrix_Chain(int *p, int len, int m[N + 1][N + 1], int s[N + 1][N + 1]);
void Print_Optimal_Parens(int s[N + 1][N + 1], int i, int j);

int main()
{
	int p[N + 1] = { 30, 35, 15, 5, 10, 20, 25 };
	int m[N + 1][N + 1] = { 0 };
	int s[N + 1][N + 1] = { 0 };
	int i, j;

	Matrix_Chain(p, N + 1, m, s);
	cout << "m value is: " << endl;

	for (i = 1; i <= N; ++i)
	{
		for (j = 1; j <= N; ++j)
			cout << m[i][j] << " ";
		cout << endl;
	}

	cout << "s value is: " << endl;

	for (i = 1; i <= N; ++i)
	{
		for (j = 1; j <= N; ++j)
		{
			cout << s[i][j] << " ";
		}

		cout << endl;

	}

	cout << "The result is:" << endl;
	Print_Optimal_Parens(s, 1, N);
	return 0;

}

void Matrix_Chain(int *p, int len, int m[N + 1][N + 1], int s[N + 1][N + 1])
{
	int i, j, k, t;

	for (i = 0; i <= N; ++i)
	{
		m[i][i] = 0;
	}

	for (t = 2; t <= N; t++)                //当前链乘矩阵的长度
	{
		for (i = 1; i <= N - t + 1; i++)   //从第一矩阵开始算起，计算长度为t的最少代价
		{
			j = i + t - 1;                  //长度为t时候的最后一个元素
			m[i][j] = MaxValue;            //初始化为最大代价

			for (k = i; k <= j - 1; k++)   //寻找最优的k值，使得分成两部分k在i与j-1之间
			{
				int temp = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];

				if (temp < m[i][j])
				{
					m[i][j] = temp;         //记录下当前的最小代价
					s[i][j] = k;           //记录当前的括号位置，即矩阵的编号
				}

			}

		}

	}

}

//s中存放着括号当前的位置
void Print_Optimal_Parens(int s[N + 1][N + 1], int i, int j)
{
	if (i == j)
	{
		cout << "A" << i;
	}
	else
	{
		cout << "(";
		Print_Optimal_Parens(s, i, s[i][j]);
		Print_Optimal_Parens(s, s[i][j] + 1, j);
		cout << ")";
	}

}

111