一、思路
1. 数学原理
矩阵相乘时不同的结合顺序,虽然结果都一样,但会导致不同的cost(标量乘法代价)。
可以从其得出规律:
- 不用考虑A1,A2,…,An的Ai的前后顺序。意思是说不用考虑这种情况
A1 x A3 x A2
(违背矩阵乘法的基本原则)。 - 只用考虑先结合谁,再结合谁的结合顺序。而且,大的结合顺序可以划分成两个更小的结合顺序,如
A1 x A3
划分成了A1
和A2 x A3
这两个子结合顺序。所以,整体的代价就是两个子结合顺序内部的代价,和将两个子部分再结合的代价。如,A1 x (A2 x A3)
分别是0
,1000
,200
。
那么如何知道该怎么结合,答案很简单,用for循环暴力试。
2. 核心公式
符号:
-
A i A_i Ai:第i个矩阵,从 A 1 A_1 A1开始。
-
m[i,j]
:表示子结合顺序内部的即 A i × A i + 1 . . . × A j A_i\times A_{i+1} ... \times A_j Ai×Ai+1...×Aj的代价,而且是 最优 代价。其中,m[i,i]
表示 A i × A i A_i \times A_i Ai×Ai,代价为0。
如:A1 x (A2 x A3)
的最优总代价就是m[1,3]。 -
k
:表示划分,左边的子结合顺序是到Ak
,右边的子结合顺序从Ak+1
开始。
i ≤ k < j i \leq k < j i≤k<j这个等于不等于是根据 m i n . . . min{...} min...那个公式来划分的。- k ≥ i k \geq i k≥i是因为当k=i时,因为可以有m[i,i]。
-
k
<
j
k < j
k<j是因为
m[k+1,j]
,k+1必须小于等于j。
-
p
:表示矩阵的维度描述,用n+1个数描述n个矩阵,用来直接计算代价。如:A1 , A2 , A3
,就是p=<10,20,50,1>
。
从0开始,A1 x A2
是p0 x p1 x p2
-
p i − 1 p k p j p_{i-1}p_{k}p_{j} pi−1pkpj:表示将两个子部分 ( A i × . . . × A k ) × ( A k + 1 × . . . × A j ) (A_i \times ... \times A_k)\times(A_{k+1} \times ... \times A_j) (Ai×...×Ak)×(Ak+1×...×Aj)再结合的代价。
如 p 0 p 1 p 3 p_0p_1p_3 p0p1p3就是 ( A 1 ) ( A 2 A 3 ) (A_1)(A_2A_3) (A1)(A2A3), p 0 p 2 p 3 p_0p_2p_3 p0p2p3就是 ( A 1 A 2 ) ( A 3 ) (A_1A_2)(A_3) (A1A2)(A3)。
3. 图示
【两个右上的三角表格】:
- m表:表示代价。n行n列。
- s表:表示k。因为m[i,i]不需要也不能划分,所以去掉最外面这条对角线,n-1行n-1列。
使用动态程序填写表m[i,j]
:
- 首先将
m[i,i]
=0,其中i=1,…,n。(对应单个的矩阵代价,如例A1 x (A2 x A3)
划分的A1代价) - 然后计算
m[1,2]
,m[2,3]
,…,m[n-1,n]
。(对应最小结合代价,两个矩阵的结合代价) - 然后
m[1,3]
,m[2,4]
,…,m[n-2,n]
,…(对应三个矩阵的结合总代价) - …一直到我们可以计算
m[1,n]
。(这就是从下到上的体现)
PS1:计算m[i,j]的k展开形状:
PS:如何从S推导出矩阵的结合顺序?
要划分Ai到Aj,就去查s[i,j],根据其值k划分为(Ai...Ak)(Ak+1...Aj)
。
例题:
p=<30,35,15,5,10,20,25>
-
m
[
1
,
2
]
=
m
[
1
,
1
]
+
m
[
2
,
2
]
+
p
0
p
1
p
2
=
0
+
0
+
(
A
1
)
(
A
2
)
=
30
×
35
×
15
=
15750
m[1,2]=m[1,1]+m[2,2]+p_0p_1p_2=0+0+(A_1)(A_2)=30\times35\times15=15750
m[1,2]=m[1,1]+m[2,2]+p0p1p2=0+0+(A1)(A2)=30×35×15=15750
k=1,s[1,2]=1。 -
m
[
1
,
3
]
m[1,3]
m[1,3]
- k=1时,即 ( A 1 ) ( A 2 A 3 ) (A_1)(A_2A_3) (A1)(A2A3): m [ 1 , 3 ] = m [ 1 , 1 ] + m [ 2 , 3 ] + p 0 p 1 p 3 = 0 + 2625 + 30 × 35 × 5 = 7875 m[1,3]=m[1,1]+m[2,3]+p_0p_1p_3=0+2625+30\times35\times5=7875 m[1,3]=m[1,1]+m[2,3]+p0p1p3=0+2625+30×35×5=7875
- k=2时,即
(
A
1
A
2
)
(
A
3
)
(A_1A_2)(A_3)
(A1A2)(A3):
m
[
1
,
3
]
=
m
[
1
,
2
]
+
m
[
3
,
3
]
+
p
0
p
2
p
3
=
15750
+
0
+
30
×
15
×
5
=
18000
m[1,3]=m[1,2]+m[3,3]+p_0p_2p_3=15750+0+30\times15\times5=18000
m[1,3]=m[1,2]+m[3,3]+p0p2p3=15750+0+30×15×5=18000
所以,选择k=1,s[1.3]=1。
4. 伪代码
// p比如<30,35,15>,表示A1(30*35)和A2(35*15)
MATRIX-CHAIN-ORDER(p):
// p用n+1个数表示n个矩阵相乘
n ← length[p] - 1
// 对表中第一条对角线初始化当i=j的情况
for i ← 1 to n:
do m[i, i] ← 0
// 控制每条红色对角线的行号区间递减,每条红色对角线列号在随行号递增的同时随对角线向上加一个数递增
for l ← 2 to n:
// 每条红色对角线的行号区间,[1,n-1],[1,n-2],...,[1,1]
do for i ← 1 to n - l + 1:
// 每条红色对角线的每行对应的列号
do j ← i + l - 1
// 先初始化正无穷
m[i, j] ← ∞
// 遍历k,i≤k<j,即[i,j-1]:找出最小的m[i,j]
for k ← i to j - 1:
// 公式
do q ← m[i, k] + m[k + 1, j] + (pi-1)(pk)(pj)
// 比较新试的k和原来哪个更小
if q < m[i, j]
then m[i, j] ← q
s[i, j] ← k
二、cpp实现
#include <iostream>
#include <vector>
using namespace std;
#define INF 0x3f3f3f3f
void matrix_chain_order(vector<int> &p, vector<vector<int>> &m, vector<vector<int>> &s)
{
// p用n+1个数表示n个矩阵相乘
int n = p.size() - 1;
// 对表中第一条对角线初始化当i=j的情况
// 表行列总数都是n个矩阵个,[1,n]
for (int i = 1; i <= n; i++)
{
m[i][i] = 0;
}
// 控制每条红色对角线的行号区间递减,每条红色对角线列号在随行号递增的同时随对角线向上加一个数递增
for (int l = 2; l <= n; l++)
{
// 每条红色对角线的行号区间,[1,n-1],[1,n-2],...,[1,1]
for (int i = 1; i <= n - l + 1; i++)
{
// 每条红色对角线的每行对应的列号
int j = i + l - 1;
// 先初始化正无穷
m[i][j] = INF;
// 遍历k,i≤k<j,即[i,j-1]:找出最小的m[i,j]
for (int k = i; k <= j - 1; k++)
{
// 公式,表示新试k的代价
int q = m[i][k] + m[k + 1][j] + p[i - 1] * p[k] * p[j];
// 比较新试的k和原来哪个更小
if (q < m[i][j])
{
m[i][j] = q;
s[i][j] = k;
}
}
}
}
}
// void print_optimal_parens(vector<vector<int>> &s, int i, int j)
// {
// if (i == j)
// {
// printf("A%d\n", i);
// }
// else
// {
// printf("(");
// print_optimal_parens(s, i, s[i, j]);
// print_optimal_parens(s, s[i, j] + 1, j);
// printf(")\n");
// }
// }
int main(void)
{
// 矩阵链,7个数表示6个矩阵
vector<int> p{30, 35, 15, 5, 10, 20, 25};
// 这里我们不从0开始,从1开始使用,[1,6]范围
// m,表示记录结果
vector<vector<int>> m(p.size(), vector<int>(p.size()));
// s,表示选择的k
vector<vector<int>> s(p.size(), vector<int>(p.size()));
// 调用
matrix_chain_order(p, m, s);
// 输出m
for (int i = 1; i < p.size(); i++)
{
for (int j = 1; j < p.size(); j++)
{
printf("%5d ", m[i][j]);
}
printf("\n");
}
printf("\n");
// 输出s
for (int i = 1; i < p.size(); i++)
{
for (int j = 1; j < p.size(); j++)
{
printf("%5d ", s[i][j]);
}
printf("\n");
}
return 0;
}