目录
1、矩阵相乘的朴素算法
2、矩阵相乘的strassen算法
3、完整测试代码c++
4、性能分析
5、参考资料
内容
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)=Θ(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 //递归求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]<<"\t";
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<<"\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 test
43 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 test
53 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算法。
分析原因:(网上总结的说法)
仔细研究后发现,采用Strassen算法作递归运算,需要创建大量的动态二维数组,其中分配堆内存空间将占用大量计算时间,从而掩盖了Strassen算法的优势。于是对Strassen算法做出改进,设定一个界限。当n<界限时,使用普通法计算矩阵,而不继续分治递归。
改进后算法优势明显,就算时间大幅下降。之后,针对不同大小的界限进行试验。在初步试验中发现,当数据规模小于1000时,下界S法的差别不大,规模大于1000以后,n取值越大,消耗时间下降。最优的界限值在32~128之间。
因为计算机每次运算时的系统环境不同(CPU占用、内存占用等),所以计算出的时间会有一定浮动。虽然这样,试验结果已经能得出结论Strassen算法比常规法优势明显。使用下界法改进后,在分治效率和动态分配内存间取舍,针对不同的数据规模稍加试验可以得到一个最优的界限。
时间复杂度就马上降下来了。。但是不要过于乐观。
从实用的观点看,Strassen算法通常不是矩阵乘法所选择的方法:
1 在Strassen算法的运行时间中,隐含的常数因子比简单的O(n^3)方法常数因子大
2 当矩阵是稀疏的时候,为稀疏矩阵设计的算法更快
3 Strassen算法不像简单方法那样子具有数值稳定性
4 在递归层次中生成的子矩阵要消耗空间。
所以矩阵乘法一般意义上还是选择的是朴素的方法,只有当矩阵变稠密,而且矩阵的阶数>20左右,才会考虑使用Strassen算法。
5、参考资料
【4】http://www.xuebuyuan.com/552410.html