快速矩阵乘法的研究——中

快速矩阵乘法的研究

渐近和不等式（Asymptotic sum inequality）

有了APA之后，我们突破了Strassen 算法的阈值，但希望进一步往下降时，总是构造不出来，遇到了瓶颈。这个时候，人们发现了一条新的道路。

同时计算多个矩阵乘法

对于 $C=AB$ 这个矩阵乘法，我们已经找不到更好的构造算法，但如果同时计算 $C_1=A_1{B}_1,C_2=A_2{B}_2$，我们会不会有好方法呢？

答案是有的， 请看下面的构造：

设：
$$\varphi =\sum _{i=1}^e\sum _{j=1}^l{a}_i{b}_j{c}_{j,i}+\sum _{i=1}^{e-1}\sum _{j=1}^{l-1}{X}_{i,j}{Y}_{i,j}Z$$

这条算式同时计算了两个不相关的矩阵乘法$C=AB,Z=XY$，一个尺寸是 $[e,1,l]$，另一个是$[1,(e-1)(l-1),1]$，这种计算我们定义为 $$[e,1,l]\oplus [1,(e-1)(l-1),1]$$
不难看出，原始需要的乘法数为 $el+(e-1)(l-1)$

现在我们来构造快速算法，令：
$$X_{i,l}=0,X_{e,j}=-\sum _{i=1}^{e-1}{X}_{i,j},Y_{i,l}=-\sum _{j=1}^{l-1}{Y}_{i,j}$$

则有：
$$F^{\prime }=\sum _{i=1}^e\sum _{j=1}^l(a_i+\lambda X_{i,j})(b_j+\lambda Y_{i,j})({\lambda }^2{c}_{j,i}+Z)-\sum _{i=1}^e{a}_i\sum _{j=1}^l{b}_{j}Z={\lambda }^2\varphi +{\lambda }^{3}G(\lambda )$$

这个式子为刚才的联合矩阵乘给出了一个 $el+1$的 Border Rank，也即：
$$[e,1,l]\oplus [1,(e-1)(l-1),1]{⊴}_3(el+1)$$，若取 $e=4,l=4$，这一次展开就减少了 32% 的乘法量，十分诱人。

但这对矩阵乘法的优化有什么作用呢，请看下面的定理

Asymptotic Sum Inequality

定理：设 $\omega$为矩阵乘法的阶，假设有一个Border Rank 为 r 的算法 $\varphi$ 能够同时计算 s 个矩阵乘法$[e_1,h_1,l_1],[e_2,h_2,l_2]...[e_s,h_s,l_s]$，也即${⨁}_{i=1}^s[e_i,h_i,l_i]{⊴}_{q}r$，那么：
$$\sum _{i=1}^s(e_i{h}_i{l}_i{)}^{\omega /3}\le r$$

定理的严格证明可自行查阅文献，这里力求简单地说明思路。

首先，我们不难看出，能够以 r 个多项式计算多个矩阵乘的算法，可以把一些项设成零，推导出一些用小于等于r个多项式计算其中一个或部分矩阵乘的算法。
$$\underset{i=1}{\overset{s}{⨁}}[e_i,h_i,l_i]{⊴}_{q}r⟶[e_i,h_i,l_i]{⊴}_{q}r$$

比如上节的 $$F^{\prime }=\sum _{i=1}^e\sum _{j=1}^l(a_i+\lambda X_{i,j})(b_j+\lambda Y_{i,j})({\lambda }^2{c}_{j,i}+Z)-\sum _{i=1}^e{a}_i\sum _{j=1}^l{b}_{j}Z={\lambda }^2\varphi +{\lambda }^{3}G(\lambda )$$

我们在 $a_i,b_j,c_{j,i}$前面都乘个0，便成了$Z=XY$的算法，类似可得$C=AB$。如此进行处理，多项式的个数只会减少，不会增加。

然后，我们进行张量积，求其n次幂，方便起见只考虑两个矩阵乘算法。
先平方(n=2)，来看看发生了什么。
$$([e_1,h_1,l_1]\oplus [e_2,h_2,l_2]{)}^2=[e_1^2,h_1^2,l_1^2]\oplus [e_2^2,h_2^2,l_2^2]\oplus [e_1{e}_2,h_1{h}_2,l_1{l}_2]\oplus [e_1{e}_2,h_1{h}_2,l_1{l}_2]$$

张量积平方后，我们同时计算了四个独立的矩阵乘法，其中有2个矩阵乘是形状一样的，如果展开较多，比如达到7个，我们能将其拼起来，按上一章所述的 Strassen 算法构造一个单独的矩阵乘。

定理的证明就是在n足够大时，找一个$\mu$，并令 $P=(C_n^{\mu }{)}^{1/\omega }$，然后把这组矩阵乘拼成一个$[P{e}_1^{\mu }{e}_2^{n-\mu },P{h}_1^{\mu }{h}_2^{n-\mu },P{l}_1^{\mu }{l}_2^{n-\mu }]$的矩阵乘，这个矩阵乘可以用不大于$((q-1)N+1{)}^2{r}^n$的乘法算出来，进一步证明出最终等式。

在定理的证明过程中，我们不难发现这个矩阵乘算法构造是一个迭代过程：
我们先以当前最好的矩阵乘算法，找到一组$n,\mu ,P$，拼出一个矩阵乘算法，然后再用这个矩阵乘算法重新确定$n,\mu ,P$，再拼，直到不能下降为止。

通过这个定理，在$e=4,l=4$时，可证明出$\omega <2.5479$，大大地前进了。

Strassen 构造

前面我们介绍了如何用多个独立矩阵乘直和去构造矩阵乘法。而 Strassen 通过 Matrix To Scalar 定理，用非独立矩阵乘构造出更好的结果，将矩阵乘的阶数降到了2.48。

Matrix To Scalar

定理：一个 $[N,N,N]$的矩阵乘法，可以近似地（APA）计算 $3N^2/4$ 个独立的实数乘积。
初看这定理会感觉莫名奇妙，本身矩阵乘法包含了 $N^3$项，去计算 $3N^2/4$个实数不是很亏么。
其实是这样的，矩阵乘法$N^3$个项并不是互相独立的，以$[2,2,2]$为例，计算出来的是：
$$a_{11}{b}_{11}{c}_{11}+a_{12}{b}_{21}{c}_{11}+a_{11}{b}_{12}{c}_{21}+a_{12}{b}_{22}{c}_{21}+a_{21}{b}_{11}{c}_{12}+a_{22}{b}_{21}{c}_{12}+a_{21}{b}_{12}{c}_{22}+a_{22}{b}_{22}{c}_{22}$$
其中 $a_{11}{b}_{11}{c}_{11}$ 和 $a_{12}{b}_{21}{c}_{11}$ 共同含有 $c_{11}$，$a_{11}{b}_{11}{c}_{11}$ 和 $a_{11}{b}_{12}{c}_{21}$ 共同含有 $a_{11}$，也即不互相独立。

这个定理就是找到一个构造方法，筛出$3N^2/4$个互相独立的乘数出来。

设定$g=[3(N+1)/2]$，然后按如下方式给$[N,N,N]$矩阵乘法中的各个乘数乘以$\lambda$：
$$\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N(x_{i,j}{\lambda }^{i^2+2ij})(y_{j,k}{\lambda }^{j^2+2j(k-g)})(z_{k,i}{\lambda }^{(k-g{)}^2+2(k-g)i})=\sum _{i+j+k=g}{x}_{i,j}{y}_{j,k}{z}_{k,i}+O(\lambda )$$

如此，构造出的等式右边是那$3N^2/4$个$x,y,z$均不相同的乘数和。

Strassen 构造

用如下式子进行一个初始的APA构造：
$$\sum _{i=1}^q(x_0+\lambda x_i)(y_0+\lambda y_i)(x_i{\lambda }^{-1})+x_0{y}_0(-\sum _{i=1}^q{z}_i{\lambda }^{-1})=\sum _{i=1}^q(x_i{y}_0{z}_i+x_0{y}_i{z}_i)$$

上面的算子用 $q+1$个乘数计算了 $2q$个乘数，但右边的式子并非矩阵乘法，也不是两个互相独立的矩阵乘法和。

我们将$x_0$和 { x 1 , x 2 , . . . , x q } \{x_1, x_2, ..., x_q\} {x1​,x2​,...,xq​}各视为一个数$X_0,X_1$，类似地 { y 1 , y 2 , . . . , y q } \{y_1, y_2, ..., y_q\} {y1​,y2​,...,yq​}和 { z 1 , z 2 , . . . , z q } \{z_1, z_2, ..., z_q\} {z1​,z2​,...,zq​}也视为一个数，上面的等式右边就成了一个$[1,2,1]$的"矩阵乘法"：
$$X_1{Y}_0{Z}_1+X_0{Y}_1{Z}_1$$

我们将 x, y, z 轮换两次，构建额外的两个式子：
$$X_1{Y}_0{Z}_1+X_1{Y}_1{Z}_0,[1,1,2]$$
$$X_0{Y}_1{Z}_1+X_1{Y}_1{Z}_0,[2,1,1]$$

这两个式子和原来那个式子作张量积，我们得到一个$[2,2,2]$的"矩阵乘法"，这个矩阵乘法需要$(q+1{)}^3$个多项式：
$X_{1,1,0}{Y}_{0,0,1}{Z}_{1,1,1}+X_{1,1,1}{Y}_{0,0,1}{Z}_{1,1,0}+X_{1,1,0}{Y}_{0,1,1}{Z}_{1,0,1}+X_{1,1,1}{Y}_{0,1,1}{Z}_{1,0,0}+X_{0,1,0}{Y}_{1,0,1}{Z}_{1,1,1}+X_{0,1,1}{Y}_{1,0,1}{Z}_{1,1,0}+X_{0,1,0}{Y}_{1,1,1}{Z}_{1,0,1}+X_{0,1,1}{Y}_{1,1,1}{Z}_{1,0,0}$

这个“矩阵乘法中”，每个项都是一个$[e,h,l],e\ast h\ast l=q^3$的矩阵乘法。

然后在这个张量积的基础上，将其作N次幂，就是一个$[2^N,2^N,2^N]$的"矩阵乘法"，对应的需要$(q+3{)}^{3N}$个多项式，这时候应用 Matrix to Scalar，将其多项式独立出来，就成了$3/4\ast 2^{2N}$个满足$e\ast h\ast l=q^{3N}$的矩阵乘法和，应用 Asymptotic Sum Inequality，可得：
$$(3/4)2^{2N}{q}^{N\omega }\le (q+1{)}^{3N}$$
令N趋于无穷大，可得：
$$2^2{q}^{\omega }<=(q+1{)}^3$$
设 q = 5，便得到了$\omega <=2.4785$

Coppersmith 和 Winograd 用另外一种方法，将矩阵乘阶数降到了 2.37，这部分内容比想象中要复杂，因此放到最后一篇介绍。

参考文献

• Matrix Multiplication Via Arithmetic Progressions (DON COPPERSMITH and SIIMUEL WINOGRAD)
• On the Complexity of Matrix Multiplication (Andrew James Stothers)