Lineare Regression im 3D



  • 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



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



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



  • 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?



  • 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.



  • 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?



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



  • 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?



  • 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

    bAx~\| 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.



  • Super, danke.



  • 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 px~<1e5\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



  • 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.



  • otze schrieb:

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

    Je nach Größe und Sparseness wird das iterativ oder direkt gelöst. Einem iterativen Löser wie CG ist das egal, ob die Matrix vollen Rang hatte. Er gibt mir eine Minimum-Norm-Lösung. Und sonst kann man das System noch regularisieren.



  • krümelkacker schrieb:

    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.

    x,y und z sind fest gegeben ohne Fehler, sprich ich habe ein überbestimmtes Gleichungssystem.

    Also ich habe zwei Objekte, die kollidieren und will die gemeinsame Kollisionsfläche bestimmen.
    Die Punkte, die diese Fläche bilden sollen, habe ich bestimmt (dunkelblaue Punkte in den Beispielbildern) und jetzt soll durch diese eine gemeinsame Ebene gelegt werden, um anschließend die Punkte auf diese zu projizieren.

    So klarer?



  • [quote="krümelkacker"]

    otze schrieb:

    Je nach Größe und Sparseness wird das iterativ oder direkt gelöst.

    Hier gings gerade um direkt gelöst. Vielleicht sollte ich in Zukunft auch immer SeppJs disclaimer dazu schreiben, damit ich, wenn auch einen monat später, interessante Antworten kriege.

    //edit ausserdem fällt CG nun wirklich nicht in die kategorie numerisch stabil. Das war ja der Ausgangspunkt der Frage.



  • icewater schrieb:

    x,y und z sind fest gegeben ohne Fehler, sprich ich habe ein überbestimmtes Gleichungssystem.

    Bei überbestimmten Gleichungen gibt es nicht unbedingt eine Lösung.

    Beispiel:
    x=2
    x=3

    Ich möchte von dir also genau wissen, was für eine Ebene du eigentlich suchst.

    icewater schrieb:

    Also ich habe zwei Objekte, die kollidieren und will die gemeinsame Kollisionsfläche bestimmen.
    Die Punkte, die diese Fläche bilden sollen, habe ich bestimmt (dunkelblaue Punkte in den Beispielbildern) und jetzt soll durch diese eine gemeinsame Ebene gelegt werden, um anschließend die Punkte auf diese zu projizieren.
    So klarer?

    Nein. Du hast nicht verraten, was das für eine Ebene sein soll. Ich weiß zwar jetzt, dass Dir klar ist, dass es bei 4 oder mehr Punkten keine Ebene geben muss, die all diese Punkte enthält, weil du in dem Zusammenhang etwas von Punkte projizieren sagtest, aber ich weiß immer noch nicht, was das jetzt für eine Ebene sein soll. Ohne Problembeschreibung keine Lösung! Was für eine Kostenfunktion soll die Ebene minimieren? Möchtest du vielleicht den mittleren quadratischen Abstand zwischen den Punkten und der Ebene minimieren? Wenn ja, dann mach dich mal bzgl. Hauptkomponentenanalyse schlau.


Anmelden zum Antworten