(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$)

  • 1
    $\begingroup$ 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. $\endgroup$ – Sycorax Sep 18 '19 at 14:19
  • 1
    $\begingroup$ 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. $\endgroup$ – Sycorax Sep 18 '19 at 14:21
  • $\begingroup$ 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. $\endgroup$ – 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)
    while True:
        key += 1
        cost = sqrerror(x,y,w)
        for ol in w:
            for il in ol:
        w = w-(eta)*np.sum((x@w.T-y)*x,axis=0) #Gradient Calculation 

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

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

    parser = argparse.ArgumentParser()

                       help='learning rate = 0.0001')

                       help='threshold value = 0.0001')

    args = parser.parse_args()

    dataset = args.data
    learning_rate = float(args.eta)
    threshold = float(args.threshold)


Works even with masked data.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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