Подгонка линии с градиентным спуском

Question

Я пытаюсь подогнать линию к паре точек, используя градиентный спуск. Я не эксперт в этом и пытался записать математический алгоритм для этого в Python. Это выполняется за пару итераций, но мои прогнозы, похоже, в какой-то момент взрываются. Вот код:

import numpy as np
import matplotlib.pyplot as plt

def mean_squared_error(n, A, b, m, c):
    e = 0
    for i in range(n):
        e += (b[i] - (m*A[i] + c)) ** 2   
    return e/n

def der_wrt_m(n,A,b,m,c):
    d = 0
    for i in range(n):
        d += (2 * (b[i] - (m*A[i] + c)) * (-A[i]))
    return d/n

def der_wrt_c(n,A,b,m,c):
    d = 0
    for i in range(n):
        d += (2 * (b[i] - (m*A[i] + c)))
    return d/n

def update(n,A,b,m,c,descent_rate):
    return descent_rate * der_wrt_m(n,A,b,m,c)), descent_rate * der_wrt_c(n,A,b,m,c))

A = np.array(((0,1),
             (1,1),
             (2,1),
             (3,1)))
x = A.T[0]
b = np.array((1,2,0,3), ndmin=2 ).T
y = b.reshape(4)

def descent(x,y):
    m = 0
    c = 0

    descent_rate = 0.00001
    iterations = 100

    n = len(x)
    plt.scatter(x, y)
    u = np.linspace(0,3,100)
    prediction = 0
    for itr in range(iterations):
        print(m,c)
        prediction = prediction + m * x + c
        m,c = update(n,x,y,m,c,descent_rate)

    plt.plot(u, u * m + c, '-')   


descent(x,y)

И это мой вывод:

0 0
19.25 -10.5
-71335.1953125 24625.9453125
5593771382944640.0 -2166081169939480.2
-2.542705027685638e+48 9.692684648057364e+47
2.40856742196228e+146 -9.202614421953049e+145
-inf inf
nan nan
nan nan
nan nan
nan nan
nan nan
nan nan
etc...

---

Обновление: значения больше не взрываются, но все равно не сходятся:

# We could also solve it using gradient descent
import numpy as np
import matplotlib.pyplot as plt

def mean_squared_error(n, A, b, m, c):
    e = 0
    for i in range(n):
        e += ((b[i] - (m * A[i] + c)) ** 2)   
    #print("mse:",e/n)
    return e/n

def der_wrt_m(n,A,b,m,c):
    d = 0
    for i in range(n):
        # d += (2 * (b[i] - (m*A[i] + c)) * (-A[i]))
        d += (A[i] * (b[i] - (m*A[i] + c)))
    #print("Dm",-2 * d/n)
    return (-2 * d/n)

def der_wrt_c(n,A,b,m,c):
    d = 0
    for i in range(n):
        d += (2 * (b[i] - (m*A[i] + c)))
    #print("Dc",d/n)
    return d/n

def update(n,A,b,m,c, descent_rate):
    return (m - descent_rate * der_wrt_m(n,A,b,m,c)),(c - descent_rate * der_wrt_c(n,A,b,m,c))

A = np.array(((0,1),
             (1,1),
             (2,1),
             (3,1)))
x = A.T[0]
b = np.array((1,2,0,3), ndmin=2 ).T
y = b.reshape(4)

def descent(x,y):
    m = 0
    c = 0

    descent_rate = 0.0001
    iterations = 10000

    n = len(x)
    plt.scatter(x, y)
    u = np.linspace(0,3,100)
    prediction = 0
    for itr in range(iterations):
        prediction = prediction + m * x + c
        m,c = update(n,x,y,m,c,descent_rate)
        loss = mean_squared_error(n, A, b, m, c)

    print(loss)
    print(m,c)
    plt.plot(u, u * m + c, '-')    

descent(x,y)

И теперь график выглядит примерно так после примерно 10000 итераций со скоростью обучения 0,0001:

[4.10833186 5.21468937]
1.503547594304175 -1.9947003678083184

В то время как наименьший квадрат соответствует примерно так:

Answer 1

В вашей функции обновления вы должны вычесть вычисленные градиенты из текущих m и c

def update(n,A,b,m,c,descent_rate):
    return m - (descent_rate * der_wrt_m(n,A,b,m,c)), c - (descent_rate * der_wrt_c(n,A,b,m,c))

Обновление: вот рабочая версия.Я получил матрицу A после получения x, y, так как это меня смущает =).Например, в ваших вычислениях градиента у вас есть выражение d += (A[i] * (b[i] - (m*A[i] + c))), но оно должно быть d += (x[i] * (b[i] - (m*x[i] + c))), поскольку x [i] дает вам один элемент, тогда как A [i] дает вам список.

Также вы забылизнак минус при расчете производной по с.Если ваше выражение (y - (m*x + c))^2), то производная по c должна быть 2 * (-1) * (y - (m*x + c)), поскольку перед c стоит минус.

# We could also solve it using gradient descent
import numpy as np
import matplotlib.pyplot as plt

def mean_squared_error(n, x, y, m, c):
    e = 0
    for i in range(n):
        e += (m*x[i]+c - y[i])**2
    e = e/n
    return e/n

def der_wrt_m(n, x, y, m, c):
    d = 0
    for i in range(n):
        d += x[i] * (y[i] - (m*x[i] + c))
    d = -2 * d/n
    return d

def der_wrt_c(n, x, y, m, c):
    d = 0
    for i in range(n):
        d += (y[i] - (m*x[i] + c))
    d = -2 * d/n
    return d


def update(n,x,y,m,c, descent_rate):
    return (m - descent_rate * der_wrt_m(n,x,y,m,c)),(c - descent_rate * der_wrt_c(n,x,y,m,c))


A = np.array(((0,1),
             (1,1),
             (2,1),
             (3,1)))
x = A.T[0]
b = np.array((1,2,0,3), ndmin=2 ).T
y = b.reshape(4)

print(x)
print(y)

def descent(x,y):
    m = 0.0
    c = 0.0

    descent_rate = 0.01
    iterations = 10000

    n = len(x)
    plt.scatter(x, y)
    u = np.linspace(0,3,100)
    prediction = 0
    for itr in range(iterations):
        prediction = prediction + m * x + c
        m,c = update(n,x,y,m,c,descent_rate)
        loss = mean_squared_error(n, x, y, m, c)
        print(loss)

    print(loss)
    print(m,c)
    plt.plot(u, u * m + c, '-')    
    plt.show()

descent(x,y)