Ich habe einen Schlüsselalgorithmus, in dem die meiste Laufzeit für die Berechnung eines dichten Matrixprodukts aufgewendet wird:
%Vor% Basierend auf diesen Dimensionen, nach den Lehren aus dem Problem der Matrixkettenmultiplikation , ist es klar, dass es optimal in einem ist Anzahl der Operationen Sinn, um die Berechnung als A*(A'*Y) zu strukturieren. Meine aktuelle Implementierung tut das, und die Leistungssteigerung durch das Erzwingen dieser Assoziativität für den Ausdruck ist spürbar.
Meine Anwendung ist in C ++ für die x86_64-Plattform geschrieben. Ich benutze die Eigen Bibliothek für lineare Algebra, mit Intels Math Kernel Library als Backend. Eigen ist in der Lage, die BLAS-Schnittstelle von IMKL zu verwenden, um die Multiplikation durchzuführen, und der Boost vom Wechsel zu Eigens nativer SSE2-Implementierung zu Intels optimierter AVX-basierter Implementierung auf meinem Sandy-Bridge-Rechner ist ebenfalls von Bedeutung.
Der Ausdruck A * (A.adjoint() * Y) (in Eigensprache geschrieben) wird jedoch in zwei allgemeine Matrix-Matrix-Produkte zerlegt (Aufrufe in die Routine xGEMM BLAS), wobei dazwischen eine temporäre Matrix erstellt wird. Ich frage mich, ob ich zu einer Implementierung gelangen kann, die schneller ist als die generische, die ich jetzt habe, wenn ich zu einer spezialisierten Implementierung gehe, um den gesamten Ausdruck auf einmal zu bewerten. Ein paar Beobachtungen, die mich dazu bringen, dies zu glauben sind:
Unter Verwendung der oben beschriebenen typischen Dimensionen passt die Eingabematrix A normalerweise nicht in den Cache. Daher wäre das spezifische Speicherzugriffsmuster, das zum Berechnen des Drei-Matrix-Produkts verwendet wird, der Schlüssel. Natürlich wäre es auch vorteilhaft, die Schaffung einer temporären Matrix für das Teilprodukt zu vermeiden.
A und seine konjugierte Transponierte haben offensichtlich eine sehr verwandte Struktur, die möglicherweise ausgenutzt werden könnte, um das Speicherzugriffsmuster für den Gesamtausdruck zu verbessern.
Gibt es Standardtechniken, um diese Art von Ausdruck cachefreundlich zu implementieren? Die meisten Optimierungsmethoden, die ich für die Matrixmultiplikation gefunden habe, sind für den Standardfall A*B , nicht für größere Ausdrücke. Ich bin mit den Mikrooptimierungsaspekten des Problems vertraut, z. B. der Übersetzung in die entsprechenden SIMD-Befehlssätze, aber ich suche nach Referenzen, um diese Struktur möglichst speicherfreundlich zu durchbrechen. p>
Bearbeiten: Aufgrund der Antworten, die bisher eingegangen sind, denke ich, dass ich oben ein wenig unklar war. Die Tatsache, dass ich C ++ / Eigen benutze, ist wirklich nur ein Implementierungsdetail aus meiner Sicht auf dieses Problem. Eigen leistet hervorragende Arbeit bei der Implementierung von Expression-Templates, aber die Auswertung dieses Problemtyps als einfacher Ausdruck wird nicht unterstützt (nur Produkte mit zwei allgemeinen dichten Matrizen).
Auf einer höheren Ebene, als die Ausdrücke von einem Compiler ausgewertet werden würden, suche ich nach einer effizienteren mathematischen Aufschlüsselung der zusammengesetzten Multiplikationsoperation, mit einer Neigung unnötige redundante Speicherzugriffe zu vermeiden aufgrund der gemeinsamen Struktur von A und seiner konjugierten Transponierten. Das Ergebnis wäre wahrscheinlich schwierig in reinem Eigen zu implementieren, so würde ich wahrscheinlich nur in einer spezialisierten Routine mit SIMD-Intrinsics implementieren.
Dies ist keine vollständige Antwort (noch - und ich bin mir nicht sicher, dass es eins wird).
Denken wir zuerst ein wenig an die Mathematik. Da die Matrixmultiplikation assoziativ ist, können wir beides tun (A * A ') Y oder A (A' * Y).
Fließkommaoperationen für (A * A ') * Y
%Vor%Fließkommaoperationen für A * (A '* Y)
%Vor%Da k viel kleiner ist als m und n, ist klar, warum der zweite Fall viel schneller ist.
Aber durch die Symmetrie könnten wir die Anzahl der Berechnungen für A * A 'um zwei reduzieren (obwohl dies mit SIMD vielleicht nicht einfach ist), um die Anzahl der Gleitkommaoperationen von (A * A') zu reduzieren. ) * Y bis
%Vor% Wir wissen, dass sowohl m als auch n größer als k sind. Wählen wir eine neue Variable für m und n mit dem Namen z und finden Sie heraus, wo Fall eins und zwei gleich sind:
Solange m und n beide mehr als zweimal k sind, hat der zweite Fall weniger Fließkommaoperationen. In Ihrem Fall sind m und n beide mehr als 100 und k weniger als 10, weshalb im zweiten Fall viel weniger Fließkommaoperationen verwendet werden.
In Bezug auf effizienten Code. Wenn der Code für eine effiziente Verwendung des Cachespeichers optimiert ist (wie MKL und Eigen sind), dann ist eine große dichte Matrixmultiplikation rechnergebunden und nicht speichergebunden, so dass Sie sich nicht um den Cache kümmern müssen. MKL ist schneller als Eigen, da MKL AVX benutzt (und vielleicht fma3 jetzt?).
Ich denke nicht, dass Sie das effizienter machen können, als Sie es bereits mit dem zweiten Fall und MKL (über Eigen) gemacht haben. Aktivieren Sie OpenMP, um maximale FLOPS zu erhalten.
Sie sollten die Effizienz berechnen, indem Sie FLOPS mit den Peak-FLOPS Ihres Prozessors vergleichen. Angenommen, Sie haben einen Sandy Bridge / Ivy Bridge Prozessor. Der Peak SP FLOPS ist
%Vor%Für doppelte Präzession durch zwei dividieren. Wenn Sie Haswell und MKL FMA verwenden, dann verdoppeln Sie den Peak FLOPS. Um die Frequenz richtig einzustellen, müssen Sie die Turbo-Boost-Werte für alle Kerne verwenden (es ist niedriger als für einen einzelnen Kern). Sie können dies überprüfen, wenn Sie Ihr System nicht übertaktet haben oder CPU-Z unter Windows oder Powertop unter Linux verwenden, wenn Sie ein übertaktetes System haben.
Verwenden Sie eine temporäre Matrix, um A '* Y zu berechnen, aber stellen Sie sicher, dass Sie eigen sagen, dass kein Aliasing stattfindet: temp.noalias() = A.adjoint()*Y . Dann berechne dein Ergebnis und sage noch einmal eigen, dass Objekte keinen Alias haben: result.noalias() = A*temp .
Es würde nur dann eine redundante Berechnung geben, wenn Sie (A*A')*Y ausführen würden, da in diesem Fall (A*A') symmetrisch ist und nur die Hälfte der Berechnung erforderlich ist. Wie Sie jedoch festgestellt haben, ist die Ausführung von A*(A'*Y) immer noch viel schneller. In diesem Fall gibt es keine redundanten Berechnungen. Ich bestätige, dass die Kosten für die temporäre Erstellung hier vernachlässigbar sind.
Ich denke, das folgende ausführen
%Vor%ist dasselbe wie das
%Vor% Wenn Ihre Matrix Y in den Cache passt, können Sie wahrscheinlich davon profitieren
oder, wenn die vorherige Schreibweise nicht erlaubt ist, so
%Vor%Dann müssen Sie nicht die Adjungierte der Matrix A tun, und speichern Sie die temporäre adjungierte Matrix für A, die viel teurer sein wird als die für Y.
Wenn Ihre Matrix Y in den Cache passt, sollte es viel schneller sein, eine Schleife über die Spalten von A für die erste Multiplikation auszuführen und dann über die Zeilen von A für die zweite Multimimplikation (mit Y.adjoint () in der Cache für die erste Multiplikation und temp.adjoint () für die zweite), aber ich denke, dass intern Eigen kümmert sich schon um diese Dinge.
Tags und Links optimization c++ eigen x86-64 matrix-multiplication