[lineare Algebra] QR-Zerlegung

[theHacker]

Moin,

ich verzweifle langsam... ich möchte gerne eine QR-Zerlegung machen, aber es will nicht so wirklich hinhauen. Da meine Infos spärlich sind, wie jetzt konkret der Algorithmus funktioniert und auch das Netz nur mathematischen Formelwust liefert, möchte ich gerne von euch mal wissen, was ich falsch mache.

Konkret gehe ich so vor:
(im Folgenden sind die Zeilen und Spalten von Matrizen und Komponenten von Vektoren von 0 ab gezählt, nicht wie mathematisch normal mit 1 ^^)

• A[sub]0[/sub] := A

• Für i von 0 bis n-1
  • Ich hole mir aus A[sub]i[/sub] den Spaltenvektor x der (i+1)-ten Zeile und ersetze die ersten i Einträge mit 0
    
    Im Beispiel der 3x3-Matrix wäre der Vektor
    • für i = 0: x[A[sub]0, 0[/sub]; A[sub]0, 1[/sub]; A[sub]0, 2[/sub]]
    • für i = 1: x[0; A[sub]1, 1[/sub]; A[sub]1, 2[/sub]]
    • für i = 2: x[0; 0; A[sub]2, 2[/sub]]
  • von diesem Vektor x wird dessen Norm an der i-ten Komponente abgezogen, somit erhalte ich den Vektor ?[sub]i[/sub]
  • Mit ?[sub]i[/sub] kann ich mir Q[sub]i[/sub] gemäß folgender Formel bauen:
    Q[sub]i[/sub] = E - 2 (?[sub]i[/sub] * (?[sub]i[/sub])[sup]T[/sup]) / ??[sub]i[/sub]?[sup]2[/sup]​
  • Mit dem Q[sub]i[/sub] bau ich mir das A[sub]i+1[/sub] als
    A[sub]i+1[/sub] = Q[sub]i[/sub] * A[sub]i[/sub]​

• Zum Schluss ermittelte ich mein ganzes Q, Q[sub]gesamt[/sub] genannt, mit Q[sub]gesamt[/sub] = Q[sub]n-1[/sub] * ... * Q[sub]1[/sub] * Q[sub]0[/sub]

• R ist dann A[sub]n-1[/sub] oder mit Hilfe von Q ergibt sich R = Q[sup]-1[/sup] * A bzw. weil weil Q orthogonal ist als R = Q[sup]T[/sup] * A

Für 2x2-Matrizen funktioniert es. Bei den 3x3-Matrizen funktioniert es nicht mehr, wenn in der ersten Spalte die unteren beiden Zeilen belegt sind.
Drum tippe ich mal auf einen Fehler beim ?.

Auch ist mein letztes ? immer der Nullvektor, was für mich nicht wirklich viel Sinn ergibt, aber nach meinem Vorgehen eigentlich ganz logisch ist: Beim letzten Schritt besteht mein Vektor x nur noch aus einer Komponente. Ich fülle dann mit Nullen auf und bekomme mein ? als ?=[0; 0; ?]. Ziehe ich nun in der letzten Komponente einmal seine Norm, die ja wieder ? ist, ab, erhalte ich logischerweise den Nullvektor als ?.
Enthält A nun negative Einträge, dann gibts keinen Nullvektor. Da hab ich mich schon gefragt, ob vielleicht irgendwie ein Betrag irgendwo fehlt.

Ich weiß irgendwie gar nicht, ob ich es jetzt kapiert hab, wie es geht oder ob ich noch was falsch habe. Stimmt mein Vorgehen und das sind nur alles Rundungsfehler oder passt da was nicht?

Ihr könnt euch ja mal angucken, was ich so gemacht hab:
https://www.thehacker.ws/dummy/qr2.html
https://www.thehacker.ws/dummy/qr3.html

Einmal eine 2x2-Matrix, wo alles klappt und einmal eine 3x3-Matrix, wo sich am Ende alles unterscheidet. Entweder ich krieg A nicht wieder raus oder aber R ist keine Dreiecksmatrix.

 

[DaPhreak]

Ich kenne den Algorithmus zwar nicht, habe ihn aber mal so wie Du es hier beschrieben hast in Matlab nachprogrammiert und da funktioniert es.

Bissl unsauber, aber zum proof of concept soll's genügen:

A0 = [1,3,14;4,5,12;2,2,9];

x = A0(1:3,1);
x(1) = x(1) - norm(x);

Q1 = eye(3) - 2*x*x'/(x'*x);

A1 = Q1*A0;

x = A1(1:3,2);
x(1) = 0;
x(2) = x(2) - norm(x);

Q2 = eye(3) - 2*x*x'/(x'*x);

A2 = Q2*A1;


Qtotal = Q2*Q1;
Qtotal*A0


ans =

    4.5826    5.8919   17.4574
   -0.0000    1.8127    9.4573
    0.0000   -0.0000   -5.1766

Und da ist jetzt Qtotal*A0 eine obere Dreiecksmatrix, wie sich das gehört. (Zur Erklärung: X(i:j,k:l) ist eine Teilmatrix von X mit Zeilen i bis j und Spalten k bis l, eye(N) ist eine Einheitsmatrix der Größe N, das Hochkomma ' steht für Transponiert).

Algorithmus stimmt also.

Ich hab's mal mit Deinen Zahlen verglichen und komm offenbar auf's selbe.
Da fiel mir dann auf: Deine Rechnung stimmt doch für die 3x3-Matrix?

Bei den 3x3-Matrizen funktioniert es nicht mehr, wenn in der ersten Spalte die unteren beiden Zeilen belegt sind.

Ja, aber mit was: 1.0848140153852E-16 und 2.3616459605893E-17, also etwa 10[sup]-16[/sup]. Das sind genau die Rundungsfehler, die man bei doubles immer bekommt.

Warum? Double sind 2 Byte, also 64 Bit. Davon entfallen 52 Bit auf die Mantisse. Die minimale Distanz zwischen 2 Zahlen ist somit 2[sup]-52[/sup]=2.2204*10[sup]-16[/sup]. Auf dieses Intervall muss der Prozessor also runden (quantisieren) und damit bekommt man diese Größenordnung als Quantisierungsrauschen.

Auch ist mein letztes ν immer der Nullvektor, was für mich nicht wirklich viel Sinn ergibt, aber nach meinem Vorgehen eigentlich ganz logisch ist

Wieso, ist doch logisch. Im ersten Schritt werden in der ersten Spalte 2 Elemente zu null gemacht. Im zweiten Schritt in der zweiten Spalte 1 Element zu null gemacht. Damit ist nach dem zweiten Schritt das Ding schon dreieckig und Du bist fertig. Es reicht wenn die äußere Schleife bis n-2 läuft. Gibt ja nur in n-1 Spalten was zu tun.

 

[theHacker]

Ich hab's mal mit Deinen Zahlen verglichen und komm offenbar auf's selbe.
Da fiel mir dann auf: Deine Rechnung stimmt doch für die 3x3-Matrix?

Guck dir mal meine Rechnung nochmal genau an.
Im untersten Abschnitt hab ich 6 Ausgaben gemacht:

• Q
• R mit Hilfe an A[sub]n-1[/sub]: das ist eine Dreiecksmatrix
• dieses R von rechts an Q multipliziert gibt eben nicht A, wie es sein sollte
• hier ist nochmal R bestimmt, diesmal als Q[sup]T[/sup] * A: dieses R ist keine Dreiecksmatrix
• nehme ich dieses R aber und multipliziere von links das Q dazu, kommt A raus, wie es sein sollte.
• A zum Vergleich

Du siehst also, ich hab 2 verschiedene R:
Das eine ist eine Dreiecksmatrix, das andere gibt mit Q zusammen wieder A.
Ich such aber eigentlich ein R, was beide Eigenschaft hat.

Nehm ich nämlich das obere R wird mir mein Gleichssystem wohl nicht richtig gelöst werden, da es ja mit A nix mehr zu tun hat.
Nehm ich das untere R, kann ich nicht mehr auflösen.

Ja, aber mit was: 1.0848140153852E-16 und 2.3616459605893E-17, also etwa 10[sup]-16[/sup]. Das sind genau die Rundungsfehler, die man bei doubles immer bekommt.

Jo klar, selbst verständlich sind 10[sup]-1x[/sup] für die Rechnung mit 0 gleichzusetzen.

Ich meinte die Unterschiede zwischen den beiden R.

Wieso, ist doch logisch. Im ersten Schritt werden in der ersten Spalte 2 Elemente zu null gemacht. Im zweiten Schritt in der zweiten Spalte 1 Element zu null gemacht. Damit ist nach dem zweiten Schritt das Ding schon dreieckig und Du bist fertig. Es reicht wenn die äußere Schleife bis n-2 läuft. Gibt ja nur in n-1 Spalten was zu tun.

Ok, als letzten Durchgang kann ich mir sparen.

 

[DaPhreak]

Achso, da liegt das Problem. Na das ist dann aber nicht die Rechnung, sondern eher die Interpretation der Ergebnisse.

• Q
• R mit Hilfe an A[sub]n-1[/sub]: das ist eine Dreiecksmatrix
• dieses R von rechts an Q multipliziert gibt eben nicht A, wie es sein sollte
• hier ist nochmal R bestimmt, diesmal als Q[sup]T[/sup] * A: dieses R ist keine Dreiecksmatrix
• nehme ich dieses R aber und multipliziere von links das Q dazu, kommt A raus, wie es sein sollte.
• A zum Vergleich

Das R was du aus A[sub]n-1[/sub] bekommst ist die "triangularisierte" Version von A. Die Triangularisierung ist erreicht worden in dem Du sukzessive mit den Q[sub]n[/sub] multipliziert hast. Damit ist das Q[sub]gesamt[/sub] genau die Matrix, die A in eine Dreiecksmatrix R überführt: R = Q[sub]gesamt[/sub] * A.

Da Q wie Du sagst orthogonal ist, lässt sich das umschreiben zu
Q[sub]gesamt[/sub][sup]T[/sup] * R = A.

Wenn Du also also eine Zerlegung der Form Q*R=A suchst, muss das entgültige Q gebildet werden, in dem Du die ganzen Q's zusammenbaust und dann nochmal transponierst.

Qtotal = Q2*Q1;
R = Qtotal * A0;

Q = Qtotal';

Q =

    0.2182    0.9457   -0.2408
    0.8729   -0.0788    0.4815
    0.4364   -0.3152   -0.8427


R =

    4.5826    5.8919   17.4574
   -0.0000    1.8127    9.4573
    0.0000   -0.0000   -5.1766

Q*R =

    1.0000    3.0000   14.0000
    4.0000    5.0000   12.0000
    2.0000    2.0000    9.0000

 

[theHacker]

Juhuu, kapiert Jetzt hab ichs raus

Danke vielmals

 

[DaPhreak]

Danke vielmals

Keine Ursache, immer gerne.