Super C ++ Optimierung der Matrixmultiplikation mit Armadillo

Ich verwende Armadillo, um sehr intensive Matrixmultiplikationen mit Seitenlängen von , wobei bis zu 20 oder sogar mehr kann. Ich verwende Armadillo mit OpenBLAS für die Matrixmultiplikation, was in parallelen Kernen sehr gute Arbeit zu leisten scheint, außer dass ich ein Problem mit dem Formalismus der Multiplikation in Armadillo zur Superoptimierung der Leistung habe.$2^n$$n$

Angenommen, ich habe eine Schleife in der folgenden Form:

arma::cx_mat stateMatrix, evolutionMatrix; //armadillo complex matrix type
for(double t = t0; t < t1; t += 1/sampleRate)
{
    ...
    stateMatrix = evolutionMatrix*stateMatrix;
    ...
}

In grundlegendem C ++ stelle ich hier das Problem fest, dass C ++ ein neues Objekt cx_matzum Speichern zuweist evolutionMatrix*stateMatrixund dann das neue Objekt stateMatrixmit kopiert operator=(). Das ist sehr, sehr ineffizient. Es ist bekannt, dass die Rückgabe komplexer Klassen großer Datentypen eine schlechte Idee ist, oder?

Ich sehe dies effizienter, wenn eine Funktion die Multiplikation in der folgenden Form ausführt:

void multiply(const cx_mat& mat1, const cx_mat& mat2, cx_mat& output)
{
    ... //multiplication of mat1 and mat2 and then store it in output
}

Auf diese Weise müssen keine großen Objekte mit einem Rückgabewert kopiert werden, und die Ausgabe muss nicht bei jeder Multiplikation neu zugewiesen werden.

Die Frage : Wie kann ich einen Kompromiss finden, bei dem ich Armadillo für die Multiplikation mit seiner schönen BLAS-Oberfläche verwenden kann, und dies effizient, ohne Matrixobjekte neu erstellen und bei jeder Operation kopieren zu müssen?

Ist das nicht ein Implementierungsproblem in Armadillo?

In grundlegendem C ++ besteht das Problem darin, dass C ++ ein neues Objekt von cx_mat zum Speichern von evolutionMatrix * stateMatrix zuweist und dann das neue Objekt mit operator = () nach stateMatrix kopiert.

Ich denke, Sie haben Recht, dass es temporäre Elemente erstellt, was zu langsam ist, aber ich denke, der Grund, warum es das tut, ist falsch.

Armadillo verwendet wie jede gute lineare C ++ - Algebra-Bibliothek Ausdrucksvorlagen, um eine verzögerte Auswertung von Ausdrücken zu implementieren. Wenn Sie ein Matrixprodukt aufschreiben wie A*B, werden keine Provisorien erstellt, macht stattdessen Armadillo ein temporäres Objekt ( x) , die Verweise auf hält Aund B, und dann wie ein Ausdruck gegeben später C = xberechnet das Matrixprodukt speichert das Ergebnis direkt in Cohne Schaffung Provisorien.

Diese Optimierung wird auch verwendet, um Ausdrücke zu verarbeiten A*B*C*D, bei denen je nach Matrixgröße bestimmte Multiplikationsordnungen effizienter sind als andere.

Ist das nicht ein Implementierungsproblem in Armadillo?

Wenn Armadillo diese Optimierung nicht durchführt, wäre dies ein Fehler in Armadillo, der den Entwicklern gemeldet werden sollte.

In Ihrem Fall gibt es jedoch ein anderes Problem, das wichtiger ist. In einem Ausdruck wie A=B*Cdem Speichern von Aenthält keine Eingabedaten, wenn Akein Alias Boder C. In Ihrem Fall würde das A = A*BSchreiben von etwas in die Ausgabematrix auch eine der Eingabematrizen ändern.

Auch angesichts Ihrer vorgeschlagenen Funktion

multiply(const cx_mat&, const cx_mat&, cx_mat&)

Wie genau würde diese Funktion im Ausdruck helfen multiply(A, B, A)? Bei den meisten normalen Implementierungen dieser Funktion würde dies zu einem Fehler führen. Es müsste selbst einen temporären Speicher verwenden, um sicherzustellen, dass die Eingabedaten nicht beschädigt werden. Ihr Vorschlag ist ziemlich genau, wie Armadillo die Matrixmultiplikation bereits implementiert, aber ich denke, es ist wahrscheinlich wichtig, Situationen zu vermeiden, wie multiply(A, B, A)durch die Zuweisung einer temporären.

Die wahrscheinlichste Erklärung dafür, warum Armadillo diese Optimierung nicht durchführt , ist, dass es falsch wäre, dies zu tun.

Schließlich gibt es eine viel einfachere Möglichkeit, das zu tun, was Sie wollen:

cx_mat *A, *Atemp, B;
for (;;) {
  *Atemp = (*A)*B;
  swap(A, Atemp);
}

Dies ist identisch mit

cx_mat A, B;
for (;;) {
  A = A*B;
}

Es wird jedoch eine temporäre Matrix anstelle einer temporären Matrix pro Iteration zugewiesen.

@BillGreene verweist auf die "Rückgabewertoptimierung" als einen Weg, um das grundlegende Problem zu umgehen, aber dies hilft tatsächlich nur für die Hälfte davon. Angenommen, Sie haben Code dieses Formulars:

struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); };

ExpensiveObject operator+ (ExpensiveObject &obj1,
                           ExpensiveObject &obj2)
{
   ExpensiveObject tmp;
   ...compute tmp based on obj1 and obj2...
   return tmp;
}

void f() {
  ExpensiveObject o1, o2, o3;
  ...initialize o1, o2...;
  o3 = o1 + o2;
}

Ein naiver Compiler wird

1. Erstellen Sie einen Slot, um das Ergebnis der Plus-Operation zu speichern (eine temporäre).
2. Betreiber anrufen +,
3. Erstellen Sie das 'tmp'-Objekt in Operator + (eine zweite temporäre).
4. tmp berechnen,
5. Kopieren Sie tmp in den Ergebnis-Slot.
6. zerstöre tmp,
7. Kopieren Sie das Ergebnisobjekt in o3
8. Zerstöre das Ergebnisobjekt

Die Rückgabewertoptimierung kann nur das 'tmp'-Objekt und den' result'-Slot vereinheitlichen, jedoch nicht die Notwendigkeit einer Kopie beseitigen. Sie müssen also noch eine temporäre Datei erstellen, den Kopiervorgang ausführen und eine temporäre Version zerstören.

Der einzige Weg, dies zu umgehen, ist, dass Operator + kein Objekt zurückgibt, sondern ein Objekt einer Zwischenklasse, das, wenn es einem zugewiesen wird ExpensiveObject, die Additions- und Kopieroperation an Ort und Stelle ausführt. Dies ist der typische Ansatz, der in Ausdrucksvorlagenbibliotheken verwendet wird .

Stackoverflow ( https://stackoverflow.com/ ) ist wahrscheinlich ein besseres Diskussionsforum für diese Frage. Hier ist jedoch eine kurze Antwort.

Ich bezweifle, dass der C ++ - Compiler Code für diesen Ausdruck generiert, wie Sie oben beschrieben haben. Alle modernen C ++ - Compiler implementieren eine Optimierung namens "Rückgabewertoptimierung" ( http://en.wikipedia.org/wiki/Return_value_optimization ). Bei der Rückgabewertoptimierung wird das Ergebnis von evolutionMatrix*stateMatrixdirekt in gespeichert stateMatrix; Es wird keine Kopie erstellt.

Es gibt offensichtlich erhebliche Verwirrung zu diesem Thema und das ist einer der Gründe, warum ich vorgeschlagen habe, dass Stackoverflow ein besseres Forum sein könnte. Es gibt dort viele C ++ "Sprachanwälte", während die meisten von uns hier lieber ihre Zeit mit CSE verbringen würden. ;-);

Ich habe das folgende einfache Beispiel basierend auf Professor Bangerths Beitrag erstellt:

#ifndef NDEBUG
#include <iostream>

using namespace std;
#endif

class ExpensiveObject  {
public:
  ExpensiveObject () {
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor called." << endl;
#endif
    v = 0;
  }
  ExpensiveObject (int i) { 
#ifndef NDEBUG
    cout << "ExpensiveObject  constructor(int) called." << endl;
#endif
    v = i; 
  }
  ExpensiveObject (const ExpensiveObject  &a) {
    v = a.v;
#ifndef NDEBUG
    cout << "ExpensiveObject  copy constructor called." << endl;
#endif
  }
  ~ExpensiveObject() {
#ifndef NDEBUG
    cout << "ExpensiveObject  destructor called." << endl;
#endif
  }
  ExpensiveObject  operator=(const ExpensiveObject  &a) {
#ifndef NDEBUG
    cout << "ExpensiveObject  assignment operator called." << endl;
#endif
    if (this != &a) {
      return ExpensiveObject (a);
    }
  }
  void print() const {
#ifndef NDEBUG
    cout << "v=" << v << endl;
#endif
  }
  int getV() const {
    return v;
  }
private:
  int v;
};

ExpensiveObject  operator+(const ExpensiveObject  &a1, const ExpensiveObject  &a2) {
#ifndef NDEBUG
  cout << "ExpensiveObject  operator+ called." << endl;
#endif
  return ExpensiveObject (a1.getV() + a2.getV());
}

int main()
{
  ExpensiveObject  a(2), b(3);
  ExpensiveObject  c = a + b;
#ifndef NDEBUG
  c.print();
#endif
}

Es sieht komplizierter aus als es tatsächlich ist, weil ich beim Kompilieren im optimierten Modus den gesamten Code für die Druckausgabe vollständig entfernen wollte. Wenn ich die mit einer Debug-Option kompilierte Version ausführe, erhalte ich die folgende Ausgabe:

ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
v=5
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.

Das erste, was zu bemerken ist, ist, dass keine Provisorien konstruiert werden - nur a, b und c. Der Standardkonstruktor und der Zuweisungsoperator werden niemals aufgerufen, da sie in diesem Beispiel nicht benötigt werden.

Professor Bangerth erwähnte Ausdrucksvorlagen. In der Tat ist diese Optimierungstechnik sehr wichtig, um eine gute Leistung in einer Matrixklassenbibliothek zu erzielen. Dies ist jedoch nur dann wichtig, wenn die Objektausdrücke komplizierter sind als nur a + b. Wenn zum Beispiel mein Test stattdessen war:

  ExpensiveObject  a(2), b(3), c(9);
  ExpensiveObject  d = a + b + c;

Ich würde die folgende Ausgabe erhalten:

ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  operator+ called.
 ExpensiveObject  constructor(int) called.
 ExpensiveObject  destructor called.
 v=14
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.
 ExpensiveObject  destructor called.

Dieser Fall zeigt den unerwünschten Aufbau eines temporären (5 Aufrufe an den Konstruktor und zwei Aufrufe von Operator +). Die ordnungsgemäße Verwendung von Ausdrucksvorlagen (ein Thema, das weit über den Rahmen dieses Forums hinausgeht) würde dies vorübergehend verhindern. (Für die Hochmotivierten finden Sie eine besonders lesbare Diskussion der Ausdrucksvorlagen in Kapitel 18 von http://www.amazon.com/C-Templates-The-Complete-Guide/dp/0201734842 ).

Der eigentliche "Beweis" dafür, was der Compiler tatsächlich tut, besteht darin, den vom Compiler ausgegebenen Assembler-Code zu untersuchen. Für das erste Beispiel ist dieser Code im kompilierten Modus erstaunlich einfach. Alle Funktionsaufrufe wurden entfernt und der Assembler-Code lädt im Wesentlichen 2 in ein Register, 3 in ein zweites und fügt sie hinzu.

$\mathcal{O}(n^{2.8})$$\mathcal{O}(n^2)$$n$

Das heißt, es sei denn , Sie eine enorme Konstante in das Kopieren entstehen - das ist eigentlich nicht so weit hergeholt ist, da die Version mit dem Kopieren ist sehr viel teurer in anderer Hinsicht: es viel mehr Speicher benötigt. Wenn Sie also am Ende von und zur Festplatte wechseln müssen, kann das Kopieren tatsächlich zum Engpass werden. Aber selbst wenn Sie nicht kopieren alles selbst, ein stark parallelisierten Algorithmus kann auch einige Kopien seiner eigenen machen. Der einzige Weg, um sicherzustellen, dass nicht zu viel Speicher in jedem Schritt verwendet wird, besteht darin, die Multiplikation in Spalten vonstateMatrix aufzuteilen , sodass jeweils nur kleine Multiplikationen durchgeführt werden. Zum Beispiel können Sie definieren

class HChunkMatrix // optimised for destructive left-multiplication updates
{
  std::vector<arma::cx_mat> colChunks; // e.g. for an m×n matrix,
                                      //  use √n chunks, each an m×√n matrix
 public:
  ...

  HChunkMatrix& operator *= (const arma::cx_mat& lhMult) {
    for (&colChunk: colChunks) {
      colChunk = lhMult * colChunk;
    }
    return *this;
  }
}

Sie sollten auch überlegen, ob Sie das überhaupt erst entwickeln müssen stateMatrix. Wenn Sie im Grunde nur eine unabhängige zeitliche Entwicklung von nStaatskets wünschen , können Sie sie genauso gut einzeln weiterentwickeln, was viel weniger speicherintensiv ist. Insbesondere wenn evolutionMatrixes spärlich ist , sollten Sie unbedingt überprüfen! Denn das ist im Grunde nur ein Hamiltonianer, nicht wahr? Hamiltonianer sind oft spärlich oder annähernd spärlich.

_{$\mathcal{O}(n^{2.38})$}

Modernes C ++ bietet eine Lösung für das Problem, indem "Verschiebungskonstruktoren" und "rWertreferenzen" verwendet werden.

Ein "Verschiebungskonstruktor" ist ein Konstruktor für eine Klasse, beispielsweise eine Matrixklasse, die eine andere Instanz derselben Klasse verwendet und die Daten von der anderen Instanz in die neue Instanz verschiebt, wobei die ursprüngliche Instanz leer bleibt. Normalerweise hat ein Matrixobjekt zwei Zahlen für die Größe und einen Zeiger auf die Daten. Wenn ein normaler Konstruktor die Daten duplizieren würde, kopiert ein Verschiebungskonstruktor nur die beiden Zahlen und den Zeiger. Dies ist also sehr schnell.

Für temporäre Variablen wird eine "rWertreferenz" verwendet, die beispielsweise als "Matrix &&" anstelle der üblichen "Matrix &" geschrieben wird. Sie würden eine Matrixmultiplikation als Rückgabe einer Matrix && deklarieren. Auf diese Weise stellt der Compiler sicher, dass ein sehr billiger Verschiebungskonstruktor verwendet wird, um das Ergebnis aus der Funktion herauszuholen, die ihn aufruft. Ein Ausdruck wie result = (a + b) * (c + d), bei dem a, b, c, d riesige Matrixobjekte sind, wird also ohne Kopieren ausgeführt.

Wenn Sie nach "rWertreferenzen und Verschiebungskonstruktoren" googeln, finden Sie Beispiele und Tutorials.

$v'*M*M*M*M*M$$M*M*M*M*M*v$

Andererseits habe OpenBLAS eine größere Sammlung architekturspezifischer Optimierungen, sodass Eigen für Sie möglicherweise ein Gewinn ist oder nicht. Leider gibt es keine Bibliothek für lineare Algebra, die so großartig ist, dass Sie nicht einmal die anderen berücksichtigen müssen, wenn Sie um "die letzten 10%" der Leistung kämpfen. Wrapper sind keine 100% ige Lösung. Die meisten (alle?) von ihnen können die Fähigkeit von Eigen nicht nutzen, Berechnungen auf diese Weise zusammenzuführen.