SVD einer 2x2Matrix

18.07.2009, 09:40

_Nikolas

Auf diesen Beitrag antworten »

SVD einer 2x2Matrix

Hallo

Ich arbeite gerade an meiner Bachelorarbeit in der Informatik und habe ein mathematisches Problem, das ich gerne lösen würde.
Es geht dabei um das Problem zwei Punktwolken (im R^2) so zu verschieben und zu drehen, dass der durchschnittliche Abstand jedes Punktes zu seinem nächsten Nachbarn in der anderen Wolke minimiert wird.

Im R^3 kann man dafür die SVD einsetzen. Seien x'_i die Punkte der ersten Wolke und p'_i die der zweiten, wobei die Mittelpunkte beider Wolken zusammen fallen. Der Punkt x'_j soll nach der Rotation moeglichst nahe bei p'_j liegen.

In einem ersten Schritt berechnet man $\begin{eqnarray*} X = \sum x'_ip'i^t \end{eqnarray*}$ und bildet seine SVD, also X = U*S*V^t Eigenschaften von U,S,V . Nach einem Satz ist dann R=U*V^t die optimale Rotationsmatrix. In der Vorlesung, in der ich die Methode (insgesamt wird das ganze dann iterative closest point, ICP genannt) gehört habe, hiess es aber, dass der Satz für rk(X)=3 gilt. (kein genau dann!).

Ich würde ihn gerne im 2D einsetzen.

Als kleines Beispiel habe ich mir zwei Punktwölkchen definiert:
x'_i=(0,3),(-3,0),(0,-3),(3,0)
p'_i=(2,2),(-2,2),(-2,-2),(2,-2)

X = (12 12 // -12 12)

Wenn ich diese Matrix hier einsetze, bekomme ich zwei Matrizen als Ergebniss, die multipliziert die Rotationsmatrix für eine Drehung um pi/4 ergeben, also exakt das, was ich erwartet habe.

Das Problem ist aber, dass ich das natürlich selber ausrechnen will. (Und das möglichst schnell, weil meine Zielplattform für die ich entwickle, ziemlich langsam ist).

Bei der Wikipedia bin ich auf einen interessanten Artikel gestoßen , der eigentlich recht gut klingt aber nicht ganz passt. Dort soll X aus beliebigen (zentrierten) Vektoren bestehen, ich habe dafür erstmal mein X gewaehlt, weil der Artikel eigentlich nur die Erstellung der SVD fuer eine Matrix erklären sollte. (Ich glaube, hier liegt mein Fehler).

Wenn ich F = X*X^t ausrechne, erhalte ich für mein Beispiel (288 0 // 0 288) und da kann man die Eigenwerte und Vektoren recht schnell rauslesen und sehen, dass das bei mir nicht zum Ziel fuehrt. (Da auch F=G)

So. Abschließend noch mal Formuliert: Wie kann ich (am besten in einer geschlossenen Form) die SVD einer 2x2-Matrix berechnen?

// Alternative Wege zur Berechnung der Rotation sind natürlich auch gerne gesehen : )

Wär schön, wenn mir da jemand weiterhelfen könnte.

Nikolas

19.07.2009, 00:14

Mazze

Auf diesen Beitrag antworten »

Zitat:

So. Abschließend noch mal Formuliert: Wie kann ich (am besten in einer geschlossenen Form) die SVD einer 2x2-Matrix berechnen?

Es gibt Zahlreiche Bibliotheken die dieses tun. Falls Du nicht die mathematische Herleitung brauchst noch die Methode selbst implementieren sollst, so wäre es doch sinnvoll auf genau eine solche zu benutzen. Für java zum Beispiel eignet sich Jama. Für C++ kenne ich da CImg, ist zwar eigentlich Bildverarbeitung, aber da in der Bildverarbeitung die Singulärwertzerlegung ebenfalls wichtig ist, gibts diese Funktion dort auch.

Zitat:

Wie kann ich (am besten in einer geschlossenen Form) die SVD einer 2x2-Matrix berechnen?

Wenn Du es effektiv mit einem Rechner machen willst, (und es selbst implementieren willst) empfehle ich dir das Buch "Numerical Recipes". Es gibt sicher andere, aber dieses ist jenes welches ich kenne. Dieses Buch solltest Du in der entsprechenden Bibliothek finden können.

Von Hand lässt sich das natürlich bei einer 2 x 2 Matrix sehr leicht tun. Hier wird die Konstruktion beschrieben.

Zitat:

Bei der Wikipedia bin ich auf einen interessanten Artikel gestoßen , der eigentlich recht gut klingt aber nicht ganz passt. Dort soll X aus beliebigen (zentrierten) Vektoren bestehen, ich habe dafür erstmal mein X gewaehlt, weil der Artikel eigentlich nur die Erstellung der SVD fuer eine Matrix erklären sollte. (Ich glaube, hier liegt mein Fehler).

In dem Artikel wird eine Approximation der SVD behandelt, ich kenne diese nicht, wollte Dich aber darauf hinweisen das dies keine exakte SVD ist.

19.07.2009, 20:48

_Nikolas

Auf diesen Beitrag antworten »

Danke schon mal, dass du dir die Zeit genommen hast. Ich schreibe in C++ und da das Programm spaeter auf einer embedded Hardware laufen soll, wäre eine eigene Implementation, die die feste groesse der Matrix ausnutzt, die beste Loesung.
Die Berechnung ueber die Eigenvektoren von AA^t und A^tA kannte ich schon, bei meiner Matrix $\begin{eqnarray*} A=\begin{pmatrix} 12 & 12 \\ -12 & 12 \end{pmatrix} \end{eqnarray*}$ komme ich da auf $\begin{eqnarray*} AA^t = \begin{pmatrix} 12 & 12 \\ -12 & 12 \end{pmatrix} \begin{pmatrix} 12 & -12 \\ 12 & 12 \end{pmatrix}=\begin{pmatrix} 288 & 0 \\ 0 & 288 \end{pmatrix} (=A^tA) \end{eqnarray*}$
Damit ist dann U=V=1, was nicht dem Ergebniss der Berechnung hier entspricht, die auf dem Code aus den Numerical Recipes basiert und auch schön auf die richtige Rotationsmatrix kommt.

Siehst du da vielleicht einen Fehler?

19.07.2009, 21:00

Mazze

Auf diesen Beitrag antworten »

Dein Fehler ist das Du U = 1 setzt. Für die Singulärwertzerlegung gilt

$\begin{eqnarray*} A = U\Sigma V^T \end{eqnarray*}$

Wenn V und U = 1 wären , dann wäre $\begin{eqnarray*} A = U\Sigma V^T = \Sigma \end{eqnarray*}$ was offensichtlich nicht stimmt.

edit:

Die C++ Bibliothek CImg ist im Übrigen open Source, d.h Du kannst Dir den Code für die SVD dort genau anschaun.

19.07.2009, 21:20

_Nikolas

Auf diesen Beitrag antworten »

Das heisst, ich soll die Eigenvektoren noch so skalieren, dass die Gleichung erfüllt ist? Joa, darauf hätte ich auch kommen können...
Ganz allgemein habe ich doch vier Eigenwerte (w,x für AA^t und y,z für A^tA) und zwei SV (s1,s2), im Allgmeinen Fall habe ich dann doch die Gleichung
$\begin{eqnarray*} \begin{pmatrix} \alpha w_1 & \beta x_1 \\ \alpha w_2 & \beta x_2 \end{pmatrix} * \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix} *\begin{pmatrix} \gamma y_1 & \gamma y_2 \\ \delta z_1 & \delta z_2 \end{pmatrix} = \begin{pmatrix} a & b \\ c & d\end{pmatrix} \end{eqnarray*}$
zu loesen, oder?

(Da sitze ich gerade dran, wenn da ein Fehler ist, denn du siehst, könntest du den noch loswerden, bisher sieht die Gleichung irgendwie nach zu vielen Unbekannten aus...)

Edit:
Wenn ich das durchrechne, komme ich z.b. auf
$\begin{eqnarray*} \alpha w_1 \sigma_1 \gamma y_1 + \alpha \delta w_1 \sigma_1 z_1 = a \end{eqnarray*}$
Damit komme ich dann auf vier Gleichungen bei vier Skalierungen der EW. Also eigenlich loesbar.
Komme ich damit zum Ziel? Wenn ich jetzt durchrechne, kann ich doch davon ausgehen, dass die Skalierungen von Null verschieden sind. (NullVektor ist nie EigenVektor), der Rest kann aber alles Null sein. (Oder gibt es einen Satz, mit dem die EW von A^tA von Null verschieden sind?)

19.07.2009, 21:27

Mazze

Auf diesen Beitrag antworten »

Die Matrix V besteht aus den Eigenvektoren von $\begin{eqnarray*} A^TA \end{eqnarray*}$ und U besteht aus den Eigenvektoren von $\begin{eqnarray*} AA^T \end{eqnarray*}$ . Ich habe Dir oben einen Link zu genau diesem gepostet, da steht auch eine Alternative zu U wenn man V hat.

Anzeige

19.07.2009, 21:29

_Nikolas

Auf diesen Beitrag antworten »

okay, da habe ich jetzt beim Aufschreiben nicht exakt drauf geachtet.
Aber wenn ich das umstelle, sollte ich damit doch alles berechnen koennen, oder habe ich etwas uebersehen?

19.07.2009, 22:03

Mazze

Auf diesen Beitrag antworten »

Ehrlich gesagt sehe nicht so ganz wo Du diese Gleichung her hast.

19.07.2009, 22:06

_Nikolas

Auf diesen Beitrag antworten »

Ich habe eigentlich nur die Matrizen von links nach rechts miteinander multipliziert. Die linke Seite der Gleichung ist dann der 0/0-Eintrag des Produkts.

19.07.2009, 22:38

Mazze

Auf diesen Beitrag antworten »

Wenn überhaupt hätte ich das Gleichungssystem

$\begin{eqnarray*} \begin{pmatrix} u_{11} & u_{12} \\ u_{21} & u_{22} \end{pmatrix} * \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix} *\begin{pmatrix} v_{11} & v_{12} \\ v_{21} & v_{22} \end{pmatrix} = \begin{pmatrix} a & b \\ c & d\end{pmatrix} \end{eqnarray*}$

unter den Bedingungen

$\begin{eqnarray*} U^TU = I \end{eqnarray*}$
$\begin{eqnarray*} V^TV = I \end{eqnarray*}$

aufgestellt. Da U und V orthogonal sind ist die skalierung der Eigenvektoren sowieso irrelevant, da sie immer auf Norm 1 skaliert werden

19.07.2009, 22:50

_Nikolas

Auf diesen Beitrag antworten »

okay. Zwei kleine Fragen:
- Müsste V nicht noch transponiert werden?
- Häh? (Im Sinne von: Was meinst du mit dem Satz zur orthogonalität von U und V?) Wenn u_i und v_i normierte Eigenvektoren von AA^t und A^tA sind, welche Freiheitsgrade habe ich noch im Gleichungssystem?

19.07.2009, 23:09

Mazze

Auf diesen Beitrag antworten »

Zitat:

üsste V nicht noch transponiert werden?

Ist es nicht egal wie man die Variablen der Matrix nennt ? Augenzwinkern

Augenzwinkern

Wenn Du es transponierst dann ändert sich nur das Du die Lösung für $\begin{eqnarray*} v_{12} \end{eqnarray*}$ mit $\begin{eqnarray*} v_{21} \end{eqnarray*}$ austauschen musst. Die Lösung selber ändert sich nicht, nur die Bezeichner.

Zitat:

Wenn u_i und v_i normierte Eigenvektoren von AA^t und A^tA sind, welche Freiheitsgrade habe ich noch im Gleichungssystem?

Die Singulärwertzerlegung als solche ist nicht eindeutig. Daher sollte man wohl mindestens einen Freiheitsgrad bekommen.

19.07.2009, 23:23

_Nikolas

Auf diesen Beitrag antworten »

Meine Frage zu den Freiheitsgraden hat sich darauf bezogen, dass in der Gleichung keine Unbekannten mehr stehen. Die normierten EV sind fest gegeben, der Rest auch, so dass ich eigentlich nichts mehr zu lösen habe.

Wenn die SVD nicht eindeutig ist, ist dann wenigsten U*V^t eindeutig? Nach dem im ersten Post erwähnten Satz würde ich mich sehr wundern, müsste U*V^t doch eindeutig festgelegt sein.

Danke für deine Hilfe smile

smile

19.07.2009, 23:32

Mazze

Auf diesen Beitrag antworten »

Hier

findest Du ein Beispiel das zeigt das die SVD nicht eindeutig ist.

Zitat:

Die normierten EV sind fest gegeben, der Rest auch, so dass ich eigentlich nichts mehr zu lösen habe.

Das würde dann genau zu dem führen was ich schon oben erwähnt habe. Wenn ich mich aber nicht ganz irre umgeht man den Weg mit der Eigenwertberechnung durch einen numerisch stabileren Weg. Näheres dazu dürfte in numerical recipes stehen.

19.07.2009, 23:37

_Nikolas

Auf diesen Beitrag antworten »

Was hast du dann mit der Gleichung gemeint? Ich dachte, dass mein ursprüngliches Problem in der Skalierung der EV liegen würde.

Wie soll ich denn die (eine) SVD der obigen Matrix berechnen?

20.07.2009, 18:06

Mazze

Auf diesen Beitrag antworten »

Zitat:

Wie soll ich denn die (eine) SVD der obigen Matrix berechnen?

Ich würde einfach V mit den normierten Eigenvektoren von $\begin{eqnarray*} A^TA \end{eqnarray*}$ füllen und U entsprechend orthogonal Ergänzen. Es ist

$\begin{eqnarray*} A^TA = \begin{pmatrix} 288 & 0 \\ 0 & 288 \end{pmatrix} \end{eqnarray*}$

Die Eigenwerte sieht man sofort, das Eigenvektorgleichungssystem ist dann

$\begin{eqnarray*} (288 I - A^TA)x = 0 \Leftrightarrow 0x = 0 \end{eqnarray*}$

Das heisst der Eigenraum zu 288 ist zweidimensional. Dort wählen wir eine orthonormale Basis, und zwar

$\begin{eqnarray*} v_1 = -e_1, v_2 = -e_2 \end{eqnarray*}$ (weil es so schön einfach ist Augenzwinkern

Augenzwinkern

)

Nun setzen wir

$\begin{eqnarray*} u_1 = \frac{Av_1}{\sqrt{288}} = \frac{1}{\sqrt{288}}\begin{pmatrix} -12 \\ 12 \end{pmatrix} \end{eqnarray*}$

und

$\begin{eqnarray*} u_2 = \frac{Av_2}{\sqrt{288}} = \frac{1}{\sqrt{288}}\begin{pmatrix} -12 \\ -12 \end{pmatrix} \end{eqnarray*}$

Mit $\begin{eqnarray*} U = (u_1,u_2) \end{eqnarray*}$

und

$\begin{eqnarray*} \Sigma = \begin{pmatrix} \sqrt{288} & 0 \\ 0 & \sqrt{288} \end{pmatrix} \end{eqnarray*}$

bekommst Du dann

$\begin{eqnarray*} A = U\Sigma V^T \end{eqnarray*}$

und Du hast genau die Zerlegung die das Javascript berechnet hat. Allerdings hätte man für V auch eine andere orthonormale Basis wählen können.

20.07.2009, 19:23

_Nikolas

Auf diesen Beitrag antworten »

Okay, das sieht interessant aus. Ich werde mich gleich mal dransetzen, wenn ich's damit hinbekomme, gibt's mindestens ne Fussnote in der Arbeit smile

smile

21.07.2009, 02:44

_Nikolas

Auf diesen Beitrag antworten »

So weit so gut. Ich habe den Algorithmus implementiert und bekomme auch für das Beispiel einen Drehwinkel von 45 grad raus.

Ein bischen sorge ich mich aber um die Berechnung der u_i. Da AtA symmetrisch ist, habe ich schon mal nur reelle Eigenwerte, nur kenne ich keinen Satz, der EW_i=0 verbietet. In diesem Fall koennte ich U nicht wie beschrieben berechnen.

Könntest du mir sagen, ob ich da in Probleme reinlaufen werde?

Vielen Dank auf jeden Fall, du hast mich da wirklich weitergebracht Freude

Freude

21.07.2009, 13:37

Mazze

Auf diesen Beitrag antworten »

Das steht auch in dem Link den ich Dir gepostet habe. Du bildest dann für alle $\begin{eqnarray*} \sigma_i \end{eqnarray*}$ die entsprechenden Vektoren $\begin{eqnarray*} \frac{1}{\sigma_i}(AV)_i \end{eqnarray*}$ und füllst die übrigen Spalten von U so auf das U eine orthogonale Matrix wird. Die Matrix V kannst Du immer bestimmen, da sie aus einer Basis der Eigenvektoren von $\begin{eqnarray*} A^TA \end{eqnarray*}$ besteht.

21.07.2009, 14:17

_Nikolas

Auf diesen Beitrag antworten »

Okay, so etwas in die Richtung hatte ich mir gedacht. (Leider wird dort aber r nie definiert, so dass man es ein bischen suchen muss).
Dann baue ich das noch ein.

Vielen Dank für deine ausführliche Unterstützung. Ich hoffe mal, dass ich die Erweiterung hinbekommen sollte smile

smile

23.07.2009, 00:18

_Nikolas

Auf diesen Beitrag antworten »

So. Die Implementation ist fertig, ich habe jetzt zufaellige Punkte in der Ebene genommen, sie um den Ursprung zentiert, rotiert und versucht die Rotation wieder uber die SVD auszulesen und es funktioniert wirklich smile

smile

Vielen Dank, das wir sicherlich ein nettes Kapitel in der Arbeit. Freude

Freude

Kennst du zufälligerweise auch noch eine zitierfähige Quelle für den Algorithmus? Eine Seite im Netz macht sich in der Literaturliste nicht allzu gut und selber herleiten will ich das ganz sicher nicht.

23.07.2009, 01:19

_Nikolas

Auf diesen Beitrag antworten »

Anscheinend darf man hier nur 15min lang bearbeiten, auch wenn noch weiterer Beitrag da ist.

Also noch ein dritter von mir:
In Sigma stehen doch die Wurzeln der Eigenwerte von M=A^tA. (Durch die ich dann später dividiere).
Warum darf ich davon ausgehen, dass diese Werte reel sind? Da M symmetrisch ist, sind die Eigenwerte reell, aber nicht unbedingt größer als Null. Bei der implementierung habe ich an der Stelle nur zwei asserts stehen, so dass die EW wirklich größer als Null sind, aber warum das der Fall ist, kann ich nicht erklären.

23.07.2009, 18:08

Mazze

Auf diesen Beitrag antworten »

Zitat:

Da M symmetrisch ist, sind die Eigenwerte reell, aber nicht unbedingt größer als Null.

Die Matrix $\begin{eqnarray*} A^TA \end{eqnarray*}$ ist positiv semidefinit. Das siehst Du daran

$\begin{eqnarray*} x^TA^TAx = (Ax)^T(Ax) = ||Ax||_2^2 \geq 0 \end{eqnarray*}$

Symmetrische, positiv semidefinite Matrizen haben nur Eigenwerte größergleich 0. Daher haut das hin Augenzwinkern

Augenzwinkern

. Wenn A invertierbar ist (und nur dann) gilt sogar das A^TA positiv definit ist, dann sind die Eigenwerte sogar echt größer als 0.

17.08.2009, 12:51

Nikolas

Auf diesen Beitrag antworten »

Der Code ist fertig und im aktiven Einsatz : )

Das Kapitel für die Bachelorarbeit nimmt langsam gestalt an, wer Interesse am Code oder der getexten Herleitung hat, kann sich gerne per PM melden.

Was jetzt nur fehlt, ist eine zitierfähige Quelle für den Algorithmus...

17.08.2009, 15:26

Mazze

Auf diesen Beitrag antworten »

Wenn Du es über die reine Eigenwertberechnung machst (also die teure Variante), dürfte Dir ein bel. lineare Algebra Buch als Quelle schon reichen. In den Numerical Recipeis ist eine numerisch bessere Methode beschrieben.

18.11.2011, 12:31

Goldnas

Auf diesen Beitrag antworten »

Zitat:

Original von Nikolas
Der Code ist fertig und im aktiven Einsatz : )

Das Kapitel für die Bachelorarbeit nimmt langsam gestalt an, wer Interesse am Code oder der getexten Herleitung hat, kann sich gerne per PM melden.

Was jetzt nur fehlt, ist eine zitierfähige Quelle für den Algorithmus...

Hätte es brauchen können....

LG Goldnas

1

Verwandte Themen

Die Beliebtesten »

Die Größten »

Die Neuesten »