整理矩阵相乘的快速算法.docx
- 文档编号:28964184
- 上传时间:2023-07-20
- 格式:DOCX
- 页数:15
- 大小:25.41KB
整理矩阵相乘的快速算法.docx
《整理矩阵相乘的快速算法.docx》由会员分享,可在线阅读,更多相关《整理矩阵相乘的快速算法.docx(15页珍藏版)》请在冰豆网上搜索。
整理矩阵相乘的快速算法
矩阵相乘的快速算法
算法介绍
矩阵相乘在进行3D变换的时候是经常用到的。
在应用中常用矩阵相乘的定义算法对其进行计算。
这个算法用到了大量的循环和相乘运算,这使得算法效率不高。
而矩阵相乘的计算效率很大程度上的影响了整个程序的运行速度,所以对矩阵相乘算法进行一些改进是必要的。
这里要介绍的矩阵算法称为斯特拉森方法,它是由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次乘法。
但因为浮点乘法的运算速度要远远慢于加/减法运算,所以新算法的整体速度有所提高。
一、两个矩阵相乘的经典算法:
若设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阶方阵
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 整理 矩阵 相乘 快速 算法