Inverse einer Matrix - zu viele Rekursionen für C#?

[Kennung Eins]

Ich arbeite mit 13x13 Matrizen und muss diverse Transformationen daran durchführen, dazu notwendig ist unter Anderem eine Inversenbildung.

Bei der Inversenbildung arbeite ich mit der Adjunkten, die wiederum die Determinante der Matrix benötigt.

http://upload.wikimedia.org/math/e/f/6/ef6391aa82a3b2316e84d7fadb987e3a.png
(de.wikipedia.org)

Programmierseitig sind das ziemlich viele Rekursionen, hauptsächlich bei der Determinantenbildung, weil da eben von 13x13 auf 2x2 runtergebrochen werden muss.

Es entsteht ein Hierarchiebaum von Rekursionsaufrufen, der ganz unten am Boden super breit ist. Ich habe keine genaue Ahnung, wie viele 2x2 Matrizen gebildet werden, aber es müssen tausende sein (die Anzahl kann man sicher irgendwie berechnen - 13! = 6227020800 ?).

Jetzt kommt meine Berechnung ungefähr bis hoch zur Stufe 8x8 und hört dann dort ohne Fehlermeldung einfach auf.
http://img243.imageshack.us/img243/6648/13x13ca0.jpg (http://imageshack.us)

Gibt es hier jemanden, der eine Ahnung hat, was da das Problem sein könnte? Zu viele Rekursionen? Der Code funktioniert problemlos mit kleineren Matrizen (<8x8).


[edit]
Grad mal getestet: Die Speicherauslastung bleibt fast Konstant beim Ausführen der Berechnungen, daran sollte es also nicht liegen.

[edit2]
Ich seh grad, er stoppt immer BEVOR er die 8x8 erreicht. Wenn sich die nötige Anzahl von Rekursionen wirklich aus n! (n-Fakultät) berechnet, dann ist 7! = 5040 und 8! = 40320 .. und dazwischen liegt ja die berühmte 32768 - vielleicht gibts da wirklich eine Einschränkung?

[Trap]

Der Stack hat in den meisten Sprachen eine feste Größe, daher kann man das bei der Speicherauslastung natürlich nicht sehen.

Mit Stackproblemen speziell in .NET hab ich allerdings keine Erfahrung. Allgemein kann man das Problem fast immer dadurch lösen, dass man manuell irgendwas aufm Heap anlegt und damit dann die Rekursion mit einer Schleife und den Daten auf dem Heap simuliert.

[DocEW]

Hab jetzt mal kurz nachgedacht... warum nimmst du nicht den Gauß-Jordan-Algorithmus (http://de.wikipedia.org/wiki/Inverse_Matrix#Gau.C3.9F-Jordan-Algorithmus)?
Sollte schneller gehen.

Das, was du machst, ist der Laplace'sche Entwicklungssatz, oder?

Achja, und dann gibt es da noch so eine Wette von irgend 'nem Mathe-Prof, daß man niemals wirklich das Inverse einer Matrix benötigt... ;)

[Kennung Eins]

Ja, ist Laplace. Den kann ich mir "vorstellen", da weiß ich was passiert. Beim G.-J.-Algorithmus müsste ich erstmal durchsteigen *G* (vor allem Programmierseitig denke ich ist das wesentlich komplizierter)

Aber wenns wirklich so kompliziert ist, muss ich wohl oder übel was anderes nehmen - Heap&Stack-Arbeiten sind doch bei C# garnicht möglich, oder?

Ich danke euch :)

[Trap]

Aber wenns wirklich so kompliziert ist, muss ich wohl oder übel was anderes nehmen - Heap&Stack-Arbeiten sind doch bei C# garnicht möglich, oder?
Heap und Stack funktionieren in C# ganz genauso wie in allen anderen Programmiersprachen. Man hat zwar einen GC für den Heap, aber das ist unwichtig.

Was ich meine ist sowas (Pseudocode):
Stack st = new Stack();
while(!st.empty()) {
args[] = st.pop();
dostuff()
if(recurse) st.push(newargs);
}

[Kennung Eins]

Oh, das ist mir noch nie über den Weg gelaufen, danke :)
Bin grad an dem anderen Algorithmus dran.

[Coda]

Achja, und dann gibt es da noch so eine Wette von irgend 'nem Mathe-Prof, daß man niemals wirklich das Inverse einer Matrix benötigt... ;)

Was ein Blödsinn.

[Trap]

Was ein Blödsinn.
Wenn es so offensichtlich Blödsinn ist, dass du keine weitere Erklärung zu schreibst, kannst du sicher einige Beispiele von Anwendungen geben in denen man wirklich das Inverse einer Matrix haben möchte.

[Kennung Eins]

Ich brauch es für die Berechnung der Pseudoinversen.

Habs übrigens mit dem Gauß-Jordan-Algorithmus hinbekommen.

[Coda]

Wenn es so offensichtlich Blödsinn ist, dass du keine weitere Erklärung zu schreibst, kannst du sicher einige Beispiele von Anwendungen geben in denen man wirklich das Inverse einer Matrix haben möchte.

Basiswechsel in eine nicht orthogonale Basis und wieder zurück?

[tokugawa]

Wenn es so offensichtlich Blödsinn ist, dass du keine weitere Erklärung zu schreibst, kannst du sicher einige Beispiele von Anwendungen geben in denen man wirklich das Inverse einer Matrix haben möchte.

Light-Space-Perspective-Shadowmapping.

[Trap]

Basistransformationen kann man einfach in beide Richtungen konstruieren, die sind dann automatisch invers ohne dass man dabei eine Matrix invertieren muss.

[Xmas]

Allerdings ist das Invertieren ab einer gewissen Anzahl schneller.

[pajofego]

Obiger algorithmus ist nur für kleine Matrizen sinnvoll, in der Regel 3x3!

Gruß

pajofego

[Trap]

Allerdings ist das Invertieren ab einer gewissen Anzahl schneller.
Wie meinst du das? Wenn man eine Basistransformation konstruiert hat, hat man die Matrix für diese Transformation. Macht man einmal in jede Richtung, dann hat man 2 zueinander inverse Matrizen, die man auf beliebige Dinge anwenden kann.

Basistransformation konstruieren geht mit linearem Gleichungssystem lösen, das ist bei vielen Dimensionen auf jeden Fall schneller, ob Invertieren bei kleinen Matrizen irgendwann schneller ist weiß ich nicht.

[tokugawa]

Also, bei Light-Space-Perspective-Shadowmapping ist mit dem ganzen Matrizen-Rumgetue auf jeden Fall der Weg über die Inverse bequemer/intuitiver, und eventuell sogar effizienter. Aber der Kern von Light-Space-Perspective-Shadowmapping ist ja auch keine wirkliche Basistransformation.

[DocEW]

Was ein Blödsinn.
Klar dürfte doch sein, was hinter dieser Wette für eine Aussage steckt: Handel dir nicht den Ärger einer Matrixinvertierung ein, es gibt immer numerisch "gutartigere" Verfahren. Außerdem benutzt man die Matrixinvertierung einfach oft unnötiger Weise, weil sie - wie tokugawa ja auch schreibt - "bequemer/intuitiver" ist. Man könnte z.B. auf die Idee kommen, ein Lineares Gleichungssystem Ax=b damit zu lösen. Invertiert man dazu A, hat man implizit die Lösung für alle rechten Seiten berechnet, obwohl man es ja nur für ein spezielles b lösen möchte.

[Coda]

Klar dürfte doch sein, was hinter dieser Wette für eine Aussage steckt

Ja, da stimm ich ja zu, aber das es gar keine Anwendung für eine Matrixinverse gibt ist nun doch wirklich nicht wahr.

[Gast]

also ne inverse dieser größe mit Gauß Jordan zu machen ist meiner meinung schon ein dicker happen (programmiertechnisch stelle ich es mir schwieriger vor, aber hab keien ahnung vom proggen)

aber die adjunkten form kannst du vielleicht vereinfachen indem du dir überlegst wie du einfacher ne determinante kriegst:
wandle die 13x13 zu ner Dreiecksmatrix. das sollte sich programm technisch noch machen lassen. die determinante ist dann recht einfach.

[Neomi]

Die Breite des Rekursionsbaumes spielt (außer für die Laufzeit) keine Rolle, da die nebeneinander liegenden Unterbäume auch nacheinander abgearbeitet werden. Für den Stackbedarf ist alleine die Tiefe des Baums relevant, und die ist unbedenklich niedrig. Der Fehler muß daher an anderer Stelle liegen.

Und nochmal dazu:

Achja, und dann gibt es da noch so eine Wette von irgend 'nem Mathe-Prof, daß man niemals wirklich das Inverse einer Matrix benötigt... ;)

Sollte es diese Wette wirklich geben, dann bestimmt ganz anders formuliert. Daß man für verschiedene Sachen die Inverse einer Matrix benötigt, sollte klar sein. Ob man die Inverse bekommt, indem man die Matrix invertiert, oder ob man sie auch einfach neu konstruieren kann, ist eine andere Sache. Ob man also eine Matrixinvertierung (die Funktion, keine konkrete Anwendung) benötigt oder ob es immer einen anderen Weg gibt, darüber würde ich keine Wetten abschließen wollen. Und streng genommen ist eine Invertierung schon dann unnötig, wenn die Möglichkeit der Neukonstruktion zwar langsamer, aber vorhanden ist.

[DocEW]

Sollte es diese Wette wirklich geben, dann bestimmt ganz anders formuliert.
Ja klar, war ja aus'm Kopf. ;)
Also dann - hier mal das Zitat aus "Numerische Mathematik I" von Deuflhard/Hohmann:
Es gab über lange Zeit die offene Wette eines wissenschaftlich hochkarätigen Kollegen, der eine hohe Summe darauf setzte, dass in praktischen Fragestellungen niemals das Problem der "Invertierung einer Matrix" unvermeidbar auftritt. Soweit bekannt, hat er seine Wette in allen Einzelfällen gewonnen.

[tokugawa]

Klar dürfte doch sein, was hinter dieser Wette für eine Aussage steckt: Handel dir nicht den Ärger einer Matrixinvertierung ein, es gibt immer numerisch "gutartigere" Verfahren. Außerdem benutzt man die Matrixinvertierung einfach oft unnötiger Weise, weil sie - wie tokugawa ja auch schreibt - "bequemer/intuitiver" ist. Man könnte z.B. auf die Idee kommen, ein Lineares Gleichungssystem Ax=b damit zu lösen. Invertiert man dazu A, hat man implizit die Lösung für alle rechten Seiten berechnet, obwohl man es ja nur für ein spezielles b lösen möchte.

Wie ich aber auch geschrieben habe glaube ich dass bei der Verzerrung des Lightspace bei LiSpSM man um eine Invertierung vernünftigerweise gar nicht herumkommt.

Einfach die inversen Matrix-Operationen nebenläufig mitzuerstellen dürfte im Endeffekt rechenaufwändiger sein.

[DaFakka]

http://de.wikipedia.org/wiki/Liste_numerischer_Verfahren

Schau mal beim Gauß-Verfahren nach dem Aufwand und überlege, was das
in deinem Fall bedeutet.
Falls eine nicht zwingend exakte Lösung erforderlich ist, so probier mal das
SOR-Verfahren.

[Kennung Eins]

Ich brauche es leider 100% korrekt, wie gesagt, für die MP-Pseudoinverse.
http://de.wikipedia.org/wiki/Pseudoinverse#Moore-Penrose_Pseudoinverse

Wie oben schon geschrieben, mit Gauß-Jordan hats problemlos geklappt :)
(auch wenn bei manchen Matrixzellen sowas wie 2,00000000001 rauskommt - Danke DOUBLE! :( )

Vielleicht probier ich die Dreiecksmatrix-Methode für die Determinantenbestimmung nochmal aus - und benutz das wieder mit der Adjunkten-Methode. Mal sehen.

[Trap]

Mit der Singulärwertzerlegung gibt es jedoch einen einfacheren Weg, die Pseudoinverse zu erhalten.
Singulärwertzerlegung gibt es fertig in LAPACK. LAPACK benötigt aber etwas Einarbeitung, ist eigentlich eine Fortran-Bibliothek und damit etwas unkonventionell aus C-Sicht. Die GNU Scientific Library hat auch ein Singulärwertzerlegen und ist in C...

Die Probleme mit den Fließkommafehlern wird man damit leider nicht los. Man kann auch das ganz mit Brüchen rechnen, ist aber deutlich langsamer. Solange der Fehler so winzig ist wie in deinem Beispiel macht es aber nichts aus, außer dass man bei der Ausgabe etwas aufpassen muss.

[Shink]

Ich brauche es leider 100% korrekt, wie gesagt, für die MP-Pseudoinverse.
http://de.wikipedia.org/wiki/Pseudoinverse#Moore-Penrose_Pseudoinverse

Wie oben schon geschrieben, mit Gauß-Jordan hats problemlos geklappt :)
(auch wenn bei manchen Matrixzellen sowas wie 2,00000000001 rauskommt - Danke DOUBLE! :( )

Vielleicht probier ich die Dreiecksmatrix-Methode für die Determinantenbestimmung nochmal aus - und benutz das wieder mit der Adjunkten-Methode. Mal sehen.
Ich würds auch mit der SVD machen (hatte damit schon mal gröber zu tun).

Ein exaktes Ergebnis wird wohl generell nicht ganz einfach werden auf einer herkömmlichen CPU - aber du kannst z.B. einfach auf eine bestimmte Genauigkeit runden (z.B. 10^-5) oder bei iterativen Verfahren so lange rechnen, bis du die erwünschte Genauigkeit erreicht hast.

[aths]

Zum Lösen linearer Gleichungssysteme nutze ich gerne Cramer, wenn ich nicht alle Lösungen brauche. Das war mal bei einem bestimmten Problem der Schattenberechnung (Projektion eines Dreiecks auf eine Ebene) der Fall.

Ja, ist Laplace. Den kann ich mir "vorstellen", da weiß ich was passiert. Beim G.-J.-Algorithmus müsste ich erstmal durchsteigen *G* (vor allem Programmierseitig denke ich ist das wesentlich komplizierter)

Aber wenns wirklich so kompliziert ist, muss ich wohl oder übel was anderes nehmen - Heap&Stack-Arbeiten sind doch bei C# garnicht möglich, oder?

Ich danke euch :)Entwicklung auf 2x2 nach Laplace habe ich auch mal geproggt, was glaube ich bis 10x10 auf jeden Fall lief. Da werden doch aber gar nicht sooo viele Untermatrizen gebildet.

- Es reicht auch 3x3, da man schon dann die Determinante direkt ausrechnen kann.

- 13 Fakultät ist sicher nicht die Zahl der Untermatrizen.

- Laplace ist nummerisch ungünstig, wenn man mit Float rechnet. Um exakte Ergebnisse zu bekommen, müsste man irgendwie mit Brüchen rechnen.

[Silpion]

Wie schon vorgeschlagen, ist SVD die Methode der Wahl, wenn es numerisch einigermaßen genau sein soll. Das Gauss-Verfahren wird schnell instabil, insbesondere bei fast singulären Matrizen.