Rešitev: Strassenov algoritem za množenje matrik (Mathematica)

Iz MaFiRaWiki

Naloga: Strassenov algoritem za množenje matrik

Uporabili bomo paket, ki omogoča lažje operiranje z matrikami in je vgrajen v Mathematici:

<< LinearAlgebra`MatrixManipulation`

1. Problem rekurzivno delimo do konca (do matrik dimenzije 1 x 1)

 Clear[strassenovoMnozenje]
 strassenovoMnozenje[{a_},{b_}]:={a*b}
 strassenovoMnozenje[A_List,B_List]/;Length[A]==Length[Transpose[A]]==Length[B]==Length[Transpose[B]] && EvenQ[Length[A]]:=
  Module[{n,A11,A12,A21,A22,B11,B12,B21,B22,P,Q,R,S,T,U,V,C11,C12,C21,C22,C},
    n=Length[A];
    A11=SubMatrix[A, {1,1},{n/2,n/2}];
    A12=SubMatrix[A, {1,n/2+1},{n/2,n/2}];
    A21=SubMatrix[A, {n/2+1,1}, {n/2,n/2}];
    A22=SubMatrix[A, {n/2+1,n/2+1}, {n/2,n/2}];
    B11=SubMatrix[B, {1,1}, {n/2,n/2}];
    B12=SubMatrix[B, {1,n/2+1}, {n/2,n/2}];
    B21=SubMatrix[B, {n/2+1,1}, {n/2,n/2}];
    B22=SubMatrix[B, {n/2+1,n/2+1}, {n/2,n/2}];
    P=strassenovoMnozenje[A11+A22, B11+B22];
    Q=strassenovoMnozenje[A21+A22, B11];
    R=strassenovoMnozenje[A11, B12-B22];
    S=strassenovoMnozenje[A22, B21-B11];
    T=strassenovoMnozenje[A11+A12, B22];
    U=strassenovoMnozenje[A21-A11, B11+B12];
    V=strassenovoMnozenje[A12-A22, B21+B22];
    C11=P+S-T+V;
    C12=R+T;
    C21=Q+S;
    C22=P+R-Q+U;
    C=BlockMatrix[{{C11,C12}, {C21,C22}}];
    Return[C]
   ]
 strassenovoMnozenje[A_List,B_List]/;Length[A]==Length[Transpose[A]]==Length[B]==Length[Transpose[B]] && OddQ[Length[A]]:=
  Print["Matrik ni mogoče zmnožiti zgolj z uporabo Strassnovega algoritma, saj nista dimenzije n x n za n=2^k za neko nenegativno 
   celo število k."];
 strassenovoMnozenje[A_List,B_List]:=Print["Ena od matrik ni kvadratna ali pa nista enakih dimenzij."]
 strassenovoMnozenje::usage="Funkcija strassenovoMnozenje[A,B] vrne produkt matrik A*B izračunan po Strassenovem algoritmu, če sta 
  obe matriki kvadratni in velja dim(A)=dim(B)=2^k za neko nenegativno celo število k. Če sta matriki kvadratni in enako veliki, 
  vendar njuna dimenzija ni enaka 2^k za neko nenegativno celo število k, potem Strassenov algoritem ne deluje in to funkcija tudi 
  izpiše. Če pa ena od matrik ni kvadratna ali pa nista enakih dimanzij, funkcija to tudi izpiše. 
  Opomba: Ta funkcija PROBLEM REKURZIVNO DELI DO KONCA (do matrik dimenzije 1x1).";

Preizkusimo na primerih:

 Clear[A1,B1]
 A1={{1,2,3}, {5,6,7}, {9,10,11}};
 B1={{17,18,19}, {21,22,23}, {25,26,27}};
 In[11]:= strassenovoMnozenje[A1,B1]
 Out[11]:= "Matrik ni mogoče zmnožiti zgolj z uporabo Strassnovega algoritma, saj nista dimenzije n x n za n=2^k za neko nenegativno  
  celo število k."

 Clear[A2,B2]
 A2={{1,2,3,4}, {5,6,7,8}, {9,10,11,12}, {13,14,15,16}};
 B2={{17,18,19,20}, {21,22,23,24}, {25,26,27,28}, {29,30,31,32}};
 In[11]:= strassenovoMnozenje[A2,B2]
 Out[11]:= {{250,260,270,280},{618,644,670,696},{986,1028,1070,1112},{1354,1412,1470,1528}}

Oglejmo si še čas, ki ga potrebujemo, da zmnožimo matrike velikih dimenzij:

 Clear[A128]
 A128=Table[Range[128], {128}];
 In[18]:= Timing[strassenovoMnozenje[A128,A128];]
 Out[18]:= {82.984 Second, Null}

 Clear[A256]
 A256=Table[Range[256], {256}];
 In[21]:= Timing[strassenovoMnozenje[A256,A256];]
 Out[21]:= {554.437 Second,Null}

2. Iskanje velikosti matrik pri katerih je (časovno) najboljše, da rekurzijo ustavimo in algoritem, ki preneha z rekurzijo že pri tej velikosti

Preizkušamo koliko časa potrebujemo, da izračunamo produkt matrik, če v funkciji spreminjamo velikost matrik pri katerih rekurzijo Strassnovega algoritma ustavimo. Izmerjeni časi (v sekundah) so podani v spodnji tabeli.


velikost matrik ob ustavitvi rekurzije \ velikost matrik, na katerih preizkusimo funkcije 128 256 521 1024 2048
rekurzivno delimo do konca (strassnovoMnozenje) 82,984 554, 437 / / /
15 0, 391 2, 047 14, 64 (veliko) /
20 0, 094 0, 437 3, 188 23, 152 /
35 0, 031 0, 156 1, 235 9, 078 78, 688
100 0, 016 0, 109 0, 813 6, 359 65, 594
200 0, 015 0, 094 0, 766 5, 797 72, 328
300 0, 016 0, 156 1, 11 8, 125 (veliko)
brez uporabe Strassnovega algoritma (A.B vgrajena funkcija) 0, 016 0, 156 1, 343 50, 141 /


Vidimo, da čas računanja za matrike velikosti 128, 256, 521 in 1024 pada vse do funkcije, kjer rekurzijo ustavimo pri matrikah velikosti 200. Torej bi bilo za te matrike rekurzijo najbolje ustaviti pri 200. Vendar pa iz tabele vidimo, da je za matrike velikosti 2048 rekurzijo najbolje ustaviti pri matrikah velikosti 100. Ker je Strassnov algoritem uporaben za matrike velikih dimenzij (to vidimo tudi iz tabele, saj so občutne razlike med običajnim množenjem, ki ga ima vgrajenega Mathematica, in Strassnovim množenjem občutne šele pri matrikah velikosti 1024 in 2048), je najbolje, da rekurzijo ustavimo pri matrikah velikosti 100.

Algoritem, ki rekurzijo ustavi že pri matrikah velikosti 100:

 Clear[strassenovoMnozenje100]
 strassenovoMnozenje100[A_List,B_List]/;Length[A]==Length[Transpose[A]]==Length[B]==Length[Transpose[B]]≤100:=A.B
 strassenovoMnozenje100[A_List,B_List]/;Length[A]==Length[Transpose[A]]==Length[B]==Length[Transpose[B]]>100 && EvenQ[Length[A]]:=
  Module[{n,A11,A12,A21,A22,B11,B12,B21,B22,P,Q,R,S,T,U,V,C11,C12,C21,C22,C},
    n=Length[A];
    A11=SubMatrix[A, {1,1},{n/2,n/2}];
    A12=SubMatrix[A, {1,n/2+1},{n/2,n/2}];
    A21=SubMatrix[A, {n/2+1,1}, {n/2,n/2}];
    A22=SubMatrix[A, {n/2+1,n/2+1}, {n/2,n/2}];
    B11=SubMatrix[B, {1,1}, {n/2,n/2}];
    B12=SubMatrix[B, {1,n/2+1}, {n/2,n/2}];
    B21=SubMatrix[B, {n/2+1,1}, {n/2,n/2}];
    B22=SubMatrix[B, {n/2+1,n/2+1}, {n/2,n/2}];
    P=strassenovoMnozenje100[A11+A22, B11+B22];
    Q=strassenovoMnozenje100[A21+A22, B11];
    R=strassenovoMnozenje100[A11, B12-B22];
    S=strassenovoMnozenje100[A22, B21-B11];
    T=strassenovoMnozenje100[A11+A12, B22];
    U=strassenovoMnozenje100[A21-A11, B11+B12];
    V=strassenovoMnozenje100[A12-A22, B21+B22];
    C11=P+S-T+V;
    C12=R+T;
    C21=Q+S;
    C22=P+R-Q+U;
    C=BlockMatrix[{{C11,C12}, {C21,C22}}];
    Return[C]
   ]
 strassenovoMnozenje100[A_List,B_List]/;Length[A]==Length[Transpose[A]]==Length[B]==Length[Transpose[B]] && OddQ[Length[A]]:=
  Print["Matrik ni mogoče zmnožiti zgolj z uporabo Strassnovega algoritma, saj nista dimenzije n x n za n=2^k za neko nenegativno 
   celo število k."];
 strassenovoMnozenje100[A_List,B_List]:=Print["Ena od matrik ni kvadratna ali pa nista enakih dimenzij."]
 strassenovoMnozenje100::usage="Funkcija strassenovoMnozenje100[A,B] vrne produkt matrik A*B izračunan po Strassenovem algoritmu, če 
  sta obe matriki kvadratni in velja dim(A)=dim(B)=2^k za neko nenegativno celo število k. Od funkcije strassenovoMnozenje se 
  razlikuje v tem, da se tu REKURZIJA USTAVI PRI MATRIKAH VELIKOSTI 100 (funkcija strassenovoMnozenje pa problem rekurzivno deli do 
  konca, do matrik dimenzije 1x1). Če sta matriki kvadratni in enako veliki, vendar njuna dimenzija ni enaka 2^k za neko nenegativno 
  celo število k, potem Strassenov algoritem ne deluje in to funkcija tudi izpiše. Če pa ena od matrik ni kvadratna ali pa nista 
  enakih dimanzij, funkcija to tudi izpiše.";

Glej tudi

Osebna orodja