Strassen矩阵乘法

实验报告一、实验目的通过算法的程序实现和执行时间测试、并与理论上的结论进行对比分析，深入理解算法时间复杂度分析中对于输入数据考虑其等价类的意义，理解算法时间复杂度渐进性态和和增长率的概念，为后续学习和实验奠定基础，同时也学习程序效率测试的基本思路。二、实验内容与实验步骤（1）了解分治的基本思想并编写算法解决Strassen矩阵乘法问题。（2）打开一台装有MyEclipse-10.7.1的PC。（3）把已经写好的代码在Java环境下运行并调试。（4）记录运行结果。三、实验环境Windows7家庭普通版，MyEclipse-10.7.1四、阐述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重写为:课程名称：算法设计与分析班级：12软工一班实验成绩：实验名称：Strassen矩阵乘法学号：1242159101批阅教师签字：实验编号：实验一姓名：陈双飞实验日期：2014年12月14日指导教师：陈平组号：实验时间：时分－时分(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)B22M3=(A21+A22)B11M4=A22(B21-B11)M5=(A11+A22)(B11+B22)M6=(A12-A22)(B21+B22)M7=(A11-A21)(B11+B12)做了这7次乘法后，再做若干次加、减法就可以得到:C11=M5+M4-M2+M6C12=M1+M2C21=M3+M4C22=M5+M1-M3-M7以上计算的正确性很容易验证。例如:C22=M5+M1-M3-M7=(A11+A22)(B11+B22)+A11(B12-B22)-(A21+A22)B11-(A11-A21)(B11+B12)=A11B11+A11B22+A22B11+A22B22+A11B12-A11B22-A21B11-A22B11-A11B11-A11B12+A21B11+A21B12=A21B12+A22B22由(2)式便知其正确性。六、问题解决我们可以得到完整的Strassen算法如下：procedureSTRASSEN(n,A,B,C);beginifn=2thenMATRIX-MULTIPLY(A，B，C)elsebegin将矩阵A和B依(1)式分块;STRASSEN(n/2,A11,B12-B22,M1);STRASSEN(n/2,A11+A12,B22,M2);STRASSEN(n/2,A21+A22,B11,M3);STRASSEN(n/2,A22,B21-B11,M4);STRASSEN(n/2,A11+A22,B11+B22,M5);STRASSEN(n/2,A12-A22,B21+B22,M6);STRASSEN(n/2,A11-A21,B11+B12,M7);;end;end;其中MATRIX-MULTIPLY(A，B，C)是按通常的矩阵乘法计算C=AB的子算法。Strassen矩阵乘积分治算法中，用了7次对于n/2阶矩阵乘积的递归调用和18次n/2阶矩阵的加减运算。由此可知，该算法的所需的计算时间T(n)满足如下的递归方程:按照解递归方程的套用公式法，其解为T(n)=O(nlog7)≈O(n2.81)。由此可见，Strassen矩阵乘法的计算时间复杂性比普通矩阵乘法有阶的改进。有人曾列举了计算2个2阶矩阵乘法的36种不同方法。但所有的方法都要做7次乘法。除非能找到一种计算2阶方阵乘积的算法，使乘法的计算次数少于7次，按上述思路才有可能进一步改进矩阵乘积的计算时间的上界。但是Hopcroft和Kerr(197l)已经证明，计算2个2×2矩阵的乘积，7次乘法是必要的。因此，要想进一步改进矩阵乘法的时间复杂性，就不能再寄希望于计算2×2矩阵的乘法次数的减少。或许应当研究3×3或5×5矩阵的更好算法。在Strassen之后又有许多算法改进了矩阵乘法的计算时间复杂性。目前最好的计算时间上界是O(n2.367)。而目前所知道的矩阵乘法的最好下界仍是它的平凡下界Ω(n2)。因此到目前为止还无法确切知道矩阵乘法的时间复杂性。关于这一研究课题还有许多工作可做。七、实验结果总结实现算法代码：/***Strassen矩阵乘法**/importjava.util.*;publicclassStrassen{publicStrassen(){A=newint[NUMBER][NUMBER];B=newint[NUMBER][NUMBER];C=newint[NUMBER][NUMBER];}/***输入矩阵函数**/publicvoidinput(inta[][]){Scannerscanner=newScanner(System.in);for(inti=0;ia.length;i++){for(intj=0;ja[i].length;j++){a[i][j]=scanner.nextInt();}}}/***输出矩阵**/publicvoidoutput(int[][]resault){for(intb[]:resault){for(inttemp:b){System.out.print(temp+);}System.out.println();}}/***矩阵乘法，此处只是定义了2*2矩阵的乘法**/publicvoidMul(int[][]first,int[][]second,int[][]resault){for(inti=0;i2;++i){for(intj=0;j2;++j){resault[i][j]=0;for(intk=0;k2;++k){resault[i][j]+=first[i][k]*second[k][j];}}}}/***矩阵的加法运算**/publicvoidAdd(int[][]first,int[][]second,int[][]resault){for(inti=0;ifirst.length;i++){for(intj=0;jfirst[i].length;j++){resault[i][j]=first[i][j]+second[i][j];}}}/***矩阵的减法运算**/publicvoidsub(int[][]first,int[][]second,int[][]resault){for(inti=0;ifirst.length;i++){for(intj=0;jfirst[i].length;j++){resault[i][j]=first[i][j]-second[i][j];}}}/***strassen矩阵算法**/publicvoidstrassen(int[][]A,int[][]B,int[][]C){//定义一些中间变量int[][]M1=newint[NUMBER][NUMBER];int[][]M2=newint[NUMBER][NUMBER];int[][]M3=newint[NUMBER][NUMBER];int[][]M4=newint[NUMBER][NUMBER];int[][]M5=newint[NUMBER][NUMBER];int[][]M6=newint[NUMBER][NUMBER];int[][]M7=newint[NUMBER][NUMBER];int[][]C11=newint[NUMBER][NUMBER];int[][]C12=newint[NUMBER][NUMBER];int[][]C21=newint[NUMBER][NUMBER];int[][]C22=newint[NUMBER][NUMBER];int[][]A11=newint[NUMBER][NUMBER];int[][]A12=newint[NUMBER][NUMBER];int[][]A21=newint[NUMBER][NUMBER];int[][]A22=newint[NUMBER][NUMBER];int[][]B11=newint[NUMBER][NUMBER];int[][]B12=newint[NUMBER][NUMBER];int[][]B21=newint[NUMBER][NUMBER];int[][]B22=newint[NUMBER][NUMBER];int[][]temp=newint[NUMBER][NUMBER];int[][]temp1=newint[NUMBER][NUMBER];if(A.length==2){Mul(A,B,C);}else{//首先将矩阵A，B分为4块for(inti=0;iA.length/2;i++){for(intj=0;jA.length/2;j++){A11[i][j]=A[i][j];A12[i][j]=A[i][j+A.length/2];A21[i][j]=A[i+A.length/2][j];A22[i][j]=A[i+A.length/2][j+A.length/2];B11[i][j]=B[i][j];B12[i][j]=B[i][j+A.length/2];B21[i][j]=B[i+A.length/2][j];B22[i][j]=B[i+A.length/2][j+A.length/2];}}//计算M1sub(B12,B22,temp);Mul(A11,temp,M1);//计算M2Add(A11,A12,temp);Mul(temp,B22,M2);//计算M3Add(A21,A22,temp);Mul(temp,B11,M3);//M4sub(B21,B11,temp);Mul(A22,temp,M4);//M5Add(A11,A22,temp1);Add(B11,B22,temp);Mul(temp1,temp,M5);//M6sub(A12,A22,temp1);Add(B21,B22,temp);Mul(temp1,temp,M6);//M7sub(A11,A2