快速矩阵乘法一
新开一个系列，学习一下平时视作黑箱略过的基础技术。如果经常用一些知名结论却不懂证明思路，恐怕也不够格自称理论计算机科学家。由于是现学现卖，内容会比较肤浅。假定读者的水平和我差不多，即学了微积分线性代数的小学二年级学生。内容以翻译为主，原文在文末“教材”中列出。中文知识不足处会使用非标准直译术语。
一、简单结论
今有二矩阵 \(x\in F^{n\times k}, y\in F^{k\times m}\), 求\(z=xy\in F^{n\times m}\), \(z_{ij}=\sum_k x_{ik} y_{kj}\). 元素假定是在域\(F\)里. 暴力算是\(O(n^3)\). 直觉上每一项都挺必要的, 但是也有冗余信息.
众所周知\(2\times 2\)分块就可以做到更好. 可以用七次乘法代替八次, 虽然不太直观:
\[\begin{align*} p_1 &= (x_{11}+x_{22})(y_{11}+y_{22})\\
p_2 &= (x_{21}+x_{22})y_{11}\\
p_3 &= x_{11}(y_{12}-y_{22})\\
p_4 &= x_{22}(y_{21}-y_{11})\\
p_5 &= (x_{11}+x_{12})y_{22}\\
p_6 &= (x_{21}-x_{11})(y_{11}+y_{12})\\
p_7 &= (x_{12}-x_{22})(y_{21}+y_{22})\\
z_{11} = x_{11}y_{11}+x_{12}y_{21} &= p_1+p_4-p_5+p_7\\
z_{12} = x_{11}y_{12}+x_{12}y_{22} &= p_3+p_5\\
z_{21} = x_{21}y_{11}+x_{22}y_{21} &= p_2+p_4\\
z_{22} = x_{21}y_{12}+x_{22}y_{22} &= p_1+p_3-p_2+p_6
\end{align*}\]
我们将\(x, y\)分成四个边长减半的小矩阵, 然后递归计算七个小矩阵乘, 再进行加减即可得到\(z\)的四个部分. 时间的递归式为\(T(n)=7T(\frac{n}{2})+O(n^2)\), 解得\(T(n)=O(n^{\log 7})=O(n^{2.807})\).
该方法[地]记作简单分块乘. 从中也可以看出, 递归时只有乘法的次数影响时间复杂度, 加减法可以忽略不计.
生活中常用的常数\(\omega\)意为两个\(n\times n\)矩阵相乘需要时间\(O(n^\omega)\). 二〇二〇年现在其值为2.3728639[玄]. 本系列试图让自己相信这个值是对的.
二、张量模型
简单分块乘算式虽繁, 细究可见从\(x, y\)到\(p\)是双线性, 从\(p\)到\(z\)是线性. 张量长于刻画多线性, 有助于深刻地分析此类计算方式.
我没学过张量, 好在学量子信息时接触过. 稍作回顾. 一个量子比特可以看作0态与1态的叠加, 用二维复向量表示 (要求模长为1), 0态和1态分别记作\(|0\rangle=(1,0)\)和\(|1\rangle=(0,1)\). 假如取两个独立的0态量子比特, 得到的状态就是\(|0\rangle\otimes|0\rangle=|00\rangle\). 这里的\(\otimes\)就是张量积, 将两个向量积成二阶张量 (2×2维). 两个量子比特的状态当然也可以是\(|01\rangle, |10\rangle, |11\rangle\). 现如我们取两个叠加态的独立量子比特放到一起，会得到\((|0\rangle+|1\rangle)\otimes(|0\rangle+|1\rangle) = |00\rangle+|01\rangle+|10\rangle+|11\rangle\) (省略归一化). 这个状态等于前面四个状态的叠加. 此例中我们看到张量积对每一部分都是线性的. 这个状态只需一次乘法, 秩是1, 两量子比特独立. 对于纠缠态秩就大于1. 纠缠态不能拆成两个独立量子比特的积，如\(|00\rangle+|11\rangle\).
下面我们把\(x, y\)拍扁, 随即需要记录拍扁后哪些项相乘. 拍扁后的下标代换为\(\kappa=(j, i), \mu=(i, k), \nu=(k, j)\). 记录乘法规则的三阶张量\(t\)定义如下. 若\(\kappa, \mu, \nu\)表示的\(i, j, k\)相容, 则\(t_{\kappa\mu\nu}\)为1, 否则为0. 从而\[z_\kappa=\sum_\mu\sum_\nu t_{\kappa\mu\nu}x_\mu y_\nu.\]
再用拍扁的下标重写简单分块乘.
\[p_\lambda = (\sum_\mu u_{\lambda\mu} x_\mu)(\sum_\nu v_{\lambda\nu} y_\nu)\]
\[z_\kappa = \sum_\lambda w_{\lambda\kappa} p_\lambda = \sum_{\lambda\mu\nu} w_{\lambda\kappa} u_{\lambda\mu} v_{\lambda\nu} x_\nu y_\nu\]
与定义式对照即得\(t_{\kappa\mu\nu}=\sum_\lambda w_{\lambda\kappa} u_{\lambda\mu} v_{\lambda\nu}\). 把\(w_\lambda, u_\lambda, v_\lambda\)看作向量, 张量\(t=\sum_\lambda w_\lambda \otimes u_\lambda \otimes v_\lambda\). 张量\(t\)的秩\(R(t)\)是\(t\)分解成独立张量积的最少个数, 不超过上式中\(\lambda\)的个数.
由上面推导可得, 已知秩为\(R\)的张量分解可以写出\(R\)次乘法的计算方式, 已知计算方式也可以写出张量分解. 其实此言稍有问题, 因为计算的时候并非只能将\(x\)的线性组合与\(y\)的线性组合相乘, 还可以算\(x, y\)的线性组合再乘. 不过这样的计算方式也可以写成两个张量分解的和, 因为\(z\)中只有\(x_\mu y_\nu, y_\mu x_\nu\)项没有\(x_\mu x_\nu, y_\mu y_\nu\)项. 总之\(t\)的秩与所需乘法次数差不多，从而也反映了计算时间。
现在我们只需研究矩阵乘张量\(t\). \(t\)只依赖于维度, 记作\(\langle n,k,m\rangle\), 秩记作\(R(n,k,m)\). 第一节实际上说的就是\(R(2,2,2)\le 7\). \(t\)的定义对于\(\kappa, \mu, \nu\)是对称的, 从而对\(i, j, k\)也是对称的. 一个直接的推论是, 如果我们调换\(n,k,m\)的次序, \(t\)只要换一下下标位置就能用同样的方法分解. 于是\(R(n,k,m)\)与\(n,k,m\)的顺序无关.
\(k\times m\times n\)维张量\(t\)与\(k'\times m'\times n'\)维张量\(t'\)的张量积\(s=t\otimes t'\)是\((kk')\times (mm')\times (nn')\)维张量， \(s_{\kappa\kappa', \mu\mu', \nu\nu'}=t_{\kappa\mu\nu} t'_{\kappa'\mu'\nu'}\). 张量积是线性的, 将\(t\)的分解与\(t'\)的分解相乘可得\(s\)的分解, 展开后共\(R(t)R(t')\)项. 故\(R(t\otimes t')\le R(t)R(t')\).
张量积可以用来构造更大的矩阵乘张量. 对于\(<n,k,m>\)与\(<n',k',m'>\)的张量积, \((\kappa\kappa', \mu\mu', \nu\nu')\)相容当且仅当\((\kappa, \mu, \nu)\)与\((\kappa', \mu', \nu')\)都相容, 故\(<n,k,m>\otimes<n',k',m'>\)与\(<nn',kk',mm'>\)是同一个张量. 将低维度的秩相乘, 即知高维度的秩的上界. 特别地, 当\(R(k,m,n)\le r\), \(R(kmn, kmn, kmn)\le\)\(R(k,m,n)R(m,n,k)R(n,k,m)\le r^3\). 于是\(\omega\le 3\log_{kmn} r\), 因为递归式将会变成\(T(N)=r^3T(N/kmn)+O(n^2)\).
三、边界秩
线性代数中, 我们面对高秩矩阵时常用奇异值分解来得到近似的低秩矩阵. 对张量也想要以容忍近似为代价降低秩.
将秩分解中的\(w, u, v\)用多项式环\(F[\epsilon]\)上的向量代替. 每个数现在都看作多项式. 若存在\(r\)组多项式向量的张量积之和是\(\epsilon^h t+O(\epsilon^{h+1})\), 记\(R_h(t)=r\). 这里的\(O(\epsilon^{h+1})\)表示只含有次数高于\(h\)的项, 也就是假装\(\epsilon\)是个无穷小量. 令\(\underline{R}(t)\)为最小的\(R_h(t)\), 称为边界秩.
边界秩可以作为秩的上界. 存在依赖于\(h\)的常数\(c_h\le {h+2 \choose 2}\), \(R(t)\le c_hR_h(t)\). 这是因为, 一旦得到\(\epsilon^h t+O(\epsilon^{h+1})\)的秩分解, 我们只考虑乘积中\(\epsilon^h\)的项就能求\(t\). 而三个多项式相乘时只有\({h+2 \choose 2}\)组系数的积对\(h\)次项有贡献. 总之\(r\)组多项式向量可以转化为\({h+2 \choose 2}r\)组向量.
现在看看靠近似降低秩的效果. 考虑\(2\times 2\)矩阵乘\(2\times 3\)矩阵. 易知\(R(2,2,3)\le 11\), 七个与第一节中\(\langle 2,2,2\rangle\)的相同，四个用来暴力计算剩下两项. 下面证明\(R_1(2,2,3)\le 10\).
\[\begin{align*} p_1 &= (x_{12}+\epsilon x_{22})y_{21}\\
p_2 &= x_{11}(y_{11}+\epsilon y_{12})\\
p_3 &= x_{12}(y_{11}+y_{21}+\epsilon y_{22})\\
p_4 &= (x_{11}+x_{12}+\epsilon x_{21})y_{11}\\
p_5 &= (x_{12}+\epsilon x_{21})(y_{11}+\epsilon y_{22})\\
\epsilon z_{11} &= \epsilon p_1+\epsilon p_2+O(\epsilon^2)\\
\epsilon z_{12} &= p_2 - p_4+ p_5+O(\epsilon^2)\\
\epsilon z_{21} &= p_1 - p_3+ p_5+O(\epsilon^2)
\end{align*}\]
\(z_{33}, z_{32}, z_{22}\)的计算与\(z_{11}, z_{12}, z_{21}\)对称, 可以通过类似定义\(p_6\)至\(p_{10}\)表出.
\(R_h\)满足\(R_{h+h'}(t\otimes t')\le R_h(t)R_{h'}(t')\), 这只需将两个分解式相乘即可看出. 于是我们有\(R_3(12,12,12)\le R_1(2,2,3)^3=1000\). 不仅如此, \(R_{3s}(12^s)\le 1000^s\). 故\(R(12^s) \le c_{3s} R_{3s}(12^s) \le 1000^s c_{3s}\).
\[\omega\le \log_{12^s} (1000^s c_{3s}) \le \log_{12} 1000 + \frac{\ln {3s+2 \choose 2}}{s\ln 12} \to 3\log_{12} 10 (s\to\infty)\]
也就是说\(\omega\le 2.780\). 于是我们超越了简单分块乘, 虽然只有一点点.
至于更强的结果，请听下回分解。
教材
[天] Markus Blaeser. Fast Matrix Multiplication. Theory of Computing, 2013.
引用
[地] Volker Strassen. Gaussian Elimination is not Optimal. Numerische Mathematik, 1969.
