还剩9页未读,继续阅读
本资源只提供10页预览,全部文档请下载后查看!喜欢就下载吧,查找使用更方便
文本内容:
矩阵相乘的快速算法算法介绍矩阵相乘在进行3D变换的时候是经常用到的在应用中常用矩阵相乘的定义算法对其进行计算这个算法用到了大量的循环和相乘运算,这使得算法效率不高而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度,所以对矩阵相乘算法进行一些改进是必要的这里要介绍的矩阵算法称为斯特拉森方法,它是由v.斯特拉森在1969年提出的一个方法我们先讨论二阶矩阵的计算方法对于二阶矩阵a11a12b11b12A=a21a22B=b21b22先计算下面7个量1x1=a11+a22*b11+b22;x2=a21+a22*b11;x3=a11*b12-b22;x4=a22*b21-b11;x5=a11+a12*b22;x6=a21-a11*b11+b12;x7=a12-a22*b21+b22;再设C=AB根据矩阵相乘的规则,C的各元素为2c11=a11*b11+a12*b21c12=a11*b12+a12*b22c21=a21*b11+a22*b21c22=a21*b12+a22*b22比较12,C的各元素可以表示为3c11=x1+x4-x5+x7c12=x3+x5c21=x2+x4c22=x1+x3-x2+x6根据以上的方法,我们就可以计算4阶矩阵了,先将4阶矩阵A和B划分成四块2阶矩阵,分别利用公式计算它们的乘积,再使用13来计算出最后结果__11__12mb11mb12A4=__21__22B4=mb21mb22其中a11a12a13a14b11b12b13b14__11=a21a22__12=a23a24mb11=b21b22mb12=b23b24a31a32a33a34b31b32b33b34__21=a41a42__22=a43a44mb21=b41b42mb22=b43b44实现//计算2X2矩阵voidMultiply2X2floatfOut_11floatfOut_12floatfOut_21floatfOut_22floatf1_11floatf1_12floatf1_21floatf1_22floatf2_11floatf2_12floatf2_21floatf2_22{constfloatx1f1_11+f1_22*f2_11+f2_22;constfloatx2f1_21+f1_22*f2_11;constfloatx3f1_11*f2_12-f2_22;constfloatx4f1_22*f2_21-f2_11;constfloatx5f1_11+f1_12*f2_22;constfloatx6f1_21-f1_11*f2_11+f2_12;constfloatx7f1_12-f1_22*f2_21+f2_22;fOut_11=x1+x4-x5+x7;fOut_12=x3+x5;fOut_21=x2+x4;fOut_22=x1-x2+x3+x6;}//计算4X4矩阵voidMultiplyCLAY__TRIXmOutconstCLAY__TRIXm1constCLAY__TRIXm2{floatfTmp
[7]
[4];//__11+__22*mb11+mb22Multiply2X2fTmp
[0]
[0]fTmp
[0]
[1]fTmp
[0]
[2]fTmp
[0]
[3]m
1._11+m
1._33m
1._12+m
1._34m
1._21+m
1._43m
1._22+m
1._44m
2._11+m
2._33m
2._12+m
2._34m
2._21+m
2._43m
2._22+m
2._44;//__21+__22*mb11Multiply2X2fTmp
[1]
[0]fTmp
[1]
[1]fTmp
[1]
[2]fTmp
[1]
[3]m
1._31+m
1._33m
1._32+m
1._34m
1._41+m
1._43m
1._42+m
1._44m
2._11m
2._12m
2._21m
2._22;//__11*mb12-mb22Multiply2X2fTmp
[2]
[0]fTmp
[2]
[1]fTmp
[2]
[2]fTmp
[2]
[3]m
1._11m
1._12m
1._21m
1._22m
2._13-m
2._33m
2._14-m
2._34m
2._23-m
2._43m
2._24-m
2._44;//__22*mb21-mb11Multiply2X2fTmp
[3]
[0]fTmp
[3]
[1]fTmp
[3]
[2]fTmp
[3]
[3]m
1._33m
1._34m
1._43m
1._44m
2._31-m
2._11m
2._32-m
2._12m
2._41-m
2._21m
2._42-m
2._22;//__11+__12*mb22Multiply2X2fTmp
[4]
[0]fTmp
[4]
[1]fTmp
[4]
[2]fTmp
[4]
[3]m
1._11+m
1._13m
1._12+m
1._14m
1._21+m
1._23m
1._22+m
1._24m
2._33m
2._34m
2._43m
2._44;//__21-__11*mb11+mb12Multiply2X2fTmp
[5]
[0]fTmp
[5]
[1]fTmp
[5]
[2]fTmp
[5]
[3]m
1._31-m
1._11m
1._32-m
1._12m
1._41-m
1._21m
1._42-m
1._22m
2._11+m
2._13m
2._12+m
2._14m
2._21+m
2._23m
2._22+m
2._24;//__12-__22*mb21+mb22Multiply2X2fTmp
[6]
[0]fTmp
[6]
[1]fTmp
[6]
[2]fTmp
[6]
[3]m
1._13-m
1._33m
1._14-m
1._34m
1._23-m
1._43m
1._24-m
1._44m
2._31+m
2._33m
2._32+m
2._34m
2._41+m
2._43m
2._42+m
2._44;//第一块mOut._11=fTmp
[0]
[0]+fTmp
[3]
[0]-fTmp
[4]
[0]+fTmp
[6]
[0];mOut._12=fTmp
[0]
[1]+fTmp
[3]
[1]-fTmp
[4]
[1]+fTmp
[6]
[1];mOut._21=fTmp
[0]
[2]+fTmp
[3]
[2]-fTmp
[4]
[2]+fTmp
[6]
[2];mOut._22=fTmp
[0]
[3]+fTmp
[3]
[3]-fTmp
[4]
[3]+fTmp
[6]
[3];//第二块mOut._13=fTmp
[2]
[0]+fTmp
[4]
[0];mOut._14=fTmp
[2]
[1]+fTmp
[4]
[1];mOut._23=fTmp
[2]
[2]+fTmp
[4]
[2];mOut._24=fTmp
[2]
[3]+fTmp
[4]
[3];//第三块mOut._31=fTmp
[1]
[0]+fTmp
[3]
[0];mOut._32=fTmp
[1]
[1]+fTmp
[3]
[1];mOut._41=fTmp
[1]
[2]+fTmp
[3]
[2];mOut._42=fTmp
[1]
[3]+fTmp
[3]
[3];//第四块mOut._33=fTmp
[0]
[0]-fTmp
[1]
[0]+fTmp
[2]
[0]+fTmp
[5]
[0];mOut._34=fTmp
[0]
[1]-fTmp
[1]
[1]+fTmp
[2]
[1]+fTmp
[5]
[1];mOut._43=fTmp
[0]
[2]-fTmp
[1]
[2]+fTmp
[2]
[2]+fTmp
[5]
[2];mOut._44=fTmp
[0]
[3]-fTmp
[1]
[3]+fTmp
[2]
[3]+fTmp
[5]
[3];}比较在标准的定义算法中我们需要进行n*n*n次乘法运算,新算法中我们需要进行7log2n次乘法,对于最常用的4阶矩阵 原算法新算法加法次数487248次加法,24次减法乘法次数6449需要额外空间16*sizeoffloat28*sizeoffloat新算法要比原算法多了24次减法运算,少了15次乘法但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高
一、两个矩阵相乘的经典算法: 若设Q=M*N其中,M是m1*n1矩阵,N是m2*n2矩阵当n1=m2时有 fori=1;i forj=1;j=n2;++j{ Q[i][j]=0; fork=1;k=n1;++k Q[i][j]+=M[i][k]*N[k][j]; } 此算法的时间复杂度是Om1*n1*n2
二、斯特拉森算法斯特拉森方法,是由v.斯特拉森在1969年提出的一个方法我们先讨论二阶矩阵的计算方法对于二阶矩阵a11a12b11b12A=a21a22B=b21b22先计算下面7个量1x1=a11+a22*b11+b22;x2=a21+a22*b11;x3=a11*b12-b22;x4=a22*b21-b11;x5=a11+a12*b22;x6=a21-a11*b11+b12;x7=a12-a22*b21+b22;再设C=AB根据矩阵相乘的规则,C的各元素为2c11=a11*b11+a12*b21c12=a11*b12+a12*b22c21=a21*b11+a22*b21c22=a21*b12+a22*b22比较12,C的各元素可以表示为3c11=x1+x4-x5+x7c12=x3+x5c21=x2+x4c22=x1+x3-x2+x6根据以上的方法,我们就可以计算4阶矩阵了,先将4阶矩阵A和B划分成四块2阶矩阵,分别利用公式计算它们的乘积,再使用13来计算出最后结果__11__12mb11mb12A4=__21__22B4=mb21mb22其中a11a12a13a14b11b12b13b14__11=a21a22__12=a23a24mb11=b21b22mb12=b23b24a31a32a33a34b31b32b33b34__21=a41a42__22=a43a44mb21=b41b42mb22=b43b44实现//计算2X2矩阵voidMultiply2X2floatfOut_11floatfOut_12floatfOut_21floatfOut_22floatf1_11floatf1_12floatf1_21floatf1_22floatf2_11floatf2_12floatf2_21floatf2_22{constfloatx1f1_11+f1_22*f2_11+f2_22;constfloatx2f1_21+f1_22*f2_11;constfloatx3f1_11*f2_12-f2_22;constfloatx4f1_22*f2_21-f2_11;constfloatx5f1_11+f1_12*f2_22;constfloatx6f1_21-f1_11*f2_11+f2_12;constfloatx7f1_12-f1_22*f2_21+f2_22;fOut_11=x1+x4-x5+x7;fOut_12=x3+x5;fOut_21=x2+x4;fOut_22=x1-x2+x3+x6;}//计算4X4矩阵voidMultiplyCLAY__TRIXmOutconstCLAY__TRIXm1constCLAY__TRIXm2{floatfTmp
[7]
[4];//__11+__22*mb11+mb22Multiply2X2fTmp
[0]
[0]fTmp
[0]
[1]fTmp
[0]
[2]fTmp
[0]
[3]m
1._11+m
1._33m
1._12+m
1._34m
1._21+m
1._43m
1._22+m
1._44m
2._11+m
2._33m
2._12+m
2._34m
2._21+m
2._43m
2._22+m
2._44;//__21+__22*mb11Multiply2X2fTmp
[1]
[0]fTmp
[1]
[1]fTmp
[1]
[2]fTmp
[1]
[3]m
1._31+m
1._33m
1._32+m
1._34m
1._41+m
1._43m
1._42+m
1._44m
2._11m
2._12m
2._21m
2._22;//__11*mb12-mb22Multiply2X2fTmp
[2]
[0]fTmp
[2]
[1]fTmp
[2]
[2]fTmp
[2]
[3]m
1._11m
1._12m
1._21m
1._22m
2._13-m
2._33m
2._14-m
2._34m
2._23-m
2._43m
2._24-m
2._44;//__22*mb21-mb11Multiply2X2fTmp
[3]
[0]fTmp
[3]
[1]fTmp
[3]
[2]fTmp
[3]
[3]m
1._33m
1._34m
1._43m
1._44m
2._31-m
2._11m
2._32-m
2._12m
2._41-m
2._21m
2._42-m
2._22;//__11+__12*mb22Multiply2X2fTmp
[4]
[0]fTmp
[4]
[1]fTmp
[4]
[2]fTmp
[4]
[3]m
1._11+m
1._13m
1._12+m
1._14m
1._21+m
1._23m
1._22+m
1._24m
2._33m
2._34m
2._43m
2._44;//__21-__11*mb11+mb12Multiply2X2fTmp
[5]
[0]fTmp
[5]
[1]fTmp
[5]
[2]fTmp
[5]
[3]m
1._31-m
1._11m
1._32-m
1._12m
1._41-m
1._21m
1._42-m
1._22m
2._11+m
2._13m
2._12+m
2._14m
2._21+m
2._23m
2._22+m
2._24;//__12-__22*mb21+mb22Multiply2X2fTmp
[6]
[0]fTmp
[6]
[1]fTmp
[6]
[2]fTmp
[6]
[3]m
1._13-m
1._33m
1._14-m
1._34m
1._23-m
1._43m
1._24-m
1._44m
2._31+m
2._33m
2._32+m
2._34m
2._41+m
2._43m
2._42+m
2._44;//第一块mOut._11=fTmp
[0]
[0]+fTmp
[3]
[0]-fTmp
[4]
[0]+fTmp
[6]
[0];mOut._12=fTmp
[0]
[1]+fTmp
[3]
[1]-fTmp
[4]
[1]+fTmp
[6]
[1];mOut._21=fTmp
[0]
[2]+fTmp
[3]
[2]-fTmp
[4]
[2]+fTmp
[6]
[2];mOut._22=fTmp
[0]
[3]+fTmp
[3]
[3]-fTmp
[4]
[3]+fTmp
[6]
[3];//第二块mOut._13=fTmp
[2]
[0]+fTmp
[4]
[0];mOut._14=fTmp
[2]
[1]+fTmp
[4]
[1];mOut._23=fTmp
[2]
[2]+fTmp
[4]
[2];mOut._24=fTmp
[2]
[3]+fTmp
[4]
[3];//第三块mOut._31=fTmp
[1]
[0]+fTmp
[3]
[0];mOut._32=fTmp
[1]
[1]+fTmp
[3]
[1];mOut._41=fTmp
[1]
[2]+fTmp
[3]
[2];mOut._42=fTmp
[1]
[3]+fTmp
[3]
[3];//第四块mOut._33=fTmp
[0]
[0]-fTmp
[1]
[0]+fTmp
[2]
[0]+fTmp
[5]
[0];mOut._34=fTmp
[0]
[1]-fTmp
[1]
[1]+fTmp
[2]
[1]+fTmp
[5]
[1];mOut._43=fTmp
[0]
[2]-fTmp
[1]
[2]+fTmp
[2]
[2]+fTmp
[5]
[2];mOut._44=fTmp
[0]
[3]-fTmp
[1]
[3]+fTmp
[2]
[3]+fTmp
[5]
[3];}比较在标准的定义算法中我们需要进行n*n*n次乘法运算,新算法中我们需要进行7log2n次乘法,对于最常用的4阶矩阵 原算法新算法加法次数487248次加法,24次减法乘法次数6449需要额外空间16*sizeoffloat28*sizeoffloat新算法要比原算法多了24次减法运算,少了15次乘法但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高
三、Strassen矩阵乘法矩阵乘法是线性代数中最常见的运算之一,它在数值计算中有广泛的应用若A和B是2个n×n的矩阵,则它们的乘积C=AB同样是一个n×n的矩阵A和B的乘积矩阵C中的元素C[ij]定义为: 若依此定义来计算A和B的乘积矩阵C,则每计算C的一个元素C[ij],需要做n个乘法和n-1次加法因此,求出矩阵C的n2个元素所需的计算时间为0n360年代末,Strassen采用了类似于在大整数乘法http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/problems/problem_set/bignumber_mul/problem.htm中用过的分治技术http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/technique/divide_and_conquer/index.htm,将计算2个n阶矩阵乘积所需的计算时间改进到Onlog7=On
2.18首先,我们还是需要假设n是2的幂将矩阵A,B和C中每一矩阵都分块成为4个大小相等的子矩阵,每个子矩阵都是n/2×n/2的方阵由此可将方程C=AB重写为: 1由此可得: C11=A11B11+A12B21 2C12=A11B12+A12B22 3C21=A21B11+A22B21 4C22=A21B12+A22B22 5如果n=2,则2个2阶方阵的乘积可以直接用2-3式计算出来,共需8次乘法和4次加法当子矩阵的阶大于2时,为求2个子矩阵的积,可以继续将子矩阵分块,直到子矩阵的阶降为2这样,就产生了一个分治http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/technique/divide_and_conquer/index.htm降阶的递归http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/technique/recursion/index.htm算法依此算法,计算2个n阶方阵的乘积转化为计算8个n/2阶方阵的乘积和4个n/2阶方阵的加法2个n/2×n/2矩阵的加法显然可以在c*n2/4时间内完成,这里c是一个常数因此,上述分治法的计算时间耗费Tn应该满足这个递归方程的解http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/complexity/chapter
6.htm仍然是Tn=On3因此,该方法并不比用原始定义直接计算更有效究其原因,乃是由于式2-5并没有减少矩阵的乘法次数而矩阵乘法耗费的时间要比矩阵加减法耗费的时间多得多要想改进矩阵乘法的计算时间复杂性,必须减少子矩阵乘法运算的次数按照上述分治法的思想http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/technique/divide_and_conquer/index.htm可以看出,要想减少乘法运算次数,关键在于计算2个2阶方阵的乘积时,能否用少于8次的乘法运算Strassen提出了一种新的算法来计算2个2阶方阵的乘积他的算法只用了7次乘法运算,但增加了加、减法的运算次数这7次乘法是: M1=A11B12-B22M2=A11+A12B22M3=A21+A22B11M4=A22B21-B11M5=A11+A22B11+B22M6=A12-A22B21+B22M7=A11-A21B11+B12做了这7次乘法后,再做若干次加、减法就可以得到: C11=M5+M4-M2+M6C12=M1+M2C21=M3+M4C22=M5+M1-M3-M7以上计算的正确性很容易验证例如: C22=M5+M1-M3-M7 =A11+A22B11+B22+A11B12-B22-A21+A22B11-A11-A21B11+B12 =A11B11+A11B22+A22B11+A22B22+A11B12 -A11B22-A21B11-A22B11-A11B11-A11B12+A21B11+A21B12 =A21B12+A22B22 由2式便知其正确性至此,我们可以得到完整的Strassen算法如下pro__dureSTRASSENnABC;begin ifn=2then__TRIX-MULTIPLYA,B,C elsebegin 将矩阵A和B依1式http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/commonalg/misc/strassen/strassen.htm\leqn1分块; STRASSENn/2A11B12-B22M1; STRASSENn/2A11+A12B22M2; STRASSENn/2A21+A22B11M3; STRASSENn/2A22B21-B11M4; STRASSENn/2A11+A22B11+B22M5; STRASSENn/2A12-A22B21+B22M6; STRASSENn/2A11-A21B11+B12M7; ; end;end;其中__TRIX-MULTIPLYA,B,C是按通常的矩阵乘法计算C=AB的子算法Strassen矩阵乘积分治算法中,用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算由此可知,该算法的所需的计算时间Tn满足如下的递归方程:按照解递归方程http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/complexity/chapter
6.htm的套用公式法http://___.comp.nus.edu.sg/~xujia/mirror/algorithm.myri__.com/algorithm/complexity/chapter6_
3.htm,其解为Tn=Onlog7≈On
2.81由此可见,Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进有人曾列举了计算2个2阶矩阵乘法的36种不同方法但所有的方法都要做7次乘法除非能找到一种计算2阶方阵乘积的算法,使乘法的计算次数少于7次,按上述思路才有可能进一步改进矩阵乘积的计算时间的上界但是Hopcroft和Kerr197l已经证明,计算2个2×2矩阵的乘积,7次乘法是必要的因此,要想进一步改进矩阵乘法的时间复杂性,就不能再寄希望于计算2×2矩阵的乘法次数的减少或许应当研究3×3或5×5矩阵的更好算法在Strassen之后又有许多算法改进了矩阵乘法的计算时间复杂性目前最好的计算时间上界是On
2.367而目前所知道的矩阵乘法的最好下界仍是它的平凡下界Ωn2因此到目前为止还无法确切知道矩阵乘法的时间复杂性关于这一研究课题还有许多工作可做。