Метод подгонки + условия Вольфа - реализация в Python3 - PullRequest
Новая
       66

Метод подгонки + условия Вольфа - реализация в Python3

0 голосов
/ 31 января 2020

Я хочу запрограммировать процедуру подгонки для алгоритма Гаусса-Ньютона + Условие Вольфа в Python 3. Уже существует аналогичный пост за 5 лет go стек . Код там написан на Python 2 и не включает условия Вольфа.

Я написал подпрограмму, которая отлично работает для тестовых данных из википедии Алгоритм Гаусса-Ньютона . 

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      
        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))

Но если попробовать эту процедуру на другом наборе данных + модель, например 

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])

с моделью f (x) = c1 * x² + c2 

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

это не работает вообще (аналогично для других примеров). Условия Вольфа никогда не заполняются полностью, и параметры не изменяются (идеальные параметры для примера: c_0 = 0,751, c_1 = 7,82).

Кто-то видит, где я допустил ошибку?

...