boost::numeric::ublas Matrizenrechnung - ist das wirklich so Benutzerunfreundlich?



  • Hallo zusammen,

    ich versuche gerade ein paar Matrizenoperationen durchzufüren (ich portiere code von Matlab nach C++).
    Ich habe mich fürs erste für boost::numeric::ublas entschieden, da boost ja meistens gut ist und ich boost sowieso schon hatte.

    Allerdings stören mich einige sachen:
    - Operatorüberladung ist nur teilweise vorhanden, für multiplikationen muss ich idr. prod(A,B) statt A*B schreiben.
    - Verschachtelte operationen von produkten führen zu kompilierungsfehlern (static_assert schlägt fehl)
    - Es gibt nichts zum invertieren einer Matrix

    Ich will folgendes implementieren:

    % This is the code which implements the discrete Kalman filter:
    
       % Prediction for state vector and covariance:
       s.x = s.A*s.x + s.B*s.u;
       s.P = s.A * s.P * s.A' + s.Q;
    
       % Compute Kalman gain factor:
       K = s.P*s.H'*inv(s.H*s.P*s.H'+s.R);
    
       % Correction based on observation:
       s.x = s.x + K*(s.z-s.H*s.x);
       s.P = s.P - K*s.H*s.P;
    

    alle großbuchstaben (P, K usw) sind matrizen, die Kleinbuchstaben sind vektoren.

    Wenn ich jetzt z.B. die erste berechnung übersetze und folgendes versuche:

    x = A * x + B * u;
    

    Dann fehlen sofort die Operatorüberladungen.

    Mit

    x = prod(A,x) + prod(B,u);
    

    Funktionierts dann sogar.

    Zweite zeile:

    P = A * P * trans(A) + Q;
    

    Fehlen natürlich die operatorüpberladungen wieder.
    Zweiter versuch:

    P = prod(prod(A, P), trans(A)) + Q;
    

    ->

    error C2338: E1::complexity==0 && E2::complexity == 0

    Das wurde übrigens durch ein BOOST_STATIC_ASSERT ausgelöst.

    Auch hier habe ich eine Lösung gefunden:

    P = prod(A, P);
    P = prod(P, trans(A));
    P += Q;
    

    Und dann bei der nächsten operation fehlt die inverse...
    Dazu habe ich das hier im Internet gefunden:

    template<class T>
    bool InvertMatrix(const matrix_expression<T>& input, matrix_expression<T>& inverse)
    {
    	typedef permutation_matrix<std::size_t> pmatrix;
    	// create a working copy of the input
    	matrix<T> A(input);
    	// create a permutation matrix for the LU-factorization
    	pmatrix pm(A.size1());
    	// perform LU-factorization
    	int res = lu_factorize(A, pm);
    	if (res != 0)
    		return false;
    
    	// create identity matrix of "inverse"
    	inverse.assign(identity_matrix<T> (A.size1()));
    	// backsubstitute to get the inverse
    	lu_substitute(A, pm, inverse);
    	return true;
    }
    

    Auch bei der verwendung davon muss ich die operationen immer schön in einzelne anweisungen stückeln damit das klappt.
    Das ganze wird dadurch extremst unübersichtlich und aufwändig zu schreiben.

    Mach ich irgendwas falsch, oder ist die Bibliothek wirklich so Benutzerunfreundlich?
    Gibt es vielleicht bessere Alternativen als boost::numeric::ublas?





  • Danke!
    Das ist in der Tat viel einfacher zu verwenden und das einbinden ist auch super einfach, da header only 🙂



  • Q schrieb:

    K = s.P*s.H'*inv(s.H*s.P*s.H'+s.R);
    

    Nie inv verwenden, wenn Du die Inverse gar nicht brauchst!

    Statt A*inv(B) kannst und solltest Du A/B schreiben.
    Statt inv(A)*B kannst und solltest Du A\B schreiben.

    Also in Deinem Fall

    K = s.P*s.H'/(s.H*s.P*s.H'+s.R);
    

    Boost.uBLAS habe ich auch hier in Benutzung. Für die nächsten Projekte würde ich aber Eigen ausprobieren wollen. das u ("micro") in uBLAS ist ernst zu nehmen.



  • Den Matlab-code mit den inversen hab ich nicht selbst geschrieben, sondern von wem anders erhalten.
    Trotzdem danke für den Tipp.
    ublas werd ich wohl nie wieder verwenden 😃 einfach grauenhaft im vergleich mit Eigen.



  • Nur so als kleiner Tipp am Rande da du gemäß deiner Gleichung offenbar ähnliche Sachen planst wie ich.

    Ich habe vor kurzem ublas gegen matlab und andere Bibliotheken gebencht.

    Boost wird mit größeren Matrizen unbrauchbar. Bsp. : bei einer 2000 x 2000 Matrix hast du zwischen MATLAB und ublas einen Faktor von 60x.

    Wenn du eine wirklich elegante und kostenlose Lösung suchst dann nimm Eigen: http://eigen.tuxfamily.org/index.php?title=Main_Page

    Scheue dich nicht vor GPL3+. Da es nur eine eine Template Header Implementierung ist entfallen viele Dinge.

    Wenn du Geld hast und Performance sehr wichtig ist dann nimm den Intel Compiler und MKL. Ist zwar sehr hässlich und unbequem .. XGEMM() like aber schneller wirst du es in C++ dann nicht hinkriegen.



  • Mit dem Faktor war die Zeit zur Bildung eines Produktes prod(A,B) gemeint.
    MATLAB:Eigen ist übrigens 1:3 ... also fast perfekt.

    Und nein ich bin kein Anfänger und habe es nicht mit DEBUG oder ähnlichem gebencht. 🙄



  • wie meisnt du das mit dem faktor 1:3? wer von beiden ist denn jetzt schneller?



  • Q schrieb:

    wie meisnt du das mit dem faktor 1:3? wer von beiden ist denn jetzt schneller?

    MATLAB ist definitiv schneller als Eigen. Nicht nur bei Produkten. Von MEX Kompilaten diverser Funktionen ganz zu schweigen.

    Alles andere wäre ein Wunder denn:
    -MAT(rix) LAB(oratory) als Programmname
    -2000€ die Basisversion, und ~6000€ je Toolbox
    -Nahezu 40 Jahre Entwicklungserfahrung
    -Feedback von Millionen professionellen Nutzern

    Das ganze natürlich unter der Voraussetzung dass man weiß wie man ordentlich MATLAB programmiert und welche Funktionen man auf welche Sorte von Matrizen anwendet. Denn auch MATLAB schützt nicht vor besch***enem Code 😉



  • Ich dachte immer Matlab wäre in java geschrieben.
    Oder gilt das nur für sachen wie die graphische oberfläche?
    Oder kriegen die das trotz java so schnell?



  • Q schrieb:

    Ich dachte immer Matlab wäre in java geschrieben.
    Oder gilt das nur für sachen wie die graphische oberfläche?
    Oder kriegen die das trotz java so schnell?

    MATLAB ist ein Mix aus vielen Sprachen je nachdem welche Elemente du verwendest. Viele Teile sind auch in der MATLAB Sprache geschrieben und sind damit ad hoc plattformunabhängig. Für Berechnungen wird soweit ich weiß ein optimiertes LAPACK:FORTRAN, FFTW:C verwendet.

    Ist auch richtig so, es macht keinen Sinn GUI Elemente in FORTRAN zu schreiben. Dadurch dass man sich nicht auf eine Sprache festlegt muss man nicht die externen Bibliotheken erst in die eigene Sprache übersetzen.



  • soweit ich weiss, verwendet MATLAB u.a. die ATLAS Bibliothek (automatically tuned linear algebra dingsbums). Aktuelle Versionen von Matlab können aber Operationen wie ein Matrixprodukt zwischen großen Matrizen noch auf die versch Kerne der CPU aufteilen.

    Was heutzutage nicht ganz unerheblich ist: Speicherzugriff. Wenn's schnell gehen soll, sollte man seine Datenstruktoren so konstruieren und Speicherzugriffe so sortieren (falls möglich), dass Cachemisses minimiert werden.



  • [quote="krümelkacker"]Aktuelle Versionen von Matlab können aber Operationen wie ein Matrixprodukt zwischen großen Matrizen noch auf die versch Kerne der CPU aufteilen.[/qoute]

    Das macht Eigen zum Glück auch.



  • Mit dem Faktor war die Zeit zur Bildung eines Produktes prod(A,B) gemeint.
    MATLAB:Eigen ist übrigens 1:3 ... also fast perfekt.

    kann ich nicht bestätigen. Zumindest kommt es ganz stark auf den Rechner an. Matlab verwendet wie bereits gesagt Atlas, das _echt_ langsamer ist. Das einzige, was Matlab macht, ist automatich multi threaded zu sein. Ist aber auch nicht unbedingt ein Vorteil, da bei vielen Numeric Anwendungen schon mit wenig Hirnschmalz ein besseres Mukltithreading über openmp zu haben ist.

    Ublas ist wirklich nicht ganz...optimal. Allerdings ist es exakt. Ich bin gerade dabei, für mein Projekt ein fast_prod mit den numeric bindings zu implementieren. Wenn ich aber ublas::axpy_prod mit cblas_gemm vergleiche (100x100 double Matrixmultiplaktion), haben im Resultat einzelne Elemente der Matrix eine Abweichung von 1.e-10. Die double multiplikation von atlas hat float-Genauigkeit. SO ist es nicht schwierig, einen Geschwindigkeitsschub zu kriegen. Für das Benchmarking von Optimierungsalgorithmen ist das allerdings total unbrauchbar. Dann lieber uBLAS.



  • otze schrieb:

    kann ich nicht bestätigen. Zumindest kommt es ganz stark auf den Rechner an. Matlab verwendet wie bereits gesagt Atlas, das _echt_ langsamer ist.

    Zeig uns doch ein Beispiel das angeblich langsamer sein soll dann zeige ich dir was du falsch machst 😉 Das "_echt_" ist immer ein Problem des Programmierers. Über Jahrzehnte machen die Leute Erfahrungen die du auf einmal über den Haufen wirfst. In Stackoverflow ärgern sich die Leute regelmäßig darüber dass MATLAB nur mit ungemein hohem Aufwand und damit verbundenen Sonderfällen und Kompromissen zu überholen ist. Und da programmierende Ingenieure zu teuer sind kauft man sich lieber gleich MATLAB.

    otze schrieb:

    Das einzige, was MATLAB macht, ist automatisch multi threaded zu sein.

    Automatisch geschieht im Computerwesen überhaupt nichts. Eine Weisheit die man als Jemand der schon einmal programmiert hat kennen sollte.

    otze schrieb:

    Ist aber auch nicht unbedingt ein Vorteil, da bei vielen Numeric Anwendungen schon mit wenig Hirnschmalz ein besseres Mukltithreading über openmp zu haben ist.

    1. OpenMP ist alles andere als langsam. Bestes Beispiel ist Intels MKL was OpenMP verwendet.

    2. Wie schon oben erwähnt... wenn du dich schon so rühmst dann bring erst einmal ein Beispiel sonst ist die Aussage einfach nicht haltbar.

    otze schrieb:

    Ublas ist wirklich nicht ganz...optimal. Allerdings ist es exakt.

    Was verstehst du unter exakt und wie kommst du darauf das MATLAB nicht exakter als ublas ist? Du kannst bei MATLAB alles steuern bis hin zu BCD oder symbolischen Berechnungen. Wenn du nummerisch arbeitest gibt es sowieso keine Exaktheit weil der IEEE753 keine Reelen Zahlen definiert.

    otze schrieb:

    Ich bin gerade dabei, für mein Projekt ein fast_prod mit den numeric bindings zu implementieren. Wenn ich aber ublas::axpy_prod mit cblas_gemm vergleiche (100x100 double Matrixmultiplaktion), haben im Resultat einzelne Elemente der Matrix eine Abweichung von 1.e-10.

    Ich kann dir jetzt schon vorhersagen dass dein Aufruf langsamer sein wird weil: "c_blas" ist 5% langsamer als pure FORTRAN Aufrufe. Das ist hinreichend bekannt und der einzige Grund warum FORTRAN noch eine Existenzberechtigung hat.

    In der Realität baut man auch sein Modell so auf dass man nicht auf nahezu double Genauigkeit spekulieren muss. Das ist einfach absurd zu denken dass bei 1e-10 noch sinnvolle NK zu suchen sind wenn man performant rechnen möchte.

    Außerdem ist es sowieso zu pauschal da du bei double nur eine Genauigkeit von 54 Bits in der Mantisse hast. 2^(-53) = 16 Nachkommastellen. Im günstigsten Fall wohlgemerkt.

    Die Abweichungen liegen auch wahrscheinlich an deiner Programmierweise. Such mal in Google nach dem Stichwort Thread Safe. Je nachdem wie du dein MATLAB oder die Bibliotheken verwendest sind diese nicht ad hoc thread safe und müssen anders konfiguriert werden. Das ist auch kein Geheimnis sondern gut in der Hilfe dokumentiert.

    otze schrieb:

    Die double multiplikation von atlas hat float-Genauigkeit.

    Ich übergehe diesen Satz mal ohne ein Kommentar abzugeben. MATLAB und NUR float Genauigkeit bei Produkten... das habe ich noch nie gehört. Meine Güte (!)

    otze schrieb:

    SO ist es nicht schwierig, einen Geschwindigkeitsschub zu kriegen. Für das Benchmarking von Optimierungsalgorithmen ist das allerdings total unbrauchbar. Dann lieber uBLAS.

    Ohje ohje... das einzige was richtig unbrauchbar ist ... nein ich lasse es.

    Wenn du es nochmal erwägst zu antworten dann informiere dich nochmal kurz über die Tatsachen und bringe ein Beispiel dann können wir diskutieren. Dein aktueller Stand ist eher Desiformation. Ein guter Rat noch: Ich würde dir raten dich nicht mit solchen Aussagen in der Öffentlichkeit zu äußern. Dann stehst du ruck zuck ohne Hose da ohne es zu merken. Noch schneller kann man sich praktisch nicht disqualifizieren. Sorry aber besser es sagt dir Jemand hier als in der Realität wo alles schon vorbei ist.

    Danke!


Log in to reply