笔记
给定两个
n
×
n
n×n
n×n正方矩阵
A
A
A和
B
B
B,这两个矩阵的乘法定义为
其中
下面是矩阵乘法的伪代码。
很显然,执行SQUARE-MATRIX-MULTIPLY需要花费
Θ
(
n
3
)
Θ(n^3)
Θ(n3)时间。然而,有一种方法可以花费更少的时间,这就是Strassen算法,它本质上也是一种分治法,它的时间复杂度为
Θ
(
n
l
g
7
)
=
O
(
n
2.81
)
Θ(n^{{\rm lg}7}) = O(n^{2.81})
Θ(nlg7)=O(n2.81)。
在介绍Strassen算法之前,我们先尝试简单的分治法来计算矩阵乘法
C
=
A
•
B
C = A•B
C=A•B。假定三个矩阵均为
n
×
n
n×n
n×n正方矩阵。并且为简化分析,假定
n
n
n为
2
2
2的幂。我们将
A
、
B
A、B
A、B和
C
C
C均分解为
4
4
4个
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)的子矩阵。
于是矩阵乘法可以表示为
上面的矩阵乘法等价于下面
4
4
4个式子。
上面每个式子都对应
2
2
2次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵乘法,以及
1
1
1次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵加法。根据以上分析,可以给出一个递归的分治算法。
现在分析这个简单分治法的时间复杂度。调用SQUARE-MATRIX-MULTIPLY-RECURSIVE计算两个
n
×
n
n×n
n×n矩阵乘法的运行时间用
T
(
n
)
T(n)
T(n)表示。对于
n
=
1
n = 1
n=1的初始情况,我们只需计算一次标量乘法,因此
T
(
1
)
=
Θ
(
1
)
T(1) = Θ(1)
T(1)=Θ(1)。当
n
>
1
n > 1
n>1时,根据上面的伪代码,
T
(
n
)
T(n)
T(n)包含
8
8
8次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵乘法的时间和
4
4
4次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵加法的时间,所以
T
(
n
)
=
8
T
(
n
/
2
)
+
Θ
(
n
2
)
T(n) = 8T(n/2) + Θ(n^2)
T(n)=8T(n/2)+Θ(n2),这里忽略了分解子矩阵的时间。于是,我们得到SQUARE-MATRIX-MULTIPLY-RECURSIVE的运行时间的递时式为
求解这个递归式得到
T
(
n
)
=
Θ
(
n
3
)
T(n) = Θ(n^3)
T(n)=Θ(n3)。可以看到,这个简单的分治法并没有带来渐近运行时间的提升。
下面介绍Strassen算法。Strassen算法同样要将每个矩阵分解为
4
4
4个
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)子矩阵。而与简单分治法不同,Strassen算法只需要递归为
7
7
7次,而不是
8
8
8次。下面直接给出Strassen算法的流程。
(1) 将输入矩阵
A
、
B
A、B
A、B以及输出矩阵
C
C
C各分解为
4
4
4个
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)子矩阵。
(2) 创建
10
10
10个
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵
S
1
,
S
2
,
…
,
S
10
S_1, S_2, …, S_{10}
S1,S2,…,S10,如下所示。由于需要进行
10
10
10次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵的加减法,所以这一步花费
Θ
(
n
2
)
Θ(n^2)
Θ(n2)时间。
(3) 用步骤(1)分解得到的子矩阵和步骤(2)中创建的
10
10
10个矩阵,递归地计算
7
7
7个矩阵乘积
P
1
,
P
2
,
…
,
P
7
P_1, P_2, …, P_7
P1,P2,…,P7,如下所示。
(4) 利用矩阵
P
1
,
P
2
,
…
,
P
7
P_1, P_2, …, P_7
P1,P2,…,P7进行加减运算,得到输出矩阵
C
C
C的子矩阵
C
11
,
C
12
,
C
21
,
C
22
C_{11}, C_{12}, C_{21}, C_{22}
C11,C12,C21,C22,如下所示。这一步需要进行
8
8
8次
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)矩阵的加减法,所以花费时为
Θ
(
n
2
)
Θ(n^2)
Θ(n2)。
由于Strassen算法只需要递归为
7
7
7次,因此它的运行时间的递归式为
求解这个递归式,可以得到Strassen算法的运行时间
T
(
n
)
=
Θ
(
n
l
g
7
)
T(n) = Θ(n^{{\rm lg}7})
T(n)=Θ(nlg7)。
练习
4.2-1 使用Strassen算法计算如下矩阵乘法:
给出计算过程。
解
(1) 分解输入矩阵
(2) 计算矩阵
S
1
,
S
2
,
…
,
S
10
S_1, S_2, …, S_{10}
S1,S2,…,S10
(3) 计算矩阵
P
1
,
P
2
,
…
,
P
7
P_1, P_2, …, P_7
P1,P2,…,P7
(4) 计算输出矩阵的
4
4
4个子矩阵
最终结果为
4.2-2 为Strassen算法编写伪代码。
解
这里还是假设了矩阵的宽高
n
n
n为
2
2
2的幂。下面给出伪代码。
4.2-3 如何修改Strassen算法,使之适应矩阵规模
n
n
n不是
2
2
2的幂的情况?证明:算法的运行时间为
Θ
(
n
l
g
7
)
Θ(n_{{\rm lg}7})
Θ(nlg7)。
解
为了保证算法的通用性,需要考虑矩阵的宽高
n
n
n不为
2
2
2的幂的情况。分两种情况讨论。
(1)
n
n
n为偶数
这种情况下
n
×
n
n×n
n×n矩阵可以直接分解为
4
4
4个
(
n
/
2
)
×
(
n
/
2
)
(n/2)×(n/2)
(n/2)×(n/2)的子矩阵,因此可以直接应用Strassen算法。为了计算矩阵乘法
C
n
×
n
=
A
n
×
n
•
B
n
×
n
C_{n×n} = A_{n×n}•B_{n×n}
Cn×n=An×n•Bn×n,令
m
=
n
/
2
m = n/2
m=n/2,需要将矩阵分解为
这种情况下,矩阵乘法所花费的时间
T
(
n
)
=
7
T
(
n
/
2
)
+
Θ
(
n
2
)
T(n) = 7T(n/2) + Θ(n^2)
T(n)=7T(n/2)+Θ(n2)。
(2)
n
n
n为奇数
这种情况不能直接应用Strassen算法。为了计算矩阵乘法
C
n
×
n
=
A
n
×
n
•
B
n
×
n
C_{n×n} = A_{n×n}•B_{n×n}
Cn×n=An×n•Bn×n,令
m
=
n
−
1
m = n−1
m=n−1,将矩阵做如下分解
如上所示,每个
n
×
n
n×n
n×n矩阵被分解为一个
(
n
−
1
)
×
(
n
−
1
)
(n−1)×(n−1)
(n−1)×(n−1)矩阵、一个
(
n
−
1
)
×
1
(n−1)×1
(n−1)×1矩阵、一个
1
×
(
n
−
1
)
1×(n−1)
1×(n−1)矩阵和一个
1
×
1
1×1
1×1矩阵。相应地,矩阵乘法
C
n
×
n
=
A
n
×
n
•
B
n
×
n
C_{n×n} = A_{n×n}•B_{n×n}
Cn×n=An×n•Bn×n可以分解为下面
4
4
4个式子。
上面4个式子包含了8个不同规模的矩阵乘法,下面逐个进行分析。
1)
A
1
1
m
×
m
•
B
1
1
m
×
m
A11_{m×m}•B11_{m×m}
A11m×m•B11m×m:由于
m
=
n
−
1
m = n−1
m=n−1是偶数,所以这个矩阵乘法可以直接应用Strassen算法。
这一矩阵乘法所花费的时间为
T
(
n
−
1
)
=
7
T
(
(
n
−
1
)
/
2
)
+
Θ
(
(
n
−
1
)
2
)
=
7
T
(
⌊
n
/
2
⌋
)
+
Θ
(
n
2
)
T(n-1)=7T((n-1)/2)+Θ((n-1)^2)=7T(⌊n/2⌋)+Θ(n^2)
T(n−1)=7T((n−1)/2)+Θ((n−1)2)=7T(⌊n/2⌋)+Θ(n2)。
2)
A
1
2
m
×
1
•
B
2
1
1
×
m
A12_{m×1}•B21_{1×m}
A12m×1•B211×m:采用朴素算法,需要做
(
n
−
1
)
2
(n−1)^2
(n−1)2次乘法,因此运行时间为
Θ
(
n
2
)
Θ(n^2)
Θ(n2)。
3)
A
1
1
m
×
m
•
B
1
2
m
×
1
A11_{m×m}•B12_{m×1}
A11m×m•B12m×1:采用朴素算法,需要做
(
n
−
1
)
2
(n−1)^2
(n−1)2次乘法和
(
n
−
1
)
(
n
−
2
)
(n−1)(n−2)
(n−1)(n−2)次加法,因此运行时间也为
Θ
(
n
2
)
Θ(n^2)
Θ(n2)。
4)
A
1
2
m
×
1
•
B
2
2
1
×
1
A12_{m×1}•B22_{1×1}
A12m×1•B221×1:采用朴素算法,需要做
(
n
−
1
)
(n−1)
(n−1)次乘法,运行时间为
Θ
(
n
)
Θ(n)
Θ(n)。
5)
A
2
1
1
×
m
•
B
1
1
m
×
m
A21_{1×m}•B11_{m×m}
A211×m•B11m×m:采用朴素算法,需要做
(
n
−
1
)
2
(n−1)^2
(n−1)2次乘法,以及
(
n
−
1
)
(
n
−
2
)
(n−1)(n−2)
(n−1)(n−2)次加法,运行时间为
Θ
(
n
2
)
Θ(n^2)
Θ(n2)。
6)
A
2
2
1
×
1
•
B
2
1
1
×
m
A22_{1×1}•B21_{1×m}
A221×1•B211×m:采用朴素算法,需要做
(
n
−
1
)
(n−1)
(n−1)次乘法,运行时间为
Θ
(
n
)
Θ(n)
Θ(n)。
7)
A
2
1
1
×
m
•
B
1
2
m
×
1
A21_{1×m}•B12_{m×1}
A211×m•B12m×1:采用朴素算法,需要做
(
n
−
1
)
(n−1)
(n−1)次乘法,以及
(
n
−
2
)
(n−2)
(n−2)次加法,运行时间为
Θ
(
n
)
Θ(n)
Θ(n)。
8)
A
2
2
1
×
1
•
B
2
2
1
×
1
A22_{1×1}•B22_{1×1}
A221×1•B221×1:这仅仅是两个元素的相乘,只花费
Θ
(
1
)
Θ(1)
Θ(1)时间。
根据以上分析,除去
A
1
1
m
×
m
•
B
1
1
m
×
m
A11_{m×m}•B11_{m×m}
A11m×m•B11m×m之外,其他
7
7
7个矩阵乘法加起来的运行时间为
Θ
(
n
2
)
Θ(n^2)
Θ(n2)。因此,当
n
n
n为奇数时,
n
×
n
n×n
n×n矩阵乘法的运行时间为
T
(
n
)
=
7
T
(
⌊
n
/
2
⌋
)
+
Θ
(
n
2
)
T(n)=7T(⌊n/2⌋)+Θ(n^2)
T(n)=7T(⌊n/2⌋)+Θ(n2)
综合以上两种情况,无论
n
n
n为奇数还是偶数,矩阵乘法的运行时间都为
T
(
n
)
=
7
T
(
⌊
n
/
2
⌋
)
+
Θ
(
n
2
)
T(n)=7T(⌊n/2⌋)+Θ(n^2)
T(n)=7T(⌊n/2⌋)+Θ(n2)。忽略其中的⌊ ⌋符号,这与之前分析的Strassen算法的运行时间是一样的。
下面给出具备通用性的Strassen算法的伪代码。
4.2-4 如果可以用
k
k
k次乘法操作(假定乘法的交换律不成立)完成两个
3
×
3
3×3
3×3矩阵相乘,那么你可以在
o
(
n
l
g
7
)
o(n^{{\rm lg}7})
o(nlg7)时间内完成
n
×
n
n×n
n×n矩阵相乘,满足这一条件的最大
k
k
k是多少?此算法的运行时间是怎样的?
解
仍然采用Strassen算法。我们现在分析该算法运行时间的递归式,不过在这里需要以
T
(
3
)
T(3)
T(3)作为边界条件,递归式如下所示。
如果我们画出递归树,该递归树一共有
l
g
(
n
/
3
)
lg(n/3)
lg(n/3)层。叶结点对应子问题
T
(
3
)
T(3)
T(3)。由于每层的结点数是上一层的
7
7
7倍,因此第
i
i
i层包含
7
i
7^i
7i个结点。因此,叶结点一共有
7
l
g
(
n
/
3
)
=
7
l
g
n
−
l
g
3
=
7
l
g
n
/
7
l
g
3
=
n
l
g
7
/
7
l
g
3
7^{{\rm lg}(n/3)} =7^{{\rm lg}n-{\rm lg}3}=7^{{\rm lg}n}/7^{{\rm lg}3} =n^{{\rm lg}7}/7^{{\rm lg}3}
7lg(n/3)=7lgn−lg3=7lgn/7lg3=nlg7/7lg3个。因此所有叶结点的代价之和为
(
n
l
g
7
/
7
l
g
3
)
•
T
(
3
)
=
k
•
(
n
l
g
7
/
7
l
g
3
)
(n^{{\rm lg}7}/7^{{\rm lg}3})•T(3)=k•(n^{{\rm lg}7}/7^{{\rm lg}3})
(nlg7/7lg3)•T(3)=k•(nlg7/7lg3)。
如果要在
o
(
n
l
g
7
)
o(n_{{\rm lg}7})
o(nlg7)时间内完成
n
×
n
n×n
n×n矩阵相乘,那么必然有
k
•
(
n
l
g
7
/
7
l
g
3
)
<
n
l
g
7
k•(n^{{\rm lg}7}/7^{{\rm lg}3})<n^{{\rm lg}7}
k•(nlg7/7lg3)<nlg7,于是得到
k
<
7
l
g
3
≈
21.85
k<7^{{\rm lg}3}≈21.85
k<7lg3≈21.85。所以
k
k
k的最大值为
21
21
21。
4.2-5 V.Pan发现一种方法,可以用
132464
132 464
132464次乘法操作完成
68
×
68
68×68
68×68的矩阵相乘,发现另一种方法,可以用
143640
143 640
143640次乘法操作完成
70
×
70
70×70
70×70的矩阵相乘,还发现一种方法,可以用
155424
155 424
155424次乘法操作完成
72
×
72
72×72
72×72的矩阵相乘。当用于矩阵相乘的分治算法时,上述哪种方法会得到最佳的渐近运行时间?与Strassen算法相比,性能如何?
解
对于采用分治法的矩阵乘法算法来说,其运行时间都为
Θ
(
n
d
)
Θ(n^d)
Θ(nd),其中
d
d
d为一个正常数。现在分析题目所给的
3
3
3种方法,其渐近运行时间中的
d
d
d分别为多少。为方便起见,假设
3
3
3种方法的运行时间分别为
T
1
(
n
)
=
n
d
1
,
T
2
(
n
)
=
n
d
2
T_1(n)=n^{d_1},T_2(n)=n^{d_2}
T1(n)=nd1,T2(n)=nd2和
T
3
(
n
)
=
n
d
3
T_3(n)=n^{d_3}
T3(n)=nd3。
用
132464
132 464
132464次乘法操作完成
68
×
68
68×68
68×68的矩阵相乘,于是有
T
1
(
68
)
=
6
8
d
1
=
132464
T_1 (68)=68^{d_1}=132464
T1(68)=68d1=132464
得到
d
1
=
l
o
g
68
132464
≈
2.795128
d_1={\rm log}_{68}132464≈2.795128
d1=log68132464≈2.795128。
用
143640
143 640
143640次乘法操作完成
70
×
70
70×70
70×70的矩阵相乘,于是有
T
2
(
70
)
=
7
0
d
2
=
143640
T_2 (70)=70^{d_2}=143640
T2(70)=70d2=143640
得到
d
2
=
l
o
g
70
143640
≈
2.795122
d_2={\rm log}_{70}143640≈2.795122
d2=log70143640≈2.795122。
用
155424
155 424
155424次乘法操作完成
72
×
72
72×72
72×72的矩阵相乘,于是有
T
3
(
72
)
=
7
2
d
3
=
155424
T_3 (72)=72^{d_3}=155424
T3(72)=72d3=155424
得到
d
3
=
l
o
g
72
155424
≈
2.795147
d_3={\rm log}_{72}155424≈2.795147
d3=log72155424≈2.795147。
根据以上分析,第(2)种方法的渐近运行时间的指数
d
2
d_2
d2是最小的,所以第(2)种方法会得到最佳的渐近运行时间。
Strassen算法的渐近运行时间为
Θ
(
n
l
g
7
)
Θ(nlg7)
Θ(nlg7),
l
g
7
≈
2.807355
>
d
2
{\rm lg}7 ≈ 2.807355 > d_2
lg7≈2.807355>d2,因此上述第(2)种方法的性能是优于Strassen算法的。
4.2-6 用Strassen算法作为子过程来进行一个
k
n
×
n
kn×n
kn×n矩阵和一个
n
×
k
n
n×kn
n×kn矩阵相乘,最快需要花费多长时间?对两个输入矩阵规模互换的情况,回答相同的问题。
解
两个矩阵
A
k
n
×
n
A_{kn×n}
Akn×n和
B
n
×
k
n
B_{n×kn}
Bn×kn相乘,得到矩阵
C
k
n
×
k
n
C_{kn×kn}
Ckn×kn。如果要利用Strassen算法,则需要将矩阵
A
、
B
A、B
A、B和
C
C
C按下面的方式分解
矩阵
C
C
C的任意一个子矩阵
C
i
j
=
A
i
•
B
j
C_{ij} = A_i • B_j
Cij=Ai•Bj, 这是一个
n
×
n
n×n
n×n矩阵乘法,采用Strassen算法,运行时间为
Θ
(
n
l
g
7
)
Θ(n^{{\rm lg}7})
Θ(nlg7)。一共有
k
2
k^2
k2个这样的
n
×
n
n×n
n×n矩阵乘法,所以总的运行时间为
Θ
(
k
2
•
n
l
g
7
)
Θ(k^2•n^{{\rm lg}7})
Θ(k2•nlg7)。
如果将输入矩阵的规模互换,即矩阵
A
n
×
k
n
A_{n×kn}
An×kn和
B
k
n
×
n
B_{kn×n}
Bkn×n相乘,得到矩阵
C
n
×
n
Cn×n
Cn×n,那么需要将矩阵
A
A
A和
B
B
B按下面的方式分解
矩阵
C
=
A
1
•
B
1
+
A
2
•
B
2
+
…
+
A
k
•
B
k
C = A1 • B1 + A2 • B2 + … + Ak • Bk
C=A1•B1+A2•B2+…+Ak•Bk。一共有
k
k
k个
n
×
n
n×n
n×n矩阵乘法,并且还有
(
k
−
1
)
(k−1)
(k−1)个
n
×
n
n×n
n×n矩阵加法,所以总的运行时间为
Θ
(
k
•
n
l
g
7
)
Θ(k•n^{{\rm lg}7})
Θ(k•nlg7)。
4.2-7 设计算法,仅使用三次实数乘法即可完成复数
a
+
b
i
a+bi
a+bi和
c
+
d
i
c+di
c+di相乘。算法需接收
a
、
b
、
c
a、b、c
a、b、c和
d
d
d为输入,分别生成实部
a
c
−
b
d
ac−bd
ac−bd和虚部
a
d
+
b
c
ad+bc
ad+bc。
解
借鉴Strassen算法的思想,该问题可以按以下步骤解决。
(1) 计算
P
1
、
P
2
P_1、P_2
P1、P2和
P
3
P_3
P3
P
1
=
a
d
P_1 = ad
P1=ad
P
2
=
b
c
P_2 = bc
P2=bc
P
3
=
(
a
–
b
)
(
c
+
d
)
=
a
c
–
b
d
+
a
d
–
b
c
P_3 = (a – b)(c + d) = ac – bd + ad – bc
P3=(a–b)(c+d)=ac–bd+ad–bc
(2) 计算实部和虚部
实部:
P
3
–
P
1
+
P
2
=
a
c
−
b
d
P_3 – P_1 + P_2 = ac−bd
P3–P1+P2=ac−bd
虚部:
P
1
+
P
2
=
a
d
+
b
c
P_1 + P_2 = ad+bc
P1+P2=ad+bc
该算法只需要
3
3
3次乘法即可。
本节代码链接:
https://github.com/yangtzhou2012/Introduction_to_Algorithms_3rd/tree/master/Chapter04/Section_4.2