numpy cholesky-zerlegung

[Arthur Dent]

Hallo

Ich muss Gleichungssysteme mit positiv definiter Koeffizientenmatrix, mit beliebig vielen rechten Seiten Lösen. Dazu bietet sich ja die Choleskyzerlegung an. Nun gibt es da zwar in numpy ne funktion, die mir die Zerlegung liefert, nur was mach ich jetzt damit? Stichwort: vorwärts und rückwärts Einsetzen. Muss ich mir für die Lösung selbst was schreiben, oder gibt es da auch schon irgendwas in numpy?

Gruß

Arthur Dent

[0x1cedd1ce]

Mit numpy hab ich nie gearbeitet. Aber Vorwärts- und Rückwärtseinsetzen ist ja kein großes Ding. Eine eigene Funktion bräuchte vermutlich zwei for-schleifen und dann hast du das Ergebniss.

[Arthur Dent]

Sicherlich ist das keine große Sache, Ich hab vor längerer Zeit mal im rahmen eines Numerik-Seminars selbst ein par Löser geschrieben. Nur ich hätte halt erwartet, dass bei dem durchaus beachtlichen Umfang von numpy und scipy schon ein Löser für zerlegte Matrizen dabei ist. Falls dem so wäre, fänd ich es aus auf jeden Fall eleganter, und für Dritte besser nachvollziehbar, die mitgelieferten Funktionen zu benutzen, als mir da selbst was hinzubasteln. Aus der Doku bin ich da leider nicht so richtig schlau geworden.

[0x1cedd1ce]

http://docs.scipy.org/doc/numpy/referen ... inalg.html
Da gibt es jede menge Funktionen zum lösen von Linearen Gleichungssystemen. 'tensorsolve' sollte das richtige für dich sein.

[Arthur Dent]

Da bin ich gestern Abend auch noch drauf gestoßen, Leider weiß ich nicht so genau was die Funktion eigentlich macht. ich will ja grad das Invertieren der Koeffizientenmatrizen vermeiden, weil dieser Aufwand schon in der Zerlegung steckt. So wie ich es Verstehe, invertiert die Tensorsolve-Funktion die zerlegten Matrizen nochmal oder.
Hab bei der Tensoralgebra aber echt nich so viel verstanden. Wenn dem also nicht so ist lasse ich mich gerne eines besseren belehren.

[0x1cedd1ce]

http://docs.scipy.org/doc/numpy/referen ... nalg.solve
Ich hab mir mal die Beschreibung von tensorsolve durchgelesen und du hast recht, das ist nicht was du brauchst.
solve hingegen ist das richtige. Du übergibst Matrix a und vektor b, das Ergebnis ist Vektor x; a*x = b

[Arthur Dent]

OK, du hast recht.
Ich hatte mal an anderer Stelle gelesen, dass die solve-Funktion für "dense matrices" - also vollbesetzte- eingesetzt werden sollte. Nichts desto trotz benutzt die Funktion die LU-Zerlegung aus der LAPACK-Bibliothek, damit lassen sich Gleichungssysteme mit regulärer, quadratischer Koeffizientenmatrix, für mehrerere rechte Seiten effizient lösen. Trotzdem wird hier die Symmetrie und die positive Definitheit leider nicht ausgenutzt. Außerdem bin ich mir relativ sicher, dass bei jedem Aufruf von "solve" zur Zerlegung das Gaußsche Eleminationsverfahren durchgeführt werden muss, selbst wenn die Koeffizientenmatrix schon auf Zeilen-Stufen-Form gebracht ist. Findet der Löser Nullen an der Stelle an der er grad ist bricht die Lösungsschleife zwar ab und macht beim nächsten Eintrag weiter, trotzdem muss immer über die halbe Matrix iteriert werden. Diesen Aufwand hab ich jedoch bereits mit der Cholesky-Zerlegung erledigt.
Ein weiteres Problem ist, das zur Laufzeit eine unbekannte Anzahl weiterer rechter Seiten dazu kommen könnte. Dann Müsste ich die Funktion jedes Mal (incl. des internen Zerlegugsaufwandes) wieder neu Aufrufen, obwohl sich an der Koeffizientenmatrix nichts geändert hat.
Ich danke dir wirklich für deine Mühe, allerdings scheint das auch nicht so das Richtige für meine Zwecke zu sein.
Ich hab mir jetzt mal "relativ quick and dirty" einen Lösungsalgorithmus zusammengebaut. wenn ich den n bissel schöner gemacht hab kann ich ja mal n Schnipsel posten.

vielen Dank

Gruß Arthur