算法导论-矩阵乘法-strassen算法

内容

1、矩阵相乘的朴素算法 T(n) = Θ(n3)

朴素矩阵相乘算法，思想明了，编程实现简单。时间复杂度是Θ(n^3)。伪码如下

1 for i ← 1 to n
2     do for j ← 1 to n
3         do c[i][j] ← 0
4             for k ← 1 to n
5                 do c[i][j] ← c[i][j] + a[i][k]⋅ b[k][j]

2、矩阵相乘的strassen算法 T(n)=Θ(n^log7) =Θ (n^2.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             //递归求M1
18             HalfSize  <-  N/2    
19             AResult  <-  A11+A22
20             BResult  <-  B11+B22                                                                     
21             Strassen( HalfSize, AResult, BResult, M1 );   //M1=(A11+A22)*(B11+B22)          p5=(a+d)*(e+h)    
22             //递归求M2
23             AResult  <-  A21+A22    
24             Strassen(HalfSize, AResult, B11, M2);          //M2=(A21+A22)B11                 p3=(c+d)*e
25             //递归求M3
26             BResult  <-  B12 - B22   
27             Strassen(HalfSize, A11, BResult, M3);         //M3=A11(B12-B22)                  p1=a*(f-h)
28             //递归求M4
29             BResult  <-  B21 - B11  
30             Strassen(HalfSize, A22, BResult, M4);         //M4=A22(B21-B11)                  p4=d*(g-e)
31             //递归求M5
32             AResult  <-  A11+A12    
33             Strassen(HalfSize, AResult, B22, M5);         //M5=(A11+A12)B22                  p2=(a+b)*h
34             //递归求M6
35             AResult  <-  A21-A11
36             BResult  <-  B11+B12      
37             Strassen( HalfSize, AResult, BResult, M6);     //M6=(A21-A11)(B11+B12)          p7=(c-a)(e+f)
38             //递归求M7
39             AResult  <-  A12-A22
40             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 to
52             //put them together and make a unit matrix which would describe our resulting Matrix.
53             for i  <-  0  to  N/2
54                 for j  <-  0  to  N/2
55                     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 array
133         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)*e
190         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)*h
202         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 to
232         //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 else
268 
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 dasht
281         }
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]<<"	";
296             if ((column+1)%((MatrixSize)) == 0)
297                 cout<<endl;
298         }
299 
300     }
301     cout<<endl;
302 }
303 #endif

Strassen.h

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;    //存放矩阵A
13     int** MatrixB;    //存放矩阵B
14     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<<"
请输入矩阵大小(必须是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 test
43         cout<<"朴素矩阵算法开始时钟:  "<< (startTime_For_Normal_Multipilication = clock());
44 
45         stra.MUL(MatrixA,MatrixB,MatrixC,MatrixSize);//朴素矩阵相乘算法 T(n) = O(n^3)
46 
47         cout<<"
朴素矩阵算法结束时钟: "<< (endTime_For_Normal_Multipilication = clock());
48 
49         cout<<"
矩阵运算结果... 
";
50         stra.PrintMatrix(MatrixC,MatrixSize);
51 
52   //*******************Strassen multiplication test
53         cout<<"
Strassen算法开始时钟: "<< (startTime_For_Strassen = clock());
54 
55         stra.Strassen( N, MatrixA, MatrixB, MatrixC ); //strassen矩阵相乘算法
56 
57         cout<<"
Strassen算法结束时钟: "<<(endTime_For_Strassen = clock());
58 
59 
60     cout<<"
矩阵运算结果... 
";
61     stra.PrintMatrix(MatrixC,MatrixSize);
62 
63     cout<<"矩阵大小 "<<MatrixSize;
64     cout<<"
朴素矩阵算法: "<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)<<" Clocks.."<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)/CLOCKS_PER_SEC<<" Sec";
65     cout<<"
Strassen算法:"<<(endTime_For_Strassen - startTime_For_Strassen)<<" Clocks.."<<(endTime_For_Strassen - startTime_For_Strassen)/CLOCKS_PER_SEC<<" Sec
";
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

相关阅读:
Fluent 18.0新功能之：其他
 【小白的CFD之旅】小结及预告
 【小白的CFD之旅】19 来自计算网格的困惑
 【小白的CFD之旅】18 控制方程基础
 【小白的CFD之旅】23 串行与并行
 【小白的CFD之旅】22 好网格与坏网格
 JS ES6的变量的结构赋值
 JS中some()和every()和join()和concat()和pop(),push(),shift(),unshfit()和map()和filter()
JS（原生js和jq方式）获取元素属性（自定义属性），删除属性（自定义属性）
javascript中this的指向问题
原文地址：https://www.cnblogs.com/zhoutaotao/p/3963048.html

算法导论-矩阵乘法-strassen算法

目录

1、矩阵相乘的朴素算法

2、矩阵相乘的strassen算法

3、完整测试代码c++

4、性能分析

5、参考资料

内容

1、矩阵相乘的朴素算法 T(n) = Θ(n3)

2、矩阵相乘的strassen算法 T(n)=Θ(n^log7) =Θ (n^2.81)

3、完成测试代码

4、性能分析

可以发现：可以看到使用Strassen算法时，耗时不但没有减少，反而剧烈增多，在n=512时计算时间就无法忍受，效果没有朴素矩阵算法好。网上查阅资料，现罗列如下：

5、参考资料

算法导论-矩阵乘法-strassen算法

目录

1、矩阵相乘的朴素算法

2、矩阵相乘的strassen算法

3、完整测试代码c++

4、性能分析

5、参考资料

内容

1、矩阵相乘的朴素算法 T(n) = Θ(n3)

2、矩阵相乘的strassen算法 T(n)=Θ(nlog7) =Θ (n2.81)

3、完成测试代码

4、性能分析

可以发现：可以看到使用Strassen算法时，耗时不但没有减少，反而剧烈增多，在n=512时计算时间就无法忍受，效果没有朴素矩阵算法好。网上查阅资料，现罗列如下：

5、参考资料

2、矩阵相乘的strassen算法 T(n)=Θ(n^log7) =Θ (n^2.81)