首页 > 代码库 > 算法导论-矩阵乘法-strassen算法
算法导论-矩阵乘法-strassen算法
目录
1、矩阵相乘的朴素算法
2、矩阵相乘的strassen算法
3、完整测试代码c++
4、性能分析
5、参考资料
内容
1、矩阵相乘的朴素算法 T(n) = Θ(n3)
朴素矩阵相乘算法，思想明了，编程实现简单。时间复杂度是Θ(n^3)。伪码如下
1 for i ← 1 to n2 do for j ← 1 to n3 do c[i][j] ← 04 for k ← 1 to n5 do c[i][j] ← c[i][j] + a[i][k]⋅ b[k][j]
2、矩阵相乘的strassen算法 T(n)=Θ(nlog7) =Θ (n2.81)
矩阵乘法中采用分治法，第一感觉上应该能够有效的提高算法的效率。如下图所示分治法方案，以及对该算法的效率分析。有图可知，算法效率是Θ(n^3)。算法效率并没有提高。下面介绍下矩阵分治法思想：
鉴于上面的分治法方案无法有效提高算法的效率，要想提高算法效率，由主定理方法可知必须想办法将2中递归式中的系数8减少。Strassen提出了一种将系数减少到7的分治法方案，如下图所示。
效率分析如下：
伪码如下：
1 Strassen (N,MatrixA,MatrixB,MatrixResult) 2 3 //splitting input Matrixes, into 4 submatrices each. 4 for i <- 0 to N/2 5 for j <- 0 to N/2 6 A11[i][j] <- MatrixA[i][j]; //a矩阵块 7 A12[i][j] <- MatrixA[i][j + N / 2]; //b矩阵块 8 A21[i][j] <- MatrixA[i + N / 2][j]; //c矩阵块 9 A22[i][j] <- MatrixA[i + N / 2][j + N / 2];//d矩阵块10 11 B11[i][j] <- MatrixB[i][j]; //e 矩阵块12 B12[i][j] <- MatrixB[i][j + N / 2]; //f 矩阵块13 B21[i][j] <- MatrixB[i + N / 2][j]; //g 矩阵块14 B22[i][j] <- MatrixB[i + N / 2][j + N / 2]; //h矩阵块15 //here we calculate M1..M7 matrices . 17 //递归求M118 HalfSize <- N/2 19 AResult <- A11+A2220 BResult <- B11+B22 21 Strassen( HalfSize, AResult, BResult, M1 ); //M1=(A11+A22)*(B11+B22) p5=(a+d)*(e+h) 22 //递归求M223 AResult <- A21+A22 24 Strassen(HalfSize, AResult, B11, M2); //M2=(A21+A22)B11 p3=(c+d)*e25 //递归求M326 BResult <- B12 - B22 27 Strassen(HalfSize, A11, BResult, M3); //M3=A11(B12-B22) p1=a*(f-h)28 //递归求M429 BResult <- B21 - B11 30 Strassen(HalfSize, A22, BResult, M4); //M4=A22(B21-B11) p4=d*(g-e)31 //递归求M532 AResult <- A11+A12 33 Strassen(HalfSize, AResult, B22, M5); //M5=(A11+A12)B22 p2=(a+b)*h34 //递归求M635 AResult <- A21-A1136 BResult <- B11+B12 37 Strassen( HalfSize, AResult, BResult, M6); //M6=(A21-A11)(B11+B12) p7=(c-a)(e+f)38 //递归求M739 AResult <- A12-A2240 BResult <- B21+B22 41 Strassen(HalfSize, AResult, BResult, M7); //M7=(A12-A22)(B21+B22) p6=(b-d)*(g+h)42 43 //计算结果子矩阵44 C11 <- M1 + M4 - M5 + M7;45 46 C12 <- M3 + M5;47 48 C21 <- M2 + M4;49 50 C22 <- M1 + M3 - M2 + M6;51 //at this point , we have calculated the c11..c22 matrices, and now we are going to52 //put them together and make a unit matrix which would describe our resulting Matrix.53 for i <- 0 to N/254 for j <- 0 to N/255 MatrixResult[i][j] <- C11[i][j];56 MatrixResult[i][j + N / 2] <- C12[i][j];57 MatrixResult[i + N / 2][j] <- C21[i][j];58 MatrixResult[i + N / 2][j + N / 2] <- C22[i][j];
3、完成测试代码
Strassen.h
1 #ifndef STRASSEN_HH 2 #define STRASSEN_HH 3 template<typename T> 4 class Strassen_class{ 5 public: 6 void ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize ); 7 void SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize ); 8 void MUL( T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize );//朴素算法实现 9 void FillMatrix( T** MatrixA, T** MatrixB, int length);//A,B矩阵赋值 10 void PrintMatrix(T **MatrixA,int MatrixSize);//打印矩阵 11 void Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC);//Strassen算法实现 12 }; 13 template<typename T> 14 void Strassen_class<T>::ADD(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize ) 15 { 16 for ( int i = 0; i < MatrixSize; i++) 17 { 18 for ( int j = 0; j < MatrixSize; j++) 19 { 20 MatrixResult[i][j] = MatrixA[i][j] + MatrixB[i][j]; 21 } 22 } 23 } 24 template<typename T> 25 void Strassen_class<T>::SUB(T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize ) 26 { 27 for ( int i = 0; i < MatrixSize; i++) 28 { 29 for ( int j = 0; j < MatrixSize; j++) 30 { 31 MatrixResult[i][j] = MatrixA[i][j] - MatrixB[i][j]; 32 } 33 } 34 } 35 template<typename T> 36 void Strassen_class<T>::MUL( T** MatrixA, T** MatrixB, T** MatrixResult, int MatrixSize ) 37 { 38 for (int i=0;i<MatrixSize ;i++) 39 { 40 for (int j=0;j<MatrixSize ;j++) 41 { 42 MatrixResult[i][j]=0; 43 for (int k=0;k<MatrixSize ;k++) 44 { 45 MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j]; 46 } 47 } 48 } 49 } 50 51 /* 52 c++使用二维数组，申请动态内存方法 53 申请 54 int **A; 55 A = new int *[desired_array_row]; 56 for ( int i = 0; i < desired_array_row; i++) 57 A[i] = new int [desired_column_size]; 58 59 释放 60 for ( int i = 0; i < your_array_row; i++) 61 delete [] A[i]; 62 delete[] A; 63 64 */ 65 template<typename T> 66 void Strassen_class<T>::Strassen(int N, T **MatrixA, T **MatrixB, T **MatrixC) 67 { 68 69 int HalfSize = N/2; 70 int newSize = N/2; 71 72 if ( N <= 64 ) //分治门槛，小于这个值时不再进行递归计算，而是采用常规矩阵计算方法 73 { 74 MUL(MatrixA,MatrixB,MatrixC,N); 75 } 76 else 77 { 78 T** A11; 79 T** A12; 80 T** A21; 81 T** A22; 82 83 T** B11; 84 T** B12; 85 T** B21; 86 T** B22; 87 88 T** C11; 89 T** C12; 90 T** C21; 91 T** C22; 92 93 T** M1; 94 T** M2; 95 T** M3; 96 T** M4; 97 T** M5; 98 T** M6; 99 T** M7;100 T** AResult;101 T** BResult;102 103 //making a 1 diminsional pointer based array.104 A11 = new T *[newSize];105 A12 = new T *[newSize];106 A21 = new T *[newSize];107 A22 = new T *[newSize];108 109 B11 = new T *[newSize];110 B12 = new T *[newSize];111 B21 = new T *[newSize];112 B22 = new T *[newSize];113 114 C11 = new T *[newSize];115 C12 = new T *[newSize];116 C21 = new T *[newSize];117 C22 = new T *[newSize];118 119 M1 = new T *[newSize];120 M2 = new T *[newSize];121 M3 = new T *[newSize];122 M4 = new T *[newSize];123 M5 = new T *[newSize];124 M6 = new T *[newSize];125 M7 = new T *[newSize];126 127 AResult = new T *[newSize];128 BResult = new T *[newSize];129 130 int newLength = newSize;131 132 //making that 1 diminsional pointer based array , a 2D pointer based array133 for ( int i = 0; i < newSize; i++)134 {135 A11[i] = new T[newLength];136 A12[i] = new T[newLength];137 A21[i] = new T[newLength];138 A22[i] = new T[newLength];139 140 B11[i] = new T[newLength];141 B12[i] = new T[newLength];142 B21[i] = new T[newLength];143 B22[i] = new T[newLength];144 145 C11[i] = new T[newLength];146 C12[i] = new T[newLength];147 C21[i] = new T[newLength];148 C22[i] = new T[newLength];149 150 M1[i] = new T[newLength];151 M2[i] = new T[newLength];152 M3[i] = new T[newLength];153 M4[i] = new T[newLength];154 M5[i] = new T[newLength];155 M6[i] = new T[newLength];156 M7[i] = new T[newLength];157 158 AResult[i] = new T[newLength];159 BResult[i] = new T[newLength];160 161 162 }163 //splitting input Matrixes, into 4 submatrices each.164 for (int i = 0; i < N / 2; i++)165 {166 for (int j = 0; j < N / 2; j++)167 {168 A11[i][j] = MatrixA[i][j];169 A12[i][j] = MatrixA[i][j + N / 2];170 A21[i][j] = MatrixA[i + N / 2][j];171 A22[i][j] = MatrixA[i + N / 2][j + N / 2];172 173 B11[i][j] = MatrixB[i][j];174 B12[i][j] = MatrixB[i][j + N / 2];175 B21[i][j] = MatrixB[i + N / 2][j];176 B22[i][j] = MatrixB[i + N / 2][j + N / 2];177 178 }179 }180 181 //here we calculate M1..M7 matrices .182 //M1[][]183 ADD( A11,A22,AResult, HalfSize);184 ADD( B11,B22,BResult, HalfSize); //p5=(a+d)*(e+h)185 Strassen( HalfSize, AResult, BResult, M1 ); //now that we need to multiply this , we use the strassen itself .186 187 188 //M2[][]189 ADD( A21,A22,AResult, HalfSize); //M2=(A21+A22)B11 p3=(c+d)*e190 Strassen(HalfSize, AResult, B11, M2); //Mul(AResult,B11,M2);191 192 //M3[][]193 SUB( B12,B22,BResult, HalfSize); //M3=A11(B12-B22) p1=a*(f-h)194 Strassen(HalfSize, A11, BResult, M3); //Mul(A11,BResult,M3);195 196 //M4[][]197 SUB( B21, B11, BResult, HalfSize); //M4=A22(B21-B11) p4=d*(g-e)198 Strassen(HalfSize, A22, BResult, M4); //Mul(A22,BResult,M4);199 200 //M5[][]201 ADD( A11, A12, AResult, HalfSize); //M5=(A11+A12)B22 p2=(a+b)*h202 Strassen(HalfSize, AResult, B22, M5); //Mul(AResult,B22,M5);203 204 205 //M6[][]206 SUB( A21, A11, AResult, HalfSize);207 ADD( B11, B12, BResult, HalfSize); //M6=(A21-A11)(B11+B12) p7=(c-a)(e+f)208 Strassen( HalfSize, AResult, BResult, M6); //Mul(AResult,BResult,M6);209 210 //M7[][]211 SUB(A12, A22, AResult, HalfSize);212 ADD(B21, B22, BResult, HalfSize); //M7=(A12-A22)(B21+B22) p6=(b-d)*(g+h)213 Strassen(HalfSize, AResult, BResult, M7); //Mul(AResult,BResult,M7);214 215 //C11 = M1 + M4 - M5 + M7;216 ADD( M1, M4, AResult, HalfSize);217 SUB( M7, M5, BResult, HalfSize);218 ADD( AResult, BResult, C11, HalfSize);219 220 //C12 = M3 + M5;221 ADD( M3, M5, C12, HalfSize);222 223 //C21 = M2 + M4;224 ADD( M2, M4, C21, HalfSize);225 226 //C22 = M1 + M3 - M2 + M6;227 ADD( M1, M3, AResult, HalfSize);228 SUB( M6, M2, BResult, HalfSize);229 ADD( AResult, BResult, C22, HalfSize);230 231 //at this point , we have calculated the c11..c22 matrices, and now we are going to232 //put them together and make a unit matrix which would describe our resulting Matrix.233 //组合小矩阵到一个大矩阵234 for (int i = 0; i < N/2 ; i++)235 {236 for (int j = 0 ; j < N/2 ; j++)237 {238 MatrixC[i][j] = C11[i][j];239 MatrixC[i][j + N / 2] = C12[i][j];240 MatrixC[i + N / 2][j] = C21[i][j];241 MatrixC[i + N / 2][j + N / 2] = C22[i][j];242 }243 }244 245 // 释放矩阵内存空间246 for (int i = 0; i < newLength; i++)247 {248 delete[] A11[i];delete[] A12[i];delete[] A21[i];249 delete[] A22[i];250 251 delete[] B11[i];delete[] B12[i];delete[] B21[i];252 delete[] B22[i];253 delete[] C11[i];delete[] C12[i];delete[] C21[i];254 delete[] C22[i];255 delete[] M1[i];delete[] M2[i];delete[] M3[i];delete[] M4[i];256 delete[] M5[i];delete[] M6[i];delete[] M7[i];257 delete[] AResult[i];delete[] BResult[i] ;258 }259 delete[] A11;delete[] A12;delete[] A21;delete[] A22;260 delete[] B11;delete[] B12;delete[] B21;delete[] B22;261 delete[] C11;delete[] C12;delete[] C21;delete[] C22;262 delete[] M1;delete[] M2;delete[] M3;delete[] M4;delete[] M5;263 delete[] M6;delete[] M7;264 delete[] AResult;265 delete[] BResult ;266 267 }//end of else268 269 }270 271 template<typename T>272 void Strassen_class<T>::FillMatrix( T** MatrixA, T** MatrixB, int length)273 {274 for(int row = 0; row<length; row++)275 {276 for(int column = 0; column<length; column++)277 {278 279 MatrixB[row][column] = (MatrixA[row][column] = rand() %5);280 //matrix2[row][column] = rand() % 2;//ba hazfe in khat 50% afzayeshe soorat khahim dasht281 }282 283 }284 }285 template<typename T>286 void Strassen_class<T>::PrintMatrix(T **MatrixA,int MatrixSize)287 {288 cout<<endl;289 for(int row = 0; row<MatrixSize; row++)290 {291 for(int column = 0; column<MatrixSize; column++)292 {293 294 295 cout<<MatrixA[row][column]<<"\t";296 if ((column+1)%((MatrixSize)) == 0)297 cout<<endl;298 }299 300 }301 cout<<endl;302 }303 #endif
Strassen.cpp
1 #include <iostream> 2 #include <ctime> 3 #include <Windows.h> 4 using namespace std; 5 #include "Strassen.h" 6 7 int main() 8 { 9 Strassen_class<int> stra;//定义Strassen_class类对象10 int MatrixSize = 0;11 12 int** MatrixA; //存放矩阵A13 int** MatrixB; //存放矩阵B14 int** MatrixC; //存放结果矩阵15 16 clock_t startTime_For_Normal_Multipilication ;17 clock_t endTime_For_Normal_Multipilication ;18 19 clock_t startTime_For_Strassen ;20 clock_t endTime_For_Strassen ;21 srand(time(0));22 23 cout<<"\n请输入矩阵大小(必须是2的幂指数值(例如:32,64,512,..): ";24 cin>>MatrixSize;25 26 int N = MatrixSize;//for readiblity.27 28 //申请内存29 MatrixA = new int *[MatrixSize];30 MatrixB = new int *[MatrixSize];31 MatrixC = new int *[MatrixSize];32 33 for (int i = 0; i < MatrixSize; i++)34 {35 MatrixA[i] = new int [MatrixSize];36 MatrixB[i] = new int [MatrixSize];37 MatrixC[i] = new int [MatrixSize];38 }39 40 stra.FillMatrix(MatrixA,MatrixB,MatrixSize); //矩阵赋值41 42 //*******************conventional multiplication test43 cout<<"朴素矩阵算法开始时钟: "<< (startTime_For_Normal_Multipilication = clock());44 45 stra.MUL(MatrixA,MatrixB,MatrixC,MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)46 47 cout<<"\n朴素矩阵算法结束时钟: "<< (endTime_For_Normal_Multipilication = clock());48 49 cout<<"\n矩阵运算结果... \n";50 stra.PrintMatrix(MatrixC,MatrixSize);51 52 //*******************Strassen multiplication test53 cout<<"\nStrassen算法开始时钟: "<< (startTime_For_Strassen = clock());54 55 stra.Strassen( N, MatrixA, MatrixB, MatrixC ); //strassen矩阵相乘算法56 57 cout<<"\nStrassen算法结束时钟: "<<(endTime_For_Strassen = clock());58 59 60 cout<<"\n矩阵运算结果... \n";61 stra.PrintMatrix(MatrixC,MatrixSize);62 63 cout<<"矩阵大小 "<<MatrixSize;64 cout<<"\n朴素矩阵算法: "<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)<<" Clocks.."<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)/CLOCKS_PER_SEC<<" Sec";65 cout<<"\nStrassen算法:"<<(endTime_For_Strassen - startTime_For_Strassen)<<" Clocks.."<<(endTime_For_Strassen - startTime_For_Strassen)/CLOCKS_PER_SEC<<" Sec\n";66 system("Pause");67 return 0;68 69 }
输出：
4、性能分析
| 矩阵大小 | 朴素矩阵算法（秒） | Strassen算法（秒） |
| 32 | 0.003 | 0.003 |
| 64 | 0.004 | 0.004 |
| 128 | 0.021 | 0.071 |
| 256 | 0.09 | 0.854 |
| 512 | 0.782 | 6.408 |
| 1024 | 8.908 | 52.391 |
可以发现：可以看到使用Strassen算法时，耗时不但没有减少，反而剧烈增多，在n=512时计算时间就无法忍受，效果没有朴素矩阵算法好。网上查阅资料，现罗列如下：
1）采用Strassen算法作递归运算，需要创建大量的动态二维数组，其中分配堆内存空间将占用大量计算时间，从而掩盖了Strassen算法的优势
2）于是对Strassen算法做出改进，设定一个界限。当n<界限时，使用普通法计算矩阵，而不继续分治递归。需要合理设置界限，不同环境（硬件配置）下界限不同
3）矩阵乘法一般意义上还是选择的是朴素的方法，只有当矩阵变稠密，而且矩阵的阶数很大时，才会考虑使用Strassen算法。
分析原因：（网上总结的说法）
http://blog.csdn.net/handawnc/article/details/7987107
仔细研究后发现，采用Strassen算法作递归运算，需要创建大量的动态二维数组，其中分配堆内存空间将占用大量计算时间，从而掩盖了Strassen算法的优势。于是对Strassen算法做出改进，设定一个界限。当n<界限时，使用普通法计算矩阵，而不继续分治递归。
改进后算法优势明显，就算时间大幅下降。之后，针对不同大小的界限进行试验。在初步试验中发现，当数据规模小于1000时，下界S法的差别不大，规模大于1000以后，n取值越大，消耗时间下降。最优的界限值在32～128之间。
因为计算机每次运算时的系统环境不同（CPU占用、内存占用等），所以计算出的时间会有一定浮动。虽然这样，试验结果已经能得出结论Strassen算法比常规法优势明显。使用下界法改进后，在分治效率和动态分配内存间取舍，针对不同的数据规模稍加试验可以得到一个最优的界限。
http://www.cppblog.com/sosi/archive/2010/08/30/125259.html
时间复杂度就马上降下来了。。但是不要过于乐观。
从实用的观点看，Strassen算法通常不是矩阵乘法所选择的方法：
1 在Strassen算法的运行时间中，隐含的常数因子比简单的O(n^3)方法常数因子大
2 当矩阵是稀疏的时候，为稀疏矩阵设计的算法更快
3 Strassen算法不像简单方法那样子具有数值稳定性
4 在递归层次中生成的子矩阵要消耗空间。
所以矩阵乘法一般意义上还是选择的是朴素的方法，只有当矩阵变稠密，而且矩阵的阶数>20左右，才会考虑使用Strassen算法。
5、参考资料
【1】http://blog.csdn.net/xyd0512/article/details/8220506
【2】http://blog.csdn.net/zhuangxiaobin/article/details/36476769
【3】http://blog.csdn.net/handawnc/article/details/7987107
【4】http://www.xuebuyuan.com/552410.html
【5】http://blog.csdn.net/chenhq1991/article/details/7599824
算法导论-矩阵乘法-strassen算法
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。