M'=P^T·M·P nach P auflösen

[DrStupid]

Hallo liebe Freunde der Matrizenrechnung,

Ich habe zwei Mengen von Punkten, von denen ich weiß, dass sie durch Rotationen um den Koordinantenursprung aufeinander abgebildet werden können. Dazu muss ich die zugehörige Rotationsmatrix finden. Leider haben die beiden Punktmengen nicht die gleiche Reihenfolge. Ich habe erst einmal zwei Matrizen X und X’ definiert, deren Spaltenvektoren, die Koordinaten der beiden Punktmengen enthalten. Damit gilt

X’ = R·X·P

wobei R die gesuchte Rotationsmatrix und P eine Übergangsmatix ist, die die Spaltenvektoren permutiert. Würde ich P kennen, könnte ich R aus den Punktkoordinaten berechnen:

R = X’·P^-1·X^T·(X·X^T)^-1

Aber leider ist P unbekannt. Da R orthogonal ist, komme ich mit

M := X^T·X

zumindest schon mal zu einer Gleichung, die nur noch P als Unbekannte enthält:

M’ = P^T·M·P

Aber dann bin ich mit meinem Latein am Ende. Kann man P aus M und M’ berechnen und wenn ja wie?

Gruß,
DrStupid

[quazee]

Hallo Dr. Stupid,

der Haken an der Sache ist, dass im Allgemeinen die Permutationsmatrix und damit die Rotaionsmatrix nicht eindeutig sind. Um ein kleines, wenn auch konstruiertes Beispiel zu geben: angenommen du hast vier Punkte, bei (1, 0), (-1, 0), (0, 1), (0, -1). Je nachdem wie ich meine Permutaionsmatrix wähle, ist es möglich, dass ich meine Daten um 90°, 180°, 270° oder 360° drehen muss.
Das Beispiel soll nur zeigen, dass man Asymmetrien der Daten ausnutzen muss um eindeutige Ergebnisse zu erhalten (wenn die Daten wie im Beispiel symmetrisch sind, hast du wohl Pech). In der Matrix M stehen ja die Skalarprodukte deiner Punkte:

m_{ij} = \langle x_i, x_j \rangle

Ich würde wahrscheinlich als erstes schauen ob ein Betragsquadrat auf der Diagonalen nur einmal vorkommt (z.B das Maximum), falls alle mehrmals vokommen, gibt es zwei Punkte mit gleichem Betrag, deren Skalarprodukt nur einmal vorkommt, etc. Nicht besonders schön, aber etwas besseres fällt mir grad nicht ein…

Viele Grüße

[Marco_F_a4ea4d]

Hi,

mir fällt nur die Brute-Force-Methode ein. Eine Matrix, die Zeilen oder Spalten vertauscht, ist einfach nur die Einheitsmatrix, bei der die entsprechenden Zeilen oder Spalten vertauscht sind wie im Ergebnis. Es gibt also endlich viele Möglichkeiten für P. Ich würde also ein kleines Programm schreiben, das alle Möglichkeiten durchprobiert.

Viele Grüße

Marco

[DrStupid]

mir fällt nur die Brute-Force-Methode ein.

Genau das will ich vermeiden.

Es gibt also endlich viele Möglichkeiten für P

Bei n Punkten sind es n! Möglichkeiten und das kann ziemlich schnell ziemlich viel werden.

Da probiere ich es lieber mit meinem Plan B: Minimierung der Summe der minimalen Abstände mit irgend einem geeigneten Optimierungsverfahren.

[Lutz_Lehmann]

Hi,

bestimme die Trägheitsmatrizen X’*X’^T und X*X^T, diese sind von der Reihenfolge unabhängig, denn P*P^T=I. Eine orthogonale Transformation erhältst Du nun, wenn Du beide Matrizen auf Eigenwertform bringst.

Oder einen Schritt früher: Bestimme die Singulärwertzerlegungen von X und X’. Wenn alle Singulärwerte unterschiedlich sind, dann lässt sich aus den linken Faktoren die Transformation ableiten, mit den rechten Faktoren läßt sich testen, ob der Unterschied tatsächlich nur eine Permutation war.

Gruß, Lutz

[quazee]

Hey,

ich weiß zwar nicht was eine Singulärwertzerlegung ist, aber was genau meinst du mit

diese sind von der Reihenfolge unabhängig

Es gilt nämlich nicht M’ = X’*X’^T = X*X^T = M
P vertauscht die Spalten, P^T nochmal die Zeilen, ist pi die zugrunde liegende Permutation, dann sind die Einträge von M’

m’_{i j} = m_{\pi(i) \pi(j)}

Viele Grüße

[quazee]

Hallo nochmal,

meine Güte bin ich durcheinander… hab mal meine Murksantworten gelöscht, jetzt nochmal strukturert.
Angenommen es gibt ein x ind X, dessen Norm nur einmal vorkommt. Dann kann man einfach ein entsprechendes x’ mit gleicher Norm finden und weiß, dass diese durch Drehung aufeinander abgebildet werden.
Problem: wir kennen weder Winkel noch Achse
Aber: wir wissen, dass sich ein Punkt bei einer Drehung innerhalb einer Ebene senkrecht zur Drehachse bewegt. Wird x auf x’ abgebildet, liegt der Verbindungsvektor in einer Ebene, die folgende Bedingung erfüllt

\langle x - x’, n \rangle = 0

wobei n die Drehachse ist. Nun kennen wir n nicht, sondern x-x’. Wir können also die Gleichung ebenso als Ebene auffassen, welche durch den Normalenvektor x-x’ bestimmt wird und besagt, in welcher Ebene mögliche Drehachsen liegen.
Nun hilft uns das nicht wirklich weiter, es sei denn wir haben in X einen weiteren Punkt y, welchen wir eindeutig einem Punkt y’ zuordnen können. Auch hier wird durch y-y’ eine Ebene bestimmt. Diese beiden Ebenen der Drehachsen stehen senkrecht zur Ebene, welche durch x-x’ und y-y’ aufgespannt wird (sofern x-x’ und y-y’ nicht parallel sind). Insbesondere liegt die Drehachse

n = (x-x’) \times (y-y’)

in der Schnittgeraden. Fehlt noch der Winkel. Dafür müsen wir x und x’ entlang n in die Ebene = 0 projezieren. Die Vektoren deren Winkel wir messen müssen sind dann

x-n’ \text{ und } x’-n’ \text{ mit } n’ := \langle x, \hat n \rangle \hat n, \quad | \hat n | = 1

und das macht man bekanntlich mit

\varphi = \arccos \left( \frac{\langle x-n’, x’-n’ \rangle}{|x-n’| | x’-n’|} \right)

Falls du dann immernoch die Permutationsmatrix brauchen solltes hab ich auch noch ne Idee.

Viele Grüße

[Moritz_fe1b4f]

Hallo,

Ich habe zwei Mengen von Punkten, von denen ich weiß, dass sie
durch Rotationen um den Koordinantenursprung aufeinander
abgebildet werden können. Dazu muss ich die zugehörige
Rotationsmatrix finden. Leider haben die beiden Punktmengen
nicht die gleiche Reihenfolge.

Wie viele Dimensionen hast du denn?

Wenn es nur zwei sind, wuerde ich als nicht-Mathematiker eine Art Korrelationsfunktion bauen. Dazu die Ortsvektoren nach dem Winkel sortieren, und dann die Summe der Abstandsquadrate berechnen, also

K(v) = Sum_i (A(i) - B((i + v) MOD n))^2

wobei A_i und B_i die Vektoren vor und nach der Rotation sind, n die Anzahl der Vekotren und i von 0 bis n-1 laeuft. Aus dem v, fuer das K minimal ist, kannst du dann den Rotationswinkel und damit die Rotationsmatrix bestimmen. Sind „nur“ O(n^2) Operationen; wenn du das ganze auf mehr Dimensionen verallgemeinern willst, faengt es vermutlich an, deutlich haesslicher zu werden.

Viele Gruesse,
Moritz

[Lutz_Lehmann]

ich weiß zwar nicht was eine Singulärwertzerlegung ist, aber

http://lmgtfy?q=Singulärwertzerlegung

Aber vielleicht kennst Du die Kahunen-Loeve (?) Transformation oder die Hauptachsentransformation,…

was genau meinst du mit

diese sind von der Reihenfolge unabhängig

Es gilt nämlich nicht M’ = X’*X’^T = X*X^T = M

Nein, natürlich nicht, M und M’ sind nur ähnlich, die orthogonale Ähnlichkeitstransformation ist ja gerade gesucht.

P vertauscht die Spalten, P^T nochmal die Zeilen, ist pi die
zugrunde liegende Permutation, dann sind die Einträge von M’

m’_{i j} = m_{\pi(i) \pi(j)}

Das möchtest Du bitte nochmal genauer nachprüfen. P wird gerade „wegsummiert“

m’_{i,j}

\sum_{k,l,m,n,o} r_{i,k} x_{k,l} \delta_{l,\pi(m)}\delta_{\pi(m),n} x_{o,n} r_{j,o}
=\sum_{k,l,o} r_{i,k}x_{k,l}x_{o,l}r_{j,o}
=\sum_{k,o} r_{i,k}m_{k,o}r_{j,o}

Gruß zurück, Lutz

[DrStupid]

Hallo,

Erstmal vielen Dank für die Antworten. Ich werde die verschiedenen Lösungsvorschläge mal durchprobieren. Inzwischen habe ich auch eine Lösung gefunden:

Ich erzeuge aus den Punktmengen unbekannter Reihenfolge neue Punktmengen mit bekannter Reihenfolge, z.B.

A_i = Xm + Σ[(X_j-Xm)·|X_j-Xm|ⁱ]/Σ|X_j-Xm|ⁱ

(mit Xm=Mittelwert von X) und kann dann die Rotationsmatrix gemäß

R = A’·A^T·(A·A^T)^-1

berechnen. Wenn ich wollte, dann könnte ich anschließend die Permutationsmatrix ermitteln:

P_ij = exp[-(X_i-R^-1·X’_j)²/σ²]

Aber die brauche ich jetzt gar nicht mehr.

Allerdings habe ich jetzt ein neues Problem: Was mache ich, wenn die Punktmengen unvollständig sind, wenn also einige Spaltenvektoren von X in X’ fehlen oder umgekehrt? Hier fällt mir außer Probieren (was wieder zu beliebig langen Rechenzeiten führen kann) wieder nur die Minimierung von ΣMin{(X_i-R·X’_j)²} (eventuell noch mit Herausstreichen von Ausreißern) ein und das ist auch nicht gerade trivial, weil es da unzählige lokale Minima gibt.

[DrStupid]

bestimme die Trägheitsmatrizen X’*X’^T und X*X^T, diese sind
von der Reihenfolge unabhängig, denn P*P^T=I.

Mit X’·X’^T=R·X·X^T·R^T habe ich im Moment das gleiche Problem wie mit X’^T·X’=P^T·X’^T·X’·P, weil ich nicht weiß, ob und wenn ja wie man anhand von Matrizengleichungen der Form A=C·B·C^T bzw. A=C·B·C^-1 von A und B auf C schließen kann.

Eine orthogonale
Transformation erhältst Du nun, wenn Du beide Matrizen auf
Eigenwertform bringst.

Mal sehen ob mich das Stichwort weiterbringt.

Oder einen Schritt früher: Bestimme die
Singulärwertzerlegungen von X und X’. Wenn alle Singulärwerte
unterschiedlich sind, dann lässt sich aus den linken Faktoren
die Transformation ableiten, mit den rechten Faktoren läßt
sich testen, ob der Unterschied tatsächlich nur eine
Permutation war.

Als Lösung für das ursprüngliche Problem hört sich das nicht sehr verlockend an, weil die Singulärwertzerlegungen wohl nur iterativ möglich ist (zumindest nach dem zu urteilen, was ich bis jetzt darüber gelesen habe) während ich die gesuchten Matrizen nach meiner oben beschriebenen Methode explizit berechnen kann. Wenn man auf diese Weise aber testen könnte, ob wirklich nur eine Permutation vorliegt, dann ließe sich so möglicherweise das zweite Problem lösen – nämlich die gezielte Eliminierung von nicht paarweise vorliegenden Punkten. Das werde ich mir also näher ansehen müssen. Notfalls habe ich da aber auch schon einen anderen Lösungsansatz - wenn auch keinen besonders elegenten.

[DrStupid]

Angenommen es gibt ein x ind X, dessen Norm nur einmal
vorkommt. Dann kann man einfach ein entsprechendes x’ mit
gleicher Norm finden und weiß, dass diese durch Drehung
aufeinander abgebildet werden.

Das ist zwar theoretisch möglich und ich habe zuerst selbst daran gedacht, die Punkte anhand ihrer Norm (bzw. der Norm ihres Abstands vom Mittelwert) zu sortieren, aber praktisch ist das nicht so ratsam weil reale Werte nicht beliebig genau sind. Dadurch kann es leicht passieren, dass Punkte, deren Normen sich im Rahmen der Messgenauigkeit kaum unterscheiden, vertauscht werden. Aber vielleicht ist sowas für ausgewählte Punkte praktikabel, die anhand irgend eines Kriteriums klar zu identifizieren sind.

[DrStupid]

Dazu die Ortsvektoren nach dem
Winkel sortieren, und dann die Summe der Abstandsquadrate
berechnen

Dabei gibt es das gleiche Problem wie bei der Sortierung nach der Norm: Wenn die Winkel im Rahmen der Messgenauigkeit gleich sind, kann man zur falschen Reihenfolge und somit zur falschen Transformation kommen. Dieses Problem tritt nicht auf, wenn man die Summe der minimalen Abstandsquadrate berechnet:

Σ{Min[(X_i-X’_j)²]}

und so was war tatsächlich mein Notfallplan für den Fall, dass es keine explizite Lösung gibt. Vielleicht verwende ich die nichtlineare Optimierung am Ende für den Feinschliff. Allerdings habe ich es mit Rotationen um bis zu drei Achsen zu tun und da kann ein einfaches Suchverfahren schon ziemlich mühselig werden – insbesondere, wenn es viele lokale Minima gibt.

[Lutz_Lehmann]

Hi,

um iterative Verfahren wirst Du nicht umhinkommen, und so hoch ist der Aufwand nicht. Proportional zu Anzahl der Dezimalstellen und Aufwand einer Matrixmultiplikation.

Zur Lösung von A=C·B·C^T:

Bestimme die Jordan-Normalformen A=U·S·U^T und B=V·T·V^T, mit U, V orthonormal. Da das hier 3x3-Matrizen sind, kann man diese Zerlegungen fast noch „per Hand“ bestimmen, wobei die Nullstellen des charakteristischen Polynoms doch besser mit einem numerischen Verfahren bestimmt werden.

Sind die Diagonalen der Diagonalmatrizen S und T nach Größe geordnet, so sollten sie übereinstimmen. Sind auch noch alle Diagonalelemente unterschiedlich, so muss U=C·V gelten oder C=U·V^T.

Sind die Punktmengen nicht nur rotiert, sondern auch noch verschoben worden, so sollte man vorher noch jeweils den Mittelpunkt bestimmen und von jeder Punktwolke abziehen, so dass beide um den Nullpunkt zentriert sind.

Gruß, Lutz