矩阵相乘的快速算法.docx
《矩阵相乘的快速算法.docx》由会员分享，可在线阅读，更多相关《矩阵相乘的快速算法.docx（11页珍藏版）》请在冰点文库上搜索。
矩阵相乘的快速算法
矩阵相乘的快速算法
算法介绍
矩阵相乘在进行3D变换的时候是经常用到的。
在应用中常用矩阵相乘的定义算法对其进行计算。
这个算法用到了大量的循环和相乘运算，这使得算法效率不高。
而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度，所以对矩阵相乘算法进行一些改进是必要的。
这里要介绍的矩阵算法称为斯特拉森方法，它是由v.斯特拉森在1969年提出的一个方法。
我们先讨论二阶矩阵的计算方法。
对于二阶矩阵a11a12b11b12A=a21a22B=b21b22先计算下面7个量
（1）x1=（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的各元素为
（2）c11=a11*b11+a12*b21
c12=a11*b12+a12*b22
c21=a21*b11+a22*b21
c22=a21*b12+a22*b22
比较
（1）
（2），C的各元素可以表示为（3）c11=x1+x4-x5+x7c12=x3+x5c21=x2+x4c22=x1+x3-x2+x6
根据以上的方法，我们就可以计算4阶矩阵了，先将4阶矩阵A和B划分成四块2阶矩阵，分别利用公式计算它们的乘积，再使用
（1）（3）来计算出最后结果。
ma11ma12mb11mb12A4=ma21ma22B4=mb21mb22
其中a11a12a13a14b11b12b13b14ma11=a21a22ma12=a23a24mb11=b21b22mb12=b23b24a31a32a33a34b31b32b33b34ma21=a41a42ma22=a43a44mb21=b41b42mb22=b43b44
实现//计算2X2矩阵voidMultiply2X2（float&fOut_11,float&fOut_12,float&fOut_21,float&fOut_22,floatf1_11,floatf1_12,floatf1_21,floatf1_22,floatf2_11,floatf2_12,floatf2_21,floatf2_22）{constfloatx1（（f1_11+f1_22）*（f2_11+f2_22））;
constfloatx2（（f1_21+f1_22）*f2_11）;
constfloatx3（f1_11*（f2_12-f2_22））;
constfloatx4（f1_22*（f2_21-f2_11））;
constfloatx5（（f1_11+f1_12）*f2_22）;
constfloatx6（（f1_21-f1_11）*（f2_11+f2_12））;
constfloatx7（（f1_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矩阵voidMultiply（CLAYMATRIX&mOut,constCLAYMATRIX&m1,constCLAYMATRIX&m2）{floatfTmp[7][4];
//（ma11+ma22）*（mb11+mb22）
Multiply2X2（fTmp[0][0],fTmp[0][1],fTmp[0][2],fTmp[0][3],
m1._11+m1._33,m1._12+m1._34,m1._21+m1._43,m1._22+m1._44,
m2._11+m2._33,m2._12+m2._34,m2._21+m2._43,m2._22+m2._44）;
//（ma21+ma22）*mb11
Multiply2X2（fTmp[1][0],fTmp[1][1],fTmp[1][2],fTmp[1][3],
m1._31+m1._33,m1._32+m1._34,m1._41+m1._43,m1._42+m1._44,
m2._11,m2._12,m2._21,m2._22）;
//ma11*（mb12-mb22）
Multiply2X2（fTmp[2][0],fTmp[2][1],fTmp[2][2],fTmp[2][3],
m1._11,m1._12,m1._21,m1._22,
m2._13-m2._33,m2._14-m2._34,m2._23-m2._43,m2._24-m2._44）;
//ma22*（mb21-mb11）
Multiply2X2（fTmp[3][0],fTmp[3][1],fTmp[3][2],fTmp[3][3],
m1._33,m1._34,m1._43,m1._44,
m2._31-m2._11,m2._32-m2._12,m2._41-m2._21,m2._42-m2._22）;
//（ma11+ma12）*mb22
Multiply2X2（fTmp[4][0],fTmp[4][1],fTmp[4][2],fTmp[4][3],
m1._11+m1._13,m1._12+m1._14,m1._21+m1._23,m1._22+m1._24,
m2._33,m2._34,m2._43,m2._44）;
//（ma21-ma11）*（mb11+mb12）
Multiply2X2（fTmp[5][0],fTmp[5][1],fTmp[5][2],fTmp[5][3],
m1._31-m1._11,m1._32-m1._12,m1._41-m1._21,m1._42-m1._22,
m2._11+m2._13,m2._12+m2._14,m2._21+m2._23,m2._22+m2._24）;
//（ma12-ma22）*（mb21+mb22）
Multiply2X2（fTmp[6][0],fTmp[6][1],fTmp[6][2],fTmp[6][3],
m1._13-m1._33,m1._14-m1._34,m1._23-m1._43,m1._24-m1._44,
m2._31+m2._33,m2._32+m2._34,m2._41+m2._43,m2._42+m2._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阶矩阵：
原算法新算法
加法次数4872（48次加法，24次减法）
乘法次数6449
需要额外空间16*sizeof（float）28*sizeof（float）
新算法要比原算法多了24次减法运算，少了15次乘法。
但因为浮点乘法的运算速度要远远慢于加/减法运算，所以新算法的整体速度有所提高。
一、两个矩阵相乘的经典算法:
若设Q=M*N其中，M是m1*n1矩阵，N是m2*n2矩阵。
当n1=m2时有：
for（i=1;i
for（j=1;j<=n2;++j）{
Q[i][j]=0;
for（k=1;k<=n1;++k） Q[i][j]+=M[i][k]*N[k][j];
}
此算法的时间复杂度是O（m1*n1*n2）。
二、斯特拉森算法
斯特拉森方法，是由v.斯特拉森在1969年提出的一个方法。
我们先讨论二阶矩阵的计算方法。
对于二阶矩阵
a11a12b11b12
A=a21a22B=b21b22
先计算下面7个量
（1）
x1=（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的各元素为
（2）
c11=a11*b11+a12*b21
c12=a11*b12+a12*b22
c21=a21*b11+a22*b21
c22=a21*b12+a22*b22
比较
（1）
（2），C的各元素可以表示为（3）
c11=x1+x4-x5+x7
c12=x3+x5
c21=x2+x4
c22=x1+x3-x2+x6
根据以上的方法，我们就可以计算4阶矩阵了，先将4阶矩阵A和B划分成四块2阶矩阵，分别利用公式计算它们的乘积，再使用
（1）（3）来计算出最后结果。
ma11ma12mb11mb12
A4=ma21ma22B4=mb21mb22
其中
a11a12a13a14b11b12b13b14
ma11=a21a22ma12=a23a24mb11=b21b22mb12=b23b24
a31a32a33a34b31b32b33b34
ma21=a41a42ma22=a43a44mb21=b41b42mb22=b43b44
实现
//计算2X2矩阵
voidMultiply2X2（float&fOut_11,float&fOut_12,float&fOut_21,float&fOut_22,
floatf1_11,floatf1_12,floatf1_21,floatf1_22,
floatf2_11,floatf2_12,floatf2_21,floatf2_22）
{
constfloatx1（（f1_11+f1_22）*（f2_11+f2_22））;
constfloatx2（（f1_21+f1_22）*f2_11）;
constfloatx3（f1_11*（f2_12-f2_22））;
constfloatx4（f1_22*（f2_21-f2_11））;
constfloatx5（（f1_11+f1_12）*f2_22）;
constfloatx6（（f1_21-f1_11）*（f2_11+f2_12））;
constfloatx7（（f1_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矩阵
voidMultiply（CLAYMATRIX&mOut,constCLAYMATRIX&m1,constCLAYMATRIX&m2）
{
floatfTmp[7][4];
//（ma11+ma22）*（mb11+mb22）
Multiply2X2（fTmp[0][0],fTmp[0][1],fTmp[0][2],fTmp[0][3],
m1._11+m1._33,m1._12+m1._34,m1._21+m1._43,m1._22+m1._44,
m2._11+m2._33,m2._12+m2._34,m2._21+m2._43,m2._22+m2._44）;
//（ma21+ma22）*mb11
Multiply2X2（fTmp[1][0],fTmp[1][1],fTmp[1][2],fTmp[1][3],
m1._31+m1._33,m1._32+m1._34,m1._41+m1._43,m1._42+m1._44,
m2._11,m2._12,m2._21,m2._22）;
//ma11*（mb12-mb22）
Multiply2X2（fTmp[2][0],fTmp[2][1],fTmp[2][2],fTmp[2][3],
m1._11,m1._12,m1._21,m1._22,
m2._13-m2._33,m2._14-m2._34,m2._23-m2._43,m2._24-m2._44）;
//ma22*（mb21-mb11）
Multiply2X2（fTmp[3][0],fTmp[3][1],fTmp[3][2],fTmp[3][3],
m1._33,m1._34,m1._43,m1._44,
m2._31-m2._11,m2._32-m2._12,m2._41-m2._21,m2._42-m2._22）;
//（ma11+ma12）*mb22
Multiply2X2（fTmp[4][0],fTmp[4][1],fTmp[4][2],fTmp[4][3],
m1._11+m1._13,m1._12+m1._14,m1._21+m1._23,m1._22+m1._24,
m2._33,m2._34,m2._43,m2._44）;
//（ma21-ma11）*（mb11+mb12）
Multiply2X2（fTmp[5][0],fTmp[5][1],fTmp[5][2],fTmp[5][3],
m1._31-m1._11,m1._32-m1._12,m1._41-m1._21,m1._42-m1._22,
m2._11+m2._13,m2._12+m2._14,m2._21+m2._23,m2._22+m2._24）;
//（ma12-ma22）*（mb21+mb22）
Multiply2X2（fTmp[6][0],fTmp[6][1],fTmp[6][2],fTmp[6][3],
m1._13-m1._33,m1._14-m1._34,m1._23-m1._43,m1._24-m1._44,
m2._31+m2._33,m2._32+m2._34,m2._41+m2._43,m2._42+m2._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阶矩阵：
原算法新算法
加法次数4872（48次加法，24次减法）
乘法次数6449
需要额外空间16*sizeof（float）28*sizeof（float）
新算法要比原算法多了24次减法运算，少了15次乘法。
但因为浮点乘法的运算速度要远远慢于加/减法运算，所以新算法的整体速度有所提高。
三、Strassen矩阵乘法
矩阵乘法是线性代数中最常见的运算之一，它在数值计算中有广泛的应用。
若A和B是2个n×n的矩阵，则它们的乘积C=AB同样是一个n×n的矩阵。
A和B的乘积矩阵C中的元素C[i,j]定义为:
若依此定义来计算A和B的乘积矩阵C，则每计算C的一个元素C[i,j]，需要做n个乘法和n-1次加法。
因此，求出矩阵C的n2个元素所需的计算时间为0（n3）。
60年代末，Strassen采用了类似于在大整数乘法中用过的分治技术，将计算2个n阶矩阵乘积所需的计算时间改进到O（nlog7）=O（n2.18）。
首先，我们还是需要假设n是2的幂。
将矩阵A，B和C中每一矩阵都分块成为4个大小相等的子矩阵，每个子矩阵都是n/2×n/2的方阵。
由此可将方程C=AB重写为:
（1）
由此可得:
C11=A11B11+A12B21
（2）
C12=A11B12+A12B22 （3）
C21=A21B11+A22B21 （4）
C22=A21B12+A22B22 （5）
如果n=2，则2个2阶方阵的乘积可以直接用
（2）-（3）式计算出来，共需8次乘法和4次加法。
当子矩阵的阶大于2时，为求2个子矩阵的积，可以继续将子矩阵分块，直到子矩阵的阶降为2。
这样，就产生了一个分治降阶的递归算法。
依此算法，计算2个n阶方阵的乘积转化为计算8个n/2阶方阵的乘积和4个n/2阶方阵的加法。
2个n/2×n/2矩阵的加法显然可以在c*n2/4时间完成，这里c是一个常数。
因此，上述分治法的计算时间耗费T（n）应该满足：
这个递归方程的解仍然是T（n）=O（n3）。
因此，该方法并不比用原始定义直接计算更有效。
究其原因，乃是由于式
（2）-（5）并没有减少矩阵的乘法次数。
而矩阵乘法耗费的时间要比矩阵加减法耗费的时间多得多。
要想改进矩阵乘法的计算时间复杂性，必须减少子矩阵乘法运算的次数。
按照上述分治法的思想可以看出，要想减少乘法运算次数，关键在于计算2个2阶方阵的乘积时，能否用少于8次的乘法运算。
Strassen提出了一种新的算法来计算2个2阶方阵的乘积。
他的算法只用了7次乘法运算，但增加了加、减法的运算次数。
这7次乘法是:
M1=A11（B12-B22）
M2=（A11+A12）B2
升级会员 