Eigen- und C++ 11-Typ-Inferenz scheitert für Cholesky des Matrixprodukts

Question

Ich versuche, die cholesky-Zerlegung des Produkts einer Matrix mit seiner Transponierung zu nehmen, unter Verwendung von Eigen und C++ 11 "Auto" -Typ. Das Problem kommt, wenn ich versuche zu tun

auto c = a * b 
auto cTc = c.tranpose() * c; 
auto chol = cTc.llt(); 

Ich verwende XCode 6.1, Eigen 3.2.2. Der Typfehler, den ich bekomme, ist here.

Dieses minimale Beispiel zeigt das Problem an meinem Computer. Ändern Sie den Typ c von auto zu MatrixXd, damit es funktioniert.

#include <iostream> 
#include <Eigen/Eigen> 
using namespace std; 
using namespace Eigen; 


int main(int argc, const char * argv[]) { 
    MatrixXd a = MatrixXd::Random(100, 3); 
    MatrixXd b = MatrixXd::Random(3, 100); 
    auto c = a * b; 
    auto cTc = c.transpose() * c; 
    auto chol = cTc.llt(); 
    return 0; 
} 

Gibt es eine Möglichkeit, dies zu tun, während immer noch Auto verwenden?

Als eine Nebenfrage, gibt es einen Leistungsgrund, die Matrix nicht zu bestätigen ist ein MatrixXd in jeder Phase? Die Verwendung von Auto würde es Eigen erlauben, die Antwort als irgendeinen seltsamen Template-Ausdruck zu behalten, den es sich vorstellt. Ich bin mir nicht sicher, ob das Tippen als MatrixXd Probleme verursachen würde oder nicht.

Answer 1

Lassen Sie mich zusammenfassen, was passiert und warum es falsch ist. Zunächst einmal wollen wir die auto Schlüsselwörter mit den Typen instanziiert sie einnehmen:

typedef GeneralProduct<MatrixXd,MatrixXd> Prod; 
Prod c = a * b; 
GeneralProduct<Transpose<Prod>,Prod> cTc = c.transpose() * c; 

Beachten Sie, dass Eigen Ausdruck Template-Bibliothek ist.Hier ist GeneralProduct<> ein abstrakter Typ, der das Produkt darstellt. Keine Berechnung wird ausgeführt. Deshalb, wenn Sie cTc auf eine MatrixXd als Kopie:

MatrixXd d = cTc; 

dem entspricht:

MatrixXd d = c.transpose() * c; 

dann das Produkt a*b wird zweimal durchgeführt! Also in jedem Fall ist es sehr bevorzugt, a*b innerhalb einer expliziten vorübergehend zu bewerten, und gleich für c^T*c:

MatrixXd c = a * b; 
MatrixXd cTc = c.transpose() * c; 

Die letzte Zeile:

auto chol = cTc.llt(); 

ist auch ziemlich falsch. Wenn cTc ein abstrakter Produkttyp ist, versucht es eine Cholesky-Faktorisierung zu instanziieren, die an einem abstrakten Produkttyp arbeitet, was nicht möglich ist. Nun, wenn cTc ein MatrixXd ist, dann sollten Sie den Code arbeiten, aber dies noch nicht die bevorzugte Art und Weise wie das Verfahren llt() ist eher Einzeiler Ausdruck zu implementieren wie:

VectorXd b = ...; 
VectorXd x = cTc.llt().solve(b); 

Wenn Sie ein benannte Cholesky Objekt mögen, dann verwenden lieber den Konstruktor:

LLT<MatrixXd> chol(cTc); 

oder sogar:

LLT chol (c.transpose() * c);

was äquivalent ist, es sei denn, Sie müssen c.transpose() * c in anderen Berechnungen verwenden.

Schließlich abhängig von der Größe der a und b, könnte es vorzuziehen sein zu berechnen cTc als:

MatrixXd cTc = b.transpose() * (a.transpose() * a) * b; 

In Zukunft (dh Eigen 3.3), wird Eigen Lage sein zu sehen:

auto c = a * b; 
MatrixXd cTc = c.transpose() * c; 

als ein Produkt aus vier Matrizen m0.transpose() * m1.transpose() * m2 * m3 und setzen Sie die Klammer an der richtigen Stelle. Allerdings kann es nicht wissen, m0==m3 und m1==m2, und daher, wenn der bevorzugte Weg ist, a*b in einem vorläufigen auszuwerten, dann müssen Sie es immer noch tun.

Answer 2

Das Problem ist, dass die erste Multiplikation eine Eigen::GeneralProduct anstelle einer MatrixXd zurückgibt und auto den Rückgabetyp abholt. Sie können implizit einen MatrixXd von einem Eigen::GeneralProduct erstellen, also wenn Sie den Typ explizit deklarieren, funktioniert es richtig. Siehe http://eigen.tuxfamily.org/dox/classEigen_1_1GeneralProduct.html

EDIT: Ich bin kein Experte auf das Eigenprodukt oder Leistungsmerkmale des Castings. Ich habe nur die Antwort aus der Fehlermeldung vermutet und aus der Online-Dokumentation bestätigt. Profiling ist immer die beste Option, um die Leistung verschiedener Teile Ihres Codes zu überprüfen.

Answer 3

Ich bin kein Experte bei Eigen, aber Bibliotheken wie diese geben häufig Proxy-Objekte aus Operationen zurück und verwenden dann implizite Konvertierung oder Konstruktoren, um die eigentliche Arbeit zu erzwingen. (Ausdruckvorlagen sind ein extremes Beispiel dafür.) Dies vermeidet unnötiges Kopieren der vollständigen Matrix von Daten in vielen Situationen. Unglücklicherweise ist ziemlich glücklich, einfach ein Objekt des Proxy-Typs zu erstellen, das normalerweise Benutzer der Bibliothek niemals explizit deklarieren würden. Da Sie letztlich die Zahlen berechnen müssen, gibt es keinen performativen Treffer von Casting auf MatrixXd. (Scott Meyers, in "Effective Moderne C++", gibt das zugehörige Beispiel für die Verwendung auto mit vector<bool>, wo operator[](size_t i) einen Proxy zurückgibt.)

Answer 4

NICHT auto mit Eigen Ausdrücke verwenden DO. Ich stieß auf noch mehr „dramatisch“ Probleme mit diesem vor, siehe

eigen auto type deduction in general product

und wurde von einem des Eigen Schöpfer (Gael) nicht verwenden auto mit Eigen Ausdrücken beraten.

Die Umwandlung von einem Ausdruck in einen bestimmten Typ wie MatrixXd sollte extrem schnell sein, es sei denn, Sie möchten eine faule Auswertung (da bei der Besetzung der Ergebnis ausgewertet wird).