# Linear Regression in Python using gradient descent

(cross-posted from data-science StackExchange)(someone recommended that this community is more appropriate for my problem)

I am trying to implement a simple multivariate linear regression model without using any inbuilt machine libraries. So far, I have been able to get a root mean squared error for training about $$2.93$$ and the model from the normal (closed-form) equation is able to produce a training RMSE of $$~2.3$$. I am looking for ways in which I can improve my implementation of the gradient descent algorithm. Below is my implementation:

My gradient descent method looks like this: $$\theta = \theta - [(\alpha/N) * X (X\theta - Y)]$$ where $$\theta$$ is the model parameter, $$N$$ is the number of training elements, $$X$$ is the input and $$Y$$ are the target elements. $$\alpha$$ is the step size.

def gradientDescent(self):
for i in range(self.iters):
# T = T - (\alpha/N) * X*(XT - Y)
self.theta = self.theta - (self.alpha/len(self.X)) * np.sum(self.X * (self.X @ self.theta.T - self.Y), axis=0)
return errors


I had set the $$\alpha$$ as $$0.1$$ and number of iterations as 1000. The gradient descent reaches convergence at around 700-800 iterations (checked).

My error function is like:

def error_function(self):
# Error function: (1/2N) * (XT - Y)^2 where T is theta
error_values = np.power(((self.X @ self.theta.T) - self.Y), 2)
return np.sum(error_values)/(2 * len(self.X))


I was expecting the training error from the gradient descent and the normal equations would turn out to be similar, but they have a bit of a huge difference. So, I wanted to know whether I am doing anything wrong or not.

PS I have not normalized the data, yet. Normalizing leads to a much lower RMSE (~$$0.22$$)

• Gradient descent is pretty much the worst way to estimate linear regression parameters. To actually get it to work, you'll have to spend a lot of computational time and perhaps implement tricks like learning rate scheduling and momentum. Gradient descent is only useful to know about because of its utility in neural network optimization; in nearly all other contexts, better methods are available. stats.stackexchange.com/questions/160179/… Unless your sole purpose is to learn about g.d., use another method. – Sycorax Sep 18 '19 at 14:19
• Also, in math you wrote $\frac{\alpha}{2N}$ but in code you wrote $\frac{\alpha}{N}$ for gradientDescent. The difference just means you have a different learning rate between the two expressions, but I think you would prefer to be consistent. – Sycorax Sep 18 '19 at 14:21
• Well, I am just getting started into the field and was assigned to do the same, just to test my understanding. Although, thanks for the alternative methods. I will surely look up to them, sooner. Also, $\alpha/2N$ looks like a newbish mistake. I will patch that up. – MaJoR21 Sep 18 '19 at 14:25

My approach without using any built-in libraries.

Libraries used: pandas,numpy,argparse

import numpy as np
import pandas as pd
import argparse

eta,t = 0.01, 0.00001  #optimal hyperparameters

#********Matrix Calculations***************
def mat(dataset):
ds = np.genfromtxt(dataset,delimiter=',')

x = ds[:,0:-1].reshape(-1, ds.shape[1]-1)
ones = np.ones([x.shape[0],1])
x = np.concatenate([ones, x],1)

y = ds[:,-1].reshape(-1,1)

w = np.zeros([1, x.shape[1]])

return x, y, w

#********Calculating SSE********************
def sqrerror(x, y, w):
e = np.power(((x@w.T)-y),2)
return np.sum(e)

#********Implementing Batch Linear Regression*************
def regressor(dataset,eta,t):
x, y, w = mat(dataset)
costmat,res=[],[]
key=0
while True:
currres=[]
currres.append(key)
key += 1
cost = sqrerror(x,y,w)
for ol in w:
for il in ol:
currres.append(il)
currres.append(cost)
costmat.append(cost)
res.append(currres)

if key>1:
if (costmat[-2]-costmat[-1]) <= t:
break
return res

#************Driver code*************************
if __name__ == "__main__":

parser = argparse.ArgumentParser()
help='your_dataset.csv')

help='learning rate = 0.0001')


Works even with masked data.