Givens Rotation & Householder < Maple < Mathe-Software < Mathe < Vorhilfe
|
| Aufgabe | Ich habe diese Frage in keinem Forum auf anderen Internetseiten gestellt.
Habe folgende Aufgabenstellung: QR- Zerlegung einer n x n Matrix mit Givens- Rotation und Householder- Transformation in Maple.
|
Hallo!
Meine Maple- Kenntnisse sind leider nicht die Besten.
Ein detailiertes Bsp. wäre schon eine gute Stütze für mich.
Kann mir bitte jemand bei diesem Thema helfen!
DANKE!
Dominik
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 01:20 Fr 02.11.2007 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
| |
|
| Status: |
(Frage) beantwortet  | | Datum: | 22:00 So 04.11.2007 | | Autor: | DominikW |
Hallo!
Habe noch immer keine Idee wie ich das implementieren sollte. Falls mir jemand helfen kann wär ich sehr dankbar.
Also falls jemand was weiß, bitte posten.
Danke!
|
|
|
| |
|
Hallo,
ich verstehe das so, dass du die Algorithmen schon verstehst und nur an Maple scheiterst.
Du solltest dir den Algorithmus mal skizzieren (in Pseudocode) und dir dann mal das Package LinearAlgebra oder linalg ansehen. An den Stellen, an denen es hakt, können wir dir dann ggf. weiterhelfen.
Zum Einarbeiten kannst du dir ja mal die Householder-Variante in Maple bei ![[]](/images/popup.gif) Wikipedia ansehen.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 13:40 Mo 05.11.2007 | | Autor: | DominikW |
Hallo!
Danke für deine Hilfe. Ich werde mal probieren das Ganze in Pseudocode zu verfassen. Poste dann sicher nocheinm zu diesem Thema.
lg
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 20:08 Mo 05.11.2007 | | Autor: | DominikW |
hallo!
Habe ein neuerliches Problem.
Wenn ich die Position [mm] a_{ij} [/mm] Null setzen will dann berechne ich für die Givens- Rotations Matrix p = [mm] sign(a_{jj})*\wurzel{a_{jj}^{2} + a_{ij}^{2}}
[/mm]
-> c = [mm] a_{jj} [/mm] / p und s = [mm] -a_{ij} [/mm] / p
Aber woher weiß ich, an welche Stelle der Matrix ich c, s und -s setzen muss. Liege ich da richtig, dass die Rotationsmatrix eine Einheitsmatrix ist mit der gleichen Größe wie die ursprüngliche Matrix und das nur an 4 Stellen c, s und -s gesetzt wird. (laut Wikipedia)
Und dann multipliziere ich die ursprüngliche Matrix mit der Rotationsmatrix und das so lange, bis ich eine Dreiecksmatrix habe.
Danke für eure Hilfe im Voraus!
LG Dominik
|
|
|
| |
|
Hallo,
>Aber woher weiß ich, an welche Stelle der Matrix ich c, s und -s setzen muss.
Nun, da das Verfahren numerisch sehr stabil ist, ist keine Pivotisierung nötig und du kannst stumpf alle Stellen abklappern, die am Ende 0 werden sollen. Bei einer MxN-Matrix A mit [mm] M$\ge$N [/mm] kannst du so vorgehen:
for spalte from 1 to N do
for zeile from M to spalte+1 by -1 do
Berechne transponierte Givens-Matrix GT(spalte, zeile);
Setze A := GT * A;
od;
od;
Du kannst aber auch schauen, welche Einträge von A noch nicht 0 sind und diese gezielt löschen.
> Liege ich da richtig, dass die Rotationsmatrix eine Einheitsmatrix ist mit der gleichen Größe wie die ursprüngliche Matrix und das nur an 4 Stellen c, s und -s gesetzt wird.
Ja!
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 23:27 Di 06.11.2007 | | Autor: | DominikW |
Hallo!
habe jetzt alles in Maple probiert zu programmierne. Erhalte aber immer eine Fehlermeldung und finde den Fehler nicht.
Stimmt das, wie die Werte c, s, -s in der Rotationsmatirx positioniert sind?
Kann mir bitte jemand helfen?
Hier mein Code:
rotation := proc (A::Matrix) local m, n, i, j, p, c, s;
global R; m := ColumnDimension(A);
n := RowDimension(A);
if (m = n) then
for i from 2 to n do
for j to n-1 do
p := sign(A[j, j])*sqrt(A[j, [mm] j]^2+A[i, j]^2);
[/mm]
c := A[j, j]/p;
s := -A[i, j]/p;
R := '<,>'('<|>'(1, 0, 0), '<|>'(0, 1, 0), '<|>'(0, 0, 1));
R[i, j] := s;
R[j, i] := -s;
R[i, i] := c;
R[j, j] := c;
A := multiply(A, R)
end do
end do
end if
end proc
---------------------
A ist eine nxn Matrix muss es nur für diesen Fall implementieren.
Danke im Voraus!
lg Dominik
|
|
|
| |
|
Hallo,
ich erlaube mir einfach kommentarlos eine korrigierte Fassung zu posten. Angesichts der späten Stunde lasse ich dich deine Fehler suchen ;)
| 1: | rotation := proc (origA::Matrix) local A, m, n, i, j, p, c, s;
| | 2: | global R;
| | 3: | A := origA;
| | 4: | m := LinearAlgebra[ColumnDimension](A);
| | 5: | n := LinearAlgebra[RowDimension](A); printf("%d",n);
| | 6: | if (m = n) then
| | 7: | for j from 1 to n do
| | 8: | for i from j+1 to n do
| | 9: | p := sign([j, j])*sqrt(A[j,j]^2+A[i,j]^2);
| | 10: | c := A[j, j]/p;
| | 11: | s := -A[i, j]/p;
| | 12: | R := <<1, 0, 0> | <0, 1, 0> | <0, 0, 1>>;
| | 13: | R[i, j] := -s;
| | 14: | R[j, i] := s;
| | 15: | R[i, i] := c;
| | 16: | R[j, j] := c;
| | 17: | A := R^%T.A;
| | 18: | end do ;
| | 19: | end do ;
| | 20: | return A;
| | 21: | else
| | 22: | error "Keine quadratische Matrix!";
| | 23: | end if ;
| | 24: | end proc; |
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 16:20 Mi 07.11.2007 | | Autor: | DominikW |
Hallo!
Danke für dein Programm. Habe jetzt auch beim mir den Fehler gefunden. Habe aber noch immer einen Fehler mit dem ich gar nichts anfangen kann. Folgender Fehler:
rotation(A)
Error, (in rotation) invalid left hand side in assignment
Kannst mir wer dabei helfen?
Danke!
lg Dominik
|
|
|
| |
|
Hallo,
> Danke für dein Programm. Habe jetzt auch beim mir den Fehler gefunden.
Ähem, eigentlich ist das dein Programm, nur eben um deine Fehler erleichtert.
Hättest du die Änderungen mal nachvollzogen (so war es eigentlich gedacht), dann wäre dir aufgefallen, dass ich eine weitere Variable eingeführt habe und nun zwischen origA und A unterscheide!
Das hat den Grund, dass die Eingabevariable nur lesbar ist und somit nicht veränderlich. Deshalb bekommst du einen Fehler, wenn du versuchst, ihr einen neuen Wert zuzuweisen (bei A := RT*A).
Also mach es wie ich: Kopier direkt am Anfang dein Original-A in eine lokale Variable A und schon ist der Spuk vorbei.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 15:15 Do 08.11.2007 | | Autor: | DominikW |
Hallo!
Ich habe jetzt alles mögliche probiert, aber ich schaffe es nicht das Programm fehlerfrei zu bekommen.
Mein Code ist:
with(LinearAlgebra);
origA := <<3 | 3 | 3>, <6 | 6 | 6>, <9 | 9 | 9>>;
rotation := proc(origA::Matrix) ... end;
global m,n,i,j,p,c,s, R, A;
A:= origA;
m := ColumnDimension(A);
n := RowDimension(A);
if (m = n) then
for j from 1 to n do
for i from j+1 to n-1 do
p := [mm] sign(A[j,j])*sqrt((A[j,j])^2+(A[i,j])^2);
[/mm]
c := (A[j,j])/p;
s := (-A[i,j])/p;
#R:=<<1|0|0>,<0|1|0>,<0|0|1>>;
R:= IdentityMatrix(n);
R[i,j] := -s;
R[j,i]:= s;
R[i,i]:= c;
R[j,j] := c;
A:=R^'%'T.A;
end do;
end do;
return A;
else
error "Keine quadratische Matrix!";
end if;
end proc;
Und beim Aufruf der Procedur kommt folgende Meldung:
rotation(origA);
Error, (in rotation) invalid assignment of non-zero to identity off-diagonal
Ich habe keine Ahnung wo da noch ein Fehler sein kann.
Bitte helft mir!
LG Dominik
|
|
|
| |
|
Hallo,
die Fehlermeldung sagt es doch schon:
Du versuchst, die Einträge außerhalb der Hauptdiagonalen einer Einheitsmatrix zu verändern, was standardmäßig verboten ist. Um das Problem zu umgehen, ersetze diese Zeile durch:
R:= IdentityMatrix(n, outputoptions=[shape=[], readonly=false]);
Was ich nicht verstehe, ist das "... end;" direkt hinter proc(). Ein Versehen?
Außerdem musst du die Apostrophe um das Prozentzeichen bei der transponierten Matrix entfernen:
A:=R^%T.A;
Ein weiterer Fehler: Die innere Schleife muss bis n laufen, nicht bis (n-1).
Und schließlich: Man muss die Funktion signum anstelle von sign benutzen. Warum, ist mir noch nicht ganz klar...
Gruß
Martin
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 17:20 Do 08.11.2007 | | Autor: | DominikW |
Hallo!
Möchte mich nur noch einmal bedanken!
Habe jetzt alles ausgebessert und das Programm funktioniert.
LG Dominik
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 20:52 Do 08.11.2007 | | Autor: | DominikW |
Hallo!
Habe die QR- Zerlegung jetzt mit Householder Transformation auch noch gemacht. Da hat mir Wikipedia sehr geholfen, da ein fast funktionsfähiges Programm auf der Homepage ist.
Habe jetzt aber das Problem, dass Maple zu ungenau oder zu genau rechnet und es so zu "Rundungsfehlern" kommt. An den Stellen wo Nullen stehen sollten, stehen Zahlen wie z.B. 3,18 * 10^-10.
Gibt es eine Möglichkeit eine Matrix auf 3 Kommastellen zu runden, damit diese minimalen Werte 0 werden.
Danke für eure Hilfe im Voraus!
LG Dominik
|
|
|
| |
|
Hallo,
> Habe die QR- Zerlegung jetzt mit Householder Transformation auch noch gemacht. Da hat mir Wikipedia sehr geholfen, da ein fast funktionsfähiges Programm auf der Homepage ist.
Hehe, das habe ich gesehen, aber nicht geprüft, ob es komplett ist.
> Gibt es eine Möglichkeit eine Matrix auf 3 Kommastellen zu runden, damit diese minimalen Werte 0 werden.
Klar. Es gibt die Umgebungsvariable Digits, deren Wert du vor oder innerhalb deiner Prozedur neu setzen kannst, z.B.:
Digits := 5:
Experimentier ruhig mit verschiedenen Werten, denn wenn du zu wenige Stellen nimmst, dann kann das Ergebnis am Ende deutlich anders aussehen, insbesondere mit großen Matrizen.
Alternativ kannst du die Matrix am Ende durch eine Funktion schicken, die alle Einträge unter einem Schwellwert einfach auf Null setzt.
Gruß
Martin
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 23:13 Do 08.11.2007 | | Autor: | DominikW |
hallo!
Danke für deine rasche Antwort. Jetzt funktioniert das auch.
LG Dominik
|
|
|
|
