Algorytm Strassena


Algorytm Strassena w encyklopedii

Z Wikipedii, wolnej encyklopedii Przejdź do nawigacji Przejdź do wyszukiwania

Algorytm Strassenaalgorytm wykorzystywany do mnożenia macierzy. Został nazwany na cześć swego twórcy – Volkera Strassena. Jest on szybszy od standardowego mnożenia macierzy i jest użyteczny dla dużych macierzy, ale działa wolniej od najszybszej znanej obecnie metody mnożenia – algorytmu Coppersmitha-Winograda – dla bardzo dużych macierzy.

Algorytm został opublikowany w 1969 roku. Pomimo że był tylko trochę szybszy od metody klasycznej, jako pierwszy wykazał, że standardowe podejście nie jest optymalne. Dzięki niemu rozpoczęto poszukiwania jeszcze szybszych algorytmów, aż do algorytmu Coppersmitha-Winograda opublikowanego w 1987 roku.

Spis treści

Działanie algorytmu | edytuj kod

Lewa kolumna przedstawia mnożenie macierzy o rozmiarze 2 na 2. Naiwna metoda mnożenia wymaga jednego mnożenia dla każdej jedynki z lewej kolumny. Pozostałe kolumny reprezentują kolejne z siedmiu mnożeń algorytmu. Suma kolumn daje prawidłowe pomnożenie macierzy pokazane w lewej kolumnie.

Niech A {\displaystyle A} i B {\displaystyle B} będą dwiema macierzami kwadratowymi nad pierścieniem R . {\displaystyle R.} Chcemy obliczyć macierz C {\displaystyle C} taką, że

C = A B A , B , C R 2 n × 2 n {\displaystyle \mathbf {C} =\mathbf {A} \mathbf {B} \qquad \mathbf {A} ,\mathbf {B} ,\mathbf {C} \in R^{2^{n}\times 2^{n}}}

Jeżeli macierze A {\displaystyle A} i B {\displaystyle B} nie są rozmiaru 2 n × 2 n {\displaystyle 2^{n}\times 2^{n}} to brakujące kolumny i wiersze wypełniamy zerami. Następnie dzielimy macierze A , {\displaystyle A,} B {\displaystyle B} i C {\displaystyle C} na macierze klatkowe o równych rozmiarach.

A = A 1 , 1 A 1 , 2 A 2 , 1 A 2 , 2 ,   B = B 1 , 1 B 1 , 2 B 2 , 1 B 2 , 2 ,   C = C 1 , 1 C 1 , 2 C 2 , 1 C 2 , 2 , {\displaystyle \mathbf {A} ={\begin{bmatrix}\mathbf {A} _{1,1}&\mathbf {A} _{1,2}\\\mathbf {A} _{2,1}&\mathbf {A} _{2,2}\end{bmatrix}},\ \mathbf {B} ={\begin{bmatrix}\mathbf {B} _{1,1}&\mathbf {B} _{1,2}\\\mathbf {B} _{2,1}&\mathbf {B} _{2,2}\end{bmatrix}},\ \mathbf {C} ={\begin{bmatrix}\mathbf {C} _{1,1}&\mathbf {C} _{1,2}\\\mathbf {C} _{2,1}&\mathbf {C} _{2,2}\end{bmatrix}},}

gdzie:

A i , j , B i , j , C i , j R 2 n 1 × 2 n 1 , {\displaystyle \mathbf {A} _{i,j},\mathbf {B} _{i,j},\mathbf {C} _{i,j}\in R^{2^{n-1}\times 2^{n-1}},}

czyli

C 1 , 1 = A 1 , 1 B 1 , 1 + A 1 , 2 B 2 , 1 , {\displaystyle \mathbf {C} _{1,1}=\mathbf {A} _{1,1}\mathbf {B} _{1,1}+\mathbf {A} _{1,2}\mathbf {B} _{2,1},} C 1 , 2 = A 1 , 1 B 1 , 2 + A 1 , 2 B 2 , 2 , {\displaystyle \mathbf {C} _{1,2}=\mathbf {A} _{1,1}\mathbf {B} _{1,2}+\mathbf {A} _{1,2}\mathbf {B} _{2,2},} C 2 , 1 = A 2 , 1 B 1 , 1 + A 2 , 2 B 2 , 1 , {\displaystyle \mathbf {C} _{2,1}=\mathbf {A} _{2,1}\mathbf {B} _{1,1}+\mathbf {A} _{2,2}\mathbf {B} _{2,1},} C 2 , 2 = A 2 , 1 B 1 , 2 + A 2 , 2 B 2 , 2 . {\displaystyle \mathbf {C} _{2,2}=\mathbf {A} _{2,1}\mathbf {B} _{1,2}+\mathbf {A} _{2,2}\mathbf {B} _{2,2}.}

W powyższym przypadku nie zredukowaliśmy liczby potrzebnych do wykonania mnożeń. Wciąż musimy wykonać 8 mnożeń, czyli tyle samo, jak w przypadku zwykłego mnożenia macierzy.

Definiujemy nowe macierze według poniższego wzoru:

M 1 := ( A 1 , 1 + A 2 , 2 ) ( B 1 , 1 + B 2 , 2 ) , {\displaystyle \mathbf {M} _{1}:=(\mathbf {A} _{1,1}+\mathbf {A} _{2,2})(\mathbf {B} _{1,1}+\mathbf {B} _{2,2}),} M 2 := ( A 2 , 1 + A 2 , 2 ) B 1 , 1 , {\displaystyle \mathbf {M} _{2}:=(\mathbf {A} _{2,1}+\mathbf {A} _{2,2})\mathbf {B} _{1,1},} M 3 := A 1 , 1 ( B 1 , 2 B 2 , 2 ) , {\displaystyle \mathbf {M} _{3}:=\mathbf {A} _{1,1}(\mathbf {B} _{1,2}-\mathbf {B} _{2,2}),} M 4 := A 2 , 2 ( B 2 , 1 B 1 , 1 ) , {\displaystyle \mathbf {M} _{4}:=\mathbf {A} _{2,2}(\mathbf {B} _{2,1}-\mathbf {B} _{1,1}),} M 5 := ( A 1 , 1 + A 1 , 2 ) B 2 , 2 , {\displaystyle \mathbf {M} _{5}:=(\mathbf {A} _{1,1}+\mathbf {A} _{1,2})\mathbf {B} _{2,2},} M 6 := ( A 2 , 1 A 1 , 1 ) ( B 1 , 1 + B 1 , 2 ) , {\displaystyle \mathbf {M} _{6}:=(\mathbf {A} _{2,1}-\mathbf {A} _{1,1})(\mathbf {B} _{1,1}+\mathbf {B} _{1,2}),} M 7 := ( A 1 , 2 A 2 , 2 ) ( B 2 , 1 + B 2 , 2 ) . {\displaystyle \mathbf {M} _{7}:=(\mathbf {A} _{1,2}-\mathbf {A} _{2,2})(\mathbf {B} _{2,1}+\mathbf {B} _{2,2}).}

Powyższe macierze zostaną użyte do wyrażenia macierzy C i , j {\displaystyle C_{i,j}} w zależności od Mk. Dzięki takiemu zdefiniowaniu macierzy możemy wyeliminować jedno z mnożeń w algorytmie, wyrażając C i , j {\displaystyle C_{i,j}} jako:

C 1 , 1 = M 1 + M 4 M 5 + M 7 , {\displaystyle \mathbf {C} _{1,1}=\mathbf {M} _{1}+\mathbf {M} _{4}-\mathbf {M} _{5}+\mathbf {M} _{7},} C 1 , 2 = M 3 + M 5 , {\displaystyle \mathbf {C} _{1,2}=\mathbf {M} _{3}+\mathbf {M} _{5},} C 2 , 1 = M 2 + M 4 , {\displaystyle \mathbf {C} _{2,1}=\mathbf {M} _{2}+\mathbf {M} _{4},} C 2 , 2 = M 1 M 2 + M 3 + M 6 . {\displaystyle \mathbf {C} _{2,2}=\mathbf {M} _{1}-\mathbf {M} _{2}+\mathbf {M} _{3}+\mathbf {M} _{6}.}

Iterujemy powyższy proces podziału n {\displaystyle n} razy, dopóki podmacierze nie przerodzą się w liczby (elementy pierścienia R {\displaystyle R} ).

Praktyczna implementacja algorytmu Strassena przekształca standardową metodę mnożenia macierzy do wystarczająco małych podmacierzy, dla których algorytm jest bardziej efektywny. Szczególny przypadek, dla którego algorytm Strassena jest bardziej efektywny, zależy od szczegółów jego implementacji oraz sprzętu.

Złożoność obliczeniowa | edytuj kod

Standardowe mnożenie macierzy posiada złożoność obliczeniową rzędu θ ( n 3 ) , {\displaystyle \theta (n^{3}),} gdyż potrzebnych jest 8 rekurencyjnych operacji mnożenia macierzy. W metodzie Strassena dzięki zastosowaniu innego rekurencyjnego podejścia do problemu wymagane jest 7 rekurencyjnych mnożeń oraz θ ( n 2 ) {\displaystyle \theta (n^{2})} operacji dodawania lub odejmowania wartości skalarnej. Wyznaczając złożoność obliczeniową, otrzymujemy wynik θ ( n lg 7 ) , {\displaystyle \theta (n^{\lg 7}),} czyli w przybliżeniu θ ( n 2 , 81 ) . {\displaystyle \theta (n^{2{,}81}).}

Implementacja algorytmu w języku C | edytuj kod

/*------------------------------------------------------------------------------*/ /* Program kompilować bez konsolidacji */ /* Zakładając, że nazwa pliku to Strassen.c to przy użyciu gcc kompilacja powinna wyglądać: */ /* gcc -c Strassen.c */ #include <stdio.h> #include <stdlib.h> void strassen(double **a, double **b, double **c, int tam); void sum(double **a, double **b, double **result, int tam); void subtract(double **a, double **b, double **result, int tam); double **allocate_real_matrix(int tam, int random); double **free_real_matrix(double **v, int tam); void strassen(double **a, double **b, double **c, int tam) { // najprostszy przypadek: kiedy macierz ma rozmiar 1 na 1: if (tam == 1) { c[0][0] = a[0][0] * b[0][0]; return; } // pozostałe przypadki: else { int newTam = tam/2; double **a11, **a12, **a21, **a22; double **b11, **b12, **b21, **b22; double **c11, **c12, **c21, **c22; double **p1, **p2, **p3, **p4, **p5, **p6, **p7; // alokacja pamięci: a11 = allocate_real_matrix(newTam, -1); a12 = allocate_real_matrix(newTam, -1); a21 = allocate_real_matrix(newTam, -1); a22 = allocate_real_matrix(newTam, -1); b11 = allocate_real_matrix(newTam, -1); b12 = allocate_real_matrix(newTam, -1); b21 = allocate_real_matrix(newTam, -1); b22 = allocate_real_matrix(newTam, -1); c11 = allocate_real_matrix(newTam, -1); c12 = allocate_real_matrix(newTam, -1); c21 = allocate_real_matrix(newTam, -1); c22 = allocate_real_matrix(newTam, -1); p1 = allocate_real_matrix(newTam, -1); p2 = allocate_real_matrix(newTam, -1); p3 = allocate_real_matrix(newTam, -1); p4 = allocate_real_matrix(newTam, -1); p5 = allocate_real_matrix(newTam, -1); p6 = allocate_real_matrix(newTam, -1); p7 = allocate_real_matrix(newTam, -1); double **aResult = allocate_real_matrix(newTam, -1); double **bResult = allocate_real_matrix(newTam, -1); int i, j; //dzielenie macierzy na cztery podmacierze: for (i = 0; i < newTam; i++) { for (j = 0; j < newTam; j++) { a11[i][j] = a[i][j]; a12[i][j] = a[i][j + newTam]; a21[i][j] = a[i + newTam][j]; a22[i][j] = a[i + newTam][j + newTam]; b11[i][j] = b[i][j]; b12[i][j] = b[i][j + newTam]; b21[i][j] = b[i + newTam][j]; b22[i][j] = b[i + newTam][j + newTam]; } } //obliczanie macierzy od M1 do M7: sum(a11, a22, aResult, newTam); // a11 + a22 sum(b11, b22, bResult, newTam); // b11 + b22 strassen(aResult, bResult, p1, newTam); // p1 = (a11+a22) * (b11+b22) sum(a21, a22, aResult, newTam); // a21 + a22 strassen(aResult, b11, p2, newTam); // p2 = (a21+a22) * (b11) subtract(b12, b22, bResult, newTam); // b12 - b22 strassen(a11, bResult, p3, newTam); // p3 = (a11) * (b12 - b22) subtract(b21, b11, bResult, newTam); // b21 - b11 strassen(a22, bResult, p4, newTam); // p4 = (a22) * (b21 - b11) sum(a11, a12, aResult, newTam); // a11 + a12 strassen(aResult, b22, p5, newTam); // p5 = (a11+a12) * (b22) subtract(a21, a11, aResult, newTam); // a21 - a11 sum(b11, b12, bResult, newTam); // b11 + b12 strassen(aResult, bResult, p6, newTam); // p6 = (a21-a11) * (b11+b12) subtract(a12, a22, aResult, newTam); // a12 - a22 sum(b21, b22, bResult, newTam); // b21 + b22 strassen(aResult, bResult, p7, newTam); // p7 = (a12-a22) * (b21+b22) //obliczanie c21, c21, c11 i c22: sum(p3, p5, c12, newTam); // c12 = p3 + p5 sum(p2, p4, c21, newTam); // c21 = p2 + p4 sum(p1, p4, aResult, newTam); // p1 + p4 sum(aResult, p7, bResult, newTam); // p1 + p4 + p7 subtract(bResult, p5, c11, newTam); // c11 = p1 + p4 - p5 + p7 sum(p1, p3, aResult, newTam); // p1 + p3 sum(aResult, p6, bResult, newTam); // p1 + p3 + p6 subtract(bResult, p2, c22, newTam); // c22 = p1 + p3 - p2 + p6 //grupowanie wyników w jedną macierz: for (i = 0; i < newTam ; i++) { for (j = 0 ; j < newTam ; j++) { c[i][j] = c11[i][j]; c[i][j + newTam] = c12[i][j]; c[i + newTam][j] = c21[i][j]; c[i + newTam][j + newTam] = c22[i][j]; } } //dealokacja pamięci: a11 = free_real_matrix(a11, newTam); a12 = free_real_matrix(a12, newTam); a21 = free_real_matrix(a21, newTam); a22 = free_real_matrix(a22, newTam); b11 = free_real_matrix(b11, newTam); b12 = free_real_matrix(b12, newTam); b21 = free_real_matrix(b21, newTam); b22 = free_real_matrix(b22, newTam); c11 = free_real_matrix(c11, newTam); c12 = free_real_matrix(c12, newTam); c21 = free_real_matrix(c21, newTam); c22 = free_real_matrix(c22, newTam); p1 = free_real_matrix(p1, newTam); p2 = free_real_matrix(p2, newTam); p3 = free_real_matrix(p3, newTam); p4 = free_real_matrix(p4, newTam); p5 = free_real_matrix(p5, newTam); p6 = free_real_matrix(p6, newTam); p7 = free_real_matrix(p7, newTam); aResult = free_real_matrix(aResult, newTam); bResult = free_real_matrix(bResult, newTam); } //koniec przypadku else } //koniec funkcji Strassen /*------------------------------------------------------------------------------*/ //funkcja sumująca dwie macierze void sum(double **a, double **b, double **result, int tam) { int i, j; for (i = 0; i < tam; i++) { for (j = 0; j < tam; j++) { result[i][j] = a[i][j] + b[i][j]; } } } /*------------------------------------------------------------------------------*/ //funkcja odejmująca jedną macierz od drugiej void subtract(double **a, double **b, double **result, int tam) { int i, j; for (i = 0; i < tam; i++) { for (j = 0; j < tam; j++) { result[i][j] = a[i][j] - b[i][j]; } } } /*------------------------------------------------------------------------------*/ // Ta funkcja alokuje macierz przy użyciu funkcji malloc oraz inicjalizuje ją. Jeżeli zmienna losowa // jest zerem to funkcja inicjalizuje macierz o wartościach zerowych, jeżeli jest jedynką to inicjalizuje // inicjalizuje macierz z losowymi wartościami. Jeżeli jest jakąkolwiek inną wartością to macierz jest // inicjalizowana bez żadnych wartości Zmienna tam określa długość macierzy. double **allocate_real_matrix(int tam, int random) { int i, j, n = tam, m = tam; double **v, a; // wskaźnik na wektor // alokuje jeden wektor wektoru (macierz) v = (double**) malloc(n * sizeof(double*)); if (v == NULL) { printf ("** Blad alokacji macierzy: Niewystarczajaca pamiec **"); return (NULL); } // alokacja każdego rzędu macierzy for (i = 0; i < n; i++) { v[i] = (double*) malloc(m * sizeof(double)); if (v[i] == NULL) { printf ("** Blad: Niewystarczajaca pamiec **"); free_real_matrix(v, n); return (NULL); } // inicjalizacja macierzy o wartościach zerowych if (random == 0) { for (j = 0; j < m; j++) v[i][j] = 0.0; } // inicjalizacja macierzy o wartościach losowych z przedziału od 0 do 10 else { if (random == 1) { for (j = 0; j < m; j++) { a = rand(); v[i][j] = (a - (int)a) * 10; } } } } return (v); // zwraca wskaźnik na wektor. } /*------------------------------------------------------------------------------*/ // Ta funkcja dealokuje macierz (zwalnia pamięć) double **free_real_matrix(double **v, int tam) { int i; if (v == NULL) { return (NULL); } for (i = 0; i < tam; i++) { if (v[i]) { free(v[i]); // zwalnia rząd macierzy v[i] = NULL; } } free(v); // zwalnia wskaźnik / v = NULL; return (NULL); //zwraca wskaźnik NULL / } /*------------------------------------------------------------------------------*/ 

Bibliografia | edytuj kod

  • Volker Strassen, Gaussian Elimination is not Optimal, Numer. Math. 13, s. 354–356, 1969.
  • Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, Clifford Stein. Wprowadzenie do algorytmów, wyd. siódme, Wydawnictwa Naukowo-Techniczne, s. 751–758, Warszawa 2005.
Na podstawie artykułu: "Algorytm Strassena" pochodzącego z Wikipedii
OryginałEdytujHistoria i autorzy