1

I wrote some code that performs gradient descent on a couple of data points. For some reason the curve is not converging correctly, but I have no idea why that is. I always end up with an exploding tail.

Am I doing one of the computations wrong? Am I actually getting stuck in a local minimum or is it something else?

Here is my code:

import numpy as np
import matplotlib.pyplot as plt

def estimate(weights, x, order):
    est = 0
    for i in range(order):
        est += weights[i] * x ** i
    return est

def cost_function(x, y, weights, m):
    cost = 0
    for i in range(m-1):
        cost += (((weights[i] * x ** i) - y) ** 2)
    return (np.sum(cost ** 2) / ( 2 * m ))

def descent(A, b, iterations, descent_rate, order):
    x = A.T[0]
    y = b.reshape(4)

    # features
    ones = np.vstack(np.ones(len(A)))
    x = np.vstack(A.T[0])
    x2 = np.vstack(A.T[0] ** 2)

    # Our feature matrix
    features = np.concatenate((ones,x,x2), axis = 1).T
    # Initialize our coefficients to zero
    weights = np.zeros(order + 1)
    m = len(y)

    # gradient descent
    for i in range(iterations):
        est = estimate(weights, x, order).T
        difference = est - y
        weights = weights + (-descent_rate * (1/m) * np.matmul(difference, features.T)[0])
        cost = cost_function(x, y, weights, m)
        print(cost)

    plt.scatter(x,y)
    u = np.linspace(0,3,100)
    plt.plot(u, (u ** 2) * weights[2] + u  * weights[1] + weights[0], '-')
    plt.show()

A = np.array(((0,1),
             (1,1),
             (2,1),
             (3,1)))

b = np.array((1,2,0,3), ndmin = 2 ).T

iterations = 150
descent_rate = 0.01
order = 2
descent(A, b, iterations, descent_rate, order)

I would like to avoid getting stuck in such a minimum. I have attempted setting the initial weights to random values but to no avail, sometimes it dips a bit more but then gives me the same behaviour again.

Here is the one of the plots that I am getting: enter image description here

And here is the expected result obtained by a least squares solution:

enter image description here

1 Answer 1

1

Your estimate function should be

def estimate(weights, x, order):
    est = 0
    for i in range(order+1):
        est += weights[i] * x ** i
    return est

Better yet, since the order information is already present in the size of the weights vector, remove the redundancy with:

def estimate(weights, x):
    est = 0
    for i in range(len(weights)):
        est += weights[i] * x ** i
    return est

This is what I got when using your code and running 2000 iterations: enter image description here

2
  • Thank you so much! You probably saved me hours of frustration tomorrow morning :D Final weights [ 1.4 -1.1 0.5] it gives me the same as when using numpy's polyfit function. Also, out of curiosity, how were you able to find the error so quickly?/ May 19, 2019 at 17:08
  • 1
    Glad to know! There's no trick. I've probably made most mistakes I could possibly make in code like that. With time spotting issues becomes automatic.
    – foglerit
    May 19, 2019 at 17:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.