动态规划之矩阵链乘
问题提出:(代码下载)
对于如下矩阵:
其中各矩阵A[i]下标为
计算其乘积的结果,以及我们需要计算其最小标量乘法次数。
问题分析:
首先我们需要明确的是何为标量:标量即为没有方向的量,而有方向的量即为矢量。(严谨的定义自己百度去)
那么标量乘法就变成了最基本的数字相乘。
其次对于两个矩阵相乘,需满足下示公式所示的形式:(左边矩阵的列数与右边矩阵的行数必须一致)
上述条件可从矩阵相乘的定义中看出:
在计算机,我们可以用一个二维数组来表示矩阵。
一个m行n列的矩阵与一个n行p列的矩阵相乘,会得到一个m行p列的矩阵,具体步骤如下:
其中第i行第j列位置上的数为第一个矩阵的第i行上的n个数与第二个矩阵第j列上的n个数对应相乘后所得的n个乘积之和。
下图展示了一个简单的例子:
从上图我们可以得出两个矩阵相乘城的次数,即为取得结果矩阵每个点所要做的标量乘法次数之和
而每个点的标量乘法次数为左边矩阵的列数n,结果矩阵有m*p个点,那么总数即为m*n*p;
另外我们也应该知道矩阵乘法符合乘法结合律,但不符合交换律;也就是说,我们可以为矩阵链乘设置括号,括号放置的位置不同,产生标量乘法的总量也是不一样的哦,这一点我们必须明确,对这点不清楚的,可以自己举个例子试试看。如果没有这点,那这题就没意义了。
前面科普了一大堆知识,现在进正题了:
对于矩阵链乘,我们可知道不同的括号化方案,会产生不同的计算代价;
算法导论书中指出了求完全括号化方案,所谓矩阵乘积链为完全括号化需满足如下性质:
要么是单一矩阵,要么是两个完全括号化的矩阵乘积链的积,且已加外括号。说实话这个定义让人有点摸不着头脑,待我仔细解释下:
首先:该乘积链中只有一个矩阵,即单一矩阵,该乘积链必是完全括号化,且该矩阵外部不用加括号。
再者:必须是两个括号化的矩阵链的乘积,我们可以把两个括号化的矩阵链看作两个结果矩阵(即两个单一矩阵),该两个矩阵相乘,并加上外括号把他们括住,即这种形式也叫着完全括号化的。
啰嗦了这么多,感觉简单的说法就是:单一矩阵(不加括号)是完全括号化,两个矩阵(单一矩阵)加上外括号也是完全括号化的,完全括号化的乘积链可以当作单一矩阵处理。
下面上一道课后的证明题:
对n个元素的表达式进行完全括号化时,恰好需要n-1对括号。
证明:
注:采用数学归纳法
已知仅有一个元素时,不需要括号,即0对括号。
假设k个元素时,需要k-1对括号;
当有k+1个元素时,我们可以看作k个元素的完全括号化的乘积链外加1个新加的元素。已知我们可以把完全括号化的乘积链看作单一矩阵,那么现有两个单一矩阵,那么需再加1对括号,才能再次完全括号化,而由假设知,k个元素的完全括号化需k-1对括号,那么k+1个元素的完全括号化方案就需要k-1+1对括号;
命题得证!
额,又科普了部分知识,只不过这些科普对后续理解问题都是有用的。
应用动态规划的方法的步骤:(出自书本)
1、刻画一个最优解的结构特征
2、递归地定义最优解的值
3、计算最优解的值,通常采用自底向上的方法
4、利用计算出的信息构造一个最优解的值
第一步:刻画一个最优解的结构特征
对于最开始提出的问题,我们看作为矩阵链乘设计最优的括号化方案。而设置括号的最终目的就是为了打乱计算顺序,通过结合律来找到最小计算代价,就比如下面一个问题33*4*25,是个正常人就会先计算后面两项,因为这样好计算,同样计算矩阵,我们也希望计算的少,计算的快啊,设置括号就变的尤为重要了。
对于完全括号化,每一个括号的作用是将一个链乘转化成两个完全括号化的链乘的积,前面有提到。
因此我最优解的结构特征就如下了:
对于第m个矩阵到第n个矩阵,他的最优解存在于对它的一次划分中,它的划分有n-m种情况,我们对这些情况中代价取最小,不就获得了最优解
而对于每种情况,我们需要划分后的两部分矩阵链都是最优解,乘积后才可能最小啊,(这部分大家可以用“剪切--粘贴”去反证),这样问题就转换成了两个子部分矩阵链的最优括号化方案的问题,以此不断的递归下去。
第二步:递归地定义最优解的值
我们的划分点k可从m到n-1,代价如下所示,每次的代价如下,而最优解只有一个,即它们的最小值:
其一般化的公式如下:
从这里可以看出,这是一个递归的问题,看着很类似上篇所述的钢条切割问题啊,仔细比对比对吧。
若采用递归策略,可能会有很多的重复子问题,无疑提高了计算的复杂度。也可以采取带备忘的递归来降低复杂度。
第三步:计算最优解的值,通常采用自底向上的方法
另外从公式中可以看出,上述计算结果要依赖于m,n(m<n)的所有可能组合。区间[m,n]组合要依赖于它们之间的子区间组合。所以我们采取自下而上的策略,逐渐扩大m,n之间的区间长度,最终算到m=0,n=n,从而获得最优的括号化方案,
具体方法为:matrix_chain_mutiply.cpp中的void matrix_divide_option(int *num,int** cost,int** point,int length)方法。
第四步、利用计算出的信息构造一个最优解的值
具体解析见matrix_chain_mutiply.cpp最后两个方法
代码如下:
matrix.h文件:
#include<iostream> #ifndef MATRIX_H #define MATRIX_H /* 矩阵类 row为其行数 col为该矩阵的列数, value表示一个二维数组,主要存储该矩阵各点的值 */ class matrix { public: matrix(); //默认的构造函数 matrix(int row,int col); //给定行数和列数,构造一个矩阵 virtual ~matrix(); //析构函数 matrix(const matrix& other); //拷贝构造函数 matrix operator*(const matrix& a); //矩阵的乘法 matrix& operator=(const matrix&); //重写赋值函数 matrix init(); //初始化矩阵 friend std::ostream& operator << (std::ostream&,const matrix &); //重写输出 protected: private: double **value; //主要用于存储矩阵各点的值 int col; //列数 int row; //行数 }; #endif // MATRIX_H
matrix.cpp
#include<iostream> #include "matrix.h" #include<cstdlib> #include<ctime> matrix::matrix() //默认构造函数 { col = row = 0; value = nullptr; } matrix::matrix(int row, int col) //根据给定行,列数来初始化一个矩阵,并分配相应的内存 { this->row = row; this->col = col; value = new double*[row]; for(int i=0;i<row;i++) { value[i] = new double[col]; } } matrix::matrix(const matrix& other) //复制构造函数,注意此处要用matrix&,重写复制构造函数为了防止值复制造成 { row = other.row; col = other.col; value = new double*[row]; //此处使用new分配内存 for(int i=0;i<row;i++) { value[i] = new double[col]; } for(int i=0;i<row;i++) { for(int j=0;j<col;j++) { value[i][j]=(other.value)[i][j]; } } } matrix::~matrix() { for(int i=0;i<row;i++) //因为构造函数中new分配有内存,所以在析构函数中要对内存进行释放。 { delete[] value[i]; } delete[] value; } matrix& matrix::operator=(const matrix& other) //重写了赋值函数 { if(this == &other) //防止a=a的情况出现 { return *this; } for(int i=0;i<row;i++) //先释放原有区域的内存,从此处我们可以上面的判断的重要性 { //若没有上述判断,a=a在赋值的同时,也将内存释放了,造成数据丢失 delete[] value[i]; } delete[] value; row = other.row; //将成员变量进行复制 col = other.col; value = new double*[row]; //重新分配内存 for(int i=0;i<row;i++) { value[i] = new double[col]; } for(int i=0;i<row;i++) //对二维数组进行复制 { for(int j=0;j<col;j++) { value[i][j]=(other.value)[i][j]; } } return *this; } matrix matrix::operator*(const matrix& n) //矩阵的乘法 { matrix result(row,n.col); for(int i=0;i<row;i++) { for(int j=0;j<n.col;j++) { (result.value)[i][j] = 0; for(int k=0;k<n.row;k++) { (result.value)[i][j] += value[i][k]*(n.value)[k][j]; } } } return result; } matrix matrix::init() //初始化矩阵 { matrix result(row,col); srand((unsigned)time(NULL)); for(int i=0;i<row;i++) { for(int j=0;j<col;j++) { (result.value)[i][j] = rand()%100; } } return result; } std::ostream& operator <<(std::ostream& os,const matrix& m)//将矩阵输出 { for(int i=0;i<m.row;i++) { for(int j=0;j<m.col;j++) { os << (m.value)[i][j] << '\t'; } os << std::endl; } return os; }
matrix_chain_mutiply.cpp(核心)
#include<iostream> #include "matrix.h" const int SIZE = 5; //SIZE的大小决定了矩阵的数量(SIZE-1),可自行改变 const int max_size = 2147483647; void matrix_divide_option(int *,int **,int **,int); //确定最优的括号化方案 void show(int **); //输出二维数组 void print(int**,int,int); //输出矩阵链乘的括号化方案 matrix matrix_mutiply(int**,int,int,matrix*); //获得矩阵链乘的最终结果矩阵 int main() { using namespace std; cout << "input " << SIZE << " number(press enter to end your input):" << endl; int *num = new int[SIZE]; //用于存储矩阵A[i]下标为num[i]*num[i+1]; int i=0; cin >> num[i]; //获得矩阵的行列数 while(cin.get() != '\n') { i++; cin >> num[i]; //读取数组维数; } cout << endl << "what you input is:" << endl; for(int j=0;j<SIZE;j++) { cout << num[j] << "-----"; //展示输入的是否正确 } cout << endl << endl; matrix* test = new matrix[SIZE-1]; //为矩阵数组分配内存 for(int j=1;j<SIZE;j++) { matrix m(num[j-1],num[j]); //初始化矩阵大小 test[j-1] = m.init(); //初始化该矩阵各点的值 } cout << "矩阵信息如下:" << endl; for(int j=0;j<SIZE-1;j++) { cout << test[j] << endl; //输出矩阵信息 } int **cost = new int*[SIZE-1]; //存储矩阵相乘时,各标量乘法次数,即代价表 int **point = new int*[SIZE-1]; //存储矩阵括号化方案时的分割点 for(int j=0;j<SIZE-1;j++) { cost[j] = new int[SIZE-1]; //分配内存 point[j] = new int[SIZE-1]; } matrix_divide_option(num,cost,point,SIZE-1); //获取划分方案,获得最小的代价 cout << "各区间的计算代价:" << endl; show(cost); cout << endl; cout << "区间划分点所在的位置:" << endl; show(point); cout << endl; cout << "输出完全括号化方案:" << endl; print(point,0,SIZE-2); cout << endl << endl; cout << "输出结果矩阵:" << endl; matrix result = matrix_mutiply(point,0,SIZE-2,test); cout << result; for(int j=0;j<SIZE-1;j++) { delete[] cost[j]; delete[] point[j]; } delete[] cost; delete[] point; delete[] test; delete[] num; return 0; } /* 此方法主要获得计算代价,以及括号化方案的划分点 为此问题解决的核心方法。 num存储的矩阵的下标,num[i-1]为第i个矩阵的行数,num[i]表示第i个矩阵的列数 因为数组存储时,最初始的下标为0,所以下述的定义就变的有点不好理解!仔细想想也容易明白。 cost----二维数组,cost[i][j]代表的就是第(i+1)个矩阵到第(j+1)个矩阵的链乘的最小计算代价 point----二维数组,point[i][j]代表的就是第(i+1)个矩阵到第(j+1)个矩阵的链乘的最佳划分处,假如为k,则k就是第k+1个矩阵, (i+1)到(k+1)矩阵属于左边的完全括号化矩阵,其余为右侧完全括号化矩阵 length代表的是矩阵的个数,即num的长度-1 */ void matrix_divide_option(int *num,int** cost,int** point,int length) { //首先将代价表中的值重置 //下标i到下标j之间的矩阵链乘的最小值 for(int i=0;i<length;i++) { for(int j=0;j<length;j++) { if(i>=j) { cost[i][j] = 0; //单个矩阵不存在乘积,代价0,同时也将不合理的下标组合置为0 } else { cost[i][j] = max_size; //先将各合理的组合代价设为最大,以便后续比较使用。 } point[i][j] = -1; //将所有组合的划分点置为-1,以便后续记录。 } } /* l代表矩阵链中矩阵的个数,我们这里从2个矩阵起步.由公式可知,我们要算遍所有的合理组合,同时大区间的组合 还要依赖于其子区间的合理组合,因为完全括号化的性质,最终都归根于2个矩阵的乘积,所以我们从2个矩阵起步, 逐步向上计算,最终就能计算(0,n)的组合的最优括号化的代价。 */ for(int l=2;l<=length;l++) { for(int i=0;i<=length-l;i++) //此处要注意的i+l-1不能超过数组的个数 { int j = i+l-1; //矩阵个数要为l,则j-i需满足j-i+1=l,所以有该等式 int q = 0; for(int k=i;k<j;k++) //注意k要在区间内移动,利用下述if条件,获得所有可能值中的最小值, { q = cost[i][k] + cost[k+1][j] + num[i]*num[k+1]*num[j+1]; if(q < cost[i][j]) { cost[i][j] = q; point[i][j] = k; } } } } } void show(int **p) { using namespace std; for(int i=0;i<SIZE-1;i++) { for(int j=0;j<SIZE-1;j++) { cout << p[i][j] << '\t'; } cout << endl; } } void print(int** point,int i,int j) { if(i == j) { std::cout << "A" << i+1; } else { std::cout << "("; print(point,i,point[i][j]); //输出划分点左侧的完全括号化表达式 print(point,point[i][j]+1,j); //输出划分点右侧的完全括号化表达式 std::cout << ")"; } } /*矩阵链乘的方法与括号化方案表达式输出的方法是同样的道理*/ matrix matrix_mutiply(int** point,int i,int j,matrix* m) { if(i==j) return m[i]; else { matrix m1 = matrix_mutiply(point,i,point[i][j],m); matrix m2 = matrix_mutiply(point,point[i][j]+1,j,m); return m1*m2; } }
该问题的解决代码在上述三个文件中。
同时上述代码也存在一定的问题,关于内存是否泄漏以及代码是否符合规范,本人不是非常熟悉。望各位大神,看后觉得有什么问题,直接回复,我也希望提高,毕竟code能力有待提高。
以上问题也属于典型的动态规划问题,也属于递归向迭代的转换。
总结:
应用动态规划的方法的步骤:
1、刻画一个最优解的结构特征
2、递归地定义最优解的值
3、计算最优解的值,通常采用自底向上的方法
4、利用计算出的信息构造一个最优解的值
源代码文件。。。。。。。点击这里