寻求最快速的矩阵乘实现

G-Spider

目标：
        时间最短是王道!
规则：

1.输入精度：
        double双精度型（64位），即使用matlab默认的数据精度。
2.输入值：
        A=rand(n,m),B=rand(m,p)  生成的0~1的两个随机数矩阵(矩阵行列通常在1000位左右)，依然遵循matlab的rand函数
        (线性同余伪随机数发生器，可自行实现)，计算C(n,p)=A(n,m)*B(m,p)。
3.时间测试：
        从rand开始，到矩阵乘计算结束。

附：
        不包括GUI或CUI的显示时间（此显示仅仅确认结果的正确性）。
        空间开销不计。指令集不限，程序语言不限，运行平台不限。Good luck!!

gxqcn

目标：
...
3.时间测试：
        从rand开始，到矩阵乘计算结束 ...
G-Spider 发表于 2010-10-25 23:56

我觉得应是“从rand后开始，到矩阵乘计算结束”，因为rand算法也会有所不同的，虽然影响不大。

我很少涉及浮点运算，也没针对通用矩阵进行优化的经验，所以静观高手表演。

G-Spider

2# gxqcn
郭兄这么说也有道理，对于时间的测试我在写的时候就有所顾虑，矩阵乘核涉及读取数据，计算，返回结果。读取数据的效率会严重影响计算效率，未对齐的数据可能会浪费更多读取时间，我猜想matlab在定义一个大矩阵的时候不会傻瓜似的顺序罗列...我加上rand的时间是希望在数据的预处理上动些手脚，这样可能使矩阵乘核工作的更加有效。我会征求大家意见和建议。
期待高人....

forcal

1# G-Spider
感谢楼主专门出贴讨论这个问题！
到底是高手，比我在这个帖子21＃的问题2总结的更到位：http://bbs.emath.ac.cn/thread-2709-1-1.html
一直期待国内的高手亮相，老听人说matlab的矩阵算法是最牛的，一直未全信。

做两个测试吧？也不麻烦。

1、从rand后开始，到矩阵乘计算结束。
2、从rand前开始，到矩阵乘计算结束。

forcal

我先献丑，垫个底，助助兴。

1、从rand后开始，到矩阵乘计算结束。

Forcal代码：

1. !using["math","sys"];
   
2. (:a,b,k,t0)=
   
3. oo{
   
4.   a=rand[1000,1000], b=rand[1000,1000],
   
5.   t0=clock(),
   
6.   k=a*b  //矩阵乘
   
7. },
   
8. [clock()-t0]/1000;

复制代码

耗时约2.219秒。

2、从rand前开始，到矩阵乘计算结束。

Forcal代码：

1. !using["math","sys"];
   
2. (:a,b,k,t0)=
   
3. t0=clock(),
   
4. oo{
   
5.   a=rand[1000,1000], b=rand[1000,1000],
   
6.   k=a*b  //矩阵乘
   
7. },
   
8. [clock()-t0]/1000;

复制代码

耗时约2.328秒。

我的电脑：Intel Core 2 Duo T5500 1.66G 1G内存。
只有一个核在运行（多核不会弄）。
可下载OpenFC（绿色免安装，2M）演示以上代码：http://www.forcal.net/xiazai/forcal9/openfc32w.rar

liangbch

可惜我的显卡的GPU不支持double计算，要不我写一个调用CUDA矩阵库的测试程序。

G-Spider

6# liangbch
那就先用单精度型的吧(32位)，实际上我暂时先写的float版本 ，雏形已经有了，使用SSE2系统指令集，
未对B矩阵转置，且A,B的列是4的整数倍的情形下......昨天主要完成了这个模块....

随机函数使用了c中的rand()改造版本，在种子相同的情形下，得到相同的随机数。
随机参数:

跟踪c版的rand函数：

;.......................................................

A(1024,1024)*B(1024,1024) 单精度乘法，A,B种子分别为1和2（单线程在core 2测试已确保结果正确性)，在320ms左右完成，预计加上B的转置和对非4整数倍的处理之后，在500ms内计算完成

直接把部分结果拷到文件中了

当然死算每一个乘积，几乎已到了尽头，虽然用到了矢量.我很不甘心，恐怕这种方式不可能超过matlab的矩阵乘。

刚才看到一篇斯特拉森矩阵乘法的文章我颇有顾虑:
数据从rand a,b到转成0~1之间的浮点，整个过程也只用了15ms左右
接着是计算a*b
乘法的代价过大，只有提供一种高效的算法来减少乘法的运算次数才能有质的飞跃....酝酿中.....

litaoye

矩阵乘法我只知道那个n^log(2,7)的分治方法，难道有类似fft的变换么？

forcal

乘法的代价过大，只有提供一种高效的算法来减少乘法的运算次数才能有质的飞跃....酝酿中..... G-Spider 发表于 2010-10-29 14:24

我在5＃的矩阵乘算法仅使用C/C++实现，如果你也仅使用C/C++实现，能达到或超过我这个速度，那我就帮不上忙了。

否则，使用我这个雕虫小技或许会使速度有更大的提高。

G-Spider

9# forcal
8# litaoye
看来要使用一些高级的算法才能降低时间复杂度.
看一下乘法的计算次数：
普通：0(n^3)......1024^3=1 073 741 824次乘运算
Strassen分治：0(n^2.81)............1024^2.81=  287 701 988次算运算
Don Coppersmith and Shmuel Wmograd改进：0(n^2.367).........1024^2.367= 13 346 887次算运算
........
当然在实现上的难度会越来越大............
最新的矩阵乘算法动态，可参考：
Fast and stable matrix multiplication