Gradient Descent on a logarithmic decline curve in Python

Question

I wish to run gradient descent on a logarithmic decline curve as represented by:

y = y0 - a * ln(b + x).

My y0 for this example: 800

I was trying to do this using partial derivatives with respect to a and b, but while this apparently minimizes the squared error, it doesn't converge. I know this isn't vectorized, and I might be taking the wrong approach altogether. Am I making any simple mistake, or off on this problem altogether?

import numpy as np

# constants my gradient descent model should find:
a = 4
b = 4

# function to fit on!
def function(x, a, b):
    y0 = 800
    return y0 - a * np.log(b + x)

# Generates data
def gen_data(numpoints):
    a = 4
    b = 4
    x = np.array(range(0, numpoints))
    y = function(x, a, b)
    return x, y
x, y = gen_data(600)

def grad_model(x, y, iterations):
    converged = False

    # length of dataset
    m = len(x)

    # guess   a ,  b
    theta = [0.1, 0.1]
    alpha = 0.001

    # initial error
    e = np.sum((np.square(function(x, theta[0], theta[1])) - y))

    for iteration in range(iterations):
        hypothesis = function(x, theta[0], theta[1])
        loss = hypothesis - y

        # compute partial deritaves to find slope to "fall" into
        theta0_grad = (np.mean(np.sum(-np.log(x + y)))) / (m)
        theta1_grad = (np.mean((((np.log(theta[1] + x)) / theta[0]) - (x*(np.log(theta[1] + x)) / theta[0])))) / (2*m)

        theta0 = theta[0] - (alpha * theta0_grad)
        theta1 = theta[1] - (alpha * theta1_grad)

        theta[1] = theta1
        theta[0] = theta0

        new_e = np.sum(np.square((function(x, theta[0], theta[1])) - y))
        if new_e > e:
            print "AHHHH!"
            print "Iteration: "+ str(iteration)
            break
        print theta
    return theta[0], theta[1]

Answer 1

I found a few bugs in your code. The line

e = np.sum((np.square(function(x, theta[0], theta[1])) - y))

is incorrect and should be replaced with

e = np.sum((np.square(function(x, theta[0], theta[1]) - y)))

The formula for new_e contains the same error.

Also, the gradient formulas are wrong. Your loss function is $L(a,b) = \sum_{i=1}^N y_0 - a \log(b + x_i)$, so you have to compute the partial derivatives of $L$ with respect to $a$ and $b$. (Does LaTeX really not work on stackoverflow?) A final point is that the gradient descent method has a step size restriction, so our step size must not be too large. Here is a version of your code that is working better:

import numpy as np
import matplotlib.pyplot as plt

# constants my gradient descent model should find:
a = 4.0
b = 4.0
y0 = 800.0

# function to fit on!
def function(x, a, b):
    # y0 = 800
    return y0 - a * np.log(b + x)

# Generates data
def gen_data(numpoints):
    # a = 4
    # b = 4
    x = np.array(range(0, numpoints))
    y = function(x, a, b)
    return x, y
x, y = gen_data(600)

def grad_model(x, y, iterations):
    converged = False

    # length of dataset
    m = len(x)

    # guess   a ,  b
    theta = [0.1, 0.1]
    alpha = 0.00001

    # initial error
    # e = np.sum((np.square(function(x, theta[0], theta[1])) - y))    #  This was a bug
    e = np.sum((np.square(function(x, theta[0], theta[1]) - y)))

    costs = np.zeros(iterations)

    for iteration in range(iterations):
        hypothesis = function(x, theta[0], theta[1])
        loss = hypothesis - y

        # compute partial deritaves to find slope to "fall" into
        # theta0_grad = (np.mean(np.sum(-np.log(x + y)))) / (m)
        # theta1_grad = (np.mean((((np.log(theta[1] + x)) / theta[0]) - (x*(np.log(theta[1] + x)) / theta[0])))) / (2*m)
        theta0_grad = 2*np.sum((y0 - theta[0]*np.log(theta[1] + x) - y)*(-np.log(theta[1] + x)))
        theta1_grad = 2*np.sum((y0 - theta[0]*np.log(theta[1] + x) - y)*(-theta[0]/(b + x)))

        theta0 = theta[0] - (alpha * theta0_grad)
        theta1 = theta[1] - (alpha * theta1_grad)

        theta[1] = theta1
        theta[0] = theta0

        # new_e = np.sum(np.square((function(x, theta[0], theta[1])) - y)) # This was a bug
        new_e = np.sum(np.square((function(x, theta[0], theta[1]) - y)))
        costs[iteration] = new_e
        if new_e > e:
            print "AHHHH!"
            print "Iteration: "+ str(iteration)
            # break
        print theta
    return theta[0], theta[1], costs

(theta0,theta1,costs) = grad_model(x,y,100000)
plt.semilogy(costs)