3. Lineare Regression im 3D

Lineare Regression im 3D

  • I

    Moin,

    ich bin auf der Suche nach C++ Code / einer Library, die mir zu einer gegebenen Punktewolke eine Ebene fittet.

    Im 2D ist dies kein Problem und man findet ohne Ende Beispiele. Im 3D bin ich allerdings noch nicht auf eine hilfreiche Lösung gestoßen.

    Für folgendes Einsatzgebiet brauche ich das:
    Ich habe eine Simulation entwickelt, in dieser kollidieren zwei konkave Objekte. Nun möchte ich die Größe der Kollisionsfläche zwischen den beiden bestimmen.

    Im ersten Schritt vergleiche ich die Punkte der beiden Objekte und speichere die Eckpunkte, dessen Abstand zueinander kleiner als ein der gegebene Wert d ist.
    Ich habe nun zwei Punktewolken (für jedes Objekt einen) die sehr nahe beieinander liegen und möchte durch diese eine Ebene fitten.

    Das Ganze sollte also so ähnlich sein wie hier: http://commons.wikimedia.org/wiki/File:Mincer_linreg.png
    Allerdings kann ich dazu keine passenden Code finden.

    Kennt vielleicht einer eine passende Library oder Code, der dies berechnet?

    Viele Grüße

    0
  • C

    Dieser Thread wurde von Moderator/in SeppJ aus dem Forum C++ (auch C++0x und C++11) in das Forum Rund um die Programmierung verschoben.

    Im Zweifelsfall bitte auch folgende Hinweise beachten:
    C/C++ Forum :: FAQ - Sonstiges :: Wohin mit meiner Frage?

    Dieses Posting wurde automatisch erzeugt.

    0
  • ?

    Siehe bei http://de.wikipedia.org/wiki/Lineare_Regression unter "Multiple Regression". Du willst ja die Funktion a*x1 + b*x2 + c an eine Punktwolke fitten, würde ich behaupten?
    Solange Du eine Bibliothek hast, die Matrizen transponieren, multiplizieren und invertieren/LU-zerlegen kann, ist das Berechnen der Koeffizienten ein Einzeiler, eben genau die Least-Squares-Formel aus dem Artikel.

    0
  • ?

    Du kannst ein Problem auf die folgende Form bringen:

    Suche x∈Rnx \in \mathbb{R}^n, so dass
    ∥Ax−b∥=minx~∈Rn∥Ax~−b∥\| A x - b \| = \min_{\tilde x \in \mathbb{R}^n} \| A \tilde x - b \|
    gilt, wobei A∈Rm×nA \in \mathbb{R}^{m \times n} und b∈Rmb \in \mathbb{R}^m.

    Das kannst du dann mit einem linear Least Squares Verfahren lösen (https://en.wikipedia.org/wiki/Linear_least_squares_%28mathematics%29).

    Besser als der obige Vorschlag ist es, die QR-Zerlegung zu verwenden. Löst man mit Hilfe der Normalengleichung, dann quadriert man die Kondition des Problems.

    0
  • ?

    Den QR-Pfad habe ich auch durchgeschaut...
    Im Wiki-Artikel steht geschrieben, dass man Q selten explizit braucht, kann mir jemand erklären wie das gemeint ist? Wenn ich die QR-Zerlegung verwende, dann komme ich doch auf

    Ax = QRx = b <=> x = (QR)^-1 * b = R^-1 * Q^T * b

    also für den simpelsten Anwendungsfall benötige ich doch Q?
    Ich habe in letzter Zeit Spaß dran gefunden, die ganzen numerischen Methoden mal in meine kleine Mathebib aufzunehmen, QR werde ich wohl auch mal hinzufügen. SVD brauche ich wohl auch zum entkoppeln von DGLs. Schön dass es immer spannende Dinge zu tun gibt 😉

    0
  • ?

    Hallo,

    in dem Fall der Regression lässt sich

    Ax=bA x = b

    nicht erfüllen, da wir ein überbestimmtes Gleichungssystem haben. Nur noch mal des Kontextes wegen, wie man ein linear least squares Problem mit der QR-Zerlegung löst, basier auf der Tatsache, dass

    ∥Ax−b∥=∥QRx−b∥=∥Q∗QRx−Q∗b∥=∥Rx−Q∗b∥\| A x - b \| = \| Q R x - b \| = \| Q^* Q R x - Q^* b \| = \| R x - Q^* b \|

    da die Multiplikation mit einer unitären Matrix die Norm nicht ändert (∥Q∗x∥=∥x∥\| Q^* x \| = \| x \|).

    Die Matrix RR hat die Form R=(R~0)R = \begin{pmatrix} \tilde R \\ 0 \end{pmatrix}. Das heißt man erhält da Minimum von ∥Rx−Q∗b∥\| R x - Q^* b \|, in dem man den oberen Teil des Gleichungssystemes Löst. D.h. teilt man Q∗b=(z_1z_2)Q^* b = \begin{pmatrix} z\_1 \\ z\_2 \end{pmatrix} analog zu der Aufteilung von RR auf, so ist x=R~−1z1x = {\tilde R}^{-1} z_1 . Da die Matrix R~\tilde R eine obere Dreieckmatrix ist, lässt sich dieses System durch einfaches Rückwärtseinsetzen lösen.

    Zu deiner Frage: Die Matrix QQ wird nicht explizit gebraucht, da du nur in der Lage sein musst, einen Vektor mit der Matrix QQ zu multiplizieren. D.h. alles was benötigt ist, ist eine Funktion, die Q∗bQ^* b berechnet. Das geht aber ohne QQ explizit aufzustellen.

    0
  • ?

    Jetzt fiele mir, dass ich bei der Erzeugung von R die Transformationen, die ich auf A anwende auch entsprechend auf b anwende und am Ende R und Q^Tb habe.
    Aber wenn ich die Zerlegung auf mehrere b's anwende, dann brauche ich die von Dir genannte Funktion und wenn ich die habe, habe ich doch praktisch auch die Matrix selber?
    Oder anders gesagt, wie kann die Funktion die b auf Q^T    b abbildet einfacher zu erzeugen sein als Q^T selber?

    Selbstverständlich rede ich hier über die einfache Zerlegung nicht im Sinne der überbestimmten Regressionsprobleme.

    0
  • ?

    Hallo,

    in der Regel wirst du die QR-Zerlegung mit Householder-Transformationen durchführen. Jede Transformation wird durch einen Vektor beschrieben. Indem du die Vektoren speicherst, kannst du die Matrix QQ auf einen Vektor bb Anwenden, indem du die Transformationen der Reihe nach auf den Vektor bb anwendest, ohne die Matrix QQ dabei explizit aufgestellt zu haben.

    0
  • O

    Mir ist nebenbei noch nie ein LSQ-Problem untergekommen, das nach Jacobi-Preconditioning so schlecht konditioniert war, das man die QR-Zerlegung brauchte. Deswegen befürworte ich immer die (unvollständige) Cholesky-Zerlegung

    0
  • ?

    Wie löst man denn mit unvollständiger Cholesky-Zerlegung ein LGS?

    0
  • O

    Wie löst du denn in LGS das nicht vollen Rang hat numerisch stabil?

    0
  • I

    Na gut dann werde ich mal ein bisschen programmieren.

    Kurze Frage vorweg:
    Welche der drei Raumrichtungen ich als Respone Y oder als Kovariablen x_1 und x_2 setze sollte keine Rolle spielen?

    0
  • ?

    otze schrieb:

    Wie löst du denn in LGS das nicht vollen Rang hat numerisch stabil?

    Garnicht, weil es die Lösung nicht gibt. D.h. um die Frage zu beantworten, muss man sagen, was man jetzt als Lösung betrachtet.

    Ich kenne als unvollständiger Cholesky-Zerlegung das hier. Und da sehe ich gerade nicht, was das mit dem Thema zu tun hat.

    icewater schrieb:

    Welche der drei Raumrichtungen ich als Respone Y oder als Kovariablen x_1 und x_2 setze sollte keine Rolle spielen?

    Das ist egal. Mir ist gerade nur aufgefallen, dass du von einer Parameterdarstellung der Ebene ausgehen musst. Bei der Normalendarstellung kommt kein lineares least squares mehr raus.

    0
  • I

    Ah jetzt wo du es sagt merke ich das meine bisherige Lösung totaler Quatsch ist.

    Ich suche ja praktische 9 Unbekannte und nicht nur 3.

    Also ich habe das Ganze jetzt mithilfe von Boost hiernach implementiert.
    http://www.codeproject.com/Articles/566326/Multi-Linear-Regression-in-Java?msg=4590356

    Nur mein b ist eben nicht (x,y,z)T,(x, y, z)^T,
    sondern (n_x,n_y,n_z,p_x,p_y,p_z,q_x,q_y,qz)T(n\_x, n\_y, n\_z, p\_x, p\_y, p\_z, q\_x, q\_y, q_z)^T

    somit muss X dann auch 9 statt 3 Spalten haben.
    Und wie mache ich das jetzt? Ich kann ja nicht 7 Spalten auf 1 setzen und nur die letzten beiden mit den x und y Werten füllen.
    Wie fülle ich also X, wenn Y die Z-Koordinaten enthält?

    0
  • I

    Em ja, ich habe ja jetzt die Normalenform raus. Mit welchen Matrizen würde ich denn gleich eine Parameterdarstellung rausbekommen?

    0
  • O

    ProgChild schrieb:

    otze schrieb:

    Wie löst du denn in LGS das nicht vollen Rang hat numerisch stabil?

    Garnicht, weil es die Lösung nicht gibt. D.h. um die Frage zu beantworten, muss man sagen, was man jetzt als Lösung betrachtet.

    Typischerweise wird aus theoretischen Überlegungen der Vektor aus dem span der möglichen Lösungen verwendet, welcher die minimale 2-Norm hat.

    incomplete cholsky zerlegung war der falsche Begriff. Passt dir Rank-Revealing Cholesky Factorization besser?

    0
  • ?

    Hallo,

    icewater schrieb:

    Em ja, ich habe ja jetzt die Normalenform raus. Mit welchen Matrizen würde ich denn gleich eine Parameterdarstellung rausbekommen?

    okay... Gehen wir davon aus, dass zz der an der Stelle (x,y)T(x ,y)^T gegebene Wert ist.

    Also... gesucht ist eine Gleichung der Form

    z=p+n⃗⋅(xy)z = p + \vec n \cdot \begin{pmatrix} x \\ y \end{pmatrix}.

    Das können wir als

    z=p+n_1x+n_2yz = p + n\_1 x + n\_2 y

    schreiben. Daraus ergibt sich

    z=(1xy)⋅(pn_1n_2)z = \begin{pmatrix} 1 & x & y \end{pmatrix} \cdot \begin{pmatrix} p \\ n\_1 \\ n\_2 \end{pmatrix}

    Hast du jetzt mehrere Messwerte, dann willst du, dass

    ∥b−Ax~∥\| b - A \tilde x \| minimal wird, wobei

    A=(1x_1y_11x_2y_2⋮⋮⋮)x~=(pn_1n_2)undb=(z_1z_2⋮)A = \begin{pmatrix} 1 & x\_1 & y\_1 \\ 1 & x\_2 & y\_2 \\ \vdots & \vdots & \vdots \end{pmatrix} \quad \tilde x = \begin{pmatrix} p \\ n\_1 \\ n\_2 \end{pmatrix} \quad \text{und} \quad b = \begin{pmatrix} z\_1 \\ z\_2 \\ \vdots \end{pmatrix}

    ist.

    otze schrieb:

    Typischerweise wird aus theoretischen Überlegungen der Vektor aus dem span der möglichen Lösungen verwendet, welcher die minimale 2-Norm hat.

    Okay... Das geht für den Fall, dass es eine Lösung gibt. Wenn die rechte Seite nicht im Bild der Matrix liegt, kommt man damit aber nicht weiter, oder übersehe ich was?

    otze schrieb:

    incomplete cholsky zerlegung war der falsche Begriff. Passt dir Rank-Revealing Cholesky Factorization besser?

    Zumindest erklärt es meine Verwirrung.

    0
  • I

    Super, danke.

    0
  • I

    Moin,

    ich muss hier leider nochmal nachharken, da ich mit dem Ergebnis nicht so zufrieden bin.

    Also ich habe die multiple lineare Regression (MLR), wie von ProgChild beschrieben, implementiert.
    In vielen Fällen liegt die Ebene ganz ordentlich, manchmal aber total daneben.

    Für ein Extrembeispiel hier mal zwei Screenshots.
    https://www.dropbox.com/s/hysvmj8g0ryz5hp/__mistMLR1.png
    https://www.dropbox.com/s/410na39ejuqhdqn/__mistMLR2.png
    Die dunkelblauen Punkte (rot umkreist) bilden die Punktwolke.
    Im ersten Bild in ihrer Ausgangsposition.
    Im zweiten Bild wurde die Ebene mit MLR generiert und alle Punkte auf diese projiziert
    und es ist in diesem Beispiel zu erkennen, das die Ebene "waagerecht liegt, anstatt senkrecht".
    Da ich anschließend den Flächeninhalt bestimmen möchte, ist das sehr ungünstig.

    Meine Fehleranalyse zeigt, dass ∑p⋅x~<1e−5\sum p\cdot\tilde x < 1e^{-5}. (p sind die projizierten Punkte)
    Also die Ebene liefert schon ein "guten" Ergebnis, nur liegen eben einige Punkte sehr weite drüber und eine andere sehr weit drunter, sodass diese sich zu 0 ergänzen.

    Ich habe jetzt die Ebene nicht mit MLR gefittet, sondern einfach durch 3 selbst gewählte Punkte bestimmte und es kommt das Ergebnis raus (vom Flächeninhalt) was ich für richtig halte.

    Nun ist es so, dass ich in meiner Abschlussarbeit bereist die MLR Methode ausführlich beschrieben habe und sie ungern wieder rausstreichen möchte 😃

    Hat vielleicht einer eine Idee, wie man die MLR-Methode leicht verändern könnte, um ein signifikant besseres Ergebnis zu erhalten? Mit fällt leider keine ein.

    Viele Grüße

    0
  • K

    Irgendwie bin ich mit der Problembeschreibung nicht zufrieden. "eine Ebene fitten" reicht da nicht wirklich. Was willst du jetzt wirklich? Was soll die Ebene, die du haben willst, erfüllen? Es gibt da ja mehrere Varianten, wie man die "Abwehcungen" gewichtet. Und da kommt es dann drauf an, woher die ganzen Daten kommen. Wenn z.B. x und y-Koordinate der Punkte jeweils vorgegeben ist und z mit Fehlern gemessen wird, dann muss man etwas ganz anderes machen, als wenn wenn man eine Ebene finden will, die den mittleren quadratischen Abstand zu den Punkten minimiert. Ersteres führt einem zu einem einfachen, ggf. überbestimmten Gleichungssystems und letzteres wäre über eine Eigen-Zerlegung der Kovarianzmatrix der x-/y-/z-Werte der Punkte lösbar. Das geht dann in Richtung Hauptkomponentenanalyse. Der Normalvektor ist dann der Eigenvektor zum kleinsten Eigenwert. Ein Punkt der Ebene ließe sich dann einfach über eine Mittelung der Punkte ausrechnen. Und fertig ist die Ebene.

    0
Anmelden zum Antworten