Benutzer:Gardini/Arbeitsseite/Algorithmen zur Matrizenmultiplikation

Algorithmen zur Matrizenmultiplikation beschreiben, wie Matrizenmultiplikationen rechnerisch durchgeführt werden. Da es sich hierbei um eine Operation handelt, die beim wissenschaftlichen Rechnen ständig anfällt, zugleich aber mit großem Ressourcenbedarf verbunden ist, hat ihre Optimierung große praktische Bedeutung.

Aufgrund dessen wird ihr im High Performance Computing viel Aufmerksamkeit gewidmet, die Fachgemeinde gilt als „besessen von effizienter Matrix-Matrix-Multiplikation“.

Definition und Notation

Beispiel einer Matrixmultiplikation: Zur Berechnung des Elements $c_{12}$ (rot) wird die erste Zeile von $\mathbf {A}$ mit der zweiten Spalte von $\mathbf {B}$ (jeweils gelb) multipliziert, zur Berechnung des Elements $c_{33}$ (blau) die dritte Zeile von $\mathbf {A}$ mit der dritten Spalte von $\mathbf {B}$ (jeweils grün).

Im Folgenden wird die Multiplikation zweidimensionaler Matrizen über einem Körper $\mathbb {K}$ – typischerweise ist dies der Körper der reellen Zahlen $\mathbb {R}$ – betrachtet. Die erste, linke Faktormatrix $\mathbf {A} \in \mathbb {K} ^{m\times p}$ hat $m$ Zeilen und $p$ Spalten, die zweite, rechte Faktormatrix $\mathbf {B} \in \mathbb {K} ^{p\times n}$ hat $p$ Zeilen und $n$ Spalten. Ihr Matrixprodukt $\mathbf {C} =\mathbf {A} \cdot \mathbf {B}$ hat dann $m$ Zeilen und $n$ Spalten. Ziel jedes der hier behandelten Algorithmen ist es, alle $m\cdot n$ Elemente $c_{ij}=\sum _{k=1}^{p}a_{ik}b_{kj}$ des Matrixprodukts $\mathbf {C}$ zu berechnen. Der Einfachheit halber wird diese Rechenvorschrift im Weiteren als „Multiplikation der $i$ten Zeile von $\mathbf {A}$ mit der $j$ten Spalte von $\mathbf {B}$“ bezeichnet. Im Falle reeller (nicht komplexer) Matrizen entspricht dies der Berechnung des kanonischen Skalarprodukts der entsprechenden Zeilen- und Spaltenvektoren. Anhand der Zeilen- und Spaltenzahlen von $\mathbf {A}$ und $\mathbf {B}$ lassen sich $2^{3}=8$ verschiedene Fälle der Matrizenmultiplikation unterscheiden, je nachdem, ob $m$, $p$ und $n$ klein oder groß sind:

Klassifikation der Matrizenmultiplikationen nach Formparametern
(nach Goto und van de Geijn)

$m$	$p$	$n$	Veranschaulichung	Klassifikation
groß	groß	groß		GEMM
groß	klein	groß		GEPP
groß	groß	klein		GEMP
klein	groß	groß		GEPM
klein	klein	groß	BILD FEHLT	GEBP
groß	klein	klein	BILD FEHLT	GEPB
klein	klein	groß	BILD FEHLT	GEPDOT
klein	klein	klein	BILD FEHLT	GEBB

Die Zeitkomplexität von Algorithmen wird mithilfe der O-Notation angegeben. Diese gibt wieder, in welchem Maße die asymptotische Laufzeit mit wachsender Problemgröße zunimmt. Konstante Faktoren tauchen dabei definitionsgemäß unabhängig von ihrer Größe nicht in der O-Notation auf.

Sequentielle Algorithmen

Standardalgorithmus

Um die Produktmatrix $\mathbf {C}$ aufzubauen, kann man $\mathbf {A}$ zeilenweise von oben nach unten und $\mathbf {B}$ spaltenweise von links nach rechts auslesen und dabei in jedem Schritt eine Zeile von $\mathbf {A}$ mit einer Spalte von $\mathbf {B}$ multiplizieren. Hierdurch erhält man ein iteratives Verfahren, das sich in Python wie folgt darstellen lässt:

C = zeros((m,n))                            # C als Nullmatrix initialisieren
for i in range(m)                           # Schleife über die Zeilen von C
  for j in range(n)                         # Schleife über die Spalten von C
    for k in range(p)                       # Schleife über die Spalten von A bzw. Zeilen von B
      C[i][j] = C[i][j] + A[i][k] * B[k][j] # Produkte aus Zeilen- und Spaltenelementen aufsummieren

Sowohl im Hinblick auf das Ergebnis als auch auf die Zeitkomplexität des Algorithmus lassen sich die drei ineinander geschachtelten For-Schleifen beliebig vertauschen. Dabei fallen stets $m\cdot n\cdot p$ Multiplikationen und $m\cdot n\cdot (p-1)$ Additionen an, daher ist die Anzahl der benötigten Operationen von der Ordnung ${\mathcal {O}}(m\cdot n\cdot p)$, im Falle quadratischer Matrizen also ${\mathcal {O}}(n^{3})$. Dieser Algorithmus wird in der Literatur bisweilen Brute-Force-Cauchy-Algorithmus genannt.

Berücksichtigung der Speicherhierarchie

In Datenstrukturen mit row-major order sind die Elemente einer Matrix zeilenweise abgelegt, in solchen mit column-major order spaltenweise. Auf das Element $a_{ik}$ einer als C-type array float A[m][p] angelegten Matrix $A\in \mathbb {R} ^{m\times p}$ kann etwa per A[i][k] zugegriffen werden. Der Speicherzugriff im Array A erfolgt zunächst auf dessen Element i, welches selbst ein Array ist (Zeile $i$ von $A$), und innerhalb dieses Arrays wird erst auf das Element k zugegriffen, das $a_{ik}$ enthält.

Mathematisch und komplexitätstheoretisch ist es für die Berechnung des Matrixprodukts nach dem Standardalgorithmus zwar unerheblich, in welcher Reihenfolge über die Zeilen und Spalten der Faktormatrizen iteriert wird, doch hängt die praktisch erreichbare Rechenleistung stark davon ab, ob und wie die Reihenfolge die Lokalitätseigenschaft für Speicherzugriffe ausnutzt. Sind die Einträge der Matrizen etwa zeilenweise abgespeichert (row-major order, Standard zum Beispiel in den Programmiersprachen C/C++ oder der Programmbibliothek NumPy), so kann zur Berechnung eines Elements $c_{ij}$ der Produktmatrix $\mathbf {C}$ die $i$te Zeile der linken Faktormatrix $\mathbf {A}$ je nach Größe entweder ganz oder abschnittsweise als zusammenhängender Datenblock in den Cache geladen werden, die entsprechende $j$te Spalte der rechten Faktormatrix $\mathbf {B}$ dagegen nicht, da im Speicher zwischen jedem Element einer Spalte noch $p-1$ weitere Elemente aus den anderen Spalten liegen. Im schlechtesten Fall muss ein in den Cache geladener Ausschnitt aus $\mathbf {B}$ nach nur einer Rechenoperation wieder gelöscht und ein neuer Ausschnitt aus $\mathbf {B}$ geladen werden, um die nächste Rechenoperation durchführen zu können, obwohl der vorherige Ausschnitt in einem späteren Schritt noch hätte gebraucht werden können, wofür er dann ein weiteres Mal geladen werden muss. Die geringe Effizienz der Cachezugriffe, die damit verbunden ist, verlängert die tatsächliche Laufzeit auf dem Computer beträchtlich.

Ein erster Schritt zu besserer Cacheeffizienz besteht darin, die rechte Matrix $\mathbf {B}$ zu transponieren (transpose-and-multiply), sodass sie gemeinsam mit dem entsprechenden Ausschnitt aus der linken Matrix $\mathbf {A}$ zeilenweise ausgelesen werden kann. Ist die Datenstruktur dagegen spaltenweise angelegt (column-major order, Standard zum Beispiel in den Programmiersprachen Fortran und Matlab), muss entsprechend die linke Matrix $\mathbf {A}$ transponiert werden, um spaltenweise ausgelesen werden zu können. Der zusätzliche Lese- und Schreibaufwand für die Transposition ist – insbesondere bei größeren Matrizen – laufzeitmäßig vernachlässigbar klein gegenüber dem Aufwand für die eigentliche Berechnung; allerdings wird zusätzlicher Speicherplatz für die transponierte Matrix benötigt.

Beispielcodes zu transpose-and-multiply

Matrizen in row-major order, Datenstruktur als „geplättete Matrix“ offengelegt (C++)

Matrizen in column-major order, Datenstruktur als Array von Arrays verborgen (Fortran)

/* Generische Funktion zu Cache-sensitiver Matrixmultiplikation in C++ */ template <typename value_t, typename index_t> void matmul_omp( value_t* A, // Zeiger auf erstes Element der Matrix A value_t* B, // Zeiger auf erstes Element der Matrix B value_t* C, // Zeiger auf erstes Element der Matrix C index_t rows_A, // Anzahl der Zeilen von A index_t cols_A, // Anzahl der Spalten von A index_t rows_B, // Anzahl der Zeilen von B index_t cols_B) // Anzahl der Spalten von B { /* Wohldefiniertheit des Matrixprodukts prüfen */ assert(cols_A == rows_B); /* Speicherplatz für transponierte Matrix auf dem Heap allokieren */ value_t* Bt = new value_t[cols_B*rows_B]; /* Transposition der rechten Matrix */ for (index_t i = 0; i < rows_B; i++) for (index_t j = 0; j < cols_B; j++) Bt[j*rows_B+i] = B[i*cols_B+j]; /* Matrixmultiplikation */ for (index_t i = 0; i < rows_A; i++) { for (index_t j = 0; j < cols_B; j++) { C[i*cols_B+j] = (value_t)(0); for (index_t k = 0; k < cols_A; k++) C[i*cols_B+j] += A[i*cols_A+k] * Bt[j*rows_B+k]; } } }

! Funktion zur Matrixmultiplikation nach Transposition in Fortran subroutine matmul_tam(& A, B, C, & ! Zeiger auf Matrizen A, B, C rows_A, cols_A, & ! Anzahl der Zeilen und Spalten von A rows_B, cols_B) ! Anzahl der Zeilen und Spalten von B ! Deklaration übergebener Variablen und nötiger Laufindizes integer :: rows_A, cols_A, rows_B, cols_B, i, j, k real(kind=8), dimension(rows_A, cols_A) :: A real(kind=8), dimension(rows_B, cols_B) :: B real(kind=8), dimension(rows_A, cols_B) :: C ! Deklaration der transponierten Matrix real(kind=8), dimension(cols_A, rows_A) :: At ! Transposition der linken Matrix do j = 1, cols_A do i = 1, rows_A At(j,i) = A(i,j) end do end do ! Matrixmultiplikation do j = 1, cols_B do i = 1, rows_A C(i,j) = 0.0 do k = 1, cols_A C(i,j) = C(i,j) + At(k,i) * B(k,j) end do end do end do end subroutine

Effizienzsteigerung lässt sich auch durch Unterteilung der Matrizen in Submatrizen – also zweidimensionale Datenblöcke (tiles) – und die Verwendung passender blocked algorithms erreichen. Dabei müssen sowohl die verschiedenen Level der Cache-Hierarchie berücksichtigt werden, als auch der Umstand, dass nicht nur bezüglich Daten und Befehlen cache misses auftreten können, sondern auch bei der Adressverwaltung im Translation Lookaside Buffer. Um solches cache blocking möglichst auf jeder Ebene der Speicherhierarchie zu verwirklichen, können die Algorithmen in einem geschichteten Ansatz so gestaltet werden, dass auf jedem Cache-Level ein möglichst großer Block von entweder $\mathbf {A}$, $\mathbf {B}$ oder $\mathbf {C}$ vorrätig gehalten wird, welcher das jeweilige Cache-Level größtenteils ausfüllt, während kleinere Blöcke der beiden anderen Matrizen nach und nach aus dem nächstlangsameren Cache-Level geladen und schließlich – auf L1-Ebene angelangt – verrechnet werden.

Parallele Algorithmen

Für sehr große Matrizen – wie sie beispielsweise im Rahmen numerischer Verfahren zur Behandlung von Vielteilchensystemen regelmäßig anfallen – dauern sequentielle Implementierungen wie die bisher dargestellten viel zu lange. Um solche Matrixmultiplikationen in vertretbarer Zeit ausführen zu können, bedarf es einer effizienten Parallelisierung.

Shared-Memory-Architekturen

Auf Mehrkernprozessoren lassen sich die Iterationen sowohl des Transpositions- als auch des Multiplikationsschritts auf mehrere Threads aufteilen. Da jedes Element des Matrixprodukts unabhängig von den anderen berechnet werden kann, lässt sich bei gemeinsamem Arbeitsspeicher Multithreading etwa mit OpenMP bereits mithilfe weniger Compileranweisungen umsetzen. Effiziente Algorithmen müssen jedoch auch hier das Fassungsvermögen des Cache berücksichtigen. Die Daten können nach einem ein- oder zweidimensionalen Schema auf die Threads verteilt werden, wobei die blockweise Verteilung der Daten als zweidimensionale Submatrizen oft bessere Performance bringt. Eine Alternative zu blocked algorithms, die aufgabenweise Daten auf Threads verteilen und wieder einsammeln, sind hier algorithms-by-blocks, die sequentielle Teilaufgaben dynamisch zur Laufzeit an die Threads verteilen, um eine bessere Lastverteilung zu erreichen.

Distributed-Memory-Architekturen

Im Falle von Hochleistungsrechnern, die aus einem großen Prozessornetzwerk ohne gemeinsamen Arbeitsspeicher bestehen, bedarf es eigenständiger Algorithmen, die den Datenaustausch zwischen den verschiedenen Knoten des Netzwerks organisieren. Wichtige Faktoren dabei sind Speicherbedarf und Kommunikationsoverhead. Um letzteren möglichst klein zu halten, sollten Rechen- und Kommunikationsphasen überlappen und die bei einem Aufruf übermittelten Datenmengen möglichst groß sein, um insgesamt weniger Kommunikationsaufrufe zu haben. Zahlreiche parallele Algorithmen zur Matrixmultiplikation mit Multiprocessing wurden entwickelt, wobei sich grob zwei Kategorien unterscheiden lassen: Algorithmen wie der Cannon-Algorithmus, die in jedem Schritt Datenblöcke über Punkt-zu-Punkt-Kommunikation weiterschieben, und solche wie SUMMA, die in jedem Schritt Kollektivkommunikationen durchführen.

Literatur

• C. Addison, Y. Ren, M. van Waveren: OpenMP Issues Arising in the Development of Parallel BLAS and LAPACK libraries. Third European Workshop on OpenMP, Barcelona, 2001. doi:10.1155/2003/278167
• Paolo D’Alberto, Alexandru Nicolau: Adaptive Strassen and ATLAS's DGEMM: a fast square-matrix multiply for modern high-performance systems. In: Proceedings of the Eighth International Conference on High-Performance Computing in Asia-Pacific Region (HPCASIA’05), 2005. doi:10.1109/HPCASIA.2005.18
• Lynn Elliot Cannon: A cellular computer to implement the Kalman Filter Algorithm. Dissertation, Montana State University, 1969.
• Anthony Danalis, Ki-Yong Kim, Lori Pollock, Martin Swany: Transformations to Parallel Codes for Communication-Computation Overlap. In: Proceedings of the ACM/IEEE Conference on Supercomputing, Seattle, 2005. doi:10.1109/SC.2005.75
• Ernie Chan, Enrique S. Quintana-Ortí, Gregorio Quintana-Ortí, Robert A. van de Geijn: SuperMatrix Out-of-Order Scheduling of Matrix Operations for SMP and Multi-Core Architectures. In: Proceedings of the nineteenth annual ACM symposium on Parallel algorithms and architectures (SPAA ’07), 2007, S. 116–125. doi:10.1145/1248377.1248397
• Jack J. Dongarra, Jeremy du Croz, Sven Hammarling, Iain Duff: A Set of Level 3 Basic Linear Algebra Subprograms. In: ACM Transactions on Mathematical Software, Band 16, Nr. 1, 1990, S. 1–17. doi:10.1145/77626.79170
• Jack J. Dongarra, Jean-François Pineau, Yves Robert, Zhiao Shi, Frédéric Vivien: Revisiting Matrix Product on Master-Worker Platforms. In: International Journal of Foundations of Computer Science, Band 19, Nr. 6, 2008, S. 1317–1336. doi:10.1142/S0129054108006303
• Erik Elmroth, Fred Gustavson, Isak Jonsson, Bo Kågström: Recursive Blocked Algorithms and Hybrid Data Structures for Dense Matrix Library Software. In: SIAM Review, Band 46, Nr. 1, 2004, S. 3–45. doi:10.1137/S0036144503428693
• Robert A. van de Geijn, Jerrell Watts: SUMMA: scalable universal matrix multiplication algorithm. In: Concurrency and Computation, Band 9, Nr. 4, 1997, S. 255–274. Online auf Netlib.org (Repository der University of Tennessee Knoxville und des Oak Ridge National Laboratory).
• Pawel Gepner, Victor Gamayunov, David L. Fraser: Effective Implementation of DGEMM on modern multicore CPU. In: Procedia Computer Science, Band 9, 2012, S. 126–135. doi:10.1016/j.procs.2012.04.014
• Pawel Gepner, Victor Gamayunov, David L. Fraser, Eric Houdard, Ludovic Sauge, Damien Declat, Mathieu Dubois: Evaluation of DGEMM Implementation on Intel Xeon Phi Coprocessor. In: Journal of Computers, Band 9, Nr. 7, 2014, S. 1566–1571. doi:10.4304/jcp.9.7.1566-1571
• Kazushige Goto, Robert A. van de Geijn: On Reducing TLB Misses in Matrix Multiplication. Tech. rep. CS-TR-02-55, University of Texas at Austin, 2002. Online auf der Website der Pennsylvania State University.
• Kazushige Goto, Robert A. van de Geijn: Anatomy of High-Performance Matrix Multiplication. In: ACM Transactions on Mathematical Software, Band 34, Nr. 3, 2008, S. 12:1–25. doi:10.1145/1356052.1356053
• John A. Gunnels, Greg M. Henry, Robert A. van de Geijn: A Family of High-Performance Matrix Multiplication Algorithms. In: Vassil N. Alexandrov, Jack J. Dongarra, Benjoe A. Juliano, Rene S. Renner, C. J. Kenneth Tan (Hrsg.): Proceedings of the International Conference on Computational Sciences Part I (LNCS 2073). Springer, Berlin/Heidelberg/New York 2001, ISBN 3-540-42232-3, S. 51–60.
• John A. Gunnels, Fred G. Gustavson, Greg M. Henry, Robert A. van de Geijn: A Family of High-Performance Matrix Multiplication Algorithms. In: Jack Dongarra, Kaj Madsen, Jerzy Wasniewski (Hrsg.): Applied Parallel Computing: State of the Art in Scientific Computing (LNCS 3732). Springer, Berlin/Heidelberg/New York 2006, ISBN 3-540-29067-2, S. 256–265.
• Himanshu Gupta, P. Sadayappan: Communication efficient matrix multiplication on hypercubes. In: Proceedings of the sixth annual ACM symposium on Parallel algorithms and architectures (SPAA ’94), 1994, S. 320–329. doi:10.1145/181014.181434
• José R. Herrero, Juan J. Navarro: Improving Performance of Hypermatrix Cholesky Factorization. In: Harald Kosch, László Böszörményi, Hermann Hellwagner (Hrsg.): Euro-Par 2003 Parallel Processing, LNCS 2790. Springer, Berlin/Heidelberg/New York 2003, ISBN 3-540-40788-X, S. 461–469.
• José R. Herrero: A Framework for Efficient Execution of Matrix Computations. Dissertation, Universitat Politècnica de Catalunya, 2006. Online auf der Website der Universitat Politècnica de Catalunya.
• Jianyu Huang, Tyler M. Smith, Greg M. Henry, Robert A. van de Geijn: Strassen’s Algorithm Reloaded. In: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC ’16), 2016. doi:10.1109/SC.2016.58
• Wladimir Wladimirowitsch Kljujew, Nikolai Iwanowitsch Kokowkin-Schtscherbak: О минимизации числа арифметических операций при решении линейных алгебраических систем уравнений. In: Журнал вычислительной математики и математической физики, Band 5, Nr. 1, 1965, S. 21–33. Original auf Russisch, englische Übersetzung erschienen unter V. V. Klyuyev, N. I. Kokovkin-Shcherbak: On the minimization of the number of arithmetic operations for the solution of linear algebraic systems of equations. In: USSR Computational Mathematics and Mathematical Physics, Band 5, Nr. 1, 1965, S. 25–43. doi:10.1016/0041-5553(65)90065-0
• Jakub Kurzak, Jack Dongarra: Implementing Linear Algebra Routines on Multi-Core Processors with Pipelining and a Look Ahead / LAPACK Working Note 178. In: Proceedings of the 8th international conference on Applied parallel computing: state of the art in scientific computing, 2006. doi:10.1007/978-3-540-75755-9_18
• Hyuk-Jae Lee, James P. Robertson, José A. B. Fortes: Generalized Cannon’s Algorithm for Parallel Matrix Multiplication. In: Procs. 11th. Intl. Conf. on Supercomputing (ICS ’97), Wien 1997, S. 44ff. doi:10.1145/263580.263591
• Monica S. Lam, Edward E. Rothberg, Michael E. Wolf: The Cache Performance and Optimizations of Blocked Algorithms. In: ACM SIGPLAN Notices, Band 26, Nr. 4, 1991, S. 63–74. doi:10.1145/106973.106981
• Roktaek Lim, Yeongha Lee, Raehyun Kim, Jaeyoung Choi: An implementation of matrix–matrix multiplication on the Intel KNL processor with AVX-512. In: Cluster Computing, Band 21, 2018, S. 1785–1795. doi:10.1007/s10586-018-2810-y
• Bryan Marker, Field G. Van Zee, Kazushige Goto, Gregorio Quintana-Ortí, Robert A. van de Geijn: Toward Scalable Matrix Multiply on Multithreaded Architectures. In: Anne-Marie Kermarrec, Luc Bougé, Thierry Priol (Hrsg.): Euro-Par 2007 Parallel Processing, LNCS 4641. Springer, Berlin/Heidelerg/New York 2007, ISBN 3-540-74465-7, S. 748–757.
• Gregorio Quintana-Ortí, Enrique S. Quintana-Ortí, Robert A. van de Geijn, Field G. van Zee, Ernie Chan: Programming Matrix Algorithms-by-Blocks for Thread-Level Parallelism. In: ACM Transactions on Mathematical Software, Band 36, Nr. 3, 2009, S. 14:1–26. doi:10.1145/1527286.1527288
• Bertil Schmidt, Jorge González-Domínguez, Christian Hundt, Moritz Schlarb: Parallel Programming: Concepts and Practice. Morgan Kaufmann (Elsevier), Cambridge 2018, ISBN 978-0-12-849890-3.
• Volker Strassen: Gaussian elimination is not optimal. In: Numerische Mathematik, Band 13, S. 354–356, 1969. doi:10.1007/BF02165411
• Field G. van Zee, Paolo Bientinesi, Tze Meng Low, Robert A. van de Geijn: Scalable Parallelization of FLAME Code via the Workqueuing Model. In: ACM Transactions on Mathematical Software, Band 34, Nr. 2, 2008, S. 10:1–29. doi:10.1145/1326548.1326552