Parallele Matrizenmultiplikation

Der Fox-Algorithmus, benannt nach Geoffrey C. Fox, ist ein Algorithmus zur parallel, auf p Prozessoren gleichzeitig ausgeführten Matrizenmultiplikation. In der Entwicklung des Algorithmus wurde eine Topologie genutzt, in der die Prozessoren in einem Hyperwürfel angeordnet sind. Als Speichermodell wird das Modell des verteilten Speicher genutzt. Hierbei hat jeder Prozessor seinen eigenen privaten Adressraum und die Kommunikation zwischen Prozessoren ist nachrichtenbasiert. Der Anspruch bei der Entwicklung war es einen effizienten und einfach einzusetzenden Algorithmus, zur Nutzung in wissenschaftlichen Berechnungen, zu erstellen. Des Weiteren war es ein Bestreben einen Algorithmus zu entwickeln, der auch für zukünftige Maschinen mit vielen Prozessoren skaliert. Essentiell für den Algorithmus sind Systolische Arrays. Dieses bezeichnet ein Pipe-Netzwerk, durch welches Datenströme hindurchgetaktet werden.[1]

Es werden zwei ${\displaystyle {n\times n}}$-Matrizen A und B multipliziert. Hierbei sind A und B vollbesetzte Matrizen. Die Elemente der Matrizen werden mit Hilfe der i- und j-Koordinaten in der Matrix auf die Prozessoren, die ebenfalls eine i- und j-Koordinate innerhalb des Hyperwürfel besitzen, verteilt. Somit hat zu Begin der Berechnung jeder Prozessor ein Element der Matrix A und ein Element der Matrix B.[2] Für Matrizen, bei denen gilt ${\displaystyle {n=p}}$ erhält jeder Prozessor einen Eintrag der Matrix A und einen Eintrag der Matrix B. Für Matrizen, bei denen gilt ${\displaystyle {n>p}}$ erhält jeder Prozessor eine Teilmatrix der Größe ${\displaystyle ({\frac {n}{\sqrt {p}}})\times ({\frac {n}{\sqrt {p}}})}$ von A und B.[3] Im Fall der Teilmatrizen führt jeder Prozessor eine eigene Matrizenmultiplikation der Teilmatrizen durch.

Während des Algorithmus werden die Elemente der Matrix A mittels Broadcast in den Prozessorreihen, also allen Prozessoren mit gleicher j-Koordinate, verteilt. Die Elemente der Matrix B werden an den darüber liegenden Prozessor weitergegeben. Bei Prozessoren am oberen Rand werden die Elemente an den Prozessor am unteren Ende weitergegeben.

Bei initialer Startkonfiguration hält jeder Prozessor das Element der Matrix A und B, das die gleichen Werte für die Koordinaten i und j besitzt, wie der Prozessor selbst im Hyperwürfel.

Der Algorithmus besteht aus folgenden Schritten:

1. Broadcast der Elemente auf der Diagonalen der Matrix A an die Prozessoren mit gleicher j-Koordinate. Anschließend halten alle Prozessoren in der ersten Reihe das Element ${\displaystyle A_{00}}$, alle Prozessoren in der zweiten Reihe das Element ${\displaystyle A_{11}}$. Dieses Muster setzt sich für alle Prozessoren fort.
2. Multiplikation des empfangenen Elements der Matrix A mit dem derzeit im Prozessor vorhandenen Element der Matrix B.
3. Alle Prozessoren reichen ihr derzeitiges Element der Matrix B an den Prozessor in der gleichen Spalte, eine Stelle über sich, weiter. Die Prozessoren am Rand übergeben ihr Element an die Prozessoren am unteren Rand.
4. Broadcast der Elemente auf der 'Diagonalen + 1' der Matrix A an die Prozessoren mit gleicher j-Koordinate. Anschließend halten alle Prozessoren in der ersten Reihe das Element ${\displaystyle A_{01}}$, alle Prozessoren in der zweiten Reihe das Element ${\displaystyle A_{12}}$. Dieses Muster setzt sich für alle Prozessoren fort.
5. Multiplikation des empfangenen Elements der Matrix A mit dem derzeit im Prozessor vorhandenen Element der Matrix B und anschließende Addition des Ergebnisses mit dem Ergebnis aus der vorherigen Multiplikation.

Dieses Muster setzt sich fort, bis die Elemente von B wieder in der Ausgangsstellung sind und alle Nebendiagonalen von A verwendet wurden.[4]

In Pseudocode kann der Algorithmus wie folgt implementiert werden:

FOX(A, B):
1 //Prozessoreinheit(i,j)
2 ${\displaystyle a:=A_{ij};}$
3 ${\displaystyle b:=B_{ij};}$
4 ${\displaystyle c_{ij}:=0;}$
5 for ${\displaystyle (l:=0;l<{\sqrt {p}};l++)}$ {
6    PE${\displaystyle (i,(i+l)mod{\sqrt {p}})}$ broadcasts a to PE(i,j) for ${\displaystyle j\in {0..{\sqrt {p}}}\land j\neq (i+l)mod{\sqrt {p}};}$
7    ${\displaystyle c_{ij}:=c_{ij}+a*b;}$
8    concurrently {
9        send b to PE${\displaystyle ((i+1)mod{\sqrt {p}},j);}$
10   } with {
11       receive b' from PE${\displaystyle ((i-1)mod{\sqrt {p}},j);}$
12   }
13   ${\displaystyle b:=b';}$
14 }

Sollen zwei Matrizen der Größe ${\displaystyle n\times n}$ mittels ${\displaystyle {\sqrt {q}}\times {\sqrt {q}}}$ Prozessoren multipliziert werden. So erhält jeder Prozessor zwei Teilmatrizen der Größe ${\displaystyle ({\frac {n}{\sqrt {p}}})\times ({\frac {n}{\sqrt {p}}})}$.

Das Muster Broadcast-Multiplikation-Rotation wiederholt sich ${\displaystyle {\sqrt {p}}}$-mal, bis das Ergebnis feststeht.

Insgesamt hat der Algorithmus eine Berechnungszeit von ${\displaystyle T={\sqrt {q}}*(2m^{3}t_{flop}+2m^{2}t_{comm}+({\sqrt {p}}1)t_{start})}$

Diese setzt sich zusammen aus den einzelnen Teilschritten:

1. Broadcast der Submatrizen der Matrix A: ${\displaystyle m^{2}t_{comm}+({\sqrt {p}}-1)t_{start}}$
2. Rotation der Submatrizen der Matrix B um die j-Achse: ${\displaystyle m^{2}t_{comm}}$
3. Multiplikation der Submatrizen und Addition zum vorherigen Teilergebnis: ${\displaystyle 2m^{3}t_{flop}}$

Hierbei ist ${\displaystyle t_{comm}}$ die Zeit, um eine Gleichkommazahl zwischen Prozessoren zu übertragen, ${\displaystyle t_{start}}$ die Startzeit um die Pipeline zu befüllen und ${\displaystyle t_{flop}}$ die Zeit für eine Multiplikation oder Addition von Gleichkommazahlen.

Die Analyse gilt für alle Maschinen, die der flynnsche Klassifikation MIMD zugeordnet werden können und das Speichermodell des verteilten Speichers haben.[5]

Ein Nachteil des Algorithmus ist die Skalierbarkeit für die Multiplikation von zwei ${\displaystyle n\times n}$-Matrizen. Es können höchstens ${\displaystyle n^{2}}$ Prozessoren benutzt werden, womit die Laufzeit in ${\displaystyle \Omega (n)}$ liegt, da es ${\displaystyle \Theta (n^{3})}$ Operationen im seriellen Algorithmus sind.[6]

Die Grafik zeigt die Matrix A, die rot eingefärbt wurde für eine bessere Wiedererkennung im Video zum Ablauf des DNS-Algorithmus.

Der DNS-Algorithmus, benannt nach Elizier Dekel, David Nassimi und Sartaj Sahni, wurde 1981 von diesen veröffentlicht. Der Algorithmus wurde entwickelt für die parallele Multiplikation von Matrizen auf general purpose Prozessoren. Der Algorithmus richtet sich an parallele Computer die der flynnsche Klassifikation SIMD zugeordnet werden können. Diese Computer haben gemeinsam, dass sie aus p verarbeitenden Elementen (processing elements) PE bestehen. Jedes dieser verarbeitenden Elemente kann die Standard arithmetischen und logischen Operationen durchführen.[7] Insgesamt können bis zu ${\displaystyle n^{3}}$ Prozessoren benutzt werden. Diese sind hierbei in einem dreidimensionalen Hyperwürfel angeordnet. Bei ${\displaystyle n^{3}}$ skalaren Multiplikationen folgt daraus, dass jeder Prozessor eine skalare Multiplikation durchführen muss.

Der Algorithmus wurde vor allem entwickelt, um Probleme der Graphentheorie zu lösen, beispielsweise das Finden des kürzesten Weges zwischen allen All-Pairs Shortest Path oder das bestimmen von Radius, Durchmesser und Mittelpunkt eines Graphen.[8]

Es werden zwei ${\displaystyle {n\times n}}$-Matrizen, A und B multipliziert. Diese sind vollbesetzte Matrizen. Konkret werden bei diesem Algorithmus ${\displaystyle m:={\sqrt[{3}]{p}}}$ und damit ${\displaystyle m^{3}}$ Prozessoren genutzt, die in einem Würfel angeordnet sind. Die Elemente der Matrizen werden mittels der i- und j-Koordinate den Prozessoren mit passenden Koordinaten innerhalb des dreidimensionalen Hyperwürfels zugeordnet. Die Multiplikation soll dann auf ${\displaystyle m\times m}$-Matrizen durchgeführt werden. Die Einträge der Matrizen sind wiederum Matrizen mit der Größe ${\displaystyle {\frac {m}{n}}\times {\frac {m}{n}}}$. Hierbei soll ${\displaystyle m}$ ein Teiler von ${\displaystyle n}$ sein.

Die Grafik zeigt die Matrix B, die zur besseren Wiedererkennung im Video zum Ablauf des DNS-Algorithmus gelb eingefärbt wurde.

Um Prozessoren zu adressieren, nutzen wir die Schreibweise ${\displaystyle PE(i,j,k)}$. Der initiale Zustand kann daher ausgedrückt werden durch ${\displaystyle A(i,k,0),B(0,k,j)}$. Es befinden sich die Elemente der Matrix A auf den Prozessoren auf der Vorderseite des Hyperwürfels und die Elemente der Matrix B auf der linken Seitenfläche des Hyperwürfels. Die Prozessoren der linken vorderen Kante haben somit bereits ein Element der Matrix A und ein Element der Matrix B. Dies ist die initiale Startkonfiguration.

Der Algorithmus besteht aus folgenden Schritten:

1. Broadcast der Matrizen A und B über die ${\displaystyle n^{3}}$ Prozessoren. Hierbei werden die Elemente der Matrix A (befindet sich auf der Vorderseite des Würfels) entlang der j-Dimension nach hinten gesendet. Die Elemente der Matrix B (befindet sich auf der linken Außenfläche des Würfels) werden entlang der k-Dimension nach rechts gesendet. Nach Abschluss des Broadcast hält jeder Prozessor 1 Element der Matrix A und ein Element der Matrix B.
2. Nun führt jeder Prozessor die Multiplikation der vorhandenen Elemente durch.
3. In der letzten Phase wird die Summe der Einzelergebnisse gebildet. Dies geschieht, indem von oben nach unten die Teilergebnisse gesendet und auf dem Weg addiert werden. Nach Abschluss dieser Operation befindet sich die Ergebnismatrix auf der Bodenfläche des Hyperwürfels.

Eine Implementierung in Pseudocode könnte wie folgt aussehen.

DNS(A, B):
1 store ${\displaystyle a_{ik}}$ in PE${\displaystyle (i,k,1);}$
2 store ${\displaystyle b_{kj}}$ in PE${\displaystyle (1,k,j);}$
3 PE${\displaystyle (i,k,1)}$ broadcasts ${\displaystyle a_{ik}}$ to PEs${\displaystyle (i,k,j)}$ for ${\displaystyle j\in {1..n};}$
4 PE${\displaystyle (1,k,j)}$ broadcasts ${\displaystyle b_{kj}}$ to PEs${\displaystyle (i,k,j)}$ for ${\displaystyle i\in {1..n};}$
5 compute ${\displaystyle c_{ikj}:=a_{ik}*b_{kj}}$ on PE${\displaystyle (i,k,j);}$
6 PEs${\displaystyle (i,k,j)}$ foreach ${\displaystyle (k\in {1..n})}$ compute ${\displaystyle c_{ij}:=\sum _{k=1}^{n}c_{ikj}}$ to PE${\displaystyle (i,1,j);}$

Angenommen es sollen zwei Matrizen der Größe ${\displaystyle {n\times n}}$ mittels ${\displaystyle m:={\sqrt[{3}]{p}}}$ und damit ${\displaystyle m^{3}}$ Prozessoren multipliziert werden. Es gilt die Anzahl Prozessoren ${\displaystyle p=q^{3}}$.

Insgesamt hat der Algorithmus hat eine Berechnungszeit von[9]:

  ${\displaystyle T_{p}={\frac {n^{3}}{p}}+t_{s}\log p+t_{w}{\frac {n^{2}}{p^{(}{\frac {2}{3}})}}\log p}$
  ${\displaystyle t_{s}}$ steht für die Startup Time. Diese wird benötigt um eine Nachricht im Sendeprozessor zu handhaben. Hierbei ist die Zeit mitinbegriffen, die benötigt wird um die Nachricht zu erstellen, den Routing-Algorithmus durchzuführen und eine Verbindung zwischen beiden Prozessoren zu etablieren.
  ${\displaystyle t_{w}}$ steht für die Transferzeit pro Wort. In diesem Fall beispielsweise die Dauer für die Übertragung eines Fließkommazahl.

Die Berechnungszeit setzt sich zusammen aus[9]:

  1 Die Broadcastoperation ausgeführt für beide Matrizen: ${\displaystyle t_{s}\log q+t_{w}({\frac {n}{q}})^{2}\log q}$
  2 Die Reduktion für die Ergebnismatrix C: ${\displaystyle t_{s}\log q+t_{w}({\frac {n}{q}})^{2}\log q}$
  3 Multiplikation der Teilmatrizen: ${\displaystyle ({\frac {n}{q}})^{3}}$

Hieraus ergibt sich[9]:

   ${\displaystyle T_{p}=({\frac {n}{q}})^{3}+3t_{s}\log q+3t_{w}({\frac {n}{q}})^{2}\log q}$

Durch einsetzen der obigen Bedingung ${\displaystyle p=q^{3}}$ erhält man die obige Berechnungszeit[9]:

   ${\displaystyle T_{p}={\frac {n^{3}}{p}}+t_{s}\log p+t_{w}{\frac {n^{2}}{p^{(}{\frac {2}{3}})}}\log p}$

Der Algorithmus ist für ${\displaystyle n^{3}=\Theta (p\log p)}$ oder ${\displaystyle p=\Omega ({\frac {n^{3}}{\log n}})}$ kostenoptimal.[10]