input: Write a python function to perform gradient descent for linear regression and track the cost function over iterations.
  assert gradient_descent(np.array([[1, 2], [1, 3], [1, 4], [1, 5]]), np.array([7, 6, 5, 4]), 0.01, 1000)[0][0] > 4.0
  assert len(gradient_descent(np.array([[1, 2], [1, 3], [1, 4], [1, 5]]), np.array([7, 6, 5, 4]), 0.01, 1000)[1]) == 1000
output: import numpy as np
  
  def gradient_descent(X, y, learning_rate, iterations):
      m, n = X.shape
      theta = np.zeros(n)
      cost_history = []
  
      for i in range(iterations):
          predictions = X.dot(theta)
          errors = predictions - y
          gradient = X.T.dot(errors) / m
          theta -= learning_rate * gradient
          cost = (1 / (2 * m)) * np.sum(errors ** 2)
          cost_history.append(cost)
          
      return theta, cost_history