4. Wie indifinites Least-Squares Problem lösen?

Wie indifinites Least-Squares Problem lösen?

  • O

    Hi Leute!

    Ich habe ein Least-Squares Problem

    minx∣∣Ax−b∣∣2\min_x ||Ax-b||^2
    wobei A eine dichte symmetrische indefinite nxn Matrixmit rank(A) = r < n ist.
    Ich weiß, dass ich das Problem mit der Moore-Penrose inverse lösen kann, aber das ist mit Kanonen auf Spatzen geschossen. Eine LDL^T Zerlegung habe ich aber nicht implementiert und die Paper die ich darüber lese sind alles andere als vertrauenserweckend.

    Nun ist A nicht völlig beliebig, sondern hat die folgende Struktur:

    A=(Q11T0)A=\begin{pmatrix}Q & 1 \\ 1^T & 0 \end{pmatrix}
    Q ist eine n-1 x n-1 dichte symmetrische positiv-semi-definite matrix. 1 ist jeweils der 1-vektor und 0 ist eine einfache skalare 0. Es sind also nur ein Zeile und Spalte die mich davon abhalten ein numerisch einfacheres Verfahren anzuwenden.

    Irgendwlche Ideen, wie ich das Gleichungssystem zerlegen kann?

    0
  • K

    Ich glaube, da gibt es Verfahren, die Dir auch ohne SVD eine Minimum-Norm-Lösung eindeutig berechnen können. Such doch mal nach "minimun norm least squares" oder sowas.

    Eventuell ist diese Regularisierung hier auch eine Alternative:

    minx∣∣Ax−b∣∣2+∣∣ϵx∣∣2\min_x ||Ax-b||^2 + ||\epsilon \; x||^2

    Das ist äquivalent zu

    minx∣∣A‾x−b‾∣∣2\min_x ||\overline{A}x-\overline{b}||^2

    mit

    \overline{A} = \begin{pmatrix} A \\ \epsilon \; I \end{pmatrix} \\ \overline{b} = \begin{pmatrix} b \\ 0 \end{pmatrix}

    I steht dabei für die Einheitsmatrix und 0 für einen passend langen 0-Vektor. Das System hat vollen Rang bei ϵ>0\epsilon>0 und wenn der Regularisierungsparameter klein genug ist, ist das auch ziemlich dicht an der Lösung, die dich interessiert. Ggf kann man das ein paar mal iterieren

    x_0 = 0 \\ x_{i+1} = x\_i + \arg\min\_t ||\overline{A} \, t-\overline{\left(b-A x_i\right)}||^2

    damit man beliebig dicht an die Lösung des ursprünglichen Problems rankommt, ohne das Epsilon zu klein wählen zu müssen. Eine QR-Zerlegung von A‾\overline{A} muss man ja nur einmal berechnen.

    Vielleicht hilft aber auch schon eine QR-Zerlegung von A mit Spaltenpivotisierung. Wenn dann in R die letzte Zeile komplett 0 ist, darfst du sie ignorieren.

    Edit: Ja, ich glaube, QR mit Pivotisierung löst dir schon dein Problem -- zumindest dann, wenn du den Rang von A schon kennst.

    0
  • O

    Danke schonmal für diee Antwort!

    krümelkacker schrieb:

    Ich glaube, da gibt es Verfahren, die Dir auch ohne SVD eine Minimum-Norm-Lösung eindeutig berechnen können. Such doch mal nach "minimun norm least squares" oder sowas.

    Hab ich schon. Die üblichen Lösungsstrategien sind die diversen Berechnungen der Pseudoinversen. Das heißt also entweder QR-Faktorisierung, SVD oder rank-revealing cholesky mit der Identität:

    A^TA = LL^T \enspace (\text{rank revealing cholesky})\\ x = A^{\dagger}b = (A^TA)^{\dagger}A^Tb = L^{\dagger} (L^T)^{\dagger}A^Tb = L(L^TL)^{-2}L^T A^Tb

    Die drei funktionieren, wobei das dritte deutlich(!) schneller als die beiden anderen Verfahren ist(mir ist numerische exaktheit nicht so wichtig, es muss nur fertig werden...). Die drei Lösungen sind aber der große Hammer auf das Problem, weil sie nicht die Symmetrie von A ausnutzen. Das sieht man an der dritten Methode am besten: Es wird zuerst aus der symmetrisch indefiniten matrix eine positiv semi-definite matrix gebastelt von der dann die rangzerlegung gebildet wird. Dann wird die Pseudoinverse auf die Teile der Rangzerlegung angewandt mit deren Hilfe schließlich das Gleichungssystem gelöst wird.

    Das wird mit der QR-Zerlegung auch tendentiell nicht anders. Hier erspare mir den Schritt, A^TA zu berechnen, löse aber am Ende doch das selbe Gleichungssystem wie bei 3, nur dass ich direkt A^TA = RR^T durch A =QR erhalte. Allerdings ist die QR Zerlegung mit pivoting auch deutlich langsamer als rank revealing cholesky. Einen tod stirbt man immer :(.

    Mit einer LDL^T Faktorisierung hingegen könnte ich direkt die Rangzerlegung von A berechnen:

    A=LDL^T\\ x = A^{\dagger}b = (LDL^T)^{\dagger}b = L^{\dagger} D^{-1}(L^T)^{\dagger}b = L(L^TL)^{-1} D^{-1} (L^TL)^{-1}L^Tb

    Dadurch erspare ich mir die kosten der Erzeugung von A^TA, komme etwas günstiger weg als mit QR und habe nicht das Problem dass ich am Ende ein Gleichungssystem mit wesentlich größerer Konditionierungszahl K(A) löse.

    Zu deinen Verfahren:

    Der erste Trick mit der Regularisierung ist ja eigentlich damit äquivalent, auf die Diagonale von A ε aufzuaddieren. also A+εI. Dieses Iterationsverfahren habe ich schon gesehen, aber ich glaube, dass das Verfahren nur für positiv semi-definite Matrizen geeignet ist. Denn wenn epsilon nahe an dem absolutwert eines negativen Eigenwertes ist, kann K(A+εI) beliebig explodieren. Dann konvergiert das Verfahren eventuell niemals - und im schlimmsten Fall ist A+εI wieder singulär.

    0
  • K

    otze schrieb:

    Der erste Trick mit der Regularisierung ist ja eigentlich damit äquivalent, auf die Diagonale von A ε aufzuaddieren. also A+εI.

    Nee, das stimmt nicht. So, wie ich das aufgeschrieben habe, ist das immer regulär bei eps>0, egal wie A aussieht. Wenn du von dem "ergänzten" System das Normalgleichungssystem aufstellst, dann ja ... dann sitzt das Epsilon mit auf der Diagonalen:
    (ATA+ϵ2I)x=ATb\left(A^T A + \epsilon^2 I\right) x = A^T b
    aber das möchtest du vielleicht nicht so aufstellen und lösen, weil es numerisch schlechter ist als zB. das Lösen über eine QR-Zerlegung von A‾\overline{A}.

    0
Anmelden zum Antworten