STRASSEN矩阵乘法

Strassen矩阵乘法

 

\(
(1) ~~
\begin{equation}\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix}
\begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix}
=
\begin{bmatrix} a_{11}b_{11}+a_{12}b_{21} & a_{11}b_{12}+a_{12}b_{22} \\ a_{21}b_{11}+a_{22}b_{21} & a_{21}b_{12}+a_{22}b_{22} \end{bmatrix}
\end{equation}
\)

 

\(
(2) ~~
\begin{equation}
\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}
\begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}
=
\begin{bmatrix} A_{11}B_{11}+A_{12}B_{21} & A_{11}B_{12}+A_{12}B_{22} \\ A_{21}B_{11}+A_{22}B_{21} & A_{21}B_{12}+A_{22}B_{22} \end{bmatrix}
\end{equation}
\)

 

\(
(3) ~~
\begin{equation}
\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}
\begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}
=
\begin{bmatrix} C_{11} & C_{12} \\ C_{21} & C_{22} \end{bmatrix}
\end{equation}
\)

 

\(
\begin{equation*}
P_1 := A_{11}(B_{12}-B_{22}) \\
P_2 := (A_{11}+A_{12})B_{22} \\
P_3 := (A_{21}+A_{22})B_{11} \\
P_4 := A_{22}(B_{21}-B_{11}) \\
P_5 := (A_{11}+A_{22})(B_{11}+B_{22}) \\
P_6 := (A_{12}-A_{22})(B_{21}+B_{22}) \\
P_7 := (A_{11}-A_{21})(B_{11}+B_{12})
\end{equation*}
\)

 

\(
\begin{equation*}
C_{11} := P_5+P_4-P_2+P_6 \\
C_{12} := P_1+P_2 \\
C_{21} := P_3+P_4 \\
C_{22} := P_5+P_1-P_3-P_7
\end{equation*}
\)

 

为更好的讨论算法本身，本文矩阵特指行列数相等(方阵)且都为2的幂的矩阵。矩阵乘法的定义是\((AB)_{ij} = \sum_{k=0}^n A_{ik}B_{kj}\)。使用定义直接求\(AB\)需要嵌套3层循环，如果\(AB\)是\(m\times m\)的，时间复杂度是\(O(m^3)\)。如何高效的做矩阵乘法是个重要的问题，计算机各领域都需要用到。目前计算机学者们虽未找到该算法足够“紧”的“理论下界”，但一直都在发明更优的具体算法，目前最快“斯坦福方法”代价为\(O(m^{2.373})\)，第一个快于\(O(m^3)\)的算法由德国数学家Strassen于1969年提出。

结合上述(1)和(2)找一种分治算法，(2)把方阵均分为四象限，把每个象限作为“子矩阵”。根据矩阵乘法定义，矩阵乘法可以正确转化成(2)等号右边的递归形式。分析这种算法的时间复杂度，其包含8个“子矩阵乘法”和4次加法，所有加法和其他操作的总代价等于遍历矩阵的\(O(m^2)\)。所以时间复杂度满足\(T(m)=8T(\frac{m}{2}) + O(m^2)\)，根据主定理\(T(m)=O(m^3)\)。

上述简单的分治算法并没有更优的时间复杂度，然后来看同样基于划分“子矩阵”思路的Strassen乘法，其划分方式如(3)所示，包含7次乘法和18次加减法，而(2)包含8次乘法和仅仅4次的加减法，所以其优化思想是“为降低乘法次数，可以不计代价的增多加减法次数”，因为根据主定理，只要加减法是常数次都不会影响整体的时间复杂度，但乘法次数影响时间复杂度。而之所以能用加减法次数换乘法次数，是因为矩阵乘法满足分配律，巧妙的运用分配律能降低乘法次数。由于只有7次乘法，Strassen乘法的时间复杂度满足\(T(m)=7T(\frac{m}{2}) + O(m^2)\)，根据主定理\(T(m)=O(m^{log_2{7}})=O(m^{2.807})\)。

下面是具体实现，这里偷懒用一维数组表示这种行列数相等的方阵，并且用到了Python面向对象的运算符重载。