Polynomfit 3.Grades



  • Hallo Community,

    ich verzweifle langsam an diesem Problem.
    Kann mir jemand Tips oder Inspirationen geben, wie ich einen Polynomfit dritten Grades hinbekomme. Habe versucht mich mit GnuSL durchzuschlagen, aber funktioniert bisher leider kein Stück und ist vielleicht auch ein bisschen Overkill für so einen relativ "einfache" Fit.

    Ich habe einige x und y Werte und möchte sie least-squared an ein Polynom der Form f(x) = a*x³ + b*x² + c*x + d fitten.

    Kann mir jemand helfen?

    Danke und beste Grüße



  • Ableitung der Fehlerfunktion berechnen und dann Gradientenabstieg anwenden.



  • otze schrieb:

    Ableitung der Fehlerfunktion berechnen und dann Gradientenabstieg anwenden.

    Ein solcher Polynomfit ist ein lineares least squares Problem. Das kann man auch direkt lösen.

    david_jon schrieb:

    Kann mir jemand helfen?

    Stelle ein lineares Gleichunssystem auf, also eine Gleichung pro x/y Paar. Die die Spaltenzahl der dazugehörigen Matrix wäre dann 4, weil du 4 unbekannte hast. Wenn du dies im Sinne der kleinsten Fehlerquadrate lösen willst, dann nehme ich an, dass du mehr als 4 Punktpaare hast; denn dann erhältst du ein überbestimmtes Gleichungssystem. Das System kannst du per Householder-Transformationen zu einem mit Dreiecksgestalt überführen, welches direkt lösbar ist. Oder du verwendest eine Bibliothek für lineare Algebra, die dir auch die QR-Zerlegung berechnen kann und machst das dann damit.



  • krümelkacker schrieb:

    otze schrieb:

    Ableitung der Fehlerfunktion berechnen und dann Gradientenabstieg anwenden.

    Ein solcher Polynomfit ist ein lineares least squares Problem. Das kann man auch direkt lösen.

    In der Tat. du hast recht. Dann braucht man aber auch keine Householder Trafo. Das Problem ist ja 3-Dimensional. Ein Newton-Schritt auf dem MSE reicht da.



  • otze schrieb:

    krümelkacker schrieb:

    otze schrieb:

    Ableitung der Fehlerfunktion berechnen und dann Gradientenabstieg anwenden.

    Ein solcher Polynomfit ist ein lineares least squares Problem. Das kann man auch direkt lösen.

    In der Tat. du hast recht. Dann braucht man aber auch keine Householder Trafo. Das Problem ist ja 3-Dimensional. Ein Newton-Schritt auf dem MSE reicht da.

    Es gibt 4 Unbekannte: a,b,c,d. Was du mit "Newton-Schritt auf dem MSE" meinst, ist mir nicht klar.

    Siehe
    http://de.wikipedia.org/wiki/Vandermonde-Matrix#Anwendung:Polynominterpolation
    http://en.wikipedia.org/wiki/Linear_least_squares
    (mathematics)#Computation

    Der erste Link verrät, wie das zu lösende Gleichungssystem aussieht, der zweite, wie man dies im Sinne der kleinsten Fehlerquadrate lösen kann.

    Im Wesentlichen gibt es folgende zwei Möglichkeiten:
    (1) Normalgleichungssystem aufstellen -> Cholesky-Faktorisierung -> Rücksubstitution
    (2) Orthogonale Transformation durchführen -> Rücksubstitution

    Das erste mag für ein bestimmtes Problem schneller sein, ist aber Rundungsfehleranfälliger. Wenn es einem auf Genauigkeit ankommt, dann ist der zweite Weg von Vorteil.

    Orthogonale Transformationen, Cholesky und Rücksubstitution muss man nicht neuerfinden. Das gibt's in jeder gescheiten Bibliothek für lineare Algebra. Zum Beispiel auch Eigen.



  • david_jon schrieb:

    Ich habe einige x und y Werte und möchte sie least-squared an ein Polynom der Form f(x) = a*x³ + b*x² + c*x + d fitten.

    Such mal in der Bücherei deines Vertrauens nach "Numerical Recipes in C++" - die dürften auch für solche Probleme die passende Formel bereit liegen haben 😉



  • krümelkacker schrieb:

    Es gibt 4 Unbekannte: a,b,c,d. Was du mit "Newton-Schritt auf dem MSE" meinst, ist mir nicht klar.

    Du weißt was der MSE ist. Du nennst ihn nur Least Squares. Das ist in diesem Fall, wie du selbst sagst ein quadratisches Optimierungsproblem ohne Randbedingungen.

    f(x,a,b,c,d)=ax+bx2+cx3+d
    MSE(a,b,c,d)=sum_i (y_i-f(x_i,a,b,c,d))^2
    => nach a,b,c,d 2 mal ableiten, dann newtonschritt.



  • otze schrieb:

    Du weißt was der MSE ist. Du nennst ihn nur Least Squares.

    Ich nehme an, dass du mit MSE "mean square error" meinst.

    otze schrieb:

    Das ist in diesem Fall, wie du selbst sagst ein quadratisches Optimierungsproblem ohne Randbedingungen.

    Ich sagte nie, es sei ein quadratisches Problem. Ich sagte "es ist ein lineares kleinste-Quadrate-Problem".

    otze schrieb:

    f(x,a,b,c,d)=ax+bx2+cx3+d
    MSE(a,b,c,d)=sum_i (y_i-f(x_i,a,b,c,d))^2
    => nach a,b,c,d 2 mal ableiten, dann newtonschritt.

    Was du meinst, ist glaub'ich ein Newton-Schritt auf dem Gradient von f. Da wo f ein Minimum hat, ist der Gradient von f ja 0:

    uk+1=uk(f(u))1f(u)u_{k+1} = u_{k} - \left( f''(u) \right)^{-1} * f'(u)

    Das führt zu dem bereits erwähnten Normalgleichungssystem. Das f''(u) ist hier eine von u unabhängige quadratische Matrix. Und dieses Normalgleichungssystem ist linear in u.

    kk



  • der Trick ist, dass f''(u) invertierbar ist und deswegen nicht solche Hämmer wie QR-Zerlegung oder andere Psudoinversen braucht. Es läuft am Ende aufs Selbe hinaus, aber wenn jemand kommt mit "ja, ähh ich krieg das nicht implementiert" dann wird er sicherlich auch keine QR bei der Hand haben, schon gar keine mit Pivoting.



  • Eine QR-Zerlegung ist keinen Deut komplizierter oder komplexer als eine Cholesky-Zerlegung, mit der man das Normalgleichungssystem lösen würde. Und wie gesagt: Das Normalgleichungssystem hat eine typischerweise schlechtere Kondition. Mach deinen Kram wie du willst. Aber behaupte nicht, dass das Lösen des Normalgleichungssystems für das Problem des OPs auf irgend eine Art und Weise "besser" sei.

    Das mit der orthogonalen Transformation ist GANZ EINFACH. Stellt man das Gleichungssystem mit 6 Punkten auf, wird es so aussehen:

    ? ? ? ?       ?
    ? ? ? ?   d   ?
    ? ? ? ? * c = ?
    ? ? ? ?   b   ?
    ? ? ? ?   a   ?
    ? ? ? ?       ?
    

    wobei die Fragezeichen für unterschiedliche Werte stehen. Nach den Householder-Transformationen, die sich IN-PLACE berechnen lassen, sieht das System dann so aus:

    ? ? ? ?       ?
    0 ? ? ?   d   ?
    0 0 ? ?   c   ?
    0 0 0 ? * b = ?
    -------   a   -
    0 0 0 0       ?
    0 0 0 0       ?
    

    Den oberen Teil kannst du direkt per Rücksubstitution lösen. Der untere Teil der rechten Seite verrät einem direkt den Fehler. Einfach quadrieren und aufaddieren und es ist die Fehlerquadratsumme.

    Das ist standard Pipikram ... lernt man gleich im Mathegrundstudium, wie man überbestimmte Gleichunssysteme im Sinne der kleinsten Fehlerquadrate löst. Und da bekommt man dann auch eingebleut, dass man das mit dem Normalgleichungssystem vermeiden soll, weil numerisch ungenau.

    aber wenn jemand kommt mit "ja, ähh ich krieg das nicht implementiert" dann wird er sicherlich auch keine QR bei der Hand haben, schon gar keine mit Pivoting.

    Für keinen Ansatz muss man Code neu erfinden. Cholesky findest du genauso wie QR in gescheiten Bibliotheken. Wenn man es wirklich selbst bauen will, macht es vom Schwierigkeitsgrad keinen Unterschied, ob Du eine Cholesky-Zerlegung oder diese orthogonale Transformation selbst programmierst. Ich würde sogar behaupten, dass der Ansatz mit dem Normalgleichungssystem komplizierter ist.



  • vielen Dank an alle die sich hier beteiligt haben!
    Hab die Sache verstanden und umgesetzt.

    Da es so aussieht als wären hier wirklich Experten am Werk, noch eine klitze kleine Frage 😉

    Habt ihr eine Idee wie man einen (wahrscheinlich nennt man es so) zweidimensionalen Fit zustande bekommt?

    sprich ich habe (x_soll, y_soll) = Funktion von (x_gemessen, y_gemessen) bis in die dritte Ordnung und muss einen Fit erstellen. Das wird vermutlich etwas komplizierter ?
    Braucht ihr noch mehr Informationen um mir Tips geben zu können ?

    Danke und Grüße



  • david_jon schrieb:

    Hab die Sache verstanden und umgesetzt.

    Dann solltest Du dein aktuelles Problem auch selbst lösen können, wenn das stimmt. 😉

    Also, wenn f(x,y) = ( f_1(x,y), f_2(x,y) ) linear in den Unbekannten ist, dann lässt sich das Problem genauso lösen. Einmal für f1 und einmal für f2.

    Wobei es da schon einen größeren Haufen an Unbekannten gibt. Genauergesagt: Insgesamt 20, wenn f_i so aussieht:

    fi(x,y)=(1,x,y,x2,xy,y2,x3,x2y,xy2,y3)(u1,iu2,iu10,i)f_i(x,y) = \left(1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3\right) \cdot \left( \begin{matrix} u_{1,i} \\ u_{2,i} \\ \vdots \\ u_{10,i} \\ \end{matrix} \right)

    Also, für i=1 hast du 10 Unbekannte und für i=2 nochmal 10 Unbekannte. Das wären dann zwei Gleichungssysteme mit jeweils 10 Unbekannten, die zu lösen wären.


Log in to reply