Algorytm mnożenia macierzy - Matrix multiplication algorithm

Ponieważ mnożenie macierzy jest tak centralną operacją w wielu algorytmach numerycznych , zainwestowano wiele pracy w uczynienie algorytmów mnożenia macierzy wydajnymi. Zastosowania mnożenia macierzy w problemach obliczeniowych można znaleźć w wielu dziedzinach, w tym w obliczeniach naukowych i rozpoznawaniu wzorców oraz w pozornie niepowiązanych problemach, takich jak liczenie ścieżek na wykresie . Zaprojektowano wiele różnych algorytmów mnożenia macierzy na różnych typach sprzętu, w tym w systemach równoległych i rozproszonych , gdzie praca obliczeniowa jest rozłożona na wiele procesorów (być może w sieci).

Bezpośrednie zastosowanie matematycznej definicji mnożenia macierzy daje algorytm, który zajmuje czas rzędu n ³ operacji na polach , aby pomnożyć dwie macierze n × n przez to pole ( Θ( n ³ ) w notacji big O ). Lepsze asymptotyczne ograniczenia czasu potrzebnego do mnożenia macierzy są znane od algorytmu Strassena w latach 60., ale nadal nie wiadomo, jaki jest czas optymalny (tj. jaka jest złożoność obliczeniowa mnożenia macierzy ). Od grudnia 2020 r. algorytm mnożenia macierzy o najlepszej asymptotycznej złożoności działa w czasie O( n ^2,3728596 ) , podanym przez Josha Almana i Virginię Vassilevską Williams , jednak algorytm ten jest algorytmem galaktycznym ze względu na duże stałe i nie można go praktycznie zrealizować.

Algorytm iteracyjny

Definicja mnożenia macierzy , że jeśli C = AB dla n x m macierzy A i m x s macierzy B , następnie C, to n x s matrycy z wpisami

${\ Displaystyle c_ {ij} = \ suma _ {k = 1} ^ {m} a_ {ik} b_ {kj}.}$

Na tej podstawie można skonstruować prosty algorytm, który zapętla indeksy i od 1 do n oraz j od 1 do p , obliczając powyższe za pomocą pętli zagnieżdżonej:

Algorytm ten zajmuje czas Θ( nmp ) (w notacji asymptotycznej ). Powszechnym uproszczeniem na potrzeby analizy algorytmów jest założenie, że wszystkie dane wejściowe są macierzami kwadratowymi o rozmiarze n × n , w którym to przypadku czas działania wynosi Θ( n ³ ) , tj. sześcienny o rozmiarze wymiaru.

Zachowanie w pamięci podręcznej

Ilustracja kolejności wierszy i kolumn głównych

Trzy pętle w iteracyjnym mnożeniu macierzy można dowolnie zamieniać między sobą bez wpływu na poprawność lub asymptotyczny czas działania. Jednak kolejność może mieć znaczny wpływ na praktyczną wydajność ze względu na wzorce dostępu do pamięci i wykorzystanie przez algorytm pamięci podręcznej ; to, która kolejność jest najlepsza, zależy również od tego, czy macierze są przechowywane w kolejności wiersz-główny , kolumna-główny, czy połączenie obu.

W szczególności w wyidealizowanym przypadku w pełni asocjacyjnej pamięci podręcznej składającej się z M bajtów i b bajtów na linię pamięci podręcznej (tj.m/bcache), powyższy algorytm jest suboptymalny dla A i B przechowywanych w kolejności wiersz-główny. Kiedy n >m/b, każda iteracja wewnętrznej pętli (jednoczesne przeciągnięcie przez wiersz A i kolumnę B ) powoduje chybienie pamięci podręcznej podczas uzyskiwania dostępu do elementu B . Oznacza to, że w najgorszym przypadku algorytm powoduje Θ( n ³ ) braków w pamięci podręcznej. Począwszy od 2010 r. szybkość pamięci w porównaniu z procesorami jest taka, że ​​w przypadku dużych macierzy dominują błędy w pamięci podręcznej, a nie rzeczywiste obliczenia.

Optymalnym wariantem algorytmu iteracyjnego dla A i B w układzie wiersz-główny jest wersja kafelkowa , w której macierz jest domyślnie podzielona na kwadratowe kafelki o rozmiarze √ M przez √ M :

W wyidealizowanym modelu pamięci podręcznej ten algorytm wiąże się tylko z Θ(n ³/b √ M) braki w pamięci podręcznej; dzielnik b √ M wynosi kilka rzędów wielkości na nowoczesnych maszynach, tak że rzeczywiste obliczenia dominują w czasie działania, a nie w chybieniach pamięci podręcznej.

Algorytm dziel i zwyciężaj

Alternatywą dla algorytmu iteracyjnego jest algorytm dziel i zwyciężaj do mnożenia macierzy. Opiera się to na partycjonowaniu bloków

${\ Displaystyle C = {\ zacząć {pmatrix} C_ {11} i C_ {12} \ \ C_ {21} i C_ {22} \ \ \ koniec {pmatrix}} \, A = {\ zacząć {pmatrix} A_ { 11}&A_{12}\\A_{21}&A_{22}\\\end{pmatrix}},\,B={\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_ {22}\\\koniec{pmatrix}},}$

co działa dla wszystkich macierzy kwadratowych, których wymiary są potęgami dwójki, tj. kształty wynoszą 2 ⁿ × 2 ⁿ dla niektórych n . Produkt macierzy jest teraz

${\ Displaystyle {\ zacząć {pmatrix} C_ {11} i C_ {12} \ \ C_ {21} i C_ {22} \ \ \ koniec {pmatrix}} = {\ zacząć {pmatrix} A_ {11} i A_ {12} \\A_{21}&A_{22}\\\end{pmatrix}}{\begin{pmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\\\end{pmatrix}} ={\begin{pmatrix}A_{11}B_{11}+A_{12}B_{21}&A_{11}B_{12}+A_{12}B_{22}\\A_{21}B_{11 }+A_{22}B_{21}&A_{21}B_{12}+A_{22}B_{22}\\\end{pmatrix}}}$

który składa się z ośmiu mnożeń par podmacierzy, po których następuje krok dodawania. Algorytm dziel i przejęcie oblicza mniejsze mnożenia rekurencyjnie stosując skalarne mnożenia C ₁₁ = ₁₁ b ₁₁ jak pudełka podstawy.

Złożoność tego algorytmu w funkcji n jest określona przez rekurencję

${\ Displaystyle T (1) = \ Theta (1);}$

${\ Displaystyle T (n) = 8 T (n/2) + \ Theta (n ^ {2})}$

uwzględnienie ośmiu rekurencyjnych wywołań macierzy o rozmiarze n /2 i Θ( n ² ) w celu zsumowania czterech par wynikowych macierzy z uwzględnieniem elementów. Zastosowanie twierdzenia głównego dla rekurencji typu dziel i zwyciężaj pokazuje, że ta rekurencja ma rozwiązanie Θ( n ³ ) , takie samo jak algorytm iteracyjny.

Macierze niekwadratowe

Wariant tego algorytmu, który działa dla macierzy o dowolnych kształtach i jest szybszy w praktyce, dzieli macierze na dwie zamiast na cztery podmacierze w następujący sposób. Podział macierzy oznacza teraz podzielenie jej na dwie części o jednakowej wielkości lub możliwie zbliżone do równych rozmiarów w przypadku nieparzystych wymiarów.

Zachowanie w pamięci podręcznej

Współczynnik chybienia pamięci podręcznej w rekurencyjnym mnożeniu macierzy jest taki sam, jak w przypadku kafelkowej wersji iteracyjnej, ale w przeciwieństwie do tego algorytmu, algorytm rekurencyjny nie zwraca uwagi na pamięć podręczną : nie ma parametru dostrajania wymaganego do uzyskania optymalnej wydajności pamięci podręcznej i zachowuje się dobrze w środowisko wieloprogramowe, w którym rozmiary pamięci podręcznej są efektywnie dynamiczne ze względu na inne procesy zajmujące miejsce w pamięci podręcznej. (Prosty algorytm iteracyjny również nie uwzględnia pamięci podręcznej, ale w praktyce jest znacznie wolniejszy, jeśli układ macierzy nie jest dostosowany do algorytmu.)

Liczba chybień pamięci podręcznej spowodowanych przez ten algorytm na maszynie z M wierszami idealnej pamięci podręcznej, każdy o rozmiarze b bajtów, jest ograniczona przez

${\ Displaystyle \ Theta \ lewo (m + n + p + {\ Frac {mn + np + MP} {b}} + {\ Frac {mnp} {b {\ sqrt {M}}}} \ prawej)}$

Algorytmy subsześcienne

Poprawa estymacji wykładnika ω w czasie dla złożoności obliczeniowej mnożenia macierzy .${\ Displaystyle O (n ^ {\ omega})}$

Istnieją algorytmy, które zapewniają lepsze czasy działania niż te proste. Jako pierwszy odkryto algorytm Strassena , opracowany przez Volkera Strassena w 1969 roku i często określany jako „szybkie mnożenie macierzy”. Opiera się na sposobie mnożenia dwóch macierzy 2 × 2, który wymaga tylko 7 mnożenia (zamiast 8), kosztem kilku dodatkowych operacji dodawania i odejmowania. Zastosowanie tego rekurencyjnie daje algorytm z mnożnikiem kosztu . Algorytm Strassena jest bardziej złożony, a stabilność numeryczna jest zmniejszona w porównaniu z algorytmem naiwnym, ale jest szybszy w przypadkach, gdy n > 100 lub więcej i pojawia się w kilku bibliotekach, takich jak BLAS . Jest to bardzo przydatne w przypadku dużych macierzy nad dokładnymi domenami, takimi jak pola skończone , gdzie stabilność numeryczna nie jest problemem. ${\ Displaystyle O (n ^ {\ log _ {2} 7}) \ około O (n ^ {2,807})}$

Pytaniem otwartym w informatyce teoretycznej jest to, jak dobrze algorytm Strassena można ulepszyć pod względem asymptotycznej złożoności . Mnożenie macierzy wykładnik , zazwyczaj oznaczone jest najmniejsza liczba rzeczywista, dla której każdy matrycy na pola mogą być mnożone przez siebie za pomocą działania pola. Obecnie najlepszym zawodnikiem jest Josh Alman i Virginia Vassilevska Williams . Algorytm ten, podobnie jak wszystkie inne najnowsze algorytmy w tej dziedzinie badań, jest uogólnieniem algorytmu Coppersmitha-Winograda, który został podany przez Dona Coppersmitha i Shmuela Winograda w 1990 roku. Koncepcja koncepcyjna tych algorytmów jest podobna do algorytmu Strassena: sposób jest przeznaczony do mnożenia dwóch macierzy k × k z mniej niż k ³ mnożeń, a technika ta jest stosowana rekurencyjnie. Jednak stały współczynnik ukryty przez notację Big O jest tak duży, że te algorytmy są warte zachodu tylko w przypadku macierzy, które są zbyt duże do obsługi na współczesnych komputerach. ${\ Displaystyle \ omega}$${\displaystyle n\razy n}$${\ Displaystyle n ^ {\ omega + O (1)}}$${\ Displaystyle \ omega}$${\ Displaystyle \ omega <2.3728596}$

Algorytm Freivaldsa jest prostym algorytmem Monte Carlo, który przy danych macierzach A , B i C weryfikuje w czasie Θ ( n ² ) jeśli AB = C .

Algorytmy równoległe i rozproszone

Równoległość pamięci współdzielonej

Dziel i zwyciężaj algorytm zarysowane wcześniej można parallelized dwojako dla wieloprocesorowych wspólną pamięcią . Opierają się one na fakcie, że osiem rekurencyjnych mnożeń macierzy w

${\ Displaystyle {\ zacząć {pmatrix} A_ {11} B_ {11} + A_ {12} B_ {21} i A_ {11} B_ {12} + A_ {12} B_ {22} \ \ A_ {21} B_ {11}+A_{22}B_{21}&A_{21}B_{12}+A_{22}B_{22}\\\end{pmatrix}}}$

mogą być wykonywane niezależnie od siebie, podobnie jak cztery podsumowania (chociaż algorytm musi „połączyć” mnożenia przed wykonaniem podsumowań). Wykorzystując pełną równoległość problemu, otrzymujemy algorytm, który można wyrazić w pseudokodzie typu fork-join :

W tym przypadku fork jest słowem kluczowym, które sygnalizuje, że obliczenia mogą być wykonywane równolegle z resztą wywołania funkcji, podczas gdy join czeka na zakończenie wszystkich wcześniej "rozwidlonych" obliczeń. partycja osiąga swój cel tylko przez manipulację wskaźnikiem.

Algorytm ten ma krytyczną długość ścieżki w Θ (log ² n ) kroków, co oznacza, że zajmuje to dużo czasu na maszynie idealny z nieskończonej liczby procesorów; W związku z tym, że ma maksymalny możliwy przyspieszenie w Θ ( n ³ / log ² N ) na każdej fizycznego komputera. Algorytm nie jest praktyczny ze względu na koszt komunikacji związany z przenoszeniem danych do iz tymczasowej macierzy T , ale bardziej praktyczny wariant osiąga przyspieszenie Θ( n ² ) bez użycia tymczasowej macierzy.

Mnożenie macierzy bloków. W algorytmie 2D każdy procesor odpowiada za jedną podmacierz C . W algorytmie 3D każda para podmacierzy z A i B, która jest mnożona, jest przypisywana do jednego procesora.

Algorytmy unikające komunikacji i rozproszone

W nowoczesnych architekturach z pamięcią hierarchiczną koszt ładowania i przechowywania elementów macierzy wejściowej zwykle dominuje nad kosztem arytmetyki. Na pojedynczej maszynie jest to ilość danych przesyłanych między pamięcią RAM a pamięcią podręczną, podczas gdy na maszynie wielowęzłowej z pamięcią rozproszoną jest to ilość przesyłana między węzłami; w obu przypadkach nazywa się to przepustowością komunikacji . Naiwny algorytm wykorzystujący trzy zagnieżdżone pętle wykorzystuje szerokość pasma komunikacyjnego Ω( n ³ ) .

Algorytm Cannona , znany również jako algorytm 2D , jest algorytm komunikacji unikania że ścianki każdej matrycy wprowadza się do matrycy bloku, którego elementy są podmatryce o wielkości √ M / 3 o √ M / 3 , w którym M ma rozmiar szybkiej pamięci. Algorytm naiwny jest następnie używany na macierzach blokowych, obliczając iloczyny podmacierzy w całości w pamięci szybkiej. Zmniejsza to pasmo komunikacyjne do O ( n ³ / √ M ) , która jest optymalna asymptotycznie (algorytmów wykonywania omów ( N ³ ) obliczeń).

W układzie rozproszonym z p procesorów ułożonych w siatkę √ p przez √ p 2D, do każdego procesora można przypisać jedną podmacierz wyniku, a iloczyn może być obliczony z każdym procesorem transmitującym O ( n ² / √ p ) słów, co jest asymptotycznie optymalne przy założeniu, że każdy węzeł przechowuje minimum O ( n ² / p ) elementów. Można to poprawić dzięki algorytmowi 3D, który układa procesory w siatkę kostki 3D, przypisując każdy iloczyn dwóch podmacierzy wejściowych do jednego procesora. Podmacierze wyników są następnie generowane przez wykonanie redukcji w każdym wierszu. Algorytm ten przesyła O ( n ² / p ^2/3 ) słów na procesor, co jest asymptotycznie optymalne. Wymaga to jednak powtórzenia każdego elementu macierzy wejściowej p ^1/3 razy, a zatem wymaga współczynnika p ^1/3 więcej pamięci niż jest to potrzebne do przechowywania danych wejściowych. Ten algorytm można połączyć z Strassen, aby jeszcze bardziej skrócić czas pracy. Algorytmy „2,5D” zapewniają ciągły kompromis między wykorzystaniem pamięci a przepustowością komunikacji. Na nowoczesnych rozproszonych środowiskach obliczeniowych, takich jak MapReduce , opracowano wyspecjalizowane algorytmy mnożenia.

Algorytmy dla siatek

Mnożenie macierzy zakończone w krokach 2n-1 dla dwóch macierzy n×n na siatce okablowanej krzyżowo.

Istnieje wiele algorytmów mnożenia na siatkach . W przypadku mnożenia dwóch n × n na standardowej dwuwymiarowej siatce przy użyciu algorytmu 2D Cannona , można wykonać mnożenie w 3 n -2 krokach, chociaż liczba ta jest redukowana do połowy tej liczby dla powtarzanych obliczeń. Standardowa tablica jest nieefektywna, ponieważ dane z dwóch macierzy nie docierają jednocześnie i muszą być uzupełnione zerami.

Rezultat jest jeszcze szybszy na dwuwarstwowej siatce z krzyżowym drutem, gdzie potrzebne są tylko 2 n -1 kroków. Wydajność ulega dalszej poprawie w przypadku powtarzanych obliczeń, co prowadzi do 100% wydajności. Macierz siatek połączonych krzyżowo może być postrzegana jako szczególny przypadek nieplanarnej (tj. wielowarstwowej) struktury przetwarzania.

Zobacz też

• Złożoność obliczeniowa operacji matematycznych
• Złożoność obliczeniowa mnożenia macierzy
• Algorytm CYK, §Algorytm Valianta
• Mnożenie łańcucha macierzy
• Rzadkie mnożenie macierzy-wektora

Bibliografia

Dalsza lektura

• Buttari, Alfredo; Langou, Julien; Kurzak, Jakub; Dongarra, Jack (2009). „Klasa równoległych kafelkowych algorytmów algebry liniowej dla architektur wielordzeniowych”. Obliczenia równoległe . 35 : 38–53. arXiv : 0709.1272 . doi : 10.1016/j.parco.2008.10.002 . S2CID  955 .
• Goto, Kazushige; van de Geijn, Robert A. (2008). „Anatomia wysokowydajnego mnożenia macierzy”. Transakcje ACM na oprogramowaniu matematycznym . 34 (3): 1–25. CiteSeerX  10.1.1.140.3583 . doi : 10.1145/1356052.1356053 . S2CID  9359223 .
• Van Zee, Field G.; van de Geijn, Robert A. (2015). „BLIS: Ramy do szybkiego tworzenia instancji BLAS Funkcjonalność”. Transakcje ACM na oprogramowaniu matematycznym . 41 (3): 1-33. doi : 10.1145/2764454 . S2CID  1242360 .
• Jak zoptymalizować GEMM