Fitting-Algorithmus mit Wolfe Bedingungen

[Particledust]

Guten Morgen zusammen,

ich habe einen Algorithmus in Python 3 geschrieben, in welchem die Gauß-Newton-Methode und die Wolfe Bedinungen implementiert sind. Es handelt sich hierbei um eine Fitting-Methode um die optimalen Parameter eines mathematischen Models zu finden.

Mein Code funktioniert für das angegebene Beispiel in Gauß-Newton-Methode

Code: Alles auswählen

from numpy import array, zeros, dot, sqrt
from numpy.linalg import inv

# Data from Wiki
x  = array([0.038, 0.194, 0.425, 0.626,  1.253,  2.5,    3.74])
y  = array([0.05,  0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317])
c  = array([1000, 40]).astype(float)  # initial parameters

accuracy = 10**(-8)
maxiter  = 10**3 

def model(x,c):
    return (c[0] * x) / (c[1] + x)

def residual(c):
    return y - model(x,c)

def jacobi(c):
    jac = zeros((len(x), len(c0)))
    jac[:,0] = -x/(c[1] + x)                # partial derivatives 
    jac[:,1] = (c[0] * x)/(c[1] + x)**2
    return jac

def squares(c):
    return sum(residual(c)**2)

def grad(c):
    val = jacobi(c).T.dot(residual(c))
    return sqrt(val[0]**2 + val[1]**2)

# main loop
counter = 0
for i in range(maxiter):
    r = residual(c)
    Jf  = jacobi(c)    # Jacobi and its Transpose
    Jft = Jf.T

    for j in range(1,100):
        c_old = c
        c_new = c - 1/j * dot(dot(inv(dot(Jft,Jf)),Jft),r)
        # Compute values to apply the wolfe conditions, 0.9 and 0.0001 from wiki      
        old = squares(c_old) + 0.0001*1/j*grad(c_old).dot(dot(dot(inv(dot(Jft,Jf)),Jft),r))                                                                                                   
        new = squares(c_new)

        ka1 = -grad(c_new).dot(dot(dot(inv(dot(Jft,Jf)),Jft),r))
        ka2 = -0.9 * grad(c_old).dot(dot(dot(inv(dot(Jft,Jf)),Jft),r))

        if new <= old and ka1 <= ka2:    # Armijo rule and curvature condition
            c = c_new
            break

    change   = grad_abs(c)
    counter +=1

    if change < accuracy:
        break

print('Iterations = ' +str(counter))
print('Parameter = ' +str(c))

Wenn ich diesen Algorithmus für ein anderes simples Beispiel verwende:

Code: Alles auswählen

x = array([1,2,3,4,5,6,7,8]).astype(float)
y = array([8.3, 11.0, 14.7, 19.7, 26.7, 35.2, 44.4, 55.9]).astype(float)
c = array([10, 7])   #optimal paramerets are: c_0 = 0.751, c_1=7.82

Mit dem Model f(x) = c1 * x² + c2

Code: Alles auswählen

def model(x,c):
    return c[0] * x**2 + c[1]

def jacobi(c):
    jac = zeros((len(x), len(c)))
    jac[:,0] = x**2 
    jac[:,1] = 1
    return jac

Dann funktioniert der Algorithmus überhaupt nicht mehr. Die Wolfe Bedingungen sind niemals erfüllt und die Parameter verändern sich nicht.
Die Wolfe Bedingungen sind höchstwahrscheinlich falsch, aber ich finde den Fehler nicht. Hat jemand von euch eine Idee, wo der Fehler liegt?

[Particledust]

Hat niemand eine Idee?

[__blackjack__]

@Particledust: Also der erste Code funktioniert so sicher nicht, denn `c0` und `grad_abs()` sind nicht definiert.

Auf Modulebene sollte nur Code stehen der Konstanten, Funktionen, und Klassen definiert. Das Hauptprogramm steht üblicherweise in einer Funktion die `main()` heisst. Wenn man das in eine Funktion steckt, fällt auf das einige Werte in anderen Funktionen einfach so auf ”magische” Weise aus der Umgebung verwendet werden. Alles was Funktionen und Methoden ausser Konstanten benötigen sollte als Argument(e) übergeben werden.

`i` wird nirgends verwendet, dafür wird in der Schleife manuell ein `counter` mitgezählt das die gleichen Werte annimmt wie `i`.

Der Teilausdruck ``dot(dot(inv(dot(Jft, Jf)), Jft), r)`` kommt in der Schleife insgesamt vier mal vor. Kann man dem vielleicht einen sinnvollen Namen geben?

Apropos sinnvolle Namen: Kryptische ein bis dreibuchstabige Namen die sich teilweise nicht an die üblichen Gross-/Kleinschreibungskonventionen halten, sind in der Regel nicht sinnvoll.

Hilft der Name `c_old` beim Verständnis? Falls nein, weg damit.

Zeichenketten und Werte mit `str()` und ``+`` zusammenstückeln ist kein Python sondern eher BASIC. Python hat dafür Zeichenkettenformatierung mit der `format()`-Methode und ab Python 3.6 f-Zeichenkettenliterale.

Zwischenstand:

Code: Alles auswählen

#!/usr/bin/env python3
from numpy import array, dot, sqrt, zeros
from numpy.linalg import inv


ACCURACY = 10 ** -8
MAX_ITERATION_COUNT = 10 ** 3


def model(x, c):
    return (c[0] * x) / (c[1] + x)


def residual(x, y, c):
    return y - model(x, c)


def jacobi(x, c):
    jac = zeros((len(x), len(c0)))
    jac[:, 0] = -x / (c[1] + x)  # partial derivatives
    jac[:, 1] = (c[0] * x) / (c[1] + x) ** 2
    return jac


def squares(x, y, c):
    return sum(residual(x, y, c) ** 2)


def grad(x, y, c):
    val = jacobi(x, c).T.dot(residual(x, y, c))
    return sqrt(val[0] ** 2 + val[1] ** 2)


def main():
    # Data from Wiki.
    x = array([0.038, 0.194, 0.425, 0.626, 1.253, 2.5, 3.74])
    y = array([0.05, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317])
    c = array([1000, 40]).astype(float)  # initial parameters

    for counter in range(1, MAX_ITERATION_COUNT + 1):
        r = residual(x, y, c)
        Jf = jacobi(x, c)  # Jacobi and its transposition.
        Jft = Jf.T

        for j in range(1, 100):
            maybe_there_is_some_good_name_for_it = dot(
                dot(inv(dot(Jft, Jf)), Jft), r
            )
            c_old = c
            c_new = c - 1 / j * maybe_there_is_some_good_name_for_it
            #
            # Compute values to apply the wolfe conditions, 0.9 and 0.0001 from
            # wiki.
            #
            old = squares(x, y, c_old) + 0.0001 * 1 / j * grad(
                x, y, c_old
            ).dot(maybe_there_is_some_good_name_for_it)
            new = squares(x, y, c_new)

            ka1 = -grad(x, y, c_new).dot(maybe_there_is_some_good_name_for_it)
            ka2 = -0.9 * grad(x, y, c_old).dot(
                maybe_there_is_some_good_name_for_it
            )
            #
            # Armijo rule and curvature condition.
            #
            if new <= old and ka1 <= ka2:
                c = c_new
                break

        if grad_abs(c) < ACCURACY:
            break

    print(f"Iterations = {counter}")
    print(f"Parameter = {c}")


if __name__ == "__main__":
    main()

[Particledust]

Vielen Dank für deine Verbesserungsvorschläge. Ich werde diese in Zukunft beachten.

Beim Überarbeiten des Codes hab ich den Fehler gefunden. Der Algorithmus funktioniert jetzt!