#算法學習錄#Strassen矩陣乘法

今天我們談談一個“土豪”算法——Strasen矩陣算法
之說以說它“土豪”就是因為其帶來了巨大的空間開銷刊殉。
先來考察一個問題:請用三次實數(shù)乘法計算復數(shù)a+bi和c+di相乘选泻。
由于:
a×(c+d)=ac+ad=s1 焰情;
b×(c-d)=bc-bd=s2 ;
d×(a+b)=ad+bd=s3 错维;
故有實部:s1 -s3 =ac-bd,
虛部:s2+ s3 =ad+bc敷待。
這樣,四次的乘法就變成三次乘法和悦。

Strassen矩陣乘法也是如此退疫,把A,B,C矩陣分解為n/2×n/2子矩陣,進行7次遞歸計算n/2×n/2矩陣的乘法鸽素,其運行時間的遞歸式:

T(n)= Θ(1)???????????? if n=1;

??????7T(n/2)+Θ(n^2 )??? if n>1;

令:
S1=B12-B22;
S2=A11+A12;
S3=A21+A22;
S4=B21-B22;
S5=A11+A22;
S6=B11+B22;
S7=A12-A22;
S8=B21+B22;
S9=A11-A21;
S10=B11+B12;
那么:?P1= A11·S1 = A11·(B12-B22)
P2= B22·S2 = B22·(A11+A12)
P3= B11·S3 = B11·(A21+A22)
P4= A22·S4 = A22·(B21-B22)
P5= S5·S6 = (A11+A22)·(B11+B22)
P6= S7·S8 = (A12-A22)·(B21+B22)
P7= S9·S10 = (A11-A21)·(B11+B12)

C11= P5 + P4 - P2 + P6=A11×B11+A12×B21
C12= P1 + P2=A11×B12+A12×B22
C21= P3 + P4=A21×B11+A22×B21
C22= P5 + P1 – P3 – P7=A21×B21+A22×B22

Strassen算法的具體實現(xiàn)(C語言):
int Strassen(int **A, int **B, int **Result, int Size){
?if (Size == 1){
??//直接計算C11
??Result[0][0] = A[0][0] * B[0][0];
??return 0;
?}
?int NewSize = Size / 2;
?/*分塊矩陣*/
?int **A11, **A12, **A21, **A22;
?int **B11, **B12, **B21, **B22;
?int **C11, **C12, **C21, **C22;

?int **P1, **P2, **P3, **P4, **P5, **P6, **P7;
?/*存放數(shù)組A褒繁、B(i、j)的臨時變量*/
?int **AResult, **BResult;

?A11 = new int*[NewSize];
?A12 = new int*[NewSize];
?A21 = new int*[NewSize];
?A22 = new int*[NewSize];

?B11 = new int*[NewSize];
?B12 = new int*[NewSize];
?B21 = new int*[NewSize];
?B22 = new int*[NewSize];

?C11 = new int*[NewSize];
?C12 = new int*[NewSize];
?C21 = new int*[NewSize];
?C22 = new int*[NewSize];

?P1 = new int*[NewSize];
?P2 = new int*[NewSize];
?P3 = new int*[NewSize];
?P4 = new int*[NewSize];
?P5 = new int*[NewSize];
?P6 = new int*[NewSize];
?P7 = new int*[NewSize];

?AResult = new int*[NewSize];
?BResult = new int*[NewSize];

?for (int i = 0; i < NewSize; i++)
?{
??A11[i] = new int[NewSize];
??A12[i] = new int[NewSize];
??A21[i] = new int[NewSize];
??A22[i] = new int[NewSize];

??B11[i] = new int[NewSize];
??B12[i] = new int[NewSize];
??B21[i] = new int[NewSize];
??B22[i] = new int[NewSize];

??C11[i] = new int[NewSize];
??C12[i] = new int[NewSize];
??C21[i] = new int[NewSize];
??C22[i] = new int[NewSize];

??P1[i] = new int[NewSize];
??P2[i] = new int[NewSize];
??P3[i] = new int[NewSize];
??P4[i] = new int[NewSize];
??P5[i] = new int[NewSize];
??P6[i] = new int[NewSize];
??P7[i] = new int[NewSize];

??AResult[i] = new int[NewSize];
??BResult[i] = new int[NewSize];


?}

?//對分塊矩陣賦值
?for (int i = 0; i < NewSize; i++)
?{
??for (int j = 0; j < NewSize; j++)
??{
???A11[i][j] = A[i][j];
???A12[i][j] = A[i][j + NewSize];
???A21[i][j] = A[i + NewSize][j];
???A22[i][j] = A[i + NewSize][j + NewSize];

???B11[i][j] = B[i][j];
???B12[i][j] = B[i][j + NewSize];
???B21[i][j] = B[i + NewSize][j];
???B22[i][j] = B[i + NewSize][j + NewSize];

??}
?}

?//計算P1 = A11*(B12-B22)
?Sub(B12, B22, BResult, NewSize);
?Strassen(A11, BResult, P1, NewSize);

?//計算P2 = (A11+A12)*B22
?Add(A11, A12, AResult, NewSize);
?Strassen(AResult, B22, P2, NewSize);

?//計算P3 = (A21+A22)*B11
?Add(A21, A22, AResult, NewSize);
?Strassen(AResult, B11, P3, NewSize);

?//計算P4 = A22*(B21-B11)
?Sub(B21, B11, BResult, NewSize);
?Strassen(A22, BResult, P4, NewSize);

?//計算P5 = (A11+A22)*(B11+B22)
?Add(A11, A22, AResult, NewSize);
?Add(B11, B22, BResult, NewSize);
?Strassen(AResult, BResult, P5, NewSize);

?//計算P6 = (A12-A22)*(B21+B22)
?Sub(A12, A22, AResult, NewSize);
?Add(B21, B22, BResult, NewSize);
?Strassen(AResult, BResult, P6, NewSize);

?//計算P7 = (A11-A21)*(B11+B12)
?Sub(A11, A21, AResult, NewSize);
?Add(B11, B12, BResult, NewSize);
?Strassen(AResult, BResult, P7, NewSize);

?//計算C11馍忽,C12棒坏,C21,C22
?//C11 = P5 + P4 - P2 + P6;
?Add(P5, P4, AResult, NewSize);
?Sub(AResult, P2, BResult, NewSize);
?Add(BResult, P6, C11, NewSize);

?//C12=P1+P2
?Add(P1, P2, C12, NewSize);

?//C21=P3+P4
?Add(P3, P4, C21, NewSize);

?//C22=P5+P1-P3-P7
?Add(P5, P1, C22, NewSize);
?Sub(C22, P3, C22, NewSize);
?Sub(C22, P7, C22, NewSize);

?//合并C11遭笋,C12俊抵,C21,C22
?for (int i = 0; i < NewSize; i++)
?{
??for (int j = 0; j < NewSize; j++)
??{
???Result[i][j] = C11[i][j];
???Result[i][j + NewSize] = C12[i][j];
???Result[i + NewSize][j] = C21[i][j];
???Result[i + NewSize][j + NewSize] = C22[i][j];
??}
?}

?//刪除數(shù)組坐梯,回收資源
?for (int i = 0; i < NewSize; i++){
??delete[] A11[i]; delete[] A12[i]; delete[] A21[i]; delete[] A22[i];
??delete[] B11[i]; delete[] B12[i]; delete[] B21[i]; delete[] B22[i];
??delete[] C11[i]; delete[] C12[i]; delete[] C21[i]; delete[] C22[i];
??delete[] P1[i]; delete[] P2[i]; delete[] P3[i]; delete[] P4[i]; delete[] P5[i]; delete[] P6[i]; delete[] P7[i];
??delete[] AResult[i]; delete[] BResult[i];
?}
?delete[] A11; delete[] A12; delete[] A21; delete[] A22;
?delete[] B11; delete[] B12; delete[] B21; delete[] B22;
?delete[] C11; delete[] C12; delete[] C21; delete[] C22;
?delete[] P1; delete[] P2; delete[] P3; delete[] P4; delete[] P5; delete[] P6; delete[] P7;
?delete[] AResult; delete[] BResult;
?return 0;
}


//矩陣相加
void Add(int **A, int **B, int **Q, int Size){
?for (int i = 0; i < Size; i++){
??for (int j = 0; j < Size; j++){
???Q[i][j] = A[i][j] + B[i][j];
??}
?}
}

//矩陣相減
void Sub(int**A, int**B, int **Q, int Size){
?for (int i = 0; i < Size; i++){
??for (int j = 0; j < Size; j++){
???Q[i][j] = A[i][j] - B[i][j];
??}
?}
}
演示結(jié)果:
?


與暴力求解相比:
?? for(i=0;i<m;i++)
???? for(j=0;j<m;j++){
???????? C[i][j]=0;
???? ?for(k=0;k<n;k++)
??????????? C[i][j]+=A[i][k]*B[k][j];????????????
}???
其運行時間(n^lg7,2.80<lg7<2.81)比暴力求解(n3)稍快徽诲,但其精度較低(在處理小數(shù)計算時),且消耗了大量存儲空間吵血。

最后附上源代碼:https://github.com/LRC-cheng/Algorithms_Practise.git

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末谎替,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子蹋辅,更是在濱河造成了極大的恐慌钱贯,老刑警劉巖,帶你破解...
    沈念sama閱讀 211,561評論 6 492
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件侦另,死亡現(xiàn)場離奇詭異秩命,居然都是意外死亡,警方通過查閱死者的電腦和手機褒傅,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,218評論 3 385
  • 文/潘曉璐 我一進店門弃锐,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人殿托,你說我怎么就攤上這事霹菊。” “怎么了支竹?”我有些...
    開封第一講書人閱讀 157,162評論 0 348
  • 文/不壞的土叔 我叫張陵旋廷,是天一觀的道長鸠按。 經(jīng)常有香客問我,道長饶碘,這世上最難降的妖魔是什么目尖? 我笑而不...
    開封第一講書人閱讀 56,470評論 1 283
  • 正文 為了忘掉前任,我火速辦了婚禮扎运,結(jié)果婚禮上瑟曲,老公的妹妹穿的比我還像新娘。我一直安慰自己绪囱,他們只是感情好测蹲,可當我...
    茶點故事閱讀 65,550評論 6 385
  • 文/花漫 我一把揭開白布莹捡。 她就那樣靜靜地躺著鬼吵,像睡著了一般。 火紅的嫁衣襯著肌膚如雪篮赢。 梳的紋絲不亂的頭發(fā)上齿椅,一...
    開封第一講書人閱讀 49,806評論 1 290
  • 那天,我揣著相機與錄音启泣,去河邊找鬼涣脚。 笑死,一個胖子當著我的面吹牛寥茫,可吹牛的內(nèi)容都是我干的遣蚀。 我是一名探鬼主播,決...
    沈念sama閱讀 38,951評論 3 407
  • 文/蒼蘭香墨 我猛地睜開眼纱耻,長吁一口氣:“原來是場噩夢啊……” “哼芭梯!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起弄喘,我...
    開封第一講書人閱讀 37,712評論 0 266
  • 序言:老撾萬榮一對情侶失蹤玖喘,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后蘑志,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體累奈,經(jīng)...
    沈念sama閱讀 44,166評論 1 303
  • 正文 獨居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 36,510評論 2 327
  • 正文 我和宋清朗相戀三年急但,在試婚紗的時候發(fā)現(xiàn)自己被綠了澎媒。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點故事閱讀 38,643評論 1 340
  • 序言:一個原本活蹦亂跳的男人離奇死亡波桩,死狀恐怖旱幼,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情突委,我是刑警寧澤柏卤,帶...
    沈念sama閱讀 34,306評論 4 330
  • 正文 年R本政府宣布冬三,位于F島的核電站,受9級特大地震影響缘缚,放射性物質(zhì)發(fā)生泄漏勾笆。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 39,930評論 3 313
  • 文/蒙蒙 一桥滨、第九天 我趴在偏房一處隱蔽的房頂上張望窝爪。 院中可真熱鬧,春花似錦齐媒、人聲如沸蒲每。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,745評論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽邀杏。三九已至,卻和暖如春唬血,著一層夾襖步出監(jiān)牢的瞬間望蜡,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,983評論 1 266
  • 我被黑心中介騙來泰國打工拷恨, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留脖律,地道東北人。 一個月前我還...
    沈念sama閱讀 46,351評論 2 360
  • 正文 我出身青樓腕侄,卻偏偏與公主長得像小泉,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子冕杠,可洞房花燭夜當晚...
    茶點故事閱讀 43,509評論 2 348

推薦閱讀更多精彩內(nèi)容

  • 貪心算法 貪心算法總是作出在當前看來最好的選擇微姊。也就是說貪心算法并不從整體最優(yōu)考慮,它所作出的選擇只是在某種意義上...
    fredal閱讀 9,223評論 3 52
  • 閱讀經(jīng)典——《算法導論》03 矩陣乘法是種極其耗時的運算拌汇。 以C = A ? B為例柒桑,其中A和B都是 n x n ...
    金戈大王閱讀 27,083評論 10 24
  • 背景 一年多以前我在知乎上答了有關(guān)LeetCode的問題, 分享了一些自己做題目的經(jīng)驗。 張土汪:刷leetcod...
    土汪閱讀 12,738評論 0 33
  • C++運算符重載-下篇 本章內(nèi)容:1. 運算符重載的概述2. 重載算術(shù)運算符3. 重載按位運算符和二元邏輯運算符4...
    Haley_2013閱讀 1,435評論 0 49
  • 一早醒來噪舀,頓悟昨晚電話中為何忽生反感: 你自然而然流露出的優(yōu)越感魁淳!縱然在談論失敗可能性時也透著一種先見之明的優(yōu)越感...
    蝎思君閱讀 267評論 1 4