矩阵乘法
一个矩阵的列数等于另一个矩阵的行数,两个矩阵可以相乘。例如,一个 1×2的矩阵A 和一个 2×3的矩阵B相乘,得到一个1×3的矩阵,运算次数为1×2×3=6,同理,一个4×10的矩阵和一个10×35的矩阵相乘,运算次数为4×10×35=1400。
多个矩阵相乘时,运算顺序不同,需要的运算次数也不同。例如矩阵1×3的矩阵A1、3×15的矩阵A2、15×4的矩阵A3相乘,有两种情况:
1.(A1*A2)A3,先计算A1*A2的次数1×3×15=45,得到一个1×15的矩阵,再加上新矩阵乘A3的次数1×15×4=60,最后得到45+60=105.
2.A1*(A2*A3),先计算A2*A3的次数3×15×4=180,再计算1×3×4=12,最后得到180+12=192.
我们要做的事情是找到一个最优的运算顺序,使得多个矩阵相乘时运算次数最少。
矩阵维数序列
若有一个矩阵维数序列 seq[] 存储5,10,15,30,这表示三个矩阵,其中A1为5×10,A2为10×15,A3为15×30. 所以A1的行数为seq[0],A1的列数为seq[1],同样地,A2的行数为seq[1],A2的列数为seq[2],这样矩阵Ai的行数和列数在seq[]中分别对应seq[i-1]和seq[i]. 若矩阵从A0开始,则Ai的行数和列数分别为seq[i]和seq[i+1]。
动态规划
动态规划问题有两个特点,一个是最优子结构,另一个是子问题重叠。对应多个矩阵连乘最优顺序的问题,子结构就是少数几个矩阵相乘的最优顺序。
最优子结构:即若𝐴1 ∗ 𝐴2 ∗ ⋯ ∗ 𝐴𝑛的一个最优计算顺序从第𝑘个矩阵处断开,即分为𝐴1 ∗ 𝐴2 ∗ ⋯ ∗ 𝐴𝑘 和𝐴𝑘 + 1 ∗ 𝐴𝑘 + 2 ∗ ⋯ ∗ 𝐴𝑛两个子问题,则该最优解应该包含𝐴1 ∗ 𝐴2 ∗ ⋯ ∗ 𝐴𝑘的一个 最优计算顺序和𝐴𝑘 + 1 ∗ 𝐴𝑘 + 2 ∗ ⋯ ∗ 𝐴𝑛的一个最优计算顺序。最后的结果应当是左边一堆的计算次数+右边一堆的计算次数+左右两边各自得到的新矩阵相乘的计算次数。
构造递归式:
我们需要一个二维数组 cost[i][j] 来记录从矩阵 Ai 到 Aj 的最优计算次数,k为分割的位置,k的位置应当大于等于 i 且小于 j。
第一种情况 i=j ,即从当前矩阵连乘到当前矩阵的最优顺序,一个矩阵不需要相乘,此时不存在一个合适的k。cost[i][j]=0.
第二种情况 i<j ,即从矩阵 Ai 到后面一个矩阵 Aj 的最优相乘顺序,我们要找到一个k将原来的矩阵分成两部分,一部分是Ai到Ak,另一部分是Ak+1 到Aj,cost[i][j] 应该为k值变化时,左边计算次数cost[i][k]+右边计算次数cost[k+1][j]+最后左右两边相乘的次数的最小值。得到如下递归式:
其中最后的pi*pk+1*pj+1是根据矩阵维数序列求得的,其实就是Ai的行数*Ak的行数*Aj的列数(Ak的行数和Ai的列数相同,Ak的列数和Aj的行数相同)。
动态规划解决矩阵连乘问题
根据上面的分析,我们已经知道了二维数组cost[i][j]中存放的是Ai到Aj的最少计算次数,且当i=j时数组内容为0.
j=1 | j=2 | j=3 | j=4 | |
---|---|---|---|---|
i=1 | 0 | |||
i=2 | 0 | |||
i=3 | 0 | |||
i=4 | 0 |
其中表格的下三角是空的,因为我们不考虑从A3到A1(即i=3,j=1)这种情况,只考虑从左边往右边乘,例如 i=2,j=4 表示从A2到A4的最少计算次数,或 i=1,j=3 表示从A1到A3的最少计算次数.
我们还需要一个二维数组trace[ ][ ]来记录每种情况对应的 k 的位置。
以该题为例,矩阵维数序列为30,35,15,5,10,20,25。
先来求cost[1][2],因为k的范围是i<=k<j,此时k只能取1,即从第一个矩阵后面断开。根据递归式,
cost[1][2]=cost[1][1]+cost[2][2]+seq[0]*seq[1]*seq[2]=0+0+30*35*15=15750
而对于cost[2][5],k的取值为2或3或4,共三种情况,根据递归式分别计算三种情况,取最小值存入cost[2][5] ,并将k值存入对应的trace[2][5].
代码
/*多个矩阵连乘,动态规划法求最少次数*/
#define N 100
int cost[N][N]; //记录最少次数
int trace[N][N]; //记录每种情况的分割位置k
void martixMultiply(int n, int seq[]){
int tempCost; //本次计算次数
int tempTrace; //本次划分位置k
int i, j, k, p; //p为子问题的规模
//p=1表示一个矩阵相乘,p=2表示两个矩阵相乘
int temp; //当前最少计算次数
for (i = 0; i < n; i ++ ) //n个矩阵相乘
cost[i][i] = 0; //将对角线赋为0,即i=j的情况,如上表
for (p = 1; p < n; p ++ ) {
/*第一层循环,p从1开始,即先考虑分成一个一个矩阵相乘的情况
p=2即考虑两个两个相乘,p=3即为三个矩阵相乘求最优顺序 */
for (i = 0; i < n - p; i ++ ) {
/*第二层循环,从第一个矩阵开始,到第n-p个结束,因为i+p等于j,
也就是开始的Ai加上问题的规模p就是结束的Aj*/
j = i+p;
tempCost = -1; //先把本次计算次数设为-1,后续更新
for (k = i; k<j ; k ++ ) {
/*第三层循环,对应每种情况下k的可能的取值进行循环,
k的范围为大于等于i小于j*/
temp = cost[i][k] + cost[k + 1][j] + seq[i] * seq[k + 1] * seq[j + 1];
//这里seq[]的下标要根据实际写,也可能是i-1,k,j
if(tempCost == -1 || tempCost > temp) { //更新当前情况的最少次数
tempCost = temp;
tempTrace = k; //记录分割位置
}
cost[i][j] = tempCost;
trace[i][j] = tempTrace;
}
}
return cost[0][n - 1]; //最后返回二维数组cost
}
这里以软考的代码为例,其中问题规模p是从1开始的,很多代码从2开始(因为一个矩阵不存在乘),思路是一样的。之前看了很多代码都不明白,主要是注释太少了所以不知道三层循环分别做了什么,循环的终止条件也很重要。该代码没有来得及运行,后续可能会增加一个完整代码。