其实这个过程不是很理解,书里的考究是使其能满足最优解,就是我们之前一直说的动态规划情况下满足的最优解。
以下部分内容引用:https://blog.csdn.net/c18219227162/article/details/50412333(感谢大神,只用于学习)
给定一个n个矩阵的序列(矩阵链)<A1,A2,...,An>,我们希望计算它们的乘积 A1A2...An,为了计算表达式,我们可以先用括号明确计算次序,然后利用标准的矩阵相乘算法进行计算。
例如如果有矩阵链为<A1,A2,A3,A4>,则共有5种完全括号化的矩阵乘积链。
(A1(A2(A3A4)))、(A1((A2A3)A4))、((A1A2)(A3A4))、((A1(A2A3))A4)、((A1A2)A3)A4)
对矩阵链加括号的方式会对乘积运算的代价产生巨大影响。我们先来分析两个矩阵相乘的代价。下面的伪代码的给出了两个矩阵相乘的标准算法,属性rows和columns是矩阵的行数和列数。
算法伪代码:
MATRIX-MULTIPLY(A,B)
1. if columns[A] !=row[B]
2. then error"矩阵乘法不相容"
3. else for i<-1 to rows[A]
4. do for j<- 1 to columns[B]
5. do C[i,j] <-0
6. for k<- 1 to columns[A]
7. do C[i,j]<-C[i,j]+A[i,k]*B[k,j]
8. return C
两个矩阵A和B只有相容(compatible),即A的列数等于B的行数时,才能相乘。如果A是p×q的矩阵,B是q×r的矩阵,那么乘积C是p×r的矩阵。计算C所需要时间由第8行的标量乘法的次数决定的,即pqr。以矩阵链<A1,A2,A3>为例,来说明不同的加括号方式会导致不同的计算代价。假设三个矩阵的规模分别为10×100、100×5和5×50。如果按照((A1A2)A3)的顺序计算,为计算A1A2(规模10×5),需要做10*100*5=5000次标量乘法,再与A3相乘又需要做10*5*50=2500次标量乘法,共需7500次标量乘法。
如果按照(A1(A2A3))的顺序计算,为计算A2A3(规模100×50),需100*5*50=25000次标量乘法,再与A1相乘又需10*100*50=50000次标量乘法,共需75000次标量乘法。因此第一种顺序计算要比第二种顺序计算快10倍。
矩阵链乘法问题(matrix-chain multiplication problem)可描述如下:给定n个矩阵的链<A1,A2,...,An>,矩阵Ai的规模为p(i-1)×p(i) (1<=i<=n),求完全括号化方案,使得计算乘积A1A2...An所需标量乘法次数最少。
步骤1:最优括号化方案的结构特征
动态规划的第一步是寻找最优子结构,然后就可以利用这种子结构从子问题的最优解构造出原问题的最优解。在矩阵链乘法问题中,我们假设A(i)A(i+1)...A(j)的最优括号方案的分割点在A(k)和A(k+1)之间。那么,继续对“前缀”子链A(i)A(i+1)..A(k)进行括号化时,我们应该直接采用独立求解它时所得的最优方案。
我们已经看到,一个非平凡(i≠j)的矩阵链乘法问题实例的任何解都需要划分链,而任何最优解都是由子问题实例的最优解构成的。为了构造一个矩阵链乘法问题实例的最优解,我们可以将问题划分为两个子问题(A(i)A(i+1)...A(k)和A(k+1)A(k+2)..A(j)的最优括号化问题),求出子问题实例的最优解,然后将子问题的最优解组合起来。我们必须保证在确定分割点时,已经考察了所有可能的划分点,这样就可以保证不会遗漏最优解。
步骤2:一个递归求解方案
下面用子问题的最优解来递归地定义原问题最优解的代价。对于矩阵链乘法问题,我们可以将对所有1<=i<=j<=n确定A(i)A(i+1)...A(j)的最小代价括号化方案作为子问题。令m[i,j]表示计算矩阵A(i..j)所需标量乘法次数的最小值,那么,原问题的最优解—计算A(1..n)所需的最低代价就是m[1,n]。
我们可以递归定义m[i,j]如下。对于i=j时的平凡问题,矩阵链只包含唯一的矩阵A(i..j)=A(i),因此不需要做任何标量乘法运算。所以,对所有i=1,2,...,n,m[i,i]=0。若i<j,我们利用步骤1中得到的最优子结构来计算m[i,j]。我们假设A(i)A(i+1)...A(j)的最优括号化方案的分割点在矩阵A(k)和A(k+1)之间,其中i<=k<j。那么,m[i,j]就等于计算A(i..k)和A(k+1..j)的代价加上两者相乘的代价的最小值。由于矩阵Ai的大小为p(i-1)*pi,易知A(i..k)和A(k+1..j)相乘的代价为p(i-1)p(k)p(j)次标量乘法运算。因此,我们得到
m[i,j]=m[i,k]+m[k+1,j]+ p(i-1)p(k)p(j)
此递归公式假定最优分割点k是已知的,但实际上我们是不知道。不过,k只有j-i种可能的取值,即k=i,i+1,...,j-1。由于最优分割点必在其中,我们只需检查所有可能情况,找到最优者即可。
因此,A(i)A(i+1)...A(j)的最小代价括号化方案的递归求解公式变为:
①如果i=j,m[i,j]=0
②如果i<j,m[i,j]=min{m[i,k]+m[k+1,j]+p(i-1)p(k)p(j)} i<=k<j
m[i,j]的值给出了子问题最优解的代价,但它并未提供足够的信息来构造最优解。为此,我们用s[i,j]保存最优括号化方案的分割点位置k,即使得m[i,j]=m[i,k]+[k+1,j]+p(i-1)p(k)p(j)成立的k值。
步骤3:计算最优代价
现在,我们可以很容易地基于递归公式写出一个递归算法,但递归算法是指数时间的,并不必检查若有括号化方案的暴力搜索方法更好。注意到,我们需要求解的不同子问题的数目是相对较少的:每对满足1<=i<=j<=n 的i和j对应一个唯一的子问题,共有n^2(最少)。递归算法会在递归调用树的不同分支中多次遇到同一个子问题。这种子问题重叠的性质是应用动态规划的另一标识(第一个标识是最优子结构)。
我们采用自底向上表格法代替递归算法来计算最优代价。此过程假定矩阵Ai的规模为p(i-1)×pi(i=1,2,...,n)。它的输入是一个序列p=<p0,p1,...,pn>,其长度为p.length=n+1。过程用一个辅助表m[1..n,1..n]来保存代价m[i,j],用另一个辅助表s[1..n-1,2..n](s[1,2]..s[n-1,n]这里i<j)记录最优值m[i,j]对应的分割点k。我们就可以利用表s构造最优解。
对于矩阵A(i)A(i+1)...A(j)最优括号化的子问题,我们认为其规模为链的长度j-i+1。因为j-i+1个矩阵链相乘的最优计算代价m[i,j]只依赖于那么少于j-i+1个矩阵链相乘的最优计算代价。因此,算法应该按长度递增的顺序求解矩阵链括号化问题,并按对应的顺序填写表m。
算法伪代码:
MATRIX-CHAIN-ORDER(p)
1. n<-length[p]-1
2. for i<-1 to n
3. do m[i,i]<-0
4. for l<-2 to n // l 为矩阵链的长度
5. do for i<-1 to n-l+1
6. do j<-i+l-1
7. m[i,j]<-∞
8. for k<-i to j-1
9. do q<-m[i,k]+m[k+1,j]+p(i-1)p(k)p(j)
10. if q<m[i,j]
11. then m[i,j]<-q //标量乘法运算数
12. s[i,j]<-k //分割点位置
13. return m and s
下图展示了对一个长度为6的矩阵链执行此算法的过程:
此图展示了将主对角线置于水平状态时的效果,计算子链AiAi+1...Aj的乘积的最小代价m[i,j]可在Ai行与Aj行的交汇处找到。即图中的顶点15125.
步骤4:构造最优解
虽然MATRIX_CHAIN_ORDER求出了计算矩阵链乘积所需的最少标量乘法运算次数,但它并未直接指出如何进行这种最优代价的矩阵链乘法计算。表s[i,j]记录了一个k值,指出A(i)A(i+1)...A(j)的最优括号化方案的分割点应在A(k)和A(k+1)之间。
因此,我们A(1..n)的最优计算方案中最后一次矩阵乘法运算应该是以s[1,n]为分界的A(1..s[1,n])*A(s[1,n]+1..n)。我们可以用相同的方法递归地求出更早的矩阵乘法的具体计算过程,因为s[1,s[1,n]]指出了计算A(1..s[1,n])时应进行的最后一次矩阵乘法运行;s[s[1,n]+1,n]指出了计算A(s[1,n]+1..n)时应进行的最后一次矩阵乘法运算。下面给出的递归过程可以输出<A(i),A(i+1),...,A(j)>的最优括号化方案。
算法伪代码:
PRINT-OPTIMAL-PARENS(s,i,j)
1. if i=j
2. then print "A"i
3. else print "("
4. PRINT-OPTIMAL-PARENS(s,i,s[i,j])
5. PRINT-OPTIMAL-PARENS(s,s[i,j]+1,j)
6. print ")"
C++:
matrixchain.h
#define _matrixchain_h
#include <stdlib.h>
#include<limits.h>
#include <stdio.h>
#include"pair.h"
pair matrixchainorder(int *p,int n){
int i,l,j,k,q;
int *m=(int*)malloc((n+1)*(n+1)*sizeof(int));
int *s=(int*)malloc((n+1)*(n+1)*sizeof(int));
for(i=0;i<=n;i++)
m[i*(n+1)+i]=0;
for(l=2;l<=n;l++)
for(i=1;i<=n-l+1;i++){
j=i+l-1;
m[i*(n+1)+j]=INT_MAX;
for(k=i;k<=j-1;k++){
q=m[i*(n+1)+k]+m[(k+1)*(n+1)+j]+p[i-1]*p[k]*p[j];
if(q<m[i*(n+1)+j]){
m[i*(n+1)+j]=q;
s[i*(n+1)+j]=k;
}
}
}
return make_pair(m,s);
}
void printoptimalparens(int *s,int i,int j,int n){
if(i==j)
printf("A%d",i);
else{
printf("(");
printoptimalparens(s,i,s[i*(n+1)+j],n);
printoptimalparens(s,s[i*(n+1)+j]+1,j,n);
printf(")");
}
}
pair.h
#define _pair_h //定义了一个结构体用来存放生成的表m以及s
typedef struct {
void *first;
void *second;
}pair;
pair make_pair(void *f,void *d){
pair p={f,d};
return p;
}
main.cpp
#include<stdio.h>
#include"matrixchain.h"
int main(){
int p[]={30,35,15,5,10,20,25};//p表示矩阵的行列数,如30行35列,35行15列,15行5列。。。
int *s,*m;
pair r=matrixchainorder(p,6);
m=(int*)r.first;
s=(int*)r.second;
printoptimalparens(s,1,6,6);
printf("\n%d\n",m[13]);//之所以m的下标是13的原因是,前面定义的i,j是加在一起的,[i*(n+1)+j]其中i为1,n为元素个数,j为6个矩阵链
free(s);//free()释放calloc, malloc, realloc动态分配的空间
free(m);
}
输出结果为:
C++:
Matrixchaincpp.h
#define _Matrixchaincpp_h
#include<iostream>
using namespace std;
pair<int*,int*> matrixchainorder(int *p,int n){
long i,l,j,k,q;
int *m=new int[(n+1)*(n+1)];
int *s=new int[(n+1)*(n+1)];
for(i=0;i<=n;i++)
m[i*(n+1)+i]=0;
for(l=2;l<=n;l++)
for(i=1;i<=n-l+1;i++){
j=i+l-1;
m[i*(n+1)+j]=numeric_limits<int>::max();
for(k=i;k<=j-1;k++){
q=m[i*(n+1)+k]+m[(k+1)*(n+1)+j]+p[i-1]*p[k]*p[j];
if(q<m[i*(n+1)+j]){
m[i*(n+1)+j]=q;
s[i*(n+1)+j]=k;
}
}
}
return make_pair(m,s);
}
void printoptimalparens(int *s,int i,int j,int n){
if (i==j)
cout<<'A'<<i;
else{
cout<<'(';
printoptimalparens(s,i,s[i*(n+1)+j],n);
printoptimalparens(s,s[i*(n+1)+j]+1,j,n);
cout<<')';
}
}
main.cpp
#include"Matrixchaincpp.h"
int main(){
int p[]={30,35,15,5,10,20,25};//p表示矩阵的行列数,如30行35列,35行15列,15行5列。。。
int *s,*m;
pair<int*,int*> r=matrixchainorder(p,6);//C++中默认预定义了数据结构类pair,可以直接使用
m=r.first;
s=r.second;
printoptimalparens(s,1,6,6);
cout<<endl<<m[13]<<endl;//之所以m的下标是13的原因是,前面定义的i,j是加在一起的,[i*(n+1)+j]其中i为1,n为元素个数,j为6个矩阵链
delete []s;//free()释放calloc, malloc, realloc动态分配的空间
delete []m;
}
输出结果如上图一样!
JAVA:
Matrixchain.java
package Jamin;
public class Matrixchain {
public static Pair matrixchainorder(int[] p) {//p为输入的矩阵行列数链
int n=p.length-1,i,l,j,k,q;
int [][] m=new int[n+1][n+1],s=new int[n+1][n+1];
for(i=0;i<=n;i++)
m[i][i]=0;
for(l=2;l<=n;l++)
for(i=1;i<=n-l+1;i++) {
j=i+l-1;
m[i][j]=Integer.MAX_VALUE;
for(k=i;k<=j-1;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;
}
}
}
return Pair.make(m, s);
}
public static void printoptimalparens(int[][] s,int i,int j) {
if(i==j)
System.out.printf("A%d",i);
else {
System.out.print("(");
printoptimalparens(s,i,s[i][j]);
printoptimalparens(s,s[i][j]+1,j);
System.out.print(")");
}
}
}
Pair.java
package Jamin;
public class Pair {//用于对m表和s表进行保存的结构变量
public Object first;
public Object second;
public Pair() {
first=new Object();
second=new Object();
}
public static Pair make(Object a,Object b) {
Pair p=new Pair();
p.first=a;
p.second=b;
return p;
}
}
Test.java
package Jamin;
public class Test {
public static void main(String[] args) {
// TODO Auto-generated method stub
int p[]= {30,35,15,5,10,20,25};
int[][] m,s;
Pair r=Matrixchain.matrixchainorder(p);
m=(int[][]) r.first;
s=(int[][]) r.second;
Matrixchain.printoptimalparens(s, 1, 6);
System.out.println();
System.out.println(m[1][6]);
}
}
输出结果:
((A1(A2A3))((A4A5)A6))
15125