樸素矩陣相乘算法,思想明瞭,編程實現簡單。時間複雜度是Θ(n^3)。僞碼以下html
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]
矩陣乘法中採用分治法,第一感受上應該可以有效的提升算法的效率。以下圖所示分治法方案,以及對該算法的效率分析。有圖可知,算法效率是Θ(n^3)。算法效率並無提升。下面介紹下矩陣分治法思想:ios
鑑於上面的分治法方案沒法有效提升算法的效率,要想提升算法效率,由主定理方法可知必須想辦法將2中遞歸式中的係數8減小。Strassen提出了一種將係數減小到7的分治法方案,以下圖所示。c++
效率分析以下:算法
僞碼以下:編程
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];
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.cpp ide
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 }
輸出:性能
矩陣大小 | 樸素矩陣算法(秒) | 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 |
1)採用Strassen算法做遞歸運算,須要建立大量的動態二維數組,其中分配堆內存空間將佔用大量計算時間,從而掩蓋了Strassen算法的優點測試
2)因而對Strassen算法作出改進,設定一個界限。當n<界限時,使用普通法計算矩陣,而不繼續分治遞歸。須要合理設置界限,不一樣環境(硬件配置)下界限不一樣this
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算法。
【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