schnellste lineare Systemlösung für kleine quadratische Matrizen (10x10)

Ich bin sehr daran interessiert, die lineare Systemlösung für kleine Matrizen (10x10), manchmal auch winzige Matrizen genannt, zu optimieren . Gibt es dafür eine fertige Lösung? Die Matrix kann als nicht singulär angenommen werden.

Dieser Solver soll mehr als 1 000 000 Mal in Mikrosekunden auf einer Intel-CPU ausgeführt werden. Ich spreche von der Optimierungsstufe, die in Computerspielen verwendet wird. Egal, ob ich es in Assembly- und Architektur-spezifischen Codes codiere oder die Reduzierung von Präzisions- oder Zuverlässigkeits-Kompromissen untersuche und Gleitkomma-Hacks verwende (ich verwende das Kompilierungsflag -ffast-math, kein Problem). Die Lösung kann sogar in etwa 20% der Fälle fehlschlagen!

Eigens PartialPivLu ist das schnellste in meinem aktuellen Benchmark und übertrifft LAPACK, wenn es mit -O3 und einem guten Compiler optimiert wird. Aber jetzt bin ich gerade dabei, einen benutzerdefinierten linearen Löser von Hand herzustellen. Jeder Rat wäre sehr dankbar. Ich werde meine Lösung Open Source machen und wichtige Erkenntnisse in Veröffentlichungen usw. kennenlernen.

Verwandte: Geschwindigkeit der Lösung des linearen Systems mit Blockdiagonalmatrix Was ist die schnellste Methode, um Millionen von Matrizen zu invertieren? https://stackoverflow.com/q/50909385/1489510

Wenn Sie einen Eigenmatrixtyp verwenden, bei dem die Anzahl der Zeilen und Spalten zur Kompilierungszeit in den Typ codiert wird, erhalten Sie einen Vorteil gegenüber LAPACK, bei dem die Matrixgröße nur zur Laufzeit bekannt ist. Diese zusätzlichen Informationen ermöglichen es dem Compiler, das vollständige oder teilweise Abrollen der Schleife durchzuführen, wodurch viele Verzweigungsanweisungen entfallen. Wenn Sie eine vorhandene Bibliothek verwenden möchten, anstatt Ihre eigenen Kernel zu schreiben, ist es wahrscheinlich wichtig, einen Datentyp zu haben, bei dem die Matrixgröße als C ++ - Vorlagenparameter angegeben werden kann. Die einzige andere Bibliothek, von der ich weiß, dass sie dies tut, ist Blaze . Es könnte sich also lohnen, sie mit Eigen zu vergleichen.

Wenn Sie sich für eine eigene Implementierung entscheiden, ist das, was PETSc für das Block-CSR-Format tut, möglicherweise ein nützliches Beispiel, obwohl PETSc selbst wahrscheinlich nicht das richtige Werkzeug für Ihre Vorstellungen ist. Anstatt eine Schleife zu schreiben, schreiben sie jede einzelne Operation für kleine Matrixvektormultiplikationen explizit aus (siehe diese Datei in ihrem Repository). Dies garantiert, dass es keine Verzweigungsanweisungen gibt, wie Sie sie mit einer Schleife erhalten könnten. Die Versionen des Codes mit AVX-Anweisungen sind ein gutes Beispiel für die tatsächliche Verwendung von Vektorerweiterungen. Zum Beispiel dieser Funktion verwendet die__m256dDatentyp für den gleichzeitigen Betrieb von vier Doppelgeräten gleichzeitig. Sie können eine spürbare Leistungssteigerung erzielen, indem Sie alle Operationen explizit mit Vektorerweiterungen ausschreiben, nur für die LU-Faktorisierung anstelle der Matrix-Vektor-Multiplikation. Anstatt den C-Code tatsächlich von Hand zu schreiben, sollten Sie ein Skript verwenden, um ihn zu generieren. Es könnte auch Spaß machen zu sehen, ob es einen nennenswerten Leistungsunterschied gibt, wenn Sie einige der Vorgänge neu anordnen, um das Pipelining von Anweisungen besser nutzen zu können.

Sie können auch einige Kilometer mit dem Tool STOKE sammeln , das zufällig den Raum möglicher Programmtransformationen untersucht, um eine schnellere Version zu finden.

Eine andere Idee könnte sein, einen generativen Ansatz zu verwenden (ein Programm, das ein Programm schreibt). Verfassen Sie ein (Meta-) Programm, das die Folge von C / C ++ - Anweisungen ausspuckt, um eine nicht drehbare ** LU auf einem 10x10-System auszuführen. Nehmen Sie im Grunde genommen das k / i / j-Schleifennest und reduzieren Sie es in etwa O (1000) Zeilen von Skalararithmetik. Führen Sie dann das generierte Programm in den optimierenden Compiler ein. Was ich hier für interessant halte, ist das Entfernen der Schleifen, um jede Datenabhängigkeit und jeden redundanten Unterausdruck freizulegen, und gibt dem Compiler die maximale Möglichkeit, Anweisungen neu zu ordnen, damit sie der tatsächlichen Hardware gut zugeordnet werden können (z. B. Anzahl der Ausführungseinheiten, Gefahren / Verzögerungen usw.) auf).

Wenn Sie alle Matrizen (oder nur einige davon) kennen, können Sie den Durchsatz verbessern, indem Sie anstelle von Skalarcode SIMD-Intrinsics / -Funktionen (SSE / AVX) aufrufen. Hier würden Sie die peinliche Parallelität zwischen den Instanzen ausnutzen, anstatt jede Parallelität innerhalb einer einzelnen Instanz zu verfolgen. Zum Beispiel könnten Sie 4 LUs mit doppelter Genauigkeit gleichzeitig mit AVX256-Intrinsics ausführen, indem Sie 4 Matrizen "quer" über das Register packen und dieselben Operationen ** für alle ausführen.

** Daher der Fokus auf nicht schwenkbare LU. Das Schwenken verdirbt diesen Ansatz auf zwei Arten. Erstens werden Zweige aufgrund der Pivot-Auswahl eingeführt, was bedeutet, dass Ihre Datenabhängigkeiten nicht so genau bekannt sind. Zweitens bedeutet dies, dass verschiedene SIMD- "Slots" unterschiedliche Aufgaben ausführen müssten, da Instanz A möglicherweise anders als Instanz B schwenkt. Wenn Sie also etwas davon verfolgen, würde ich vorschlagen, Ihre Matrizen vor der Berechnung statisch zu schwenken (permutieren Sie den größten Eintrag) jeder Spalte zur Diagonale).

Ihre Frage führt zu zwei unterschiedlichen Überlegungen.

Zunächst müssen Sie den richtigen Algorithmus auswählen. Daher sollte die Frage berücksichtigt werden, ob die Matrizen eine Struktur haben. Wenn beispielsweise die Matrizen symmetrisch sind, ist eine Cholesky-Zerlegung effizienter als LU. Wenn Sie nur eine begrenzte Genauigkeit benötigen, kann eine iterative Methode schneller sein.

$10 \times 10$

Insgesamt hängt die Antwort auf Ihre Frage stark von der Hardware und den Matrizen ab, die Sie berücksichtigen. Es gibt wahrscheinlich keine eindeutige Antwort und Sie müssen einige Dinge ausprobieren, um eine optimale Methode zu finden.

Ich würde blockweise Inversion versuchen.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen verwendet eine optimierte Routine, um die Inverse einer 4x4-Matrix zu berechnen. Dies ist wahrscheinlich die beste, die Sie erhalten werden. Versuchen Sie das so oft wie möglich zu verwenden.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Oben links: 8x8. Oben rechts: 8x2. Unten links: 2x8. Unten rechts: 2x2. Invertieren Sie das 8x8 mit dem optimierten 4x4-Inversionscode. Der Rest sind Matrixprodukte.

BEARBEITEN: Die Verwendung von 6x6-, 6x4-, 4x6- und 4x4-Blöcken hat sich als etwas schneller erwiesen als oben beschrieben.

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Hier sind die Ergebnisse eines Benchmark-Laufs mit einer Million Eigen::Matrix<double,10,10>::Random()Matrizen und Eigen::Matrix<double,10,1>::Random()Vektoren. Bei all meinen Tests ist meine Umkehrung immer schneller. Meine Lösungsroutine besteht darin, das Inverse zu berechnen und es dann mit einem Vektor zu multiplizieren. Manchmal ist es schneller als Eigen, manchmal nicht. Meine Benchmarking-Methode ist möglicherweise fehlerhaft (Turbo-Boost usw. wurde nicht deaktiviert). Außerdem sind Eigens Zufallsfunktionen möglicherweise nicht repräsentativ für reale Daten.

• Eigener partieller Pivot invers: 3036 Millisekunden
• Meine Umkehrung mit 8x8 oberem Block: 1638 Millisekunden
• Meine Umkehrung mit 6x6 oberem Block: 1234 Millisekunden
• Eigene partielle Pivot-Lösung: 1791 Millisekunden
• Meine Lösung mit 8x8 oberen Block: 1739 Millisekunden
• Meine Lösung mit 6x6 oberen Block: 1286 Millisekunden

Ich bin sehr interessiert zu sehen, ob jemand dies weiter optimieren kann, da ich eine Finite-Elemente-Anwendung habe, die eine Unmenge von 10x10-Matrizen invertiert (und ja, ich benötige einzelne Koeffizienten der Inversen, so dass das direkte Lösen eines linearen Systems nicht immer eine Option ist). .