gradient descent diverges extremely
Although this answer shows code with Python, the theory is exactly the same.
Two advices when running a linear regression via gradient descent from scratch:


use matrix notation for your features data (i.e. input data) so you can apply it to n-dimensional datasets (not just 1-D as in this case)
standardize your feature matrix data (although in this case you have actually 1-D, but it is a good practice which always works)


Your data plotted with matplotlib:

import numpy as np 
from numpy.random import randn, seed

x_data = 20 * randn(1000) + 100
y_data = x_data + (10 * randn(1000) + 50)
n=1000
t=1.4
sigma_R = t*0.001
min_value_t = t-sigma_R
max_value_t = t+sigma_R
y_data = min_value_t + (max_value_t - min_value_t) * np.random.rand(n,1)
x_data=np.array(range(1000))

pyplot.scatter(x_data, y_data)
pyplot.show()




Build your dataset to train the model:

final_df_DICT = {'X': x_data}
import pandas as pd
H = pd.DataFrame(final_df_DICT)
feature_matrix = np.zeros(n*2) 
feature_matrix.shape = (n, 2) 
feature_matrix[:,0] = 1 
feature_matrix[:,1] = H['X'] 
#standardize features
feature_matrix = (feature_matrix - feature_matrix.mean()) / feature_matrix.std()
target_data = y_data.reshape(len(y_data), )


Look at the feature matrix, having the first column filled with 1's and the second column having the actual x_data:


Model training with your hiperparameters:

w = [0, 0]
L=0.0001
epochs=1000
iteration = 0
cost=[]
while iteration &lt; epochs:
    pred = np.dot(feature_matrix, w)
    residuals = pred-target_data
    #we calculate the gradient for the 2 coeffs with the scalar product 
    for i in range(len(w)):
        partial = 2*np.dot(residuals, feature_matrix[:, i])
        w[i] = w[i] - L*partial

    iteration += 1
    computed_cost = np.sum(np.power((pred - target_data), 2)) / n

    cost.append(computed_cost)

print('coef: {}'.format(w))
print('cost: {}'.format(cost[-1]))


Result:


  coef: [-1.80963253e+00 -6.15189807e-06] cost: 6.466287828899486e-07


Let's plot the fitted regression model predictions over the original dataset (we trained on the whole dataset, not taking into account any validation set in this case...):

my_predictions = np.dot(feature_matrix, w)
pyplot.scatter(feature_matrix[:, 1], target_data)
pyplot.scatter(feature_matrix[:, 1], my_predictions, color='r')

pyplot.show()