矩阵乘法是种极其耗时的运算。以
为例,其中
和
都是
矩阵。根据矩阵乘法的定义,
中的每个元素需要按照如下方式计算
式(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的研究涉及到更高级的理论,我看过后就放弃了,目前还没法完全理解。我会把相关资料放在文章末尾,感兴趣的同学可以前往阅读。
最后,值得一提的是,研究者们经过多年努力,到目前为止,已经把矩阵乘法的时间复杂度降到了
,由斯坦福大学的学者于提出。未来,不知道会不会有更快速的算法出现,但可以肯定的是,人们精益求精的脚步从来不会停息。
参考资料