矩阵乘法_矩阵乘法Strassen算法

矩阵乘法是种极其耗时的运算。以

为例，其中

和

都是

矩阵。根据矩阵乘法的定义，

中的每个元素需要按照如下方式计算

式(4.8)包含一个

次的循环，因此计算

的时间复杂度为

。而

共有

个元素，因此总时间复杂度为

。

这是一个很可怕的数字，当

较大时，连计算机都无法承受如此的计算时间。不过，凭着科学家们的智慧，人们逐渐发现存在一些低于

的矩阵乘法算法。Strassen是这方面的先驱，他率先把矩阵乘法的时间复杂度降到了

。别小看这个细微的改进，当

非常大时，该算法将比平凡算法节约大量时间。

分治法

Strassen算法基于分治的思想，因此我们首先考虑一个简单的分治策略。

每个

的矩阵都可以分割为四个

的矩阵，我们将

、

和

重写如下。

于是

可以改写为

展开得

每个等式需要计算两次矩阵乘法和一次矩阵加法。用

表示两个

矩阵之间的乘法，则上式可以表示为如下的递归式

其中，

表示8次子矩阵乘法，子矩阵的规模为

。

表示4次矩阵加法的时间复杂度以及合并

矩阵的时间复杂度。

我们可以用上篇文章递归式求解——代入法、递归树与主定理中介绍的递归树法或主方法计算式(4.16)的解。结果是

，与平凡算法时间复杂度相同。可见，简单的分治策略并没有起到加速运算的效果。

Strassen算法

1969年，Volker Strassen发表文章提出一种渐进快于平凡算法的矩阵相乘算法，引起巨大轰动。在此之前，很少人敢设想一个算法能渐近快于平凡算法。矩阵乘法的渐近上界自此被改进了。

让我们回头观察前面使用分治策略的时候为什么无法提高速度。

因为分解后的问题包含了8次矩阵相乘和4次矩阵相加，就是这8次矩阵相乘导致了速度不能提升。于是我们想到能不能减少矩阵相乘的次数，哪怕代价是更多的矩阵相加。Strassen正是利用了这一点。

现在，我们来看一下Strassen算法的原理。

仍然把每个矩阵分割为4份，然后创建如下10个中间矩阵。

接着，计算7次矩阵乘法。

最后，根据这7个结果就可以计算出

矩阵。

是不是很神奇呢？话说我第一次看到这个算法的时候真的是惊呆了，10个

矩阵和7个

矩阵究竟是怎么凑出来的，简直不可思议。

我们可以把

矩阵和

矩阵展开，并代入最后的式子验算，会发现恰好是式(4.11~14)中的四个式子。也就是说，Strassen为了计算式(4.11~14)，绕了一大圈，用了更多的步骤，成功的把计算量变成了7个矩阵乘法和18个矩阵加法。虽然矩阵加法增加了好几倍，而矩阵乘法只减小了1个，但在数量级面前，18个加法仍然渐进快于1个乘法。这就是该算法的精妙之处。

同样地，我们可以写出Strassen算法的递推公式。

使用递归树或主方法得到解为

。

下图展示了平凡算法和Strassen算法的性能差异，

越大，Strassen算法节约的时间越多。

代码与实验

为了实际检验Strassen算法的性能，我分别实现了矩阵乘法的平凡算法和Strassen算法，完整代码见PlainMultiplier.java和StrassenMultiplier.java。实现并不复杂，不再赘述。但在实现的过程中，我发现当数据规模

比较小时，Strassen算法会明显慢于平凡算法。这很容易解释，因为Strassen算法的数据分解与合并、数据拷贝、以及大量的加法运算使得其常数因子很大。为了解决这个问题，我在递归调用的时候先判断当前

的规模，是否大于512，是则继续递归调用Strassen算法，否则调用平凡算法终止递归。其中，512是实验后得出的交叉点，在该点切换算法可以收到比较满意的效果。

在实践中还有一点需要注意，Strassen算法只适合用于

的方阵，且

必须是2的幂，因为只有这样，递归时对矩阵的划分才能正常进行。但这个要求太苛刻，实际中满足这一要求的情况很少。所以实际中，会采取一些措施，比如将一些行或列补0，使其成为方阵，再用Strassen算法计算。具体的操作细节我们就不深究了。

对Strassen的一些好奇

我猜很多人像我一样，会比较好奇Strassen当年是怎么想出这个算法的。这个疑问困扰了我很久，在网上很少能查到相关的资料，甚至Strassen当年的论文中也没有提到他的思路。

后来有人提到，Gauss曾经发现一个技巧，计算复数乘积

看似需要4次乘法，但其实换种计算方式

就只需要3次乘法，代价是加法从2次增加到5次。Strassen受到这一算法的启发，认为矩阵乘法也可以通过类似的方式实现。但Strassen的研究涉及到更高级的理论，我看过后就放弃了，目前还没法完全理解。我会把相关资料放在文章末尾，感兴趣的同学可以前往阅读。

最后，值得一提的是，研究者们经过多年努力，到目前为止，已经把矩阵乘法的时间复杂度降到了

，由斯坦福大学的学者于2014年提出。未来，不知道会不会有更快速的算法出现，但可以肯定的是，人们精益求精的脚步从来不会停息。

参考资料

Gaussian Elimination is not Optimal Strassen

Strassen's Algorithm Made (Somewhat) More Natural: A Pedagogical Remark Ann Q. Gates, Vladik Kreinovich

How did Strassen come up with his matrix multiplication method? StackExchange

Toward an Optimal Algorithm for Matrix Multiplication Sara Robinson

Multiplying matrices in O(n^2.373) time Virginia Vassilevska Williams