What's wrong with my gradient descent? (simple, one variable)

Question

# calculating the cost function
theta = np.array([[0], [0]])
def h(X, theta):
    return np.matmul(X, theta)

def cost(X, y, theta, m):
    return 1/(2* m) * np.sum((h(X, theta)-y)**2)

print(cost(X, y, theta, m))
# 1-variable gradient decent
cost_vals = []
def grad_decent(X, y, theta, m, alpha):
    theta1 = theta[0, 0]
    theta2 = theta[1, 0]
    for i in range(0, iter):
        cost_vals.append(cost(X, y, np.array([[theta1], [theta2]]), m))
        theta1 = theta1 - alpha/m * np.sum((h(X, theta)-y) * X[:, 0])
        theta2 = theta2 - alpha/m * np.sum((h(X, theta)-y) * X[:, 1])
        # print(f"Iteration {i}: ", theta1, ", ", theta2)
        
    return np.array([[theta1], [theta2]])

As one can see I'm doing gradient descent with only one feature, i.e. I want to fit a line of the form theta1 + theta2 * x to my data. I know that for such small feature sizes one would definitely use the normal equation, but I wanted to implement it with gradient descent. My learning rate alpha is 0.01 and I'm following the machine learning course of andrew ng on coursera. However, instead of converging, my cost values increase so there is some mistake in the algorithm which I can't find. X is of the shape 97x2 where the first column is just ones since there is only one feature but I want to fit a polynomial with two coefficients.

Answer 1

You set your model parameters using np.array([[0], [0]]), meaning that your parameter array is of shape (2, 1). Changing this to np.array([0, 0]) (with shape (2,)) decreases the loss over time as expected using a learning rate $\alpha$ of 0.01. Looking more closely at the difference of using the two, I see that your way returns way higher losses (exactly 100x). Using a smaller $\alpha$, e.g. 0.0001, with your implementation also works:

Therefore the main reason of your implementation is your choice of $\alpha$, which is too high, preventing the model to reach the optimum.

EDIT: The issues is twofold, the first is indeed the shape or your theta array. The second has to do with the fact that you are using the theta array in your for loop which always contains the original theta parameters, as you are not storing and using the updated theta1 and theta2. The following code seems to give a good result of $\theta_1 \approx -3.2$ and $\theta_2 \approx 1.12$ using 1000 iterations with $\alpha=0.01$

import numpy as np

# load and transform data
data = np.genfromtxt("data.txt", delimiter=",")
X = np.stack((np.ones(data.shape[0]), data[:,0]), axis=-1)
y = data[:,1]

m = X.shape[0]
iter = 1000
 
# calculating the cost function
theta = np.array([[0.3], [0.5]])
def h(X, theta):
    return np.matmul(X, theta)
 
def cost(X, y, theta, m):
    return 1/(2* m) * np.sum((h(X, theta)-y)**2)
 
print(cost(X, y, theta, m))
# 1-variable gradient decent
cost_vals = []
def grad_decent(X, y, theta, m, alpha):
    theta1 = theta[0, 0]
    theta2 = theta[1, 0]
    for i in range(0, iter):
        theta = np.array([theta1, theta2])
        cost_vals.append(cost(X, y, np.array([theta1, theta2]), m))
        theta1 = theta1 - alpha/m * np.sum((h(X, theta)-y) * X[:, 0])
        theta2 = theta2 - alpha/m * np.sum((h(X, theta)-y) * X[:, 1])
        # print(f"Iteration {i}: ", theta1, ", ", theta2)
        
    return np.array([[theta1], [theta2]])

grad_decent(X, y, theta, m, 0.01)
# 1740.4098044406999
# array([[-3.20070025],
#        [ 1.12320526]])