MastersForum

Quoted

Original von DRDK_Fragman
Ich studiere Informatik. Naja, es handelt ich um Matrizen die bei der Lösung der Euler-Lagrange Gleichungen von Funktionalen auftreten. Das Anwendungsgebiet ist das Berechnen sogenannter Optischer Flußfelder in Bildsequenzen.
Es handelt sich um Gleichungen mit vielen Unbekannten (eine pro Pixel in den Bilder der Sequenzen), also so in der Größenordnung von 16384 (128x128 Bilder) oder 65536 (256x256 Bilder). Das bedeutet Pentadiagonalmatrizen haben deutlich weniger als 10% nicht null Einträge.

Ich benutze derzeit für manche Sachen die newmat Bibliothek, aber die kann zwar Bandmatrizen, aber nur solche wo das Band in der Mitte verläuft. Meine haben aber die Struktur das in der Mitte 3 Diagonalen sind und und dann außen noch mal (auf der Hälfte der Seiten) (bin nicht sicher ob man die dann überhaupt Pentadiagonal nennt). Die Methoden die verwende sind mit Differentialgleichungen recht eng verwand.

Ich schreibe die Arbeit beim Lehrstuhl für Angewandte Mathematik. Die verwenden selbst Multigrid Verfahren die sie selbst in C implemtiert haben, aber da die nicht sonderlich modular sind ist es mir zu Blöd die von Hand so anzupassen das man sie einsetzen kann. Möglicherweise ist SOR in meinem Fall auch ausreichend, das konvergiert für den geringen Aufwand ja auch schon ganz gut. Aber wenn das nicht reicht werde ich mich nochmal nach geeigneten Bibliotheken umsehen.

Hm ok, du hast also eine Tridiagonalmatrix und außen nochmal 2 Bänder... das ist was anderes, was man nicht mehr Pentadiagonalmatrizen nennt.

Wenn die Matrix nicht nur sym. pos. semi-definit sondern sym. pos. definit wäre, dann wäre wohl eine normale Cholesky-Zerlegung am besten, da hast du einen Aufwand von O(n) jeweils für Speicher und Flops.

Man kann Cholesky theoretisch auch auf semi-definite Matrizen anwenden, aber üblich ist es nicht.

Du solltest sowieso mal sagen, was der Löser für dein LGS machen soll, wenn A nur noch semi-definit und nicht mehr definit ist. Im eigentlichen Sinne existiert ja dann keine Lösung x=A^-1 b mehr, also was genau berechnest du dann? Eine Pseudoinverse?

Vielleicht ist deine Matrix auch gar nicht semi-definit, du hast bisher nur so schlechte Solver verwendet, dass numerische Instabilitäten sie semi-definit erscheinen haben lassen.

Das Verfahren der Wahl wäre wohl nach Cholesky ein CG-Verfahren als no-Brainer. Der nicht-triviale Kern bei nur semi-definiten Matrizen (es werden Vektoren außer dem Nullvektor von der Matrix auf Null abgebildet, daher der Rangabfall) stört nicht, wenn der Startvektor orthogonal zum Kern und die rechte Seite immer im Bild liegen.

CG solltest du also am ehesten ins Auge fassen.

Vorher solltest aber erst einmal mehr über deine Matrix erfahren. hast du mal bei einer kleineren die EW ausgerechnet und sie geplottet? ich schätze mal die fallen irgendwie exponentiell ab und du bekommst durch numerische Ungenauigkeiten welche Null hin.
Je nach der Größe der Matrix und der Kardinalität der Menge, aus welcher die Matrixeinträge sein können kann sich natürlich die Wahrscheinlichkeit erhöhen, überhaupt einmal auf einen pathologischen Fall zu treffen, d.h. verschwindende EW zu haben und damit semi-definit zu sein.

Vielleicht postest du mal eine allgemeine 16x16 Matrix hier, oder besser in einem neuen Thread. Dann solltest du sicher von ein paar Leuten Tipps bekommen, ob die Matrix überhaupt a) wirklich nur semi-definit ist bzw. b) ob sich durch geeignete Verfahren verhindern lässt, dass eine definite Matrix semi-definit erscheint und c) wie man letztlich die Inverse am besten berechnet bzw. das LGS löst.

Wenn du nur auf eine Lirary verlinken willst, dann würde ich heute abend mal ein CG-Verfahren probieren.

Quoted

Original von DRDK_Fragman
Ich benutze derzeit für manche Sachen die newmat Bibliothek, aber die kann zwar Bandmatrizen, aber nur solche wo das Band in der Mitte verläuft.

Wenn du noch nicht zu viel gemacht hast, dann würde ich definitiv GNU Scientific Library oder Boost (uBLAS) als Bibliotheken benutzen.

@ZwergOrca & Mari:
Man sollte dann ja auch mit Bibliotheken arbeiten. Benutzt man Boost, kann man ein CG-Verfahren praktisch in 5-10 Zeilen Code runterschreiben und ist dabei noch effizient.
C++ heißt ja nicht, dass man das Rad jedes mal neu erfinden muss.

Quoted

Original von Coold0wn

Quoted

Original von -=)GWC(RaMsEs
und im übrgien lernt man sehr leicht Java wenn man c++ kann, aber andersrum wird das schon schwieriger ...

ja weil man das schwierige hinter sich hat
c++ is doch viel härter als java oder?

Am Anfang wirst du keinen großen Unterschied merken.

Quoted

Original von Coold0wn
also compilieren - kein Problem (mehr , thx leute =) )
aber wie lass ich das dann laufen?
oO

Wenn du eine gute IDE hast, sollte das Programm nach dem compilieren automatisch starten.
Ansonsten wird unter Windows durch das linken und compilieren eine exe erstellt - auf die doppelt draufklicken.

Eventuell brauchst du das System("PAUSE"); vom Anfang des Threads vor return 0; damit das Fenster offen bleibt und du überhaupt etwas siehst.

Quoted

Original von -=)GWC(RaMsEs

Quoted

Original von Coold0wn

Quoted

Original von -=)GWC(RaMsEs
und im übrgien lernt man sehr leicht Java wenn man c++ kann, aber andersrum wird das schon schwieriger ...

ja weil man das schwierige hinter sich hat
c++ is doch viel härter als java oder?

kommt drauf an.
java nimmt dir halt sehr viele dinge wie zum beispiel speichermanagement ab. c++ hat eben kein garabge collector. Zusätzlich ist, wenn man mit pointern rbeitet es alles etwas schieriger zu handeln. das ist in java deutlich leichter. deswegen ist wenn man c++ mal kann java nicht mehr wild, aber anders rum schon.

Man muss sich in C++ nicht notwendigerweise drum kümmern, kann aber. Das ist die stärke: wenig muss, viel kann.

Ansonsten wird ja immer wieder gesagt: "Know your templates" und STL (Standard Template Library). Wenn man ganz ohne programmiert, sieht so ein C++ Code natürlich archaisch aus, aber das macht ja auch kaum jemand.

@Worf:

Also danke erstmal das du dich da so reinbegibst. Ich habe Donnerstag einen Vortrag über ein anderes Thema, das heißt bis dahin werde ich nichts mehr machen. Ab Donnerstag muss ich erstmal zusehen das ich Qt unter OS X mit Eclipse ans Laufen kriege (da muss ich nämlich irgendwie deren Precompiler in die Compiler Toolchain reinbekommen). Das brauche ich um erstmal meinen Algorithmus zu debuggen (weil sich das bei der Problemgröße so schwer ohne Visualisierung machen lässt). Wenn dann mein Code läuft, kann ich mich mit Performance befassen. Dann werde ich meine Matrizen noch mal genau checken. Es könnte gut sein das sie sogar positiv definit sind. Aber wie gesagt: Möglicherweise langt für meine Zwecke auch SOR. Wenn ich mehr Zeit habe werde ich das natürlich optimieren, aber ich habe nur bis Ende Juli und letztlich geht es ja nicht um Gleichungslösungsverfahren. Das ist hier nur Mittel zum Zweck. Außerdem kann es passieren das mein Algorithmus zur Berechnung des Flußfeldes nicht genau genug arbeitet, in dem Fall könnte es sein das ich anstelle eines linearen Gleichungssystems ein nonlineares lösen müsste, dann würde ich mit Lagged-Nonlinearity-Method arbeiten oder mit Warping, sodass meine Gleichungssysteme dann möglicherweise wieder ganz anders aussehen als jetzt. Die Fazit ist also: Es ist noch viel zu früh um großartig zu optimieren.

Zu den genannten Libraries: Ich hatte ursprünglich MTL 4 benutzt, die arbeitet ja unter der Haube angeblich auch mit BLAS, war aber um ein vielfaches langsamer und dazu hässlicher (weil mächtiger). Wenn du Boost oder GNU SL enpfiehlst schaue ich mir die noch mal an, aber erst wenn der Rest läuft. Sonst stehe ich da nachher mit ner numerisch ultra optimierten Methode da die aber einfach nicht funktioniert weil ich nicht genug Zeit genommen habe die Algorithmen geeignet zu kombinieren...

Du hast recht, der bessere Algorithmus schlägt eine bessere Implementierung um längen.

Genauigkeit ist allerdings sowohl eine Funktion des Problems (Kondition) sowie des verwendeten Algorithmus. Mit dem falschen Algorithmus (z.B. Gauss) wird die Genauigkeit der Lösung auch sehr schlecht sein.

Die Lagged-Nonlinearity-Method sagt mir nichts, Warping auch nicht. Entweder es sind sehr spezielle Verfahren für eine mir unbekannte Anwendung, oder ich kenne sie einfach nicht. -_-

Das wichtigste ist sicher erst einmal, Ergebnisse hinzubekommen. Allerdings kann ich nicht zustimmen, dass du dich erst ganz am Ende um Libraries kümmern solltest. Die Librarys, Boost beispielsweise, stellen alle ein Interface für Daten zur Verfügung, z.B. Matrizen sind jetzt nicht neu, es ist einfach nur eine gute Implementierung. Aber andere Dinge in Verbindung mit Boost sind neu, und so ist alles gut aufeinander abgestimmt.
Hat man eine Sparse- oder Bandmatrix, dann sind diese Typen sehr gut implementiert.

Das ist schon sehr nett und damit lässt sich auch schnell arbeiten, weil man schnell Ergebnisse erzielen kann.

Prototypst du auch in C++ oder in Matlab oder Mathematica oder etwas ähnlichem?

Wenn du nochmal über deine Problemstellung reden willst und vielleicht Anregungen suchst, dann mach doch einfach im Masters ein neues Thema auf.

schau hier: Ausgelagert: Optic Flow Computation